%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Solution to CSE Your Homework Assignment Project 6    %
%  More Models of Infection: It's Epidemic               %
%  solutiona.m Dianne P. O'Leary   11/03                 %
%                                                        %
%  This program solves Problems 1, 2, and 3.             %
%                                                        %
% We model the spread of an infection through a          %
% population, assuming that a person once infected       %
% cannot be reinfected.                                  %
%                                                        %
% We model two parts of the population:                  %
%    I(t) = proportion of population that is infected    %
%           at time t                                    %
%    S(t) = proportion of population that, at time t,    %
%           has never been infected                      %
% These quantities satisfy 0 <= I(t) <= 1 for t >= 0     %
%                          0 <= S(t) <= 1 for t >= 0     %
% The third part, the proportion of the population that  %
% was once infected but has now recovered, can be        %
% derived from these two:                                %
%    R(t) = 1 - I(t) - S(t).                             %
%                                                        %
% The model is                                           %
%    dI/dt =   mytrans * I * S - I / myrecov             %
%    dS/dt = - mytrans * I * S                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global mytrans myrecov myinitinfec mythreshhold

% Plotting symbols
pl(1).s = 'b';
pl(2).s = 'r-.';
pl(3).s = 'k:';
pl(4).s = 'g--';
pl(5).s = 'm';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set the parameters for the model:                      %
%                                                        %
%   mytrans = transmission rate for the infection        %
%   myrecov = # days contagious                          %
%   myinitinfec = proportion of population initially     %
%                 infected                               %
%   mythreshhold = tolerance; when infected or           %
%                  susceptible drop below this value,    %
%                  we conclude that the epidemic ended.  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mytrans = .8;  
myrecov = 4;  
myinitinfec = .005;
mythreshhold = 1.e-5;

nt = 3000;      % maximum number of days in simulation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construct the initial conditions.                      %
% The vector y0 has three components:                    %
%     infected, susceptible, recovered.                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

y0 = epi1initial(0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Problem 1:                                             % 
% Solve the ODE model, stopping when the number of       %
% infected or the number of susceptible drops below a    %
% threshhold.                                            %
% Plot the results.                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

strng = sprintf('Solution from Ordinary Differential Equation Model')
options=odeset('Events',@epi1end);
tspan = 0:.25:nt;
[to,yo] = ode23s(@epi1,tspan,y0,options);

figure(1)
i = length(to);
plot(to(1:i),yo(1:i,1),pl(2).s,to(1:i),yo(1:i,2),pl(3).s, ...
     to(1:i),yo(1:i,3),pl(1).s,'Linewidth',1.5)
axis([0 to(i) 0 1])
legend('Infected','Susceptible','Recovered')
xlabel('time')
ylabel('proportion of population')
title(strng)
print -depsc prob1.eps
sprintf('Proportion of population infected: %f',yo(i,3))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Problem 2:                                             % 
% Solve the DAE model, stopping when the number of       %
% infected or the number of susceptible drops below a    %
% threshhold.                                            %
% Plot the results.                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

strng = sprintf('Solution from Differential Algebraic Equation Model')
options=odeset('Events',@epi1end,'Mass',[1 0 0;0 1 0;1 1 1]);
[ta,ya] = ode23s(@epi1dae,tspan,y0,options);

figure(2)
i = length(ta);
tfinal = ta(i);
plot(ta(1:i),ya(1:i,1),pl(2).s,ta(1:i),ya(1:i,2),pl(3).s, ...
     ta(1:i),ya(1:i,3),pl(1).s,'LineWidth',1.5)
axis([0 tfinal 0 1])
legend('Infected','Susceptible','Recovered')
xlabel('time')
ylabel('proportion of population')
title(strng)
print -depsc prob2.eps
sprintf('Proportion of population infected: %f',ya(i,3))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Problem 3:                                             % 
% Solve the DDE model, stopping when the number of       %
% infected or the number of susceptible drops below a    %
% threshhold.                                            %
% Plot the results.                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

strng = sprintf('Solution from Delay Differential Equation Model')
options=ddeset('Jumps',[0 myrecov],'Events',@epi1end);
sol = dde23(@epi1delay,[myrecov],@epi1initial,[0,nt],options);
ntt = sol.xe;
% Note that deval returns the solution transposed
td = 0:.25:ntt;
yd = deval(sol,td)';

figure(3)
i = length(td);
plot(td,yd(:,1),pl(2).s,td,yd(:,2),pl(3).s, ...
     td,yd(:,3),pl(1).s,'LineWidth',1.5)
axis([0 tfinal 0 1])
legend('Infected','Susceptible','Recovered')
xlabel('time')
ylabel('proportion of population')
title(strng)
print -depsc prob3.eps
sprintf('Proportion of population infected: %f',yd(i,3))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check conservation: the sum I+S+R should always be 1.  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:length(ta),
  yya = sum(ya(i,:));
end
for i=1:length(to),
  yyo = sum(yo(i,:));
end
for i=1:length(td),
  yyd = sum(yd(i,:));
end
disp('ODE conservation (should be zero):')
disp(max(abs(1-yyo)))
disp('DAE conservation (should be zero):')
disp(max(abs(1-yya)))
disp('DDE conservation (should be zero):')
disp(max(abs(1-yyd)))