% Generate p5ab hairpin figures for manuscript.

clear;

% Load p5ab RNA hairpin analysis.
load p5ab-analysis.mat;

%% PLOT FORCE TRAJECTORY AND HISTOGRAM
disp('Plotting force trajectory and histogram...');

figure(1)
clf;

% Plot force trajectory.
% Determine trajectory length.

x_ddag = 12.57; % good dividing surface guess (determined by manual inspection)
x_t = x_t_1kHz;
dt = dt_1kHz;
T = length(x_t); 

% Construct vector of data observation times.

tpoints = (1:T) * dt;
subplot('position', [0.1 0.2 0.7 0.7]);
plot(tpoints, x_t, 'k.');
hold on;
plot([0 T]*dt, x_ddag*[1 1], 'r-');
axis([0 T*dt min(x_t) max(x_t)]);
%set(gca, 'YTick', []);
set(gca, 'YTick', [11 12 13 14]);
ylabel('x / pN');
xlabel('t / s');

% Plot force histogram to right of trajectory.
subplot('position', [0.8 0.2 0.15 0.7]);
nbins = 100;
[Nhist, xhist] = hist(x_t, nbins);
%barh(xhist, Nhist);
plot(Nhist, xhist, 'k-', 'markersize', 1);
hold on;
fill([Nhist ones(1,nbins)], [xhist fliplr(xhist)], [0.5 0.5 0.5]);
axis([0 max(Nhist) min(x_t) max(x_t)]);
plot([0 max(Nhist)], x_ddag * [1 1], 'r-');
set(gca, 'XTick', []);
set(gca, 'YTick', []);
xlabel('P(x)');

% Export figure
exportfig(gcf, 'p5ab-1kHz-trajectory.eps', 'width', 3.4, 'height', 1.5, 'color', 'cmyk', 'fontmode', 'fixed', 'fontsize', 6, 'markermode', 'fixed', 'markersize', 1.0); 
system('epstopdf p5ab-1kHz-trajectory.eps');

%% PLOT FORCE HISTOGRAM, RATE ESTIMATES FOR TST AND RF, AND EQUILIBRIUM PROBABILITIES AT 1 KHZ AND 50 KHZ

figure(2);
clf;

% Plot force histogram.
%subplot(3,1,1);
haxes = axes('position', [0.18 0.82 0.80 0.18]);
plot(xhist, Nhist, 'k-');
hold on;
fill([xhist fliplr(xhist)], [Nhist zeros(1,nbins)], [0.5 0.5 0.5]);
xlabel('x / pN');
ylabel('P(x)');
set(gca, 'YTick', []);
axis([x_ddag_min x_ddag_max 0 max(Nhist)]);
box on;

% Plot rate estimate as a function of dividing surface position.
%subplot(3,1,2);
haxes = axes('position', [0.18 0.55 0.80 0.18]);
hold on;
%plot(x_ddag_points, kim_it_1kHz(:,end), 'k--'); 
%plot(x_ddag_points, kim_it_50kHz(:,end), 'k-');
plot(x_ddag_points, kim_it_1kHz(:,60), 'k--'); 
plot(x_ddag_points, kim_it_50kHz(:,60*50), 'k-');
plot(x_ddag_points, kTS_i_1kHz, 'r--'); 
plot(x_ddag_points, kTS_i_50kHz, 'r-'); 
plot(x_ddag_points, kRF_it_1kHz(:,3), 'g--'); 
plot(x_ddag_points, kRF_it_50kHz(:,150), 'g-'); 
set(gca, 'XTick', []);
ylabel('k / s^{-1}')
axis([x_ddag_min x_ddag_max 0 max(kTS_i_50kHz)]);
legend('1 kHz', '50 kHz');
%text(12.45, 250, '1 kHz');
box on;

