% Perform analysis for p5ab RNA hairpin data and save it to a mat file.

clear;

% Load p5ab hairpin data at 50 kHz.
load fiber1.mat;

% Parameters
x_ddag = 12.57; % good dividing surface guess (determined by manual inspection)
%x_ddag = 13.0; % bad dividing surface guess (for debugging)
tmax = 60; % maximum number of samples to evaluate reactive flux and implied rates

% Select observed data from data.
trace_index = 1; % trace index to use
x_t_50kHz = trace{trace_index};
dt_50kHz = 1 / 50000; % s - sampling resolution of this trace

% Filter data to 1 kHz.
nfilter = 50;
x_t_1kHz = x_t_50kHz(nfilter:nfilter:end);
dt_1kHz = dt_50kHz * nfilter;

% Compute rate estimates as a function of dividing surface position.
disp('Computing rate estimates as a function of dividing surface position...');

%x_ddag_min = 11.5;
%x_ddag_max = 13.5;
x_ddag_min = x_ddag - 0.5; 
x_ddag_max = x_ddag + 0.5;
ndividing = 51;
tmax = 60;
x_ddag_points = linspace(x_ddag_min, x_ddag_max, ndividing);

tmax_1kHz = 60; % 60 ms
tmax_50kHz = 60*50; % 60 ms

PA_i_1kHz = zeros(ndividing,1);
PB_i_1kHz = zeros(ndividing,1);
kTS_i_1kHz = zeros(ndividing,1); % transition state theory estimate from crossing frequency
kRF_it_1kHz = zeros(ndividing,tmax_1kHz); % reactive flux rate estimates
kim_it_1kHz = zeros(ndividing,tmax_1kHz); % implied rates
kint_it_1kHz = zeros(ndividing,tmax_1kHz); % correlation function integral

PA_i_50kHz = zeros(ndividing,1);
PB_i_50kHz = zeros(ndividing,1);
kTS_i_50kHz = zeros(ndividing,1); % transition state theory estimate from crossing frequency
kRF_it_50kHz = zeros(ndividing,tmax_50kHz); % reactive flux rate estimates
kim_it_50kHz = zeros(ndividing,tmax_50kHz); % implied rates
kint_it_50kHz = zeros(ndividing,tmax_50kHz); % correlation function integral

for i = 1:ndividing
  i

  x_ddag = x_ddag_points(i);

  % 1 kHz data
  x_t = x_t_1kHz;
  dt = dt_1kHz;
  h_A = double(x_t <= x_ddag);
  h_B = double(x_t > x_ddag);
  PA_i_1kHz(i) = mean(h_A);
  PB_i_1kHz(i) = mean(h_B);
  kTS_i_1kHz(i) = compute_crossing_rate(x_t, x_ddag, dt);
  kRF_it_1kHz(i,:) = compute_reactive_flux(x_t, x_ddag, dt, tmax_1kHz);
  kim_it_1kHz(i,:) = compute_implied_rate(x_t, x_ddag, dt, tmax_1kHz);
  kint_it_1kHz(i,:) = compute_correlation_integral_rate(x_t, x_ddag, dt, tmax_1kHz);

  % 50 kHz data
  x_t = x_t_50kHz;
  dt = dt_50kHz;
  h_A = double(x_t <= x_ddag);
  h_B = double(x_t > x_ddag);
  PA_i_50kHz(i) = mean(h_A);
  PB_i_50kHz(i) = mean(h_B);
  kTS_i_50kHz(i) = compute_crossing_rate(x_t, x_ddag, dt);
  kRF_it_50kHz(i,:) = compute_reactive_flux(x_t, x_ddag, dt, tmax_50kHz);
  kim_it_50kHz(i,:) = compute_implied_rate(x_t, x_ddag, dt, tmax_50kHz);
  kint_it_50kHz(i,:) = compute_correlation_integral_rate(x_t, x_ddag, dt, tmax_50kHz);

  % Save analyzed data.
  save p5ab-analysis.mat
end

