function annealingTSP()
%  Deze functie-subroutine lost het handelsreizigersprobleem
%  (traveling salesman problem [tsp]) op met simulated annealing.
%  
%  Stel S_0, S_1, ...,  S_N is een N+1 tal steden en stel dat
%  er een tabel beschikbaar is die tussen ieder tweetal van
%  deze steden de afstand tussen die steden geeft (via de weg).
%  Een rondreis tussen deze steden doet iedere stad (in een of
%  andere volgorde) exact een keer aan: als pi een permutatie
%  is van 1,...,N dan is de reis van S_0 naar S_pi(1) 
%  vervolgens naar S_pi(2), dan naar S_pi(3), enzovoort tot 
%  naar S_pi(N) en tenslotte terug naar S_0 een rondreis. 
%  Iedere permuatie defineert een rondreis. De lengte L(pi) van
%  de rondreis die gedefieerd is door de permutatie pi is de
%  totale afstand die wordt afgelegd in die rondreis (de som
%  van de afstanden tussen de steden die achtereenvolgens in de
%  reis bezocht worden. Voor een gekozen rondreis kan deze
%  lengte gemakkelijk middels de beschikbare tabel uitgerekend
%  worden).
% 
%  We proberen de permutatie pi te vinden waarvoor L(pi)
%  minimaal is (dit is het traveling salesman problem).
%
%  Met Simulated Annealing wordt (middels een Markov proces)                
%  een rij (x_n) van "kandidaat" oplossingen pi gegenereerd.
%  De eerste rondreis x_0 is een random gekozen permutatie.
%  De gevonden "oplossing" x_laatst heeft een zeer grote kans 
%  om de correcte oplossing te zijn (de lengte van de gevonden
%  rondreis zal zeer waarschijnlijk de kleinste lengte zijn of er zeert 
%  dich bij in de buurt liggen). Omdat het aantal permutaties zeer
%  groot is er overigens geen manier om uit te rekenen hoe groot 
%  dat die kans daadwerkelijk is. 
%
%  In deze routine worden random N+1 punten p_0, p_1, ..., p_N
%  gekozen in het eenheidsvierkant (het vierkant stelt een
%  landkaart voor, en ieder punt p_j een stad. De afstand
%  tussen een tweetal punten is de afstand tussen de
%  corresponderende steden. Alleen maar om gemakkelijk 
%  voorbeelden te kunnen genereren reizen hier in
%  "vogelvlucht"). Afstanden worden gemeten in "vogelvlucht".
%  (Voor het algoritme gaat het om de afstand tussen twee
%  punten en is het niet van belang hoe we aan die afstanden
%  komen: het algoritme heeft alleen een afstandstabel nodig).
%  
%  In Figuur 2 wordt in stap n de beste oplossing tot die stap
%  getekend.


%  Gerard Sleijpen
%  Copyright (c) 28-9-2010

help annealingTSP

close all
clear all

global FORPRINT FORSLIDES
FORSLIDES=1;
global AANTALCHECKS
AANTALCHECKS=0;

N=31;
% Kies random punten in het eenheidsvierkant
X=rand(1,N);
Y=rand(1,N);
% Bereken  de afstandstabel voor deze  verzameling punten
S=afstandsmatrix(X,Y);

% Kies random een begin rondreis (dwz een permutatie)
P=randperm(N-1)+1;
% P=2:N;
P=[1,P,1];

% Bereken de lengte van de rondreis
L=lengtereis(X,Y,P); L0=L;
plotreis(X,Y,P,L,1)
drawnow

T=2;
x=P; x0=P; Lx=L; Lx0=Lx;
k=0; m=0;
PLOT=1; ok=1; figure(2)
tic

% stop als de temperatuur T lager is dan 0.005
while (T>0.005)

  %%% vind een "buur" rondreis
  m=1-m; 
  I=floor((N-1)*rand(1,2)-eps)+2;
  if m
    if I(1)<I(2), I=I(1):I(2); else, I=I(2):I(1); end
  end
  y=x; y(I)=flipdim(y(I),2);

  Ly=lengtereis(X,Y,y);
 
  %%% simulated annealing
  if Ly<Lx
    %%% kies de betere
    x=y; Lx=Ly;
  else
    %%% kies ook eventueel een slechtere
    p=exp((Lx-Ly)/T);
    if rand(1,1)<p, x=y; Lx=Ly; end
  end
 
  %%% bewaar het beste resultaat tot nu toe
  if  Lx<Lx0, 
     x0=x; Lx0=Lx; ok=1; AS=AANTALCHECKS;
  end
  
  %%% verlaag de temperatuur na 3000 trials
  if k>3000; k=0; T=T*0.98; end

  %%% plot de beste oplossing tot dusver
  if  PLOT  & (ok | k==0)
      plotreis(X,Y,x0,Lx0), ok=0; 
      xlabel(sprintf('Temperatuur=%g',T)), drawnow
  end
  k=k+1;
