Skip to content

Commit 1c0f997

Browse files
FS-CodeBaseFS-CodeBase
authored
Added scripts and functions
Scripts and functions to solve model, perform Sobol analysis, and make plots
1 parent ade9256 commit 1c0f997

15 files changed

+1717
-0
lines changed

‎barplot_sobol_indices.m

+34
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
function [BPLOT,LEG] = barplot_sobol_indices(FirstOrderIdx,TotalOrderIdx,...
2+
prms_str,leg_bool,axis_lab_bool)
3+
LEG = 0;
4+
C = linspecer(2);
5+
6+
SobolMat = [FirstOrderIdx.indices;TotalOrderIdx.indices];
7+
BPLOT = bar(SobolMat','grouped','FaceColor','flat'); hold on
8+
for c_idx = 1:size(SobolMat',2)
9+
BPLOT(c_idx).CData = repmat(C(c_idx,:),[numel(FirstOrderIdx.indices),1]);
10+
end
11+
ylim([0 1])
12+
yticks(0:.1:1)
13+
xticks(1:numel(prms_str))
14+
xticklabels(prms_str)
15+
16+
% Plot error bars
17+
ERRS = [FirstOrderIdx.error;TotalOrderIdx.error];
18+
errorbar([BPLOT(1).XEndPoints;BPLOT(2).XEndPoints],...
19+
SobolMat,ERRS,'k','linestyle','none','LineWidth',1);
20+
hold off
21+
if leg_bool
22+
LEG = legend({'FO','TE'},'Location',...
23+
'northwest',...
24+
'Interpreter','latex','FontSize',12,...
25+
'NumColumns',2);
26+
end
27+
if axis_lab_bool
28+
ylabel('Fraction of Variance','Interpreter','latex','FontSize',13)
29+
xlabel('Model Factor','Interpreter','latex','FontSize',13)
30+
end
31+
grid on
32+
set(gcf,'Position',[10 10 600 320])
33+
set(gca,'TickLabelInterpreter','latex','FontSize',13)
34+
ytickformat('%0.02f');

‎compute_sobol_indices_doubling_time.m

+63
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,63 @@
1+
% % SCRIPT: compute_sobol_indices_doubling_time
2+
% % AUTHOR: Fabian Santiago
3+
% % EMAIL: [email protected]
4+
% % DATE: 11/21/2020
5+
6+
% Make a folder to store DT Sobol indices if one does not exist
7+
if ~exist('sobol_indices_dt', 'dir')
8+
mkdir('sobol_indices_dt')
9+
end
10+
11+
% Load parameter information for slobal sensitivity analysis
12+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13+
prms_info = fun_model_parameter_ranges;
14+
15+
% Load model solutions for each contact scenario
16+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
17+
model_sols_file_id = dir('model_sols/');
18+
model_sols_file_id(1:2) = [];
19+
20+
% Set sub-sampling size for bootstrapping and number of resamples
21+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
22+
sub_sample_size = 1500*(sum(prms_info(:,1))+2);
23+
n_resamples = 2000; % Number of times to perform re-sampling
24+
y_sol_idx = 1; % index of model solution of interest
25+
26+
% Compute Sensitivity Indices
27+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
28+
for case_solution_idx = 1:numel(model_sols_file_id)
29+
tic
30+
% Load model solutions YA, YB, YA_Bj, and
31+
% parameter information prms_info and prms_str
32+
load(['model_sols/',model_sols_file_id(case_solution_idx).name])
33+
34+
% Display contact scenario being considered
35+
disp(['contact scenario: ',num2str(class_cap)]);
36+
37+
% Compute first and total-order Sobol indices
38+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
39+
[FirstOrderIdx,TotalOrderIdx,Stats] = ...
40+
fun_sobol_indices_by_measure(...
41+
y_sol_idx,...
42+
sub_sample_size,...
43+
n_resamples,...
44+
YA_cell,...
45+
YB_cell,...
46+
YA_Bj_cell,...
47+
prms_info(:,1)...
48+
);
49+
save(... % File name:
50+
['sobol_indices_dt/Sobol_dt_'...
51+
'sol',num2str(y_sol_idx),...
52+
'sub',num2str(sub_sample_size),...
53+
'rep',num2str(n_resamples),...
54+
'CC_',num2str(class_cap),...
55+
'Vu',num2str(Vu),'Vd',num2str(Vd),...
56+
'Vg',num2str(Vg),'Vf',num2str(Vf),'.mat'],...
57+
... % Variables to save:
58+
'FirstOrderIdx',...
59+
'TotalOrderIdx',...
60+
'Stats',...
61+
'class_cap');
62+
disp(['cs: ',num2str(class_cap),', time: ',num2str(toc/60),'min']);
63+
end

‎compute_sobol_indices_over_time.m

+54
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,54 @@
1+
% SCRIPT: compute_sobol_indices_doubling_time
2+
% AUTHOR: Fabian Santiago
3+
% EMAIL: [email protected]
4+
% DATE: 8/27/2021
5+
6+
% Make a folder to store Sobol indices if one does not exist
7+
if ~exist('sobol_indices_in_time', 'dir')
8+
mkdir('sobol_indices_in_time')
9+
end
10+
11+
% Load parameter information for global sensitivity analysis
12+
parameters_info = fun_model_parameter_ranges;
13+
14+
% Load model solutions for each contact scenario
15+
model_sols_file_id = dir('model_sols/');
16+
model_sols_file_id(1:2) = [];
17+
18+
% Set sub-sampling size for bootstrapping and number of resamples
19+
SubSampN = 1500*(sum(parameters_info(:,1))+2); % Sub-samples to compute Si from resampled solutions
20+
n_resamples = 2000; % Number of times to perform re-sampling
21+
y_sol_idx = 1;
22+
23+
% Compute Sensitivity Indices
24+
for case_solution_idx = 1:numel(model_sols_file_id)
25+
tic
26+
% Load model solutions YA, YB, YA_Bj, and
27+
% parameter information prms_info and prms_str
28+
load(['model_sols/',model_sols_file_id(case_solution_idx).name])
29+
30+
% Display contact scenario being considered
31+
disp(['contact scenario: ',num2str(class_cap)]);
32+
33+
% Compute first and total-order Sobol indices
34+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35+
[FirstOrderIdx,TotalOrderIdx,Stats] = ...
36+
fun_sobol_indices_by_solution_over_time(...
37+
y_sol_idx,...
38+
SubSampN,...
39+
n_resamples,...
40+
YA_cell,...
41+
YB_cell,...
42+
YA_Bj_cell,...
43+
parameters_info(:,1));
44+
45+
save(['sobol_indices_in_time/SobolT_',...
46+
'sol',num2str(y_sol_idx),...
47+
'sub',num2str(SubSampN),...
48+
'rep',num2str(n_resamples),...
49+
'CC_',num2str(class_cap),...
50+
'Vu',num2str(Vu),'Vd',num2str(Vd),...
51+
'Vg',num2str(Vg),'Vf',num2str(Vf),'.mat'],...
52+
'FirstOrderIdx','TotalOrderIdx','Stats','class_cap');
53+
disp(['cs: ',num2str(class_cap),', time: ',num2str(toc/60),'min']);
54+
end

‎fun_compute_model_solutions.m

+159
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,159 @@
1+
function fun_compute_model_solutions(CASES,reduced_infect,SubSampN)
2+
% FUNCTION: FUN_COMPUTE_MODEL_SOLUTIONS
3+
% AUTHOR: Fabian Santiago
4+
% EMAIL: [email protected]
5+
% DATE: 11/27/2021
6+
%
7+
% DESCRIPTION:
8+
% Computes model solutions for computation of the FIRST ORDER
9+
% and TOTAL ORDER sensitivity indices of parameters COV-19
10+
% Model indices following (Saltelli et al. 2008 notation).
11+
%
12+
% INPUTS:
13+
% CASES: A vector of class caps: 25, 50, 100 or 50000 (no cap)
14+
% reduced_infect: Undergraduate vaccination scenarios vector (U and D)
15+
16+
% Check for invalid class-caps
17+
for ccap_idx = unique(CASES(:,1)')
18+
if not(any(ccap_idx==[25,50,100,50000]))
19+
error(['Acceptable class caps: 25, 50, 100, or 50000',...
20+
' (no class cap)'])
21+
end
22+
end
23+
24+
% Check for invalid vacc entry (values outside of the interval [0,100])
25+
if any([[CASES(:,2)' CASES(:,3)']<0 [CASES(:,2)' CASES(:,3)']>100])
26+
error('Percentage of initally vaccinated must be a value in [0,100]')
27+
end
28+
29+
% Make a folder to store model solutions if one does not exist
30+
if ~exist('model_sols','dir')
31+
mkdir('model_sols')
32+
end
33+
34+
% Load parameters of interest and their ranges for Latin Hypercube Sampling
35+
parameters_info = fun_model_parameter_ranges;
36+
37+
% Load contact matrices and number of students by classification
38+
[ConMat,Con_livingMat,...
39+
nUVec,nDVec,nGVec,nFVec] = fun_initialize_contact_matrices(reduced_infect);
40+
41+
% Model solution length in days
42+
DaysRefined = 30; % Number of days model is solved daily
43+
MaxDays = 7*15; % End of semester (15 week semester)
44+
45+
% Refine solutions during the first 30 days, then every 3 days
46+
tspan = 24*[0:DaysRefined (DaysRefined+3):3:MaxDays];
47+
converttohrs = 24; % 24 hrs/day
48+
49+
% Number of times to sample the parameter space (ps)
50+
% SubSampN (user input)
51+
52+
% Generate Parameter Value Matrices: A, B
53+
A = fun_lhs_sampling_of_parameters(parameters_info,SubSampN);
54+
B = fun_lhs_sampling_of_parameters(parameters_info,SubSampN);
55+
56+
% Pre-allocate space for matrices B_Aj and A_Bj. Notations indicates
57+
% matrix A or B from above but the jth column is replaced with that of the
58+
% other matrix.
59+
B_Aj_cell = cell(size(parameters_info,1),1);
60+
A_Bj_cell = cell(size(parameters_info,1),1);
61+
62+
% Define a matrix C_i formed by all columns of B except the ith column,
63+
% which is taken from A
64+
for param = 1:size(parameters_info,1)
65+
% Only parameters whose sensitivity is being computed
66+
if logical(parameters_info(param,1))
67+
ABjtmp = A;
68+
BAjtmp = B;
69+
70+
% Matrix, where column J comes from matrix B and
71+
% all other k − 1 columns come from matrix A
72+
ABjtmp(:,param) = B(:,param);
73+
74+
% Matrix, where column J comes from matrix A and
75+
% all other k − 1 columns come from matrix B
76+
BAjtmp(:,param) = A(:,param);
77+
78+
% Save A_Bj and B_Aj
79+
A_Bj_cell{param} = ABjtmp;
80+
B_Aj_cell{param} = BAjtmp;
81+
end
82+
end
83+
84+
% Compute model solutions
85+
tic
86+
for ccap_idx = 1:size(CASES,1)
87+
TOC = toc;
88+
class_cap = CASES(ccap_idx,1);
89+
switch class_cap
90+
case 25
91+
con_idx = 1;
92+
case 50
93+
con_idx = 2;
94+
case 100
95+
con_idx = 3;
96+
case 50000
97+
con_idx = 4;
98+
end
99+
if class_cap ~= 50000
100+
disp(['Class Cap: ',num2str(class_cap)])
101+
else
102+
disp('Class Cap: None')
103+
end
104+
% Percent Faculty and Graduate Students Vaccinated
105+
Vu = CASES(ccap_idx,2); Vd = Vu;
106+
Vf = CASES(ccap_idx,3); Vg = Vf;
107+
108+
disp(['und vax: ',num2str(Vu),' & fac vax: ',num2str(Vf)])
109+
% Percent of undergraduates vaccinated
110+
111+
% Load a particular contact scenario
112+
[Con,Con_living,nu,nd,ng,nf] = deal(...
113+
ConMat(:,:,con_idx),...
114+
Con_livingMat,...
115+
nUVec(con_idx),...
116+
nDVec(con_idx),...
117+
nGVec(con_idx),...
118+
nFVec(con_idx)...
119+
);
120+
% Compute YA
121+
YA_cell = fun_solve_covid_model(A,tspan,converttohrs,Con,Con_living,...
122+
nu,nd,ng,nf,Vu,Vd,Vg,Vf);
123+
124+
% Compute YB
125+
YB_cell = fun_solve_covid_model(B,tspan,converttohrs,Con,Con_living,...
126+
nu,nd,ng,nf,Vu,Vd,Vg,Vf);
127+
128+
129+
% Compute YA_Bj
130+
YA_Bj_cell = cell(1,numel(parameters_info(:,1)));
131+
for param = find(parameters_info(:,1))'
132+
A_Bj_tmp = A_Bj_cell{param};
133+
YA_Bj_tmp = fun_solve_covid_model(...
134+
A_Bj_tmp,tspan,converttohrs,...
135+
Con,Con_living,nu,nd,ng,nf,...
136+
Vu,Vd,Vg,Vf...
137+
);
138+
YA_Bj_cell{param} = YA_Bj_tmp;
139+
end
140+
141+
% Save the solutions using case specific information
142+
save(... % File name:
143+
['./model_sols/model_sols_dailyMCn',...
144+
num2str(SubSampN),...
145+
'Vu',num2str(Vu),'Vd',num2str(Vd),...
146+
'Vg',num2str(Vg),'Vf',num2str(Vf),...
147+
'CC_',num2str(class_cap),'.mat'],...
148+
'YA_cell','YB_cell','YA_Bj_cell',...
149+
... % Variables:
150+
'parameters_info',...
151+
'Vu','Vd','Vg','Vf',...
152+
'class_cap',...
153+
'-v7.3');
154+
155+
disp(['und vax: ',num2str(Vu),...
156+
' & fac vax: ',num2str(Vf),...
157+
' & elapsed Time: ',num2str(round(toc/60,2)),' mins',...
158+
' & simulation Time: ',num2str(round((toc-TOC)/60,2)),' mins'])
159+
end

0 commit comments

Comments
 (0)