% A two-patch SIS program for MBI-NIMBioS-CAMBAM PROGRAM
% Euler's method is used for solving the SDEs
% nsamp sample paths are calculated 
% s1,s2,i1,i2 are population sizes in two patches
% gam1,gam2,di12,di21,ds12,ds21,bet1,bet2 are constants
% s10,s20,i10,i20 are the initial population sizes
% Problem-dependent statements are marked with a %*** (15 statements)
% Accuracy generally increases as h decreases
% icase=1 corresponds to the deterministic problem and so nsamp=1 
clf
clear
%diary outputgroupsde.m
for icase=1:2
clear tt
clear s1
clear s2
clear i1
clear i2
nsamp=1000;                                         %***
tmax=25;                                            %*** 
nt=500;                                             %***                    
s10=25;                                             %***
s20=20;                                             %***
i10=30;                                             %***
i20=25;                                             %*** 
gam1=.05;                                           %***
gam2=.4;                                            %***
di12=.3;                                            %***
di21=.05;                                           %***
ds12=0;                                             %***  
ds21=.25;                                           %***
bet1=.5;                                            %***
bet2=.01;                                           %***
if(icase==1) nsamp=1; end
h=tmax/nt;
hs=sqrt(h);
randn('state',3);  %initiates the random number generator
ss1=0;
ss12=0;
ss2=0;
ss22=0;
si1=0;
si12=0;
si2=0;
si22=0;
for jj=1:nsamp 
   s1=s10; 
   s2=s20;
   i1=i10;
   i2=i20;  
   s1p(1)=s1;
   s2p(1)=s2;
   i1p(1)=i1;
   i2p(1)=i2;
   r=randn(nt+1,8);
   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
   nt1=s1+i1;
   nt2=s2+i2;
   fs1=-s1*i1*bet1/nt1+gam1*i1+ds21*s2-ds12*s1;
   fi1=s1*i1*bet1/nt1-gam1*i1-di12*i1+di21*i2;
   fs2=-s2*i2*bet2/nt2+gam2*i2-ds21*s2+ds12*s1;
   fi2=s2*i2*bet2/nt2-gam2*i2-di21*i2+di12*i1;
   hlp1=-sqrt(s1*i1*bet1/nt1)*r(n,1)+sqrt(gam1*i1)*r(n,2);
   hlp2=sqrt(ds21*s2)*r(n,3)-sqrt(ds12*s1)*r(n,4);
   hlp3=sqrt(di21*i2)*r(n,5)-sqrt(di12*i1)*r(n,6);
   hlp4=-sqrt(s2*i2*bet2/nt2)*r(n,7)+sqrt(gam2*i2)*r(n,8);
   gs1= hlp1+hlp2;
   gi1= -hlp1+hlp3;
   gs2= hlp4-hlp2;
   gi2= -hlp4-hlp3;
   if(icase==1) gs1=0; end
   if(icase==1) gs2=0; end
   if(icase==1) gi1=0; end
   if(icase==1) gi2=0; end
   s1=s1+h*fs1+hs*gs1;
   i1=i1+h*fi1+hs*gi1;
   s2=s2+h*fs2+hs*gs2;
   i2=i2+h*fi2+hs*gi2;
% This is Euler's approximation to the SDE
   if(jj==nsamp) s1p(n+1)=s1; end
   if(jj==nsamp) i1p(n+1)=i1; end
   if(jj==nsamp) s2p(n+1)=s2; end
   if(jj==nsamp) i2p(n+1)=i2; end 
%if (s1 < 0) 
%      i1=i1+s1;
%      s1=0;
%end
%if (i1 < 0)
%      s1=s1+i1;
%      i1=0;
%end
%if (s2 < 0) 
%      i2=i2+s2;
%      s2=0;
%end
%if (i2 < 0)
%      s2=s2+i2;
%      i2=0;
%end 
if (t > tmax)
      ss1=s1+ss1;
      ss12=s1*s1+ss12;
      ss2=s2+ss2;
      ss22=s2*s2+ss22;
      si1=i1+si1;
      si12=i1*i1+si12;
      si2=i2+si2;
      si22=i2*i2+si22;
      chk=1;
end 
end % end of while (chk==0) loop
end % end of for jj=1:nsamp loop
      ss1=ss1/nsamp;
      ss12=ss12/nsamp-ss1*ss1;
      ss2=ss2/nsamp;
      ss22=ss22/nsamp-ss2*ss2;
      si1=si1/nsamp;
      si12=si12/nsamp-si1*si1;
      si2=si2/nsamp;
      si22=si22/nsamp-si2*si2;
      ntot=ss1+ss2+si1+si2;
disp('     ')
if(icase==1) disp('    ODE Calculational Results'); end
if(icase==2) disp('    SDE Calculational Results (Means and Standard Deviations)'); end
disp('          icase       nsamp       h            tmax')
disp((sprintf(' %12.0f %12.0f %12.5f %12.2f',icase,nsamp,h,tmax)));
disp('       S1(mean)      I1(mean)     S2(mean)     I2(mean)')       
disp((sprintf(' %12.6f %12.6f %12.6f %12.6f', ss1, si1, ss2, si2)));
disp('       S1(std)       I1(std)     S2(std)      I2(std)') 
disp((sprintf(' %12.6f %12.6f %12.6f %12.6f', sqrt(ss12), sqrt(si12), sqrt(ss22), sqrt(si22))));
disp('       ntot ') 
disp((sprintf(' %12.6f %12.6f %12.6f %12.6f', ntot)));
  subplot(2,2,2*icase-1)
  plot(tt,s1p,'k-',tt,i1p,'r--','linewidth',2)
  xlabel(' Time t')
  ylabel('S_1 and I_1')
  legend('S_1', 'I_1');
  axis([0,tmax,0,60])
  if(icase==1) title('ODE'); end
  if (icase==2) title('SDE'); end
  hold off
  subplot(2,2,2*icase) 
  plot(tt,s2p,'k-',tt,i2p,'r--','linewidth',2)
  xlabel('Time t')   
  ylabel('S_2 and I_2')
  legend('S_2', 'I_2')
  axis([0,tmax,0,60])
  if(icase==1) title('ODE'); end
  if(icase==2) title('SDE'); end
  hold off  
end % end of for icase=1:2 loop
hold off
%diary off