Skip to content

SysBiOThe / CRA-Matlab

Go to a project
Toggle navigation
Toggle navigation pinning
Commit 99e1a80f by Chiara Antonini
Options

update covid19

1 parent ba01dbde
Inline Side-by-side
Showing 150 changed files with 449 additions and 80 deletions
CRA_Updated_Version/CRA_Toolbox_bash/Examples/ode_covid19_v3.m
......@@ -8,13 +8,13 @@ b1=p(3);
b2=p(4);
b3=p(5);
N=10^5;
Norm=10^5;
be=be/(N);
b0=b0/(N);
b1=b1/(N);
b2=b2/(N);
b3=b3/(N);
be=be/(Norm);
b0=b0/(Norm);
b1=b1/(Norm);
b2=b2/(Norm);
b3=b3/(Norm);
FracSevere = p(6);
FracCritical = p(7);
......@@ -54,7 +54,7 @@ g3=(1/TimeICUDeath)-u;
s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3)) * (t<Tlock(4)) + s00(3)*(t>=Tlock(4));
s1 = (t<Tlock(2)) + s11(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s11(2)*(t>=Tlock(3))* (t<Tlock(4)) + s00(3)*(t>=Tlock(4));
s1 = (t<Tlock(2)) + s11(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s11(2)*(t>=Tlock(3))* (t<Tlock(4)) + s11(3)*(t>=Tlock(4));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
%s3 = (t<Tlock(1)) + s33 * (t>=Tlock(1));
......
CRA_Updated_Version/CRA_Toolbox_bash/compute_MIRI.m
......@@ -34,10 +34,16 @@ function compute_MIRI(model,variable_name,current_func,tail_size,current_tm,Nr,N
nominal_parameters=model.nominal_parameters;
%num_parameters=length(nominal_parameters);
nominal_parameters_name=model.nominal_parameters_name;
derived_parameters=model.derived_parameters;
derived_parameters_name=model.derived_parameters_name;
num_parameters=length(nominal_parameters)+length(derived_parameters);
parameters_name=[nominal_parameters_name derived_parameters_name];
if(isfield(model,'derived_parameters'))
derived_parameters=model.derived_parameters;
derived_parameters_name=model.derived_parameters_name;
num_parameters=length(nominal_parameters)+length(derived_parameters);
parameters_name=[nominal_parameters_name derived_parameters_name];
else
num_parameters=length(nominal_parameters);
parameters_name=[nominal_parameters_name];
end
obs_name=model.observables_name;
......@@ -131,16 +137,20 @@ try
%[all_xbin_max,all_xbin_min,all_ks_max,all_ks_min]=obj_MIRI.estimate_parampdfs(handles1.nominal_p, handles1.perturbation, Ncloud, handles.XiMax, handles.XiMin);
[MIRI_param,CloudCondMin, CloudCondMax]=obj_MIRI.MIRI_computation(all_ks_max,all_ks_min,all_xbin_max,all_xbin_min);
obj_MIRI_derived=MIRI();
[all_xbin_max_d,all_xbin_min_d,all_ks_max_d,all_ks_min_d]=obj_MIRI_derived.estimate_parampdfs(derived_parameters, AllDerivedParam{k,1}, XiMax, XiMin);
%[all_xbin_max,all_xbin_min,all_ks_max,all_ks_min]=obj_MIRI.estimate_parampdfs(handles1.nominal_p, handles1.perturbation, Ncloud, handles.XiMax, handles.XiMin);
[MIRI_param_d,CloudCondMin_d, CloudCondMax_d]=obj_MIRI_derived.MIRI_computation(all_ks_max_d,all_ks_min_d,all_xbin_max_d,all_xbin_min_d);
if(isfield(model,'derived_parameters'))
obj_MIRI_derived=MIRI();
[all_xbin_max_d,all_xbin_min_d,all_ks_max_d,all_ks_min_d]=obj_MIRI_derived.estimate_parampdfs(derived_parameters, AllDerivedParam{k,1}, XiMax, XiMin);
%[all_xbin_max,all_xbin_min,all_ks_max,all_ks_min]=obj_MIRI.estimate_parampdfs(handles1.nominal_p, handles1.perturbation, Ncloud, handles.XiMax, handles.XiMin);
[MIRI_param_d,CloudCondMin_d, CloudCondMax_d]=obj_MIRI_derived.MIRI_computation(all_ks_max_d,all_ks_min_d,all_xbin_max_d,all_xbin_min_d);
MIRI_param=[MIRI_param MIRI_param_d];
MIRI_param=[MIRI_param MIRI_param_d];
CloudCondMax=[CloudCondMax CloudCondMax_d];
CloudCondMin=[CloudCondMin CloudCondMin_d];
end
MIRIT(k,:)=MIRI_param;
CloudCondMax=[CloudCondMax CloudCondMax_d];
CloudCondMin=[CloudCondMin CloudCondMin_d];
CloudCondMaxT=[CloudCondMaxT;CloudCondMax];
CloudCondMinT=[CloudCondMinT;CloudCondMin];
mean_CloudMax=mean(CloudCondMaxT);
......@@ -153,11 +163,13 @@ try
ks_param{2,k}(ip,:)=all_ks_min(ip,:);
end
for ip=1:length(derived_parameters)
xbin_param_d{1,k}(ip,:)=all_xbin_max_d(ip,:);
xbin_param_d{2,k}(ip,:)=all_xbin_min_d(ip,:);
ks_param_d{1,k}(ip,:)=all_ks_max_d(ip,:);
ks_param_d{2,k}(ip,:)=all_ks_min_d(ip,:);
if(isfield(model,'derived_parameters'))
for ip=1:length(derived_parameters)
xbin_param_d{1,k}(ip,:)=all_xbin_max_d(ip,:);
xbin_param_d{2,k}(ip,:)=all_xbin_min_d(ip,:);
ks_param_d{1,k}(ip,:)=all_ks_max_d(ip,:);
ks_param_d{2,k}(ip,:)=all_ks_min_d(ip,:);
end
end
end
......@@ -196,8 +208,10 @@ try
save(fullfile(folder,variable_name,'pdf_eval_func.mat'),'xbinT','ks_y_T');
disp('saving probability density functions of all parameters')
save(fullfile(folder,variable_name,'pdf_param.mat'),'xbin_param','ks_param');
save(fullfile(folder,variable_name,'pdf_param_derived.mat'),'xbin_param_d','ks_param_d');
if(isfield(model,'derived_parameters'))
save(fullfile(folder,variable_name,'pdf_param_derived.mat'),'xbin_param_d','ks_param_d');
end
catch ME
......
CRA_Updated_Version/CRA_Toolbox_bash/compute_R0.m
function R0=compute_R0(p)%,Tlock2,s02,s12)
function R0=compute_R0(param)%,Tlock2,s02,s12)
be = p(1);
b0=p(2);
b1=p(3);
b2=p(4);
b3=p(5);
be = param(1);
b0=param(2);
b1=param(3);
b2=param(4);
b3=param(5);
N=10^5;
be=be/(N);
b0=b0/(N);
b1=b1/(N);
b2=b2/(N);
b3=b3/(N);
Norm=10^5; %882000
%Norm=882000;
be=be/(Norm);
b0=b0/(Norm);
b1=b1/(Norm);
b2=b2/(Norm);
b3=b3/(Norm);
FracSevere = p(6);
FracCritical = p(7);
FracSevere = param(6);
FracCritical = param(7);
FracMild = 1-FracSevere-FracCritical;
FracAsym = p(8);
FracAsym = param(8);
IncubPeriod = p(9);
DurMildInf = p(10);
IncubPeriod = param(9);
DurMildInf = param(10);
DurAsym = p(11);
DurAsym = param(11);
DurHosp = p(12);
TimeICUDeath = p(13);
ProbDeath=p(14);
PresymPeriod=p(15)*IncubPeriod;
DurHosp = param(12);
TimeICUDeath = param(13);
ProbDeath=param(14);
PresymPeriod=param(15)*IncubPeriod;
%PresymPeriod=param(15);
CFR=ProbDeath*FracCritical/100;
......@@ -54,4 +55,4 @@ g3=(1/TimeICUDeath)-u;
%s0 = (t<Tlock1) + s01 * (t>=Tlock1); %* (t<Tlock2) + s01(2)*(t>=Tlock2);
%s1 = (t<Tlock1) + s11 * (t>=Tlock1); %* (t<Tlock2) + s11(2)*(t>=Tlock2);
R0=N*(((be)/a1)+f*((b0)/g0)+(1-f)*(((b1)/(p1+g1))+(p1/(p1+g1))*(b2/(p2+g2)+ (p2/(p2+g2))*(b3/(u+g3)))));
\ No newline at end of file
R0=Norm*(((be)/a1)+f*((b0)/g0)+(1-f)*(((b1)/(p1+g1))+(p1/(p1+g1))*(b2/(p2+g2)+ (p2/(p2+g2))*(b3/(u+g3)))));
\ No newline at end of file
CRA_Updated_Version/CRA_Toolbox_bash/main_PulseGeneratorNetwork.m
......@@ -66,11 +66,11 @@ u=470;
num_observables=2;
observables_name={'R1','Y'};
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,'input',u,'total_proteins',0);
model=struct('name',model_name,'odesolver',ode_solver,'time',time_axis,'stop',stop_time,'step',step_size,'nominal_parameters',nominal_parameters,'nominal_parameters_name',{parameters_name},'num_observables',num_observables,'observables_name',{observables_name},'initial_conditions',x0,'input',u,'total_proteins',0);
Nr=10;
Nr=2;
folder='PulseGeneratorNetwork_results';
folder='2014sortedprovaparfor_PulseGeneratorNetwork_results';
LBpi=0.1;
UBpi=10;
......@@ -80,7 +80,8 @@ variable_name='Y';
current_func=Area(); %current evaluation function
tail_size=1000; %number of samples for the lower and upper tail
step_size=0.1; %this parameter needs to be defined only when current_tm=tmp_sum(step_size)
current_tm=tmp_sum(0.1);
%current_tm=tmp_sum(0.1);
current_tm=sorted();
%model simulation for each sample of the Latin Hypercube
disp('Starting model simulation with perturbed parameters');
......@@ -91,9 +92,11 @@ disp('All done! Model simulation completed!');
%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)
AllDerivedParam=cell(Nr,1);
compute_MIRI(model,variable_name,current_func,tail_size,current_tm,Nr,Ns,AllResults,AllPerturbations,AllDerivedParam,folder)
catch ME
break
%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');
......
CRA_Updated_Version/CRA_Toolbox_bash/main_covid19.m
......@@ -65,7 +65,7 @@ time_axis=[0:step_size:stop_time]';
% S0 = N-E00;
% x0 = [S0 E00 0 0 0 0 0 0 0];
%
N = 882000; %population Lombardia
N = 882000; %population Umbria
E00 = (1/N)*10^5;
S0 = 10^5-E00;
......@@ -80,9 +80,9 @@ x0 = [S0 E00 0 0 0 0 0 0 0];
%dopo 19 giorni (10 Marzo): total lockdown
Tlock=[3 14 17 19];
s0 = [0.7 0.3 0.2]; %reduction of 30%, 70%, 80%
s1=[0.6 0.25 0.15]; %reduction of 40%, 75%, 85%
s2=0.4; %reduction of 60%
s00 = [0.7 0.3 0.2]; %reduction of 30%, 70%, 80%
s11=[0.6 0.25 0.2]; %reduction of 40%, 75%, 80%
s22=0.4; %reduction of 60%
load('Examples/moda_parametri_inizio9step_LINLOG_100000_UMBRIA_I2I3DE_3may.mat');
......@@ -91,7 +91,7 @@ nominal_parameters=moda_parametri_inizio9step_LINLOG;
%s11=moda_parametri_inizio14step_LOG(17);
%nominal_parameters=[be b0 b1 b2 b3 a0 a1 f g0 g1 p1 g2 p2 g3 u];
nominal_parameters_name={'be', 'b0', 'b1', 'b2','b3', 'FracSevere', 'FracCritical', 'FracAsym', 'IncubPeriod', 'DurMildInf', 'DurAsym', 'DurHosp', 'TimeICUDeath', 'ProbDeath', 'PresymPeriod'};
nominal_parameters_name={'be', 'b0', 'b1', 'b2','b3', 'FracSevere', 'FracCritical', 'FracAsym', 'IncubPeriod', 'DurMildInf', 'DurAsym', 'DurHosp', 'TimeICUDeath', 'ProbDeath', 'PresymPercentage'};
derived_parameters=1;
derived_parameters_name={'R0'};
......@@ -99,27 +99,42 @@ derived_parameters_name={'R0'};
%Tlock1=23;
num_observables=9;
observables_name={'S','E','P_S','A','M','H','ICU','R','D'};
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,'s00',s0,'s11',s1,'s22',s2);
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,'s00',s00,'s11',s11,'s22',s22);
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.5 0.25 0.25 0.25 0.25 0.125 0.125 0.25 0.5 0.4311 0.6988 0.5556 0.4 0.5 0.7321];
UBpi=[2 4 4 4 4 3 3 2 2 2.8741 3.5532 2.7778 4 2 1.3903];
Ns=10000;
%LBpi=[0.5 0.25 0.25 0.25 0.25 0.125 0.125 0.25 0.5 0.4311 0.6988 0.5556 0.4 0.5 0.7321];
%UBpi=[2 4 4 4 4 3 3 2 2 2.8741 3.5532 2.7778 4 2 1.3];
variable_name='D';
%LBpi=[0.5 0.25 0.25 0.25 0.25 0.125 0.125 0.25 0.5 0.4311 0.6988 0.5556 0.4 0.5 0.7321];
%UBpi=[2 4 4 4 4 3 3 2 2 2.8741 3.5532 2.7778 4 2 1.3];
%load('intervalli_param_perCRA_stessi_delle_bande_stretti.mat')
%LBpi=interv_low;
%UBpi=interv_high;
%LBpi=[0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5];
%UBpi=[2 2 2 2 2 2 2 2 2 2 2 2 2 2 1.3235];
LBpi=0.5;
UBpi=1.3;
Ns=1000;
variable_name='I2';
current_func=Area(); %current evaluation function
tail_size=1000; %number of samples for the lower and upper tail
tail_size=100; %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_results_3may';
folder='parfor_p_2_altraprova';
%model simulation for each sample of the Latin Hypercube
disp('Starting model simulation with perturbed parameters');
......
CRA_Updated_Version/CRA_Toolbox_bash/main_covid19.m~ 0 → 100644
%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';
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];
%
N = 882000; %population Umbria
E00 = (1/N)*10^5;
S0 = 10^5-E00;
x0 = [S0 E00 0 0 0 0 0 0 0];
%Tlock=23;
%inizio epidemia: 21 Febbraio
%dopo 3 giorni (24 Febbraio) DPI
%dopo 14 giorni (5 Marzo): chiusura scuole
%dopo 17 giorni (8 Marzo): social distancing + no events
%dopo 19 giorni (10 Marzo): total lockdown
Tlock=[3 14 17 19];
s00 = [0.7 0.3 0.2]; %reduction of 30%, 70%, 80%
s1=[0.6 0.25 0.15]; %reduction of 40%, 75%, 85%
s2=0.4; %reduction of 60%
load('Examples/moda_parametri_inizio9step_LINLOG_100000_UMBRIA_I2I3DE_3may.mat');
nominal_parameters=moda_parametri_inizio9step_LINLOG;
%s01=moda_parametri_inizio14step_LOG(16);
%s11=moda_parametri_inizio14step_LOG(17);
%nominal_parameters=[be b0 b1 b2 b3 a0 a1 f g0 g1 p1 g2 p2 g3 u];
nominal_parameters_name={'be', 'b0', 'b1', 'b2','b3', 'FracSevere', 'FracCritical', 'FracAsym', 'IncubPeriod', 'DurMildInf', 'DurAsym', 'DurHosp', 'TimeICUDeath', 'ProbDeath', 'PresymPercentage'};
derived_parameters=1;
derived_parameters_name={'R0'};
%Tlock1=23;
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,'s00',s0,'s11',s1,'s22',s2);
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.5 0.25 0.25 0.25 0.25 0.125 0.125 0.25 0.5 0.4311 0.6988 0.5556 0.4 0.5 0.7321];
%UBpi=[2 4 4 4 4 3 3 2 2 2.8741 3.5532 2.7778 4 2 1.3];
%LBpi=[0.5 0.25 0.25 0.25 0.25 0.125 0.125 0.25 0.5 0.4311 0.6988 0.5556 0.4 0.5 0.7321];
%UBpi=[2 4 4 4 4 3 3 2 2 2.8741 3.5532 2.7778 4 2 1.3];
%load('intervalli_param_perCRA_stessi_delle_bande_stretti.mat')
%LBpi=interv_low;
%UBpi=interv_high;
%LBpi=[0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5];
%UBpi=[2 2 2 2 2 2 2 2 2 2 2 2 2 2 1.3235];
LBpi=0.5;
UBpi=1.3;
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='2014parfor_p_2_10';
%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
%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;
\ No newline at end of file
CRA_Updated_Version/CRA_Toolbox_bash/plotpdf_param.m
......@@ -20,6 +20,7 @@ derived_parameters_name=model.derived_parameters_name;
plot(xbin_param{1,k}(i,:),ks_param{1,k}(i,:),'b',xbin_param{2,k}(i,:),ks_param{2,k}(i,:),'r');
hold on
end
nominal_parameters_name{1,i}
xlabel(nominal_parameters_name{1,i},'Interpreter','none');
ylabel('pdf');
legend('Conditional upper density','Conditional lower density');
......
CRA_Updated_Version/CRA_Toolbox_bash/start_simulation.m
......@@ -34,8 +34,8 @@ for k=1:Nr
psim_all=cell(Ns,1);
parfor i=1:Ns
pi=perturbed_p(i,:);
%modelsim=@(t,x)feval(model.name,t,x,pi,model.input,model.total_proteins);
modelsim=@(t,x)feval(model.name,t,x,pi,model.Tlock1,model.s01,model.s11);
modelsim=@(t,x)feval(model.name,t,x,pi,model.input,model.total_proteins);
%modelsim=@(t,x)feval(model.name,t,x,pi,model.Tlock1,model.s01,model.s11);
[tsimi,psimi]=feval(model.odesolver,modelsim,model.time,model.initial_conditions);
% for t=1:model.num_observables
% psim_all{i,t}=psimi(:,t);
......
CRA_Updated_Version/CRA_Toolbox_bash/start_simulation_v2_apr2020.m
......@@ -28,35 +28,47 @@ for k=1:Nr
waitbar(k/Nr,h,['Realization number ',num2str(k)]);
perturbation=zeros(Ns,num_parameters);
for i=1:length(LBpi)
hypercube=LatinHypercube(LBpi(i),UBpi(i),Ns);
perturbation_i=hypercube.generate_sample(1);
perturbation(:,i)=perturbation_i;
end
hypercube=LatinHypercube(LBpi,UBpi,Ns);
perturbation=hypercube.generate_sample(num_parameters);
% perturbation=zeros(Ns,num_parameters);
% for b=1:num_parameters
% %LBpi(i)
% %UBpi(i)
% b
% hypercube=LatinHypercube(LBpi(b),UBpi(b),Ns);
% perturbation_i=hypercube.generate_sample(1);
% perturbation(:,b)=perturbation_i;
% end
perturbed_p=repmat(nominal_parameters,length(perturbation),1).*perturbation;
size(perturbed_p)
AllPerturbations{k,1}=perturbation;
AllPerturbations{k,1}=perturbation;
disp(strcat('Generated a Latin Hypercube of size',{' '},sprintf('%.0f',size(perturbation,1)),'x',sprintf('%.0f',size(perturbation,2))));
%psim_all=cell(Ns,model.num_observables);
psim_all=cell(Ns,1);
R0_all=zeros(Ns,1);
parfor i=1:Ns
lastwarn('')
parfor i=1:Ns %parfor
pi=perturbed_p(i,:);
%modelsim=@(t,x)feval(model.name,t,x,pi,model.input,model.total_proteins);
modelsim=@(t,x)feval(model.name,t,x,pi,model.Tlock,model.s00,model.s11,model.s22);
%modelsim=@(t,x)feval(model.name,t,x,pi,model.Tlock1,model.s01,model.s11);
[tsimi,psimi]=feval(model.odesolver,modelsim,model.time,model.initial_conditions);
R0i=compute_R0(pi);
%plot(tsimi, psimi(:,6))
%hold on;
% for t=1:model.num_observables
% psim_all{i,t}=psimi(:,t);
%end
msgstr=lastwarn;
msgstr=lastwarn
if isequal(msgstr,'')
psim_all{i,1}=psimi;
R0_all(i,1)=R0i;
else
psim_all{i,1}=NaN(length(model.time),model.num_observables);
R0_all{i,1}=NaN;
R0_all(i,1)=NaN;
end
lastwarn('')
end
......
CRA_Updated_Version/CRA_Toolbox_bash/untitled.m 0 → 100644
%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=70;
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];
%N = 882000; %population Umbria
%E00 = (1/N)*10^5;
%S0 = 10^5-E00;
%x0 = [S0 E00 0 0 0 0 0 0 0];
Tlock=23;
%inizio epidemia: 21 Febbraio
%dopo 3 giorni (24 Febbraio) DPI
%dopo 14 giorni (5 Marzo): chiusura scuole
%dopo 17 giorni (8 Marzo): social distancing + no events
%dopo 19 giorni (10 Marzo): total lockdown
%Tlock=[3 14 17 19];
%s0 = [0.7 0.3 0.2]; %reduction of 30%, 70%, 80%
%s1=[0.6 0.25 0.15]; %reduction of 40%, 75%, 85%
%s2=0.4; %reduction of 60%
load('Examples/moda_parametri_inizio14step_LOG_100000_UMBRIA_I2I3DE.mat');
nominal_parameters=moda_parametri_inizio14step_LOG(1:15);
s01=moda_parametri_inizio14step_LOG(16);
s11=moda_parametri_inizio14step_LOG(17);
%nominal_parameters=[be b0 b1 b2 b3 a0 a1 f g0 g1 p1 g2 p2 g3 u];
nominal_parameters_name={'be', 'b0', 'b1', 'b2','b3', 'FracSevere', 'FracCritical', 'FracAsym', 'IncubPeriod', 'DurMildInf', 'DurAsym', 'DurHosp', 'TimeICUDeath', 'ProbDeath', 'PresymPeriod'};
derived_parameters=1;
derived_parameters_name={'R0'};
%Tlock1=23;
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,'Tlock1',Tlock,'s01',s01,'s11',s11);
Nr=5;
%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.5 0.9 0.9 0.9 0.9 0.9 0.9 0.9 1 0.9 0.9 0.9 0.9 0.9 0.9];
%UBpi=[2 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1 1.1];
%LBpi=[0.5 0.25 0.25 0.25 0.25 0.125 0.125 0.25 0.5 0.4311 0.6988 0.5556 0.4 0.5 0.7321];
%UBpi=[2 4 4 4 4 3 3 2 2 2.8741 3.5532 2.7778 4 2 1.3];
%LBpi=[0.5 0.25 0.25 0.25 0.25 0.125 0.125 0.25 0.5 0.4311 0.6988 0.5556 0.4 0.5 0.7321];
%UBpi=[2 4 4 4 4 3 3 2 2 2.8741 3.5532 2.7778 4 2 1.3];
%load('intervalli_param_perCRA_stessi_delle_bande_stretti.mat')
%LBpi=interv_low;
%UBpi=interv_high;
%LBpi=[0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5];
%UBpi=[2 2 2 2 2 2 2 2 2 2 2 2 2 2 1.3235];
Ns=1000;
variable_name='I2';
current_func=Area(); %current evaluation function
tail_size=100; %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='2014provaparfor_v2';
%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
%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;
\ No newline at end of file
Styling with Markdown is supported
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!