compute_MIRI.m
5.76 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
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
%Function executed when the button COMPUTE MIRI in the second GUI (gui_MIRI) is
%pushed. It takes in input the handles of both GUIs and it computes MIRI of
%the selected protein. It plots boxplot of MIRI values and the pdf of the
%evaluation function of the chosen protein. 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(handles,hObject,handles1)
MIRIT=zeros(handles1.Nr,handles1.num_param);
CloudCondMaxT=[];
CloudCondMinT=[];
xbinT=[];
ks_y_T=[];
try
xbin_param=cell(2,handles1.Nr);
ks_param=cell(2,handles1.Nr);
h = waitbar(0,'Please wait...');
for k=1:handles1.Nr
waitbar(k/handles1.Nr,h,['Realization number ',num2str(k)]);
Results_reshape=zeros(handles1.Nsample,length(handles1.time_results));
for j=1:handles1.Nsample
Results_reshape(j,:)=reshape(handles1.AllResults{k,handles1.num_protein}{j,1},[1 length(handles1.time_results)]);
end
%creation of an object TimeBehavior where results of the chosen
%protein are stored
handles.Results=TimeBehavior(handles1.time_results,Results_reshape);
%values of the evaluation function are computed and stored in the
%TimeBehavior object Results
handles.Results.currentEvalFunc=handles.current_func;
handles.Results.computeEvalFunc();
guidata(hObject,handles);
samples=handles.Results.evalFuncValues;
% an evaluation function may have identical values when parameters
% do not affect its value. For example, if a protein 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 protein.');
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];
%method for extracting the upper and lower tail of the evaluation
%function
handles.Results.currentTailMethod=handles.current_tm;
[XiMax,XiMin]=handles.Results.compute_tail(handles1.AllPerturbations{k,1},handles.tail_size);
%[XiMax,XiMin]=handles.Results.compute_tail(handles1.perturbation,handles.tail_size);
handles.XiMax=XiMax;
handles.XiMin=XiMin;
guidata(hObject,handles);
NcloudMax=size(handles.XiMax,1);
NcloudMin=size(handles.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(handles1.nominal_p, handles1.AllPerturbations{k,1}, handles.XiMax, handles.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);
MIRIT(k,:)=MIRI_param;
CloudCondMaxT=[CloudCondMaxT;CloudCondMax];
CloudCondMinT=[CloudCondMinT;CloudCondMin];
mean_CloudMax=mean(CloudCondMaxT);
mean_CloudMin=mean(CloudCondMinT);
for ip=1:length(handles1.nominal_p)
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
end
close (h)
%MIRI boxplot are displayed in gui_MIRI
boxplot(handles.axes1,MIRIT,'labels',handles1.param_name,'labelorientation','inline');
ylabel('MIRI');
labelSize = 13; %size of the label
set(findobj(handles.axes1,'Type','text'),'FontSize',labelSize);
%saving results
mkdir(handles1.folder,handles1.chosen_protein);
save(fullfile(handles1.folder,handles1.chosen_protein,'MIRIT.mat'),'MIRIT');
save(fullfile(handles1.folder,handles1.chosen_protein,'meanCloudCondMax.mat'),'mean_CloudMax');
save(fullfile(handles1.folder,handles1.chosen_protein,'meanCloudCondMin.mat'),'mean_CloudMin');
save(fullfile(handles1.folder,handles1.chosen_protein,'pdf_eval_func.mat'),'xbinT','ks_y_T');
save(fullfile(handles1.folder,handles1.chosen_protein,'pdf_param.mat'),'xbin_param','ks_param');
%the evaluation function pdf is displayed in gui_MIRI
plotpdf_evalfunc(handles,hObject,handles1);
catch ME
if (strcmp(ME.identifier,'MATLAB:badsubscript'))
msg='The pdf is a Dirac delta function: change evaluation function or protein';
errordlg(msg);
close(h)
else
errordlg(ME.message);
close(h)
end
end
end