Upload new file

CRA_Updated_Version/CRA_Toolbox_GUI/compute_MIRI.m

+%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
\ No newline at end of file