def chi1d(h,energies = [0.],t=0.0001,delta=0.01,q=0.001,nk=1000,U=None,adaptive=True,ks=None):
  hkgen = h.get_hk_gen() # get the generating hamiltonian
  n = len(h.geometry.x) # initialice response
  m = np.zeros((len(energies),n*n),dtype=np.complex) # initialice
  if not adaptive: # full algorithm  
    if ks==None: ks = np.linspace(0.,1.,nk)  # create klist
    for k in ks:
#      print "Doing k=",k
      hk = hkgen(k) # fist point
      e1,ev1 = lg.eigh(hk) # get eigenvalues
      hk = hkgen(k+q) # second point
      e2,ev2 = lg.eigh(hk) # get eigenvalues
      ct = calculate_xychi(ev1.T,e1,ev2.T,e2,energies,t,delta) # contribution
      m += ct
    m = m/nk # normalice by the number of kpoints  
    ms = [m[i,:].reshape(n,n) for i in range(len(m))] # convert to matrices
  if adaptive: # adaptive algorithm 
    from integration import integrate_matrix
    def get_chik(k): # function which returns a matrix
      """ Get response at a cetain energy"""
      hk = hkgen(k) # first point
      e1,ev1 = lg.eigh(hk) # get eigenvalues
      hk = hkgen(k+q) # second point
      e2,ev2 = lg.eigh(hk) # get eigenvalues
      ct = calculate_xychi(ev1.T,e1,ev2.T,e2,energies,t,delta) # contribution
      return ct # return response
    ms = []
    m = integrate_matrix(get_chik,xlim=[0.,1.],eps=.1,only_imag=False) # add energy
    ms = [m[i,:].reshape(n,n) for i in range(len(m))] # convert to matrices

  if U==None: # raw calculation
    return ms 
  else: # RPA calculation
    return rpachi(ms,U=U)

def create_subspace(M, k):
    [size, images] = M.shape
    # calculate the mean
    mean = np.dot(M, np.ones((images, 1), dtype=np.int)) / images
    if images > size:
        covariance = np.dot((M - mean), (M - mean).T)
        [eigenvectors, eigenvalues] = la.eigh(covariance)

    # this should usually be the case since the number of pixels in a picture is probably
    # greater that the number of input pictures so instead of creating a huge Covariance
    # matrix which can be very large we instead calculate the eigenvectors of NxN matrix
    # and then use this to calculate the N eigenvectors of the DxD sized matrix
    else:
        L = np.dot((M - mean).T, (M - mean))
        [eigenvalues, eigenvectors] = la.eigh(L)
        eigenvectors = np.dot((M - mean), eigenvectors)
    # wow python no scoping in loops, it's kinda hard to take you serious as a language sometimes

    # to make the eigenvectors unit length or orthonormal
    for i in range(images):
        eigenvectors[:, i] = eigenvectors[:, i] / la.norm(eigenvectors[:, i])

    sorted_order = np.argsort(eigenvalues)
    sorted_order = np.flipud(sorted_order)

    eigenvalues = eigenvalues[sorted_order]
    eigenvectors = eigenvectors[:, sorted_order]

    principle_eigenvalues = eigenvalues[0:k]
    principle_eigenvectors = eigenvectors[:, 0:k]

    return principle_eigenvalues, principle_eigenvectors, mean

 def solve_kernel(self, regparam):
     self.regparam = regparam
     K1, K2 = self.K1, self.K2
     Y = self.Y.reshape((K1.shape[0], K2.shape[0]), order='F')
     #assert self.Y.shape == (self.K1.shape[0], self.K2.shape[0]), 'Y.shape!=(K1.shape[0],K2.shape[0]). Y.shape=='+str(Y.shape)+', K1.shape=='+str(self.K1.shape)+', K2.shape=='+str(self.K2.shape)
     if not self.trained:
         self.trained = True
         evals1, V  = la.eigh(K1)
         evals1 = mat(evals1).T
         V = mat(V)
         self.evals1 = evals1
         self.V = V
         
         evals2, U = la.eigh(K2)
         evals2 = mat(evals2).T
         U = mat(U)
         self.evals2 = evals2
         self.U = U
         self.VTYU = V.T * self.Y * U
     
     newevals = 1. / (self.evals1 * self.evals2.T + regparam)
     
     self.A = multiply(self.VTYU, newevals)
     self.A = self.V * self.A * self.U.T
     self.model = KernelPairwiseModel(self.A)

    def get_msig_pos( self, sctx, eps_app_eng, *args, **kw ):
        '''
        get biggest positive principle stress
        @param sctx:
        @param eps_app_eng:
        '''
        sig_eng, D_mtx = self.get_corr_pred( sctx, eps_app_eng, 0, 0, 0 )
        ms_vct = zeros( 3 )
        shape = sig_eng.shape[0]
        if shape == 3:
            s_mtx = self.map_sig_eng_to_mtx( sig_eng )
            m_sig, m_vct = linalg.eigh( s_mtx )

            # @todo: - this must be written in a more readable way
            # 
            if m_sig[-1] > 0:
                # multiply biggest positive stress with its vector                
                ms_vct[:2] = m_sig[-1] * m_vct[-1]
        elif shape == 6:
            s_mtx = self.map_sig_eng_to_mtx( sig_eng )
            m_sig = linalg.eigh( s_mtx )
            if m_sig[0][-1] > 0:
                # multiply biggest positive stress with its vector                
                ms_vct = m_sig[0][-1] * m_sig[1][-1]
        return ms_vct

