import numpy as np
from math import sqrt, exp

def tri2full(H_nn, UL='L'):
    """Fill in values of hermitian matrix.

    Fill values in lower or upper triangle of H_nn based on the opposite
    triangle, such that the resulting matrix is symmetric/hermitian.

    UL='U' will copy (conjugated) values from upper triangle into the
    lower triangle.

    UL='L' will copy (conjugated) values from lower triangle into the
    upper triangle.
    """
    N, tmp = H_nn.shape
    assert N == tmp, 'Matrix must be square'
    #assert np.isreal(H_nn.diagonal()).all(), 'Diagonal should be real'
    if UL != 'L':
        H_nn = H_nn.T

    for n in range(N - 1):
        H_nn[n, n + 1:] = H_nn[n + 1:, n].conj()

def dagger(matrix):
    return np.conj(matrix.T)

def rotate_matrix(h, u):
    return np.dot(u.T.conj(), np.dot(h, u))

def get_subspace(matrix, index):
    """Get the subspace spanned by the basis function listed in index"""
    assert matrix.ndim == 2 and matrix.shape[0] == matrix.shape[1]
    return matrix.take(index, 0).take(index, 1)   

permute_matrix = get_subspace

def normalize(matrix, S=None):
    """Normalize column vectors.

    ::

      <matrix[:,i]| S |matrix[:,i]> = 1

    """
    for col in matrix.T:
        if S is None:
            col /= np.linalg.norm(col)
        else:
            col /= np.sqrt(np.dot(col.conj(), np.dot(S, col)))

def subdiagonalize(h_ii, s_ii, index_j):
    nb = h_ii.shape[0]
    nb_sub = len(index_j)
    h_sub_jj = get_subspace(h_ii, index_j)
    s_sub_jj = get_subspace(s_ii, index_j)
    e_j, v_jj = np.linalg.eig(np.linalg.solve(s_sub_jj, h_sub_jj))
    normalize(v_jj, s_sub_jj) # normalize: <v_j|s|v_j> = 1
    permute_list = np.argsort(e_j.real)
    e_j = np.take(e_j, permute_list)
    v_jj = np.take(v_jj, permute_list, axis=1)
    
    #setup transformation matrix
    c_ii = np.identity(nb, complex)
    for i in xrange(nb_sub):
        for j in xrange(nb_sub):
            c_ii[index_j[i], index_j[j]] = v_jj[i, j]

    h1_ii = rotate_matrix(h_ii, c_ii)
    s1_ii = rotate_matrix(s_ii, c_ii)

    return h1_ii, s1_ii, c_ii, e_j

def cutcoupling(h, s, index_n):
    for i in index_n:
        s[:, i] = 0.0
        s[i, :] = 0.0
        s[i, i] = 1.0
        Ei = h[i, i]
        h[:, i] = 0.0
        h[i, :] = 0.0
        h[i, i] = Ei

def fermidistribution(energy, kt):
    #fermi level is fixed to zero
    return 1.0 / (1.0 + np.exp(energy / kt) )

def fliplr(a):
    length=len(a)
    b = [0] * length
    for i in range(length):
        b[i] = a[length - i - 1]
    return b

def plot_path(energy):
    import pylab
    pylab.plot(np.real(energy), np.imag(energy), 'b--o')
    pylab.show()
    

