% Example illustrating how rates can be measured.

clear;
clf;

% Transition matrix (row-stochastic)
Tij = [0.9 0.1;
       0.1 0.9;]; 

Tij = [0.95 0.05;
       0.05 0.95;]; 

T = 5000; % trajectory length
T = 50000; % trajectory length
nstates = 2;
sigmas = [0.1 0.2 0.3 0.4]; % sigma values to use
nsigmas = length(sigmas);

% Simulate trajectory.
state = zeros(1,T,'int8');
state(1) = 1;
for t = 2:T
  state(t) = randsample(1:nstates, 1, true, Tij(state(t-1),:));
end  

%% Generate trajectories
tpoints = 1:T;
for sigma_index = 1:nsigmas
  sigma = sigmas(sigma_index);

  % Simulate observables.
  obs = double(state) + sigma * randn(1,T);

  % Quantize using dividing surface.
  HA = (obs < 1.5);
  
  % Compute lifetime distribution of HA.
  counts = double(HA);
  for t0 = 1:T-1
    if (counts(t0)>0) && (counts(t0+1)>0)
      counts(t0+1) = counts(t0+1) + counts(t0);
    end
  end
  lifetime_counts = hist(counts, 0:max(counts));
  lifetime_counts = lifetime_counts(2:end);
  
  % Compute normalized time-correlation function for HA.
  tmax = floor(3/(1-Tij(1,1)));
  Ct = zeros(1,tmax);
  Ct_true = zeros(1,tmax);
  PA = mean(HA);
  for t = 1:tmax
    Ct(t) = (mean(HA(1:T-t) & HA(1+t:T)) - PA^2) / (PA - PA^2);
    Tijt = Tij^t;
    Ct_true(t) = (0.5 * Tijt(1,1) - 0.5^2) / (0.5 - 0.5^2);
  end
  Ct = [1 Ct];
  Ct_true = [1 Ct_true];
  ylabel('C(t)');
  xlabel('t');
  
  % Plot
  %subplot('position', [0.1 0.4 0.7 0.25]);
  subplot(nsigmas,4,4*(sigma_index-1)+1);
  plot(tpoints, obs, 'k.');
  hold on;
  plot([0 T], 1.5*[1 1], 'r-');
  axis([0 T 0 3]);
  set(gca, 'YTick', []);
  ylabel('x(t)');

  %subplot('position', [0.8 0.4 0.1 0.25]);
  subplot(nsigmas,4,4*(sigma_index-1)+2);
  nbins = 40;
  [Nhist, xhist] = hist(obs, nbins);
  barh(xhist, Nhist);
  hold on;
  axis([0 max(Nhist) 0 3]);
  plot([0 max(Nhist)], [1.5 1.5], 'r-');
  set(gca, 'XTick', []);
  set(gca, 'YTick', []);
  
%  subplot('position', [0.1 0.08 0.3 0.25]);
%  plot(0:tmax, Ct_true, 'r-');
%  hold on
%  plot(0:tmax, Ct, 'k.');
%  axis([0 tmax -0.1 1]);
%  title('correlation function');
%  ylabel('C(t)');
  
  % Plot lifetime fit
  %subplot('position', [0.5 0.08 0.3 0.25]);
  subplot(nsigmas,4,4*(sigma_index-1)+3);
  tmaxlife = length(lifetime_counts);
  N = sum(lifetime_counts);
  %lifetime_true = N * poisspdf(1:tmaxlife, Tij(1,2));
  lifetime_true = Tij(1,1).^(1:tmaxlife) * Tij(1,2);
  lifetime_true = N * lifetime_true / sum(lifetime_true);
  plot(1:tmaxlife, log(lifetime_counts), 'k.');
  hold on
  plot(1:tmaxlife, log(lifetime_true), 'r-');
  set(gca,'YTick', []);
  xlabel('t');
  ylabel('log P(t)'); 
  title('lifetimes');

  % Compute normalized time-correlation function for HA.
  tmax = 30;
  Ct = zeros(1,tmax);
  Ct_true = zeros(1,tmax);
  PA = mean(HA);
  for t = 1:tmax
    Ct(t) = (mean(HA(1:T-t) & HA(1+t:T)) - PA^2) / (PA - PA^2);
  end
  Ct = [1 Ct];

  % Compute reactive flux correlation function
  kRF = - (Ct(2:end) - Ct(1:end-1)) / PA;
  
  % Compute implied rate constant.
  tval = [0:tmax];
  kim = - tval.^(-1) .* log(Ct);

  % Plot reactive flux and implied rate constants
  subplot(nsigmas,4,4*(sigma_index-1)+4);
  %  subplot('position', [0.1 0.08 0.3 0.25]);
  %title('correlation function');
  %ylabel('C(t)');
  plot(tval(2:end) - 0.5, kRF, 'k.');
  hold on;
  plot(tval, kim, 'r.');
  xlabel('t');
  ylabel('k(t)')
  %set(gca, 'XTick', [0 10 20 30 40 50 60]);
  %set(gca, 'YTick', [0 25 50 75 100 125 150]);
  legend('k_{rf}', 'k_{im}');
  legend boxoff;
  %axis square;
  %title('1 kHz');
  true_rate = -log(1-Tij(1,2)-Tij(2,1));
  plot([0 tmax], true_rate*[1 1], 'r-');

end



% Export figure
exportfig(gcf, 'lifetime-fitting.eps', 'width', 2.5, 'height', 3.5, 'color', 'cmyk'); 