% Simple 1D model system illustrating single large perturbation followed by many steps of Metropolis MC propagation in a harmonic potential.
% JDC 2011-08-12

% Define potential.
potential = @(x) (1/2)*(x.^2);
dx_plt=0.1;
xplt=-4:dx_plt:4;
% Define run parameters.
nsteps_ncmc = 200000; % number of NCMC steps
dx_ncmc = 10.0; % stddev of outer NCMC trial moves
nsteps_mmc = 50; % number of Metropolis MC steps per NCMC step
dx_mmc = 1; % stddev of inner MMC trial moves

% Initialize sampler.
x_current = 0.0; % sampler location
xhist = zeros(nsteps_ncmc,1); % xhist(i) is position of sampler at iteration i
naccepted_ncmc = 0; % number of accepted NCMC moves

% Sample.
xhist=[];
count=0; first_time=1;
for ncmc_iteration = 1:nsteps_ncmc
  % Make perturbation, accumulating work.
  x_proposed = x_current + dx_ncmc*randn();  
  work = potential(x_proposed) - potential(x_current);

  % Take propagation steps.
  naccepted_mmc = 0; % number of inner MMC steps accepted
  for mmc_step = 1:nsteps_mmc
    x_proposed_mmc = x_proposed + dx_mmc*randn();
    du = potential(x_proposed_mmc) - potential(x_proposed);
    if (du < 0.0) || (rand() < exp(-du))
      % Accept.
      x_proposed = x_proposed_mmc;
      naccepted_mmc = naccepted_mmc + 1;
    else
      % Reject.
    end
  end

  % Accept or reject NCMC.
  if (work < 0.0) || (rand() < exp(-work))
    % Accept.
    x_current = x_proposed;
    naccepted_ncmc = naccepted_ncmc + 1;
  else
    % Reject.
  end

  % Store.
  count=count+1;
  xhist(count) = x_current;
  %xhist(ncmc_iteration) = x_current;
  
  % Display.
  if(mod(ncmc_iteration,5000)==1)
    disp(sprintf('Iteration %d | work %8.3f | accepted %8.1f%% of inner MMC steps', ncmc_iteration, work, (naccepted_mmc/nsteps_mmc)*100.0));
    disp(sprintf('Accepted %.3f%% of NCMC moves.', (naccepted_ncmc/ncmc_iteration)*100.0));
    disp(sprintf('Accepted %.3f%% of NCMC moves.', (naccepted_ncmc/nsteps_ncmc)*100.0));
    clf        ; 
    nbins = 50 ;  % number of histogram bins
    N_tmp = histc(xhist, xplt)          ;
    x=xplt;  N_tmp = histc(xhist, xplt) ;
    if(first_time==1)
        first_time=0                    ;
        N=N_tmp                         ;
    else,       N = N + N_tmp; ,    end
    
    clear xhist; count=0;
    
    hold on;
    
    Px = exp(-potential(x)); Px = Px * max(N) / max(Px);
    Px=Px/(sum(Px))/dx_plt;
    N_n=N/sum(N)/dx_plt;
    plot(x,N_n,'k-',x,Px,'r-')
    set(gca,'YTick', []);
    xlabel('x');
    ylabel('P(x)');
    legend('observed', 'expected');
    drawnow;
  end
  
end

% Display.
disp(sprintf('Accepted %.3f%% of NCMC moves.', (naccepted_ncmc/nsteps_ncmc)*100.0));
clf;
nbins = 50;  % number of histogram bins
[N,x] = hist(xhist, nbins);
hold on;
Px = exp(-potential(x)); Px = Px * max(N) / max(Px);
plot(x,N,'k-',x,Px,'r-')
set(gca,'YTick', []);
xlabel('x');
ylabel('P(x)');
legend('observed', 'expected');

