% Generate all plots for alanine dipeptide system.

%% Temperature mixing statistics.
%figure;
%clf;
%Tij = textread('swap-statistics.txt');
%imagesc(Tij);
%axis square;
%colorbar;

scale = 1.3; % scale factor for exportfig

% Read temperatures.
temperatures = textread('temperatures.txt');
temperatures

%% Autocorrelation functions.

temperature_indices = [1 8 15 30 40]; % temperature indices to plot
%colors = 'rgbcmyk'; % colors to cycle through for plots
colors = [1.0 0.0 0.0;
	  0.0 1.0 0.0;
	  0.0 0.0 1.0;
	  1.0 1.0 0.0;
	  1.0 0.0 1.0;
	  0.0 1.0 1.0;
	  1.0 1.0 1.0;];
ncolors = size(colors,1);
figure(1);
clf; 
hold on;
box on;

% Reweighted correlation functions.
data = textread('correlation-reweighted.txt');
reweighted_temperatures = data(:,1);
C_lt = data(:,2:end);
[L,T] = size(C_lt);
data = textread('correlation-reweighted-error.txt');
dC_lt = data(:,2:end);
tvec = (0:(T-1)) * 0.5; % times at which correlation function is computed
handles = [];
for i = 1:length(temperature_indices)
  k = temperature_indices(i);
  l = 2*(k-1)+1;
  color = colors(mod(i, ncolors)+1,:);
  %plot(tvec, C_lt(l,:), color);
  %errorbar(tvec, C_lt(l,:), 2*dC_lt(l,:), color); % 95% confidence interval
  h = fill([tvec fliplr(tvec)], [C_lt(l,:)+2*dC_lt(l,:) fliplr(C_lt(l,:)-2*dC_lt(l,:))], 'k');
  set(h, 'FaceColor', color * 0.2 + [1 1 1] * 0.8);
  set(h, 'EdgeColor', color * 0.2 + [1 1 1] * 0.8);
  hreweight = plot(tvec, C_lt(l,:), '-', 'LineWidth', 2);
  set(hreweight, 'color', color);
  handles = [handles hreweight];  
end
xlabel('\tau / ps');
ylabel('C(\tau;\beta)');
axis([0 10 0 1])

% Individual temperature estimates.
data = textread('correlation.txt');
temperatures = data(:,1);
C_kt = data(:,2:end);
data = textread('correlation-error.txt');
dC_kt = data(:,2:end);
dC_kt(isnan(dC_kt)) = 0.0;
[K,T] = size(C_kt);
tvec = (0:(T-1)) * 0.5; % times at which correlation function is computed
for i = 1:length(temperature_indices)
  k = temperature_indices(i);
  color = colors(mod(i, ncolors)+1,:);
  %plot(tvec, C_kt(k,:), color);
  hsingle = errorbar(tvec, C_kt(k,:), 2*dC_kt(k,:), '.', 'MarkerSize', 6); % 95% confidence interval
  set(hsingle, 'color', color);
end
xlabel('\tau / ps');
ylabel('C(\tau;\beta)');
%legend([hsingle, hreweight], 'single temp', 'reweighting');
hlegend = legend(handles, '300 K', '334 K', '375 K', '493 K', '600 K');
set(hlegend, 'location', 'southwest');
%legend([hsingle, hreweight], 'single temp', 'reweighting');
axis([0 10 0 1])

% Generate PDF.
filename = 'alanine-dipeptide-correlation-functions.eps';
%print('-depsc', filename);
exportfig(gcf, filename, 'color', 'cmyk', 'width', scale*3.5, 'height', scale*2.5);
system(sprintf('epstopdf %s', filename));


%% Plot bias and variance.

% bias
figure(2);
clf;
X = textread('reweighting-bias.txt'); X = X(:,2:end); Y = textread('singletemp-bias.txt'); Y = Y(:,2:end);    
min_bias = min(min([Y X]));
max_bias = max(max([Y X]));
max_bias = max([max_bias abs(min_bias)]);
clim = [-max_bias +max_bias];
imagesc([tvec tvec], temperatures, [Y X], clim); axis square; colorbar;
title('bias');
xlabel('single temp | reweighting');
ylabel('temperature / K');
set(gca,'XTick',[]); 
%set(gca,'YTick',[]);
colormap(jet);
colormap(flipud([linspace(1,1,32)' linspace(0,1,32)' linspace(0,1,32)'; linspace(1,0,32)' linspace(1,0,32)' linspace(1,1,32)']));

% Generate PDF.
filename = 'alanine-dipeptide-bias.eps';
%print('-depsc', filename);
exportfig(gcf, filename, 'color', 'cmyk', 'width', scale*3.5 * 0.5, 'height', scale*3.5);
system(sprintf('epstopdf %s', filename));

% bias (full-column)
figure(2);
clf;
X = textread('reweighting-bias.txt'); X = X(:,2:end); Y = textread('singletemp-bias.txt'); Y = Y(:,2:end);    
min_bias = min(min([Y X]));
max_bias = max(max([Y X]));
max_bias = max([max_bias abs(min_bias)]);
clim = [-max_bias +max_bias];
imagesc([tvec tvec], temperatures, [Y X], clim); axis square; colorbar;
title('bias');
xlabel('single temp | reweighting');
ylabel('temperature / K');
set(gca,'XTick',[]); 
%set(gca,'YTick',[]);
colormap(jet);
colormap(flipud([linspace(1,1,32)' linspace(0,1,32)' linspace(0,1,32)'; linspace(1,0,32)' linspace(1,0,32)' linspace(1,1,32)']));