% Plot rate estimate as a function of dividing surface position.
haxes = axes('position', [0.18 0.35 0.80 0.18]);
%subplot(3,1,2);
hold on;
%plot(x_ddag_points, kim_it_1kHz(:,end), 'k--'); 
%plot(x_ddag_points, kim_it_50kHz(:,end), 'k-');
plot(x_ddag_points, kim_it_1kHz(:,60), 'k--'); 
plot(x_ddag_points, kim_it_50kHz(:,60*50), 'k-');
plot(x_ddag_points, kTS_i_1kHz, 'r--'); 
plot(x_ddag_points, kTS_i_50kHz, 'r-'); 
plot(x_ddag_points, kRF_it_1kHz(:,3), 'g--'); 
plot(x_ddag_points, kRF_it_50kHz(:,150), 'g-'); 
set(gca, 'XTick', []);
ylabel('k / s^{-1}')
axis([x_ddag_min x_ddag_max 0 200]);
%text(12.45, 5000, '50 kHz');
box on;

% Plot pi1 and pi2 as a funtion of dividing surface position.
%subplot(3,1,3);
haxes = axes('position', [0.18 0.15 0.80 0.18]);
hold on;
plot(x_ddag_points, PA_i_1kHz, 'm-');
plot(x_ddag_points, PB_i_1kHz, 'c-');
%plot(x_ddag_points, PA_i_50kHz, 'r--');
%plot(x_ddag_points, PB_i_50kHz, 'r-');
axis([x_ddag_min x_ddag_max 0 1]);
xlabel('x^+ / pN');
ylabel('\pi')
legend('\pi_A', '\pi_B');
box on;

% Export figure
exportfig(gcf, 'p5ab-dividing.eps', 'width', 3.4, 'height', 3.8, 'color', 'cmyk', 'fontmode', 'fixed', 'fontsize', 6, 'markermode', 'fixed', 'markersize', 3); 
system('epstopdf p5ab-dividing.eps');

%% COMPARE RF CORRELATION FUNCTION WITH IMPLIED RATES AT 1 KHZ AND 50 KHZ

figure(3);
clf;

% 1 kHz
subplot(1,2,1);
%hmain = axes('position', [0.20 0.15 0.80 0.80]);
tval = dt_1kHz*[1:tmax_1kHz]*1000; % ms
PA = PA_i_1kHz(26);
plot(tval, kRF_it_1kHz(26,:), 'k.'); 
hold on;
plot(tval, kim_it_1kHz(26,:), 'r.');
xlabel('t / ms');
ylabel('k / s^{-1}')
set(gca, 'XTick', [0 10 20 30 40 50 60]);
set(gca, 'YTick', [0 25 50 75 100 125 150]);
legend('k_{RF}(t)', 'k_{im}(t)');
legend boxoff;
axis square;
axis([0 60 0 max(kRF_it_1kHz(26,:))]);
title('1 kHz');

% 50 kHz
subplot(1,2,2);
%hmain = axes('position', [0.20 0.15 0.80 0.80]);
tval = dt_50kHz*[1:tmax_50kHz]*1000; % ms
plot(tval, kRF_it_50kHz(26,:), 'k.');
hold on;
plot(tval, kim_it_50kHz(26,:), 'r.');
xlabel('t / ms');
ylabel('k / s^{-1}')
set(gca, 'XTick', [0 1 2 3 4 5 6]);
set(gca, 'YTick', [0 500 1000 1500 2000]);
%legend('k_{RF}', 'k_{im}');
%legend boxoff;
axis square;
axis([0 6 0 1500]);
title('50 kHz');

% 50 kHz inset
%hmain = axes('position', [0.70 0.48 0.18 0.18]);
hmain = axes('position', [0.64 0.48 0.10 0.18]);
tval = dt_50kHz*[1:tmax_50kHz]*1000; % ms
plot(tval, kRF_it_50kHz(26,:), 'k.');
hold on;
plot(tval, kim_it_50kHz(26,:), 'r.');
%xlabel('t / ms');
%ylabel('k(t) / s^{-1}')
set(gca, 'XTick', [3 3.5 4]);
set(gca, 'YTick', [30 35 40]);
box off
axis square;
axis([3 4 30 40]);

%hmain = axes('position', [0.70 0.48 0.18 0.18]);
hmain = axes('position', [0.79 0.48 0.10 0.18]);
tval = dt_50kHz*[1:tmax_50kHz]*1000; % ms
plot(tval, kRF_it_50kHz(26,:), 'k.');
hold on;
plot(tval, kim_it_50kHz(26,:), 'r.');
%xlabel('t / ms');
%ylabel('k(t) / s^{-1}')
set(gca, 'XTick', [59 59.5 60]);
set(gca, 'YTick', [30 35 40]);
box off
axis square;
axis([59 60 30 40]);

