def get_left_channels(self, energy, nchan=1):
        self.initialize()
        g_s_ii = self.greenfunction.retarded(energy)
        lambda_l_ii = self.selfenergies[0].get_lambda(energy)
        lambda_r_ii = self.selfenergies[1].get_lambda(energy)

        if self.greenfunction.S is not None:
            s_mm = self.greenfunction.S
            s_s_i, s_s_ii = linalg.eig(s_mm)
            s_s_i = np.abs(s_s_i)
            s_s_sqrt_i = np.sqrt(s_s_i)  # sqrt of eigenvalues
            s_s_sqrt_ii = np.dot(s_s_ii * s_s_sqrt_i, dagger(s_s_ii))
            s_s_isqrt_ii = np.dot(s_s_ii / s_s_sqrt_i, dagger(s_s_ii))

        lambdab_r_ii = np.dot(np.dot(s_s_isqrt_ii, lambda_r_ii), s_s_isqrt_ii)
        a_l_ii = np.dot(np.dot(g_s_ii, lambda_l_ii), dagger(g_s_ii))
        ab_l_ii = np.dot(np.dot(s_s_sqrt_ii, a_l_ii), s_s_sqrt_ii)
        lambda_i, u_ii = linalg.eig(ab_l_ii)
        ut_ii = np.sqrt(lambda_i / (2.0 * np.pi)) * u_ii
        m_ii = 2 * np.pi * np.dot(np.dot(dagger(ut_ii), lambdab_r_ii), ut_ii)
        T_i, c_in = linalg.eig(m_ii)
        T_i = np.abs(T_i)

        channels = np.argsort(-T_i)[:nchan]
        c_in = np.take(c_in, channels, axis=1)
        T_n = np.take(T_i, channels)
        v_in = np.dot(np.dot(s_s_isqrt_ii, ut_ii), c_in)

        return T_n, v_in

