function main
% This is MATLAB code for ODE yaws simulations
% Authors Presley Kimball, Jacob Levenson, Amy Moore, Jan Rychtar, Dewey Taylor
% Email questions to Jan Rychtar rychtarj@vcu.edu

% Version from December 7, 2021
%-----------
% Written on Matlab 2020a (version 9.8.0.1417392 (R2020a) Update 4')
% with optimization toolbox

% The program has several major parts, each starts with %%
%
% The main parts are as follows
%
% (Part 0) Technical preliminaries (Users do not have to touch)
%
% (Part 1) Setting up parameters (names and values)
%     This is a large part further divided into
%      1a) List of all parameters and their base values
%      1b) Technical part to set the parameters as empty variables to make them global.
%      1c) Set auxiliary variables as functions of parameters
%      1d) Actually set the parameter values to their base values
%
% (Part 2) Setting up compartments
%
% (Part 3) Setting up the equations for the system
%
% (Part 4) Actual code to run
%
% (Part 5) Codes of the functions that depend on the model
%       This is the actual core of the code. However, this part generally
%       does not have to be touched by a casual user that is interested
%       only in modifying a few parameter values
%
% This includes code for
% - analytical solutions of the equilibria
% - analytical solution of the basic reproduction number
% - solving when R0=1
% - calibration
% - uncertainty analysis
% - validation
%
% (Part 6) Codes of functions that do not depend on model
%     This is a large part that essentially should not be modified by a
%     casual user
%     Functions in this part should be working for all sorts of models
%      - Resetting Parameters Functions
%      - Sensitivity analysis Functions
%      - Setting up system of equations to use for ODE solvers
%      - numerically solving the system of equations
%      - checking for consistency
%
% (Part 7) Some auxiliary functions for better plotting, easy saving etc


% The important textual display is saved in the file yaws_code_output.txt


%% (Part 0) Technical preliminaries

% Close all figures
close all;

% Clear all variables
clear all;

% Set the randomg number generator for reproducibility
rng(0)

% Check if the folder for figures exists and create if not
if not(isfolder('Figures'))
    mkdir('Figures')
end

if license('test','optimization_toolbox')~=1
    disp_Diary('Optimization toolbox is not installed')
    disp_Diary('Calibration and uncertainty functions may not run properly')
    pause
end


% Set the file for textual output and write the date and time we run the code
% After multiple run, the last output is at the bottom of this text file
diary 'yaws_code_output.txt'
disp_Diary('--------------- New run  ----------------')
disp_Diary(datetime)

%% (Part 1) Setting up parameters (names and values)

% Part 1a - list of all parameters, base values and ranges

% Set the main parameters
% Each line is for one parameter
% The first column is the name that will be used in Matlab
% The second column is the label that will be used for plotting
% The third is the base value
% The fourth is the minimal and the fifth is the maximal value (for
% sensitivity analysis)

Params = ...
    {'Lambda', '$\Lambda$',  27.2/12000, 0.001,   0.003;  % Birth rate
    'mu',      '$\mu$',      1/(65*12),  0.0008,  0.0017; % Mortality rate
    'beta',    '$\beta$',    0.014273,   0.01,    0.02;   % Transmission rate 
    'sigma',   '$\sigma$',   30/21,      1/3,     10/3;   % Incubation rate
    'lambda1', '$\lambda_1$',1/3,        1/6,     1/2;    % Healing rate of primary yaws
    'lambda2', '$\lambda_2$',1/3,        1/6,     1/2;    % Healing rate of secondary yaws
    'rho1',    '$\rho_1$',   1/1.5,      1/24,    1;      % Primary latency period length
    'rho2',    '$\rho_2$',   1/25,       1/60,    1/12;   % Secondary latency period length
    'pL2Y2',   '$p_{L_2Y_2}$', 0.9999,   0.9998,  1;      % Probability of relapse from secondary latent to secondary yaws
    'pY1Y2',   '$p_{Y_1Y_2}$', 0.12,     0.09,    0.15;   % Probability of primary yaws to secondary yaws
    'pY2Y3',   '$p_{Y_2Y_3}$', 0.0001,   0,       0.0002; % Probability of secondary yaws to tertiary yaws
    };


% Part 1b - Set the parameters as global variables

% This enables them to be passed from a function to function without actually calling them
% It may not be the cleanest way to do it from coding perspective, but
% increases legibility of the code in some places

% For functionality, we have to manual assign of all the variables here
% once as empty
Lambda = [];
mu = [];
epsilon = 0;  % This one hardoced and is not going to change
beta = [];
sigma = [];
lambda1 = [];
lambda2 = [];
rho1 = [];
rho2 = [];
pL2Y2 = [];
pY1L1 = [];
pY1Y2 = [];
pY2L2 = [];
pY2Y3 = [];
% The ones below are for treatment rates
cE = [];
cY1 = [];
cY2 = [];
cY3 = [];
cL1 = [];
cL2 = [];
tauM = 0; % Just initialize the treatment rate as 0



% Part 1c - Set the auxiliary variables and make them global parameters

% (they are not parameters but are calculated from them easily and
% they simplify the model)
% Each line is for one aux variable
% The first column is the name that will be used in Matlab
% The second column is the label that will be used for plotting
% The third is the formula

AuxVars = {'tauE' , '$\tau_E$', 'cE .* tauM';
    'tau1' , '$\tau_{Y_1}$', 'cY1 .* tauM';
    'tau2' , '$\tau__{Y_2}$', 'cY2 .* tauM';
    'tau3' , '$\tau_{Y_3}$', 'cY3 .* tauM';
    'tauL1', '$\tau_{L_1}$', 'cL1 .* tauM';
    'tauL2', '$\tau_{L_2}$', 'cL2 .* tauM';
    'pY1L1', '$p_{Y_1L_1}$', '1-pY1Y2';
    'pY2L2', '$p_{Y_2L_2}$', '1-pY2Y3';
    'N', '$N$', 'Lambda./mu' };

% Set all auxiliary variables as global
tauE = [];
tau1 = [];
tau2 = [];
tau3 = [];
tauL1 = [];
tauL2 = [];
N=[];

% Part 1d - Actually set the values of the parameters

% Get the number of parameters, discard the number of columns
[numParams,~] = size(Params); %

% Get the number of auxiliary variables, discard the number of columns
[numAuxVars,~] = size(AuxVars); %

% Get the parameter names
Names = Params(:,1);

% Get the parameter labels
Labels = Params(:,2);

% Get the parameter base values
BaseValues = Params(:,3);

% Actually set the parameters
ResetParameters

% Tentatively set treatment as total community treatment
TCT;




%% (Part 2) Set up the compartments

Comps = {'S', '$S$';    % Susceptible
    'E',  '$E$';   % Exposed
    'Y1', '$Y_1$'; % Primary yaws
    'Y2', '$Y_2$'; % Secondary yaws
    'Y3', '$Y_3$'; % Tertiray yaws
    'L1', '$L_1$'; % Primarily Latently infected
    'L2', '$L_2$'; % Secondarily Latently infected
    };
% Set the compartments as empty global
S = [];
E = [];
Y1 = [];
Y2 = [];
Y3 = [];
L1 = [];
L2 = [];

dS = [];
dE = [];
dY1 = [];
dY2 = [];
dY3 = [];
dL1 = [];
dL2 = [];


% Get the number of compartments, discard the number of columns
[numComps,~] = size(Comps); %

