%SCRIPT FOR RUNNING THE CRA TOOLBOX WITHOUT USING THE GUI
%The model is written as a Matlab function.

%%%% THE FIRST THING TO DO FOR RUNNING THE SCRIPT IS TO ADD THE FOLDER
%%%% CLASSES AND THE FOLDER Examples TO THE PATH %%%%

%Parameters to be set by the user are:
% 1-model_name: name of the model in .xml format. A struct is created to
% store the following characteristics of an ODE model: nominal values and
% names of the parameters, initial conditions and names of variables, input
% values.
% 2-stop_time: final time point for the model simulation
% 3-step_size: time interval for the vector of time points
% 4- ode_solver: type of ode_solver for simulating the model (possible
% choices are: ode45, ode15s, ode23t, sundials, stochastic, explicit tau
% and implicit tau
% 5- Nr: number of independent realizations to perform 
% 6- LBpi: lower boundary of the Latin Hypercube Sampling for perturbation
% of the model parameter space
% 7- UBpi: upper boundary of the Latin Hypercube Sampling
% 8- Ns: number of samples of the Latin Hypercube (parameters n of the
% lhsdesign function)
% 9- variable_name: variable of the model to set as reference node to be
% measured
% 10- current_func: is the type of evaluation function. Currently, it is
% possible to choose among
% three evaluation functions: area under the curve, maximum value and time
% of maximum for the time behavior of the selected variables. The user can
% also define his evaluation function in a .m file by extending the
% abstract class EvaluationFunction
% 11- tail_size: number of samples to include in the upper and lower tail
% when computing the probability density function of the evaluation
% function
% 12- current_tm: is the method for computing the tails of the pdf of the
% evaluation function. Right now, it is possible to choose between two
% methods: sorted() which sorts the values of the evaluation function and
% selects the first and last samples according to tail_size; tmp_sum()
% computes the tails by selecting the upper and lower quartile. When using
% tmp_sum(), the parameter step_size needs also to be specified. Step_size
% is the step for computing the lower and upper quartile of the pdf in an
% iterative way, i.e. when the upper and lower tails do not have the number of
% samples specified by the user, the threshold is increased of a quantity
% equal to step_size and the calculation of the tails is repeated. 
% The user can also define his own method for the tails computation by
% extending the abstract class TailMethod().
% 13- folder: name of the folder to create where the results are saved

tic;

model_name='ode_covid19_v2';

stop_time=80;
step_size=1;
ode_solver='ode15s';

%time axis for model simulation
time_axis=[0:step_size:stop_time]';

%parameters and initial conditions of the model

N = 882000;         
InitInf=1;

E00 = InitInf;
S0 = N-E00;
x0 = [S0 E00 0 0 0 0 0 0 0];

load('Examples/moda_parametri_inizio14step_LOG_100000_UMBRIA_I2I3DE.mat');
nominal_input_parameters=moda_parametri_inizio14step_LOG;

be = nominal_input_parameters(1);
b0=nominal_input_parameters(2);
b1=nominal_input_parameters(3);
b2=nominal_input_parameters(4);
b3=nominal_input_parameters(5);
        
be=be/(N);
b0=b0/(N);
b1=b1/(N);
b2=b2/(N);
b3=b3/(N);
        
FracSevere = nominal_input_parameters(6);
FracCritical = nominal_input_parameters(7);
FracMild = 1-FracSevere-FracCritical;
FracAsym = nominal_input_parameters(8);
               
IncubPeriod = nominal_input_parameters(9);
DurMildInf = nominal_input_parameters(10);

DurAsym = nominal_input_parameters(11);
   
DurHosp = nominal_input_parameters(12);
TimeICUDeath = nominal_input_parameters(13);
ProbDeath=nominal_input_parameters(14);

PresymPeriod=nominal_input_parameters(15);
s01=nominal_input_parameters(16);
s11=nominal_input_parameters(17);
            
CFR=ProbDeath*FracCritical/100; 
            
a1=1/PresymPeriod; %presymptomatic period of transmission
a0=1/(IncubPeriod-PresymPeriod); %true latent period  

f=FracAsym;
    
g0=1/DurAsym;

g1=(1/DurMildInf)*FracMild;
p1=(1/DurMildInf)-g1;
    
p2=(1/DurHosp)*(FracCritical/(FracSevere+FracCritical));
g2=(1/DurHosp)-p2;
    
u=(1/TimeICUDeath)*(CFR/FracCritical);
    
g3=(1/TimeICUDeath)-u;

nominal_parameters=[be b0 b1 b2 b3 a0 a1 f g0 g1 p1 g2 p2 g3 u]; 
nominal_parameters=ones(1,15);
parameters_name={'be', 'b0', 'b1', 'b2','b3', 'a0', 'a1', 'f', 'g0', 'g1', 'p1', 'g2', 'p2', 'g3','u'};

Tlock1=23;

num_observables=6;
observables_name={'E1','I0','I1','I2','I3','D'};

model=struct('name',model_name,'odesolver',ode_solver,'time',time_axis,'stop',stop_time,'step',step_size,'nominal_parameters',nominal_parameters,'parameters_name',{parameters_name},'num_observables',num_observables,'observables_name',{observables_name},'initial_conditions',x0,'Tlock1',Tlock1,'s01',s01,'s11',s11);

Nr=3;

%define specific boundaries for each parameter
LBpi=[0.001 0.001 0.001 0.001 0.001 1e-5 1e-5 0.01 2 1 1 1 1 1 1];
UBpi=[3 3 3 3 3 0.5 0.4 0.7 8 20 20 10 30 60 4];
Ns=1000;

variable_name='I0';
current_func=Maximum(); %current evaluation function
tail_size=10;  %number of samples for the lower and upper tail
step_size=0.001;  %this parameter needs to be defined only when current_tm=tmp_sum(step_size)
current_tm=sorted();

folder='covid19_I0_results';

%model simulation for each sample of the Latin Hypercube
disp('Starting model simulation with perturbed parameters');
[AllResults,AllPerturbations]=start_simulation_v2_apr2020(model,Nr,LBpi,UBpi,Ns);
disp('All done! Model simulation completed!');

%computation of the MIRI for the chosen evaluation function and model
%variable
disp('Starting computation of the MIRI for each parameter...');
try
  compute_MIRI(model,variable_name,current_func,tail_size,current_tm,Nr,Ns,AllResults,AllPerturbations,folder)
catch ME
    %break  
    return
end
%plot and save probability density function of the evaluation function
disp('Plot of the probability density function of the chosen evaluation function');
plotpdf_evalfunc(folder,variable_name);

%plot and save conditional probability density functions of the parameters 
disp('Plot of the parameter probability density functions');
plotpdf_param(folder,variable_name,model,Nr);

toc;