def FLQTEB(engine,app):
    nmatrix=len(engine.generators['h'].table)
    if app.path!=None:
        result=zeros((app.path.rank,nmatrix+1))
        key=app.path.mesh.keys()[0]
        if len(app.path.mesh[key].shape)==1:
            result[:,0]=app.path.mesh[key]
        else:
            result[:,0]=array(xrange(app.path.rank[key]))
        for i,parameter in enumerate(list(app.path.mesh[key])):
            result[i,1:]=phase(eig(engine.evolution(t=app.ts.mesh['t'],**{key:parameter}))[0])/app.ts.volume['t']
    else:
        result=zeros((2,nmatrix+1))
        result[:,0]=array(xrange(2))
        result[0,1:]=angle(eig(engine.evolution(t=app.ts.mesh['t']))[0])/app.ts.volume['t']
        result[1,1:]=result[0,1:]
    if app.save_data:
        savetxt(engine.dout+'/'+engine.name.full+'_EB.dat',result)
    if app.plot:
        plt.title(engine.name.full+'_EB')
        plt.plot(result[:,0],result[:,1:])
        if app.show:
            plt.show()
        else:
            plt.savefig(engine.dout+'/'+engine.name.full+'_EB.png')

def SlowDownFactor(temporalnet):    
    """Returns a factor S that indicates how much slower (S>1) or faster (S<1)
    a diffusion process in the temporal network evolves on a second-order model 
    compared to a first-order model. This value captures the effect of order
    correlations on dynamical processes.
    """
    g2 = temporalnet.iGraphSecondOrder().components(mode="STRONG").giant()
    g2n = temporalnet.iGraphSecondOrderNull().components(mode="STRONG").giant()
    
    A2 = np.matrix(list(g2.get_adjacency()))
    T2 = np.zeros(shape=(len(g2.vs), len(g2.vs)))
    D2 = np.diag(g2.strength(mode='out', weights=g2.es["weight"]))
    
    for i in range(len(g2.vs)):
        for j in range(len(g2.vs)):
            T2[i,j] = A2[i,j]/D2[i,i]
    
    A2n = np.matrix(list(g2n.get_adjacency()))
    T2n = np.zeros(shape=(len(g2n.vs), len(g2n.vs)))
    D2n = np.diag(g2n.strength(mode='out', weights=g2n.es["weight"]))
    
    for i in range(len(g2n.vs)):
        for j in range(len(g2n.vs)):
            T2n[i,j] = A2n[i,j]/D2n[i,i]
    
    w2, v2 = spl.eig(T2, left=True, right=False)
    w2n, v2n = spl.eig(T2n, left=True, right=False)
    
    return np.log(np.abs(w2n[1]))/np.log(np.abs(w2[1]))

def test_streamline_tensors():
    # Small streamline
    streamline = [[1, 2, 3], [4, 5, 3], [5, 6, 3]]
    # Non-default eigenvalues:
    evals = [0.0012, 0.0006, 0.0004]
    streamline_tensors = life.streamline_tensors(streamline, evals=evals)
    npt.assert_array_almost_equal(streamline_tensors[0],
                                  np.array([[0.0009, 0.0003, 0.],
                                            [0.0003, 0.0009, 0.],
                                            [0., 0., 0.0004]]))

    # Get the eigenvalues/eigenvectors:
    eigvals, eigvecs = la.eig(streamline_tensors[0])
    eigvecs = eigvecs[np.argsort(eigvals)[::-1]]
    eigvals = eigvals[np.argsort(eigvals)[::-1]]

    npt.assert_array_almost_equal(eigvals,
                                  np.array([0.0012, 0.0006, 0.0004]))

    npt.assert_array_almost_equal(eigvecs[0],
                                  np.array([0.70710678, -0.70710678, 0.]))
    # Another small streamline
    streamline = [[1, 0, 0], [2, 0, 0], [3, 0, 0]]
    streamline_tensors = life.streamline_tensors(streamline, evals=evals)

    for t in streamline_tensors:
        eigvals, eigvecs = la.eig(t)
        eigvecs = eigvecs[np.argsort(eigvals)[::-1]]
        # This one has no rotations - all tensors are simply the canonical:
        npt.assert_almost_equal(np.rad2deg(np.arccos(
            np.dot(eigvecs[0], [1, 0, 0]))), 0)
        npt.assert_almost_equal(np.rad2deg(np.arccos(
            np.dot(eigvecs[1], [0, 1, 0]))), 0)
        npt.assert_almost_equal(np.rad2deg(np.arccos(
            np.dot(eigvecs[2], [0, 0, 1]))), 0)