def CholDecomp(amatrix):
    # Routine from "An iterative algorithm to produce
    # a positive definite correlation matrix from an
    # approximate correlation matrix" Iman and
    # Davenport, 1982 (Sandia report SAND81-1376)
    EigVal = linalg.eig(amatrix)[0]
    EigVec = linalg.eig(amatrix)[1]
    epsilon = np.array(EigVal <= 0, dtype="i") * 0.001
    if np.sum(epsilon) > 0:
        for i in range(10):
            EigVal = EigVal + epsilon
            EigValMat = np.diag(EigVal)
            amatrix = np.dot(EigVec, EigValMat)
            amatrix = np.dot(amatrix, transpose(EigVec))
            EigVal = linalg.eig(amatrix)[0]
            EigVec = linalg.eig(amatrix)[1]
            epsilon = np.array(EigVal <= 0, dtype="i") * 0.001
            if np.sum(epsilon) == 0:
                break
    decompmatrix = linalg.cholesky(amatrix)
    # need to scale decompmatrix so that diagonals equal 1
    # simply setting them to one is preferred--see Method A in
    # Iman and Davenport
    step = NVar + 1
    decompmatrix.flat[::step] = 1.0
    return decompmatrix

    def compute_score(self):
        self.scores = []
    
        # We now have a dictionary
        # start with a 'row' of all zeroes
        adjacency = []
        adjacency = adjacency + [0]*(len(self.user_dict) - len(adjacency))
        # Adjacency Matrix
        A = np.zeros( shape=(len(self.user_dict), len(self.user_dict)) )
        # keep track of A's rows
        outer_count = 0
        for mentioning_user in self.user_dict:
            inner_count = 0
            for mentioned_user in self.user_dict:
                if( mentioned_user in self.user_dict[mentioning_user]['mentioned'] ):
                    adjacency[inner_count] = 1
                else:
                    adjacency[inner_count] = 0
                inner_count += 1
            # print adjacency
            A[outer_count] = adjacency
            outer_count += 1

        self.scores = [np.dot(A, np.transpose(A)), np.dot(np.transpose(A), A)]
        dictList = []
        print "Hub:"
        w, v = LA.eig(np.dot(A, np.transpose(A)))
        i = np.real_if_close(w).argmax()
        principal = v[:,i]
        print self.user_dict.keys()[principal.argmax()]
        print "Authority:"
        w, v = LA.eig(np.dot(A, np.transpose(A)))
        i = np.real_if_close(w).argmax()
        principal = v[:,i]
        print self.user_dict.keys()[principal.argmax()]

    def _predict(self,k1,k2,y,gamma):
        la,Qa = LA.eig(k1)
        lb,Qb = LA.eig(k2)

        la = la.flatten()
        lb = lb.flatten()
        la = np.diag(la)
        lb = np.diag(lb)

        # http://stackoverflow.com/questions/17035767/kronecker-product-in-python-and-matlab
        diagLa = np.diag(la)
        diagLa = diagLa.reshape((len(diagLa),1))
        diagLbTrans = np.diag(lb).transpose()
        diagLbTrans = diagLbTrans.reshape((1,len(diagLbTrans)))

        l = sparse.kron( diagLbTrans,diagLa ).toarray()
        inverse = l / (l+gamma)

        m1 = Qa.transpose().dot(y).dot(Qb)
        m2 = m1 * inverse

        ypred = Qa.dot(m2).dot( Qb.transpose() )
        ypred = ypred.real

        return ypred

    def fit(self):
        from numpy.linalg import eig
        import scipy as sp
        import scipy.stats as stats

        n_samps, n_feats = self.shape()

        #Cross-Correlation and covariance matrix
        #Eigenvalues
        lcorr = eig (np.corrcoef(self._data.T))[0][::-1]
        lcov  = eig (np.cov     (self._data.T))[0][::-1]

        ems = []
        for i in range(len(self._far)):
            n_ems = 0
            pf    = self._far[i]
            for j in range(n_feats):
                sigma_sqr = (2*lcov[j]/n_samps) + (2*lcorr[j]/n_samps) + (2/n_samps) * lcov[j] * lcorr[j]
                sigma = sp.sqrt(sigma_sqr)

                print(sigma)
                # stats.norm.ppf not valid with sigma
                # using the module of the complex number : abs(sigma)
                tau = -stats.norm.ppf(pf, 0, abs(sigma))
                if (lcorr[j]-lcov[j]) > tau: 
                    n_ems += 1

            ems.append(n_ems)

        self.vd_  = ems

        return self.vd_

    def turn(self):
        return
        adp = self.adp['cart_int']
        adp = matrix([[float(adp[0]), float(adp[3]), float(adp[4])],
                      [float(adp[3]), float(adp[1]), float(adp[5])],
                      [float(adp[4]), float(adp[5]), float(adp[2])]])
        w, v = eig(adp)

        keep = w.tolist().index(min(w))
        vectors = [array((w[i] * v[:, i]).flatten().tolist()[0]) for i in range(3)]
        # print(vectors)

        value = 0
        for i in range(3):
            if not i == keep:
                value += w[i]
        value *= 0.5
        for i in range(3):
            if not i == keep:
                w[i] = value
        v = [array((w[i] * v[:, i]).flatten().tolist()[0]) for i in range(3)]
        # print vectors
        adp = matrix([[v[0][0], v[1][0], v[2][0]],
                      [v[0][1], v[1][1], v[2][1]],
                      [v[0][2], v[1][2], v[2][2]]])
        adp = (adp + adp.T) / 2
        # print adp

        w, v = eig(adp)

        keep = w.tolist().index(min(w))
        vectors = [array((w[i] * v[:, i]).flatten().tolist()[0]) for i in range(3)]

