%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,AllDerivedParam,folder)

nominal_parameters=model.nominal_parameters;
%num_parameters=length(nominal_parameters);
nominal_parameters_name=model.nominal_parameters_name;

if(isfield(model,'derived_parameters'))
   derived_parameters=model.derived_parameters;
   derived_parameters_name=model.derived_parameters_name;
   num_parameters=length(nominal_parameters)+length(derived_parameters);
   parameters_name=[nominal_parameters_name derived_parameters_name];
else
   num_parameters=length(nominal_parameters);
   parameters_name=[nominal_parameters_name];
end

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);
  
   xbin_param_d=cell(2,Nr);
   ks_param_d=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);
        
        if(isfield(model,'derived_parameters'))
           obj_MIRI_derived=MIRI();
           [all_xbin_max_d,all_xbin_min_d,all_ks_max_d,all_ks_min_d]=obj_MIRI_derived.estimate_parampdfs(derived_parameters, AllDerivedParam{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_d,CloudCondMin_d, CloudCondMax_d]=obj_MIRI_derived.MIRI_computation(all_ks_max_d,all_ks_min_d,all_xbin_max_d,all_xbin_min_d);
        
        
           MIRI_param=[MIRI_param MIRI_param_d];
           CloudCondMax=[CloudCondMax CloudCondMax_d];
           CloudCondMin=[CloudCondMin CloudCondMin_d];
        end
        
        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   
        
        if(isfield(model,'derived_parameters'))
            for ip=1:length(derived_parameters)
               xbin_param_d{1,k}(ip,:)=all_xbin_max_d(ip,:);
               xbin_param_d{2,k}(ip,:)=all_xbin_min_d(ip,:);
               ks_param_d{1,k}(ip,:)=all_ks_max_d(ip,:);
               ks_param_d{2,k}(ip,:)=all_ks_min_d(ip,:);
            end
        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:num_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');
   saveas(gcf, fullfile(folder,variable_name,'MIRI'),'fig');
      
   
   %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');
   
   if(isfield(model,'derived_parameters'))
      save(fullfile(folder,variable_name,'pdf_param_derived.mat'),'xbin_param_d','ks_param_d');
   end
      
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