%
%  CFD_path.m
%
%               Run coupled finite differences on the model
%
%                   A -> B -> C  
%
%               with rates k1 = theta, k2 = 1.
%
%  usage:	CFD_path(tend,initial,theta,epsilon)
%
%               tend = max time for the program to run
%
%               initial = initial condition, ex: [30, 0, 0]
%
%               theta = value of parameter at which we are computing the
%                       derivative.
%
%               epsilon = epsilon used in centered finite difference
%
%  example:	CFD_path(3,[300, 0, 0] , 0.5, 0.5*1/20)
%



function [x_plus,x_minus] = CFD_path(tend,initial,theta,epsilon)


%rate constants
k = [theta 1];

%maximum number of steps to allow.
maxi = 10^4;

% reaction vectors (there are three for each of the 2 reactions now)
zeta = [-1 -1  0    0  0  0; 
		 1  1  0   -1 -1  0; 
		 0  0  0    1  1  0;
        -1  0 -1    0  0  0;
         1  0  1   -1  0 -1;
         0  0  0    1  0  1];

%initialize.  T represents time and x is the state of the enlarged system.
%x(1:3) gives the system with parameter theta + epsilon/2 whereas x(4:6)
%gives the system with parameter theta - epsilon/2.
T = 0;
x = zeros([6,1]);
x(1:3) = initial;
x(4:6) = initial;

%lambda1 are intensities for x(1:3).  lambda2 are intensities for x(4:6)
lambda1 = zeros([2,1]);
lambda2 = zeros([2,1]);

%As are intensities for actual representation.  The first 3
%components represnets the first reaction, the second 3 represent
% the second reaction.

A = zeros([6,1]);

%Set Poisson times and next Poisson Jumps.  One for each of the 6 actual
%reactions.  Also set t_ks.
Tk = zeros([6,1]);
Pk = zeros([6,1]);
t  = zeros([6,1]);

%set the first unit exponentials for the 6 Poisson processes.
for j = 1:6
    r = rand;
    Pk(j) = log(1/r);
end


for i=1:maxi-1,    
    
    %update intensity functions
    lambda1(1) = (k(1)+epsilon/2)*x(1);
    lambda1(2) = k(2)*x(2);
    
    lambda2(1) = (k(1)-epsilon/2)*x(4);
    lambda2(2) = k(2)*x(5);
 
    %Now get the intensity functions for the 6 actual channels.
    A(1) = min(lambda1(1),lambda2(1));
    A(2) = lambda1(1) - A(1);
    A(3) = lambda2(1) - A(1);
    
    A(4) = min(lambda1(2),lambda2(2));
    A(5) = lambda1(2) - A(4);
    A(6) = lambda2(2) - A(4);
    
    %find which reaction happens next.
    for j = 1:6
        if A(j) > 0
            t(j) = (Pk(j) - Tk(j))/A(j);
        else
            t(j) = inf;
        end
    end
    
    %find the time of the next reaction.
    loc = 1;
    for ell = 2:6
        if t(loc)>t(ell)
            loc = ell;
        end
    end
    
    %if all intensities are zero, we're done.  quit.
    if t(loc) == inf;
        break
    end
    
    %update time    
    T = T + t(loc);
    
    %update the state of the system.
    x = x + zeta(:,loc);
    
    %update the one Poisson process that fired.
    r = rand;
    Pk(loc) = Pk(loc) + log(1/r);

    %update local times.
    Tk = Tk + A*t(loc);
    
    if T >= tend
        break
    end
    
end 

x_plus = x(1:3);
x_minus = x(4:6);



%diff = (x(2) - x(4))/h;

%CPUTime = etime(clock,bt)


end %The program

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