function [cbest,tbest]=SA_TSP(Tmax,Tmin)
// Calling example: SA_TSP(10000,0.001) - user specified data (mouse clicking)
//[x,n]=gendata(); // user selected coordinates 
mat=csvRead('eil51.tsp',' ')';   //read data from file
x=mat(2:3,:);
n=size(x,2);
t=init(n);     //initial tour
c=cost(x,t);   //evaluation of the initial tour 
tbest=t;
cbest=c;
T=Tmax;
while T>Tmin
    tnou=transform(t);   // tour transformation
    cnou=cost(x,tnou);
    disp(cnou,'cost=');
    a=accept(c,cnou,T);
    if a==1
        t=tnou;
        c=cnou;
    end
    if c<cbest
        cbest=c;    
        tbest=t;
    end;    
    T=T*0.9995;   // new temperature
end    
//plot the tour
clf
plot(x(1,:),x(2,:),'ko');
for i=1:n-1
    plot([x(1,tbest(i)),x(1,tbest(i+1))],[x(2,tbest(i)),x(2,tbest(i+1))],'r*-');
end
plot([x(1,tbest(n)),x(1,tbest(1))],[x(2,tbest(n)),x(2,tbest(1))],'r*-');
endfunction

// user selected coordinates
function [x,n]=gendata()
   // data input
   clf;
   plot2d(0,0,0,rect=[0,0,1,1]);
   xgrid(3);
   p=locate(-1,1);   
   n=size(p,2); 
   x = p;
endfunction  

// tour initialization
function t=init(n)
    for i=1:n
        t(i)=i;
    end   
//    ninv=grand(1,1,"uin",1,n);
//    for i=1:ninv
//        i1=grand(1,1,"uin",1,n);
//        i2=grand(1,1,"uin",1,n);
//        aux=t(i1); t(i1)=t(i2); t(i2)=aux;
//    end     
    t=grand(1,'prm',t); 
endfunction 

//cost computation
function c=cost(x,t)
    n=length(t);
    c=0;
    for i=1:n-1
        c=c+sqrt((x(1,t(i))-x(1,t(i+1)))^2+(x(2,t(i))-x(2,t(i+1)))^2);
    end    
    c=c+sqrt((x(1,t(n))-x(1,t(1)))^2+(x(2,t(n))-x(2,t(1)))^2);    
endfunction    

// 2-opt transformation
function tnou=transform(t)
n=length(t);
i1=grand(1,1,"uin",1,n); i2=grand(1,1,"uin",1,n);
if i1>i2
    aux=i1; i1=i2; i2=aux;
end
tnou=t
for i=1:i1-1
    tnou(i)=t(i);
end
for i=i1:i2
    tnou(i)=t(i1+i2-i);
end
for i=i2+1:n
    tnou(i)=t(i);
end    
endfunction

// acceptance function
// a=1 - the new tour is accepted
// a=0 - the new tour is not accepted
function a=accept(c,cnou,T)
if cnou<c
    a=1;
else
    p=exp(-(cnou-c)/T);
    if grand(1,1,"unf",0,1)<p
        a=1;
    else
        a=0;
    end
end    
endfunction

[cbest,tbest]=SA_TSP(5000,0.001)