def function_integral(function, calcutype):
    #return the integral of the 'function' on 'intrange'    
    #the function can be a value or a matrix, arg1,arg2 are the possible
    #parameters of the function
    
    intctrl = function.intctrl
    if calcutype == 'eqInt':
        intrange = intctrl.eqintpath
        tol = intctrl.eqinttol
        if hasattr(function.intctrl, 'eqpath_radius'):
            radius = function.intctrl.eqpath_radius
        else:
            radius = -1
        if hasattr(function.intctrl, 'eqpath_origin'):
            origin = function.intctrl.eqpath_origin
        else:
            origin = 1000
    elif calcutype == 'neInt':
        intrange = intctrl.neintpath
        tol = intctrl.neinttol
        radius = -1
        origin = 1000
    elif calcutype == 'locInt':
        intrange = intctrl.locintpath
        tol = intctrl.locinttol
        if hasattr(function.intctrl, 'locpath_radius'):
            radius = function.intctrl.locpath_radius
        else:
            radius = -1
        if hasattr(function.intctrl, 'locpath_origin'):
            origin = function.intctrl.locpath_origin
        else:
            origin = 1000
    trace = 0    
    a = 0.
    b = 1.

    #Initialize with 13 function evaluations.
    c = (a + b) / 2
    h = (b - a) / 2
    realmin = 2e-17

    s = [.942882415695480, sqrt(2.0/3),
         .641853342345781, 1/sqrt(5.0), .236383199662150]
    s1 = [0] * len(s)
    s2 = [0] * len(s)
    for i in range(len(s)):
        s1[i] = c - s[i] * h
        s2[i] = c + fliplr(s)[i] * h
    x0 = [a] + s1 + [c] + s2 + [b]

    s0 = [.0158271919734802, .094273840218850, .155071987336585,
          .188821573960182,  .199773405226859, .224926465333340]
    w0 = s0 + [.242611071901408] + fliplr(s0)
    w1 = [1, 0, 0, 0, 5, 0, 0, 0, 5, 0, 0, 0, 1]
    w2 = [77, 0, 432, 0, 625, 0, 672, 0, 625, 0, 432, 0, 77]
    for i in range(len(w1)):
        w1[i] = w1[i] / 6.0
        w2[i] = w2[i] / 1470.0
                                                        
    dZ = [intrange[:len(intrange) - 1], intrange[1:]]
    hmin = [0] * len(dZ[1])

    path_type = []
    for i in range(len(intrange) - 1):
        rs = np.abs(dZ[0][i] - origin)
        re = np.abs(dZ[1][i] - origin)
        if abs(rs - radius) < 1.0e-8 and abs(re - radius) < 1.0e-8:
            path_type.append('half_circle')
        else:
            path_type.append('line')
   
    for i in range(len(dZ[1])):
        if path_type[i] == 'half_circle':
            dZ[0][i] = 0
            dZ[1][i] = np.pi
    for i in range(len(dZ[1])):
        dZ[1][i] = dZ[1][i] - dZ[0][i]
        hmin[i] = realmin / 1024 * abs(dZ[1][i])


    temp = np.array([[1] * 13, x0]).transpose()
    
    Zx = np.dot(temp, np.array(dZ))
      
    Zxx = []
    for i in range(len(intrange) - 1):
        for j in range(13):
            Zxx.append(Zx[j][i])

    ns = 0
    ne = 12
    if path_type[0] == 'line':
        yns = function.calgfunc(Zxx[ns], calcutype)
    elif path_type[0] == 'half_circle':
        energy = origin + radius * np.exp((np.pi - Zxx[ns + i]) * 1.j)
        yns = -1.j * radius * np.exp(-1.j* Zxx[ns +i])* function.calgfunc(energy, calcutype)        
    fcnt = 0
    

    for n in range(len(intrange)-1):
        # below evaluate the integral and adjust the tolerance
        Q1pQ0 = yns * (w1[0] - w0[0])
        Q2pQ0 = yns * (w2[0] - w0[0])
        fcnt = fcnt + 12
        for i in range(1,12):
            if path_type[n] == 'line':
                yne = function.calgfunc(Zxx[ns + i], calcutype)
            elif path_type[n] == 'half_circle':
                energy = origin + radius * np.exp((np.pi -Zxx[ns + i]) * 1.j)
                yne = -1.j * radius * np.exp(-1.j * Zxx[ns + i])* function.calgfunc(energy, calcutype)
            Q1pQ0 += yne * (w1[i] - w0[i])
            Q2pQ0 += yne * (w2[i] - w0[i])

        # Increase the tolerance if refinement appears to be effective
        r = np.abs(Q2pQ0) / (np.abs(Q1pQ0) + np.abs(realmin))
        dim = np.product(r.shape)
        r = np.sum(r) / dim
        if r > 0 and r < 1:
            thistol = tol / r
        else:
            thistol = tol
        if path_type[n] == 'line':
            yne = function.calgfunc(Zxx[ne], calcutype)
        elif path_type[n] == 'half_circle':
            energy = origin + radius * np.exp((np.pi -Zxx[ne]) * 1.j)
            yne = -1.j * radius * np.exp(-1.j * Zxx[ne])* function.calgfunc(energy, calcutype)
        #Call the recursive core integrator
       
        Qk, xpk, wpk, fcnt, warn = quadlstep(function, Zxx[ns],
                                            Zxx[ne], yns, yne,
                                            thistol, trace, fcnt,
                                            hmin[n], calcutype, path_type[n],
                                            origin, radius)
        if n == 0:
            Q = np.copy(Qk)
            Xp = xpk[:]
            Wp = wpk[:]
        else:
            Q += Qk
            Xp = Xp[:-1] + xpk
            Wp = Wp[:-1] + [Wp[-1] + wpk[0]] + wpk[1:]
        if warn == 1:
            print 'warning: Minimum step size reached,singularity possible'
        elif warn == 2:
            print 'warning: Maximum function count excced; singularity likely'
        elif warn == 3:
            print 'warning: Infinite or Not-a-Number function value encountered'
        else:
            pass
        
        ns += 13
        ne += 13
        yns = np.copy(yne)
      
    return Q,Xp,Wp,fcnt

