# salesman problem

n=20  # number of cities 
L=10  # side of square where cities are
x=L*runif(n) # x-coordinates of cities
y=L*runif(n) # y-coordinates of cities
path=seq(1:n) # path of a salesman

plot(x,y,xlim=c(0,L),ylim=c(0,L),type="b")

dist=function(path){
  #d=0;
    #for (i in 2:length(path)){
    # x1=x[path[i-1]]
    # y1=y[path[i-1]]
    # x2=x[path[i]]
    # y2=y[path[i]]
    # d=d+sqrt((x1-x2)^2+(y1-y2)^2)
    #}
  k=length(path)
  xx=x[path]; yy=y[path]
  d=sum(sqrt((xx[1:(k-1)]-xx[2:k])^2+(yy[1:(k-1)]-yy[2:k])^2))
  d;
 }

tmin=0.01 # admissible minimum of temperature
alpha=0.999  # multiplicative coefficient for fast simulated annealing T[i+1]=alpha*T[i]
T=tini=500    # starting temperature
E=dist(path)

while(T >= tmin){
  pair=sample.int(n,2);
  nmin=min(pair)
  nmax=max(pair)
  newpath=path
  newpath[nmin:nmax]=rev(path[nmin:nmax])
  Enew=dist(newpath)
  if (T*log(runif(1)) < E-Enew){E=Enew; path=newpath}
  T=T*alpha
  }

plot(x[path],y[path],xlim=c(0,L),ylim=c(0,L),type="b")
dist(path)

prop=function(path){
  pair=sample.int(n,2)
  nmin=min(pair)
  nmax=max(pair)
  proppath=path
  proppath[nmin:nmax]=rev(path[nmin:nmax])
  proppath
}

install.packages("ConsPlan")