#-----------------------------#
#                             #
#    SDEs from Theorem 4.2    #
#      w2 = 1, EW only        #
#                             #
#-----------------------------#


require(doParallel)
require(doRNG)

ncores = detectCores()                   
ncores                                        # 20 on my pc
ncores = 16
mycl = makeCluster(ncores)
registerDoParallel(mycl)

# if you are using Microsoft R Open with Intel Math Kernel Library:
# set the number of threads to 1 for each cluster process ( => no 
# automatic parallelization of matrix products )

clusterEvalQ( cl = mycl , {setMKLthreads(1)} )


 
set.seed(123456)



###   Start Input Variables   ###

npath = 10000                                 # for 10K: 38 seconds on my pc

L = 3
M = 2                                         # L x M Gitter

u = 4
vep = -1.0                                    # varepsilon, hopping parameter

T = 4    
dt = 0.01        

bc = "none"                                   # boundary conditions
#bc = "pbc"

###   End Input Variables   ###



tk = seq(0,T,by=dt)
nt = length(tk)                       

LM   = L*M

site = matrix(0,nrow=LM,ncol=2)
k = 1
for(jy in 1:M)
{
  for(jx in 1:L)
  {
    site[k,] = c(jx,jy)                     
    k = k+1
  }
}

eps = matrix(0,LM,LM) 
eps0 = matrix(0,LM,LM)                                    

for(k in 1:LM)
{
  for(ell in 1:LM)
  {
    abstand = sum( (site[k,]-site[ell,])^2 )
    if( abstand==1  )
    {
      eps[k,ell] = vep 
      eps0[k,ell] = vep
    }
    
    # periodic bc's
    d = site[k,]-site[ell,]
    dx = d[1]
    dy = d[2]
    if( abs(dx) == L-1 & dy==0 )
    {
      eps[k,ell] = vep 
    }
    if( dx==0 & abs(dy) == M-1 )
    {
      eps[k,ell] = vep 
    }
    # end periodic bc's
    
  }
}

if(bc=="none")
{
  epsh = eps0
}
if(bc=="pbc")
{
  epsh = eps
}



#---------------------#
#        SDE          #
#---------------------#

SolveSDEperCoreE2W = function(npath,nt,LM,dt,epsh,u)
{
  
  LMS = 2*LM
  sqrtdt = sqrt(dt)
  absu = abs(u)
  epsu = sign(u)
  Id   = diag(rep(1,LM))
  Null  = diag(rep(0,LM))
  Id2   = diag(rep(1,LMS))
  Null2  = diag(rep(0,LMS))
  epsepsh = rbind( cbind( epsh , Null ) , cbind( Null , epsh ) )
  
  # Notation: rho -> G,  F -> F
  
  G0 = Null2                                    
  F0 = Null2
  
  Weps = matrix(0,nrow=nt,ncol=npath)
  Wu   = matrix(0,nrow=nt,ncol=npath)
  intW = matrix(0,nrow=nt,ncol=npath)
  
  nNA = rep(0,nt)
  iValid = rep(TRUE,npath)
  
  for( i in 1:npath)
  {
    G = G0
    F = F0
    
    # solve SDE
    for( k in 2:nt )
    {
      xi = rnorm(LM)
      dy = diag( sqrtdt*xi )
      dBy = rbind( cbind(Null,dy) , cbind(epsu*dy,Null) ) 
      
      # Calculate DeltaG and DeltaF #

      Guu = G[1:LM,1:LM]
      Gud = G[1:LM,(LM+1):LMS]
      Fuu = F[1:LM,1:LM]
      Fud = F[1:LM,(LM+1):LMS]
      Gdu = epsu*Gud
      Fdu = epsu*Fud
      Gdd = Guu
      Fdd = Fuu
      tiF = rbind( cbind( Fdd , Fdu ) , cbind( Fud , Fuu ) )
      DGud = diag(diag(Gud))
      DGdu = diag(diag(Gdu))
      temp1 = cbind( Null , DGud )
      temp2 = cbind( epsu*DGud , Null )
      DG = rbind( temp1 , temp2 )
      dh = -epsepsh*dt + absu*dt * DG/2 - sqrt(absu) * dBy  
      dht = t(dh)
      Gdht = G %*% dht
      Mat = Gdht %*% (Id2+F) 
      DeltaG = ( (Id2-F) %*% dh %*% (Id2+tiF) - Gdht %*% G )/2 
      DeltaF = ( t(Mat) - Mat )/2

      # END Calculate DeltaG and DeltaF #
      
      Gnew = G + DeltaG 
      Fnew = F + DeltaF
      
      # Calculate energies and check path #

      Guu = Gnew[1:LM,1:LM]
      Gud = Gnew[1:LM,(LM+1):LMS]
      Fuu = Fnew[1:LM,1:LM]
      Fud = Fnew[1:LM,(LM+1):LMS]
      
      weps = ( sum( diag( epsh%*%Guu ) ) )
      wu   = -absu/4 * sum( diag(Gud)^2 )
      nrgdensity = ( abs(weps) + abs(wu) )/LM
      
      # KEEP PATH ?
      if( is.na(nrgdensity) || nrgdensity > 10 )
      {
        nNA[k] = nNA[k] + 1
        iValid[i] = FALSE
        break              # ignore path, leave k-loop, start new MC-path
      }
      # END Calculate energies and check path # 
      
      # Update G and F and store energies #
      G = Gnew
      F = Fnew
      Weps[k,i] = weps
      Wu[k,i]   = wu
      
    } # next k
    
  } # next path
  
  result = list(Weps,Wu,iValid,sum(nNA))
  
}  #  END SolveSDEperCoreE2W  #


