# A simulation algorithm for the Ising model
# A Gibbs sampler first--I've fudged a little since I use random
# selection of pixels rather than moving through them systematically
ising.gs=function(nx = 50, ny=20, beta = .2, istop = 100)
{
	x <- matrix(0,nrow=nx+2,ncol=ny+2)
        x[2:(nx+1),2:(ny+1)]=matrix(sample(c(-1,1),nx*ny,replace=T),nrow=nx,ncol=ny)
        x.init=x[2:(nx+1),2:(ny+1)]
        iestop=istop*nx*ny
	for(i1 in 1:iestop) {
           i=sample(1:nx,1)+1
           j=sample(1:ny,1)+1
           r=exp(2*beta*(x[i,j-1]+x[i,j+1]+x[i-1,j]+x[i+1,j]))
           r=r/(1+r)
           if(runif(1)<r) x[i,j]=1 else x[i,j]=-1
	}
	gridout <- list(x[2:(nx+1),2:(ny+1)],x.init)
        names(gridout)=c("Final Grid","Start Grid")
	return(gridout)
}
ising.mh=function(nx = 50, ny=20, beta = .2, istop = 100)
{
	x <- matrix(0,nrow=nx+2,ncol=ny+2)
        x[2:(nx+1),2:(ny+1)]=matrix(sample(c(-1,1),nx*ny,replace=T),nrow=nx,ncol=ny)
        x.init=x[2:(nx+1),2:(ny+1)]
        iestop=istop*nx*ny
	for(i1 in 1:iestop) {
           i=sample(1:nx,1)+1
           j=sample(1:ny,1)+1
# The only other choice for a move is -x, and the proposal distribution
# is simply the random selection of a new pixel, which is symmetric by
# default, so the jump probability simplifies to the form below:
           r=exp(-2*beta*x[i,j]*(x[i,j-1]+x[i,j+1]+x[i-1,j]+x[i+1,j]))
           if(runif(1)<r) x[i,j]=-x[i,j]
	}
	gridout <- list(x[2:(nx+1),2:(ny+1)],x.init)
        names(gridout)=c("Final Grid","Start Grid")
	return(gridout)
}
# Here are commands for the MH model output.  Be sure to load the 
# library lattice before proceeding
# The contour command doesn't work, but it sure looks interesting
y=ising.mh()
contour(1:50,1:50,y$'Start Grid')
levelplot(y$'Start Grid',col.regions=c("black","white"),cuts=1)
win.graph()
levelplot(y$'Final Grid',col.regions=c("black","white"),cuts=1)
y=ising.mh(beta=1.5)
levelplot(y$'Start Grid',col.regions=c("black","white"),cuts=1)
levelplot(y$'Final Grid',col.regions=c("black","white"),cuts=1)
y=ising.mh(beta=-2)
levelplot(y$'Start Grid',col.regions=c("black","white"),cuts=1)
levelplot(y$'Final Grid',col.regions=c("black","white"),cuts=1)