% Generate PDF.
filename = 'alanine-dipeptide-bias-large.eps';
%print('-depsc', filename);
exportfig(gcf, filename, 'color', 'cmyk', 'width', scale*7.5 * 0.5, 'height', scale*3.5);
system(sprintf('epstopdf %s', filename));

% variance
figure(3);
clf;
X = textread('reweighting-stddev.txt'); X = X(:,2:end); Y = textread('singletemp-stddev.txt'); Y = Y(:,2:end);
max_variance = max(max([Y.^2 X.^2]));
clim = [0 max_variance];
imagesc([tvec tvec], temperatures, [Y.^2 X.^2], clim); axis square; colorbar;
title('variance');
xlabel('single temp | reweighting');
set(gca,'XTick',[]); 
set(gca,'YTick',[]);
%ylabel('temperature / K');
colormap(flipud([linspace(1,1,64)' linspace(0,1,64)' linspace(0,1,64)']));

% Generate PDF.
filename = 'alanine-dipeptide-variance.eps';
%print('-depsc', filename);
exportfig(gcf, filename, 'color', 'cmyk', 'width', scale*3.5 * 0.5, 'height', scale*3.5);
system(sprintf('epstopdf %s', filename));

% variance
figure(3);
clf;
X = textread('reweighting-stddev.txt'); X = X(:,2:end); Y = textread('singletemp-stddev.txt'); Y = Y(:,2:end);
max_variance = max(max([Y.^2 X.^2]));
clim = [0 max_variance];
imagesc([tvec tvec], temperatures, [Y.^2 X.^2], clim); axis square; colorbar;
title('variance');
xlabel('single temp | reweighting');
set(gca,'XTick',[]); 
set(gca,'YTick',[]);
%ylabel('temperature / K');
colormap(flipud([linspace(1,1,64)' linspace(0,1,64)' linspace(0,1,64)']));

% Generate PDF.
filename = 'alanine-dipeptide-variance-large.eps';
%print('-depsc', filename);
exportfig(gcf, filename, 'color', 'cmyk', 'width', scale*7.5 * 0.5, 'height', scale*3.5);
system(sprintf('epstopdf %s', filename));

%% Autocorrelation at fixed tau as a function of temperature.

tau = 21; % time index of autocorrelation function

figure(4);
clf; 
hold on;
box on;

% Reweighted correlation functions.
data = textread('correlation-reweighted.txt');
reweighted_temperatures = data(:,1);
C_lt = data(:,2:end);
[L,T] = size(C_lt);
data = textread('correlation-reweighted-error.txt');
dC_lt = data(:,2:end);

h = fill([reweighted_temperatures' fliplr(reweighted_temperatures')], [C_lt(:,tau)'+2*dC_lt(:,tau)' fliplr(C_lt(:,tau)'-2*dC_lt(:,tau)')], 'k');
color = [0 0 1];
set(h, 'FaceColor', color * 0.2 + [1 1 1] * 0.8);
hreweight = plot(reweighted_temperatures, C_lt(:,tau), '-', 'LineWidth', 2);
set(hreweight, 'color', color);
  
xlabel('temperature / K');
ylabel('C(\tau;\beta) for \tau = 10 ps');
axis([300 600 0 1])

% Individual temperature estimates.
data = textread('correlation.txt');
temperatures = data(:,1);
C_kt = data(:,2:end);
data = textread('correlation-error.txt');
dC_kt = data(:,2:end);
dC_kt(isnan(dC_kt)) = 0.0;
[K,T] = size(C_kt);

color = 'k';
%plot(tvec, C_kt(k,:), color);
hsingle = errorbar(temperatures, C_kt(:,tau), 2*dC_kt(:,tau), '.', 'MarkerSize', 10); % 95% confidence interval
set(hsingle, 'color', color);
legend([hsingle, hreweight], 'single temp', 'reweighting');

% Generate PDF.
filename = 'alanine-dipeptide-temperature-dependence.eps';
%print('-depsc', filename);
exportfig(gcf, filename, 'color', 'cmyk', 'width', scale*3.5, 'height', scale*2.5);
system(sprintf('epstopdf %s', filename));

%% 2D PMF

figure(5);
clf;
box on;
data = textread('2d-pmf.txt');
phi_bins = data(1,:);
psi_bins = phi_bins;
pmf = data(2:end,:);
data = textread('2d-pmf-error.txt');
dpmf = data(2:end,:);
clim = [0 10];
imagesc(phi_bins, psi_bins, pmf', clim);
xlabel('\phi / degrees');
ylabel('\psi / degrees');
axis square;
colorbar;
cmap = hot;
%cmap = [cmap; 1 1 1]; % add white as top color
colormap(cmap);
set(gca,'YDir','normal')

% Generate PDF.
filename = 'alanine-dipeptide-2d-pmf.eps';
%print('-depsc', filename);
exportfig(gcf, filename, 'color', 'cmyk', 'width', scale*3.5, 'height', scale*2.5);
system(sprintf('epstopdf %s', filename));

%% Contributing temperature weights.

figure(6);
clf;
box on;
data = textread('contributing-temperatures.txt');
temperatures = data(1,:);
ln_weights = data(2:end,:);
log10_weights = ln_weights * log10(exp(1)); % log base 10
clim = [-10 0];
imagesc(temperatures, temperatures, log10_weights, clim);
xlabel('contributing temperature / K');
ylabel('target temperature / K');
axis square;
colorbar;
colormap(flipud(hot));

% Generate PDF.
filename = 'alanine-dipeptide-log-contributions.eps';
%print('-depsc', filename);
exportfig(gcf, filename, 'color', 'cmyk', 'width', scale*3.5, 'height', scale*2.5);
system(sprintf('epstopdf %s', filename));

