% Mathematical Modeling [MATH 4630]
% P.Parisi
% v6.0  2May2020

%% Message to the Professor

% Dear Dr. Clough,

% The code is broken into 2 sections for 2 topics I focused on.
% It really should be two separate .m files... but to operate well,
% Please use the 'run section' command for each section to run independently.

% Thank you for a great class and help on this project!

%       -Phil



%% Topic 1. What If?
%Curve Varieties (Mountain, Hill, Meadow)

clc, close all, clear all, format compact

%%%%%%%%%%%%%%%%%%% PARAMETERS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Curve Parameters
totalcases = 332000; %total cases of the virus (used to scale values of the normal dist. curves)
Days = 180; %6month period
ICU_capacity = 33000; %110,000 ICU's nationally - 77,000 avg filled up. 
ICU_bedstay = 8; %how long each infected stays in ICU


%%%%%%%%%%%%%%%%%%%% MOUNTAIN CURVE %%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Mountain Curve
mu1 = Days/6; %position of center curve
sigma1 = 8; %how 'fat' the curve is
f1 = @(x) 1/(sigma1*sqrt(2*pi))*exp(-((x-mu1)^2)/(2*sigma1^2)); %normal dist curve, used for new cases shape


%Loop Parameters
CumulativeInfect1(1) = 0;
ICUBeds_alive1(1:Days+15) = 0;


%Loop over the Days
for i = 1:Days
    
     
        %How many COVID Patients come thru today?
        DailyInfect1(i) = totalcases*f1(i); %daily new covid patients
        CumulativeInfect1(i+1) = CumulativeInfect1(i) + DailyInfect1(i); %cumulative COVIDs that entered
     

        %Calculate Deaths
        if DailyInfect1(i) + ICUBeds_alive1(i) > ICU_capacity
           deaths1(i) = DailyInfect1(i) + ICUBeds_alive1(i) - ICU_capacity;
          lives1(i) = DailyInfect1(i) - deaths1(i);
        
        else
        deaths1(i) = 0;
        lives1(i) = DailyInfect1(i) - deaths1(i);
     
        end
     
        %Assign alive patients to ICU beds
        ICUBeds_alive1(i:i+ICU_bedstay-1) = ICUBeds_alive1(i:i+ICU_bedstay-1) + lives1(i);
    
           
end



%%%%%%%%%%%%%%%%%%%% HILL CURVE %%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Hill Curve
mu2 = mu1+30;
sigma2 = sigma1*2;
f2 = @(x) 1/(sigma2*sqrt(2*pi))*exp(-((x-mu2)^2)/(2*sigma2^2));


%Loop Parameters
CumulativeInfect2(1) = 0;
ICUBeds_alive2(1:Days+15) = 0;

%Loop over the Days
for i = 1:Days
    
     
        %How many COVID Patients come thru today?
        DailyInfect2(i) = totalcases*f2(i); %daily new covid patients
        CumulativeInfect2(i+1) = CumulativeInfect2(i) + DailyInfect2(i); %cumulative COVIDs that entered


        %Calculate Deaths
        if DailyInfect2(i) + ICUBeds_alive2(i) > ICU_capacity
           deaths2(i) = DailyInfect2(i) + ICUBeds_alive2(i) - ICU_capacity;
           lives2(i) = DailyInfect2(i) - deaths2(i);
        
        else
        deaths2(i) = 0;
        lives2(i) = DailyInfect2(i) - deaths2(i);
     
        end
     
        %Assign alive patients to ICU beds
        ICUBeds_alive2(i:i+ICU_bedstay-1) = ICUBeds_alive2(i:i+ICU_bedstay-1) + lives2(i);
    
           
end


%%%%%%%%%%%%%%%%%%%% MEADOW CURVE %%%%%%%%%%%%%%%%%%%%%%%%%%
     
%Meadow Curve
mu3 = mu2+30;
sigma = sigma2*2;
f = @(x) 1/(sigma*sqrt(2*pi))*exp(-((x-mu3)^2)/(2*sigma^2));


%Loop Parameters
CumulativeInfect(1) = 0;
ICUbeds(1:Days+15) = 0;

%Loop over the Days
for i = 1:Days
    
     
        %How many COVID Patients come thru today?
        DailyInfect(i) = totalcases*f(i); %daily new covid patients
        CumulativeInfect(i+1) = CumulativeInfect(i) + DailyInfect(i); %cumulative COVIDs that entered
     
        
        %Calculate Deaths
        if DailyInfect(i) + ICUbeds(i) > ICU_capacity
           deaths3(i) = DailyInfect(i) + ICUbeds(i) - ICU_capacity;
           lives3(i) = DailyInfect(i) - deaths3(i);
       
        else
        deaths3(i) = 0;
        lives3(i) = DailyInfect(i) - deaths3(i);
     
        end
     
        %Assign alive patients to ICU beds
        ICUbeds(i:i+ICU_bedstay-1) = ICUbeds(i:i+ICU_bedstay-1) + lives3(i);
    
           
