% Gillespie simulation for a birth-death process
% reaction 1 death S -> I at rate alpha
% reaction 2 birth S -> 2S at rate beta
function birth_death_stochastic
% initialize the states
global Nde  alpha beta% N = total number of particles

cmykflg=1;
%if cmykflg=1, this code will make cmyk figure files
%cmykflg = 1;

N = 20; %initial population size
Nde = 30;  %this is the number of p_k's to simulate, needed because of birth
K = 10000;  % number of trials
alpha = 1;  % without loss of generality
beta = 0.5;  % must be smaller than alpha to guarantee eventual extinction

%set the rate constants for two reactions:
c(1) = alpha;  %S -> I  % I tracks the number of deaths
c(2) = beta; %S -> 2S

%Specify the change matrix

C = [-1,1;1,0];  %What happens to s and i when a reaction occurs:

%initialize the state space

s = N* ones(K,1);
i = zeros(K,1); %start each population with no deaths yet
  
S = s ; %This will keep track of the   trajectories
I = i ;
T = zeros(K,1);  %This will track the transition times
j = 1;
rk = [];

ndx = find(s>0);  %update only the s's that are not zero
Nt = length(ndx); %=K at first
Kt = 5000; % max number of time steps to prevent very long runs without termination
while (Nt>0&j<Kt)
    j = j+1;
    h(:,1) = c(1)*s;  %reaction rate 1
    h(:,2) = c(2)*s; % reaction rate 2
    hc = cumsum(h')'; % the cumulative sum of h
    H = sum(h')';
   
    rn = rand(Nt,2);
    
    T(:,j) = T(:,j-1);  % add the current time to T
    
    for k = 1:Nt
     % update only those for which s >0
     
    T(ndx(k),j) = - log(rn(k,1))/H(ndx(k)) +T(ndx(k),j-1); % time of next reaction
    rk  = min(find(rn(k,2) <=hc(ndx(k),:)/H(ndx(k)) )); % this determines which reaction occurs
    s(ndx(k)) = s(ndx(k)) + C(rk,1); % update s and i
    i(ndx(k)) = i(ndx(k)) + C(rk,2);
    
    end 
    % save the values of the first trajectory
    S(:,j) = s ;
    I(:,j) = i ;
    ndx = find(s>0);  %check to see which trajectories are not extinct yet
    if (isempty(ndx)==1)
        Nt = 0;
    else
        Nt = length(ndx); %Nt is the number of trajectories not yet extinct
    
    
    end
end
nj = j; 
  kj = 1;; % use the first trial to plot a sample trajectory 
 tt = reshape([T(kj,:);T(kj,:)],2*nj,1);  % make a plot with stepwise decay
 nn = reshape([S(kj, :);S(kj, :)],  2*nj,1); 
 %plot one decay trajectory
 figure(1)
 plot(tt(2:end),nn(1:end-1,1),tt,N*exp(-(alpha-beta)*tt),'--','linewidth',2)
 xlabel('t', 'fontsize', 20)
 ylabel('Particle Number', 'fontsize', 20)
 set(gca, 'YTick', 0:4:20);
 
 if(cmykflg==1)
print('../../figs_c/chapt_1/b_d_traj','-deps','-cmyk')
 end

% create the pdf for extinction times from the data
 
 
 [NN,TT]=hist(T(:,end),50);  % histogram of the extinction times with 50 boxes
 tt = sort(TT);
   
 dt = mean(TT(2:end)-TT(1:end-1));  %timestep increment on the histogram
 
 % NN/(K*dt) is the approximate pdf for the extinction times from data
 % it is normalized to have total integral = 1
 
 

% simulate the de's
p0 = zeros(Nde+1,1); %initialize p_k 
p0(N+1) = 1;
 
 tstep = 0.1;
t_end =  max(T(:,end))
tspan = [0:tstep:t_end];   
[T,S] = ode23s(@deRHS,tspan,p0); 
 
figure(3)
for j = 1:Nde
plot(T,S(:,j))
hold on
end
hold off
 
figure(2)
plot(TT,NN/(K*dt),'*' ,'linewidth',4)

 hold on
  plot( T,alpha*S(:,2),'--','linewidth',2)
 hold off
  xlabel('t','fontsize',20)
ylabel('dp_0/dt','fontsize',20)
 if(cmykflg==1)
print('../../figs_c/chapt_1/b_d_ext_hist','-deps','-cmyk')
 end

 
 % create a histogram for number of births
 Nb = max(I(:,end)-N)
 HNb = zeros(Nb+1,1);
 for j = 1:K
     jdx = I(j,end)-N+1;
     HNb(jdx) = HNb(jdx)+1;
 end
 xnb = [0:Nb];
 figure(5)
 plot(xnb,HNb,'*')
 xlabel('Number of Births','fontsize',20)
 ylabel('Number of Trials','fontsize',20)
 
function p_prime=deRHS(t,p)  % right hand side for ode system
global Nde alpha beta
 
K = [0:Nde]'; 
 
p_prime = -(alpha*K+beta*K).*p+alpha*[K(2:end);0].*[p(2:end);0] ...
    +beta*[0;K(1:end-1)].*[0;p(1:end-1)];