def eig_linearised(Z, modes):
    """Solves a linearised approximation to the eigenvalue problem from
    the impedance calculated at some fixed frequency.

    The equation :math:`L = -s^2 S` is solved for `s`

    Parameters
    ----------
    Z : EfieImpedanceMatrixLoopStar
        The impedance matrix calculated in a loop-star basis
    modes : ndarray (int)
        A list or array of the mode numbers required

    Returns
    -------
    s_mode : ndarray, complex
        The resonant frequencies of the modes (in Hz)
        The complex pole `s` corresponding to the mode's eigenfrequency
    j_mode : ndarray, complex
        Columns of this matrix contain the corresponding modal currents
    """

    modes = np.asarray(modes)

    L = Z.matrices['L']
    S = Z.matrices['S']

    try:
        # Try to find the loop and star parts of the matrix (all relevant
        # matrices and vectors follow the same decomposition)
        loop, star = loop_star_indices(L)
    except AttributeError:
        loop = [[], []]
        star = [slice(None), slice(None)]

    if len(loop[0]) > 0 and len(loop[1]) > 0:
        L_conv = la.solve(L[loop[0], loop[1]],
                          L[loop[0], star[1]])
        L_red = (L[star[0], star[1]] -
                 np.dot(L[star[0], loop[1]], L_conv))

        # find eigenvalues, and star part of eigenvectors
        w, v_s = la.eig(S[star[0], star[1]], -L_red)

        vr = np.empty((L.shape[0], len(w)), np.complex128)
        vr[star[1]] = v_s
        vr[loop[1]] = -np.dot(L_conv, v_s)
    else:
        # Matrix does not have loop-star decomposition, so use the whole thing
        # TODO: implement some filtering to eliminate null-space solutions?
        w, vr = la.eig(S, -L)

    w_freq = np.sqrt(w)
    # make sure real part is negative
    w_freq = np.where(w_freq.real > 0, -w_freq, w_freq)

    w_selected = np.ma.masked_array(w_freq, abs(w_freq.real) > abs(w_freq.imag))
    which_modes = np.argsort(abs(w_selected.imag))[modes]

    return w_freq[which_modes], vr[:, which_modes]

    def algorithm3(self, e0, M):

        # Compute global maximum expansion rate
        Vtest = linspace(0,0.5,1000)
        muplot = []
        for Vi in Vtest:
            J = self.Jac(0,[Vi,0])
            muplot.append( 0.5*max(real(eig(J+J.T)[0])) )
        index = argmax(muplot)
        mustar = muplot[index]
        Vstar = Vtest[index]
        
        self.d1 = [e0]
        self.d2 = [e0]
        self.theta = [0]
        c = []
        for i in range(len(self.T)-1):
            # compute maximal expansion rate c_i in a neigbourhood of V[i] using global vector field bound M
            if abs(self.V[i] - Vstar) <= self.d1[i] + M*(self.T[i+1]-self.T[i]):
                c.append(mustar)
            elif self.V[i] + self.d1[i] + M*(self.T[i+1]-self.T[i]) < Vstar:
                J = Jac(0,[self.V[i] + self.d1[i] + M*(self.T[i+1]-self.T[i]),0])
                c.append(0.5*max(real(eig(J+J.T)[0])))
            else:
                J = Jac(0,[self.V[i] - self.d1[i] - M*(self.T[i+1]-self.T[i]),0])
                c.append(0.5*max(real(eig(J+J.T)[0])))            
            
            # compute diameter of ball based on bound on expansion rate in neighbourhood of current state
            self.d1.append(exp(c[i]*(self.T[i+1]-self.T[i]))*self.d1[i]+self.tolerance)
            self.d2.append(exp(c[i]*(self.T[i+1]-self.T[i]))*self.d2[i]+self.tolerance)
            self.theta.append(0)

    def coupled_modes(self, w, ignore_damping=False, **kwargs):
        M, B, C = self.linearised_matrices(w, **kwargs)
        if ignore_damping:
            wn, vn = linalg.eig(C, M)
            order = np.argsort(abs(wn))
            wn = np.sqrt(abs(wn[order]))
            vn = vn[:, order]
        else:
            AA = r_[c_[zeros_like(C), C], c_[C, B]]
            BB = r_[c_[C, zeros_like(C)], c_[zeros_like(C), -M]]
            wn, vn = linalg.eig(AA, BB)
            order = np.argsort(abs(wn))
            wn = abs(wn[order])
            # Mode shapes are the first half of the rows; the second
            # half of the rows should be same multiplied by eigenvalues.
            vn = vn[:M.shape[0], order]

            # We expect all the modes to be complex conjugate; return
            # every other one.
            # First: make sure all are the same sign
            norm_vn = vn / vn[np.argmax(abs(vn), axis=0), range(vn.shape[1])]
            assert (np.allclose(wn[::2], wn[1::2], rtol=1e-4) and
                    np.allclose(norm_vn[:, ::2], norm_vn[:, 1::2].conj(), atol=1e-2)), \
                "Expect conjugate modes"
            wn = wn[::2]
            vn = norm_vn[:, ::2]
        return wn, vn

