function BoardingSchoolModelComparison

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   This code fits the data from the English boarding school to an SIR
%   model and an SEIR model. The code then does a model comparison between
%   the two using an AIC. 
%
%   Inputs:
%   1. data for the problem
%   2. tspan is the time mesh we plot the model on
%
%   Ouputs:
%   1. Plot of model fitted to data
%   2. Plot of Residuals
%   3. Comparison between the models
%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
close all;

%% Inputs

%% Data (In practice this should be loaded from an external file)

tdata=3:1:14;
Idata=[25,75,227,296,258,236,192,126,71,28,11,7]';
tspan=linspace(3,14,1000);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Fitting SIR
%   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha=1/3;
beta=.0025; %This is set as an initial guess for the force of infection.

param=[alpha,beta];


%% Initial Conditions for SIR model
I0=25;
S0=738;
R0=0;

y0=[S0,I0,R0];


%% SIR ODE System
    function dydt = SIR(t,y,alpha,beta)
        
        %Conversion of y into the variables we are using for the problem. I do this
        %to make the debugging easier and to make the code more readable.
        
        S=y(1);
        I=y(2);
        R=y(3);
        
        
        Sdot=-beta*S*I;
        Idot=beta*S*I-alpha*I;
        Rdot=alpha*I;
        
        %Putting the derivatives into a single row vector dydt.
        dydt(1)=Sdot;
        dydt(2)=Idot;
        dydt(3)=Rdot;
        
        dydt=dydt';
        
        
    end

%% Computing L^2 Error in SIR Model
    function error=L2SIRError(param)
        alpha=param(1);
        beta=param(2);
        [~,y] = ode45(@(t,y) SIR(t,y,alpha,beta),tdata,y0);
        
        Imodel=y(:,2);
        error=sum((Idata-Imodel).^2);
    end


%% Finding the best fit for the model

options = optimset('Display','iter','PlotFcns',@optimplotfval);

[param,fval]=fminsearch(@L2SIRError,param,options);
alpha=param(1);
beta=param(2);
SIRError=L2SIRError(param);

%% Constructing model fit:

[~,y] = ode45(@(t,y) SIR(t,y,alpha,beta), tspan, y0);

ISIR=y(:,2);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Fitting SEIR
%   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha=1/3;
beta=.0025; %This is set as an initial guess for the force of infection.
eta=.3; %This is set as a guess for rate of exposure
E0=2;

param=[alpha,beta,eta,E0];


%% Initial Conditions for SEIR model
I0=25;
R0=0;




%% SEIR ODE System
    function dydt = SEIR(t,y,alpha,beta,eta)
        
        %Conversion of y into the variables we are using for the problem. I do this
        %to make the debugging easier and to make the code more readable.
        
        S=y(1);
        E=y(2);
        I=y(3);
        R=y(4);
        
        
        Sdot=-beta*S*I;
        Edot=beta*S*I-eta*E;
        Idot=eta*E-alpha*I;
        Rdot=alpha*I;
        
        %Putting the derivatives into a single row vector dydt.
        dydt(1)=Sdot;
        dydt(2)=Edot;
        dydt(3)=Idot;
        dydt(4)=Rdot;
        
        dydt=dydt';
        
        
    end

%% Computing L^2 Error in SIR Model
    function error=L2SEIRError(param)
        alpha=param(1);
        beta=param(2);
        eta=param(3);
        E0=param(4);
        
        S0=738-E0;
        
        y0=[S0,E0,I0,R0];
        [~,y] = ode45(@(t,y) SEIR(t,y,alpha,beta,eta),tdata,y0);
        
        Imodel=y(:,3);
        error=sum((Idata-Imodel).^2);
    end


%% Finding the best fit for the model

options = optimset('Display','iter','PlotFcns',@optimplotfval);

[param,fval]=fminsearch(@L2SEIRError,param,options);
alpha=param(1);
beta=param(2);
eta=param(3);
E0=param(4);

S0=738-E0;

y0=[S0,E0,I0,R0];

%% Constructing Model Fit
[~,y] = ode45(@(t,y) SEIR(t,y,alpha,beta,eta), tspan, y0);

ISEIR=y(:,3);
SEIRError=L2SEIRError(param);


%% Plotting Model Fits
figure('units','normalized','outerposition',[0 0 1 1]);
plot(tdata,Idata,'ro',tspan,ISIR,tspan,ISEIR,'linewidth',2)

xlabel('$$t$$','interpreter','latex')
ylabel('Infected Population','interpreter','latex')

ax=gca;
ax.FontSize = 40;

set(gcf,'color','w');

legend('Data','SIR Model','SEIR Model','interpreter','latex','fontsize',40)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   AIC Comparison
%
%   We need to use the AICc since we do not have eneough data points
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[n,~]=size(Idata);
kSIR=3;
kSEIR=5;

AICSIR=n*log(SIRError/n)+2*kSIR+2*kSIR*(kSIR+1)/(n-kSIR-1)
AICSEIR=n*log(SEIRError/n)+2*kSEIR+2*kSEIR*(kSEIR+1)/(n-kSEIR-1)

d1=AICSEIR-AICSIR;
wSIR=1/(1+exp(-d1/2))
wSEIR=exp(-d1/2)/(1+exp(-d1/2))
end