## Gibbs sampler for the Ising model

f.gibbs <- function(x,beta)
{
  ## Purpose: Gibbs sampler for the Ising model. It alternates between
  ##          updating values at positions (i,j) with i+j even and
  ##          updating values at positions where i + j odd.
  ## -------------------------------------------------------------------------
  ## Arguments: x=configuration (n times n matrix with n even and values +1/-1
  ##            beta=inverse temperature
  ## -------------------------------------------------------------------------
  ## Author: Hans-Ruedi Kuensch, Date: 22 Nov 2001, 12:10
  beta <- 2*beta
  n <- nrow(x)
  nq <- 0.25*n*n # Number of indices divided by 4.
  ger <- 2*(1:(n/2)) # even indices
  ung <- ger-1  # odd indices 
  z <- rep(0,n) # set x to zero outside the {1,2,...,n}^2
  nb <- rbind(x[-1,],z) + rbind(z,x[-n,]) + cbind(x[,-1],z) + cbind(z,x[,-n])
                                        # sum of neighbors
  nb1 <- nb[ung,ung] # i,j both odd
  pr <- exp(beta*nb1) # Compute conditional probabilities given the neighbors
  pr <- pr/(pr+1)
  x[ung,ung] <- 2*(pr > runif(nq))-1 # update of x
  nb1 <- nb[ger,ger] # update for i,j both even, same thing
  pr <- exp(beta*nb1)
  pr <- pr/(pr+1)
  x[ger,ger] <- 2*(pr > runif(nq))-1
  nb <- rbind(x[-1,],z) + rbind(z,x[-n,]) + cbind(x[,-1],z) + cbind(z,x[,-n])
                  # recomputing the sum of neighbors
  nb1 <- nb[ung,ger] # update for i odd, j even
  pr <- exp(beta*nb1)
  pr <- pr/(pr+1)
  x[ung,ger] <- 2*(pr > runif(nq))-1
  nb1 <- nb[ger,ung] # update for i even, j odd
  pr <- exp(beta*nb1)
  pr <- pr/(pr+1)
  x[ger,ung] <- 2*(pr > runif(nq))-1
  x
}

par(mfrow=c(2,2))
n <- 32
beta <- 0.25
# phase transition occurs for beta>0.5*log(1+sqrt(2))=0.4407
x <- matrix(2*(runif(n*n) < 0.5)-1,nrow=n,ncol=n) # i.i.d. starting values
image(1:n,1:n,z=x,col=gray(0:1),xaxs="i",yaxs="i",axes=F,xlab="",ylab="")
for (i in (1:1000)) x <- f.gibbs(x,beta)
image(1:n,1:n,z=x,col=gray(0:1),xaxs="i",yaxs="i",axes=F,xlab="",ylab="")
      # visualization after a fixed number of iterations 

# behavior of mean magnetisation for a large number of iterations
m <- rep(0,10000)

for (i in (1:10000)){
  x <- f.gibbs(x,beta)
  m[i] <- mean(x)
}
par(mfrow=c(2,1))
plot(m[10*(1:1000)],type="l",xlab="Iteration",ylab="Magnetisation")
hist(m,breaks=30) #distribution of magnetisation 
acf(m,lag.max=1000) # autocorrelations, shows the degree of dependence
                    # for values separated by 1,2,... iterations
summary(m)
sqrt(var(m))

#Generating a nice figure (also for other courses)
n <- 128 #be patient
pdf("Ising.pdf")
par(mfrow=c(2,2),pty="s")
x <- matrix(2*(runif(n*n) < 0.5)-1,nrow=n,ncol=n) #Startwert i.i.d.
image(1:n,1:n,z=x,col=gray(0:1),xaxs="i",yaxs="i",pty="s",
      axes=F,xlab="",ylab="")
beta <- 0.25
for (i in (1:1000)) x <- f.gibbs(x,beta)
image(1:n,1:n,z=x,col=gray(0:1),xaxs="i",yaxs="i",pty="s",
      axes=F,xlab="",ylab="")
beta <- 0.35
for (i in (1:1000)) x <- f.gibbs(x,beta)
image(1:n,1:n,z=x,col=gray(0:1),xaxs="i",yaxs="i",pty="s",
      axes=F,xlab="",ylab="")
beta <- 0.44
for (i in (1:1000)) x <- f.gibbs(x,beta)
image(1:n,1:n,z=x,col=gray(0:1),xaxs="i",yaxs="i",pty="s",
      axes=F,xlab="",ylab="")
dev.off()
par(pty="m")