%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_v3_OK';

stop_time=300;
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
init_italia
load('Examples/ITALIA/italia_moda_nr4.mat');
nominal_parameters=[italia_moda_nr4(1,1:19) 1 1 italia_moda_nr4(1,20:21)];
%nominal_parameters=italia_moda_nr4;

nominal_parameters_name={'be', 'b0', 'b1', 'b2','b3', 'FracSevere', 'FracCritical', 'FracAsym', 'IncubPeriod', 'DurMildInf', 'DurAsym', 'DurHosp', 'TimeICUDeath', 'ProbDeath', 'PresymPeriod','s01','s02','s03','s04','s05','s06','s1','s2'};
s00=[nominal_parameters(16) nominal_parameters(17) nominal_parameters(18) nominal_parameters(19) nominal_parameters(20) nominal_parameters(21)];
s11=nominal_parameters(22);
s22=nominal_parameters(23);

derived_parameters=1;
derived_parameters_name={'R0'};


num_observables=9;
observables_name={'S','E0','E1','I0','I1','I2','I3','R','D'};

model=struct('name',model_name,'odesolver',ode_solver,'time',time_axis,'stop',stop_time,'step',step_size,'nominal_parameters',nominal_parameters,'nominal_parameters_name',{nominal_parameters_name},'derived_parameters',derived_parameters,'derived_parameters_name',{derived_parameters_name},'num_observables',num_observables,'observables_name',{observables_name},'initial_conditions',x0,'Tlock',Tlock,'region_name',region_name);

Nr=10;

%define specific boundaries for each parameter
%LBpi=[0.5 0.25 0.25 0.25 0.25 0.125 0.125 0.25 1 0.4311 0.6988 0.1504 0.9213 0.5 0.2995];
%UBpi=[2 4 4 4 4 32 16 10 2.81 2.8741 6.9881 1.5042 9.2128 2 1.1980];

%LBpi=[0.125 0.25 0.25 0.25 0.25 0.125 0.125 0.25 1 0.5 0.6988 0.5556 0.4 0.25 0.25];
%UBpi=[1.02 4 4 4 4 1.1 1.1 1.5 3 2 1.5 2.7778 4 1.5 1.02];

%LBpi=[0.125 0.25 0.25 0.25 0.25 0.125 0.125 0.25 0.5 0.5 0.6988 0.5556 0.4 0.25 0.25];
%UBpi=[1.02 4 4 4 4 1.1 1.1 1.5 2 2 1.5 2.7778 4 1.5 1.02];


load('Examples/ITALIA/intervalli_figura_bande_90perc_CON_PRCTILE_ITALIA_4Nr.mat')
%LBpi=parameters{1,1}(1,:)./nominal_parameters;
%UBpi=parameters{1,1}(2,:)./nominal_parameters;

LBpi=[parameters{1,1}(1,1:19) 0.1 0.1 parameters{1,1}(1,20:21)]./nominal_parameters;
UBpi=[parameters{1,1}(2,1:19) 0.9 0.9 parameters{1,1}(2,20:21)]./nominal_parameters;

%LBpi=[0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.6396 0.7018 0.5403 0.7539 0.6150 0.6177 0.5640 0.6903 0.6754 0.5958 0.5023 0.4850 0.1 0.1 0.6887 0.5429];
%UBpi=[2 2 2 2 2 1.5 1.38 1.2792 1.2281 1.6210 1.2565 1.5376 3.0885 1.0152 1.2426 1.3509 1.3903 1.5068 1.4549 0.9 0.9 2.4105 2.1716];   

%LBpi=[0.99 0.5 0.5 0.5 0.5 0.5 0.5 0.6396 0.7018 0.5403 0.7539 0.6150 0.6177 0.5640 0.99 0.6754 0.5958 0.5023 0.4850 0.6887 0.5429];
%UBpi=[1.05 2 2 2 2 1.5 1.38 1.2792 1.2281 1.6210 1.2565 1.5376 3.0885 1.0152 1.05 1.3509 1.3903 1.5068 1.4549 2.4105 2.1716];   


Ns=10000;

variable_name='I2';
current_func=Area(); %current evaluation function
tail_size=1000;  %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_area_10Nr_1000sample_time300_ITALIA_4Nr_90perc_predizione';

%model simulation for each sample of the Latin Hypercube
disp('Starting model simulation with perturbed parameters');
[AllResults,AllPerturbations,AllDerivedParam]=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,AllDerivedParam,folder)
catch ME 
    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;