nn <- function(m, c) {
    nr <- nrow(m)
    nc <- ncol(m)
    nn <- matrix(0, nr, nc)
    nn[1:(nr)-1,] <- nn[1:(nr)-1,] + (m[2:nr,] == c)
    nn[2:nr,] <- nn[2:nr,] + (m[1:(nr-1),] == c)
    nn[,1:(nc)-1] <- nn[,1:(nc)-1] + (m[,2:nc] == c)
    nn[,2:nc] <- nn[,2:nc] + (m[,1:(nc-1)] == c)
    nn
}

simGroup <- function(m, l2, l1, beta, which) {
    pp2 <- l2 * exp(beta * nn(m, 2))
    pp1 <- l1 * exp(beta * nn(m, 1))
    pp <- pp2 / (pp2 + pp1)
    ifelse(runif(sum(which)) < pp[which], 2, 1)
}


simImgV <- function(m, img, N, beta, p) {
    white <- outer(1:nrow(m), 1:ncol(m), FUN=`+`) %% 2 == 1
    black <- ! white
    if (is.null(img)) {
        l2 <- 1
        l1 <- 1
    }
    else {
        l2 <- ifelse(img == 2, p, 1 - p)
        l1 <- ifelse(img == 1, p, 1 - p)
    }
    for (i in 1:N) {
	m[white] <- simGroup(m, l2, l1, beta, white)
	m[black] <- simGroup(m, l2, l1, beta, black)
    }
    m
}

# load the data
load("mrf.Rda")

# Draws from the sampler taken 50 apart
image(img)
m <- img
for (i in 1:40) {
    m <- simImgV(m, img ,50, 0.9 , .7)
    image(m)
}

# Posterior means based on draws taken 50 apart
image(img)
mm<-img
m<-img
N<-0
for (i in 1:40) {
    m <- simImgV(m, img, 50, 0.9, .7)
    mm <- (N / (N + 1)) * mm + (1 / (N + 1)) * m; N <- N + 1
    image(mm)
}