% Make plots from data.

clear;
load model_system_data.mat

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

nshowiterations = 5000;

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

ntemps_index_to_plot = 4; % index to plot
temperature_k = data.gibbs{ntemps_index_to_plot}.temperature_k;
beta_k = data.gibbs{ntemps_index_to_plot}.beta_k;
ntemps = ntemps_to_try(ntemps_index_to_plot);

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

figure(1)
clf;

ny = 6;
nx = 1;

% Plot potential.
%subplot(ny,nx,1);
subplot('position', [0.10 0.25+5*0.7/6 0.85 0.7/6]);
plot(xpoints, Upoints, 'k-');
axis([xmin xmax min(Upoints) max(Upoints)]);
%xlabel('x');
ylabel('U(x)');
set(gca,'XTick',[]);
box on;

%subplot(ny,nx,2);
subplot('position', [0.10 0.25+4*0.7/6 0.85 0.7/6]);
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,'YTick',[]);
box on;

% Neighbor update.

simulation_data = data.neighbor{ntemps_index_to_plot};
x_t = simulation_data.x_t;
k_t = simulation_data.k_t;
g_x = simulation_data.g_x;
g_k = simulation_data.g_k;
tau_x = simulation_data.tau_x;
tau_k = simulation_data.tau_k;
% DEBUG
[tau_k, dtau_k] = autocorrelation_time_with_error(k_t);
[tau_x, dtau_x] = autocorrelation_time_with_error(x_t);

%subplot(ny,nx,3);
subplot('position', [0.10 0.15+3*0.7/6. 0.85 0.7/6]);
plot(k_t, 'r-');
axis([0 nshowiterations 0 ntemps+1]);
title(sprintf('neighbor (\\tau_k = %.0f\\pm%.0f, \\tau_x = %.1f\\pm%.1f)', tau_k, dtau_k, tau_x, dtau_x));
set(gca,'XTick',[]);
%xlabel('iteration');
ylabel('k');
set(gca,'YTick',[1 16]);
box on;

%subplot(ny,nx,4);
subplot('position', [0.10 0.15+2*0.7/6 0.85 0.7/6]);
plot(x_t, 'k.');
axis([0 nshowiterations xmin xmax]);
%title(sprintf('neighbor, g_x = %.1f', g_x));
set(gca,'XTick',[0 1000 2000 3000 4000 5000]);
set(gca,'YTick',[-1 0 +1]);
xlabel('iteration');
ylabel('x');
box on;

% Gibbs update.

simulation_data = data.gibbs{ntemps_index_to_plot};
x_t = simulation_data.x_t;
k_t = simulation_data.k_t;
g_x = simulation_data.g_x;
g_k = simulation_data.g_k;
tau_x = simulation_data.tau_x;
tau_k = simulation_data.tau_k;
% DEBUG
[tau_k, dtau_k] = autocorrelation_time_with_error(k_t);
[tau_x, dtau_x] = autocorrelation_time_with_error(x_t);

%subplot(ny,nx,5);
subplot('position', [0.10 0.05+1*0.7/6 0.85 0.7/6]);
plot(k_t, 'r-');
axis([0 nshowiterations 0 ntemps+1]);
title(sprintf('independence (\\tau_k = %.3f\\pm%.3f, \\tau_x = %.1f\\pm%.1f)', tau_k, dtau_k, tau_x, dtau_x));
set(gca,'XTick',[]);
%xlabel('iteration');
ylabel('k');
set(gca,'YTick',[1 16]);
box on;

%subplot(ny,nx,6);
subplot('position', [0.10 0.05+0*0.7/6 0.85 0.7/6]);
plot(x_t, 'k.');
axis([0 nshowiterations xmin xmax]);
%title(sprintf('Gibbs, g_x = %.1f', g_x));
set(gca,'XTick',[0 1000 2000 3000 4000 5000]);
set(gca,'YTick',[-1 0 +1]);
xlabel('iteration');
ylabel('x');
box on;

filename = 'model-example-trajectory.eps';
exportfig(gcf, filename, 'width', 3.5, 'height', 7.0, 'color', 'cmyk');
system(sprintf('epstopdf %s', filename));