def quadlstep(f, Za, Zb, fa, fb, tol, trace, fcnt, hmin, calcutype,
                                                   path_type, origin, radius):
    #Gaussian-Lobatto and Kronrod method
    #QUADLSTEP Recursive core routine for integral
    #input parameters:
    #      f      ----------   function, here we just use the module calgfunc
    #                          to return the value, if wanna use it for
    #                          another one, change it
    #     Za, Zb  ----------   the start and end point of the integral
    #     fa, fb  ----------   the function value on Za and Zb
    #     fcnt    ----------   the number of the funtion recalled till now
    #output parameters:
    #      Q      ----------   integral
    #     Xp      ----------   selected points
    #     Wp      ----------   weight
    #    fcnt     ----------   the number of the function recalled till now

    maxfcnt = 10000

    # Evaluate integrand five times in interior of subintrval [a,b]
    Zh = (Zb - Za) / 2.0
    if abs(Zh) < hmin:
        # Minimun step size reached; singularity possible
        Q = Zh * (fa + fb)
        if path_type == 'line':
            Xp = [Za, Zb]
        elif path_type == 'half_circle':
            Xp = [origin + radius * np.exp((np.pi - Za) * 1.j),
                                      origin + radius * np.exp((np.pi - Zb) * 1.j)]
        Wp = [Zh, Zh]
        warn = 1
        return Q, Xp, Wp, fcnt, warn
    fcnt += 5
    if fcnt > maxfcnt:
        #Maximum function count exceed; singularity likely
        Q = Zh * (fa + fb)
        if path_type == 'line':
            Xp = [Za, Zb]
        elif path_type == 'half_circle':
            Xp = [origin + radius * np.exp((np.pi - Za) * 1.j),
                                      origin + radius * np.exp((np.pi - Zb) * 1.j)]
        Wp = [Zh, Zh]
        warn = 2
        return Q, Xp, Wp, fcnt, warn
    x = [0.18350341907227,   0.55278640450004,   1.0,
         1.44721359549996,   1.81649658092773];
    Zx = [0] * len(x)
    y = [0] * len(x)
    for i in range(len(x)):
        x[i] *= 0.5
        Zx[i] = Za + (Zb - Za) * x[i]
        if path_type == 'line':
            y[i] = f.calgfunc(Zx[i], calcutype)
        elif path_type == 'half_circle':
            energy = origin + radius * np.exp((np.pi - Zx[i]) * 1.j)
            y[i] = f.calgfunc(energy, calcutype)
    #Four point Lobatto quadrature
    s1 = [1.0, 0.0, 5.0, 0.0, 5.0, 0.0, 1.0]
    s2 = [77.0, 432.0, 625.0, 672.0, 625.0, 432.0, 77.0]
    Wk = [0] * 7
    Wp = [0] * 7
    for i in range(7):
        Wk[i] = (Zh / 6.0) * s1[i]
        Wp[i] = (Zh / 1470.0) * s2[i]
    if path_type == 'line':
        Xp = [Za] + Zx + [Zb]
    elif path_type == 'half_circle':
        Xp = [Za] + Zx + [Zb] 
        for i in range(7):
            factor = -1.j * radius * np.exp(1.j * (np.pi - Xp[i]))
            Wk[i] *= factor
            Wp[i] *= factor
            Xp[i] = origin + radius * np.exp((np.pi - Xp[i]) * 1.j)
    Qk = fa * Wk[0] + fb * Wk[6]
    Q = fa * Wp[0] + fb * Wp[6]
    for i in range(1, 6):
        Qk += y[i-1] * Wk[i]
        Q  += y[i-1] * Wp[i]
    if np.isinf(np.max(np.abs(Q))):
        Q = Zh * (fa + fb)
        if path_type == 'line':
            Xp = [Za, Zb]
        elif path_type == 'half_circle':
            Xp = [origin + radius * np.exp((np.pi - Za) * 1.j),
                                      origin + radius * np.exp((np.pi - Zb) * 1.j)]
        Wp = [Zh, Zh]
        warn = 3
        return Qk, Xp, Wp, fcnt, warn
    else:
        pass
    if trace:
        print fcnt, real(Za), imag(Za), abs(Zh)
    #Check accurancy of integral over this subinterval
    XXk = [Xp[0], Xp[2], Xp[4], Xp[6]]
    WWk = [Wk[0], Wk[2], Wk[4], Wk[6]]
    YYk = [fa, y[1], y[3], fb]
    if np.max(np.abs(Qk - Q)) <= tol:
        warn = 0
        return Q, XXk, WWk, fcnt, warn
    #Subdivide into six subintevals
    else:
        Q, Xk, Wk, fcnt, warn = quadlstep(f, Za, Zx[1], fa, YYk[1],
                                           tol, trace, fcnt, hmin,
                                               calcutype, path_type,
                                                origin, radius)

        Qk, xkk, wkk, fcnt, warnk = quadlstep(f, Zx[1],
                          Zx[3], YYk[1], YYk[2], tol, trace, fcnt, hmin,
                                             calcutype, path_type,
                                             origin, radius)
        Q += Qk
        Xk = Xk[:-1] + xkk
        Wk = Wk[:-1] + [Wk[-1] + wkk[0]] + wkk[1:]
        warn = max(warn, warnk)
        
        Qk, xkk, wkk, fcnt, warnk = quadlstep(f, Zx[3], Zb, YYk[2], fb,
                                           tol, trace, fcnt, hmin,
                                                   calcutype, path_type,
                                                   origin, radius)
        Q += Qk
        Xk = Xk[:-1] + xkk
        Wk = Wk[:-1] + [Wk[-1] + wkk[0]] + wkk[1:]
        warn = max(warn, warnk)
    return Q, Xk, Wk, fcnt, warn

def mytextread0(filename):
    num = 0
    df = file(filename)
    df.seek(0)
    for line in df:
        if num == 0:
            dim = line.strip().split(' ')
            row = int(dim[0])
            col = int(dim[1])
            mat = np.empty([row, col])
        else:
            data = line.strip().split(' ')
            if len(data) == 0 or len(data)== 1:
                break
            else:
                for i in range(len(data)):
                    mat[num - 1, i] = float(data[i])
        num += 1
    return mat

def mytextread1(filename):
    num = 0
    df = file(filename)
    df.seek(0)
    data = []
    for line in df:
        tmp = line.strip()
        if len(tmp) != 0: 
            data.append(float(tmp))
        else:
            break
    dim = int(sqrt(len(data)))
    mat = np.empty([dim, dim])
    for i in range(dim):
        for j in range(dim):
            mat[i, j] = data[num]
            num += 1
    return mat

def mytextwrite1(filename, mat):
    num = 0
    df = open(filename,'w')
    df.seek(0)
    dim = mat.shape[0]
    if dim != mat.shape[1]:
        print 'matwirte, matrix is not square'
    for i in range(dim):
        for j in range(dim):
            df.write('%20.20e\n'% mat[i, j])
    df.close()