Skip to content

SysBiOThe / CRA-Matlab

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

covid19

1 parent f1e72144
Inline Side-by-side
Showing 14 changed files with 742 additions and 0 deletions
CRA_Updated_Version/CRA_Toolbox_bash/COVID19/compute_MIRI.m 0 → 100644
%This function computes MIRI of the parameter vector for the selected output variable.
%It plots boxplot of MIRI values and the pdf of the
%evaluation function of the chosen variable.
%It takes in input the following parameters:
%1. model is the struct where all the features of the model are saved;
%2. variable_name is the name of the output node of interest;
%3. current_func is the chosen evaluation function to compute;
%4. tail_size is the number of samples to include in the upper and lower
%tails of the conditional parameter pdfs;
%5. current_tm is the method chosen for computing the tails of the
%parameter pdfs;
%6. Nr is the number of independent realizations;
%7. Ns is the number of samples generated with the Latin Hypercube;
%8. AllResults and AllPerturbations are the cell arrays in output from the
%function start_simulation;
%9. folder is the name of the folder where the results are saved;
%It saves the following variables:
% 1) MIRIT is the array containing MIRI values of all realizations
% 2) mean_CloudMax is the array containing the mean of values in
% CloudCondMaxT. CloudCondMaxT is the array of the mode of parameter values
% that give rise to the upper tail of the evaluation function
% 3) mean_CloudMin is the array containing the mean of values in
% CloudCondMinT. CloudCondMinT is the array of the mode of parameter values
% that give rise to the lower tail of the evaluation function
% 4) pdf_eval_func contains values of the evaluation function in all
% realizations
% 5) pdf_param contains values of the conditional pdf of parameters in all
% realizations
function compute_MIRI(model,variable_name,current_func,tail_size,current_tm,Nr,Ns,AllResults,AllPerturbations,AllDerivedParam,folder)
nominal_parameters=model.nominal_parameters;
%num_parameters=length(nominal_parameters);
nominal_parameters_name=model.nominal_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;
%extract index of variable_name in the ODE model
index=strfind(obs_name,variable_name);
variable_number=find(not(cellfun('isempty',index)));
MIRIT=zeros(Nr,num_parameters);
CloudCondMaxT=[];
CloudCondMinT=[];
xbinT=[];
ks_y_T=[];
try
%create a directory to save all results
mkdir(folder,variable_name);
xbin_param=cell(2,Nr);
ks_param=cell(2,Nr);
xbin_param_d=cell(2,Nr);
ks_param_d=cell(2,Nr);
h = waitbar(0,'Please wait...');
for k=1:Nr
waitbar(k/Nr,h,['Realization number ',num2str(k)]);
Results_reshape=zeros(Ns,length(model.time));
for j=1:Ns
Results_reshape(j,:)=reshape(AllResults{k,variable_number}{j,1},[1 length(model.time)]);
end
%creation of an object TimeBehavior where results of the chosen
%variable are stored
Results=TimeBehavior(model.time,Results_reshape);
%values of the evaluation function are computed and stored in the
%TimeBehavior object Results
Results.currentEvalFunc=current_func;
Results.computeEvalFunc();
samples=Results.evalFuncValues;
% an evaluation function may have identical values when parameters
% do not affect its value. For example, if a variable has a
% descending temporal behavior, the maximum and time of maximum
% will always be the same, independently of parameter perturbation
if (length(unique(samples))==1)
ME=MException('Matlab:EvalFuncValuesNull','Array of evaluation function contains identical values. Change evaluation function or variable.');
throw(ME)
close(h)
end
%creation of an object pdfEstimator for estimating the pdf of the
%evaluation function
pdf_obj=pdfEstimator();
BinEdges=[min(samples):(max(samples)-min(samples))/length(samples):max(samples)];
[ks_y,xbin]=pdf_obj.evaluate_pdf(samples,BinEdges);
xbinT=[xbinT;xbin];
ks_y_T=[ks_y_T;ks_y];
%in addition to the pdf it is possible to plot the histogram of the
%evaluation function
fh = figure('Visible','off');
hist(samples);
xlabel(variable_name);
set(fh,'CreateFcn','set(fh,''Visible'',''on'')')
saveas(fh,fullfile(folder,variable_name,strcat('hist_',variable_name,'_Nr',num2str(k))),'jpeg');
close(fh)
%method for extracting the upper and lower tail of the evaluation
%function
Results.currentTailMethod=current_tm;
[XiMax,XiMin]=Results.compute_tail(AllPerturbations{k,1},tail_size);
[XiMaxD,XiMinD]=Results.compute_tail(AllDerivedParam{k,1},tail_size);
%[XiMax,XiMin]=handles.Results.compute_tail(handles1.perturbation,handles.tail_size);
NcloudMax=size(XiMax,1);
NcloudMin=size(XiMin,1);
Ncloud=min([NcloudMax NcloudMin]);
%creation of an object MIRI for computing MIRI and conditional pdfs
%of parameters
obj_MIRI=MIRI();
[all_xbin_max,all_xbin_min,all_ks_max,all_ks_min]=obj_MIRI.estimate_parampdfs(nominal_parameters, AllPerturbations{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,CloudCondMin, CloudCondMax]=obj_MIRI.MIRI_computation(all_ks_max,all_ks_min,all_xbin_max,all_xbin_min);
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}, XiMaxD, XiMinD);
%[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];
CloudCondMax=[CloudCondMax CloudCondMax_d];
CloudCondMin=[CloudCondMin CloudCondMin_d];
end
MIRIT(k,:)=MIRI_param;
CloudCondMaxT=[CloudCondMaxT;CloudCondMax];
CloudCondMinT=[CloudCondMinT;CloudCondMin];
mean_CloudMax=mean(CloudCondMaxT);
mean_CloudMin=mean(CloudCondMinT);
for ip=1:length(nominal_parameters)
xbin_param{1,k}(ip,:)=all_xbin_max(ip,:);
xbin_param{2,k}(ip,:)=all_xbin_min(ip,:);
ks_param{1,k}(ip,:)=all_ks_max(ip,:);
ks_param{2,k}(ip,:)=all_ks_min(ip,:);
end
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
close (h)
%MIRI boxplot (or bar if the realization is only one) are displayed in gui_MIRI
figure
if Nr==1
bar(MIRIT);
%set(gca,'xticklabel',parameters_name);
xticklabel_rotate([1:num_parameters],90,parameters_name);
ylabel('MIRI');
else
boxplot(MIRIT,'labels',parameters_name,'labelorientation','inline');
ylabel('MIRI');
end
saveas(gcf, fullfile(folder,variable_name,'MIRI'),'jpeg');
saveas(gcf, fullfile(folder,variable_name,'MIRI'),'fig');
%saving results
disp('saving results...');
save(fullfile(folder,variable_name,'MIRIT.mat'),'MIRIT');
disp(strcat('Array of MIRI has size',{' '},sprintf('%.0f',size(MIRIT,1)),'x',sprintf('%.0f',size(MIRIT,2))));
disp('saving mode of the conditional upper pdf of the parameter vector');
save(fullfile(folder,variable_name,'meanCloudCondMax.mat'),'mean_CloudMax');
disp(strcat('The mode vector of the upper pdf has size',{' '},sprintf('%.0f',size(mean_CloudMax,1)),'x',sprintf('%.0f',size(mean_CloudMax,2))));
disp('saving mode of the conditional lower pdf of the parameter vector');
save(fullfile(folder,variable_name,'meanCloudCondMin.mat'),'mean_CloudMin');
disp(strcat('The mode vector of the upper pdf has size',{' '},sprintf('%.0f',size(mean_CloudMin,1)),'x',sprintf('%.0f',size(mean_CloudMin,2))));
disp('saving the probability density function of the evaluation function')
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');
if(isfield(model,'derived_parameters'))
save(fullfile(folder,variable_name,'pdf_param_derived.mat'),'xbin_param_d','ks_param_d');
end
catch ME
if (strcmp(ME.identifier,'MATLAB:badsubscript'))
msg='The pdf is a Dirac delta function: change evaluation function or variable';
errordlg(msg);
close(h)
else
errordlg(ME.message);
close(h)
end
end
CRA_Updated_Version/CRA_Toolbox_bash/COVID19/compute_R0.m 0 → 100644
function R0=compute_R0(param)%,Tlock2,s02,s12)
be = param(1);
b0=param(2);
b1=param(3);
b2=param(4);
b3=param(5);
Norm=10^5; %882000
%Norm=882000;
be=be/(Norm);
b0=b0/(Norm);
b1=b1/(Norm);
b2=b2/(Norm);
b3=b3/(Norm);
FracSevere = param(6);
FracCritical = param(7);
FracMild = 1-FracSevere-FracCritical;
FracAsym = param(8);
IncubPeriod = param(9);
DurMildInf = param(10);
DurAsym = param(11);
DurHosp = param(12);
TimeICUDeath = param(13);
ProbDeath=param(14);
PresymPeriod=param(15)*IncubPeriod;
%PresymPeriod=param(15);
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;
%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=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 \ No newline at end of file
CRA_Updated_Version/CRA_Toolbox_bash/COVID19/compute_intervention.m 0 → 100644
function [s0,s1,s2]=compute_intervention(t,Tlock,s00,s11,s22,region_name)
if(strcmp(region_name,'Italy'))
s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(2)) + s00(2)*(t>=Tlock(2)) * (t<Tlock(3)) + s00(3)*(t>=Tlock(3))*(t<Tlock(4)) + s00(4)*(t>=Tlock(4));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
if(strcmp(region_name,'Italy_pred'))
s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(2)) + s00(2)*(t>=Tlock(2)) * (t<Tlock(3)) + s00(3)*(t>=Tlock(3))*(t<Tlock(4)) + s00(4)*(t>=Tlock(4))*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6))*(t<Tlock(7))+ s00(7)*(t>=Tlock(7))*(t<Tlock(8))+ s00(8)*(t>=Tlock(8))*(t<Tlock(9))+ s00(9)*(t>=Tlock(9));
s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
if(strcmp(region_name,'Umbria'))
s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(2)) + s00(2)*(t>=Tlock(2));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
if(strcmp(region_name,'Umbria_pred'))
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))*(t<Tlock(5)) + s00(4)*(t>=Tlock(5))*(t<Tlock(6)) + s00(5)*(t>=Tlock(6))*(t<Tlock(7)) + s00(6)*(t>=Tlock(7))*(t<Tlock(8)) + s00(7)*(t>=Tlock(8));
s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
end
end
CRA_Updated_Version/CRA_Toolbox_bash/COVID19/italy/init_italia.m 0 → 100644
region_name='Italy_pred';
%define initial conditions
N = 60e6; %population Italy
%define time
%time = 0:1:110;
E00 = (1/N)*10^5;
S0 = 10^5-E00;
x0 = [S0 E00 0 0 0 0 0 0 0];
%epidemic start: 27 January
%after 24 days (20 February): first covid-19 patient
%after 25 days (21 February): red areas in Lombardia e Veneto
%after 28 days (24 February): school closure in Emilia-Romagna, Friuli-Venezia Giulia,
%Liguria, Lombardia, Marche, Piemonte, Veneto
%after 25 days (21 February) PPE
%after 38 days (5 March): school closure in the rest of Italy
%after 41 days (8 March): lockdown northern of italy + social distancing
%after 43 days (10 March): total lockdown
%after 98 days (4 May 2020): first reopenings
%after 112 days (18 May 2020): second lockdown release
%after 140 days (15 June 2020): apertura discoteche
%after 203 days (17 August 2020): chiusura discoteche
%after 231 days (14 September 2020): apertura scuole
Tlock=[25 28 38 43 98 112 140 203 231];
CRA_Updated_Version/CRA_Toolbox_bash/COVID19/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_OK';
stop_time=250;
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('italy/italia_moda_nr4.mat');
moda_parametri_inizio9step_LINLOG=italia_moda_nr4;
nominal_parameters=[moda_parametri_inizio9step_LINLOG(1,1:19) 1 1 1 1 1 moda_parametri_inizio9step_LINLOG(1,20:21)];
%nominal_parameters=moda_parametri_inizio9step_LINLOG(3,:);
nominal_parameters_name={'be', 'b0', 'b1', 'b2','b3', 'FracSevere', 'FracCritical', 'FracAsym', 'IncubPeriod', 'DurMildInf', 'DurAsym', 'DurHosp', 'TimeICUDeath', 'ProbDeath', 'PresymPeriod','s01','s02','s03','s04','s05','s06','s07','s08','s09','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;
load('italy/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 0.1 0.1 0.1 parameters{1,1}(1,20:21)]./nominal_parameters;
UBpi=[parameters{1,1}(2,1:19) 0.9 0.9 0.9 0.9 0.9 parameters{1,1}(2,20:21)]./nominal_parameters;
Ns=10000;
variable_name='I3';
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='ITALIA_3Nr_Area_time250_discoteche_scuole'; %where results are saved
%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;
\ No newline at end of file \ No newline at end of file
CRA_Updated_Version/CRA_Toolbox_bash/COVID19/ode_covid19_v3_OK.m 0 → 100644
%ODE model to describe the spread and clinical progression of COVID-19
function dx=ode_covid19_v3_OK(t,x,p,Tlock,region_name) %s00,s11,s22,region_name)%,Tlock2,s02,s12)
be = p(1);
b0=p(2);
b1=p(3);
b2=p(4);
b3=p(5);
Norm=10^5;
be=be/(Norm);
b0=b0/(Norm);
b1=b1/(Norm);
b2=b2/(Norm);
b3=b3/(Norm);
FracSevere = p(6);
FracCritical = p(7);
FracMild = 1-FracSevere-FracCritical;
FracAsym = p(8);
IncubPeriod = p(9);
DurMildInf = p(10);
DurAsym = p(11);
DurHosp = p(12);
TimeICUDeath = p(13);
ProbDeath=p(14);
PresymPeriod=p(15)*IncubPeriod;
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;
s00=[p(16) p(17) p(18) p(19) p(20) p(21) p(22) p(23) p(24)];
s11=p(25);
s22=p(26);
[s0,s1,s2]=compute_intervention(t,Tlock,s00,s11,s22,region_name);
%susceptibles individuals (not infected)
S = x(1);
%exposed individuals, infected but no symptoms and no transmission
E0 = x(2);
%exposed individuals, infected but no symptoms and can transmit
E1 = x(3);
%infected individuals but asymptomatic infection
I0 = x(4);
%infected individuals, mild infection
I1 = x(5);
%infected individuals, severe infection
I2 = x(6);
%infected individuals, critical infection
I3 = x(7);
%recovered
R = x(8);
%dead
D = x(9);
%dS = - (be *(1/(1+s0^n))* E1 + b0*(1/(1+s0^n))*I0 + b1*(1/(1+s1^n))*I1 + b2*(1/(1+s2^n))*I2 + b3*(1/(1+s2^n))*I3)*S; % + k_e0s * E0;
%dE0 = (be *(1/(1+s0^n))* E1 + b0*(1/(1+s0^n))*I0 + b1*(1/(1+s1^n))*I1 + b2*(1/(1+s2^n))*I2 + b3*(1/(1+s2^n))*I3)*S -a0*E0; % - k_e0s * E0;
dS = - ((be * s0)* E1 + (b0*s0)*I0 + (b1*s1)*I1 + (b2*s2)*I2 + (b3*s2)*I3)*S; % + k_e0s * E0;
dE0 = ((be *s0)* E1 + (b0*s0)*I0 + (b1*s1)*I1 + (b2*s2)*I2 + (b3*s2)*I3)*S -a0*E0; % - k_e0s * E0;
dE1 = a0 * E0 - a1*E1; % -g4*E1;
dI0 = f * a1 * E1 - g0*I0;
dI1 = (1-f) * a1 * E1 - g1*I1 - p1*I1;
dI2 = p1 * I1 - g2*I2 - p2*I2; % -u2*I2;
dI3 = p2 * I2 - g3 *I3 - u * I3 ;
dR = g0*I0 + g1*I1 + g2*I2 + g3*I3; % +g4 * E1;
dD = u * I3; % + u2 * I2;
dx=[dS;dE0;dE1;dI0;dI1;dI2;dI3;dR;dD];
CRA_Updated_Version/CRA_Toolbox_bash/COVID19/plotpdf_evalfunc.m 0 → 100644
%Function that plots and saves the figure of the estimated pdf of the evaluation function
%for the chosen variable.
function plotpdf_evalfunc(folder,variable_name)
try
load(fullfile(folder,variable_name,'pdf_eval_func.mat'));
figure
for i=1:size(xbinT,1)
plot(xbinT(i,:),ks_y_T(i,:));
hold on
end
xlabel(variable_name,'Interpreter','none');
ylabel('pdf of evaluation function');
saveas(gcf, fullfile(folder,variable_name,variable_name),'jpeg');
saveas(gcf, fullfile(folder,variable_name,variable_name),'fig');
catch ME
errordlg(ME.message)
end
end
\ No newline at end of file \ No newline at end of file
CRA_Updated_Version/CRA_Toolbox_bash/COVID19/plotpdf_param.m 0 → 100644
%Function that plots and saves the figures of the conditional pdfs of parameters that are
%used for computing MIRI. It takes in input the struct of the model, the
%folder where the results are saved, the output node of interest and the
%number of independent realizations
function plotpdf_param(folder,variable_name,model,Nr)
nominal_parameters=model.nominal_parameters;
nominal_parameters_name=model.nominal_parameters_name;
derived_parameters=model.derived_parameters;
derived_parameters_name=model.derived_parameters_name;
try
load(fullfile(folder,variable_name,'pdf_param.mat'));
for i=1:length(nominal_parameters)
figure
for k=1:Nr
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');
savefig(fullfile(folder,variable_name,nominal_parameters_name{1,i}));
end
load(fullfile(folder,variable_name,'pdf_param_derived.mat'));
for i=1:length(derived_parameters)
figure
for k=1:Nr
plot(xbin_param_d{1,k}(i,:),ks_param_d{1,k}(i,:),'b',xbin_param_d{2,k}(i,:),ks_param_d{2,k}(i,:),'r');
hold on
end
xlabel(derived_parameters_name{1,i},'Interpreter','none');
ylabel('pdf');
legend('Conditional upper density','Conditional lower density');
savefig(fullfile(folder,variable_name,derived_parameters_name{1,i}));
end
catch ME
errordlg(ME.message)
end
end
\ No newline at end of file \ No newline at end of file
CRA_Updated_Version/CRA_Toolbox_bash/COVID19/start_simulation_v2_apr2020.m 0 → 100644
%Function for simulating the ODE model using the perturbed parameter
%vector. It takes in input:
%1. model which is the struct containing all the characheristics of the
%model
%2. Nr is the number of independent realizations to perform;
%3. LBpi and UBpi are the lower and upper boundaries of the Latin
%Hypercube. They are arrays since each parameter may have its own range of
%variation;
%4. Ns is the number of samples to generate for creating the Latin
%Hypercube;
%It returns in output the array AllResults containing all the output
%variables simulated using the parameter vectors in AllPerturbations
function [AllResults, AllPerturbations,AllDerivedParam]=start_simulation_v2_apr2020(model,Nr,LBpi,UBpi,Ns)
nominal_parameters=model.nominal_parameters;
num_parameters=length(nominal_parameters);
AllResults=cell(Nr,model.num_observables);
AllPerturbations=cell(Nr,1);
AllDerivedParam=cell(Nr,1);
h = waitbar(0,'Please wait...');
parpool(4)
for k=1:Nr
waitbar(k/Nr,h,['Realization number ',num2str(k)]);
%hypercube=LatinHypercube(LBpi,UBpi,Ns);
%perturbation=hypercube.generate_sample(num_parameters);
perturbation=zeros(Ns,num_parameters);
for b=1:num_parameters
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;
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);
lastwarn('')
parfor l=1:Ns %parfor
pi=perturbed_p(l,:);
%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.region_name);
%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;
if isequal(msgstr,'')
psim_all{l,1}=psimi;
R0_all(l,1)=R0i;
else
psim_all{l,1}=NaN(length(model.time),model.num_observables);
R0_all(l,1)=NaN;
end
lastwarn('')
end
AllDerivedParam{k,1}=R0_all;
for j=1:model.num_observables
for w=1:Ns
AllResults{k,j}{w,1}=psim_all{w,1}(:,j);
end
end
%AllResults{k,1:end}=psim_all;
disp(strcat('Completed one realization. The array of results has size',{' '},sprintf('%.0f',size(psim_all,1)),'x',sprintf('%.0f',size(psim_all,2))));
end
delete(gcp) %close parallel pool
close (h) %close the waitbar
end
\ No newline at end of file \ No newline at end of file
CRA_Updated_Version/CRA_Toolbox_bash/COVID19/umbria/init_umbria.m 0 → 100644
region_name='Umbria_pred';
%define time
%time = 0:1:80;
%define initial conditions
N = 882000; %population Umbria
E00 = (1/N)*10^5;
S0 = 10^5-E00;
x0 = [S0 E00 0 0 0 0 0 0 0];
%epidemic start: 21 February
%after 3 days (24 February): PPE
%after 13 days (5 March): school closure
%after 16 days (8 March): social distancing + no events
%after 18 days (10 March): total lockdown
%after 73 days (4 May): first reopenings
%after 87 days (18 May): second lockdown release
%after 118 days (19 June): reopening discos and clubs
%after 176 days (17 August): closing discos and clubs
%after 204 days (14 September): school reopening
Tlock=[3 13 18 73 87 118 176 204];
Styling with Markdown is supported
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!