# a simple implementation of 2D  Ising model
# 
#   version 2 deals with 8-pt neighborhood (for simplicity) and also has 1-color boundary
# 

 pix <- function(i,j,col){ 
    if (col == 1) clr = "white"
    if (col == -1) clr = "red"
    rect(dx*(i-1), dx*(j-1), dx*i, dx*j , col= clr)
}

 N <- 32;
 par(mfrow=c(1,2))
 plot(c(0,1),c(0,1), col="white") 
 dx = 1/N

 # start with all white 

 X <- sample(c(-1,1),N*N, replace=TRUE)
 dim(X) <- c(N,N)

 X[1,] <- 1;   X[N,] <- 1; X[,1] <- 1;  X[,N] <- 1

 # X <- matrix(1,N,N)


 for (i in 1:N)
   for (j in 1:N){
   pix(i,j,X[i,j])
 }

 
 T <- 100000        	# simulation length
 Et <- rep(0,T/10)	# this is used to estimate the partition "likelihood" 
 E <- 0    # 2*N*(N-1)          # initially, E = 2*N*(N-1) for all same color

 beta <- 0.3

 for (t in 1:T){
   ij <- floor(runif(2)*(N-2)) + 2
   i <- ij[1]
   j <- ij[2]
   
   sumall <- sum(X[(i-1):(i+1), (j-1):(j+1)] ) - X[i,j]

   Eold <- X[i,j] * sumall
   flip <- -X[i,j] 
   Enew <- -Eold
   dE <- Enew - Eold

    ## the rule for Metropolis sampler: if energy increases, do it. 
    #  If energy decreases, do it with probability P = e(-beta*delta E) 

    p <- exp(beta*min(0, dE))
    rnd <- runif(1)
    if (rnd < p){
      X[i,j] <- flip
      E <- E + dE       # total energy
      pix(i,j,X[i,j])
    }
    if((t %% 10) == 0) Et[t %/% 10] <- E

 }

  plot(-Et)