% Export figure
filename = 'p5ab-1kHz-50kHz.eps';
scale = 1.0;
exportfig(gcf, filename, 'width', scale*3.65, 'height', scale*3.65, 'color', 'cmyk', 'fontmode', 'fixed', 'fontsize', 6, 'markermode', 'fixed', 'markersize', 3); 
system(sprintf('epstopdf %s', filename));

% DEBUG
return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% OLD STUFF


figure(2);
clf;
plot([0:tmax]*dt, Ct, 'k.');
axis([0 tmax*dt -0.1 1]);
title('correlation function');
ylabel('<\delta h_A(0) \delta h_A(t)> / <\delta h_A^2>');
xlabel('t / s');

figure(3);
clf;
tmaxlife = length(lifetime_counts);
N = sum(lifetime_counts);
%lifetime_true = N * poisspdf(1:tmaxlife, Tij(1,2));
plot([1:tmaxlife] * dt, log(lifetime_counts), 'k.');
set(gca,'YTick', []);
xlabel('t / s');
ylabel('log P(t)'); 
title('lifetimes');

% Compute normalized time-correlation function for dh_A.
Cdt = zeros(1,tmax);
Cdt_true = zeros(1,tmax);
PA = mean(h_A);
dh_A = h_A - PA;
for t = 1:tmax
  Cdt(t) = mean(dh_A(1:T-t) .* dh_A(1+t:T)) / mean(dh_A.*dh_A);
end
Cdt = [1 Cdt];

% Compute reactive flux correlation function
kRF = - (Ct(2:tmax) - Ct(1:tmax-1)) / dt;

% Compute implied rate constant.
tval = [0:tmax] * dt;
kim = - tval.^(-1) .* log(Ct);






%% ESTIMATE STATIONARY PROBABILITIES

figure(5);
%hmain = axes('position', [0.20 0.15 0.80 0.80]);

% Quantize using dividing surface.
h1 = (x_t <= x_ddag);
h2 = (x_t > x_ddag);

% Compute correlation functions between sets.
c11_t = zeros(1,tmax+1);
c12_t = zeros(1,tmax+1);
c22_t = zeros(1,tmax+1);
for t = 0:tmax
  c11_t(t+1) = mean(h1(1:T-t) .* h1(1+t:T));
  c12_t(t+1) = 0.5 * mean(h1(1:T-t) .* h2(1+t:T) + h2(1:T-t) .* h1(1+t:T));
  c22_t(t+1) = mean(h2(1:T-t) .* h2(1+t:T));
end
c21_t = c12_t;
pi1 = mean(h1);
pi2 = mean(h2);

% Fit correlation function.
k = kim(end); % rate estimate
t = dt*[0:tmax]; % time vector
expmodel = @(a,t) a*exp(-k*t) + pi1^2;
a11 = nlinfit(t, c11_t, expmodel, (pi1 - pi1^2));

expmodel = @(a,t) a*exp(-k*t) + pi1*pi2;
a12 = nlinfit(t, c12_t, expmodel, -pi1*pi2);

expmodel = @(a,t) a*exp(-k*t) + pi2^2;
a22 = nlinfit(t, c22_t, expmodel, (pi2 - pi2^2));

% Plot correlation functions.
clf;
% c11
subplot(2,2,1); 
hold on;
plot(dt*[0:tmax]*1000, c11_t, 'k.');
plot(dt*[0 tmax]*1000, pi1*[1 1], 'b-');
plot(dt*[0 tmax]*1000, pi1^2*[1 1], 'g-');
plot(t*1000, a11*exp(-k*t) + pi1^2, 'r-');
plot(t*1000, (pi1-pi1^2)*exp(-k*t) + pi1^2, 'm-');
xlabel('t / ms');
ylabel('c_{11}(t)')
set(gca, 'XTick', [0 10 20 30 40 50 60]);
legend('c_{11}(t)', '\pi_1', '\pi_1^2', 'a_{11} e^{-k t} + \pi_1^2');
legend boxoff;
axis([0 dt*tmax*1000 0 1]);
title(sprintf('a_{11} = %.5f (f_1 = %f)', a11, a11 / (pi1 - pi1^2)));

