function[y] = covidEconFunc_pareto(sds)
numIntervals = length(sds);
%  August 18/20 V2 with added CB cost of day in bed
%  SEIR model, adapted from Andersion, Coombs and Colijn etc.
%  costs in $k until final summary

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

% control variables
sd_index = 1;
theta = 0.2; % sd effectivesness = f(x) = x^theta,
Nvac = 20400; % Number of vaccines that can be administered in one day
vac_available = 360; % time before vaccine plus time to vaccinate all
time_to_vac = ceil(Npop/Nvac);
timedays = vac_available  + time_to_vac;


% Parameter values
alpha = 0.17;      %rate I to R
R0 = 3.8;  % default is 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.021;         % prob. needing hospital if infected (ICU or regular)
pu = 0.0047;         % prob. needing ICU (Davies paper)
pb = 1-ph - pu;          % prob. NOT needing hospital or icu if infected
phd = 0.17;          % prob. dying in hospital
pud = 0.27;         % prob. dying in ICU
pwd = 0.8 ;         % prob. dying while waiting for hospital
pd = ph*phd + pu*pud; % prob. dying when infected (disregarding ICU or hospital overflow)

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

% costs (in thousands of dollars)
VL = 7000;          % value of life
CDH =  1;			% cost per day in hospital
CNU = 8.0;			% setup ICU costs
CDU = 2.0;          % cost per day in ICU
CB = 0.1;

% Initial infected and exposed
IEzero = 200;



% Initilize populations

% Broad populations
S = Npop-2*IEzero;
E = IEzero;
I =  IEzero;
Q = 0;
V = 0;
% Q subgroups
B = 0;
H = 0;
U = 0;
W = 0;
R = 0;
D = 0;

% Reset old costs - get new percent economic shutdown
medicalCosts_i = 0;
distancingCosts_i = 0;
deathCosts_i = 0;
maxInfected_i = 0;

% Epidemic loop
xi = sds(sd_index);
betat = beta*(1- xi^theta); % effective infection rate
for t = 1:timedays
    
    if t >= vac_available
       vaccinated = Nvac;
    else
        vaccinated = 0;
    end
    
    if mod(t,ceil((timedays)/(numIntervals))) == 0 && t ~= timedays     
       sd_index = sd_index + 1;
       xi = sds(sd_index);
       betat = beta*(1- xi^theta); % effective infection rate
    end
    %% Update SEIQ populations
   
    
    St = S; % Old values
    Et = E;
    It = I;
    Qt = Q;
    Vt = V;
    
    % Save infected pop
    Is(t) = It;
    Qs(t) = Qt;
    Ss(t) = St;
    % Update populations
    S = max(0,St - betat*It*St/Npop - vaccinated);
    E = Et +  betat*It*St/Npop - gamma*Et;
    I = It + gamma*Et - alpha*It;
    Q = Qt + alpha*It;
    V = Vt + min(vaccinated,St-betat*It*St/Npop);
    % 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*Qnew; % Needs hospital and not ICU
    ICUdemand = Ut*(1-ku) + pu*Qnew;               % Needs hospital and ICU
    
    
    Overflow_wait(t) = max(0,Hdemand - Hmax); % save overflow populations
    Overflow_dead(t) = max(0,ICUdemand - Umax);
    
    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 + V
    %         TotalQ = S + E + I + V + B + H + W + U + R + D
    
    %% Calculate costs
    distancingCosts_i = distancingCosts_i + GD*xi;
    medicalCosts_i = medicalCosts_i + (CDH*H) + (CDU*U) + (CNU*Unew) + CB*B;
    if I > maxInfected_i
        maxInfected_i = I;
    end    
end % End time loop without vaccine
deathCosts_i =  VL*D;
% End calculate costs


% Save final costs
max_infected = maxInfected_i;
%% For pareto opt. do not include deaths as part of the cost
costs = distancingCosts_i + medicalCosts_i;


y = [costs;D];
%y = y(1); % Comment for 1 variable optimization


% if costs < 4.5e7 || costs > 5.5e7
%     y = y.^2;;
% end
    
end