%% Part 3 set up the equations of the system

    function equations = getEquations
        % Set the equations
        
        % This returns the right hand sides of the equations
        % It is done via a function so that the equations are set when
        % needed. They can be used for symbolic solutions as well as for
        % numerical solutions or even for differential equations
        
        equations = { ...
            Lambda+tauE*E+(tau1 +(1-(pY1Y2+pY1L1))*lambda1)*Y1+(tau2+(1-(pY2Y3+pY2L2))*lambda2)*Y2+tau3*Y3+tauL1*L1+tauL2*L2-(epsilon+beta*(Y1/(S+E+Y1+Y2+Y3+L1+L2)+Y2/(S+E+Y1+Y2+Y3+L1+L2))+mu)*S; %dS/dt
            (epsilon+beta*(Y1/(S+E+Y1+Y2+Y3+L1+L2)+Y2/(S+E+Y1+Y2+Y3+L1+L2)))*S-(sigma+tauE+mu)*E; %dE/dt
            sigma*E-(lambda1+tau1+mu)*Y1; %dY1/dt
            pY1Y2*lambda1*Y1+rho1*L1+pL2Y2*rho2*L2-(lambda2+tau2+mu)*Y2; %dY2/dt
            pY2Y3*lambda2*Y2+(1-pL2Y2)*rho2*L2-(tau3+mu)*Y3; %dY3/dt
            pY1L1*lambda1*Y1-(rho1+tauL1+mu)*L1; %dL1/dt
            pY2L2*lambda2*Y2-(rho2+tauL2+mu)*L2; %dL2/dt
            };
    end


%% Part 4  Actual code to run starts here

% Check for errors in the code/math
% Run with an argument larger than 1 for more through check
% It will generate random parameters specified number of times and
% crosschecks the formulas and analytical solutions and R0 etc
% checkConsistency(1)






        
% Calibrate, i.e. find beta, and display the outcomes
calibration('display')

% Calculate the value of R0
disp(['R0 is ' num2str(R0)])

%pause

% Produce the validation graphs
Samples = 100;
validation2



dependence_on_cL2


%Compute Critical values for tau for TCT and TTT
%Create a plot of value of tau versus value of R0
plotR0onTau


% Do the global stability
globalStability(1000)


% Do the uncertainty analysis
% Set to 100 for good and 1000 for better
% The function uncertainty saves the results, so make sure it runs again
% with more samples by renaming the saved file
Samples = 1000; % How many samples for sensitivity analysis

% Actually set the parameters
ResetParameters


% Note that it generally takes a long time to run the uncertainty analysis
% Thus, the procedure loads previously saved results and one should rename
% or delete the file YawsUncertainty.mat if one want to see new pictures.
uncertainty
return


