# Path: nrgmap/utils.py
# def H2G(h,w,tp='r',geta=1e-2,sigma=None):
#     '''
#     Get Green's function g from Hamiltonian h.
# 
#     Parameters:
#         :h: An array of hamiltonian.
#         :w: The energy(frequency).
#         :tp: The type of Green's function.
# 
#             * 'r': retarded Green's function.(default)
#             * 'a': advanced Green's function.
#             * 'matsu': finite temperature Green's function.
#         :geta: Smearing factor. default is 1e-2.
#         :sigma: Additional self energy.
# 
#     Return:
#         A(Array of) Green's function.
#     '''
#     if tp=='r':
#         z=w+1j*geta
#     elif tp=='a':
#         z=w-1j*geta
#     elif tp=='matsu':
#         z=1j*w
#     if sigma!=None:
#         h=h+sigma
#     if ndim(h)>0:
#         return inv(z*identity(h.shape[-1])-h)
#     else:
#         return 1./(z-h)
# 
# def s2vec(s):
#     '''
#     Transform a 2 x 2 matrix to a 4 dimensional vector, corresponding to s0,sx,sy,sz component.
# 
#     Parameters:
#         :s: The matrix.
#     Return:
#         len-4 vector indicating the pauli components.
#     '''
#     res=array([trace(s),trace(dot(sx,s)),trace(dot(sy,s)),trace(dot(sz,s))])/2
#     return res
# 
# Path: nrgmap/chain.py
# class Chain(object):
#     '''
#     Wilson chain class.
#     
#     Attributes:
#         :elist/tlist: 2D array/4D array, A list of on-site energies and coupling terms.
#         :nband: integer, the number of bath sites.(readonly)
#         :nsite: integer, the number of bands.(readonly)
#     '''
#     def __init__(self,elist,tlist):
#         if ndim(elist)!=3 or ndim(tlist)!=3:
#             raise ValueError()
#         self.elist=complex128(elist)
#         self.tlist=complex128(tlist)
# 
#     @property
#     def nsite(self):
#         '''The number of sites.'''
#         return len(self.elist)
# 
#     @property
#     def nband(self):
#         '''number of bands.'''
#         return self.tlist.shape[-1]
# 
#     def save(self,file_prefix):
#         '''
#         save a Chain instance to files.
# 
#         Parameters:
#             :file_prefix: A string as a prefix to store datas of a chain.
# 
#         **Note:**
#         The data is stored in 3 files,
# 
#         1. <file_prefix>.info.dat, the true shape of elist.
#             e.g. for multi-band, it's (# of sites, degree of freedom, degree of freedom).
#                  for single-band, it's (# of sites).
#         2. <file_prefix>.el.dat, an array of on-site energies, in the format(2-band,2-site as an example)
#             | E0[0,0].real E0[0,0].imag
#             | E0[0,1].real E0[0,1].imag
#             | E0[1,0].real E0[1,0].imag
#             | E0[0,1].real E0[0,1].imag
#             | E1[0,0].real E1[0,0].imag
#             |          ... ...
#         3. <file_prefix>.tl.dat, an array of hopping terms with the first element t0 the coupling strength with impurity.
#             Its data format is similar to `el.dat`.
#         '''
#         savetxt(file_prefix+'.info.dat',array(self.elist.shape))
#         savetxt(file_prefix+'.el.dat',self.elist.ravel().view('float64'))
#         savetxt(file_prefix+'.tl.dat',self.tlist.ravel().view('float64'))
# 
#     def get_H0(self,e0):
#         '''
#         Hamiltonian without interaction terms.
# 
#         Return:
#             sparse matrix, the hamiltonian.
#         '''
#         tlist=self.tlist
#         elist=self.elist
#         N=self.nsite+1
#         B=ndarray([N,N],dtype='O')
#         elist=concatenate([e0[newaxis,...],elist],axis=0)
#         tlistH=swapaxes(tlist,1,2).conj()
#         offset=[-1,0,1]
#         return construct_tridmat([tlist,elist,tlistH],offset).toarray()


from numpy import *
from scipy.sparse import block_diag
from scipy.linalg import eigvalsh,norm,qr
from scipy.integrate import trapz
from matplotlib.pyplot import *
from matplotlib import cm
from .utils import H2G,s2vec
from .tridiagonalize import *
from .chain import Chain
from .lib.futils import hybri_chain
import pdb
            if normalize_method=='mpqr':
                data,offset=tridiagonalize_mpqr(H,q=qq.reshape([-1,model.nband]),m=nsite)
            else:
                data,offset=tridiagonalize_qr(H,q=qq.reshape([-1,model.nband]),m=nsite)
            chain=Chain(asarray(data[1]),concatenate([t0[newaxis],data[0]],axis=0))
        cl.append(chain)
    return cl

def check_spec(chains,rhofunc,wlist,mode='eval',smearing=1.):
    '''
    Check mapping quality for wilson chain.

    Parameters:
        :chains: list, a list of chain after mapping.
        :rhofunc: function, Hybridization function.
        :wlist: 1D array, the frequency space.
        :mode: str, Choose the checking method,
            * 'eval' -> check eigenvalues.
            * 'pauli' -> check pauli components, it is only valid for 2 band system.

        :smearing: The smearing factor.
    '''
    nband=chains[0].nband
    nz=len(chains)
    print('Recovering Spectrum ...')
    dlv=hybri_chain(tlist=array([chain.tlist for chain in chains]),elist=array([chain.elist for chain in chains]),wlist=wlist,smearing=smearing*1./nz)
    if nband==1:
        dlv=dlv[:,0,0]
        nplt=1
    elif nband==2 and mode=='pauli':