def propagateRLHamiltonian(t, k, omega, delta, epsilon, U, n):  
    t=np.array(t)    

    Energy1, V1 = LA.eig(RamanLatHamiltonian(0.0, 0.0, 0.0, 0.0, U, n))
    sort=np.argsort(Energy1)
    V1sorted=V1[:,sort]
    psi0=V1sorted[:,0]
   # psi0[np.divide(3*n,2)]=1.0+0.0*1j
    H = RamanLatHamiltonian(k, omega, delta ,epsilon,U,n)
    Energy, V = LA.eig(H)

    V = V + 1j*0.0
    Vinv = np.conjugate(np.transpose(V))

    # np.outer(t, Energy).flatten() creates a matrix for all t
    U = np.diag(np.exp(-1j*np.outer(t, Energy).flatten()))  

    a = np.dot(Vinv, psi0)
    # This repeats a so that the shape is consitent with U
    aa = np.outer(np.ones(t.size),a).flatten()
                      
    # Have to add the transpose to make shapes match 
    b = np.dot(U, aa)                                     
    # Same block diagonal trick for eigenvector matrix
    VV = sLA.block_diag(*([V]*t.size))                          
    psi = np.dot(VV, b)
    
    pops=np.absolute(psi)**2.0                     
    # Since you want the first value, need to take every 3rd row 
    # and extract the values you want from the diagonal
    latPops=np.sum(pops.reshape(t.size,n,3)[:,np.divide(n,2)-1:np.divide(n,2)+2,:],axis=2).flatten() 
    #populations in the -2k_L, 0, and +2k_L lattice sites, summed over spin sites,in time step blocks
    spinPops=np.sum(pops.reshape(t.size,n,3),axis=1).flatten() 
    #populations in each spin state, summed over lattice sites, in time step blocks 
    return spinPops

 def fit(self, x):
     self.matrix = x
     x = np.cov(x)
     ev = eig(x)[0]
     self.eg = eig(x)[1]
     self.ev = []
     for i in range(ev.shape[0]):
         self.ev.append([ev[i], i])
     self.ev[::-1].sort()

def linear_algebra():
    """ Use the `numpy.linalg` library to do Linear Algebra 
        For a reference on math, see 'Linear Algebra explained in four pages'
        http://minireference.com/static/tutorials/linear_algebra_in_4_pages.pdf
    """

    ### Setup two vectors
    x = np.array([1, 2, 3, 4])
    y = np.array([5, 6, 7, 8])

    ### Vector Operations include addition, subtraction, scaling, norm (length),
    # dot product, and cross product
    print np.vdot(x, y)  # Dot product of two vectors


    ### Setup two arrays / matrices
    a = np.array([[1, 2],
                  [3, 9]])
    b = np.array([[2, 4],
                  [5, 6]])


    ### Dot Product of two arrays
    print np.dot(a, b)


    ### Solving system of equations (i.e. 2 different equations with x and y)
    print LA.solve(a, b)


    ### Inverse of a matrix undoes the effects of the Matrix
    # The matrix multipled by the inverse matrix returns the 
    # 'identity matrix' (ones on the diagonal and zeroes everywhere else); 
    # identity matrix is useful for getting rid of the matrix in some equation
    print LA.inv(a)  # return inverse of the matrix
    print "\n"


    ### Determinant of a matrix is a special way to combine the entries of a
    # matrix that serves to check if matrix is invertible (!=0) or not (=0)
    print LA.det(a)  # returns the determinant of the array
    print "\n"  # e.g. 3, means that it is invertible


    ### Eigenvectors is a special set of input vectors for which the action of
    # the matrix is described as simple 'scaling'.  When a matrix is multiplied
    # by one of its eigenvectors, the output is the same eigenvector multipled
    # by a constant (that constant is the 'eigenvalue' of the matrix)
    print LA.eigvals(a)  # comput the eigenvalues of a general matrix
    print "\n"
    print LA.eigvalsh(a)  # Comput the eigenvalues of a Hermitian or real symmetric matrix
    print "\n"
    print LA.eig(a)  # return the eigenvalues for a square matrix
    print "\n"
    print LA.eigh(a)  # return the eigenvalues or eigenvectors of a Hermitian or symmetric matrix
    print "\n"

