% an agent based simulation of SIRS epidemic with OU process movement in one dimension 
clear

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

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

% set some parameters
L = 100; % length of 1D domain
dt = 0.01; %the time step  

Nstep = 3000;  % the number of time steps to take
D = 1;  % the diffusion coefficient for the particles
kappa  = 0; % OU S restoring force
 
scD = sqrt(2*D*dt);  % variance for diffusion
%there are three possible reactions
rxn = 3;  % reaction (transmission) rate I+S -> 2I
% set rxn = 0 to see pure diffusion
Rxn = 0.1; %size of the reaction neighborhood
rmv = 0.5; % recovery rate I -> R
rmv = 0; % for no recovery
% set rvm = 0 to see Fisher's equation in an agent based model
rs = 0.; % R -> S rate 
% set rs >0 to see S->I->R->S disease dynamics
 
%assign colors to the different particle types
col(1) = 'g';
col(2) = 'r';
col(3) = 'b';
 
    clear X
    clear S
   S0 = 0.75;  %This is the density of the initial susceptible population

% this means that there are a total of p*L particles 
N = fix(S0*L); % total number of particles
 
%For the OU process, distribute ``homes" regularly
%  assign home positions
xhome = [1:N]'*L/N;
%initial positions
x = xhome;
t = 0;
%there are three possible states
% s=1; susceptible
% s=2 ; infected
% s=3; recovered
s = ones(N,1); %initially all are susceptible

% except one infected at the center of the domain
 
s(fix(N/2)) = 2;
X = x; % save the positions
T (1) = t;
S = s; %Save the states;
for n = 1:Nstep

%there are two steps to take
%1) move the particles by OU process  

x = x -kappa .*(x-xhome)*dt+ scD*randn(N,1);
x = abs(x); % reflecting boundaries at x=0;
x = L - abs(x-L); % reflecting boundaries at x=L;

% calculate reaction probabilities

for j = 1:N
% reaction probability for susceptibles:    
    if(s(j)==1)
        % calculate distances between susceptible and infected particles
        dist = zeros(1,N);
    for k = 1:N
        if (s(k) ==2)
        dist(k) =  (x(j)-x(k)).^2;
        end
    end
    DltR = sum((dist < Rxn^2)&(dist>0)); % find number of particles within distance Rxn
    rxp(j) = rxn*DltR/(2*Rxn);
    end
    if (s(j)==2)
        rxp(j) = rmv;% the removal rate
    end
    if (s(j) ==3)
        rxp(j) = rs; % the recovery rate
    end
end
 
rxpp = find(rxp*dt>rand(1,N)); %These are the particles that change state
s(rxpp) = s(rxpp) + 1;
rxpp = find(s==4);
s(rxpp) = 1;
% one cycle of the process is complete
t = t + dt;
T(n+1) = t;
X=[X,x];
S = [S,s];

end

figure(1)
for k = 1:N
    tp1 = find(S(k,:)==1);
    tp2 = find(S(k,:) == 2);
    tp2m= min(tp2);
 plot(T(tp1),X(k,tp1), strcat(col(1) ));
 plot(T(tp2),X(k,tp2), strcat(col(2) ),T(tp2m),X(k,tp2m),'b*');
 hold on
end
axis([0 T(end) 0 L])
xlabel('time','fontsize',20)
ylabel('x','fontsize',20)
hold off

if(cmykflg==1) 
print('../../figs_c/chapt_6/oned_si_agnt','-deps','-cmyk')
end

nI = sum(S==2);


nIx = find(nI<N);
sl = sum(T(nIx).*nI(nIx))/sum(T(nIx).^2);
figure(2)
plot(T,nI/N,T(nIx),sl*T(nIx)/N,'--','linewidth',2)
xlabel('time','fontsize',20)
ylabel('% of Infected','fontsize',20)
axis([0 T(end) 0 1])

 

