% Gillespie simulation for the SIS process
% reaction 1 infection S + I -> 2I
% reaction 2 recovery I -> S
 
function stochastic_sis
global N alpha beta 

set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7);

 
 
    
N = 50; %initial susceptible population size
% initial infected population = 1;
 
K = 2000;  % number of trials
alpha = 1;  % without loss of generality
beta = 35;   

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

%Specify the change matrix

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

%initialize the state space

s = zeros(K,1);
i = N*ones(K,1); %start each population with Ninfected

  
S = s ; %This will keep track of the   trajectories
I = i ;
T = zeros(K,1);  %This will track the transition times
j = 1; % number of reaction steps
rk = [];

ndx = find(i>0);  %update only the s's that are not zero
Nt = length(ndx); %=K at first
Kt = 15000; % 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.*i;  %reaction rate 1
    h(:,2) = c(2)*i; % reaction rate 2
     
    hc = cumsum(h')'; % the cumulative sum of h
    H = sum(h')';
   
    rn = rand(Nt,2); %find random number for the number of not yet terminated trajectories
    
    T(:,j) = T(:,j-1);  % add the current time to T
    
    for k = 1:Nt
     % update only those for which i >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  trajectories
    S(:,j) = s ;
    I(:,j) = i ;
    ndx = find(i>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
 Nt  % this should be zero to make sure everything has gone extinct
 kj = 1; % use the first several trials to plot   sample trajectories
  
%  some trajectories
 figure(1)
  for k = 1:kj
 plot(T(k,:),I(k,:)  ,'linewidth',2)
 hold on
 end
 xlabel('t', 'fontsize', 20)
 ylabel('I', 'fontsize', 20)
title('Sample Trajectory','fontsize',20)
 hold off
 
% now process the data 
 
 % 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
 figure(2 )
 
 NNne0=find(NN>0);
 plot(TT(NNne0), NN(NNne0)/(K*dt),'*')
 xlabel('recovery time','fontsize',20)
 ylabel('pdf of recovery times','fontsize',20)
 hold on

% now find p_k(t)
% simulate the de's
p0 = zeros(N+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:N+1
plot(T,S(:,j))
hold on
end
xlabel('t','fontsize',20)
ylabel('p_n(t)','fontsize',20)

hold off
 
figure(2)
 plot( T,beta*S(:,2),'--','linewidth',2)
  xlabel('t','fontsize',20)
ylabel('dp_0/dt','fontsize',20)

 hold off
 

 
function p_prime=deRHS(t,p)  % right hand side for ode system
global N alpha beta  
 
K = [0:N]'; 
 
p_prime = -(alpha*K.*(N-K)+beta*K).*p...
+alpha*[0;K(1:end-1)].*(N-[0;K(1:end-1)]).*[0;p(1:end-1)]...
    +beta*[K(2:end);0] .*[p(2:end);0];
     

 