% c12
subplot(2,2,2); 
hold on;
plot(dt*[0:tmax]*1000, c12_t, 'k.');
plot(dt*[0 tmax]*1000, 0*[1 1], 'b-');
plot(dt*[0 tmax]*1000, pi1*pi2*[1 1], 'g-');
plot(t*1000, a12*exp(-k*t) + pi1*pi2, 'r-');
plot(t*1000, -pi1*pi2*exp(-k*t) + pi1*pi2, 'm-');
xlabel('t / ms');
ylabel('c_{12}(t)')
set(gca, 'XTick', [0 10 20 30 40 50 60]);
legend('c_{12}(t)', '0', '\pi_1 \pi_2', 'a_{12} e^{-k t} + \pi_1*\pi_2');
legend boxoff;
axis([0 dt*tmax*1000 0 1]);
title(sprintf('a_{12} = %.5f', a12));


% c22
subplot(2,2,4);
hold on;
plot(dt*[0:tmax]*1000, c22_t, 'k.');
plot(dt*[0 tmax]*1000, pi2*[1 1], 'b-');
plot(dt*[0 tmax]*1000, pi2^2*[1 1], 'g-');
plot(t*1000, a22*exp(-k*t) + pi2^2, 'r-');
plot(t*1000, (pi2-pi2^2)*exp(-k*t) + pi2^2, 'm-');
xlabel('t / ms');
ylabel('c_{22}(t)')
set(gca, 'XTick', [0 10 20 30 40 50 60]);
legend('c_{22}(t)', '\pi_2', '\pi_2^2', 'a_{22} e^{-k t} + \pi_2^2');
legend boxoff;
axis([0 dt*tmax*1000 0 1]);
title(sprintf('a_{22} = %.5f (f_2 = %f)', a22, a22 / (pi2 - pi2^2)));

% Export figure
exportfig(gcf, 'p5ab-1kHz-stationary.eps', 'width', 5, 'height', 3.75, 'color', 'cmyk', 'fontmode', 'fixed', 'fontsize', 6, 'markermode', 'fixed', 'markersize', 3); 
system('epstopdf p5ab-1kHz-stationary.eps');

%% ESTIMATE STATIONARY PROBABILITIES USING SINGLE-TIME ESTIMATE OF c11

figure(6);
%hmain = axes('position', [0.20 0.15 0.80 0.80]);

% Quantize using dividing surface.
h1 = (x_t <= x_ddag);
h2 = (x_t > x_ddag);

% Compute correlation functions between sets.
c11_t = zeros(1,tmax+1);
c12_t = zeros(1,tmax+1);
c22_t = zeros(1,tmax+1);
for t = 0:tmax
  c11_t(t+1) = mean(h1(1:T-t) .* h1(1+t:T));
  c12_t(t+1) = 0.5 * mean(h1(1:T-t) .* h2(1+t:T) + h2(1:T-t) .* h1(1+t:T));
  c22_t(t+1) = mean(h2(1:T-t) .* h2(1+t:T));
end
c21_t = c12_t;
pi1 = mean(h1);
pi2 = mean(h2);

% Fit correlation function.
t = dt*[0:tmax]; % time vector
a11_t = (c11_t - pi1*pi2) ./ exp(-k*t);
a12_t = (c12_t - pi1*pi2) ./ exp(-k*t);
a22_t = (c22_t - pi1*pi2) ./ exp(-k*t);

% Plot a11.
clf;
% a11
subplot(2,2,1); 
hold on;
plot(dt*[0:tmax]*1000, a11_t, 'k.');
xlabel('t / ms');
ylabel('a_{11}(t)')
set(gca, 'XTick', [0 10 20 30 40 50 60]);
title('a_{11}(t)');

% a12
subplot(2,2,2); 
hold on;
plot(dt*[0:tmax]*1000, a12_t, 'k.');
xlabel('t / ms');
ylabel('a_{12}(t)')
set(gca, 'XTick', [0 10 20 30 40 50 60]);
title('a_{12}(t)');

