
clear
clf

% Monte Carlo sampling of a probability density function
% using the simulated annealing approach

% Heiner Igel, 2002

% for simplicity we use the "peaks" function. The corresponding
% section can simply be replaced with the probability from a
% proper inverse problem (e.g. hypocenter location, gravity inversion, etc. )

% Some global parameters
np=1000000;				% number of tries to take
ns=100;				% number of models to keep
T0=1.;
alpha=0.995;
T=T0;

% Let's define a pdf

n=100;				% number of dimension in x and y
pdf=peaks(n);
pdf=sqrt(pdf.*pdf);
pdf=pdf/(sum(sum(pdf)));
pdf=pdf/(max(max(pdf)));

surf(pdf), pause

% Random walk

% find an initial vector x

xcur=round([rand rand]*100);

% Let s get moving

is=1;
np=0;
xa(1,1)=xcur(1);
xa(1,2)=xcur(2);

imagesc(pdf'),

while is<=ns,

	np=np+1;   
   % make a random choice for the next move
   xnew=round([rand rand]*100);
   xnew(find(xnew==0))=xnew(find(xnew==0))+1;
   
   % compare probabilities
   Pcur=pdf(xcur(1),xcur(2));
   Pnew=pdf(xnew(1),xnew(2));
   
   T=alpha*T;
   
   if Pnew >= Pcur, 
      xcur=xnew;
      is=is+1;
      disp(sprintf(' Made the %i-th move to [%i,%i] ', is,xcur(1),xcur(2))) 
      
      % display move graphically
      
		imagesc(pdf'), hold on, plot(xnew(1),xnew(2),'k+'), hold off
      drawnow
      
   end
   
   if Pnew < Pcur, 
      deltaP=Pcur-Pnew;
      P=exp(-deltaP/T);
      test=rand;
      if test<=P, 
      xcur=xnew;
      is=is+1;       
        
      disp(sprintf(' Made the %i-th move to [%i,%i] P= %g ', is,xcur(1),xcur(2),P)) 
      imagesc(pdf'), hold on, plot(xnew(1),xnew(2),'k+'), hold off
      drawnow
      
   end
end

   Pa(np)=Pnew;
   xa(np,1)=xcur(1);
   xa(np,2)=xcur(2);

   
end


hold off

% Final sampling

[y,i]=max(Pa);  		% index of the best model

subplot(121)
imagesc(pdf'), hold on
xlabel('x')
ylabel('y')
plot(xa(1:ns,1),xa(1:ns,2),'w+')
plot(xa(i,1),xa(i,2),'w+','LineWidth',4)
hold off

% Show as a function of the iterations 
% a. the best model
% b. the current model

for i=1:np,
   best(i)=max(Pa(1:i));
end



subplot(122)
plot(1:np,best(1:np),'r-',1:np,Pa(1:np),'b-')
xlabel(' # of points visited ')
ylabel(' Current model (blue) and best model (red) ') 