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))