% a22
subplot(2,2,3); 
hold on;
plot(dt*[0:tmax]*1000, a22_t, 'k.');
xlabel('t / ms');
ylabel('a_{22}(t)')
set(gca, 'XTick', [0 10 20 30 40 50 60]);
title('a_{22}(t)');

% Export figure
exportfig(gcf, 'p5ab-1kHz-a11t.eps', 'width', 5, 'height', 3.75, 'color', 'cmyk', 'fontmode', 'fixed', 'fontsize', 6, 'markermode', 'fixed', 'markersize', 3); 
system('epstopdf p5ab-1kHz-a11t.eps');

%% ESTIMATE AMPLITUDE OF PRIMARY DECAY PROCESS AS A FUNCTION OF DIVIDING SURFACE POSITION

figure(7);
%hmain = axes('position', [0.20 0.15 0.80 0.80]);

x_ddag_min = 11.5;
x_ddag_max = 13.5;
ndividing = 50;
x_ddag_points = linspace(x_ddag_min, x_ddag_max, ndividing);

pi1_i = zeros(1,ndividing);
pi2_i = zeros(1,ndividing);
k_i = zeros(1,ndividing);
a11_i = zeros(1,ndividing);
a12_i = zeros(1,ndividing);
a22_i = zeros(1,ndividing);

for i = 1:ndividing
  i

  % Quantize using dividing surface.
  h1 = (x_t <= x_ddag_points(i));
  h2 = (x_t > x_ddag_points(i));
  
  % Compute correlation functions between sets.
  c11_t = zeros(1,tmax+1);
  c12_t = zeros(1,tmax+1);
  c22_t = zeros(1,tmax+1);
  for t = 0:tmax
    c11_t(t+1) = mean(h1(1:T-t) .* h1(1+t:T));
    c12_t(t+1) = 0.5 * mean(h1(1:T-t) .* h2(1+t:T) + h2(1:T-t) .* h1(1+t:T));
    c22_t(t+1) = mean(h2(1:T-t) .* h2(1+t:T));
  end
  c21_t = c12_t;
  pi1 = mean(h1);
  pi2 = mean(h2);
  
  % Compute normalized time-correlation function for dh_A.
  Ct = zeros(1,tmax);
  dh1 = h1 - pi1;
  for t = 1:tmax
    Ct(t) = mean(dh1(1:T-t) .* dh1(1+t:T)) / mean(dh1.*dh1);
  end
  Ct = [1 Ct];
  
  % Compute implied rate constant.
  tval = [0:tmax] * dt;
  kim = - tval.^(-1) .* log(Ct);

  % Fit correlation function.
  k = kim(end); % rate estimate
  t = dt*[0:tmax]; % time vector
  expmodel = @(a,t) a*exp(-k*t) + pi1^2;
  a11 = nlinfit(t, c11_t, expmodel, (pi1 - pi1^2));
  
  expmodel = @(a,t) a*exp(-k*t) + pi1*pi2;
  a12 = nlinfit(t, c12_t, expmodel, -pi1*pi2);
  
  expmodel = @(a,t) a*exp(-k*t) + pi2^2;
  a22 = nlinfit(t, c22_t, expmodel, (pi2 - pi2^2));

  % Store data.
  pi1_i(i) = pi1;
  pi2_i(i) = pi2;
  k_i(i) = k;
  kTS_i(i) = kim(2);
  a11_i(i) = a11;
  a12_i(i) = a12;
  a22_i(i) = a22;
end

figure(8);
clf;

% Plot histogram.
%subplot(4,1,1);
haxes = axes('position', [0.18 0.82 0.80 0.18]);
plot(xhist, Nhist, 'k-');
hold on;
fill([xhist fliplr(xhist)], [Nhist zeros(1,nbins)], [0.5 0.5 0.5]);
xlabel('x / pN');
ylabel('P(x)');
set(gca, 'YTick', []);
axis([x_ddag_min x_ddag_max 0 max(Nhist)]);
box on;

