Skip to content

SysBiOThe / CRA-Matlab

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

update robustness covid

1 parent 06a93842
Inline Side-by-side
Showing 32 changed files with 206 additions and 92 deletions
CRA_Updated_Version/CRA_Toolbox_bash/Examples/ode_covid19_v2.m
......@@ -3,20 +3,70 @@
function dx=ode_covid19_v2(t,x,p,Tlock1,s01,s11)%,Tlock2,s02,s12)
be = p(1);
b0 = p(2);
b1 = p(3);
b2 = p(4);
b3 = p(5);
a0 = p(6);
a1 = p(7);
f = p(8);
g0 = p(9);
g1 = p(10);
p1 = p(11);
g2 = p(12);
p2 = p(13);
g3 = p(14);
u = p(15);
b0=p(2);
b1=p(3);
b2=p(4);
b3=p(5);
N=882000;
be=be/(N);
b0=b0/(N);
b1=b1/(N);
b2=b2/(N);
b3=b3/(N);
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);
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;
% be = p(1);
% b0 = p(2);
% b1 = p(3);
% b2 = p(4);
% b3 = p(5);
% a0 = p(6);
% a1 = p(7);
% f = p(8);
% g0 = p(9);
% g1 = p(10);
% p1 = p(11);
% g2 = p(12);
% p2 = p(13);
% g3 = p(14);
% u = p(15);
%seas_amp = p(16);
%seas_phase = p(17);
%s0 = p(18);
......@@ -42,6 +92,7 @@ u = p(15);
s0 = (t<Tlock1) + s01 * (t>=Tlock1); %* (t<Tlock2) + s01(2)*(t>=Tlock2);
s1 = (t<Tlock1) + s11 * (t>=Tlock1); %* (t<Tlock2) + s11(2)*(t>=Tlock2);
%susceptibles individuals (not infected)
S = x(1);
%exposed individuals, infected but no symptoms and no transmission
......@@ -63,7 +114,7 @@ D = x(9);
dS = - (be *s0* E1 + b0*s0*I0 + b1*s1*I1 + b2*I2 + b3*I3)*S; % + k_e0s * E0;
dE0 = (be *s0* E1 + b0*s0*I0 + b1*s1*I1 + b2*I2 + b3*I3)*S -a0*E0; % - k_e0s * E0;
dE0 = (be *s0* E1 + b0*s0*I0 + b1*s1*I1 + b2*I2 + b3*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;
......
CRA_Updated_Version/CRA_Toolbox_bash/compute_MIRI.m
......@@ -29,11 +29,16 @@
function compute_MIRI(model,variable_name,current_func,tail_size,current_tm,Nr,Ns,AllResults,AllPerturbations,folder)
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);
parameters_name=model.parameters_name;
%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];
obs_name=model.observables_name;
......@@ -57,6 +62,9 @@ try
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
......@@ -123,7 +131,16 @@ 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);
MIRI_param=[MIRI_param MIRI_param_d];
MIRIT(k,:)=MIRI_param;
CloudCondMax=[CloudCondMax CloudCondMax_d];
CloudCondMin=[CloudCondMin CloudCondMin_d];
CloudCondMaxT=[CloudCondMaxT;CloudCondMax];
CloudCondMinT=[CloudCondMinT;CloudCondMin];
mean_CloudMax=mean(CloudCondMaxT);
......@@ -135,6 +152,14 @@ try
ks_param{1,k}(ip,:)=all_ks_max(ip,:);
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,:);
end
end
close (h)
......@@ -144,13 +169,14 @@ try
if Nr==1
bar(MIRIT);
%set(gca,'xticklabel',parameters_name);
xticklabel_rotate([1:length(nominal_parameters)],90,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
......@@ -170,6 +196,7 @@ 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');
catch ME
......
CRA_Updated_Version/CRA_Toolbox_bash/compute_R0.m 0 → 100644
function R0=compute_R0(p,Tlock1,s01,s11)%,Tlock2,s02,s12)
be = p(1);
b0=p(2);
b1=p(3);
b2=p(4);
b3=p(5);
N=882000;
be=be/(N);
b0=b0/(N);
b1=b1/(N);
b2=b2/(N);
b3=b3/(N);
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);
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=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
CRA_Updated_Version/CRA_Toolbox_bash/main_covid19.m
......@@ -66,93 +66,48 @@ 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'};
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=6;
observables_name={'E1','I0','I1','I2','I3','D'};
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,'parameters_name',{parameters_name},'num_observables',num_observables,'observables_name',{observables_name},'initial_conditions',x0,'Tlock1',Tlock1,'s01',s01,'s11',s11);
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',Tlock1,'s01',s01,'s11',s11);
Nr=3;
Nr=10;
%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;
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];
Ns=10000;
variable_name='I0';
current_func=Maximum(); %current evaluation function
tail_size=10; %number of samples for the lower and upper tail
variable_name='I3';
current_func=TimeOfMaximum(); %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_I0_results';
folder='covid19_TimeMax_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);
[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,folder)
compute_MIRI(model,variable_name,current_func,tail_size,current_tm,Nr,Ns,AllResults,AllPerturbations,AllDerivedParam,folder)
catch ME
%break
return
......
CRA_Updated_Version/CRA_Toolbox_bash/plotpdf_evalfunc.m
......@@ -13,7 +13,7 @@ function plotpdf_evalfunc(folder,variable_name)
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
......
CRA_Updated_Version/CRA_Toolbox_bash/plotpdf_param.m
......@@ -7,8 +7,10 @@
function plotpdf_param(folder,variable_name,model,Nr)
nominal_parameters=model.nominal_parameters;
parameters_name=model.parameters_name;
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'));
......@@ -18,10 +20,24 @@ parameters_name=model.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
xlabel(parameters_name{1,i},'Interpreter','none');
xlabel(nominal_parameters_name{1,i},'Interpreter','none');
ylabel('pdf');
legend('Conditional upper density','Conditional lower density');
savefig(fullfile(folder,variable_name,parameters_name{1,i}));
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)
......
CRA_Updated_Version/CRA_Toolbox_bash/start_simulation_v2_apr2020.m
......@@ -11,7 +11,7 @@
%It returns in output the array AllResults containing all the output
%variables simulated using the parameter vectors in AllPerturbations
function [AllResults, AllPerturbations]=start_simulation_v2_apr2020(model,Nr,LBpi,UBpi,Ns)
function [AllResults, AllPerturbations,AllDerivedParam]=start_simulation_v2_apr2020(model,Nr,LBpi,UBpi,Ns)
nominal_parameters=model.nominal_parameters;
num_parameters=length(nominal_parameters);
......@@ -20,6 +20,8 @@ AllResults=cell(Nr,model.num_observables);
AllPerturbations=cell(Nr,1);
AllDerivedParam=cell(Nr,1);
h = waitbar(0,'Please wait...');
for k=1:Nr
......@@ -38,22 +40,28 @@ for k=1:Nr
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
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);
[tsimi,psimi]=feval(model.odesolver,modelsim,model.time,model.initial_conditions);
R0i=compute_R0(pi,model.Tlock1,model.s01,model.s11);
% for t=1:model.num_observables
% psim_all{i,t}=psimi(:,t);
%end
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;
end
lastwarn('')
end
end
AllDerivedParam{k,1}=R0_all;
for j=1:model.num_observables
for w=1:Ns
......
Styling with Markdown is supported
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!