def plotEqqFqqA(streams, Q_t, alpha, p = 0):
    """
     p = plot e_qq and f_qq (YES/NO)
     flag = record all data for cov_mat and eigenvalues (YES/NO)
    """
    # N = number of timesteps + 1 for initial Q_0
    N = len(Q_t) 

    # Calculate F_qq #  (deviation fron orthogonality)
    f_qq = zeros((N,1))                                                                
    index = 0
    for q_t_i in Q_t:       
        X = dot(q_t_i.T , q_t_i) 
        FQQ = X - eye(X.shape[0])  
        f_qq[index, 0] = 10 * log10(trace(dot(FQQ.T, FQQ)))
        index += 1

    # Calculate E_qq (deviation from eigenvector subspace)
    e_qq = zeros((N-1,1))
    g_qq = zeros((N-1,1))
    cov_mat = zeros((streams.shape[1],streams.shape[1]))    
    for i in range(streams.shape[0]):
        
        data = streams[i,:]
        data = data.reshape(data.shape[0],1) # store as column vector 
        cov_mat = alpha * cov_mat + dot(data , data.T)
        W , V = eig(cov_mat)
                    
        # sort eigenVectors in according to deccending eigenvalue
        eig_idx = W.argsort() # Get sort index
        eig_idx = eig_idx[::-1] # Reverse order (default is accending)
        
        # v_r = highest r eigen vectors accoring to thier eigenvalue.
        V_r = V[:, eig_idx[:Q_t[i+1].shape[1]]]
        # Hopefuly have sorted correctly now 
        # Currently V_r is [1st 2nd, 3rd 4th] highest eigenvector 
        # according to eigenvalue     
        
        Y = dot(V_r , V_r.T) - dot(Q_t[i+1] , Q_t[i+1].T)  
        e_qq[i, 0] = 10 * log10(trace(dot(Y.T , Y)))
        
        # Calculate angle between projection matrixes
        A = dot(dot(dot(V_r.T , Q_t[i+1]) , Q_t[i+1].T) , V_r) 
        eigVal , eigVec = eig(A)
        angle = arccos(sqrt(max(eigVal)))        
        g_qq[i,0] = angle
        
    if p != 0:    
        figure()
        plot(f_qq)
        title('Deviation from orthonormality')    
        figure()
        plot(e_qq)
        title('Deviation of true tracked subspace') 
    
    return e_qq, f_qq, g_qq     

def Godunov_linear_solv(A , q_l , q_r , mode):

	dim = np.size(q_l)
	#Distinguish between 1dim case and system.
	#Case of a system: 
	if(dim > 1):
		eigenvalue , eigenvector = LA.eig(A)
		r = eigenvector
		eigenvalue , l = LA.eig(A.T)
		wavespeed_Godunov = eigenvalue * t_step / x_step	#Vector!
		wavespeed_LF = t_step / (x_step*2)					#Skalar!
		
		U = np.empty((x.size,t.size))
		#Sets the Q_i up for the first time step with the initial data. 
		#Q[j,i] is a matrix and contains the values for component i at x[j] at each time step
		q = initial_values( q_l , q_r , eigenvector , x )
		#iterating over time
		for j in range(np.size(t)) :
			#Values for the animation are saved in U
			U[:,j] = q[:,2]									#Change here to animate other components

			#Godunov
			if( mode == 1):
				q = update_Godunov(wavespeed_Godunov, q , x , l , r)
			
			#Lax-Friedrich		
			if( mode == 2):
				q = update_LF(wavespeed_LF, q , x , l , r , A)
				
			#Lax-Wendroff
			if( mode == 3):
				q = update_LW(wavespeed_LF, q , x , l , r , A)
	
	#1dim case with the wavespeed A, that gets passed instead of a matrix:	
	
	else:
		U = np.empty((x.size,t.size))
		q = np.zeros( np.size(x) )
		for i in range(np.size(q)) :
			q[i] = initial_values_1dim( x[i] )

		for j in range(np.size(t)) :
			U[:,j] = q
			qtemp = q
			if ( A > 0):
				for i in range(np.size(q)):
					if( i == 0 ):				qtemp[i] = q[i] - ((A * t_step / x_step) * (q[i] - q[np.size(q) - 1]))
					else:						qtemp[i] = q[i] - ((A * t_step / x_step) * (q[i] - q[i-1]))
			else:
				for i in range(np.size(q)):
					if( i == np.size(q) - 1):	qtemp[i] = q[i] - ((A * t_step / x_step) * (q[0] - q[i]))
					else:						qtemp[i] = q[i] - ((A * t_step / x_step) * (q[i+1] - q[i]))
			q = qtemp
	
	return U

def sorted_eigenvalues_vectors(matrix, hermitian=False):
    # i-th column(!) of v is eigenvector to i-th eigenvalue in w
    if hermitian:
        w,V = la.eigh(matrix)
    else:
        w,V = la.eig(matrix)
    w,V = la.eig(matrix)
    order = w.argsort()
    w = w[order]
    V = V[:,order]
    return w,V