% Plot rate estimate as a function of dividing surface position.
%subplot(4,1,2); 
haxes = axes('position', [0.18 0.55 0.80 0.18]);
hold on;
plot(x_ddag_points, k_i, 'k-'); % reactive-flux estimate of rate
%plot(x_ddag_points, kTS_i, 'r-'); % transition-state-theory rate (***)
fill([x_ddag_points fliplr(x_ddag_points)], [max(k_i)*ones(1,length(x_ddag_points)) fliplr(k_i)], [0.5 0.5 0.5]);
%xlabel('$x^\ddag$ / pN');
set(gca, 'XTick', []);
ylabel('k_{im} / s^{-1}')
%axis([x_ddag_min x_ddag_max 0 max(k_i)]);
box on;

% Plot a11 as a function of dividing surface position.
%subplot(4,1,3); 
haxes = axes('position', [0.18 0.35 0.80 0.18]);
hold on;
plot(x_ddag_points, a11_i ./ (pi1_i .* pi2_i), 'k-');
fill([x_ddag_points fliplr(x_ddag_points)], [a11_i ./ (pi1_i .* pi2_i) zeros(1,length(x_ddag_points))], [0.5 0.5 0.5]);
axis([x_ddag_min x_ddag_max 0 1]);
%xlabel('$x^\ddag$ / pN');
set(gca, 'XTick', []);
ylabel('a_{11} / (\pi_1 \pi_2)')
box on;

% Plot pi1 and pi2 as a funtion of dividing surface position.
%subplot(4,1,4);
haxes = axes('position', [0.18 0.15 0.80 0.18]);
hold on;
plot(x_ddag_points, pi1_i, 'r-');
plot(x_ddag_points, pi2_i, 'b-');
axis([x_ddag_min x_ddag_max 0 1]);
xlabel('x_{TS} / pN');
ylabel('\pi')
%legend('\pi_1', '\pi_2');
box on;

% Export figure
exportfig(gcf, 'p5ab-1kHz-dividing.eps', 'width', 3.4, 'height', 3.8, 'color', 'cmyk', 'fontmode', 'fixed', 'fontsize', 6, 'markermode', 'fixed', 'markersize', 3); 
system('epstopdf p5ab-1kHz-dividing.eps');

%%%%%%%%%%%%%%%%%%%

figure(9);
clf;

% Plot histogram.
subplot(3,1,1);
%haxes = axes('position', [0.18 0.80 0.80 0.18]);
plot(xhist, Nhist, 'k-');
hold on;
fill([xhist fliplr(xhist)], [Nhist zeros(1,nbins)], [0.5 0.5 0.5]);
xlabel('x / pN');
ylabel('P(x)');
set(gca, 'YTick', []);
axis([x_ddag_min x_ddag_max 0 max(Nhist)]);
box on;

% Plot rate estimate as a function of dividing surface position.
subplot(3,1,2); 
%haxes = axes('position', [0.18 0.55 0.80 0.18]);
hold on;
plot(x_ddag_points, kTS_i, 'r-'); % transition-state-theory rate (***)
%fill([x_ddag_points fliplr(x_ddag_points)], [max(k_i)*ones(1,length(x_ddag_points)) fliplr(k_i)], [0.5 0.5 0.5]);
%xlabel('x^\ddag / pN');
set(gca, 'XTick', []);
ylabel('k_{TST} / s^{-1}')
%axis([x_ddag_min x_ddag_max 0 max(k_i)]);
box on;

% Plot rate estimate as a function of dividing surface position.
subplot(3,1,3); 
%haxes = axes('position', [0.18 0.55 0.80 0.18]);
hold on;
plot(x_ddag_points, k_i, 'k-'); % reactive-flux estimate of rate
%fill([x_ddag_points fliplr(x_ddag_points)], [max(k_i)*ones(1,length(x_ddag_points)) fliplr(k_i)], [0.5 0.5 0.5]);
%xlabel('x^\ddag / pN');
set(gca, 'XTick', []);
ylabel('k_{RF} / s^{-1}')
xlabel('x_{TS} / pN');
%axis([x_ddag_min x_ddag_max 0 max(k_i)]);
box on;



% Export figure
exportfig(gcf, 'p5ab-1kHz-dividing-2.eps', 'width', 3, 'height', 3, 'color', 'cmyk', 'fontmode', 'fixed', 'fontsize', 6, 'markermode', 'fixed', 'markersize', 3); 
system('epstopdf p5ab-1kHz-dividing-2.eps');





