% Compare volumes between replica-exchange and no swaps.

clear;

nequil = 150;

volumes_allswap = textread('data-allswap/volumes.out');
volumes_noswap = textread('data-noswap/volumes.out');

% Determine number of states.
nstates = size(volumes_allswap,2);
nsamples = size(volumes_allswap,1);

% Determine statistical inefficiencies.
for i = 1:nstates
  i

  [t, g, Neff] = detect_equilibration(volumes_allswap(:,i));
  g_i_allswap(i) = g;
  V_i_allswap(i) = mean(volumes_allswap(t:end,i));
  dV_i_allswap(i) = std(volumes_allswap(t:end,i)) / sqrt(Neff);

  [t, g, Neff] = detect_equilibration(volumes_noswap(:,i));
  g_i_noswap(i) = g;
  V_i_noswap(i) = mean(volumes_noswap(t:end,i));
  dV_i_noswap(i) = std(volumes_noswap(t:end,i)) / sqrt(Neff); 
end

clf;

errorbar(V_i_allswap, 2*dV_i_allswap, 'r.');
hold on;
errorbar(V_i_noswap, 2*dV_i_noswap, 'g.');
legend('repex', 'standard');
xlabel('alchemical intermediate');
ylabel('average box volume (nm^3)');
title('Hamiltonian exchange vs independent simulations (95% CI)');

exportfig(gcf, 'average-volume-comparison.eps', 'width', 4, 'height', 3, 'color', 'cmyk');
system('epstopdf average-volume-comparison.eps');