end
toc

% Geef de statistiek resulaten van het proces
fprintf('Oplossing gevonden in %0.3g trials, uit %0.3g mogelijkheden\n',...
   AANTALCHECKS, factorial(N))
fprintf('Na %0.3g trials was er geen verbetering meer.\n\n',AS)

% Plot de gevonden "kortste" rondreis
plotreis(X,Y,x0,Lx0)
title(sprintf('Oorspronkelijke lengte=%g, huidige Lengte=%g',...
   L0,lengtereis(X,Y,x0)))


return
%===============================================
function plotreis(X,Y,P,L,ok)

    clf
    MyPlot(X,Y,'bo6'); hold on
    MyPlot(X(1),Y(1),'ro6');
    plot(X(P),Y(P),'b-')
    title(sprintf('Rondreis langs %d steden. Lengte reis=%g',size(X,2),L))

return
%===============================================
function L=lengtereis(X,Y,P);

global AANTALCHECKS
AANTALCHECKS=AANTALCHECKS+1;

  L=0;
  for j=1:size(X,2)
    p=P(j); q=P(j+1);
    L=L+sqrt((X(q)-X(p))^2+(Y(q)-Y(p))^2);
  end
return
%===============================================
function S=afstandsmatrix(X,Y)
  N=size(X,2);
  S=zeros(N,N);
  for i=2:N
    S(i,i)=0;
    for j=1:i-1
      S(i,j)=sqrt((X(i)-X(j))^2+(Y(i)-Y(j))^2);
      S(j,i)=S(i,j);
    end
  end
return


%===============================================
%===============================================
%===============================================

function varargout=MyPlot(X,Y,str1,str2)
% Routine om dikke punten te plotten
%
% str1 voor de line specificaties, str2 voor de markers
% str1=cml, c=color, m=linetype, l=linewidth  (default 2)
% str2=cml, c=color, m=markertype, l=markersize (default 6)

   global FORSLIDES

   if isempty(FORSLIDES), FORSLIDES=0; end

   if nargin < 3, str1='b-'; end
   if nargin < 4, str2='';  end

   if ~FORSLIDES
     if (isempty(str1) & isempty(str2)), p=plot(X,Y); 
     else
       if ~isempty(str1)
         [c,l,z]=PlotSign(str1); str1=[c,l];
         p=plot(X,Y,str1); 
       end
       if ~isempty(str2)
         [c2,l2,z2]=PlotSign(str2); str2=[c2,l2];
         p=plot(X,Y,str2); 
       end
     end
     if nargout==1, varargout(1)={p}; end
     return
   end

   cmd1=''; cmd2='';
   if ~isempty(str1)
      [c,l,z]=PlotSign(str1,2); str1=[c,l];
      cmd1=',''LineWidth'',z,''LineStyle'',l,''Color'',c';
   end
   if ~isempty(str2)
     [c2,l2,z2]=PlotSign(str2,6); str2=[c2,l2];
     cmd2=',''Marker'',l2,''MarkerFaceColor'',c2';
     cmd2=[cmd2,',''MarkerEdgeColor'',c2,''MarkerSize'',z2'];
   end
   cmd=['p=line(X,Y',cmd1,cmd2,');'];
% cmd, pause
   eval(cmd)
   if nargout==1, varargout(1)={p}; end
return

%=============================================
function [c,l,s]=PlotSign(str,s)

  global FORSLIDES

  if nargin < 2
    s=0.5; if nargin < 1, c='b'; l='-'; return, end
  end

  c=str(1); l=str(2:length(str)); 
  ts='0123456789'; K=[];
  for s1=ts, K=[K,findstr(l,s1)]; end
  if ~isempty(K)
    k1=min(K); s=str2num(l(k1:max(K))); l=l(1:k1-1);
  end
  if c=='g', c=[0,0.4,0]; if FORSLIDES==2, c=[0,0.8,0]; end, end
  if c=='w', c=[0.95,0.95,0.95]; end

return



