function annealing()
%  Deze functie-subroutine visualiseert de werking van 
%  het Simulated Annealing proces voor het vinden van 
%  het minimum (en bijbehorend minimalizerend punt) van 
%  een reeel-waardige functie.
%  
%  Een functie F (zie Fig.1) wordt gedefinieerd op een verza- 
%  meling van n equidistante gekozen punten t_k tussen 0 en 1.
%  Gevonden moet worden 'n t_k waarin F zijn minimale waarde 
%  aanneemt en de waarde w van dit minimum.
%  (F is een combinatie van random gekozen sinussen.)
%
%  Met Simulated Annealing wordt (middels een Markov proces)                
%  een rij (x_n) van "kandidaat" oplossingen t_k gegenereerd.
%
%  In Figuur 1 wordt de grafiek van F geplot plus
%  in iedere 10-de stap de "kandidaat" oplossing in die stap:
%  als in de n-de stap x_n=t_k die "kandidaat" oplossing is, 
%  dan markeert het rode bolletje het punt (x_n,F(x_n)) 
%
%  In Fig.2 en 3 is voor iedere t_k de kans p_k uitgezet dat
%  het proces uiteindelijk het punt t_k kiest als oplossing.
%  Het rode bolletje in deze twee figuren markeert de laatste              
%  x_laatst=x_n, het groene bolletje markeert de beste
%  x_best=x_n (indien de beste beter is dan de laatste: 
%            dus als F(x_best)<F(x_laatst)).
%  Fig.2 en 3 laten hetzelfde zien alleen op een andere schaal.

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


% Print bovenstaande toelichting
help annealing



close all
clear all

global FORPRINT FORSLIDES
FORSLIDES=1;
AANTALCHECKS=0;


n=1000; N=1:n; 
% Kies het aantal buurpunten
L=100;

% bepaal de te minimaliseren functie F
F=rand(1,n); alpha=1.2;
S=0;
for j=N
  S=S+sin(2*j*pi*N/n)*(F(j)/j^alpha);
end
F=S;
% F is gedefineerd op n punten t_k tussen 0 en 1
N=N/n;
% Plot F
plot(N,F,'*')

% Kies de begintemperatuur T
T=1;

% Kies een begin punt x_0
x=round(n/5); x0=x;
i=1; k=0;

%  Voer het simulated annealing algorithme uit tot dat T voldoende klein is
while (T>0.002)
  % selecteer een random y uit de L/2 buurpunten links of rechts
  z=[x-L:x-1,x0,x+1:x+L];
  y=z(randi(size(z,2)));
  % maak er een periodiek probleem van
  if y<1, y=y+n; end
  if y>n, y=y-n; end

  % Hoeveel potentiele oplossingen hebben we al gecheckt?
  AANTALCHECKS=AANTALCHECKS+1;
  
  % De simulated annealing stap
  if F(y)<F(x)
    x=y;
  else
    p=exp((F(x)-F(y))/T);
    if rand(1,1)<p, x=y; end
  end

  % onthou de beste oplossing tot zo ver.
  if F(x)<F(x0), x0=x; end
  
  % verklein T iedere 10-de stap
  % plot de gevonden oplossing iedere 10-de stap
  if k==10; k=0; T=T*0.96;
    clf, plot(N,F,'*'), hold on
    MyPlot(N(x),F(x),'ro6')
    xlabel(sprintf('T=%g',T))
    pause(0.1)
  end
  k=k+1;
end

%% Resultaten plus statistiek van het proces
J=find(F<=F(x0));
fprintf('\nDe minimale waarde van F (gevonden met simulated annealing)')
fprintf('\n                 is %g in t=%g\n',F(x0),N(x0))
[f,j]=min(F);
fprintf('Echte minimale waarde (gevonden met \"brute force\")')
fprintf('\n                 is %g in t=%g.\n',f,N(j))
if size(J,2)==1
  fprintf('De oplossing is uniek en is gedetecteerd.\n')
else
  fprintf('Behoort bij de %g%% kleinste.\n',100*size(J,2)/n)
end
fprintf('Gevonden in %0.3g trials, uit %0.3g mogelijkheden.\n\n',...
   AANTALCHECKS, n)


% Bepaal en plot de kansverdeling van de dominante eigenvector
% voor de eindtemperatuur T
FF=exp(-F/T); S=sum(FF); FF=FF/S;

% Plot de kansverdeling
figure(2)
MyPlot(N(x0),FF(x0),'go6'); hold on
MyPlot(N(x),FF(x),'ro6'); 
plot(N,FF,'b.'); 
ylabel('p_k, waarbij p_k de kans is dat t_k gevonden wordt')
title(sprintf('De dominante eigenvector bij T=%g',T))
xlabel('t_k')

figure(3)
plot(N,log10(FF),'b.'), hold on
MyPlot(N(x0),log10(FF(x0)),'go6')
MyPlot(N(x),log10(FF(x)),'ro6')
ylabel('log_{10}(p_k), waarbij p_k de kans is dat t_k gevonden wordt')
title(sprintf('De dominante eigenvector bij T=%g',T))
xlabel('t_k')

% Return de beste gevonden waard
x=x0;

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


