Upload new file

CRA_Updated_Version/CRA_Toolbox_bash/compute_MIRI.m

+%This function computes MIRI of the parameter vector for the selected output variable. 
+%It plots boxplot of MIRI values and the pdf of the
+%evaluation function of the chosen variable.
+%It takes in input the following parameters:
+%1. model is the struct where all the features of the model are saved;
+%2. variable_name is the name of the output node of interest;
+%3. current_func is the chosen evaluation function to compute;
+%4. tail_size is the number of samples to include in the upper and lower
+%tails of the conditional parameter pdfs;
+%5. current_tm is the method chosen for computing the tails of the
+%parameter pdfs;
+%6. Nr is the number of independent realizations;
+%7. Ns is the number of samples generated with the Latin Hypercube;
+%8. AllResults and AllPerturbations are the cell arrays in output from the
+%function start_simulation;
+%9. folder is the name of the folder where the results are saved;
+%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(model,variable_name,current_func,tail_size,current_tm,Nr,Ns,AllResults,AllPerturbations,folder)
+
+nominal_parameters=model.nominal_parameters;
+num_parameters=length(nominal_parameters);
+parameters_name=model.parameters_name;
+obs_name=model.observables_name;
+
+
+%extract index of variable_name in the ODE model
+index=strfind(obs_name,variable_name);
+variable_number=find(not(cellfun('isempty',index)));
+
+MIRIT=zeros(Nr,num_parameters);
+CloudCondMaxT=[];
+CloudCondMinT=[];
+xbinT=[];
+ks_y_T=[];
+
+
+try
+   
+   %create a directory to save all results
+   mkdir(folder,variable_name);
+   
+    
+   xbin_param=cell(2,Nr);
+   ks_param=cell(2,Nr);
+  
+   h = waitbar(0,'Please wait...');
+   
+   for k=1:Nr 
+        waitbar(k/Nr,h,['Realization number ',num2str(k)]);
+        
+        Results_reshape=zeros(Ns,length(model.time));
+        for j=1:Ns
+            Results_reshape(j,:)=reshape(AllResults{k,variable_number}{j,1},[1 length(model.time)]);
+        end
+        
+        %creation of an object TimeBehavior where results of the chosen
+        %variable are stored
+        Results=TimeBehavior(model.time,Results_reshape);
+        
+        %values of the evaluation function are computed and stored in the
+        %TimeBehavior object Results
+        Results.currentEvalFunc=current_func;
+        Results.computeEvalFunc();  
+        
+        samples=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(variable_name);
+        set(fh,'CreateFcn','set(fh,''Visible'',''on'')')
+        saveas(fh,fullfile(folder,variable_name,strcat('hist_',variable_name,'_Nr',num2str(k))),'jpeg');
+        close(fh)
+        
+        %method for extracting the upper and lower tail of the evaluation
+        %function
+        Results.currentTailMethod=current_tm;
+        
+        [XiMax,XiMin]=Results.compute_tail(AllPerturbations{k,1},tail_size);
+        %[XiMax,XiMin]=handles.Results.compute_tail(handles1.perturbation,handles.tail_size);
+        
+        NcloudMax=size(XiMax,1);
+        NcloudMin=size(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(nominal_parameters, AllPerturbations{k,1}, XiMax, 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(nominal_parameters)
+            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
+   figure
+   if Nr==1
+       bar(MIRIT);
+       %set(gca,'xticklabel',parameters_name);
+       xticklabel_rotate([1:length(nominal_parameters)],90,parameters_name);
+       ylabel('MIRI');
+   else
+       boxplot(MIRIT,'labels',parameters_name,'labelorientation','inline');
+       ylabel('MIRI');
+   end
+   saveas(gcf, fullfile(folder,variable_name,'MIRI'),'jpeg');
+      
+   
+   %saving results
+   disp('saving results...');
+   save(fullfile(folder,variable_name,'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(folder,variable_name,'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(folder,variable_name,'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(folder,variable_name,'pdf_eval_func.mat'),'xbinT','ks_y_T');
+   disp('saving probability density functions of all parameters')
+   save(fullfile(folder,variable_name,'pdf_param.mat'),'xbin_param','ks_param');
+   
+      
+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