#!/usr/bin/env python
import os, time
from scipy import *
from pylab import *

def noninteracting_G0(beta, N):
    om = arange(1,2*N,2)*pi/beta
    return (om, 2j*om*(1-sqrt(1+1/om**2)) )  # Bethe lattice, half-filling

def InverseFourier(Giom, om, tau):
    "Inverse Fourier computes G(tau) from G(iom)"
    beta = pi/om[0]
    Gtau = zeros(len(tau), dtype=float)
    for it,t in enumerate(tau):
        dsum=0
        for iw,w in enumerate(om):
            dsum += cos(w*t)*Giom[iw].real + sin(w*t)*(Giom[iw].imag+1./w)
        Gtau[it] = 2*dsum/beta - 0.5
    return Gtau

def g0_2D(tau, G0):
    L = len(tau)-1
    print 'L=', L
    g0 = zeros((L,L), dtype=float)
    for i in range(L):
        for j in range(L):
            if i>=j:
                g0[i,j] = -G0[i-j]
            else:
                g0[i,j] = G0[L+i-j]
    return g0

def CleanUpdate(g, g0, vn):
    L = len(vn)
    A = zeros((L,L), dtype=float)
    for i,s in enumerate([1,-1]):
        # a = e^V-1
        a = [exp(vn[l]*s)-1 for l in range(L)]
        print 'a=', a
        # Matrix A = 1 + (1-g)(e^V-1) == (1+a)*I - g*a
        for l1 in range(L):
            for l2 in range(L): A[l1,l2] = -g0[l1,l2]*a[l2]
            A[l1,l1] += 1+a[l1]
        print 'A=', A

        g[i] = solve(A,g0)
    

if __name__=='__main__':
    beta = 16.       # inverse temperature
    N = 1500         # number of matsubara points
    L = 64           # number of time slices
    U = 2.           # interaction U
    #######
    nwarm0 = 100     # how many sweeps to warm up
    nmeas0 = 10      # how often to meassure (in sweeps)
    nsteps = 100000  # number of all MC steps
    ncout  = 10000   # how often to print info

    nwarmup = int(nwarm0*L)
    measure = int(nmeas0*L)

    (om, G0start) = noninteracting_G0(beta, N)

    seed(1)
    S_ising = sign(rand(L)-0.5)

    G0iom = G0start
    tau = linspace(0,beta,L+1)
    G0tau = InverseFourier(G0iom, om, tau)

    xlam = arccosh(exp(0.5*(tau[1]-tau[0])*U)) # lambda
    vn = array(S_ising)*xlam

    # g0 in 2D representation
    g0 = g0_2D(tau, G0tau)
    g = zeros((2,L,L), dtype=float)
    CleanUpdate(g, g0, vn)
    print 'g='
    print g
    
    #plot(tau[:-1], g0[:,0], label='1')
    #plot(tau[:-1], g0[0,:], label='2')
    #legend(loc='best')
    #show()