%
%  Likelihood_path.m
%
%               Run Gillespie on the model
%
%                   A -> B -> C  
%
%               with rates k1 = theta, k2 = 1 so as to compute the
%               derivative of E C(tend) with respect to theta.
%
%  usage:	Likelihood_path(tend,initial,theta)
%
%               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.
%
%  example:	Likelihood_path(3,[300, 0, 0],0.5)
%

function [x,M] = Likelihood_path(tend,initial,theta)

%rate constants
k = [theta 1];

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

%reaction vectors
zeta = [-1  0; 
         1 -1;
         0  1];

%"M" stores the weights used to compute the derivatives.
M = 0;

%initialize. T represents time and x is the state of the system.
T = 0;
x = initial;

%intensity functions
lambda = zeros([2,1]);

%main loop
for i=1:maxi-1,
    
    %update intensity functions
    lambda(1) = k(1)*x(1);
    lambda(2) = k(2)*x(2);
    
    lambda0 = lambda(1) + lambda(2);
    %calculate when the next reaction will take place.
    r = rand;
    
    if lambda0 == 0
        break
    else
        delta = log(1/r)/lambda0;
    end
    
    %update time
    T = T + delta;
    
    %find which reaction fired.
    r = rand;
    
    if r < lambda(1)/lambda0;
        loc = 1;
    else
        loc = 2;
    end
    
    %update x
    x = x + zeta(:,loc);
    
    %update M
    if loc == 1
        M = M + (1/theta)*(1 - lambda(1)*delta);
    else
        M = M - (1/theta)*lambda(1)*delta;
    end
    
    
    if T + delta > tend
        break
    end
    
end

end %The program

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