AddResEW = function(list1,list2)
{
  res1 = cbind(list1[[1]],list2[[1]])
  res2 = cbind(list1[[2]],list2[[2]])
  res3 = c(list1[[3]],list2[[3]])
  res4 = list1[[4]] + list2[[4]] 
  res = list(res1,res2,res3,res4)
  return(res)
}


npercore = as.integer(npath/ncores)
npath = npercore*ncores

calctime = Sys.time()
res = foreach(i=1:ncores, .combine='AddResEW') %dorng% SolveSDEperCoreE2W(npercore,nt,LM,dt,epsh,u)
Sys.time() - calctime



Weps = res[[1]]
Wu   = res[[2]]
iValid = res[[3]] 
sumNA = res[[4]]
sumNA

##  Notation:  E1W := EW(w1=1), E2W := EW(w2=1)  ##

E2Weps  = rep(0,nt)
E2Wu    = rep(0,nt)
E2W     = rep(0,nt)
minint2W = rep(0,nt)
intW = matrix(0,nrow=nt,ncol=npath)
for( i in 1:npath)
{
  intW[,i] = cumsum( (Weps[,i] + Wu[,i]) * dt )
}


##  Plain MC: Energies EWeps, EWu, EW  ##

for(k in 2:nt)
{
  minint2W[k-1] = min(intW[k-1,iValid])
  Zk = exp( -( intW[k-1,iValid] - minint2W[k-1] ) )
  
  E2Weps[k] = sum( Weps[k,iValid]/LM * Zk ) / sum(Zk)
  E2Wu[k] = sum( Wu[k,iValid]/LM * Zk ) / sum(Zk)
}
E2W = E2Weps + E2Wu


# picture: 

info = "Exact Diagonalization (green) vs. Monte Carlo (blue)"
plot(tk, E2W, col="blue", main=info, font.main=10, xlab="beta", ylab="< H >_beta" )

# the quantity EH is to be calculated with 2D-Fermi-Hubbard-Exact-Diagonalization.txt, then:
points(tk,EH,col="green",cex=0.5)


# atomistic limit:
points(tk, -u/4 *tanh(tk*u/4), cex=0.5, col="cyan" ) 





stopCluster(mycl)