%
%  NRM_LV_timetrack.m	Run NRM on a Lotka Volterra model.
%      
%  usage:	NRM_LV_timetrack(tend)
%
%               tend = max time for the program to run = 20
%
%  example:	NRM_LV_timetrack(20)
%

function [] = NRM_LV_timetrack(tend)

%maximum number of steps to run before quitting.
maxi = tend*10^5;

%rate constants.
k = [10, 0.01, 0.01, 0.01, 10];

%reaction vectors (columns);
zeta= [1  0  -1  -1  0;
       0  1   0   0 -1];      


%initialize.  Assume start with 300 animals each.
initial = [300,300]';

%T is real time.
T = zeros([maxi,1]);

%state of the system.
x = zeros([2,maxi]);

%initialize the state of the system
x(:,1) = initial;

%holder for intensities.
lambda = zeros([5,1]);


%Set Poisson times and next Poisson Jumps.  Also set t's.
Tk = zeros([5,1]);
Pk = zeros([5,1]);
t = zeros([5,1]);

%set the first exponential holding times of the Poisson processes.
for i = 1:5
    r = rand;
    Pk(i) = log(1/r);
end


%main loop
for i =1:maxi-1
    
    %set intensity functions.
    lambda(1) = k(1)*x(1,i);
    lambda(2) = k(2)*x(1,i)*x(2,i);
    lambda(3) = k(3)*x(1,i)*x(2,i);
    lambda(4) = k(4)*x(1,i);
    lambda(5) = k(5)*x(2,i);
    
    %set t's, the times until each of the reactions would go if left alone.
    % This is equivalent to finding how long each "exponential alarm clock"
    % would require.
    for ell = 1:5
        if lambda(ell) > 0
            t(ell) = (Pk(ell) - Tk(ell))/lambda(ell);
            if t(ell)<0
                t
            end
        else
            t(ell) = inf;
        end
    end
    
    %loc gives the index of the minimum and t(loc) is the time until the
    %first alarm clock would go off.
    loc = 1;
    for ind = 2:5
        if t(loc)>t(ind)
            loc = ind;
        end
    end
    
    %if the time takes us past the end time, we should quit.
    if T(i) + t(loc)> tend
        T(i+1) = tend;
        x(:,i+1) = x(:,i);
        break
    end
    
    %update the state of the system and the time.
    x(:,i+1) = x(:,i) + zeta(:,loc);
    T(i+1) = T(i) + t(loc);
    
    %update the internal times of the Poisson processes.
    for m = 1:5
        Tk(m) = Tk(m) + lambda(m)*t(loc);
    end
    
    %We used one exponential random variable (the one from the Poisson
    %process that fired).  We need to set the next firing time of that
    %Poisson process.
    r = rand;
    Pk(loc) = Pk(loc) + log(1/r);
    
end

count = i;

%plot the timecourse of the state of the system.

figure
plot(T(1:count),x(1,(1:count)),T(1:count),x(2,(1:count)))
     legend('Rabbits','Foxes')

     %     title('Full version')

%Now for the exact solution
% 
% Num = 1000000;
% Delta = tend/Num;
% 
% ex = zeros([2,Num+1]);
% ex(:,1) = initial;
% T = zeros([Num+1,1]);
% for i = 1:Num
%     lambda(1) = k(1)*ex(1,i);
%     lambda(2) = k(2)*ex(1,i)*ex(2,i);
%     lambda(3) = k(3)*ex(1,i)*ex(2,i);
%     lambda(4) = k(4)*ex(1,i);
%     lambda(5) = k(5)*ex(2,i);
%     
%     ex(:,i+1) = ex(:,i) + zeta*lambda*Delta;
%     T(i+1) = T(i) + Delta;
% end
% 
% figure
% plot(T,ex(1,:),T,ex(2,:))
% legend('Rabbits','Foxes')
end %The program