%% Show efficiency vs number of temperatures.

figure(2);
clf;

ny = 2;
nx = 1;

subplot(ny,nx,1);
tau_k_neighbor = zeros(1,length(ntemps_to_try));
dtau_k_neighbor = zeros(1,length(ntemps_to_try));
tau_k_gibbs = zeros(1,length(ntemps_to_try));
dtau_k_gibbs = zeros(1,length(ntemps_to_try));
for i = 1:length(ntemps_to_try)
  tau_k_neighbor(i) = data.neighbor{i}.tau_k;
  tau_k_gibbs(i) = data.gibbs{i}.tau_k;

  % Compute error bars.
  k_t = data.neighbor{i}.k_t;
  [tau_k, dtau_k] = autocorrelation_time_with_error(k_t);  
  tau_k_neighbor(i) = tau_k;
  dtau_k_neighbor(i) = dtau_k;
  k_t = data.gibbs{i}.k_t;
  [tau_k, dtau_k] = autocorrelation_time_with_error(k_t);  
  tau_k_gibbs(i) = tau_k;
  dtau_k_gibbs(i) = dtau_k;
end
semilogx(ntemps_to_try, tau_k_neighbor, 'k.-', ntemps_to_try, tau_k_gibbs, 'r*-', 'markersize', 6);
hold on;
for i = 1:length(ntemps_to_try)
  semilogx(ntemps_to_try(i) * [1 1], tau_k_neighbor(i)+dtau_k_neighbor(i)*[-1 +1], 'k-');
  semilogx(ntemps_to_try(i) * [1 1], tau_k_gibbs(i)+dtau_k_gibbs(i)*[-1 +1], 'r-');
end
legend('neighbor', 'independence', 'location', 'northwest');
set(gca,'XTick',ntemps_to_try);
ylabel('\tau_k');
xlabel('number of temperatures');
oldaxis = axis;
axis([1.5 max(ntemps_to_try)*1.5 0 oldaxis(4)]);
box on;

subplot(ny,nx,2);
tau_x_neighbor = zeros(1,length(ntemps_to_try));
dtau_x_neighbor = zeros(1,length(ntemps_to_try));
tau_x_gibbs = zeros(1,length(ntemps_to_try));
dtau_x_gibbs = zeros(1,length(ntemps_to_try));
for i = 1:length(ntemps_to_try)
  tau_x_neighbor(i) = data.neighbor{i}.tau_x;
  tau_x_gibbs(i) = data.gibbs{i}.tau_x;

  % Compute error bars.
  x_t = data.neighbor{i}.x_t;
  [tau_x, dtau_x] = autocorrelation_time_with_error(x_t);  
  tau_x_neighbor(i) = tau_x;
  dtau_x_neighbor(i) = dtau_x;
  x_t = data.gibbs{i}.x_t;
  [tau_x, dtau_x] = autocorrelation_time_with_error(x_t);  
  tau_x_gibbs(i) = tau_x;
  dtau_x_gibbs(i) = dtau_x;
end
semilogx(ntemps_to_try, tau_x_neighbor, 'k.-', ntemps_to_try, tau_x_gibbs, 'r*-', 'markersize', 6);
hold on;
for i = 1:length(ntemps_to_try)
  semilogx(ntemps_to_try(i) * [1 1], tau_x_neighbor(i)+dtau_x_neighbor(i)*[-1 +1], 'k-');
  semilogx(ntemps_to_try(i) * [1 1], tau_x_gibbs(i)+dtau_x_gibbs(i)*[-1 +1], 'r-');
end
%legend('neighbor', 'Gibbs', 'location', 'northwest');
set(gca,'XTick',ntemps_to_try);
ylabel('\tau_x');
xlabel('number of temperatures');
oldaxis = axis;
axis([1.5 max(ntemps_to_try)*1.5 0 oldaxis(4)]);
box on;

filename = 'model-g-vs-ntemps.eps';
exportfig(gcf, filename, 'width', 3.5, 'height', 3.5, 'color', 'cmyk');
system(sprintf('epstopdf %s', filename));


