% Run model system simulation and analysis.

% PARAMETERS
xmin = -1.75;
xmax = +1.75;

kT = 1.0; % thermal energy
beta = 1.0 / kT;

ntemps = 15;
%temperature_k = linspace(1.0, 10.0, ntemps);
temperature_k = 1.0*(10.0/1.0).^linspace(0,1,ntemps);
beta_k = temperature_k .^ (-1);

% Define potential
%K = 10;
%potential = @(x) K * (x-1).^2 .* (x+1).^2 + 0.1*x; % asymmetric double well

% Plot potential and stationary distribution at multiple temperatures.
npoints = 1000;
dx = (xmax - xmin) / (npoints-1);
xpoints = linspace(xmin, xmax, npoints);
Upoints = potential(xpoints);

% Clear figure
clf;

% Plot potential.
subplot(2,4,1);
plot(xpoints, Upoints, 'k-');
axis([xmin xmax min(Upoints) max(Upoints)]);
xlabel('x');
ylabel('U(x)');
box on;

subplot(2,4,5);
hold on;
Pmax = 0.0;
for k = 1:ntemps
  beta = beta_k(k);
  Peq = exp(-beta*Upoints);
  Peq = Peq / sum(Peq) * dx;
  plot(xpoints, Peq);
  Pmax = max(Pmax, max(Peq));
end
axis([xmin xmax 0 Pmax]);
xlabel('x');
ylabel('\pi(x)');
set(gca,'XTick',[]);
box on;

% Neighbor update.

% Run expanded ensemble simulation.
[x_t, k_t] = expanded_ensemble_simulation(beta_k, 'neighbor');

% Compute statistical inefficiencies.
g_k = statistical_inefficiency(k_t);
g_x = statistical_inefficiency(x_t);

subplot(2,4,2);
plot(k_t, 'r-');
axis([0 nshowiterations 0 ntemps+1]);
title(sprintf('neighbor, g_k = %.1f', g_k));
xlabel('iteration');
ylabel('temperature index');
box on;

subplot(2,4,6);
plot(x_t, 'k.');
axis([0 nshowiterations xmin xmax]);
title(sprintf('neighbor, g_x = %.1f', g_x));
xlabel('iteration');
ylabel('x');
box on;

% Run expanded ensemble simulation.
[x_t, k_t] = expanded_ensemble_simulation(beta_k, 'gibbs');

% Compute statistical inefficiencies.
g_k = statistical_inefficiency(k_t);
g_x = statistical_inefficiency(x_t);

subplot(2,4,3);
plot(k_t, 'r-');
axis([0 nshowiterations 0 ntemps+1]);
title(sprintf('gibbs, g_k = %.1f', g_k));
xlabel('iteration');
ylabel('temperature index');
box on;

subplot(2,4,7);
plot(x_t, 'k.');
axis([0 nshowiterations xmin xmax]);
title(sprintf('gibbs, g_x = %.1f', g_x));
xlabel('iteration');
ylabel('x');
box on;

% Collect data for different numbers of temperatures.
%ntemps_to_try = 5:5:50;
ntemps_to_try = [2 4 8 16 32 64];
g_k_gibbs = zeros(1,length(ntemps_to_try));
g_x_gibbs = zeros(1,length(ntemps_to_try));
g_k_neighbor = zeros(1,length(ntemps_to_try));
g_x_neighbor = zeros(1,length(ntemps_to_try));
for i = 1:length(ntemps_to_try)
  ntemps = ntemps_to_try(i);
  temperature_k = 1.0*(10.0/1.0).^linspace(0,1,ntemps);
  beta_k = temperature_k .^ (-1);

  % Run expanded ensemble simulation.
  [x_t, k_t] = expanded_ensemble_simulation(beta_k, 'gibbs');
  
  % Compute statistical inefficiencies.
  g_k_gibbs(i) = statistical_inefficiency(k_t);
  g_x_gibbs(i) = statistical_inefficiency(x_t);

  % Run expanded ensemble simulation.
  [x_t, k_t] = expanded_ensemble_simulation(beta_k, 'neighbor');
  
  % Compute statistical inefficiencies.
  g_k_neighbor(i) = statistical_inefficiency(k_t);
  g_x_neighbor(i) = statistical_inefficiency(x_t);
end

subplot(2,4,4);
plot(ntemps_to_try,g_k_neighbor,ntemps_to_try,g_k_gibbs);
legend('neighbor', 'gibbs');
ylabel('g_k');
xlabel('ntemps');

subplot(2,4,8);
plot(ntemps_to_try,g_x_neighbor,ntemps_to_try,g_x_gibbs);
legend('neighbor', 'gibbs');
ylabel('g_x');
xlabel('ntemps');

filename = 'model-expanded-ensemble.eps';
exportfig(gcf, filename, 'width', 7.5, 'height', 3.5, 'color', 'cmyk');
system(sprintf('epstopdf %s', filename));

save model_workspace.mat
