% Generate plot for model system.

clear;
figure(4);
clf;

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

Tij = [0.95 0.05;
       0.05 0.95;]; 

true_rate = -log(1-Tij(1,2)-Tij(2,1)); % true rate

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

% Simulate discrete state sequence.
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);

  % Plot trajectory.
  subplot('position', [0.1 1.0 - 0.7*(sigma_index/nsigmas) 0.27 0.7/nsigmas]);
  %subplot(nsigmas,4,4*(sigma_index-1)+1);
  plot(tpoints(1:5000), obs(1:5000), 'k.');
  hold on;
  plot([0 5000], 1.5*[1 1], 'r-');
  axis([0 5000 0 3]);
  set(gca, 'YTick', [1 2 ]);
  set(gca, 'XTick', []);
  ylabel('x(t)');
  if (sigma_index == 4)
    xlabel('t');
    set(gca, 'XTick', [0 1000 2000 3000 4000 5000]);
  end

  % Plot histogram next to trajectory.
  subplot('position', [0.37 1.0 - 0.7*(sigma_index/nsigmas) 0.1 0.7/nsigmas]);
  %subplot(nsigmas,4,4*(sigma_index-1)+2);
  nbins = 51;
  [Nhist, xhist] = hist(obs, nbins);
  Nhist = Nhist / max(Nhist);
  %barh(xhist, Nhist);
  plot(Nhist, xhist, 'k-');
  hold on;
  fill([Nhist zeros(1,nbins)], [xhist fliplr(xhist)], [0.5 0.5 0.5]);
  axis([0 max(Nhist) 0 3]);
  plot([0 max(Nhist)], [1.5 1.5], 'r-');
  set(gca, 'XTick', []);
  set(gca, 'YTick', []);
  if (sigma_index == 4)
    xlabel('P(x)');
  end
  text(0.6, 2.6, sprintf('\\sigma = %.1f', sigma));

  % Plot crossing rate as a function of dividing surface.
  x_ddag = 1.5;
  x_ddag_min = x_ddag - 0.25;
  x_ddag_max = x_ddag + 0.25;
  ndividing = 51;
  tmax = 20;
  x_ddag_points = linspace(x_ddag_min, x_ddag_max, ndividing);

  PA_i = zeros(ndividing,1);
  PB_i = zeros(ndividing,1);
  kTS_i = zeros(ndividing,1); % transition state theory estimate from crossing frequency
  kRF_it = zeros(ndividing,tmax); % reactive flux rate estimates
  kim_it = zeros(ndividing,tmax); % implied rates
  kint_it = zeros(ndividing,tmax); % correlation function integral rate

  for i = 1:ndividing
    i

    x_ddag = x_ddag_points(i);
    x_t = obs;
    dt = 1;
    h_A = double(x_t <= x_ddag);
    h_B = double(x_t > x_ddag);
    PA_i(i) = mean(h_A);
    PB_i(i) = mean(h_B);
    kTS_i(i) = compute_crossing_rate(x_t, x_ddag, dt);
    kRF_it(i,:) = compute_reactive_flux(x_t, x_ddag, dt, tmax);
    kim_it(i,:) = compute_implied_rate(x_t, x_ddag, dt, tmax);
    kint_it(i,:) = compute_correlation_integral_rate(x_t, x_ddag, dt, tmax);
  end

  subplot('position', [0.52 1.0 - 0.7*(sigma_index/nsigmas) 0.2 0.7/nsigmas]);
  %subplot(nsigmas,4,4*(sigma_index-1)+3);
  hold on;
  plot([x_ddag_min x_ddag_max], true_rate*[1 1], 'b:', 'LineWidth', 2); % plot true rate
  plot(x_ddag_points, kim_it(:,end), 'k-'); 
  plot(x_ddag_points, kTS_i, 'r-'); 
  plot(x_ddag_points, kRF_it(:,2), 'g-'); 
  %plot(x_ddag_points, kint_it(:,end), 'c-', 'LineWidth', 2); 
  set(gca, 'XTick', []);
  ylabel('k / s^{-1}')
  ymax = max([max(kTS_i) max(kim_it(:,end)) max(kRF_it(:,2))]) + 0.2;
  ymax = 0.3;
  axis([x_ddag_min x_ddag_max 0 ymax]);
  box on;
  if (sigma_index == 4)
    xlabel('x^+');
    set(gca, 'XTick', [1.3 1.4 1.5 1.6 1.7]);
  end

  % Plot reactive flux and implied rate constants
  %subplot(nsigmas,4,4*(sigma_index-1)+4);
  subplot('position', [0.77 1.0 - 0.7*(sigma_index/nsigmas) 0.2 0.7/nsigmas]);
  %title('correlation function');
  hold on;
  plot([0 tmax], true_rate*[1 1], 'b:', 'LineWidth', 2); % plot true rate in red
  plot(1:tmax, kRF_it(26,:), 'g-');
  plot(1:tmax, kim_it(26,:), 'k-');
  plot([0 tmax], kTS_i(26) * [1 1], 'r-');
  %plot(1:tmax, kint_it(26,:), 'c-', 'LineWidth', 2); 
  set(gca, 'XTick', []);
  ylabel('k / s^{-1}')
  %set(gca, 'XTick', [0 10 20 30 40 50 60]);
  %set(gca, 'YTick', [0 25 50 75 100 125 150]);
  axis([0 tmax 0 ymax]);
  box on;
  if (sigma_index == 1)
    legend('k', 'k_{RF}', 'k_{im}', 'k_{crossing}');
  end
  if (sigma_index == 4)
    xlabel('t');
    set(gca, 'XTick', [0 5 10 15 20]);
  end
end

% Export figure
scale = 1.5;
exportfig(gcf, 'model-rate-estimation.eps', 'width', scale*7.5, 'height', scale*6.0, 'color', 'cmyk'); 
system('epstopdf model-rate-estimation.eps');