def get_red_rho_A(Sx,Sz,b,N):
    '''
    Form the reduced ground state density matrix.
    Only defined when N >= 2.
    '''
    
    if N < 2:
        raise Exception("N must be greater or equal to 2.")
    
    H = get_tran_ising_H(Sx,Sz,b,N)
    E,V = eigh(H.toarray())
    
    rho_A_B = np.outer(V[:,0],V[:,0])
    l,basis = eigh(Sz.toarray())
    rho_A_k = rho_A_B
    for k in range(N-1,0,-1):
        rho_A_m = np.zeros([D**k,D**k])
        for m in range(D):
            # Rewrite basis states in the full N space.
            basis_m = get_full_matrix(dok_matrix(basis)[m].transpose(),k,k)
            basis_m = np.kron(np.eye(D),basis_m.toarray())
            
            # Trace the density matrix over the k-th particle.
            rho_A_m += np.dot(np.dot(np.transpose(basis_m),rho_A_k),basis_m)
        rho_A_k = rho_A_m
    rho_A = rho_A_k
    
    return rho_A

def HEG2(n,V):
    """\
    Does the diagonalization in the discrete variable representation
    """
    X = zeros((n,n),'d')
    for i in range(n):
        if i > 0:
            X[i,i-1] = X[i-1,i] = sqrt(i)/sqrt(2)
    # Eq 1 from HEG
    lam,T = eigh(X)

    KEho = zeros((n,n),'d')
    for i in range(n):
        KEho[i,i] = 0.25*(2*i+1)
        if i > 1:
            KEho[i,i-2] = KEho[i-2,i] = -0.25*sqrt(i-1)*sqrt(i)

    KEx = matmul(transpose(T),matmul(KEho,T))
    #KEx = matmul(T,matmul(KEho,transpose(T)))

    # Form the potential matrix
    # Eq 2 from HEG
    Vx = diag([V(li) for li in lam])

    Hx = KEx + Vx
    print "x\n",lam[:5]
    #from scipy.special.orthogonal import h_roots
    #print h_roots(n)
    
    matprint(KEx,label="T")
    matprint(Vx,label="V")
    matprint(Hx,label="H")
    
    E,U = eigh(Hx)
    return lam,E,U

