beta.param<-0.95

set.seed(100)
require(PottsUtils)
neighbors <- getNeighbors(mask=matrix(1, 128, 128), neiStruc=c(2,2,0,0))
blocks <- getBlocks(mask=matrix(1, 128, 128), nblock=2)
GB.sample<-BlocksGibbs(n=500, nvertex=128^2, ncolor=2, neighbors=neighbors, blocks=blocks, beta=beta.param)
pdf('Ising-BlockGibbs.pdf',paper='USr',width=9,height=9)
par(pty='s')
for(j in 1:500){
	image(0:128,0:128,matrix(GB.sample[,j],128,128),col=c('red','black'),xlab='x',ylab='y')
	title(substitute(beta==be,list(be=beta.param)))
}
dev.off()

set.seed(100)
edges <- getEdges(mask=matrix(1, 128, 128), neiStruc=c(2,2,0,0))
SW.sample<-SW(n=500, nvertex=128^2, ncolor=2, edges, beta=beta.param)

pdf('Ising-SW.pdf',paper='USr',width=9,height=9)
par(pty='s')
for(j in 1:500){
	image(0:128,0:128,matrix(SW.sample[,j],128,128),col=c('red','black'),xlab='x',ylab='y')
	title(substitute(beta==be,list(be=beta.param)))
}
dev.off()

count.match<-function(Xm){

	Nv<-nrow(Xm)
	count<-0
	for(i in 1:Nv){
		for(j in 1:Nv){
			imin<-i-1
			iplus<-i+1
			jmin<-j-1
			jplus<-j+1
			if(i == 1){imin<-Nv}
			if(i == Nv){iplus<-1}
			if(j == 1){jmin<-Nv}
			if(j == Nv){jplus<-1}
			count<-count+(Xm[i,j]==Xm[i,jmin])
			count<-count+(Xm[i,j]==Xm[i,jplus])
			count<-count+(Xm[i,j]==Xm[imin,j])
			count<-count+(Xm[i,j]==Xm[iplus,j])
		}
	}
	return(count)
	
}

SW.counts<-GB.counts<-rep(0,500)
for(j in 1:500){
	SW.counts[j]<-count.match(matrix(SW.sample[,j],128,128))
	GB.counts[j]<-count.match(matrix(GB.sample[,j],128,128))
}
plot(1:500,GB.counts,pch=19,ylim=range(55000,60000),xlab='Iteration',ylab='Score')
points(1:500,SW.counts,pch=3,col='blue')