def plot_ritz(A, n, iters):
    ''' Plot the relative error of the Ritz values of `A'.
    '''
    Amul = A.dot
    b = np.random.rand(A.shape[0])
    Q = np.empty((len(b), iters+1), dtype = np.complex128)
    H = np.zeros((iters+1, iters), dtype = np.complex128)
    Q[:, 0] = b / la.norm(b)
    eigvals = np.sort(abs(la.eig(A)[0]))[::-1]
    eigvals = eigvals[:n]
    abs_err = np.zeros((iters,n))

    for j in xrange(iters):
        Q[:, j+1] = Amul(Q[:, j])
        for i in xrange(j+1):
            H[i,j] = np.vdot(Q[:,i].conjugate(), (Q[:, j+1]))
            Q[:,j+1] = Q[:,j+1] - H[i,j] * (Q[:,i])

        H[j+1, j] = np.sqrt(np.vdot(Q[:, j+1], Q[:, j+1].conjugate()))
        Q[:,j+1] = Q[:,j+1] / H[j+1, j]

        if j < n:
            rit = np.zeros(n, dtype = np.complex128)
            rit[:j+1] = np.sort(la.eig(H[:j+1, :j+1])[0])[::-1]
            abs_err[j,:] = abs(eigvals - rit) / abs(eigvals)
        else:
            rit = np.sort(la.eig(H[:j+1,:j+1])[0])[::-1]
            rit = rit[:n]
            abs_err[j,:] = abs(eigvals - rit) / abs(eigvals)

    for i in xrange(n):
        plt.semilogy(abs_err[:,i])
    plt.show()

