
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=1000;				% number of tries to take
ns=100;				% number of models to keep
T0=1.;
alpha=0.99;

% 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)));
ppdf=pdf*0+1;
ppdf=ppdf/(sum(sum(ppdf)));


Ta=input(' Please enter T_a ');
T=Ta;

% Random walk

% find an initial vector x



% Let s get moving

xcur=round([rand rand]*100);
is=1;
xa(1,1)=xcur(1);
xa(1,2)=xcur(2);

% Random walk (keep the best)

ii=1;
it=0;

best_model=xcur;
best_prob=0;

while it < np,
   it=it+1;
   ii=-1*ii;
   
   % change T
   T=alpha*T;
   
   % calculate current probability
      S=-T0*log(pdf(xcur(1),xcur(2))/ppdf(xcur(1),xcur(2)));
      smTc=ppdf(xcur(1),xcur(2))*exp( -S/T);
      
   % calculate current probability
   if ii>0,
      xnew=xcur;
      xnew(1)=ceil(rand*n); 
      S=-T0*log(pdf(xnew(1),xnew(2))/ppdf(xnew(1),xnew(2)));
      smTn=ppdf(xnew(1),xnew(2))*exp( -S/T);
      else
      xnew=xcur;
      xnew(2)=ceil(rand*n);if xnew(2)==0,xnew(2)==1;end
      S=-T0*log(pdf(xnew(1),xnew(2))/ppdf(xnew(1),xnew(2)));
      smTn=ppdf(xnew(1),xnew(2))*exp( -S/T);
   end
   
   % decide whether to accept the new model with Metropolis
   if smTn > smTc, 
      xcur=xnew;
      if smTn > best_prob; best_prob=smTn; best_model=xnew; end
      disp(sprintf(' Move accepted : %g > %g at T = %g',smTn,smTc,T));
      imagesc(pdf'), hold on,plot(xcur(1),xcur(2),'w+'),
      plot(best_model(1),best_model(2),'r+')
      drawnow,pause, hold off
   else
      x=rand;
      if x < smTn/smTc,
         xcur=xnew;
         if smTn > best_prob; best_prob=smTn; best_model=xnew; end
         disp(sprintf(' Move accepted : %g < %g (%g%%) at T = %g ',smTn,smTc,smTn/smTc*100,T));
         imagesc(pdf'), hold on,plot(xcur(1),xcur(2),'w+'),
         plot(best_model(1),best_model(2),'r+')
			drawnow,pause, hold off
      end
   end
   
   
end

   