traveling salesman problem
来週までに、inversionとtranslationを行う関数を作っておいてください。

# initialisation, etc.
n.city=20
location=matrix(runif(n.city*2,min=0,max=100),nrow=n.city)
route=sample(n.city) 
calc.Dist<-function(location,route) {
  n.city=length(route)
  total.Dist=0
  route=c(route,route[1])
  for (i_city in 1:n.city){
    total.Dist=total.Dist+dist(location[c(route[i_city],route[i_city+1]),])
  }
  return(total.Dist)
}

# plotting results
plot(rbind(location[route,],location[route[1],]),type='o',pch=20,cex=2.5, col='red',
  xlab='location @X',ylab='location @Y',main='Initial route')
plot(rbind(location[best.route,],location[best.route[1],]),type='o',pch=20, 
  col='magenta',cex=2.5,xlab='location @X',ylab='location @Y',main='Optimized route')

#switching 
TSP.switch<-function(route) {
  two.cities=sample(route,2,replace=F)
  route[two.cities]=route[rev(two.cities)]
  return(route)
}

# inversion
TSP.inv<-function(route,inv.length) { 
  inv.begin=sample(route,1)
  inv.end=min(inv.begin+inv.length-1,length(route))
  route[inv.begin:inv.end]=rev(route[inv.begin:inv.end])
  return(route)
}

# translation
TSP.trans<-function(route,tr.length) {
  trP1=sample(route,1)
  tr.vec=route[trP1:min(length(route),trP1+tr.length-1)]
  temp.vec=route[-(trP1:min(length(route),trP1+tr.length-1))]
  trP2=sample(1:length(temp.vec),1)
  if (trP2==length(temp.vec)) {
    new=c(temp.vec,tr.vec)
  } else {
    new=c(temp.vec[1:trP2],tr.vec,temp.vec[(trP2+1):length(temp.vec)])
  }
return(new)
}

# demo
TSP.demo<-function(n.city=20, maxItr=1000) {
  location=matrix(runif(n.city*2,min=0,max=100),nrow=n.city)
  route=sample(n.city) 
  ## param. for simulated annealing - change values if necessary 
  C=1;eta=0.99;TEMP=50; 
  ##
  optDist=calc.Dist(location,route)
  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,sample(2:(round(n.city/2)),1))}
    new.dist=calc.Dist(location,new.route)
    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=TEMP*eta
  }
  par(mfrow=c(1,2)) 
  plot(rbind(location,location[1,]),type='o',pch=20,cex=2.5, col='red',
    xlab='location @X',ylab='location @Y',main='Initial route')
  plot(location[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)
}