def f(X):
    M = inv(X + .000001*np.eye(X.shape[0]))
    #return np.trace(M.dot(X))
    w, v = eig(M.dot(X))
    w_M, _ = eig(M)
    w_X, _ = eig(X)
    w.sort()
    w_M.sort()
    w_X.sort()

    print w[-5] - w_X[-5] * w_M[4]
    return w.max()

def get_bond_fc_with_sem(crds, fcmatrix, nat1, nat2, scalef, bondavg):

    crd1 = crds[3*nat1-3:3*nat1]
    crd2 = crds[3*nat2-3:3*nat2]
    disbohr = calc_bond(crd1, crd2) #unit is bohr
    dis = disbohr * B_TO_A #Transfer bohr to angstrom

    vec12 = array(crd2) - array(crd1) #vec12 is vec2 - vec1
    vec12 = [i/(disbohr) for i in vec12]
    vec12 = array(vec12)
 
    #bond force constant matrix, size 3 * 3
    bfcmatrix = array([[float(0) for x in range(3)] for x in range(3)])
   
    #1. First way to chose the matrix-----------------
    for i in range(0, 3):
      for j in range(0, 3):
        bfcmatrix[i][j] = -fcmatrix[3*(nat1-1)+i][3*(nat2-1)+j]
    eigval, eigvector = eig(bfcmatrix)
    fc = 0.0
    for i in range(0, 3):
      ev = eigvector[:,i]
      fc = fc + eigval[i] * abs(dot(ev, vec12))
    fcfinal1 = fc * HB2_TO_KCAL_MOL_A2 * 0.5

    if bondavg == 1:
      #2. Second way to chose the matrix-----------------
      for i in range(0, 3):
        for j in range(0, 3):
          bfcmatrix[i][j] = -fcmatrix[3*(nat2-1)+i][3*(nat1-1)+j]
      eigval, eigvector = eig(bfcmatrix)
      fc = 0.0
      for i in range(0, 3):
        ev = eigvector[:,i]
        fc = fc + eigval[i] * abs(dot(ev, vec12))
      fcfinal2 = fc * HB2_TO_KCAL_MOL_A2 * 0.5

      #Hatree/(Bohr^2) to kcal/(mol*angstrom^2)
      #Times 0.5 factor since AMBER use k(r-r0)^2 but not 1/2*k*(r-r0)^2

      fcfinal = average([fcfinal1, fcfinal2])
      stdv = std([fcfinal1, fcfinal2])
      fcfinal = fcfinal * scalef * scalef
      stdv = stdv * scalef * scalef
      return dis, fcfinal, stdv

    elif bondavg == 0:

      fcfinal = fcfinal1 * scalef * scalef
      return dis, fcfinal

def _parallelAnalysis(ff, n):

	""" Select the number of components for PCA using parallel analysis.
	
	Parameters
	----------
	ff : array_like
		Flat field data as numpy array. Each flat field is a single row 
		of this matrix, different rows are different observations.

	n : int
		Number of repetitions for parallel analysis.

	Return value
	------------
	V : array_like
		Eigen values.

	numPC : int
		Number of components for PCA.

	"""
	# Disable a warning:
	simplefilter("ignore", ComplexWarning)
	stdEFF = std(ff, axis=1, ddof=1)

	kpTrk = zeros((ff.shape[1], n), dtype=float32)
	stdMat = tile(stdEFF,(ff.shape[1], 1)).T

	for i in range(0, n):
		
		sample = stdMat * (randn(ff.shape[0], ff.shape[1])).astype(float32)		
		D, V = eig(cov(sample, rowvar=False))
		kpTrk[:,i] = sort(D).astype(float32)

	mean_ff_EFF = mean(ff,axis=1)
	
	F = ff - tile(mean_ff_EFF, (ff.shape[1], 1)).T
	D, V = eig(cov(F, rowvar=False))

	# Sort eigenvalues from smallest to largest:
	idx = D.argsort()   
	D = D[idx]
	V = V[:,idx]

	sel = zeros(ff.shape[1], dtype=float32)	
	sel[D > (mean(kpTrk, axis=1) + 2*std(kpTrk, axis=1, ddof=1))] = 1
	numPC = sum(sel).astype(int_)

	return (V, numPC)

