Compute Hjorth Fractal Dimension of a time series X, kmax
     is an HFD parameter
487300
def hfd(X, Kmax):
    
    L = []
    x = []
    N = len(X)
    for k in range(1, Kmax):
        Lk = []
        for m in range(0, k):
            Lmk = 0
            for i in range(1, int(numpy.floor((N - m) / k))):
                Lmk += abs(X[m + i * k] - X[m + i * k - k])
            Lmk = Lmk * (N - 1) / numpy.floor((N - m) / float(k)) / k
            Lk.append(Lmk)
        L.append(numpy.log(numpy.mean(Lk)))
        x.append([numpy.log(float(1) / k), 1])

    (p, r1, r2, s) = numpy.linalg.lstsq(x, L)
    return p[0]
1164754
def K(self, X, X2, target):
        """"""
        if X2 is None: X2 = X
       # i1 = X[:,1]
       # i2 = X2[:,1]
       # X = X[:,0].reshape(-1,1)
       # X2 = X2[:,0].reshape(-1,1)
        dist = np.abs(X - X2.T)
        
        ly=1/self.lengthscaleY
        lu=np.sqrt(3)/self.lengthscaleU
        #ly=self.lengthscaleY
        #lu=self.lengthscaleU

        k1 = np.exp(-ly*dist)*(2*lu+ly)/(lu+ly)**2
        k2 = (np.exp(-lu*dist)*(ly-2*lu+lu*ly*dist-lu**2*dist) + np.exp(-ly*dist)*(2*lu-ly) ) / (ly-lu)**2 
        k3 = np.exp(-lu*dist) * ( (1+lu*dist)/(lu+ly) + (lu)/(lu+ly)**2 )

        np.add(self.varianceU*self.varianceY*(k1+k2+k3), target, target)
968974
def _K(self, R):
        """"""
        return self._ns * self._N / R / self._sin_alpha
1057752
def __f2d(frac_coords, v):
        
        # frac_coords = frac_coords % 1
        return np.array([int(frac_coords[0] * v.shape[0]),
                         int(frac_coords[1] * v.shape[1]),
                         int(frac_coords[2] * v.shape[2])])
1075066
def tau_davies(Z, fmax, t, y, dy, normalization='standard', dH=1, dK=3):
    """"""
    N = len(t)
    NH = N - dH  # DOF for null hypothesis
    NK = N - dK  # DOF for periodic hypothesis
    Dt = _weighted_var(t, dy)
    Teff = np.sqrt(4 * np.pi * Dt)
    W = fmax * Teff
    if normalization == 'psd':
        # 'psd' normalization is same as Baluev's z
        return W * np.exp(-Z) * np.sqrt(Z)
    elif normalization == 'standard':
        # 'standard' normalization is Z = 2/NH * z_1
        return (_gamma(NH) * W * (1 - Z) ** (0.5 * (NK - 1))
                * np.sqrt(0.5 * NH * Z))
    elif normalization == 'model':
        # 'model' normalization is Z = 2/NK * z_2
        return (_gamma(NK) * W * (1 + Z) ** (-0.5 * NK)
                * np.sqrt(0.5 * NK * Z))
    elif normalization == 'log':
        # 'log' normalization is Z = 2/NK * z_3
        return (_gamma(NK) * W * np.exp(-0.5 * Z * (NK - 0.5))
                * np.sqrt(NK * np.sinh(0.5 * Z)))
    else:
        raise NotImplementedError("normalization={0}".format(normalization))
1172202
def compute_k(self, memory_antecedent):
    
    if self.shared_kv:
      raise ValueError("compute_k cannot be called with shared_kv")
    ret = mtf.einsum(
        [memory_antecedent, self.wk], reduced_dims=[self.memory_input_dim])
    if self.combine_dims:
      ret = mtf.replace_dimensions(ret, ret.shape.dims[-1], self.k_dims)
    return ret
