% A program for exit calculations for MBI-NIMBioS-CAMBAM PROGRAM
% Euler's method is used for solving the SDEs
% nsamp sample paths are calculated 
% The proportion exiting up to time t are calculated
% Exits occur when either population size is less than 1
% y1 and y2 are the two populations (different species)
% b1,d1,b2,d2 are the per capita birth and death rates 
% y10 y20 are the initial population sizes
% Problem-dependent statements are marked with a %*** (9 statements)
% Accuracy generally increases as h decreases
% icase=1 corresponds to the deterministic problem and so nsamp=1 
clf
clear
diary outputsde.m
for icase=1:2
clear tt
clear yp1
clear yp2
nsamp=1000;                                         %***
tmax=100;                                           %*** 
nt=1000;                                            %*** 
y10=15;                                             %***
y20=15;                                             %*** 
if(icase==1) nsamp=1; end
h=tmax/nt;
hs=sqrt(h);
randn('state',2);  %initiates the random number generator
te1=zeros(nsamp,1);
te2=zeros(nsamp,1);
te3=zeros(nsamp,1);
jj1=0;
jj2=0;
jj3=0;
for jj=1:nsamp 
   y1=y10; 
   y2=y20;
   yp1(1)=y1;
   yp2(1)=y2;
   r=randn(nt+1,2);
   nchk1=0;
   nchk2=0;
   n=0;
   t=0;
   chk=0;
   tt(1)=0;
   while (chk==0)
   n=n+1;
   t=t+h;
   if(jj==nsamp) tt(n+1)=t; end
   b1=.84;                                          %***
   d1=.40+.01*y1+.022*y2;                           %***
   b2=.90;                                          %***
   d2=.75+.0067*y2+.005*y1;                         %***
   f1=b1*y1-d1*y1;
   f2=b2*y2-d2*y2;
   g1=sqrt(b1*y1+d1*y1);
   g2=sqrt(b2*y2+d2*y2); 
   if(icase==1) g1=0; end
   if(icase==1) g2=0; end
   y1=y1+h*f1+hs*g1*r(n,1);
   y2=y2+h*f2+hs*g2*r(n,2);
   if(jj==nsamp) yp1(n+1)=y1; end
   if(jj==nsamp) yp2(n+1)=y2; end
% This is Euler's approximation to the SDE
if (y1 < 1) 
      chk=1;
      jj1=jj1+1;
      te1(jj1)=t; 
end
if (y2 < 1)
      chk=1;
      jj2=jj2+1;
      te2(jj2)=t;
end 
if (t > tmax)
      chk=1;
      jj3=jj3+1;
      te3(jj3)=t;
      chk=1; 
end 
end % end of while (chk==0) loop
end % end of for jj=1:nsamp loop
tp=0; tp1=0; tp2=0; tp3=0;
if(jj1 ~= 0) tp1=sum(te1)/jj1; end
if(jj2 ~= 0) tp2=sum(te2)/jj2; end
if(jj3 ~= 0) tp3=sum(te3)/jj3; end   
if(jj1+jj2 ~= 0) tp=(sum(te1)+sum(te2))/(jj1+jj2); end
p1=jj1/nsamp;
p2=jj2/nsamp;
p3=jj3/nsamp;
% tp1 and tp2 are mean exit times for populations 1 and 2
% tp3 is mean time for sample paths not exiting
% tp is mean exit time for paths that exit
% p1, p2 are proportions exiting for populations 1 and 2
% p3 is proportion not exiting in time tmax
disp('     ')
if(icase==1) disp('    Deterministic Calculational Results'); end
if(icase==2) disp('    Stochastic Calculation Results'); end
disp('          icase       nsamp       h            tmax')
disp((sprintf(' %12.0f %12.0f %12.5f %12.2f',icase,nsamp,h,tmax)));
disp('        tp1        p1')       
disp((sprintf(' %12.6f %12.6f', tp1, p1)));
disp('        tp2        p2') 
disp((sprintf(' %12.6f %12.6f', tp2, p2)));
disp('        tp3        p3')
disp((sprintf(' %12.6f %12.6f', tp3, p3)));
disp('        tp         p1+p2')
disp((sprintf(' %12.6f %12.6f', tp, p1+p2)));
  subplot(2,2,2*icase-1)
  set(gca,'fontsize',18,'linewidth',1.5);
  plot(yp1,yp2,'k-')
  xlabel('Pop. y1')
  ylabel('Pop. y2')
  if(icase==1) title('Deterministic'); end
  %if(icase==2) title('Stochastic'); end
  hold on
  subplot(2,2,2*icase)
  set(gca,'fontsize',18,'linewidth',1.5);  
  plot(tt,yp1,'r-',tt,yp2,'k-')
  xlabel('Time t')   
  ylabel('Pops. 1 and 2')
  if(icase==1) title('Deterministic'); end
  %if(icase==2) title('Stochastic'); end
  hold on  
end % end of for icase=1:2 loop
hold off
diary off