clear; clc;
%  SEIR model, adapted from Andersion, Coombs and Colijn etc.
%  27.4.2020	 matches Excel spreadsheet calculations
%
%  Econ2 - simple constant Xt
%  revised 4.5.20 to add costs for infections also
%  updated 15.5.20 for Salje values  Xt as %
%  15.5.20 revised to include H overflow deaths
%  costs in $k until final summary

%% Initialize Parameters
% control variables
timedays = 450;
xruns = 50;  % xruns-1 is max % shutdown
theta = 0.2; % sd effectivesness = f(x) = x^theta,
xstep = 1.0;
xvals = 0:xstep:(xruns-1); % percent of economy shutdown

% Parameter values
alpha = 0.17;      %rate I to R
R0 = 3.3;
beta = alpha *R0;   %infection rate
gamma = 0.2;        %rate E to I
kb = 1/14;          
kh = 1/10;
ku = 1/15;
kw = 1/2;

% transition probabilities
ph = 0.026;         % prob. needing hospital if infected (ICU or regular)
pb = 1-ph;          % prob. NOT needing hospital if infected 
pu = 0.182;         % prob. needing ICU given hospital is needed
phd = 0.2;          % prob. dying in hospital
pud = 0.22;         % prob. dying in ICU
pwd = 0.8 ;         % prob. dying while waiting for hospital
pd = ph*(1-pu)*phd + ph*pu*pud; % prob. dying when infected (disregarding ICU or hospital overflow)

% Hospital Capacities
Hmax = 5600;        % Maximum H beds
Umax = 300;         % Maximum ICU beds

% economic and population variables
Npop = 5100000;		% pop. BC
GD = 800000;        % GDP BC per day

% costs
VL = 7000;          % value of life
VLvec = 1000:2000:9000; % vector of value of life values
CDH =  1.9;			% cost per day in hospital
CNU = 8.0;			% setup ICU costs
CDU = 2.0;          % cost per day in ICU

 % Initial infected and exposed
IEzero = 1000;

%% Vary VL Loop
vlj = length(VLvec);
for j = 1:length(VLvec)
    VL = VLvec(j);

%% Percent Shutdown Loop
for i = 1:xruns
    % Initilize populations
    
    % Broad populations
    S = Npop-2*IEzero;
    E = IEzero;
    I =  IEzero;
    Q = 0;
    % Q subgroups
    B = 0;
    H = 0;
    U = 0;
    W = 0;
    R = 0;
    D = 0;    

    % Reset old costs - get new percent economic shutdown
    xi = xvals(i)/100;
    medicalCosts_i = 0;    
    distancingCosts_i = 0;
    deathCosts_i = 0;
    maxInfected_i = 0;
    % Epidemic loop
    for t = 1:timedays
       %% Update SEIQ populations       
        betat = beta*(1- xi^theta); % effective infection rate        
        
        St = S; % Old values
        Et = E;
        It = I;
        Qt = Q;
        
        % Update populations
        S = St - betat*It*St/Npop;
        E = Et +  betat*It*St/Npop - gamma*Et;
        I = It + gamma*Et - alpha*It;        
        Q = Qt + alpha*It;        
        % Total = S + E + I + Q             % Sanity check - uncomment population should remain constant                
        
        %% Decompose Q into B,H,U,W,R,D
        % Save old populations for calculations
        Bt = B;
        Ht = H;
        Ut = U;
        Wt = W;
        Rt = R;
        Dt = D;
        
        % Get new Q patients to be sorted into bins
        Qnew = Q - Qt;
        
        % Hospital and ICU bed demand
        Hdemand = Ht*(1-kh) + Wt*(1-kw) + ph*(1-pu)*Qnew; % Needs hospital and not ICU
        ICUdemand = Ut*(1-ku) + ph*pu*Qnew;               % Needs hospital and ICU 
        
        B = Bt*(1-kb) + pb*Qnew;
        H = min(Hmax,Hdemand);
        W = max(0,Hdemand - Hmax); % Those who can't find hostpital space must wait
        U = min(Umax,ICUdemand);             
        R = Rt + kh*(1-phd)*Ht + kw*(1-pwd)*Wt + ku*(1-pud)*Ut + kb*Bt; % Those who don't die recover
        D = Dt + kh*phd*Ht + kw*pwd*Wt + ku*pud*Ut + max(0,ICUdemand - Umax);                
        Unew = max(0,U - Ut); % Used to calculate ICU setup costs
        
%         % Sanity checks - make sure you don't loose track of anyone - totals should stay constant
%         Total = S + E + I + Q;
%         TotalQ = S + E + I + B + H + W + U + R + D
        
        
       %% Calculate costs
       distancingCosts_i = distancingCosts_i + GD*xi;
       medicalCosts_i = (CDH*H) + (CDU*U) + (CNU*Unew);  
       if I > maxInfected_i
           maxInfected_i = I;
       end
       
    end % End time loop
    DeadAfterT = pd*(E + I) + phd*H + pud*U;
    deathCosts_i =  VL*(D + DeadAfterT);    
    % End calculate costs
    
    % Save final costs 
    maxInfected(i) = maxInfected_i;
    distancingCosts(i) = distancingCosts_i;
    medicalCosts(i) = medicalCosts_i;
    deathCosts(i) = deathCosts_i;
    costs(i) = distancingCosts_i + medicalCosts_i + deathCosts_i;               
end
%% Summarys
% Critical values
xc = ((R0-1)/R0)^(1/theta)
VLc = xc*timedays*GD/(pd*Npop)

% Cost summary
[MinCosts,minIdx] = min(costs)
sdMinCost = xvals(minIdx)

allcosts(j,:) = costs;

end % End VL loop

%% Plots
h_costs = figure('Name', 'Shutdown costs', 'Units', 'normalized');
plot(xvals,allcosts, 'LineWidth', 3)
set(gca,'FontSize',14)
set(findall(gcf,'type','text'),'FontSize',14)
xlabel('percent shutdown')
ylabel('costs (thousands of dollars CAD)')
%title('shutdown cost with varying VL')
grid
VLstring = arrayfun( @(VLvec) [' VL = ',num2str(VLvec)], VLvec, 'Uni',0);
legend(VLstring)

   
    