1082753
def calc_K_xx_wo_variance(self,X):
        """"""
        K_xx = np.ones([X.shape[0],X.shape[0]]) #ones now as a product occurs over each dimension
        for i,x in enumerate(X):
            for j,x2 in enumerate(X):
                for il,l in enumerate(self.lengthscale):
                    idx = il*2 #each pair of input dimensions describe the limits on one actual dimension in the data
                    K_xx[i,j] *= self.k_xx(x[idx],x2[idx],x[idx+1],x2[idx+1],l)
        return K_xx
384844
def fit(self, X):
        
        self.centers_, self.labels_, self.sse_arr_, self.n_iter_ = \
              _kmeans(X, self.n_clusters, self.max_iter, self.n_trials, self.tol)
1233241
def _hotellings_fourier(self):
        
        sigma_2 = self._catalog_object.sigma**2
        # compute residual
        df = self._df #degrees of freedom
        R = self._A - np.dot(self._X, self._Bhat)
        R = np.matrix(R)
        # residual covariance matrix
        Sigma_R = (R.H*R) * (1./df)
        # noise covariance matrix
        #raise TypeError("fourier domain hypothesis tests not implemented yet")
        Zft = np.matrix(self._basis_object._Z)
        Sigma_S = Zft.H*np.matrix(np.diag(sigma_2))*Zft
        # add covariances (sum of multivariate normals is normal)
        Sigma_Z = np.array(Sigma_R + Sigma_S)

        T_2_list = []
        p_value_list = []
        for i in np.arange(0,self._Bhat.shape[0]):

            Bstar       = self._Bhat[i,np.arange(0,np.shape(self._A)[1])]
            lstar       = float(np.shape(Bstar)[0])
            cx          = self._Cx[i,i]
            Einv        = np.linalg.pinv(Sigma_Z)
            Zs          = Bstar/np.sqrt(cx)
            last_part   = np.array(np.matrix(Zs)*np.matrix(Einv)*np.matrix(Zs).H).real
            T_2         = ((df - lstar + 1.)/(df*lstar))*last_part
            p_value     = 1. - stats.f.cdf(T_2, 2.*lstar, 2.*(df - lstar + 1.))
            p_value_list.append(p_value)
            T_2_list.append(T_2)
        # save pvalue results in a table
        results = [['Comparison','Hotellings T^2', "p-value"]]
        for i in np.arange(0,len(self._col_names)):
            results.append([self._col_names[i], T_2_list[i], p_value_list[i]])
        self._results = results
128702
def fRD( a, M):
    """"""
    f = (lal.C_SI**3.0 / (2.0*lal.PI*lal.G_SI*M*lal.MSUN_SI)) * (1.5251 - 1.1568*(1.0-a)**0.1292)
    return f
155363
def poisson_se_multiclust(data, k, n_runs=10, **se_params):
    
    clusters = []
    norm_data = cell_normalize(data)
    if sparse.issparse(data):
        log_data = data.log1p()
        log_norm = norm_data.log1p()
    else:
        log_data = np.log1p(data)
        log_norm = np.log1p(norm_data)
    tsvd_50 = TruncatedSVD(50)
    tsvd_k = TruncatedSVD(k)
    km = KMeans(k)
    tsvd1 = tsvd_50.fit_transform(data.T)
    tsvd2 = tsvd_k.fit_transform(data.T)
    tsvd3 = tsvd_50.fit_transform(log_data.T)
    tsvd4 = tsvd_k.fit_transform(log_data.T)
    tsvd5 = tsvd_50.fit_transform(norm_data.T)
    tsvd6 = tsvd_k.fit_transform(norm_data.T)
    tsvd7 = tsvd_50.fit_transform(log_norm.T)
    tsvd8 = tsvd_k.fit_transform(log_norm.T)
    tsvd_results = [tsvd1, tsvd2, tsvd3, tsvd4, tsvd5, tsvd6, tsvd7, tsvd8]
    clusters = []
    for t in tsvd_results:
        clust = km.fit_predict(t)
        clusters.append(clust)
    clusterings = np.vstack(clusters)
    consensus = CE.cluster_ensembles(clusterings, verbose=False, N_clusters_max=k)
    init_m, init_w = nmf_init(data, consensus, k, 'basic')
    M, W, ll = poisson_estimate_state(data, k, init_means=init_m, init_weights=init_w, **se_params)
    return M, W, ll
