# Illustrating the idea of the MC on tables with fixed row and col sums

x0 = matrix(c(8,1,4,2,7,1,3,1,1),nrow = 3)

#> x0
#     [,1] [,2] [,3]
#[1,]    8    2    3
#[2,]    1    7    1
#[3,]    4    1    1
                                      
# x = x0
move = function(x, print=T) {
	m = dim(x)[1]
	n = dim(x)[2]
	temp = sample(m,2)
	i1 = min(temp)
	i2 = max(temp)
	temp = sample(n,2)
	j1 = min(temp)
	j2 = max(temp)
	delta1 = matrix(c(1,-1,-1,1),nrow=2)
	delta2 = matrix(c(-1,1,1,-1),nrow=2)
	change = F
	if(runif(1) < 0.5){
    if(min(x[i1,j2],x[i2,j1]) > 0){
        change = T
        x[c(i1,i2),c(j1,j2)]=x[c(i1,i2),c(j1,j2)]+delta1
    }
	}	else{
    if(min(x[i1,j1],x[i2,j2]) > 0){
        change = T
        x[c(i1,i2),c(j1,j2)]=x[c(i1,i2),c(j1,j2)]+delta2
    }
	}
	if(print){
    cat(paste("(i1,i2)=(",i1,",",i2,") (j1,j2)=(",j1,",",j2,")",sep=""))
    if(!change)cat("  no change")
    cat("\n")
    print(x)
	}
	return(x)
}


x = x0
x0
# repeat this several times to see how it works:
x = move(x)   

nit = 10                         # <-- change to, e.g., 10000
results = array(0, dim=c(3,3,nit))  # <-- just to set up placeholders
x = x0                          # or could try x = 100*x0
results[,,1] = x
# results
for(t in 2:nit){   
  x = move(x,print=F)
  results[,,t] = x
}

corner = results[1,1,]  
# ^^ can change this to collect different elements
                        #  # Or, equivalently:
                        #  f = function(x){x[1,1]}  
                        #  corner = apply(results, 3, f)