#-----------------------------#
#                             #
#    SDEs from Theorem 4.2    #
#      w1 = 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: 17 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          #
#---------------------#

SolveSDEperCoreEW = function(npath,nt,LM,dt,epsh,u)
{
  
  sqrtdt = sqrt(dt)
  absu = abs(u)
  epsu = sign(u)
  Id   = diag(rep(1,LM))
  Null  = diag(rep(0,LM))
  
  # Notation: rhouu -> G,  Fuu -> F
  
  G0 = Null                                    
  F0 = Null
  
  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 )
    {
      phi = rnorm(LM)
      dx = diag( sqrtdt*phi )
      
      # Calculate DeltaG and DeltaF #

      DG = diag(diag(G))
      dh = -epsh*dt + absu*dt * DG/2 - sqrt(absu) * dx  
      dhG = dh %*% G
      Mat = (Id-F) %*% dhG 
      DeltaG = ( (Id-F) %*% dh %*% (Id+F) - G %*% dhG )/2 
      DeltaF = ( Mat - t(Mat) )/2

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

      weps = ( sum( diag( epsh%*%Gnew ) ) )
      wu   = -absu/4 * sum( diag(Gnew)^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 SolveSDEperCoreEW  #


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% SolveSDEperCoreEW(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)  #

E1Weps  = rep(0,nt)
E1Wu    = rep(0,nt)
E1W     = rep(0,nt)
minint1W = 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)
{
  minint1W[k-1] = min(intW[k-1,iValid])
  Zk = exp( -( intW[k-1,iValid] - minint1W[k-1] ) )
  
  E1Weps[k] = sum( Weps[k,iValid]/LM * Zk ) / sum(Zk)
  E1Wu[k] = sum( Wu[k,iValid]/LM * Zk ) / sum(Zk)
}
E1W = E1Weps + E1Wu


# picture: 

info = "Exact Diagonalization (green) vs. Monte Carlo (blue)"
plot(tk, E1W, 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)