## com100e of october 2021
#  revisiting magical squares and so on  
#  trying first 4 x 4 magical squares,
#  summing to mean(1:16)*4 = 34 in rows and columns
#  perhaps also diagonals
#  the programme also works for other sizes,
#  but will be slower for bigger squares 

#  see the two Nils FocuStat Blog Posts, 
#  on sudoko solving and magical squares 

kk <- 5
magic <- 0.5*kk*(kk^2+1) # 34 for magical 4 x 4 squares

rr <- 0*(1:kk)
cc <- 0*(1:kk) 

R <- function(x)
{
for (i in 1:kk)
{
rr[i] <- abs( sum(x[i, ]) - magic )
}
sum(rr)
}

C <- function(x)
{
for (j in 1:kk)
{
cc[j] <- abs( sum(x[ ,j]) - magic )
}
sum(cc)
}

D <- function(x)
{
oi1 <- sum(diag(x))
oi2 <- sum(diag(t(x))) 
abs(oi1 - magic) + abs(oi2 - magic)
}

gaudiset <- 1:(kk^2)
aux <- sample(gaudiset,kk^2)
xx <- array(aux,dim=c(kk,kk)) # i start my mcmc here 

Q <- function(x)
{
R(x) + C(x) + D(x) 
}

sim <- 10*10^4 # make it bigger, for bigger squares 
uu <- runif(sim)
ok <- 0*(1:sim)
probs <- 0*(1:sim) 
Qcheck <- 0*(1:sim) 
Qcheck[1] <- Q(xx) 

xxsim <- array(0,dim=c(sim,kk,kk))
xxsim[1, , ] <- xx

rAsim <- sample(1:kk,sim,replace=T)
cAsim <- sample(1:kk,sim,replace=T)
rBsim <- sample(1:kk,sim,replace=T)
cBsim <- sample(1:kk,sim,replace=T)

lambda0 <- 1.234
# can play with smaller and bigger values
# and a good values depends on the k x k size 

for (ss in 2:sim)
{
old <- xxsim[ss-1, , ]
# 
prop <- old
prop[rAsim[ss],cAsim[ss]] <- old[rBsim[ss],cBsim[ss]]
prop[rBsim[ss],cBsim[ss]] <- old[rAsim[ss],cAsim[ss]]
# 
diff <- Q(prop) - Q(old)
# pure metropolis: fixed lambda0,
# but i can make it slowly grow with time: 
lambda <- lambda0 + 0.00*log(ss) 
probs[ss] <- min(1,exp(-lambda*diff))
# 
ok[ss] <- 1*(uu[ss] <= probs[ss])
xxsim[ss, , ] <- (1-ok[ss])*old + ok[ss]*prop
Qcheck[ss] <- Q(xxsim[ss, , ])
print(c(ss,Qcheck[ss]))
}

print(range(Qcheck))

starthere <- 1001 # after burn-in

par(bty="l")
par(cex.lab=1.4142) 

matplot(starthere:sim,Qcheck[starthere:sim],
   type="l",ylim=c(0,max(Qcheck[starthere:sim])),
   xlab="long magical chain is running",ylab="Q score",lwd=2)

ssgood <- (1:sim)[Qcheck==min(Qcheck)]
firsttime <- min(ssgood)
print(firsttime)

magicsolution <- xxsim[firsttime, , ]

print("here it is")
print(magicsolution) 
print(c(R(magicsolution),C(magicsolution),D(magicsolution)))

for (i in 1:kk)
{
print(sum(magicsolution[i, ]))
print(sum(magicsolution[ ,i]))
}

print(sum(diag(magicsolution)))
print(sum(diag(t(magicsolution))))


#  magic 4 x 4, after 2140 steps, all sums 34: 
#    1   12   16    5
#   10   13    7    4
#    8    3    9   14
#   15    6    2   11

#  magic 5 x 5, after 27802 steps, all sums 65: 
#   24   12    1   22    6
#   10    9   16    5   25
#    2   18   21   13   11
#   14   19    4    8   20
#   15    7   23   17    3

#  magic 6 x 6, after 23241 steps, all sums 111: 
#   24   17   36    7    9   18
#    6   20   10   26   15   34
#   16    4   31   27   14   19
#   25   12    3   13   30   28
#    8   35    2   33   22   11
#   32   23   29    5   21    1

#  magic 7 x 7, after 74012 steps, all sums 175: 
#  22   34   20   23   24   37   15
#  31   16   46   43   17    8   14
#  39   47    9    2   29   48    1
#  36    6   21   26   41    7   38
#  32   10   12   35   28   40   18
#   4   49   25   19    3   30   45
#  11   13   42   27   33    5   44

#  magic 8 x 8, after 134901 steps, all sums 260: 
#  35   47   39   64   62   10    1    2
#  58   16   12   15   34   36   49   40
#   6   30    5   44   63   33   25   54
#  41   13   21   57   22   37   60    9
#  27   29   48    8   19   46   59   24
#  52   26   32    7   18   55   17   53
#   3   56   61   14   11   20   45   50
#  38   43   42   51   31   23    4   28