% solves the traveling salesman problem
% with random position of points 

n=50; %no. of points
T=2; %temperature
noitemax=800000;

xp=rand([n+1,1]); xp(n+1)=xp(1);
yp=rand([n+1,1]); yp(n+1)=yp(1);

 
x=zeros(n+1,1);
y=zeros(n+1,1);
x1=zeros(n,1);
y1=zeros(n,1);
x2=zeros(n,1);
y2=zeros(n,1);
 
path=1:n;

 
x1(1:n)=xp(1:n); x2(1:n)=xp(2:n+1);
y1(1:n)=yp(1:n); y2(1:n)=yp(2:n+1);
cost=sum(sqrt((x1-x2).^2+(y1-y2).^2))
evolcost=zeros(noitemax,1);
evolT=zeros(noitemax,1);
noitevec=1:noitemax;
gam=1.3; %cooling exponent
 

    
for (noite=1:noitemax)
    %T=40/floor(log(noite+1)+1);
    if (mod(noite,1000)==1) 
        
        T=1*(10000^gam)/((10000^gam)+(noite^gam)); 
       
        %clf reset;
        x(1:n)=xp(path);
        y(1:n)=yp(path);
        x(n+1)=x(1); y(n+1)=y(1);
        plot([1:n+1],x,'-ro','MarkerFaceColor','black', 'MarkerSize',3);
        plot(x,y,'-ro','MarkerFaceColor','black', 'MarkerSize',3);
        title('Salesman''s route','FontSize',10);

        hold on;        
        axis equal; xlim([0,1]); ylim([0,1]);
        drawnow;
        if (noite==1) waitforbuttonpress;
        end
         
        hold off;
 
    end
    
  
    d1=randi(n); d2=randi(n-1);
    if (d2==d1) 
        d2=n;
    end
    ind1=min(d1,d2); 
    ind2=max(d1,d2);
     
 
    
    %reverse the path between ind1 and ind2
    pathn=path;
    pathn(ind1:ind2)=path(ind2:-1:ind1);
    xn=xp(pathn);
    yn=yp(pathn);
    xn(n+1)=xn(1); yn(n+1)=yn(1);
    x1(1:n)=xn(1:n); x2(1:n)=xn(2:n+1);
    y1(1:n)=yn(1:n); y2(1:n)=yn(2:n+1);
    cost2=sum(sqrt((x1-x2).^2+(y1-y2).^2));
    deltacost=cost2-cost;
    prob=exp(-deltacost/T);
    
    if ((deltacost<0)||(rand<prob))
        path=pathn;
        cost=cost2;
    end 
        
 
    if (mod(noite,10000)==0)
        fprintf('%d: T=%f, cost=%f\n', noite,T,cost);
    end
    evolcost(noite)=cost;
    evolT(noite)=T;
end
figure;
 
fprintf('cost=%f\n', cost);
plot(noitevec,evolcost, 'LineWidth',1.5, 'Color', 'red' );
hold on;
 
plot(noitevec,evolT, 'LineWidth',1.5, 'Color', 'blue' );
xlabel('time', 'Interpreter', 'Latex')
ylabel('cost', 'Interpreter', 'Latex');
title('Evolution of cost','FontSize',10);