def pca(data, base_num=1):
    N, dim = data.shape

    data_m = data.mean(0)
    data_new = data - data_m

    # データ数 > 次元数
    if N > dim:
        # データ行列の共分散行列
        cov_mat = sp.dot(data_new.T, data_new) / float(N)
        # 固有値・固有ベクトルを計算
        l, vm = linalg.eig(cov_mat)
        # 固有値が大きい順に並び替え
        axis = vm[:, l.argsort()[-min(base_num, dim) :][::-1]].T

    # 次元数 > データ数
    else:
        base_num = min(base_num, N)
        cov_mat = sp.dot(data_new, data_new.T) / float(N)
        l, v = linalg.eig(cov_mat)
        # 固有値と固有ベクトルを並び替え
        idx = l.argsort()[::-1]
        l = l[idx]
        v = vm[:, idx]
        # 固有ベクトルを変換
        vm = sp.dot(data_m.T, v[:, :base_num])
        # （主成分の）基底を計算
        axis = sp.zeros([base_num, dim], dtype=sp.float64)
        for ii in range(base_num):
            if l[ii] <= 0:
                break
            axis[ii] = vm[:, ii] / linalg.norm(vm[:, ii])

    return axis

def FindMaximumQAQ(A, vertices, tetra):
  lambdas = []
  Q = np.zeros((4,4))
  for i in range(4):
    Q[:,i] = vertices[tetra[i],:]
  print "Q", Q
  # Full problem:
  A_ = Q.T.dot(A).dot(Q) 
  B_ = Q.T.dot(Q)  
  e, V = eig(A_, B_)
  alpha = np.real(V[:,np.argmax(e)])
  if np.all(alpha >= 0.) or np.all(alpha <= 0.):
    lambdas.append(np.max(np.real(e)))
  # Only three qs: 
  for comb in combinations(range(4), 3):
    A__ = np.array([[A_[i,j] for j in comb] for i in comb])
    B__ = np.array([[B_[i,j] for j in comb] for i in comb])
    e, V = eig(A__, B__)
    alpha = np.real(V[:,np.argmax(e)])
    if np.all(alpha >= 0.) or np.all(alpha <= 0.):
      lambdas.append(np.max(np.real(e)))
  # Only two qs: 
  for comb in combinations(range(4), 2):
    A__ = np.array([[A_[i,j] for j in comb] for i in comb])
    B__ = np.array([[B_[i,j] for j in comb] for i in comb])
    e, V = eig(A__, B__)
    alpha = np.real(V[:,np.argmax(e)])
    if np.all(alpha >= 0.) or np.all(alpha <= 0.):
      lambdas.append(np.max(np.real(e)))
  # Only one q: 
  for i in range(4):
    lambdas.append((Q[:,i]).T.dot(A).dot(Q[:,i]))
  print lambdas
  return np.max(np.array(lambdas))

    def test_create(self):
        basis_set = SphericalGTOSet()
        xyz = (0.0, 0.0, 0.0)
        r0 = 10.0
        for n in range(-10, 10):
            z = 2.0**n
            basis_set.add_one_basis(0, 0, xyz, z)
            for L in [0]:
                basis_set.add_basis(L, (0.0, 0.0, +r0), z)
                basis_set.add_basis(L, (0.0, 0.0, -r0), z)
                basis_set.add_basis(L, (0.0, +r0, 0.0), z)
                basis_set.add_basis(L, (0.0, -r0, 0.0), z)
                basis_set.add_basis(L, (+r0, 0.0, 0.0), z)
                basis_set.add_basis(L, (-r0, 0.0, 0.0), z)

        smat = basis_set.s_mat()
        for e in sorted(abs(la.eig(smat)[0]))[0:10]:
            print e

        print "-----"
        hmat = basis_set.t_mat() + basis_set.v_mat(1.0, xyz)
        zmat = basis_set.xyz_mat((0, 0, 1))
        for e in  sorted(la.eig(hmat, smat)[0].real)[0:10]:
            print e
        
        """

def regular_svd_by_pca(K, k=0):
    
    K_size = K.shape;
    if( K_size[0] < K_size[1] ):
        K_squared = np.dot(K, K.T);        
        tsp, tUp = la.eig(K_squared);
    else:
        K_squared = np.dot(K.T, K);
        tsp, tVp = la.eig(K_squared);
    # As la.eig returns complex number, use its absolute value.
    tsp = abs(tsp);
    tsp = np.sqrt(tsp);
    n_pos_sigs = sum(tsp > 0);
    tSp = np.diag(map(lambda s: 1.0/s, tsp[0:n_pos_sigs]));
    if( K_size[0] < K_size[1] ):
        tVp = np.dot(K.T, tUp);
        tVp[:, 0:n_pos_sigs] = np.dot(tVp[:, 0:n_pos_sigs], tSp);
    else:
        tUp = np.dot(K, tVp);
        tUp[:, 0:n_pos_sigs] = np.dot(tUp[:, 0:n_pos_sigs], tSp);
    if( 0 < k and k < min(K_size) ):
        tUp = tUp[:, 0:k];
        tVp = tVp[:, 0:k];
        tsp = tsp[0:k];
    return tUp, tsp, tVp;

def tridiag_eigs():
    # Most of this code is just constructing
    # tridiagonal matrices and calling functions
    # they have already written.
    m = 1000
    k = 100
    A = np.zeros((m, m))
    a = rand(m)
    b = rand(m-1)
    np.fill_diagonal(A, a)
    np.fill_diagonal(A[1:], b)
    np.fill_diagonal(A[:,1:], b)
    Amul = lambda u: tri_mul(a, b, u)
    alpha, beta = lanczos(rand(m), Amul, k)
    H = np.zeros((alpha.size, alpha.size))
    np.fill_diagonal(H, alpha)
    np.fill_diagonal(H[1:], beta)
    np.fill_diagonal(H[:,1:], beta)
    H_eigs = eig(H, right=False)
    H_eigs.sort()
    H_eigs = H_eigs[::-1]
    print H_eigs[:10]
    A = np.zeros((m, m))
    np.fill_diagonal(A, a)
    np.fill_diagonal(A[1:], b)
    np.fill_diagonal(A[:,1:], b)
    A_eigs = eig(A, right=False)
    A_eigs.sort()
    A_eigs = A_eigs[::-1]
    print A_eigs[:10]

def spatialFilter(Ra,Rb):
	R = Ra + Rb
	E,U = la.eig(R)

	# CSP requires the eigenvalues E and eigenvector U be sorted in descending order
	ord = np.argsort(E)
	ord = ord[::-1] # argsort gives ascending order, flip to get descending
	E = E[ord]
	U = U[:,ord]

	# Find the whitening transformation matrix
	P = np.dot(np.sqrt(la.inv(np.diag(E))),np.transpose(U))

	# The mean covariance matrices may now be transformed
	Sa = np.dot(P,np.dot(Ra,np.transpose(P)))
	Sb = np.dot(P,np.dot(Rb,np.transpose(P)))
	
	# Find and sort the generalized eigenvalues and eigenvector	# Find and sort the generalized eigenvalues and eigenvector

	E1,U1 = la.eig(Sa,Sb)
	ord1 = np.argsort(E1)
	ord1 = ord1[::-1]
	E1 = E1[ord1]
	U1 = U1[:,ord1]

	# The projection matrix (the spatial filter) may now be obtained
	SFa = np.dot(np.transpose(U1),P)
	return SFa.astype(np.float32)

 def __init__(self,matrix,error,overlap=None,overlap_err=None):
   if overlap_err is None:
     self.func = lambda mat:lin.eig(mat,overlap)[0]
     self.resample = lambda: gaussian_matrix_resample(matrix,error)
   else:
     self.func = lambda mats:lin.eig(mats[0],mats[1])[0]
     self.resample = lambda: (gaussian_matrix_resample(matrix,error),
                                gaussian_matrix_resample(overlap,overlap_err))

def eig(A,B):
	"""
	To ensure matlab compatibility, we need to
	swap matrices A and B around !!
	"""
	(XX1,XX2) = LIN.eig(A,B)
	(XX3,XX4) = LIN.eig(B,A)
	return (mat(XX4),mat(XX1))

    def compute_method(self,  parametro= None):
        ''' eigenvalues and eigen vectors '''
        #first era method
#         self.__eraresOut.A, self.__eraresOut.B, self.__eraresOut.C = mr.compute_ERA_model(np.array(self._signalOut), 2)
        # second, eigenvalues and eigenvectors
        self.__eraresOut.lambdaValues, self.__eraresOut.lambdaVector = linalg.eig(self.__eraresOut.A)
#         self.__eraresRef.A, self.__eraresRefB, self.__eraresRefC = mr.compute_ERA_model(np.array(self._signalRef), 2)
        # second, eigenvalues and eigenvectors
        self.__eraresRef.lambdaValues, self.__eraresRef.lamdaVector = linalg.eig(self.__eraresRef.A)

def designCSP(dataA, dataB, nb):

    # return v, a, d
    n_channels = dataA.shape[0]
    q = dataA.shape[1]

    cA = np.zeros([dataA.shape[0], n_channels, n_channels])
    cB = np.zeros([dataB.shape[0], n_channels, n_channels])

    # Compute the covariance matrix of each epoch of the same class (A and B)
    for i in range(dataA.shape[0]):
        # cA[i,...] = np.cov(dataA[i,:,:])
        c = np.dot(dataA[i, :, :], dataA[i, :, :].transpose())
        cA[i, ...] = c / (np.trace(c) * q)
        # cA[i,...] = c

    # compute the mean of the covariance matrices of each epoch
    cA_mean = cA.mean(0)

    for i in range(dataB.shape[0]):
        # cB[i,...] = np.cov(dataB[i,:,:])
        c = np.dot(dataB[i, :, :], dataB[i, :, :].transpose())
        cB[i, ...] = c / (np.trace(c) * q)
        # cB[i,...] = c

    # compute the mean of the covariance matrices of each epoch
    cB_mean = cB.mean(0)

    lamb, v = lg.eig(cA_mean + cB_mean)  # eigvalue and eigvector decomposition

    lamb = lamb.real  # return only real part of eigen vector

    # returns the index of array lamb in crescent order
    index = np.argsort(lamb)

    # reverse the order, now index has the positon of lamb in descendent order
    index = index[::-1]

    lamb = lamb[index]  # sort the eingenvalues in descendent order

    # the same goes for the eigenvectors along axis y
    v = v.take(index, axis=1)

    # whitening matrix computation
    Q = np.dot(np.diag(1 / np.sqrt(lamb)), v.transpose())

    # eig decomposition of whiten cov matrix
    D, V = lg.eig(np.dot(Q, np.dot(cA_mean, Q.transpose())))

    W_full = np.dot(V.transpose(), Q)

    # select only the neighbours defined in NB; get the first 3 eigenvectors
    W = W_full[:nb, :]
    W = np.vstack((W, W_full[-nb:, :]))  # get the three last eigenvectors

    return W

def printEigen(A, F):
    print 'Pole Locations:'
    (w_A, v_A) = la.eig(A)
    (w_F, v_F) = la.eig(F)
    for i in range(0, len(w_A)):
        print 'w_A = ', w_A[i]
        print 'v_A = \n', v_A[:,i]
    for i in range(0, len(w_F)):
        print 'w_F = ', w_F[i]
        print 'v_F = \n', v_F[:,i]
    return

def test_aligned_mem():
    """Check linalg works with non-aligned memory"""
    # Allocate 804 bytes of memory (allocated on boundary)
    a = arange(804, dtype=np.uint8)

    # Create an array with boundary offset 4
    z = np.frombuffer(a.data, offset=4, count=100, dtype=float)
    z.shape = 10, 10

    eig(z, overwrite_a=True)
    eig(z.T, overwrite_a=True)

def mysqrtm(m):
    m = 0.5 * (m.H + m)
    ls, vs = la.eigh(m)
    vs = np.matrix(vs)
    try:
        ls = [math.sqrt(max(l.real, 0)) for l in ls]
    except ValueError:
        print m.H - m
        print la.eig(m)
        raise ValueError
    return vs * np.diag(ls) * vs.H