compute_MIRI.m
7.58 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
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
%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 variable. It plots boxplot of MIRI values and the pdf of the
%evaluation function of the chosen variable. 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
%create a directory to save all results
mkdir(handles1.folder,handles1.chosen_variable);
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_variable}{j,1},[1 length(handles1.time_results)]);
end
%creation of an object TimeBehavior where results of the chosen
%variable 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 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(handles1.chosen_variable);
set(fh,'CreateFcn','set(fh,''Visible'',''on'')')
saveas(fh,fullfile(handles1.folder,handles1.chosen_variable,strcat('hist_',handles1.chosen_variable,'_Nr',num2str(k))),'jpeg');
close(fh)
%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 (or bar if the realization is only one) are displayed in gui_MIRI
if handles1.Nr==1
bar(handles.axes1,MIRIT);
%set(handles.axes1,'XTickLabel',handles1.param_name);
% xticks=get(handles.axes1,'XTickLabel');
% xticklabel_rotate(xticks,90);
ylabel('MIRI');
else
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);
end
fh2 = figure('Visible','off');
h_new2=copyobj(handles.axes1, fh2);
set(h_new2, 'Units', 'Normalized');
set(h_new2,'OuterPosition',[.1, .1, .85, .85]);
set(gcf,'Visible','off','CreateFcn','set(gcf,''Visible'',''on'')')
saveas(fh2, fullfile(handles1.folder,handles1.chosen_variable,'MIRI'),'jpeg');
close(fh2);
%saving results
disp('saving results...');
save(fullfile(handles1.folder,handles1.chosen_variable,'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(handles1.folder,handles1.chosen_variable,'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(handles1.folder,handles1.chosen_variable,'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(handles1.folder,handles1.chosen_variable,'pdf_eval_func.mat'),'xbinT','ks_y_T');
disp('saving probability density functions of all parameters')
save(fullfile(handles1.folder,handles1.chosen_variable,'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 variable';
errordlg(msg);
close(h)
else
errordlg(ME.message);
close(h)
end
end
end