694568
def lm_ffinal(freqs, damping_times, modes):
    

    f_max = {}
    for lmn in modes:
        l, m, nmodes = int(lmn[0]), int(lmn[1]), int(lmn[2])
        for n in range(nmodes):
            f_max['%d%d%d' %(l,m,n)] = qnm_freq_decay(freqs['%d%d%d' %(l,m,n)],
                                      damping_times['%d%d%d' %(l,m,n)], 1./1000)

    f_final = max(f_max.values())
    if f_final > max_freq:
        f_final = max_freq

    return f_final
461528
def _MStep(self, l, K, x, z_h1, pi_h1, p_h1, p_h, delta_i_j_s):
        
        for k in range(K):
            normconst = 0
            if l == 0: # only need to compute pi at first MM iteration
                pi_h1[k] = np.sum(z_h1.T[k]) / len(z_h1)
            for j in range(self.m):
                omega_k_j = EMMMixPLAggregator.omega(k, j, z_h1, x) # numerator
                denom_sum = 0
                for i in range(self.n):
                    sum1 = 0
                    for t in range(len(x[i])):
                        sum2 = 0
                        sum3 = 0
                        for s in range(t, self.m):
                            sum2 += p_h[k][EMMMixPLAggregator.c(x[i], s)]
                        for s in range(t, self.m + 1):
                            sum3 += delta_i_j_s[i][j][s]
                        sum1 += z_h1[i][k] * (sum2 ** -1) * sum3
                    denom_sum += sum1
                p_h1[k][j] = omega_k_j / denom_sum
                normconst += p_h1[k][j]
            for j in range(self.m):
                p_h1[k][j] /= normconst
401210
def ks_synth(freq):
  
  ks_mem = (sum(lz.sinusoid(x * freq) for x in [1, 3, 9]) +
            lz.white_noise() + lz.Stream(-1, 1)) / 5
  return lz.karplus_strong(freq, memory=ks_mem)
923293
def to_evenly_discretized_mfd(self):
        
        return EvenlyDiscretizedMFD(self.mmin + self.bin_width / 2.,
                                    self.bin_width,
                                    self.occurrence_rate.tolist())
31838
def timeseries(X, **kwargs):
    """"""
    pl.figure(figsize=(2*rcParams['figure.figsize'][0], rcParams['figure.figsize'][1]),
              subplotpars=sppars(left=0.12, right=0.98, bottom=0.13))
    timeseries_subplot(X, **kwargs)
1248223
def Q_srk(self, k):
        
        ind = self.index[self.Q_time_var_index, k]
        Q = self.Q[:, :, ind]

        if (Q.shape[0] == 1):  # 1-D case handle simplier. No storage
            # of the result, just compute it each time.
            square_root = np.sqrt(Q)
        else:
            if self.svd_each_time:

                (U, S, Vh) = sp.linalg.svd(Q, full_matrices=False,
                                           compute_uv=True,
                                           overwrite_a=False,
                                           check_finite=True)

                square_root = U * np.sqrt(S)
            else:

                if ind in self.Q_square_root:
                    square_root = self.Q_square_root[ind]
                else:
                    (U, S, Vh) = sp.linalg.svd(Q, full_matrices=False,
                                               compute_uv=True,
                                               overwrite_a=False,
                                               check_finite=True)

                    square_root = U * np.sqrt(S)

                    self.Q_square_root[ind] = square_root

        return square_root