function BoardingSchool

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   This code fits the data from the English boarding school to an SIR
%   model. The SIR model is simulated using ODE45 using varying values of
%   beta and a fixed value of alpha.
%
%   Inputs:
%   1. alpha
%   2. data for the problem
%   3. tspan is the time mesh we plot the model on
%
%   Ouputs:
%   1. Plot of model fitted to data
%   2. Plot of Residuals
%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
close all;

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

param=[alpha,beta];

%% 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);

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

y0=[S0,I0,R0];

%% Plotting Initial Conditions for model
[~,y] = ode45(@(t,y) SIR(t,y,alpha,beta), tspan, y0);

Iguess=y(:,2);

figure('units','normalized','outerposition',[0 0 1 1]);
plot(tdata,Idata,'ro',tspan,Iguess,'b','linewidth',2)

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

ax=gca;
ax.FontSize = 40;

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

legend('Data of Infected Population','Model Prediction','interpreter','latex','fontsize',40)


%% 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 Model
    function error=L2Error(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(@L2Error,param,options);
alpha=param(1);
beta=param(2);

%% Plotting Model Fit
[~,y] = ode45(@(t,y) SIR(t,y,alpha,beta), tspan, y0);
[~,yres]=ode45(@(t,y) SIR(t,y,alpha,beta), tdata, y0);

Imodel=y(:,2);
Ires=yres(:,2);

figure('units','normalized','outerposition',[0 0 1 1]);
subplot(1,2,1)
plot(tdata,Idata,'ro',tspan,Imodel,'b','linewidth',2)

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

ax=gca;
ax.FontSize = 40;

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

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

%% Plotting Residuals
subplot(1,2,2);
plot(tdata,Ires-Idata,'ro','linewidth',2)

xlabel('$$t$$','interpreter','latex')
ylabel('Residuals','interpreter','latex')

ax=gca;
ax.FontSize = 40;

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