def HEG(n,V):
    X = zeros((n,n),'d')
    for i in range(n):
        if i > 0:
            X[i,i-1] = X[i-1,i] = sqrt(i)/sqrt(2)
    # Eq 1 from HEG
    lam,T = eigh(X)
    print lam

    # Form the potential matrix
    # Eq 2 from HEG
    Vx = [V(li) for li in lam]
    Vho = matmul(T,matmul(diag(Vx),transpose(T)))

    KEho = zeros((n,n),'d')
    for i in range(n):
        KEho[i,i] = 0.25*(2*i+1)
        if i > 1:
            KEho[i,i-2] = KEho[i-2,i] = -0.25*sqrt(i-1)*sqrt(i)
    Hho = KEho+Vho
    Hx = matmul(transpose(T),matmul(Hho,T))
    E,U = eigh(Hho)

    # The eigenvectors are in terms of the HO eigenvectors, so
    # we have to multiply by X before returning
    return lam,E,matmul(transpose(T),U)

    def set(self,matrix):
        '''
        Set the basis.

        Parameters
        ----------
        matrix : callable
            The function to get the single particle matrix.
        '''
        Eup,Uup,Edw,Udw=[],[],[],[]
        for k in [()] if self.BZ is None else self.BZ.mesh('k'):
            m=matrix(k)
            es,us=nl.eigh(m[:m.shape[0]//2,:m.shape[0]//2])
            Edw.append(es)
            Udw.append(us)
            es,us=nl.eigh(m[m.shape[0]//2:,m.shape[0]//2:])
            Eup.append(es)
            Uup.append(us)
        Eup,Uup=np.asarray(Eup),np.asarray(Uup).transpose((1,0,2))
        Edw,Udw=np.asarray(Edw),np.asarray(Udw).transpose((1,0,2))
        if self.polarization=='up':
            self._E1_=Edw
            self._E2_=Eup
            self._U1_=Udw
            self._U2_=Uup
        else:
            self._E1_=Eup
            self._E2_=Edw
            self._U1_=Uup
            self._U2_=Udw

    def __calcEigChan(self,A1,G2,Left,channels=10):
        # Calculate Eigenchannels using recipe from PRB
        # For right eigenchannels, A1=A2, G2=G1 !!!
        if isinstance(A1,MM.SpectralMatrix):
            ev, U = LA.eigh(MM.mm(A1.L,A1.R))
        else:
            ev, U = LA.eigh(A1)

        # This small trick will remove all zero contribution vectors
        # and will diagonalize the tt matrix in the subspace where there
        # are values.
        idx = (ev > 0).nonzero()[0]
        ev = N.sqrt(ev[idx] / ( 2 * N.pi ))
        ev.shape = (1, -1)
        Utilde = ev * U[:, idx]
        
        nuo,nuoL,nuoR = self.nuo, self.nuoL, self.nuoR
        if Left:
            tt=MM.mm(MM.dagger(Utilde[nuo-nuoR:nuo,:]),2*N.pi*G2,Utilde[nuo-nuoR:nuo,:])
        else:
            tt=MM.mm(MM.dagger(Utilde[:nuoL,:]),2*N.pi*G2,Utilde[:nuoL,:])

        # Diagonalize (note that this is on a reduced tt matrix (no 0 contributing columns)
        evF, UF = LA.eigh(tt)
        EC = MM.mm(Utilde, UF[:,-channels:]).T
        return EC[::-1, :], evF[::-1] # reverse eigenvalues

    def diag_H(self):
        '''Diagonalize the Hamiltonians of spin a and spin b.
        '''
        e_a, w_a = nl.eigh(self.Ha)
        e_b, w_b = nl.eigh(self.Hb)
        e_gs = np.sum(e_a[:self.Na]) + np.sum(e_b[:self.Nb]) + self.C
        tmp_a = (w_a[:, :self.Na].dot(w_a.conj().T[:self.Na, :])).diagonal()
        print "tmpa:", tmp_a
        tmp_b = (w_b[:, :self.Nb].dot(w_b.conj().T[:self.Nb, :])).diagonal()
        print "tmpb:", tmp_b
        def vec_equal(v1, v2, Min = 10e-5):
            a = abs(v1-v2)
            a = (a < Min)*1.
            one = np.ones(len(v1))
            return np.array_equal(a, one)

        not_conv = True
        if vec_equal(tmp_a, self.lattice_a) and vec_equal(tmp_b,\
                self.lattice_b):
            not_conv = False
            return e_gs, not_conv

        self.lattice_a = tmp_a.copy() # update new density
        print "lattice a", self.lattice_a
        self.lattice_b = tmp_b.copy()
        print "lattice b", self.lattice_b
        return e_gs, not_conv #ground state energy

 def get_chik(k): # function which returns a matrix
   """ Get response at a cetain energy"""
   hk = hkgen(k) # first point
   e1,ev1 = lg.eigh(hk) # get eigenvalues
   hk = hkgen(k+q) # second point
   e2,ev2 = lg.eigh(hk) # get eigenvalues
   ct = calculate_xychi(ev1.T,e1,ev2.T,e2,energies,t,delta) # contribution
   return ct # return response

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 energy(s,n = None, scaleEnergy = 1):
  values = []
  if n == None:
    for q in numpy.arange(0,1.+qstep,qstep):
      evalues, evectors = linalg.eigh(ham(q,s,scaleEnergy))
      values.append(evalues[:5])
  else:
    for q in numpy.arange(0,1.+qstep,qstep):
      evalues, evectors = linalg.eigh(ham(q,s,scaleEnergy))
      values.append(evalues[n])
  return values

def test_eig_vs_eigh_above_560():
        # gh-6896
        N = 560

        A = np.arange(N*N).reshape(N, N)
        A = A + A.T

        w1 = np.sort(linalg.eig(A)[0])
        w2 = np.sort(linalg.eigh(A, UPLO='U')[0])
        w3 = np.sort(linalg.eigh(A, UPLO='L')[0])
        assert_array_almost_equal(w1, w2)
        assert_array_almost_equal(w1, w3)

def pca(X, d_prime):
    d,n = X.shape
    # mu: vector promedio
    mu = X.mean(axis=0)
    # Restamos la media 
    for i in range(n):
        X[i] -= mu 
    A = X.copy()

    if d>200 and n<3*d:
        if d_prime > n:
            d_prime = n
        # C: Matriz de covarianzas
        C_prime = 1.0/d * np.dot(A.T,A)
        #Delta=eigenvalues B=eigenvectors
        D_prime,B_prime = la.eigh(C_prime)
        #print "B prime: ", B_prime.shape, "- delta: ",  D_prime.shape

        for i in xrange(n):
            B_prime[:,i] = B_prime[:,i]/np.linalg.norm(B_prime[:,i])

        B = np.dot(A, B_prime)
        D = d/n * D_prime
        #print "B complete: ", B.shape, "- delta: ",  D.shape
        # Ordenamos los vectores propios, primero los que más varianza recogen 
        order = np.argsort(D, axis=0)[::-1] 
        # Ordenamos los vectores propios & los valores propios
        B = B[:,order]
        D = D[order]

    else:
        C = 1.0/n * np.dot(A,A.T)
        D,B = la.eigh(C) 
        # Ordenamos los vectores propios, primero los que más varianza recogen 
        order = np.argsort(D)[::-1] # sorting the eigenvalues
        # Ordenamos los vectores propios & los valores propios
        B = B[:,order]
        D = D[order]

    # B_dprime (d'xn)
    #print "B: ", B.shape, " - ", B[:,:d_prime].shape
    #print "D: ", D.shape
    #print "X: ", X.shape
    #print "d': ",d_prime
    #print "mu: ", mu.shape
    #Proyectamos los datos en d'
    B_dprime = B[:,:d_prime]
    y = np.dot(B_dprime.T,X)
    #print y[0]
    #print 
    #print
    #return ['B_dprime':B_dprime,D,B,mu,X]
    return {'B':B, 'B_dprime':B_dprime,'mu':mu,'y':y}, d_prime

def PCA(data_mat, p):
    'reduce the data dimensionality to p, this function will substract the \
    mean from the original data'
    d, N=data_mat.shape
    m=matlib.mean(data_mat, 1)
    data_mat-=m
    if d<N:
        AAT=data_mat*data_mat.T
        w, v=linalg.eigh(AAT)
        return v[:,-p:], m
    else:
        ATA=data_mat.T*data_mat
        w, v=linalg.eigh(ATA)
        return data_mat*v[:, -p:], m

def SVD(mat):
	matT = mat.transpose()
	matmatT = mat.dot(matT)
	matTmat = matT.dot(mat)
	egnvalU, egnvecU = LA.eigh(matmatT)
	egnvalV, egnvecV = LA.eigh(matTmat)
	V = np.fliplr(egnvecV)
	VT = V.transpose()
	egnvalV = egnvalV[::-1]
	S = np.zeros(mat.shape)
	for i in range(min(mat.shape)):
		S[i][i] = math.sqrt(egnvalV[i])
	U = np.dot(np.dot(mat,np.transpose(VT)),LA.pinv(S))
	return U, S, VT

 def get_msig_pm( self, sctx, eps_app_eng, *args, **kw ):
     sig_eng, D_mtx = self.get_corr_pred( sctx, eps_app_eng, 0, 0, 0 )
     t_field = zeros( 9 )
     shape = sig_eng.shape[0]
     if shape == 3:
         s_mtx = self.map_sig_eng_to_mtx( sig_eng )
         m_sig = linalg.eigh( s_mtx )
         if m_sig[0][-1] > 0:
             t_field[0] = m_sig[0][-1]#biggest positive stress
     elif shape == 6:
         s_mtx = self.map_sig_eng_to_mtx( sig_eng )
         m_sig = linalg.eigh( s_mtx )
         if m_sig[0][-1] > 0:
             t_field[0] = m_sig[0][-1]
     return t_field

    def calculate_edf(self, useibl=True):
        """Calculate the coefficients b_il in the expansion of the EDF.

        ``|phi_l> = sum_i b_il |f^u_i>``, in terms of ``|f^u_i> = P^u|f_i>``.

        To use the infinite band limit set useibl=True.
        N is the total number of bands to use.
        """
        
        for k, L in enumerate(self.L_k):
            if L==0:
                assert L!=0, 'L_k=0 for k=%i. Not implemented' % k
        
        self.Vo_kni = [V_ni[:M] for V_ni, M in zip(self.V_kni, self.M_k)]
        
        self.Fo_kii = np.asarray([np.dot(dagger(Vo_ni), Vo_ni) 
                                  for Vo_ni in self.Vo_kni])
        
        if useibl:
            self.Fu_kii = self.s_lcao_kii - self.Fo_kii
        else:
            self.Vu_kni = [V_ni[M:self.N] 
                           for V_ni, M in zip(self.V_kni, self.M_k)]
            self.Fu_kii = np.asarray([np.dot(dagger(Vu_ni), Vu_ni) 
                                     for Vu_ni in self.Vu_kni])
        self.b_kil = [] 
        for Fu_ii, L in zip(self.Fu_kii, self.L_k):
            b_i, b_ii = la.eigh(Fu_ii)
            ls = b_i.real.argsort()[-L:] 
            b_il = b_ii[:, ls] #pick out the eigenvec with largest eigenvals.
            normalize2(b_il, Fu_ii) #normalize the EDF: <phi_l|phi_l> = 1
            self.b_kil.append(b_il)

def _pca1 (X, verbose=False):
    """
    Simple principal component decomposition (PCA)
    X as Npix by Nt matrix
    X should be normalized and centered beforehand

    returns:
    - EV (Nt by Nesq:esq>0), matrix of PC 'signals'
     (eigenvalues of temporal covariance matrix). Signals are in columns
    - esq, vector of eigenvalues
    """
    print "Please don't use this, it's not ready"
    #return
    n_data, n_dimension = X.shape # (m x n)
    Y = X - X.mean(axis=0)[np.newaxis,:] # remove mean
    #C = dot(Y, Y.T) # (n x n)  covariance matrix
    C = dot(Y.T, Y)
    print C.shape
    es, EV = eigh(C)  # eigenvalues, eigenvectors

    ## take non-negative eigenvalues
    non_neg, = where(es>=0)
    neg = where(es<0)
    if len(neg)>0:
        if verbose:
            print "pca1: Warning, C have %d negative eigenvalues" %len(neg)
        es = es[non_neg]
        EV = EV[:,non_neg]
    #tmp = dot(Y.T, EV).T
    #V1 = tmp[::-1]
    #S1 = sqrt(es)[::-1]
    return EV