%% (Part 5) Essential Functions that depend on the model

    function TTT(coverage)
        % Resets treatments probabilities for targeted treatments
        % This is the WORST CASE SCENARIO treating only the active yaws
        if nargin == 0
            % if the function is called without an argument,
            % assume the coverage is 1
            coverage = 1;
        end
        
        cE = coverage;
        cY1 = coverage;
        cY2 = coverage;
        cY3 = coverage;
        cL1 = coverage; %coverage;
        cL2 = 0.1; %0.01;
        CalculateAuxiliaryVariables;
    end

    function TCT(coverage)
        % Resets treatment probabilities for total community treatment
        if nargin == 0
            % if the function is called without an argument,
            % assume the coverage is 1
            coverage = 1;
        end
        
        cE = coverage;
        cY1 = coverage;
        cY2 = coverage;
        cY3 = coverage;
        cL1 = coverage;
        cL2 = coverage;
        CalculateAuxiliaryVariables;
    end


    function sols = getAnalyticalSolutions
        % Returns analytical solutions to the equilibrium
        % equations
        % This function is coded in full generality for nonzero epsilon
        
        % Define constants for simplification
        vE = sigma + tauE + mu;
        vY1 = lambda1 + tau1 + mu;
        vY2 = lambda2 + tau2 + mu;
        vY3 = mu + tau3;
        vL1 = rho1 + tauL1 + mu;
        vL2 = rho2 + tauL2 +mu;
        
        kY1 = sigma / vY1;
        kL1 = (pY1L1 * lambda1 * sigma) / (vY1 *vL1);
        kY2 = ((pY1Y2 * lambda1 * kY1 + rho1 * kL1)* vL2) / (vY2*vL2 - pL2Y2 * rho2 * pY2L2 * lambda2);
        kL2 = (pY2L2 * lambda2 * kY2) / vL2;
        kY3 = (pY2Y3 * lambda2 * kY2 + (1-pL2Y2)* rho2 *kL2) / vY3;
        kN = 1+kY1 +kY2 +kY3  + kL1 + kL2;
        
        w0 = -epsilon;
        w1= vE+ kN*epsilon -vE*R0;
        w2 = vE*kN*R0;
        
        % Solve the quadratic for E/N, pick only the positive solution
        EoverN = (-w1+sqrt((w1)^2-(4*w2*w0)))/(2*w2); %Positive solutions
        % Get E by mutliplying E/N by N
        ansE = EoverN * N;
        
        if ansE ~= 0
            % If endemic equilibrium
            aS = vE * ansE/ (epsilon+beta*(kY1+kY2)*EoverN);
        else
            % Disease-free equilibrium
            aS = N;
        end
        
        % Calculate the rest of the compartments as constant times E
        aE = ansE;
        aY1 = kY1*ansE;
        aY2 = kY2*ansE;
        aY3 = kY3*ansE;
        aL1 = kL1*ansE;
        aL2 = kL2*ansE;
        
        % Return the analytical solutions
        sols = [aS, aE, aY1, aY2, aY3, aL1, aL2];
    end

    function output = R0
        % Returns the basic reproduction number for the system
        
        % Calculate the outgoing rates for the simplification
        vE = sigma + tauE + mu;
        vY1 = lambda1+tau1+mu;
        vY2 = lambda2+tau2+mu;
        vY3 = mu + tau3;
        vL1 = rho1+tauL1+mu;
        vL2 = rho2+tauL2 + mu;
        
        % Return R0
        output = (beta.*sigma)./(vE.*vY1) + (beta.*(lambda1.*pY1L1.*rho1.*sigma.*vL2 + ...
            lambda1.*pY1Y2.*sigma.*vL1.*vL2))./(vE.*vL1.*(vL2.*vY1.*vY2 - lambda2.*pL2Y2.*pY2L2.*rho2.*vY1));
        
    end


    function output = calibration(input)
        % Calibrates the model by finding the values of beta, pY1L1 and
        % pY2L2 to fit data from Lihir Island
        
        % "input" is an optional parameters. If any is given, it will
        % display the results of the calibration.
        
        % output is the error after the calibration
        
        function output = getPrevalences(params)
            % Calculates prevalences of S, active and latent in the
            % population,
            % params give beta
            % output is the prevalences
            
            % Set paramters as given
            beta =  params(1);
            CalculateAuxiliaryVariables;
            
            % Get analytical solutions
            sol = getAnalyticalSolutions;
            % Output the prevalences
            output = [sol(1)+sol(2), sol(3)+sol(4)+sol(5), sol(6)+sol(7)]/sum(sol);
        end
        
        % Do the optimization only if toolbox is installed
        if license('test','optimization_toolbox')==1
            % Data from observations Mitja et al (2018)
            realData = [0.797, 0.024, 0.179];
            
            % Set the function to minimize (it is the norm distance of our estimated
            % prevalences from the real data)
            fun = @(params) norm((realData - getPrevalences(params))./realData);
            % Set the lower bound constraint (rates and probabilities are non-negative)
            lb = [0];
            % Set the upper bound constraint (beta  probabilities have to add to no more than 1)
            ub = [1];
            
            % Set initial guess
            init = [0];
            
            % There are no linear or nonlinear constraints, so set those arguments to [].
            A = [];
            b = [];
            Aeq = [];
            beq = [];
            nonlcon = [];
            % Suppress the display of minimization results
            options = optimoptions('fmincon','Display','off');
            % Do the actual minimization subject to the constraints specified above
            x = fmincon(fun,init,A,b,Aeq,beq,lb,ub, nonlcon, options);
        else
            disp('Optimization toolbox not installed')
            disp('Using hardcoded values')
            x = beta;
        end
        
        % Unpack the results
        beta = x;
        
        % Calculate the error of the approximation
        output = fun(beta);
        
        if nargin>0
            disp(['Calibration complete. Optimal values are ' ...
                ' beta = ' num2str(beta) ]);
            disp(['Real data are ' num2str(realData)])
            disp(['Our estimated values are ' num2str(getPrevalences(x))])
            disp(['Error is ' num2str(output)])
        end
        
    end

    function results = solve_for_tau
        % Numerically solves where  R0=1 and outputs the solution
        
        % Store what the original value of tauM is
        tauOriginal = tauM;
        
        % Check whether there is a solution or not
        tauM = 0; % sets treatment to 0
        CalculateAuxiliaryVariables % this sets all treatment rates appropriate values
        
        if R0<1
            results = 0; % R0<1 even without treatment, so the solution is 0
        else
            
            % Create a symbolic parameter for tau
            TAU = sym('TAU');
            
            % Set tau as a symbolic variable
            tauM = TAU;
            CalculateAuxiliaryVariables
            
            % Set the equation as R0=1.
            eqn = (R0 == 1);
            
            % Numerically solve the equation
            % There will be 5 solutions, only one is positive and the maximal
            % is the right solution
            
            tauSolve = vpasolve(eqn, tauM);
            
            results = tauSolve(tauSolve>0);
            
        end
        % Restore the original value of tau
        tauM = tauOriginal;
        CalculateAuxiliaryVariables;
    end

    function results = solve_for_tau_TTT
        % Numerically solves where  R0=1 and outputs the solution under TTT
        
        TTT;
        results = solve_for_tau;
        
    end



    function dependence_on_cL2
        
        tauM = 0;
        % Get init conditions based on calibrated model (i.e. close to
        % current true data)
        ResetParameters;
        CalculateAuxiliaryVariables;
        InitConditions  = getAnalyticalSolutions;
        
        all_cL2_values = linspace(0.01,1,100);
        for aux1 = 1:length(all_cL2_values)
            
            % Start with a year of TCT treatment
            tspan = [0, 12]; % 1 year
            tauM = 1/6;
            TCT;
            
            % Numerically solve the ODE
            [t,y] = ode45(@(t,y) ODEmodel(t,y), tspan, InitConditions);
            
            % Get prevalence of all infections as a function of time
            I_TCT = 1-y(:,1)./sum(y,2); %(y:,1) is S, sum(y,2) is a sum of all (as a function of a time)
            
            % Get the times needed for elimination
            decrease = 1000;
            
            % Now we will simulate TTT after two rounds of TCT
            % See where we stopped after a year (two rounds of TCT)
            Init_afterTCT = y(end,:);
            
            % Reset the time to run again for 40 years
            tspan = [0, 480];
            
            % Set the treatment rate to reflect the rate from the
            % validation
            tauM = 1/6;
            % Set the value of cL2
            cL2 = all_cL2_values(aux1);
            CalculateAuxiliaryVariables;
            
            % Numerically solve the ODE
            [t,y] = ode45(@(t,y) ODEmodel(t,y), tspan, Init_afterTCT);
            
            % Get prevalence of all infections as a function of time
            I_TTT = 1-y(:,1)./sum(y,2); %(y:,1) is S, sum(y,2) is a sum of all (as a function of a time)
            % Get the truly initial conditions
            I_TTT(1) = 1-InitConditions(1)/(sum(InitConditions));
            
            
            % Get how long it will take to decrease the prevalence given number of times
            if I_TTT(end)<I_TTT(1)/decrease
                % If it got low enough, get the first time it went below
                % Add 12 for the first 12 months of TCT
                tElimTTT(aux1) = min(t(I_TTT<I_TTT(1)/decrease))+12;
            else
                % It did not get low enough, add 2 extra years
                tElimTTT(aux1) = 480+12+24;
            end
            
        end
        
        
        % Plot the figure of how TTT elim time depends on tau2
        figure
        plot(all_cL2_values*100, tElimTTT/12, 'k-')
        ylabel('Time (in years)')
        xlabel('Percentage of $L_2$ treated in TTT')
        FormatFigure
        SaveFigure(['tElimTTT_cL2'])
        
    end



    function uncertainty
        % Runs uncertainty analysis
        % Does PRCC uncertainty and sensitivity analysis
        % Based on Marino, S., Hogue, I. B., Ray, C. J., and Kirschner, D. E. (2008).
        % A methodology for performing global uncertainty and sensitivity analysis in systems
        % biology. Journal of Theoretical Biology, 254(1):178–196
        
        
        % Simulates a bunch of parameter values and calculates R0, how long
        % it will take to decrease the disease prevalence and
        % what will be the disease prevalence in 10 years
        
        paramsUn = []; % For uncertainty analysis
        
        % Check whether we really have to run the analysis again or if we
        % already did it and want just to see the plots again
        
        TTT; % This is just to set cL2 to an appropriate value
        savedFileName = ['YawsUncertainty' num2str(cL2) '.mat']
        
        if isfile(savedFileName) == 1
            % Already saved the results once
            disp('Already run uncertainty, just loading the results')
            disp('If you want to redo, delete the file YawsUncertainty.mat')
            load(savedFileName,'R0Un','paramsUn', 'SolsUn', 'tElimTCT', 'tElimTTT', 'R0treatment', 'tEst', 'error' )
            
        else
            % No results saved
            disp('Running uncertainty analysis, this may take a while')
            
            tauM = 0;
            % Get init conditions based on calibrated model (i.e. close to
            % current true data)
            ResetParameters;
            CalculateAuxiliaryVariables;
            InitConditions  = getAnalyticalSolutions;
            
            % For the number of experiments, generate parameters randomly and
            % see what the outcomes are
            for aux1 = 1:Samples
                disp(['Running simulation ' num2str(aux1) ' of ' num2str(Samples) ])
                if aux1 == 1
                    % Use the base values for the first run
                    % We will just set R0 as it is and then the values will not
                    % have to be generated again
                    R0Un(aux1) = R0;
                    error(aux1) = 0;
                else
                    % Tentatively set R0 to 0 so that we can discard parameters
                    % that produce "bad" R0 and too large error
                    R0Un(aux1) = 0;
                    error(aux1) = 1;
                end
                
                % Keep assigning parameter values randomly until 3>R0>1 and
                % fitting error is small
                while error(aux1)>0.005 %(R0Un(aux1)<=1)||(R0Un(aux1)>3)||
                    % Assign parameter values at random (but within the bounds)
                    for aux = 1:numParams
                        eval([Names{aux} '=' num2str((Params{aux,5}-Params{aux,4})*rand+Params{aux,4}) ';'])
                    end
                    % Set tau = 0 for  no treatment
                    tauM = 0;
                    
                    % Calibrate to find beta and qY1 and qY2, store fittting
                    % error for next round of "while cycle"
                    error(aux1) = calibration;
                    % Calculate the auxiliary variables
                    CalculateAuxiliaryVariables
                    
                    % Evaluate R0
                    R0Un(aux1) = R0;
                end
                % Now we have random parameters such that 1<R0<3 and good fit
                % to data
                % Store the parameters
                for aux = 1:numParams
                    % Go through all the parameters
                    % Evaluate the command that assigns the parameter with a given name the proper base value
                    eval(['paramsUn(aux1, aux) = ' Names{aux} ';'])
                end
                
                % Store endemic equilibria for these parameters
                SolsUn(aux1, :) = getAnalyticalSolutions;
                SolsUn(aux1, :) = SolsUn(aux1, :)/sum(SolsUn(aux1,:));
                
                %% See how many cases we will have in 10years of TCT treatment
                tspan = [0, 120]; % 12 years is 120 months
                tauM = 1/6;
                TCT;
                
                %InitConditions
                
                %pause
                % Numerically solve the ODE
                [t,y] = ode45(@(t,y) ODEmodel(t,y), tspan, InitConditions);
                
                
                
                % Get prevalence of all infections as a function of time
                I_TCT = 1-y(:,1)./sum(y,2); %(y:,1) is S, sum(y,2) is a sum of all (as a function of a time)
                
                % Get the times needed for elimination
                decrease = 1000;
                if I_TCT(end)<I_TCT(1)/decrease
                    % If it got low enough, get the first time it went below
                    % Add 12 for the first 12 months of TCT
                    tElimTCT(aux1) = min(t(I_TCT<I_TCT(1)/decrease));
                else
                    % It did not get low enough, add 2 extra years
                    tElimTCT(aux1) = 120+24;
                end
                
                % Now we will simulate TTT after two rounds of TCT
                % See where we stopped after a year (two rounds of TCT)
                Init_afterTCT = y(find(t>12, 1),:);
                
                % Reset the time to run again for 40 years
                tspan = [0, 480];
                
                % Set the treatment rate to reflect the rate from the
                % validation
                tauM = 1/6; %tauM = 1/3;1/3
                TTT;
                
                
                
                % Numerically solve the ODE
                [t,y] = ode45(@(t,y) ODEmodel(t,y), tspan, Init_afterTCT);
                
                % Get prevalence of all infections as a function of time
                I_TTT = 1-y(:,1)./sum(y,2); %(y:,1) is S, sum(y,2) is a sum of all (as a function of a time)
                % Get the truly initial conditions
                I_TTT(1) = 1-InitConditions(1)/(sum(InitConditions));
                
                if aux1 == 1
                    figure
                    plot(t,I_TTT)
                    xlabel('Time (in months)')
                    ylabel('Proportion of infected')
                    FormatFigure;
                end
                
                % Get how long it will take to decrease the prevalence given number of times
                if I_TTT(end)<I_TTT(1)/decrease
                    % If it got low enough, get the first time it went below
                    % Add 12 for the first 12 months of TCT
                    tElimTTT(aux1) = min(t(I_TTT<I_TTT(1)/decrease))+12;
                else
                    % It did not get low enough, add 2 extra years
                    tElimTTT(aux1) = 480+12+24;
                end
                
                % Get R0 with the treatment
                R0treatment(aux1) = R0;
                tEst(aux1) = -log(10)/log(R0treatment(aux1));
            end
            
            % We have just finished the analysis, save the results so that
            % we do not have to rerun again
            % save it only if the analyses was run with many samples
            if Samples > 20
                save(savedFileName,'R0Un','paramsUn', 'SolsUn', 'tElimTCT', 'tElimTTT', 'R0treatment', 'tEst', 'error' )
            end
            
        end
        
        % Process and display the results of the uncertainty analysis
        disp(['Time to eliminate by TCT with base values ' num2str(tElimTCT(1)/12) ' years'])
        disp(['Probability to eliminate in 5 years ' num2str(mean((tElimTCT<60)))])
        disp(['Probability to eliminate in 10 years ' num2str(mean((tElimTCT<120)))])
        
        
        disp(['Time to eliminate by TTT with base values ' num2str(tElimTTT(1)/12) ' years'])
        disp(['Probability to eliminate in 10 years ' num2str(mean((tElimTTT<120)))])
        disp(['Probability to eliminate in 15 years ' num2str(mean((tElimTTT<15*12)))])
        disp(['Probability to eliminate in 20 years ' num2str(mean((tElimTTT<240)))])
        disp(['Probability to eliminate in 25 years ' num2str(mean((tElimTTT<300)))])
        
        % Investigate what is going on with sigma (nothing significant
        % found)
        %         for aux = 1:numParams
        %             figure
        %             sc = scatter(paramsUn(:,4), paramsUn(:,aux), ...
        %                 'MarkerFaceColor', [0.5 0.5 0.5], ...
        %                 'MarkerEdgeColor', [0.5 0.5 0.5]);
        %             alpha(sc, 1/Samples)
        %             sc.SizeData = 5;
        %             xlabel(Labels{4})
        %             ylabel(Labels{aux})
        %             % Set limits on y axis based on parameter limits
        %             ylim([Params{aux,4} Params{aux,5}])
        %             FormatFigure
        %             SaveFigure(['Correlation_sigma_' Names{aux}])
        %
        %
        %         end
        
        %         % See how well we can estimate time to elimination from R0
        %         % (NOT AS WELL AS WE HOPED FOR)
        %         figure
        %         sc = scatter(tElimTTT/12, tEst, ...
        %             'MarkerFaceColor', [0.5 0.5 0.5], ...
        %             'MarkerEdgeColor', [0.5 0.5 0.5]);
        %         alpha(sc, 1/Samples)
        %         sc.SizeData = 5;
        %         xlabel('Time to elimination (in years)')
        %         ylabel('$-(\log(R_0))^{-1}$')
        %         FormatFigure
        %         SaveFigure('Correlation between times')
        
        figure
        %edges = [10:1:31];
        % histogram(tElimTTT/12, edges, 'Normalization', 'probability');
        histogram(tElimTTT/12,'Normalization', 'probability');
        xlabel('Time of TTT application (in years)')
        ylabel('Probability of elimination')
        FormatFigure
        SaveFigure(['tElimTTT_hist_Marino'])
        
        % Plot the bars with PRCC indices
        prcc = PRCClocal(paramsUn, tElimTTT', 'tElimTTT')
        
        figure
        histogram(tElimTCT/12,  'Normalization', 'probability');
        xlabel('Time of TCT application (in years)')
        ylabel('Probability of elimination')
        FormatFigure
        SaveFigure(['tElimTCT_hist_Marino'])
        
        % Plot the bars with PRCC indices
        prcc = PRCClocal(paramsUn, tElimTCT', 'tElimTCT')
        
        pause
        
        
        figure
        histogram(R0treatment,  'Normalization', 'probability');
        xlabel('$R_0$ with TTT treatment')
        ylabel('Probability')
        FormatFigure
        SaveFigure(['R0Treatment_hist_Marino'])
        
        % Plot the bars with PRCC indices
        prcc = PRCClocal(paramsUn, R0treatment', 'R0treatment')
        
        
        
        figure
        histogram(R0Un,  'Normalization', 'probability');
        xlabel('$R_0$ prior treatment')
        ylabel('Probability')
        FormatFigure
        SaveFigure(['R0_hist_Marino'])
        
        % Plot the bars with PRCC indices
        prcc = PRCClocal(paramsUn, R0Un', 'R0')
        
        %         figure
        %         histogram(error,  'Normalization', 'probability');
        %         xlabel('Fit error')
        %         ylabel('Probability')
        %         FormatFigure
        %         SaveFigure(['error_hist_Marino'])
        %
        %         % Plot the bars with PRCC indices
        %         prcc = PRCClocal(paramsUn, error', 'error')
        %
        
        
        
        
        
        
        
        %         % Plot correlation between sigma and compartments
        %         for i1 = 1:numComps
        %             figure
        %             sc = scatter(paramsUn(:,4), SolsUn(:,i1),  ...
        %                         'MarkerFaceColor', [0.5 0.5 0.5], ...
        %                         'MarkerEdgeColor', [0.5 0.5 0.5]);
        %             alpha(sc, 1/Samples)
        %             xlabel(Labels{4})
        %             ylabel(Comps{i1,2})
        %
        %             FormatFigure
        %             SaveFigure(['Correlations_sigma_' Comps{i1,1}])
        %
        %             cor = corrcoef(paramsUn(:,4),SolsUn(:,i1))
        %             disp(['Correlation between sigma and ' Comps{i1,1} ' is ' num2str(cor(1,2))])
        %         end
        
        % Plot correlations between compartments
        for i1 = 1:numComps
            for i2 = 1:numComps
                figure
                if i1==i2
                    histogram(SolsUn(:,i1),  'Normalization', 'probability');
                    xlabel(Comps{i1,2})
                    ylabel('Probability')
                else
                    sc = scatter(SolsUn(:,i1), SolsUn(:,i2),  ...
                        'MarkerFaceColor', [0.5 0.5 0.5], ...
                        'MarkerEdgeColor', [0.5 0.5 0.5]);
                    alpha(sc, 1/Samples)
                    xlabel(Comps{i1,2})
                    ylabel(Comps{i2,2})
                end
                FormatFigure
                SaveFigure(['Correlations' num2str(i1) num2str(i2)])
            end
        end
        
        % Plot parameter distribution
        for i1 = 1:numParams
            figure
            histogram(paramsUn(:,i1), 'Normalization', 'probability');
            xlabel(Labels{i1})
            ylabel('Probability')
            FormatFigure
            SaveFigure([Names{i1} '_hist'])
        end
    end


    function plotR0onTau
        % Plots how R0 depends on tau and graphically solves the equation
        % R0=1
        
        % Store what tau was
        tauOrig = tauM;
        tauMax = 0.005;
        % Reset tau for plotting
        tauM = linspace(0,tauMax, 101);
        
        % Calculate R0 during TTT
        TTT; % Sets probabilities c to TTT setting
        R0TTT = R0; % Stores the R0 values
        tau0TTT = solve_for_tau; % Solves for tau0 during TTT
        dispTTT = ['Tau for TTT is ', num2str(sym2poly(tau0TTT)), ' which corresponds to a round of TTT every ', num2str(1/(12*sym2poly(tau0TTT))), ' years.'];
        disp(dispTTT)
        
        % Calculate R0 during TCT
        TCT; % Sets probabilities c to TTT setting
        R0TCT = R0; % Stores the R0 values
        tau0TCT = solve_for_tau; % Solves for tau0 during TCT
        dispTCT = ['Tau for TCT is ', num2str(sym2poly(tau0TCT)), ' which corresponds to a round of TCT every ', num2str(1/(12*sym2poly(tau0TCT))), ' years.'];
        disp(dispTCT)
        
        % Plot R0 on tau
        figure
        pTTT = plot(tauM, R0TTT, 'k--', 'DisplayName', 'TTT');
        hold on
        pTCT = plot(tauM, R0TCT, 'k-', 'DisplayName', 'TCT');
        % Plot horizontal line where R0 is 1
        plot(tauM, ones(length(tauM),1), ':', 'color', [0.5 0.5 0.5])
        % Plot vertical lines where tau = tau 0
        plot([tau0TTT, tau0TTT], [0,1], ':', 'color', [0.5 0.5 0.5])
        plot([tau0TCT, tau0TCT], [0,1], ':', 'color', [0.5 0.5 0.5])
        
        % Set x and y labels
        xlabel('Treatment rate $\tau$')
        ylabel('Basic reproduction number')
        
        % Show the legend
        lgd = legend([pTTT pTCT], 'Location', 'NorthEast', 'Interpreter', 'LaTeX');
        
        % Specify x axis labels
        xticks([0 sym2poly(tau0TCT) 0.001 sym2poly(tau0TTT) 0.002 0.003  0.004 0.005 ])
        xticklabels({'','$\tau_{0,TCT}$', '', '$\tau_{0,TTT}$',  '', '0.003', '0.004', '0.005'})
        
        yticks([0 0.5 1 1.253 1.5])
        xlim([min(tauM) max(tauM)])
        
        % Format and save the figure
        FormatFigure
        SaveFigure('R0ontau')
        
        % Reset parameters back to what they were
        ResetParameters
        CalculateAuxiliaryVariables
    end


    function check_TTT
        % This checks whether analytical solutions for TTT are correct
        
        tauM = 0; % sets treatment to 0
        CalculateAuxiliaryVariables % this sets all treatment rates appropriate values
        
        %Defining constants for simplification
        r1 = (lambda1*(rho2+mu))/(rho1+mu);
        r2= pY1L1*rho1+pY1Y2*(rho1+mu);
        r3= (rho2+mu)*(lambda2+mu)-pY2L2*lambda2*pL2Y2*rho2;
        r4 = (sigma + mu)*(lambda1 + mu);
        
        r5 = r3*r4*R0;
        t2= (rho2+mu)*(sigma+mu)*cY1*cY2;
        t1 = (rho2+mu)*cY2*r4+(sigma+mu)*cY1*r3-beta*sigma*(rho2+mu)*cY2;
        t0=r3*r4-r5;
        
        %Calculate Tau critical
        display((-t1+sqrt(t1^2-4*t2*t0))/(2*t2));
        %display((-t1-sqrt(t1^2-4*t2*t0))/(2*t2));
        
        % Reset the parameters
        ResetParameters
        CalculateAuxiliaryVariables
        
    end

    function validation
        % Fits real data to the curves and produces the figure
        ResetParameters
        %calibration
        CalculateAuxiliaryVariables
        % Set the data
        % Number of cases of active yaws (Table 1 from
        % https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5920722/ )
        month = [0 6 12 18 24 30 36 42]';
        cases = [323/13490 44/13166 34/13204 38/15977 24/11792 52/14935 53/14765 63/13601]';
        latent = [18.9 14.7 7.3 10.8 5.4 7.9 4.7 6.6]'/100;
        
        % Plot the results
        % It will be done in two figures at the same time, one for active
        % and one for latent
        
        % Start the active infections
        activeFigure = figure;
        axActive  = gca;
        
        % plot active infections real data
        plot(axActive, month,cases, '.', 'MarkerSize', 15);
        hold on
        
        % Start the latent infections
        latentFigure = figure;
        axLatent = gca;
        plot(axLatent, month-3,latent, '.', 'MarkerSize', 15);
        hold on
        
        % Plot our model prediction
        % Plots the solutions of the ODE system (specified by the equations listed in the function getEquations)
        
        % Specify the time of interest, first x months of total community
        % treatment
        tspan = [0 6];
        
        % Specify initial conditions as endemic equilibrium
        InitConditions = getAnalyticalSolutions/N;
        
        % One round of TCT, i.e. the rate is 1/time
        tauM = 1/max(tspan);
        TCT;
        
        
        % Numerically solve the ODE.
        % ode45 is the most universal matlab solver; use another one only if this one causes some issues (unlikely)
        % The input ODEmodel is a function that has to be defined separately
        [t,sols] = ode45(@(t,sols) ODEmodel(t,sols), tspan, InitConditions);
        
        % Unpack the solutions (give it human names)
        % sols is an array: the first component is time (stored in t), the second component are the compartments.
        for comp = 1:numComps
            % Evaluate the value of the compartment
            eval([Comps{comp} ' =  sols(:,comp) ;']);
            % The above is basically like writing S = sols(:,1), i.e. the first component of the ODE solution
            % but for all compartments automatically and in one line
        end
        
        % Plot the active infections in the proper axis (figure)
        plot(axActive, t, (Y1+Y2)./(S+E+Y1+Y2+Y3+L1+L2), 'k');
        
        % Plot the latent infections in the proper axis (figure)
        plot(axLatent, t,(L1+L2)./(S+E+L1+L2), 'k')
        
        % Specify the time of interest, up to 24 months of TTT
        tspan = [max(tspan) 24];
        
        % Start where we ended
        InitConditions = sols(end,:);
        
        % 3 rounds of TTT plus additional non-strategic treatment between
        % rounds
        tauM = 1/6; %tauM = 1/3;
        TTT;
        
        
        % Numerically solve the ODE.
        % ode45 is the most universal matlab solver; use another one only if this one causes some issues (unlikely)
        % The input ODEmodel is a function that has to be defined separately
        [t,sols] = ode45(@(t,sols) ODEmodel(t,sols), tspan, InitConditions);
        
        % Unpack the solutions (give it human names)
        % sols is an array: the first component is time (stored in t), the second component are the compartments.
        for comp = 1:numComps
            % Evaluate the value of the compartment
            eval([Comps{comp} ' =  sols(:,comp) ;']);
            % The above is basically like writing S = sols(:,1), i.e. the first component of the ODE solution
            % but for all compartments automatically and in one line
        end
        
        plot(axActive, t, (Y1+Y2)./(S+E+Y1+Y2+Y3+L1+L2), 'k--');
        plot(axLatent, t, (L1+L2)./(S+E+L1+L2), 'k--');
        
        % Specify the time of interest
        tspan = [24 42];
        
        % Start where we ended
        InitConditions = sols(end,:);
        
        %Set treatment rate - Non-strategic treatment
        tauM = 1/24;
        TTT;
        
        %Numerically solve the ODE.
        %ode45 is the most universal matlab solver; use another one only if this one causes some issues (unlikely)
        %The input ODEmodel is a function that has to be defined separately
        [t,sols] = ode45(@(t,sols) ODEmodel(t,sols), tspan, InitConditions);
        
        %Unpack the solutions (give it human names)
        %sols is an array: the first component is time (stored in t), the second component are the compartments.
        for comp = 1:numComps
            %Evaluate the value of the compartment
            eval([Comps{comp} ' =  sols(:,comp) ;']);
            %The above is basically like writing S = sols(:,1), i.e. the first component of the ODE solution
            %but for all compartments automatically and in one line
        end
        
        plot(axActive, t, (Y1+Y2)./(S+E+Y1+Y2+Y3+L1+L2), 'k:');
        plot(axLatent, t,(L1+L2)./(S+E+L1+L2), 'k:');
        
        %Set axis ticks, labels, and bounds
        xticks(axLatent, month)
        xlabel(axLatent, 'Month')
        ylabel(axLatent, 'Prevalence of latent yaws')
        xlim([-6,42])
        ylim([0, 0.2])
        
        % Format the figure
        FormatFigure(latentFigure)
        % Save the figure
        SaveFigure('Validation_Latent_Yaws_TCT_TTT_lowerTTT')
        close(latentFigure)
        
        %Set axis ticks, labels, and bounds
        xticks(axActive, month)
        xlabel(axActive, 'Month')
        ylabel(axActive, 'Prevalence of active yaws')
        xlim([0,42])
        ylim([0, 0.03])
        
        % Format the figure
        FormatFigure(activeFigure)
        % Save the figure
        SaveFigure('Validation_Active_Yaws_TCT_TTT_lowerTTT')
        
        % Display output for the caption of the figure
        disp('Validation output to the caption')
        tauM = 0;
        disp(['$R_0$ Predicted by the model is ' num2str(R0)]);
    end

    function validation2
        % Fits real data to the curves and produces the figure with some
        % sort of uncertainty
        ResetParameters
        %calibration
        CalculateAuxiliaryVariables
        % Set the data
        % Number of cases of active yaws (Table 1 from
        % https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5920722/ )
        month = [0 6 12 18 24 30 36 42]';
        cases = [323/13490 44/13166 34/13204 38/15977 24/11792 52/14935 53/14765 63/13601]';
        latent = [18.9 14.7 7.3 10.8 5.4 7.9 4.7 6.6]'/100;
        
        myColor = [0 0 0]; % Sets the original color as black
        
        % Plot the results
        % It will be done in two figures at the same time, one for active
        % and one for latent
        
        % Start the active infections
        activeFigure = figure;
        axActive  = gca;
        
        % plot active infections real data
        plot(axActive, month,cases, '.', 'MarkerSize', 15);
        hold on
        
        % Start the latent infections
        latentFigure = figure;
        axLatent = gca;
        plot(axLatent, month-3,latent, '.', 'MarkerSize', 15);
        hold on
        
        tauM = 0;
        % Get init conditions based on calibrated model (i.e. close to
        % current true data)
        ResetParameters;
        CalculateAuxiliaryVariables;
        InitConditions  = getAnalyticalSolutions;
        
        % For the number of experiments, generate parameters randomly and
        % see what the outcomes are
        for aux1 = 1:Samples+1
            disp(['Running simulation ' num2str(aux1) ' of ' num2str(Samples) ])
            if aux1 == 1
                % Use the base values for the first run
                % We will just set R0 as it is and then the values will not
                % have to be generated again
                R0Un(aux1) = R0;
                error(aux1) = 0;
            else
                % Tentatively set R0 to 0 so that we can discard parameters
                % that produce "bad" R0 and too large error
                R0Un(aux1) = 0;
                error(aux1) = 1;
            end
            
            % Keep assigning parameter values randomly until 3>R0>1 and
            % fitting error is small
            while error(aux1)>0.9 %(R0Un(aux1)<=1)||(R0Un(aux1)>3)||
                % Assign parameter values at random (but within the bounds)
                for aux = 1:numParams
                    eval([Names{aux} '=' num2str((Params{aux,5}-Params{aux,4})*rand+Params{aux,4}) ';'])
                end
                % Set tau = 0 for  no treatment
                tauM = 0;
                
                % Calibrate to find beta and qY1 and qY2, store fittting
                % error for next round of "while cycle"
                error(aux1) = calibration;
                % Calculate the auxiliary variables
                CalculateAuxiliaryVariables
                
                % Evaluate R0
                R0Un(aux1) = R0;
            end
            % Now we have random parameters such that 1<R0<3 and good fit
            % to data
            
            % For the last run, resume main parameters and draw black
            if aux1 == Samples+1
                ResetParameters;
                CalculateAuxiliaryVariables;
                InitConditions  = getAnalyticalSolutions;
                myColor = [0 0 0];
            end
            
            % Plot our model prediction
            % Plots the solutions of the ODE system (specified by the equations listed in the function getEquations)
            
            % Specify the time of interest, first x months of total community
            % treatment
            tspan = [0 6];
            
            % Specify initial conditions as endemic equilibrium
            InitConditions = getAnalyticalSolutions/N;
            
            % One round of TCT, i.e. the rate is 1/time
            tauM = 1/max(tspan);
            TCT;
            
            
            % Numerically solve the ODE.
            % ode45 is the most universal matlab solver; use another one only if this one causes some issues (unlikely)
            % The input ODEmodel is a function that has to be defined separately
            [t,sols] = ode45(@(t,sols) ODEmodel(t,sols), tspan, InitConditions);
            
            % Unpack the solutions (give it human names)
            % sols is an array: the first component is time (stored in t), the second component are the compartments.
            for comp = 1:numComps
                % Evaluate the value of the compartment
                eval([Comps{comp} ' =  sols(:,comp) ;']);
                % The above is basically like writing S = sols(:,1), i.e. the first component of the ODE solution
                % but for all compartments automatically and in one line
            end
            
            % Plot the active infections in the proper axis (figure)
            plot(axActive, t, (Y1+Y2)./(S+E+Y1+Y2+Y3+L1+L2), '-', 'color', myColor);
            
            % Plot the latent infections in the proper axis (figure)
            plot(axLatent, t,(L1+L2)./(S+E+L1+L2), '-', 'color', myColor)
            
            % Specify the time of interest, up to 24 months of TTT
            tspan = [max(tspan) 24];
            
            % Start where we ended
            InitConditions = sols(end,:);
            
            % 3 rounds of TTT plus additional non-strategic treatment between
            % rounds
            tauM = 1/6; %tauM = 1/3;
            TTT;
            
            
            % Numerically solve the ODE.
            % ode45 is the most universal matlab solver; use another one only if this one causes some issues (unlikely)
            % The input ODEmodel is a function that has to be defined separately
            [t,sols] = ode45(@(t,sols) ODEmodel(t,sols), tspan, InitConditions);
            
            % Unpack the solutions (give it human names)
            % sols is an array: the first component is time (stored in t), the second component are the compartments.
            for comp = 1:numComps
                % Evaluate the value of the compartment
                eval([Comps{comp} ' =  sols(:,comp) ;']);
                % The above is basically like writing S = sols(:,1), i.e. the first component of the ODE solution
                % but for all compartments automatically and in one line
            end
            
            plot(axActive, t, (Y1+Y2)./(S+E+Y1+Y2+Y3+L1+L2), '--', 'color', myColor);
            plot(axLatent, t, (L1+L2)./(S+E+L1+L2), '--', 'color', myColor);
            
            % Specify the time of interest
            tspan = [24 42];
            
            % Start where we ended
            InitConditions = sols(end,:);
            
            %Set treatment rate - Non-strategic treatment
            tauM = 1/24;
            TTT;
            
            %Numerically solve the ODE.
            %ode45 is the most universal matlab solver; use another one only if this one causes some issues (unlikely)
            %The input ODEmodel is a function that has to be defined separately
            [t,sols] = ode45(@(t,sols) ODEmodel(t,sols), tspan, InitConditions);
            
            %Unpack the solutions (give it human names)
            %sols is an array: the first component is time (stored in t), the second component are the compartments.
            for comp = 1:numComps
                %Evaluate the value of the compartment
                eval([Comps{comp} ' =  sols(:,comp) ;']);
                %The above is basically like writing S = sols(:,1), i.e. the first component of the ODE solution
                %but for all compartments automatically and in one line
            end
            
            plot(axActive, t, (Y1+Y2)./(S+E+Y1+Y2+Y3+L1+L2), ':', 'color', myColor);
            plot(axLatent, t,(L1+L2)./(S+E+L1+L2), ':', 'color', myColor);
            
            % Set the next colors as gray
            myColor = [0.9 0.9 0.9 0.2];
        end
        
        % Redraw the dots
        plot(axActive, month,cases, '.', 'MarkerSize', 15, 'color', [0 0 1]);
        plot(axLatent, month-3,latent, '.', 'MarkerSize', 15, 'color', [0 0 1]);
        %Set axis ticks, labels, and bounds
        xticks(axLatent, month)
        xlabel(axLatent, 'Month')
        ylabel(axLatent, 'Prevalence of latent yaws')
        xlim([-6,42])
        ylim([0, 0.2])
        
        % Format the figure
        FormatFigure(latentFigure)
        % Save the figure
        SaveFigure('Validation_Latent_Yaws_TCT_TTT_lowerTTT')
        close(latentFigure)
        
        %Set axis ticks, labels, and bounds
        xticks(axActive, month)
        xlabel(axActive, 'Month')
        ylabel(axActive, 'Prevalence of active yaws')
        xlim([0,42])
        ylim([0, 0.03])
        
        % Format the figure
        FormatFigure(activeFigure)
        % Save the figure
        SaveFigure('Validation_Active_Yaws_TCT_TTT_lowerTTT')
        
        % Display output for the caption of the figure
        disp('Validation output to the caption')
        tauM = 0;
        disp(['$R_0$ Predicted by the model is ' num2str(R0)]);
    end



% Part 6 -- Essential functions that do not depend on the model

    function ResetParameters
        % Assigns the parameter base values
        for aux = 1:numParams
            % Go through all the parameters
            % Evaluate the command that assigns the parameter with a given name the proper base value
            eval([Names{aux} '=' num2str(BaseValues{aux}) ';'])
        end
    end

    function CalculateAuxiliaryVariables
        % Calculates the auxiliary variables
        for aux = 1:numAuxVars
            eval([AuxVars{aux,1} '=' AuxVars{aux,3} ';']);
        end
    end



    function error = CheckSolution(solution)
        % Checks whether the supplied solutions indeed solves the system of equations = 0
        
        % Takes the array "solution" in the order of Compartments
        % Evaluates the equations with this potential solution
        % Returns the values of the equations
        % If the values are near 0, it is likely that the supplied vector
        % is a solution of the equations
        
        % Assign the possible solutions to the compartments
        for comp = 1:numComps
            eval([Comps{comp} ' = solution(comp);']);
        end
        
        % Evaluate the equations with these values
        error = getEquations;
    end



    function sols = numSolve
        % Numerically solves the system of equilibria equations
        
        % Initialize the list of unknowns for the equations
        dummm = sym('dummm');
        unknowns = [dummm]; % This makes vars the list of appropriate type. FIX IT WITH CLEANER PROGRAMMING
        
        % Set up the unknowns as symbolic and put them into variables
        for comp = 1:numComps
            % Make the compartment a symbolic variable
            eval([Comps{comp} '= sym('' ' Comps{comp} ' '');']);
            % The above is basically like writing S = sym('S') but for all
            % compartments automatically and in one line
            
            % Set the variable as the unknown in the equations
            eval(['unknowns(comp) = ' Comps{comp} ';' ]);
        end
        
        % Get the system of equations
        EqsRHS = getEquations;
        
        % Assign the system of equations
        for eq = 1:length(EqsRHS)
            eqns(eq) = EqsRHS{eq} == 0;
        end
        
        % Numerically solve the equations and store it in "struct" variable
        sols = vpasolve(eqns, unknowns);
        
        % Interpret the solutions
        numSols = []; % initialize the variable for number of solutions
        eval(['numSols = size(sols.' Comps{1} ');']) % gets the number of solutions
        
        disp_Diary(['There are ' num2str(numSols) ' sets of solutions']);
        
        numericalSolutions = []; % Just initialize the array storing all solutions
        % numericalSolutions(i,j) will be the ith set, compartment j
        
        for solSet = 1:numSols
            disp_Diary(['In the order ' Comps{:} ', the solution set number ' num2str(solSet) ' is'])
            
            % Actually get the solutions from the structured array
            for comp = 1:numComps
                eval(['numericalSolutions(solSet, comp) = sols.' Comps{comp} '(solSet);']);
            end
            % Print out the solutions
            disp_Diary(num2str(numericalSolutions(solSet,:) ));
        end
        
        
    end

    function dXdt = ODEmodel(t, X)
        % Returns the equations of the ODE model for the plotting of the
        % ODEs
        
        % Unpack the variable X
        for comp = 1:numComps
            % Evaluate the value of each compartment
            eval([Comps{comp} '=  X(comp) ;']);
            % The above is basically like writing S = X(1); etc but for all
            % compartments automatically and in one line
        end
        
        % Plug the compartment values into the equations to get the value of the derivatives
        % Have to convert the cell array to a (column) vector
        dXdt = cell2mat(getEquations);
    end


    function output = checkConsistency(num_of_tests)
        % Checks for errors
        
        % will generate num_of_tests sets of parameters (more or less
        % randomly around the predefined base values) and displays warning
        % if anything suspicious is found.
        
        % Checks include
        % - confirming the accuracy of analytical solutions
        % - checking that R0>1 corresponds to endemic equilibrium (by
        % comparing it to numerical solutions
        % ADD MORE AS NEEDED
        
        disp('Checking for consistency. It may take some time')
        
        % Initialize the number of times R0>1
        Endemic = 0;
        % Initialize the number of time endemic is unstable
        Unstable = 0;
        
        for  dummy_variable= 1:num_of_tests
            disp(['Test ' num2str(dummy_variable) ' of ' num2str(num_of_tests) '.' ])
            
            % Assign parameter values at random (but within the specified range)
            for aux = 1:numParams
                % Go through all the parameters
                % Evaluate the command that assigns the parameter with a given name the proper base value
                %eval([Names{aux} '=' num2str(BaseValues{aux}*2*rand) ';'])
                eval([Names{aux} '=' num2str((Params{aux,5}-Params{aux,4})*rand+Params{aux,4}) ';'])
            end
            
            % Calculate the auxiliary variables
            CalculateAuxiliaryVariables
            
            % Calculate Analytical Solutions
            ansols = getAnalyticalSolutions
            
            % Check whether analytical solutions are correct
            tol = 10^(-10);  % tolerance level for correctness of analytical solutions
            if max(abs(cell2mat(CheckSolution(ansols))))>tol
                % if the max error for analytical solutions is bigger then
                % tolerance, display the warning
                disp('Analytical solutions may be wrong.')
                displayParameters
            end
            
            % Check whether R0<1 corresponds to no endemic equilibrium
            % Get R0
            Reff = R0;
            
            % Update the number of times R0>1
            if Reff>1
                Endemic = Endemic + 1;
            end
            
            % Get numerical solutions
            numsols = numSolve;
            
            if (((Reff<1) && (any(numsols.E>0))) && (epsilon ==0))
                % if R0<1 but there is an endemic equilibrium (positive exposed)
                disp('R0<1 but endemic equilibrium exists')
                displayParameters
            end
            
            if ((Reff>1) && (all(numsols.E<=0)))
                % if R0>1 but there is no  endemic equilibrium (no positive exposed)
                disp('R0>1 but endemic equilibrium does not exist')
                displayParameters
            end
            
            n_of_numsols = [];
            eval(['n_of_numsols = size(numsols.' Comps{1} ',1);']);
            
            % Get the index where numerical and analytical should agree
            solIdx = 0; % Initialize the index
            for aux = 1:n_of_numsols % For every numerical solution
                if abs(numsols.E(aux) - ansols(2))<tol  % if numerical value for E agrees with analytical
                    % Set the index
                    solIdx = aux;
                end
            end
            
            % Check whether the solutions agree for all compartments
            if solIdx == 0 % did not agree even for one compartment
                disp('Numerical and analytical solutions do not agree, Idx = 0')
                displayParameters
            else
                % Initialize the agreement array
                agree = [];
                for comp = 1:numComps
                    % Check the agreement at compartment comp
                    eval(['agree(comp) = abs(numsols.' Comps{comp} '(solIdx)-ansols(comp))<tol;'])
                end
                if ~all(agree)  % If there was any disagreement
                    disp('Numerical and analytical solutions do not agree')
                    displayParameters
                end
            end
            
            % Check the global stability of the endemic equilibria
            if Reff>1 % only if the endemic equilibrium is possible
                disp('Checking global stability of endemic equilibria')
                % Set initial conditions
                InitConditions  = [1 1 0 0 0 0 0];
                % Set thenumber of months to run the ODEs
                tspan = [0, 12000];
                % Numerically solve the ODE
                [t,y] = ode45(@(t,y) ODEmodel(t,y), tspan, InitConditions);
                % Get the difference betweennumerical and analytical
                % solution
                difference = max(abs(ansols - y(end,:)));
                if difference >10^(-3)
                    disp('Equilibrium may not be stable')
                    Unstable = Unstable +1;
                    displayParameters
                end
                
            end
        end
        disp(['Out of ' num2str(Endemic) ' endemic equilibria, ' num2str(Unstable) ' were unstable'])
        
        disp('Consistency check complete. If no warning was displayed, no inconsistency was found')
        
        
        % Reset parameters back to base values
        ResetParameters
        CalculateAuxiliaryVariables
    end



    function prcc=PRCClocal(values, Y, FileName);
        % PRCC procedure
        % input:  Y, response function that we investigate how it depends
        % on the parameters
        % output - the prcc coefficients
        
        % Based on Marino, S., Hogue, I. B., Ray, C. J., and Kirschner, D. E. (2008).
        % A methodology for performing global uncertainty and sensitivity analysis in systems
        % biology. Journal of Theoretical Biology, 254(1):178–196
        % Part of the code for PRCC taken from
        % http://malthus.micro.med.umich.edu/lab/usanalysis.html
        
        % Initiate the LHS matrix that will carry all the parameters
        LHSmatrix= values;
        
        for i=1:numParams  % Loop for the parameters
            
            % Remove the ith parameter
            LHStemp=LHSmatrix;
            LHStemp(:,i)=[];
            Ztemp=LHStemp;
            LHStemp=[];
            
            % Get the correlation coefficient
            [rho_aux,~]=partialcorr([LHSmatrix(:,i),Y],Ztemp,'type','Spearman');
            
            % Record the correlation coefficient
            prcc(1,i)=rho_aux(1,2);
            
        end
        
        % Sort the correlation coefficients for better plotting
        [sortedPRCC,I] = sort(prcc, 'descend');
        
        %             % Plot PRCC coefficients as a horizontal bar plot
        %             figure;
        %             h=barh(sortedPRCC);
        %             % Set the y limits by the number of parameters
        %             ylim([0.5 numParams+0.5]);
        %             % Ticks at integers and labeled by the parameter names/labels
        %             ax = gca;
        %             ax.YTick = 1:numParams;
        %             ax.YTickLabel = Labels(I);
        %             ylabel('Parameters');
        %             % Set the x axis between -1 and 1
        %             xlim([-1 1]);
        %             xlabel('PRCC');
        %             % Make the color gray
        %             h.FaceColor = [0.5 0.5 0.5];
        %             % Format the figure
        %             FormatFigure;
        %             SaveFigure([FileName '_bar_Marino'])
        %
        % PLot only the important coefficients
        figure
        % Determine which ones are important
        importantPRCC = sortedPRCC(abs(sortedPRCC)>0.01);
        h=barh(importantPRCC);
        ylim([0.5 length(importantPRCC)+0.5]);
        % Ticks at integers and labeled by the parameter names/labels
        ax = gca;
        ax.YTick = 1:length(importantPRCC);
        ax.YTickLabel = Labels(I(abs(sortedPRCC)>0.01));
        ylabel('Parameters');
        % Set the x axis between -1 and 1
        xlim([-1 1]);
        xlabel('PRCC');
        % Make the color gray
        h.FaceColor = [0.5 0.5 0.5];
        
        FormatFigure;
        SaveFigure([FileName '_bar_Marino_important'])
    end


    function out = displayParameters
        % Displays parameter values for cross-checking
        disp('Parameter values are')
        for aux = 1:numParams
            % Go through all the parameters and display their values
            eval([Names{aux}])
        end
        pause
        
    end

    function globalStability(numOfRepetition)
        % Checks whether the endemic equilibrium is globally stable
        % Defines a Lyapunov function and checks that the derivative is
        % negative by randomly generating a bunch of compartments and
        % pluggin it in
        
        ResetParameters;
        
        % Store what the original value of tauM is
        tauOriginal = tauM;
        
        % Check whether there is a solution or not
        tauM = 0; % sets treatment to 0
        
        TCT(0);
        
        
        % Get analytical solutions for endemic equilibrium
        asols = getAnalyticalSolutions;
        % Unpack the solutions
        aS = asols(1);
        aE = asols(2);
        aY1 = asols(3);
        aY2 = asols(4);
        aY3 = asols(5);
        aL1 = asols(6);
        aL2 = asols(7);
        
        for aux = 1:numOfRepetition
            % Randomly assign values to all compartments
            for comp = 1:numComps % for all compartments
                % Assign a value to the compartment
                eval([Comps{comp,1} '=' num2str(rand(1)*10) ';'])
            end
            % Evaluate the derivatices
            ders = cell2mat(getEquations);
            % Unpack the derivatives
            for comp = 1:numComps % for all compartments
                % Assign a value to the compartment
                eval(['d' Comps{comp,1} '=' num2str(ders(comp)) ';'])
            end
            % Evaluate the derivative to Lyapunov function
            dLyap(aux) = (1-aS/S)*dS + (1-aE/E)*dE +(1-aY1/Y1)*dY1 +...
                (1-aY2/Y2)*dY2 +(1-aY3/Y3)*dY3 +(1-aL1/L1)*dL1 + ...
                (1-aL2/L2)*dL2;

        end
        % Check that no value is positive
        if max(dLyap)>0
            disp('SOMETHING IS WRING WITH GLOBAL STABILITY')
        else
            disp('Endemic equilibrium seems globally stable')
        end
        
        % Reset tauM to whatever it was
        tauM = tauOriginal;
        
    end



%% (Part 7) Auxiliary Functions for easy plotting, etc


    function FormatFigure(h)
        % Formats figures so that all figures are uniformly formatted
        
        if nargin == 0 % if the function is called without an argument, assume we want to format current figure
            h = gcf;
        end
        
        % Get axes handle
        haxes=get(h, 'CurrentAxes');
        
        % Get figure number
        fnbr = get(h,'Number');
        [fig_x,fig_y] = getPosition(fnbr);
        
        % Set dimensions of the figure
        set(h,'Units','Inches','Position',[fig_x, fig_y, 3, 2.5]);
        
        % Set default font size
        set(haxes,'FontSize',10)
        % Set interpreter to LaTeX
        set(0,'defaultTextInterpreter','latex')
        
        % Set the interpreter for latex
        set(haxes, 'TickLabelInterpreter', 'latex');
        
        % Find all line objects
        hline= findobj(haxes,'Type','line');
        % Set line properties (relatively fat lines)
        set(hline, 'LineWidth'   , 1);
        
        % Set axes properties
        set(haxes, ...
            'Box'         , 'on'     , ...  %puts a box around the axis
            'TickDir'     , 'in'     , ...
            'TickLength'  , [.02 .01] , ...
            'XMinorTick'  , 'off'      , ...
            'YMinorTick'  , 'off'     , ...
            'LineWidth'   , 1);
        
        function [xpos,ypos] = getPosition(n)
            % gets the x y position of the figure to plot
            n = mod(n-1,15)+1;
            ypos = floor((n-1)/5)*5;
            xpos = (mod(n-1,5))*3;
        end
        
    end

    function SaveFigure(filename)
        % Saves a current graphics file to file named filename.pdf
        
        % Get the handle
        h=gcf;
        set(h,'Units','Inches');
        pos = get(h,'Position');
        
        % Cut all the white margins
        set(h,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
        
        % Actually save
        print(h,['Figures/' filename, '.pdf'],'-dpdf','-r0')
        
        % This outputs on the screen what figure was saved into what file
        disp_Diary(['Saved a figure ' num2str(get(gcf,'Number')) ' into a file ' ...
            filename '.pdf'])
    end

    function disp_Diary(string)
        % Displays text on the screen and also logs it in the diary
        diary on
        disp(string)
        diary off
    end

end