//==================================================================================
//
//  METROPOLIS ADJUSTED LANGEVIN ALGORITHM
//  ---------------------------------------
//
//  FIGURES 
//  0 = comparison between theoretical and empirical distributions
//  1 = trajectory x^n as a function of the step index n
//  2 = distribution of the increments x^{n+1}-x^n
//  3 = average of the observable under consideration as a function of the iteration step
//
//  WORK TO DO
//  1 -- compute the rejection rate as a function of alpha (for h = 1 and Beta = 2)
//  2 -- what is the optimal value to have the fastest sampling? 
//       For this question, it is worth considering the following observables: 
//       potential energy V(q) and position q
//  3 -- what happens as h is increased? as Beta is decreased?
//  4 -- implement another Metropolis move, for instance MALA but with a noise given by 
//       a uniform distribution between [-1/2, 1/2] rather than a Gaussian of variance 1
//
//  Last modified: Paris, July 3rd, 2017 (Gabriel STOLTZ)
// 
//==================================================================================


//-----------------------
//   POTENTIAL USED
//-----------------------

function y = V(x)
  y = h*((x-1).^2).*(x+1).^2;
endfunction

//-- derivative of the potential --
function y = nablaV(x)
  y = 2*h*( (x-1).*(x+1).^2 + (x-1).^2*(x+1) );
endfunction

//------------------------
//  VARIABLES TO DECLARE
//------------------------

h = 1;
Beta = 2;
Niter = 10000;
alpha = 0.3;
Rplot = 3;
Rmax = 5;
Nhisto = 30;
Ncomp = 100;

//------------------
// INITIALIZATIONS
//-----------------

x = zeros(Niter);
X = x;
x(1) = 1;
V_old = V(x(1));
V_new = V_old; 
E = zeros(Niter);
E(1) = V_new;

//-------------------------------------------
//  GENERATION OF THE STATISTICAL ENSEMBLE 
//-------------------------------------------

rejection = 0;
for i=2:Niter

  V_old = V(x(i-1));
 
  //--------------------
  //   PROPOSITION    
  //--------------------
  noise = rand(1,'normal');
  x_proposed = x(i-1) - nablaV(x(i-1))*alpha + sqrt(2*alpha/Beta)*noise;
  X(i) = x_proposed;
  V_new = V(x_proposed);
  pxy = exp(-0.5*noise^2); //-- this is the probability T(q^n,q^n+1)
  Noise = ( x(i-1) - x_proposed + nablaV(x_proposed)*alpha )/sqrt(2*alpha/Beta);
  pyx = exp(-0.5*Noise^2); //-- this is T(q^n+1,q^n)

  //--------------------------
  //   ACCEPTANCE / REJECTION
  //---------------------------
  p = exp(-Beta*(V_new - V_old))*pyx/pxy;
  if (p >= 1)
    x(i) = x_proposed;
  else 
    if (rand(1,'uniform') < p)
      x(i) = x_proposed;
    else
      x(i) = x(i-1);
      rejection = rejection + 1;
    end
  end

  //-- update of the observables
  E(i) = E(i-1) + x(i); //V(x(i));

end

//--------------------------------------------------------------
//  COMPARISON BETWEEN EMPIRICAL AND THEORETICAL DISTRIBUTIONS 
//--------------------------------------------------------------

//-- un-normalized Boltzmann distribution 
function z = mu(x)
  z = exp(-Beta*V(x));
endfunction

//-- plot of the empirical distribution 
histplot(Nhisto,x);

//-- plot of the theoretical distribution 
C = intg(-Rmax,Rmax,mu);
z = -Rplot:0.05:Rplot;
plot(z,mu(z)/C,'r');

//-- plot of the visited positions 
set("figure_style","new");
f = get("current_figure");
set("current_figure",1);
plot(x);

//-- difference between x^n+1 and x^n
set("figure_style","new");
f = get("current_figure");
set("current_figure",2);
histplot(Ncomp,x(2:Niter)-x(1:Niter-1)); 

//-- convergence of the observable
energie = E' ./ (1:Niter);
set("figure_style","new");
f = get("current_figure");
set("current_figure",4);
plot(energie);

//----------------------------
//     REJECTION RATE
//----------------------------

rejection = rejection/Niter