# Path: nrgmap/utils.py
# def H2G(h,w,tp='r',geta=1e-2,sigma=None):
#     '''
#     Get Green's function g from Hamiltonian h.
# 
#     Parameters:
#         :h: An array of hamiltonian.
#         :w: The energy(frequency).
#         :tp: The type of Green's function.
# 
#             * 'r': retarded Green's function.(default)
#             * 'a': advanced Green's function.
#             * 'matsu': finite temperature Green's function.
#         :geta: Smearing factor. default is 1e-2.
#         :sigma: Additional self energy.
# 
#     Return:
#         A(Array of) Green's function.
#     '''
#     if tp=='r':
#         z=w+1j*geta
#     elif tp=='a':
#         z=w-1j*geta
#     elif tp=='matsu':
#         z=1j*w
#     if sigma!=None:
#         h=h+sigma
#     if ndim(h)>0:
#         return inv(z*identity(h.shape[-1])-h)
#     else:
#         return 1./(z-h)
# 
# def s2vec(s):
#     '''
#     Transform a 2 x 2 matrix to a 4 dimensional vector, corresponding to s0,sx,sy,sz component.
# 
#     Parameters:
#         :s: The matrix.
#     Return:
#         len-4 vector indicating the pauli components.
#     '''
#     res=array([trace(s),trace(dot(sx,s)),trace(dot(sy,s)),trace(dot(sz,s))])/2
#     return res
# 
# Path: nrgmap/chain.py
# class Chain(object):
#     '''
#     Wilson chain class.
#     
#     Attributes:
#         :elist/tlist: 2D array/4D array, A list of on-site energies and coupling terms.
#         :nband: integer, the number of bath sites.(readonly)
#         :nsite: integer, the number of bands.(readonly)
#     '''
#     def __init__(self,elist,tlist):
#         if ndim(elist)!=3 or ndim(tlist)!=3:
#             raise ValueError()
#         self.elist=complex128(elist)
#         self.tlist=complex128(tlist)
# 
#     @property
#     def nsite(self):
#         '''The number of sites.'''
#         return len(self.elist)
# 
#     @property
#     def nband(self):
#         '''number of bands.'''
#         return self.tlist.shape[-1]
# 
#     def save(self,file_prefix):
#         '''
#         save a Chain instance to files.
# 
#         Parameters:
#             :file_prefix: A string as a prefix to store datas of a chain.
# 
#         **Note:**
#         The data is stored in 3 files,
# 
#         1. <file_prefix>.info.dat, the true shape of elist.
#             e.g. for multi-band, it's (# of sites, degree of freedom, degree of freedom).
#                  for single-band, it's (# of sites).
#         2. <file_prefix>.el.dat, an array of on-site energies, in the format(2-band,2-site as an example)
#             | E0[0,0].real E0[0,0].imag
#             | E0[0,1].real E0[0,1].imag
#             | E0[1,0].real E0[1,0].imag
#             | E0[0,1].real E0[0,1].imag
#             | E1[0,0].real E1[0,0].imag
#             |          ... ...
#         3. <file_prefix>.tl.dat, an array of hopping terms with the first element t0 the coupling strength with impurity.
#             Its data format is similar to `el.dat`.
#         '''
#         savetxt(file_prefix+'.info.dat',array(self.elist.shape))
#         savetxt(file_prefix+'.el.dat',self.elist.ravel().view('float64'))
#         savetxt(file_prefix+'.tl.dat',self.tlist.ravel().view('float64'))
# 
#     def get_H0(self,e0):
#         '''
#         Hamiltonian without interaction terms.
# 
#         Return:
#             sparse matrix, the hamiltonian.
#         '''
#         tlist=self.tlist
#         elist=self.elist
#         N=self.nsite+1
#         B=ndarray([N,N],dtype='O')
#         elist=concatenate([e0[newaxis,...],elist],axis=0)
#         tlistH=swapaxes(tlist,1,2).conj()
#         offset=[-1,0,1]
#         return construct_tridmat([tlist,elist,tlistH],offset).toarray()


# Path: nrgmap/chainmapper.py
from numpy import *
from scipy.sparse import block_diag
from scipy.linalg import eigvalsh,norm,qr
from scipy.integrate import trapz
from matplotlib.pyplot import *
from matplotlib import cm
from .utils import H2G,s2vec
from .tridiagonalize import *
from .chain import Chain
from .lib.futils import hybri_chain
import pdb

            if normalize_method=='mpqr':
                data,offset=tridiagonalize_mpqr(H,q=qq.reshape([-1,model.nband]),m=nsite)
            else:
                data,offset=tridiagonalize_qr(H,q=qq.reshape([-1,model.nband]),m=nsite)
            chain=Chain(asarray(data[1]),concatenate([t0[newaxis],data[0]],axis=0))
        cl.append(chain)
    return cl

def check_spec(chains,rhofunc,wlist,mode='eval',smearing=1.):
    '''
    Check mapping quality for wilson chain.

    Parameters:
        :chains: list, a list of chain after mapping.
        :rhofunc: function, Hybridization function.
        :wlist: 1D array, the frequency space.
        :mode: str, Choose the checking method,
            * 'eval' -> check eigenvalues.
            * 'pauli' -> check pauli components, it is only valid for 2 band system.

        :smearing: The smearing factor.
    '''
    nband=chains[0].nband
    nz=len(chains)
    print('Recovering Spectrum ...')
    dlv=hybri_chain(tlist=array([chain.tlist for chain in chains]),elist=array([chain.elist for chain in chains]),wlist=wlist,smearing=smearing*1./nz)
    if nband==1:
        dlv=dlv[:,0,0]
        nplt=1
    elif nband==2 and mode=='pauli':

        dlv=array([s2vec(d) for d in dlv])