def create(l,k,p): #esta funcion crea las quasienergias del operador de floquet, las quasienergias estan separadas en 
    x,y=LA.eig(R(l,J_y)) # un bloque de paridad positiva y otro de paridad negativa, al final regresa un vector que 
    x_sort=np.sort(x) #contiene la distribucion NNS de las quasienergias para ambos bloques si se desean solo 
    y_sort=y[:,x.argsort()] # las quasienergias hacer return dist_eig
    new_flo=chabas(flo(l,k,p),y_sort)
    pos_new_flo=pos_block(l,new_flo)
    neg_new_flo=neg_block(l,new_flo)
    eig_vals_pos,eig_vecs_pos=LA.eig(pos_new_flo)
    eig_vals_neg,eig_vecs_neg=LA.eig(neg_new_flo)
    r_pos,eig_ene_pos1=cart2pol(np.real(eig_vals_pos),np.imag(eig_vals_pos))
    r_neg,eig_ene_neg1=cart2pol(np.real(eig_vals_neg),np.imag(eig_vals_neg))
    eig_ene_pos=np.sort(eig_ene_pos1)
    eig_ene_neg=np.sort(eig_ene_neg1)
    eig_ene=np.append(eig_ene_pos,eig_ene_neg)
    return eig_ene

 def calcPCA(self, data):
     data -= np.mean(data, axis=0)
     # data = data / np.std(data, axis=0)
     c = np.cov(data, rowvar=0)
     values, vectors = la.eig(c)
     featureVector = vectors[:, [values.tolist().index(x) for x in np.sort(values)[::-1]]]
     return (np.matrix(featureVector) * np.matrix(data.T)).T

def myPCA(X,num_comp):
    nd = np.shape(X)
    n = nd[0]
    d = nd[1]
    # Verificar que el numero de comp. princ. sea menor a la dimension de X
    if (num_comp > d):
        print('Error: numero de componentes principales es mayor a la dimension de X')
    elif (num_comp <= 0):
        print('Error: numero de comp. principales menor o igual a 0')
    else:
        A = np.asmatrix(X)
        # Debemos restar de cada columna la media de los datos
        for i in range(d):
            A[0:,i] = np.asarray(A[0:,i]) - np.mean(A[0:,i])
        
        B = A.transpose()
        #Calculamos la matriz de covarainza
        M = B*A
        # Diagonaliza W son lo val propios ordenados y v la matriz de vect. propios corresp.
        w, v = LA.eig(M)
        # Define la matrix para la reduccion de dim
        T = v[:,0:num_comp]
        # transforma los datos
        T_X = T.transpose()*B
        return {'eigen_val':w, 'eigen_vec':v, 'T':T, 'TX':(T_X.transpose())}

def test_jordbloc():
    """Simple test of jordbloc. Do we recover the correct eigenvalues?"""
    evalue = 5.911  # arbitrary value
    a = rogues.jordbloc(10, evalue)
    w, v = nl.eig(a)
    b = evalue * np.ones(10)
    npt.assert_array_equal(w, b)

    def construct_rankarray(self):
        """Constructs the rank array"""
        from numpy import ones, argmax
        from numpy import linalg
        # M = (1 - m)*A + m*S
        # A = (1 - alpha)*B + alpha*C
        # B = normalised gain matrix
        # C = normalised intrinsic value matrix
        # S = gets rid of sub stochasticity for rows of all 0
        self.a_matrix = (1 - self.alpha)*self.g_matrix + \
            self.alpha*self.intrinsicvalue_array
        self.s_matrix = (1.0/self.size)*ones((self.size, self.size))
        mval = 0.15
        # Basic PageRank algorithm
        self.m_matrix = (1 - mval)*self.a_matrix + mval*self.s_matrix
        # Calculate eigenvalues, eigenvectors as usual
        [eigval, eigvec] = linalg.eig(self.m_matrix)

        maxeigindex = argmax(eigval)
        # Store value for downstream checking
        self.maxeig = eigval[maxeigindex].real
        # Cuts array into the eigenvector corrosponding to the eigenvalue above
        self.rank_array = eigvec[:, maxeigindex]
        # This is the 1-dimensional array composed of rankings (normalised)
        self.rank_array = (1/sum(self.rank_array))*self.rank_array
        # Remove the useless imaginary +0j
        self.rank_array = self.rank_array.real