Traveling Salesman Problem
TSP.inv <- function(x, n){
x.length = length(x)
start.idx = sample(1:x.length,1)
end.idx = min(x.length, (start.idx + n -1))
x[start.idx:end.idx]=x[end.idx:start.idx]
return(x)
}
TSP.switch <- function(x){
x.length = length(x)
switch.idx = sample(1:x.length,2)
x[switch.idx] = x[rev(switch.idx)]
return(x)
}
TSP.trans <-function(x){
x.length = length(x)
trans.idx = sample(1:x.length,2)
new.x = x[-trans.idx[1]]
new.x = append(new.x, x[trans.idx[1]],after = (trans.idx[2]-1))
return(new.x)
}
TSP.init<-function(n_city=10) {
return(matrix(runif(n_city*2,1,100),nrow=n_city,ncol=2))
}
cities = TSP.init(10)
TSP.totDist<-function(route, cities) {
route=c(route,route[1]);
n.cities=nrow(cities);
totDist=0;
for (i.cities in 1:n.cities) {
totDist = totDist + dist(cities[route[i.cities:(i.cities + 1)],])
}
return(totDist)
}
TSP.totDist(route,cities)
TSP.demo<-function(n.city=20, maxItr=1000, parm.set = c(1,0.99, 50)) {
loc=TSP.init(n.city);
route=1:n.city
## param. for simulated annealing - change values if necessary
C=parm.set[1];
eta=parm.set[2];
TEMP=parm.set[3];
##
optDist=TSP.totDist(route,loc)
optHist=matrix(0,nrow=maxItr,ncol=(length(route)+1))
optHist[1,]=c(route,optDist)
for (i_loop in 2:maxItr) {
rand.op=sample(c('inv','sw','trans'),1,prob=c(0.75,0.125,0.125))
if (rand.op=='inv') {
new.route=TSP.inv(route,sample(2:n.city,1))
} else if (rand.op=='sw') {
new.route=TSP.switch(route)
} else {
new.route=TSP.trans(route)
}
new.dist=TSP.totDist(new.route,loc)
delta=new.dist-optDist
prob=1/(1+exp(delta/(C*TEMP)))
if (runif(1) < prob) {
route=new.route;optDist=new.dist;
}
optHist[i_loop,]=c(route,optDist);
TEMP=min(TEMP*eta,0.1)
}
par(mfrow=c(1,2))
plot(rbind(loc,loc[1,]),type='o',pch=20,cex=2.5, col='red',
xlab='location @X',ylab='location @Y',main='Initial route')
plot(loc[optHist[1000,c(1:n.city,1)],],type='o',pch=20, col='magenta',
cex=2.5,xlab='location @X',ylab='location @Y',main='Optimized route')
return(optHist)
}
# regression with SimAnneal.
reg.ERR<-function(b,X,y){
yhat<-X%*%b
return(sum((y-yhat)^2))
}
##
SA.reg <- function(maxItr,parm.set=c(C=1,eta=0.995,TEMP=100,sigma=100)){
set.seed(20)
nData = 100
x=matrix(rnorm(9*nData,mean=10,sd=2),ncol=9);X=cbind(rep(1,nData),x)
y=X%*%c(10,2,5,-3,-5,0,0,0,0,0)+rnorm(nData,mean=0,sd=2);
lm.sum = summary(lm(y~X[,2:10]))
C = parm.set[1]
eta = parm.set[2]
TEMP= parm.set[3]
sigma = parm.set[4]
b = rnorm(10,0,100)
optERR=reg.ERR(b,X,y)
optHist=matrix(0,nrow=maxItr,ncol=(10+1))
optHist[1,]=c(b,optERR)
for (i_loop in 2:maxItr) {
new.b = b + rnorm(10,0,sigma*TEMP)
new.ERR=Fun.reg(new.b, X,y)
delta=new.ERR-optERR
prob=1/(1+exp(delta/(C*TEMP)))
if (runif(1) < prob) {
b=new.b;
optERR=new.ERR;
}
optHist[i_loop,]=c(b,optERR);
TEMP=max(TEMP*eta,0.01)
}
print(optHist[maxItr,])
print(coef(lm.sum)[,1])
return(list(optHist=optHist, gt = coef(lm.sum)))
}
res=SA.reg(1e5,c(80,0.999,1,1))