%
%  CRN_main.m
%
%           Use common random numbers to compute the 
%           sensitivity of the expected number of C molecules with respect 
%           to theta for the model
%
%                   A -> B -> C  
%
%               with rates k1 = theta, k2 = 1.  
%
%  usage:	CRN_main(n,tend,theta,epsilon)
%
%               n = number of times to run CMC_path to get stats.
%
%               tend = terminal time
%
%               theta = value of parameter at which we are computing the
%                       derivative.
%
%               epsilon = epsilon used in centered finite difference
%
%  output:  m_C, m_d, v_d
%
%               m_C = estimated value of E C(tend).
%               m_d = estimated value of the derivative of E C(tend)
%               v_d = variance of estimate for m_d.
%
%  example:	CRN_main(5000,3,0.5,0.5*1/20)
%

function [m_C,m_d,v_d] = CRN_main(n,tend,theta,epsilon)

%get the necessary seeds for the random number generator.
rng('shuffle')
R = floor(rand([n,1])*sum(100*clock));

%SEED_CLOCK = sum(100*clock);
%rand('state',SEED_CLOCK)
%R = rand([n,1])*sum(100*clock);

%initial condition
initial = [300; 
            0; 
            0];
%"data_plus" will store the output with theta + epsilon/2 at time tend.  
%Each column will represent a realization.
data_plus = zeros([3,n]);

%"data_minus" will store the output with theta - epsilon/2 at time tend.  
%Each column will represent a realization.
data_minus = zeros([3,n]);

%"derivative" will give the estimated derivative achieved by taking a
%difference and dividing by epsilon.
derivative = zeros([n,1]);

for j = 1:n
    %rand('state',R(j));
    rng(R(j));
    [data_plus(:,j)] = CRN_path(tend,initial,theta + epsilon/2);
   
    %rand('state',R(j));
    rng(R(j));
    [data_minus(:,j)] = CRN_path(tend,initial,theta - epsilon/2);
    
    derivative(j) = (data_plus(3,j) - data_minus(3,j))/epsilon;
    
end


m_C = mean(data_plus(3,:) + data_minus(3,:))/2;
m_d = mean(derivative);
v_d = var(derivative)/n;


end %The program

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