end



%%%%%%%%%%%%%%%%%%%%%% SUBPLOT 1 %%%%%%%%%%%%%%%%%%%%%%%     
figure(1)

%Daily Cases over Time Plot
subplot(1,2,1)
plot(1:Days,DailyInfect1/1000,'b--','LineWidth',2)
hold on
plot(1:Days,DailyInfect2/1000,'r-.','LineWidth',2)
hold on
plot(1:Days,DailyInfect/1000,'m:','LineWidth',2)
title('Daily Cases over Time')
ylabel('Daily Cases [1000s]'),xlabel('Day')
legend('Mountain','Hill','Meadow'), grid on
xlim([1 180]), ylim([0 max(DailyInfect1)*1.1/1000 ])
% text2place1 = {'\mu = 30'; '\sigma = 8'};
% text2place2 = {'\mu = 60'; '\sigma = 16'};
% text2place3 = {'\mu = 180'; '\sigma = 32'};
% gtext(text2place1)
% gtext(text2place2)
% gtext(text2place3)

%Cumulative Cases over Time Plot
subplot(1,2,2)
plot(1:Days+1,CumulativeInfect1/1000,'bd')
hold on
plot(1:Days+1,CumulativeInfect2/1000,'rV')
hold on
plot(1:Days+1,CumulativeInfect/1000,'ms')
title('Total Cases over Time')
ylabel('Total Cases [1000s]'),xlabel('Day')
legend('Mountain','Hill','Meadow','Location','Southeast'), grid on
xlim([1 180]), ylim([0 max(CumulativeInfect1)*1.1/1000 ])
text(140,totalcases*1.05/1000,num2str(totalcases))


%%%%%%%%%%%%%%%%%%%%% SUBPLOT 2 %%%%%%%%%%%%%%%%%%%%%
figure(2)

%Daily Deaths
subplot(3,1,3)
plot(1:Days,deaths1/1000,'b')
hold on
plot(1:Days,deaths2/1000,'r')
hold on
plot(1:Days,deaths3/1000,'m')
title('Daily Deaths over Time')
ylabel('Daily Deaths [1000s]'),xlabel('Day')
legend('Mountain','Hill','Meadow','Location','Northeast'), grid on
xlim([1 180])

%Daily Lives
subplot(3,1,2)
plot(1:Days,lives1/1000,'b')
hold on
plot(1:Days,lives2/1000,'r')
hold on
plot(1:Days,lives3/1000,'m')
title('Daily Lives admitted to ICU over Time')
ylabel('Daily Lives [1000s]'),xlabel('Day')
legend('Mountain','Hill','Meadow','Location','Northeast'), grid on
xlim([1 180])

%ICU Capacity
subplot(3,1,1)
plot([1 Days], [ICU_capacity/1000 ICU_capacity/1000],'k--','LineWidth',2)
hold on
plot(1:Days+15,ICUBeds_alive1/1000,'b')
hold on
plot(1:Days+15,ICUBeds_alive2/1000,'r')
hold on
plot(1:Days+15,ICUbeds/1000,'m')
title('ICU Occupancy over Time')
ylabel('ICU Beds Full [1000s]'),xlabel('Day')
legend('ICU Capacity','Mountain','Hill','Meadow','Location','Northeast'), grid on
xlim([1 180])


%%%%%%%%%%%%%%%%%%% COMMAND WINDOW OUTPUTS %%%%%%%%%%%%%%%%%%

%Calculate Total Deaths for each Curve
DeathsMountain = ceil(sum(deaths1));
DeathsHill = ceil(sum(deaths2));
DeathsMeadow = ceil(sum(deaths3));

    %Assemble death counts into a cute lil' table
    table1 = table(DeathsMountain,DeathsHill,DeathsMeadow);
    disp(table1)

%Display Max Values from Daily Cases graph
PeakDeathsMountain = ceil(max(DailyInfect1))
PeakDeathsHill = ceil(max(DailyInfect2))
PeakDeathsMeadow = ceil(max(DailyInfect))









%% Topic 2. Plan Ahead.
% Determining hopsital beds and ICUs needed for a pandemic
% randomness introduced to replicate a realistic scenario

clc, close all, clear all, format compact

% We Loop over DAYS and the PATIENTS every day
% Thus, the modeled is scaled down to run with 10,000 totalcases

%Hill Curve Model
Days = 100;
totalcases = 5000000;
mu = Days/2; %mean --> center of curve
sigma = 16; %std dev. --> 'width' of curve
f = @(x) 1/(sigma*sqrt(2*pi))*exp(-((x-mu)^2)/(2*sigma^2)); %curve shape model, normal dist.

%Infection Parameters
alpha = .7; % 70% of those with infection are hospitalized
beta = .4; % 40% of those hopsitalized go to ICU
delta = .95; %survival rate in ICU
gamma = .995; %survival rate in hospital bed

