%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                        %
%  A solution to CSE Your Homework Assignment Project 5  %
%  Models of Infection: Person to Person                 %
%  Dianne P. O'Leary   09/03                             %
%                                                        %
%  solution.m                                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  This project investigates a model for the spread of   %
%  infection in a hospital ward.                         %
%                                                        %
%  The parameters for the model:                         %
%                                                        %
%   n,m:     dimensions of the grid of hospital beds     %
%   k:       number of days of contagion                 %
%   tau:     transmission rate for the infection         %
%   delta:   mobility rate of the population             %
%   nu:      vaccination rate                            %
%                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following global variables are passed to           %
% function zeroinfect, for Problem 5.                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global target ntimes n k tau delta   

movieon = 1;
movieoff = 0;

% Plotting symbols:

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Problem 1                                              %
%   Run a Monte Carlo simulation with no mobility or     %
%   vaccination.                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1)
n = 10;
mn = n*n;
k = 4;
tau = .2;
delta = 0;
nu = 0;

[nrecovered1,infected,recovered,susceptible,vaccinated] = ...
          modelmc(n,k,tau,delta,nu,movieon);

t = length(infected);
plot(1:t,infected/mn,pl(2).s,1:t,susceptible/mn,pl(3).s, ...
     1:t,recovered/mn,pl(1).s,'LineWidth',1.5)
legend('Infected','Susceptible','Recovered')
xlabel('day')
ylabel('Proportion of individuals')
title(sprintf('Disease Status with tau = %f',tau))
print -depsc figure1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Problem 2                                              %
%   Add mobility of the patients to the model.           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(2)
delta = .01;

[nrecovered2,infected,recovered,susceptible,vaccinated] = ...
          modelmc(n,k,tau,delta,nu,movieoff);

t = length(infected);
plot(1:t,infected/mn,pl(2).s,1:t,susceptible/mn,pl(3).s, ...
     1:t,recovered/mn,pl(1).s,'LineWidth',1.5)
legend('Infected','Susceptible','Recovered')
xlabel('day')
ylabel('Proportion of individuals')
title(sprintf('Disease Status with tau = %f, delta = %f' ...
       ,tau, delta))
print -depsc figure2

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Problem 3                                              %
%   Add vaccination to the model.                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(3)
nu = .1;

[nrecovered3,infected,recovered,susceptible,vaccinated] = ...
          modelmc(n,k,tau,delta,nu,movieoff);

t = length(infected);
plot(1:t,infected/mn,pl(2).s,1:t,susceptible/mn,pl(3).s, ...
     1:t,recovered/mn,pl(1).s,1:t,vaccinated/mn,pl(4).s,'LineWidth',1.5)
legend('Infected','Susceptible','Recovered','Vaccinated')
xlabel('day')
ylabel('Proportion of individuals')
title(sprintf('Disease Status with tau = %f, delta = %f, nu= %f' ...
       ,tau, delta,nu))
print -depsc figure3

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Problem 4                                              %
%   Investigate the variance of the estimate of          %
%   the number of patients eventually infected           %
%   (nrecovered) by taking 1000 trials using various     %
%   nu values.                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(4)

jj = 1;
for nu=0:.1:.3,

   for i=1:1000,
      recovered(i) = modelmc(n,k,tau,delta,nu,movieoff)/mn;
   end

   subplot(2,2,jj)
   hist(recovered*100)
   xlabel('Percent infected')
   ylabel('Number of trials')
   title(sprintf('Histogram of infection rate for nu = %f',nu))

   disp(sprintf('For nu= %f, Mean = %f,  Variance = %f', ...
         nu,mean(recovered),var(recovered)))
   jj = jj + 1;
end
print -depsc figure4

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Problem 5                                              %
%   Compute the vaccination rate that confines the       %
%   expected infection level to 20 per cent of the       %
%   population.                                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

target = .2;
ntimes = 100;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   We can try to use the zerofinder fzero, but the      %
%   results are unreliable.                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nu = fzero(@zeroinfect,.1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Problem 6                                              %
%   Investigate a Markov Chain model for a 2-day         %
%   infection on a 3-person population.                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

data = [
 2 1 (1-tau)^2
 3 1 tau*tau
 5 1 tau*(1-tau)
 8 1 tau*(1-tau)
13 2 tau*(1-tau)
14 2 tau*tau
15 2 tau*(1-tau)
17 2 (1-tau)^2
 4 3 1
18 4 1
 6 5 tau
11 5 1-tau
 7 6 1
18 7 1
 9 8 tau
12 8 1-tau
10 9 1
18 10 1
19 11 1
16 12 1
12 13 1
 4 14 1
11 15 1
16 16 1
17 17 1
18 18 1
19 19 1];

a = sparse(data(:,1),data(:,2),data(:,3));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute a stationary vector for the matrix.            %
% Since the matrix is reducible, there are many          %
% stationary vectors, so we isolate the one of interest  %
% by using the power method started with a vector z      %
% corresponding to starting in the first state.          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

z = zeros(19,1);
znew = z;
znew(1)=1;

while norm(z-znew) > eps,
   z = znew;
   znew = a*z;
end

prob(1) = znew(17);
prob(2) = znew(16)+znew(19);
prob(3) = znew(18);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compare with results of Monte Carlo trials.            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

i=0;
probmc = zeros(1,3);
while (norm(probmc-prob,'inf')>.005)

   i=i+1;
   nrecovered(i) = modelmc1d(3,2,tau);

   probmc(1) = sum(nrecovered==1)/i;
   probmc(2) = sum(nrecovered==2)/i;
   probmc(3) = sum(nrecovered==3)/i;

end

ntrials = i;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compare with analytical results.                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

probtr(1) = (1-tau)^4;
probtr(2) = 2*(tau*(1-tau)^2 + tau*(1-tau)^3);
probtr(3) = tau^2 + 2*tau^2*(1-tau) + (1-tau)^2*tau^2;

disp('The 3-bed case')
for i=1:3,
  disp(sprintf('probability of %d infection predicted to be',i))
  disp(sprintf( ...
    ' %f by Markov chain, %f by Monte Carlo, and %f analytically',...
                prob(i),probmc(i),probtr(i)))
end

disp(sprintf(...
    'It took %d trials of Monte Carlo to get 2 digits of accuracy',...
    ntrials))
