#Modified from: http://xianblog.wordpress.com/2012/03/17/simulated-annealing-for-sudokus-2
#Calculates the error from a proposed solution
target=function(s){
tar=sum(apply(s,1,duplicated)+apply(s,2,duplicated))
for (r in 1:9){
bloa=(1:3)+3*(r-1)%%3
blob=(1:3)+3*trunc((r-1)/3)
tar=tar+sum(duplicated(as.vector(s[bloa,blob])))
}
return(tar)
}
cost=function(i,j,s){
#constraint violations in cell (i,j)
cos=sum(s[,j]==s[i,j])+sum(s[i,]==s[i,j])
boxa=3*trunc((i-1)/3)+1;
boxb=3*trunc((j-1)/3)+1;
cos+sum(s[boxa:(boxa+2),boxb:(boxb+2)]==s[i,j])
}
#Creates a starter suggestion substituting ZEROS by sampled numbers!
entry=function(to_solve){
s=to_solve
pop=NULL
for (i in 1:9)
pop=c(pop,rep(i,9-sum(to_solve==i)))
s[s==0]=sample(pop)
return(s)
}
# generates a step following simulated annealing strategy!
move=function(tau,s,con){
pen=(1:81)
for (i in pen[con==0]){
a=((i-1)%%9)+1
b=trunc((i-1)/9)+1
pen[i]=cost(a,b,s)
}
#It seems that Evolutionary Algorithms or Tabu Search would work better
#as it is statying a lot of time in the same solution!
#Maybe simply storing the best from the past and trying go back after stuck!
wi=sample((1:81)[con==0],2,prob=exp(pen[(1:81)[con==0]]))
prop=s
prop[wi[1]]=s[wi[2]]
prop[wi[2]]=s[wi[1]]
if (target(prop) < target(s))
s=prop
else if (runif(1)<exp((target(s)-target(prop)))/tau)
s=prop
return(s)
}
run_annealing=function(num_it, num_try, tau, to_solve, temp_degradation=0.99, draw=FALSE, point_color="black") {
s=entry(to_solve)
for(t in 1:num_it) {
for(v in 1:num_try) s=move(tau,s,to_solve)
tau=tau*temp_degradation
if (draw) points(t,target(s),col=point_color, cex=.3,pch=19)
if(target(s)==0) break()
}
}
#Too easy one
s=matrix(0,ncol=9,nrow=9)
s[1,c(6,8)]=c(6,4)
s[2,c(1:3,8)]=c(2,7,9,5)
s[3,c(2,4,9)]=c(5,8,2)
s[4,3:4]=c(2,6)
s[6,c(3,5,7:9)]=c(1,9,6,7,3)
s[7,c(1,3:4,7)]=c(8,5,2,4)
s[8,c(1,8:9)]=c(3,8,5)
s[9,c(1,7,9)]=c(6,9,1)
plot(c(0,0),col="white",xlim=c(1,100),ylim=c(0,target(s)),xlab="iterations",ylab="penalty")
run_annealing(100,100,100,s,0.99,TRUE,"blue")
run_annealing(100,100,50,s,0.99,TRUE,"red")
run_annealing(100,100,100,s,0.999,TRUE,"green")
run_annealing(100,100,50,s,0.999,TRUE,"black")