start_simulation_v2_apr2020.m
2.85 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
%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