%Loop Parameters
CumulativeInfect(1) = 0;
CumulativeInfect_hospital(1) = 0;
DailyInfect_hospital(1:Days) = 0;
ICUbeds(1:Days+16) = 0;
Hospitalbeds(1:Days+16) = 0;
deaths(1:Days) = 0;
lives(1:Days) = 0;


disp('This takes about 15 seconds to run...')

%Loop over the Days
for i = 1:Days
    
    %How many infected this day?
    DailyInfect(i) = ceil(totalcases*f(i)*(1 + (rand(1)*.4 - .2)));  %daily new covid patients, +/- 20% of nominal
    CumulativeInfect(i+1) = CumulativeInfect(i) + DailyInfect(i); %cumulative infected
     
    
    %Loop over every infected person
    for k = 1:DailyInfect(i)
    
    
    if rand(1) < alpha % 70% show symptoms and go to hospital
        
        %count how many of infected actually go to the hospital
        DailyInfect_hospital(i) = DailyInfect_hospital(i) + 1;
        
        if rand(1) < beta % 40% go to ICUs
            
            %had symptoms, went to the ICU, how long is stay?
            ICU_bedstay = ceil(rand(1)*9 + 7); %8-16 day stay duration
            
            %put patient into ICU beds!
            ICUbeds(i:i+ICU_bedstay-1) = ICUbeds(i:i+ICU_bedstay-1) + 1;
            
            %does patient survive?
            if rand(1) > delta % 5% of patients, even after an ICU stay, will die
                deaths(i) = deaths(i) + 1;
            else % 95% survive after ICU stay
                lives(i) = lives(i) + 1;
            end
            
        else % 60% go to regular Hospital Bed
            
            %had symptoms, went to hospital bed, how long is stay?
            Hospital_bedstay = ceil(rand(1)*7 + 1); %2-8 day stay duration

            %put patient into hospital beds!
            Hospitalbeds(i:i+Hospital_bedstay-1) = Hospitalbeds(i:i+Hospital_bedstay-1) + 1;
            
            %does patient survive? much higher than in ICU!!
            if rand(1) > gamma % 0.5% of patients die suddenly when in hospital beds - surprisingly!!
                deaths(i) = deaths(i) + 1;
            else % 99.5% survive from hospital bed stay
                lives(i) = lives(i) + 1;
            end
            
        end
            
    else %Nothing happens to those that didn't go to hospital. They live. Asymptomatic.
    end
    
    end
    
    %Calculate cumulative hopsital visits
    CumulativeInfect_hospital(i+1) = CumulativeInfect_hospital(i) + DailyInfect_hospital(i);
    
end


%%%%%%%%%%%%%%%%%%% COMMAND WINDOW OUTPUTS %%%%%%%%%%%%%%%%%%%%%%%

disp('Finished!')
HospitalBedsNeeded = max(Hospitalbeds)
ICUsNeeded = max(ICUbeds)
TotalDeaths = sum(deaths)

%%%%%%%%%%%%%%%%%%% SUBPLOT 1: Infection Cruve %%%%%%%%%%%%%%%%%%%
figure(1)

%Daily Infect curve
subplot(1,2,1)
plot(1:Days, DailyInfect/1000,'b'), grid on
hold on
plot(1:Days, DailyInfect_hospital/1000,'m')
title('Daily Infected over Time')
ylabel('New Cases per Day [1000s]'), xlabel('Days')
legend('Infected','Hospitalized')
xlim([0 Days])

%Cumulative Infect curve
subplot(1,2,2)
plot(1:Days+1, CumulativeInfect/1000000,'bd'), grid on
hold on
plot(1:Days+1, CumulativeInfect_hospital/1000000,'ms')
title('Cumulative Infections over Time')
ylabel('Cumulative Cases per Day [Millions]'), xlabel('Days')
legend('Infected','Hospitalized','Location','Northwest')
xlim([0 Days])

%%%%%%%%%%%%%%%%%%%% SUBPLOT 2: Deaths & Beds %%%%%%%%%%%%%%%%%%%%%
figure(2)

%Daily Deaths over Time
subplot(3,1,1)
plot(1:Days,deaths,'m'), title('Daily Deaths over Time')
ylabel('Daily Deaths'), xlabel('Days'), grid on
xlim([0 Days])

%ICU bed occupancy
subplot(3,1,2)
plot(1:Days+16,ICUbeds/1000,'m'), title('ICU Occupancy over Time')
ylabel('ICUs Full [1000s]'), xlabel('Days'), grid on
xlim([0 Days])

%Hospital bed occupancy 
subplot(3,1,3)
plot(1:Days+16,Hospitalbeds/1000,'m'), title('Hospital Bed Occupancy over Time')
ylabel('Hosp. Beds Full [1000s]'), xlabel('Days'), grid on 
xlim([0 Days])