index to the first n value for a give m within the spherical 
    coefficients vector. Used by sc_to_fc
424296
def mindx(m, nmax, mmax):
    """"""

    ind = 0
    NN = nmax + 1

    if np.abs(m) > mmax:
        raise Exception("|m| cannot be larger than mmax")

    if (m != 0):
        ind = NN
        ii = 1
        for i in xrange(1, np.abs(m)):
            ind = ind + 2 * (NN - i)
            ii = i + 1

        if m > 0:
            ind = ind + NN - ii

    return ind
108202
def obfn_dfd(self):
        
        """

        Ef = sl.inner(self.Zf, self.obfn_fvarf(), axis=self.cri.axisM) \
          - self.Sf
        return (np.linalg.norm(self.W * sl.irfftn(Ef, self.cri.Nv,
                                                  self.cri.axisN))**2) / 2.0
949178
def sc_to_fc(spvec, nmax, mmax, nrows, ncols):
    """"""

    fdata = np.zeros([int(nrows), ncols], dtype=np.complex128)

    for k in xrange(0, int(ncols / 2)):
        if k < mmax:
            kk = k
            ind = mindx(kk, nmax, mmax)
            vec = spvec[ind:ind + nmax - np.abs(kk) + 1]
            fdata[:, kk] = fcvec_m_sc(vec, kk, nmax, nrows)

            kk = -(k + 1)
            ind = mindx(kk, nmax, mmax)
            vec = spvec[ind:ind + nmax - np.abs(kk) + 1]
            fdata[:, kk] = fcvec_m_sc(vec, kk, nmax, nrows)

        if k == mmax:
            kk = k
            ind = mindx(kk, nmax, mmax)
            vec = spvec[ind:ind + nmax - np.abs(kk) + 1]
            fdata[:, kk] = fcvec_m_sc(vec, kk, nmax, nrows)

    return fdata
1052128
def eye(n, d=None):
    
    c = _matrix.matrix()
    c.tt = _vector.vector()
    if d is None:
        n0 = _np.asanyarray(n, dtype=_np.int32)
        c.tt.d = n0.size
    else:
        n0 = _np.asanyarray([n] * d, dtype=_np.int32)
        c.tt.d = d
    c.n = n0.copy()
    c.m = n0.copy()
    c.tt.n = (c.n) * (c.m)
    c.tt.r = _np.ones((c.tt.d + 1,), dtype=_np.int32)
    c.tt.get_ps()
    c.tt.alloc_core()
    for i in xrange(c.tt.d):
        c.tt.core[
        c.tt.ps[i] -
        1:c.tt.ps[
              i +
              1] -
          1] = _np.eye(
            c.n[i]).flatten()
    return c
1295614
private boolean smo(double epsgr) {
        SupportVector i = svmax;
        int ii = gmaxindex;
        double old_alpha_i = i.alpha[ii];

        if (i.kcache == null) {
            gram(i);
        }

        SupportVector j = svmin;        
        int jj = gminindex;        
        double old_alpha_j = j.alpha[jj];
        
        // Second order working set selection.
        double best = 0.0;
        double gi = ii == 0 ? -i.g[0] : i.g[1];
        for (SupportVector v : sv) {
            double curv = i.k + v.k - 2 * kernel.k(i.x, v.x);
            if (curv <= 0.0) curv = TAU;
            
            double gj = -v.g[0];
            if (v.alpha[0] > 0.0 && gj < gi) {
                double gain = -Math.sqr(gi - gj) / curv;
                if (gain < best) {
                    best = gain;
                    j = v;
                    jj = 0;
                    old_alpha_j = j.alpha[0];
                }
            }
            
            gj = v.g[1];
            if (v.alpha[1] < v.C && gj < gi) {
                double gain = -Math.sqr(gi - gj) / curv;
                if (gain < best) {
                    best = gain;
                    j = v;
                    jj = 1;
                    old_alpha_j = j.alpha[1];
                }                
            }
        }
        
        if (j.kcache == null) {
            gram(j);
        }
        
        // Determine curvature
        double curv = i.k + j.k - 2 * kernel.k(i.x, j.x);
        if (curv <= 0.0) curv = TAU;

        if (ii != jj) {
            double delta = (-i.g[ii] - j.g[jj]) / curv;
            double diff = i.alpha[ii] - j.alpha[jj];
            i.alpha[ii] += delta;
            j.alpha[jj] += delta;

            if (diff > 0.0) {
                // Region III
                if (j.alpha[jj] < 0.0) {
                    j.alpha[jj] = 0.0;
                    i.alpha[ii] = diff;
                }
            } else {
                // Region IV
                if (i.alpha[ii] < 0.0) {
                    i.alpha[ii] = 0.0;
                    j.alpha[jj] = -diff;
                }
            }
            
            if (diff > i.C - j.C) {
                // Region I
                if (i.alpha[ii] > i.C) {
                    i.alpha[ii] = i.C;
                    j.alpha[jj] = i.C - diff;
                }
            } else {
                // Region II
                if (j.alpha[jj] > j.C) {
                    j.alpha[jj] = j.C;
                    i.alpha[ii] = j.C + diff;
                }
            }
        } else {

            double delta = (i.g[ii] - j.g[jj]) / curv;
            double sum = i.alpha[ii] + j.alpha[jj];
            i.alpha[ii] -= delta;
            j.alpha[jj] += delta;

            if (sum > i.C) {
                if (i.alpha[ii] > i.C) {
                    i.alpha[ii] = i.C;
                    j.alpha[jj] = sum - i.C;
                }
            } else {
                if (j.alpha[jj] < 0) {
                    j.alpha[jj] = 0;
                    i.alpha[ii] = sum;
                }
            }

            if (sum > j.C) {
                if (j.alpha[jj] > j.C) {
                    j.alpha[jj] = j.C;
                    i.alpha[ii] = sum - j.C;
                }
            } else {
                if (i.alpha[ii] < 0) {
                    i.alpha[ii] = 0.0;
                    j.alpha[jj] = sum;
                }
            }
        }
        
        double delta_alpha_i = i.alpha[ii] - old_alpha_i;
        double delta_alpha_j = j.alpha[jj] - old_alpha_j;

        int si = 2 * ii - 1;
        int sj = 2 * jj - 1;
        for (int k = 0; k < sv.size(); k++) {
            SupportVector v = sv.get(k);
            v.g[0] -= si * i.kcache.get(k) * delta_alpha_i + sj * j.kcache.get(k) * delta_alpha_j;
            v.g[1] += si * i.kcache.get(k) * delta_alpha_i + sj * j.kcache.get(k) * delta_alpha_j;
        }

        // optimality test
        minmax();
        b = -(gmax + gmin) / 2;
        if (gmax - gmin < epsgr) {
            return false;
        }

        return true;
    }
1233496
def matrix_chain_mult(M):
    
    opt, arg = matrix_mult_opt_order(M)
    return _apply_order(M, arg, 0, len(M)-1)
834511
def _n_parameters(self):
        """"""
        ndim = self.means_.shape[1]
        if self.covariance_type == 'full':
            cov_params = self.n_components * ndim * (ndim + 1) / 2.
        elif self.covariance_type == 'diag':
            cov_params = self.n_components * ndim
        elif self.covariance_type == 'tied':
            cov_params = ndim * (ndim + 1) / 2.
        elif self.covariance_type == 'spherical':
            cov_params = self.n_components
        mean_params = ndim * self.n_components
        return int(cov_params + mean_params + self.n_components - 1)
755042
def spin1x_from_xi1_phi_a_phi_s(xi1, phi_a, phi_s):
    
    phi1 = phi1_from_phi_a_phi_s(phi_a, phi_s)
    return xi1 * numpy.cos(phi1)
882706
def get_factors_iterative1(n):
    

    todo, res = [(n, 2, [])], []
    while todo:
        n, i, combi = todo.pop()
        while i * i <= n:
            if n % i == 0:
                res += combi + [i, n//i],
                todo.append((n//i, i, combi+[i])),
            i += 1
    return res
664698
private void moveAtFirst(int n_from) {
		int eT = elements[0];
		int eF = elements[n_from];
		positions[eT] = -1;
		positions[eF] = 0;
		elements[0] = eF;
		elements[n_from] = -1;
	}
621921
def _vmomentsurfaceIntegrand(vR,vT,R,az,df,n,m,sigmaR1,sigmaT1,t,initvmoment):
    """"""
    o= Orbit([R,vR*sigmaR1,vT*sigmaT1,az])
    return vR**n*vT**m*df(o,t)/initvmoment
2344
def solver(AA, N_max, symNx = 2, throw_out_modes=False):
    

    # Find all integer component n_vectors which lie within sphere of radius N_max
    # Here we have assumed that the potential is symmetric x->-x, y->-y, z->-z
    # This can be relaxed by changing symN to 1
    # Additionally due to time reversal symmetry S_n = -S_-n so we only consider
    # "half" of the n-vector-space

    angs = unroll_angles(AA.T[3:].T,np.ones(3))

    symNz = 2
    NNx = range(-N_max, N_max+1, symNx)
    NNy = range(-N_max, N_max+1, symNz)
    NNz = range(-N_max, N_max+1, symNz)
    n_vectors = np.array([[i,j,k] for (i,j,k) in product(NNx,NNy,NNz)
                          if(not(i==0 and j==0 and k==0)            # exclude zero vector
                             and (k>0                               # northern hemisphere
                                  or (k==0 and j>0)                 # half of x-y plane
                                  or (k==0 and j==0 and i>0))       # half of x axis
                             and np.sqrt(i*i+j*j+k*k)<=N_max)])     # inside sphere

    xxx = check_each_direction(n_vectors,angs)

    if(throw_out_modes):
        n_vectors = np.delete(n_vectors,check_each_direction(n_vectors,angs),axis=0)

    n = len(n_vectors)+3
    b = np.zeros(shape=(n, ))
    a = np.zeros(shape=(n,n))

    a[:3,:3]=len(AA)*np.identity(3)

    for i in AA:
        a[:3,3:]+=2.*n_vectors.T[:3]*np.cos(np.dot(n_vectors,i[3:]))
        a[3:,3:]+=4.*np.dot(n_vectors,n_vectors.T)*np.outer(np.cos(np.dot(n_vectors,i[3:])),np.cos(np.dot(n_vectors,i[3:])))
        b[:3]+=i[:3]
        b[3:]+=2.*np.dot(n_vectors,i[:3])*np.cos(np.dot(n_vectors,i[3:]))

    a[3:,:3]=a[:3,3:].T

    return np.array(solve(a,b)), n_vectors
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
78715
def one_hot(cls, n, l):
        
        v = [0.0] * l
        v[n] = 1.0
        return Vector(v)
364509
def ones(n, d=None):
    
    c = _vector.vector()
    if d is None:
        c.n = _np.array(n, dtype=_np.int32)
        c.d = c.n.size
    else:
        c.n = _np.array([n] * d, dtype=_np.int32)
        c.d = d
    c.r = _np.ones((c.d + 1,), dtype=_np.int32)
    c.get_ps()
    c.core = _np.ones(c.ps[c.d] - 1)
    return c
1267650
def kmm_init(X, m = 10):
    

    # compute the distances
    XXT = np.dot(X, X.T)
    D = (-2.*XXT + np.diag(XXT)[:,np.newaxis] + np.diag(XXT)[np.newaxis,:])

    # select the first point
    s = np.random.permutation(X.shape[0])[0]
    inducing = [s]
    prob = D[s]/D[s].sum()

    for z in range(m-1):
        s = np.random.multinomial(1, prob.flatten()).argmax()
        inducing.append(s)
        prob = D[s]/D[s].sum()

    inducing = np.array(inducing)
    return X[inducing]
703558
def smpar(C):
    """"""
    m2 = C['m2'].real
    Lambda = C['Lambda'].real
    v = (sqrt(2 * m2 / Lambda) + 3 * m2**(3 / 2) /
         (sqrt(2) * Lambda**(5 / 2)) * C['phi'])
    GF = 1 / (sqrt(2) * v**2)  # TODO
    Mh2 = 2 * m2 * (1 - m2 / Lambda * (3 * C['phi'] - 4 * Lambda * C['phiBox'] +
                                       Lambda * C['phiD']))
    eps = C['phiWB'] * (v**2)
    gb = (C['g'] / (1 - C['phiW'] * (v**2))).real
    gpb = (C['gp'] / (1 - C['phiB'] * (v**2))).real
    gsb = (C['gs'] / (1 - C['phiG'] * (v**2))).real
    MW = gb * v / 2
    ZG0 = 1 + C['phiD'] * (v**2) / 4
    MZ = (sqrt(gb**2 + gpb**2) / 2 * v
          * (1 + eps * gb * gpb / (gb**2 + gpb**2)) * ZG0)
    Mnup = -(v**2) * C['llphiphi']
    Mep = v / sqrt(2) * (C['Ge'] - C['ephi'] * (v**2) / 2)
    Mup = v / sqrt(2) * (C['Gu'] - C['uphi'] * (v**2) / 2)
    Mdp = v / sqrt(2) * (C['Gd'] - C['dphi'] * (v**2) / 2)
    UeL, Me, UeR = ckmutil.diag.msvd(Mep)
    UuL, Mu, UuR = ckmutil.diag.msvd(Mup)
    UdL, Md, UdR = ckmutil.diag.msvd(Mdp)
    UnuL, Mnu = ckmutil.diag.mtakfac(Mnup)
    eb = (gb * gpb / sqrt(gb**2 + gpb**2) *
          (1 - eps * gb * gpb / (gb**2 + gpb**2)))
    K = UuL.conj().T @ UdL
    # U = UeL.conj().T @ UnuL
    sm = {}
    sm['GF'] = GF
    sm['alpha_e'] = eb**2 / (4 * pi)
    sm['alpha_s'] = gsb**2 / (4 * pi)
    sm['Vub'] = abs(K[0, 2])
    sm['Vcb'] = abs(K[1, 2])
    sm['Vus'] = abs(K[0, 1])
    sm['delta'] = phase(-K[0, 0] * K[0, 2].conj()
                        / (K[1, 0] * K[1, 2].conj()))
    # sm['U'] = Uu
    sm['m_W'] = MW
    sm['m_Z'] = MZ
    sm['m_h'] = sqrt(abs(Mh2))
    sm['m_u'] = Mu[0]
    sm['m_c'] = Mu[1]
    sm['m_t'] = Mu[2]
    sm['m_d'] = Md[0]
    sm['m_s'] = Md[1]
    sm['m_b'] = Md[2]
    sm['m_e'] = Me[0]
    sm['m_mu'] = Me[1]
    sm['m_tau'] = Me[2]
    return  {k: v.real for k, v in sm.items()}