function  matrix_analysis_v3(res_data, minSSD, project_name, funct, params, X, fields, ...
                smooth_x_intervals, num_rnd_sets , ...
                SSD_matrix, SSD_levels, SSD_level_toshow, SSD_levels_colors, DPM_path)
% Analysis of parameter correlation matrices
% v.2 - vectorized version
%               761 sec - before optimization
%               < 30 sec after. Vectorized code is really fast!!!
%
% 3/21/08 - uses function with swapped order of parameters
%                - uses X with column-oriented structure
% 3/21/08 - adapted for LineShapeKin analysis: 
%               -  don't do data presentation, only pull out equivalent data sets
%
% 3/22/08 - working version after bug fixes


global p1name  p2name p3name
global X_column

% %--- set global printing options --
% set(0,'DefaultFigurePaperOrientation','portrait');
% 
%     %============= First plot: data presentation ==================
%     %-- Generate ideal Y data {'experimental') --
%     Y=feval(funct,params,X);
%     %-- Generate smooth curve for plotting ---
%     X_taucp=X(:,X_column);
%     x_smooth=[];
%     fields_smooth=[];
%     x_smooth_segment=[ min(X_taucp) : (max(X_taucp)-min(X_taucp))/smooth_x_intervals : max(X_taucp) ]';
%     field_segment=ones(length(x_smooth_segment), 1);
%     for i=1:length(fields)
%         x_smooth=[ x_smooth x_smooth_segment];
%         fields_smooth=[fields_smooth field_segment*fields(i)/fields(1)];
%     end
%     x_smooth=[x_smooth  fields_smooth];
%     y_smooth=feval(funct, params, x_smooth);
% 
%     HF1=figure('Name',project_name);
%     set(HF1,'Position',[200,0,700,700]);
%     set(HF1,'FileName',sprintf('%s.fig',project_name));
%     subplot(3,2,1);
%     title(sprintf('%s ; Y+%.2f%% noise ',project_name,noiseRMSDvalue*100*SSD_level_toshow ));
%     hold on;
%     plot(X_taucp,Y,'o'); xlabel('X'); ylabel('Y');
%     hold on;
%     plot(x_smooth(:,X_column),y_smooth,'r');
% 

% %========== find minSSD - sqrt (sum of squared) ==================
% %    if SSD < minSSD - such a  parameter set will be considered acceptable
% % 1. Generate randomized a Y sets
% % 2. Calculate SSD to ideal Y data
% % 3. Average SSD to obtain minSSD corresponding to particular
% % noiseRMSDvalue
% % 4. Repeat that for 2x and 3x random noise value
% disp(sprintf('Finding minSSD  %.2f%% noise RMSD in Y-values',100*noiseRMSDvalue));
% rndGauss = normrnd(0,1,length(Y),num_rnd_sets);
% minSSD=0;
% noise = Y*noiseRMSDvalue;  %--- calculate absolute noise RMSD values for particular Ys
% for i=1:num_rnd_sets
%     %-- apply noise to Y data
%     Ynoisy= rndGauss(:,i) .* noise + Y;
%     Ynoisy_toshow= rndGauss(:,i) .* noise*SSD_level_toshow + Y;
%     plot(X_taucp,Ynoisy_toshow,'.','MarkerSize',1);
%     %--- calculate SSD to ideal data 
%     SSD=sqrt(sum((Ynoisy-Y).^2));
%     minSSD=minSSD+SSD/num_rnd_sets;
% end
% disp(sprintf('minSSD(%.2f%%Y) = %.2e',100*noiseRMSDvalue,minSSD));
% noisepercent=noiseRMSDvalue*100;%-- to set relative values for contour levels
% 
HF1=figure; hold on
%-- display experimental data
start_index=1;
for dataset=1:length(res_data.indices)
    
    x=res_data.Xs{dataset};
    y=res_data.Ynormalized{dataset};
    %-- extract composite X-data for line shape simulation
    next_start_index=start_index+res_data.lengths(dataset);
    x_composite=res_data.Xcomposite(start_index : next_start_index-1,:);
   
    y_fit=model(params.param,x_composite);

    
        plot(x,y_fit,'k','LineWidth',1);

        plot(x,y,'bo',x,y,'b:','LineWidth',1); hold on;

    start_index=next_start_index;
 

end


%=================== Generate parameter arrays =======================
SSD_matrix_size=size(SSD_matrix);
for i=1:length(params.param)
    param_grid(i).range=params.range(i,:);
    param_grid(i).array = [ param_grid(i).range(1): params.grid_step(i) : param_grid(i).range(2)];
    param_grid(i).length = length( param_grid(i).array );   
    %--- check if that was correct SSD matrix loaded:
    if param_grid(i).length~=SSD_matrix_size(i)
        error('SSD_matrix_dimensions do not match size of your parameter grid!!!');
    end
end

%================ Scan parameter grid and generate SSD matrix ==============

switch length(params.param)
    case 2
        error('I did not do 2D. Recommendation: make you 2-parameter function quasi-3D.');
       
    %=================================== 3D ===========================================        
    case 3
            
        
         %==== creating report =======
        file=fopen(fullfile(DPM_path,sprintf('%s_params.txt',project_name)),'wt');
        fprintf(file,'Function: %s\n\n', funct);
        fprintf(file,'Parameter values:\n');
        fprintf(file,'%10s\t%10s\t%10s\n',p1name,p2name,p3name);
        fprintf(file,'%10.3e\t%10.3e\t%10.3e\n\n ', params.param(1),params.param(2),params.param(3));
        fprintf(file,'Grid scale:\n%10s\t%10s\t%10s\n\n ', params.funct{1},params.funct{2},params.funct{3});
        fprintf(file,'Grid steps:\n%10.3e\t%10.3e\t%10.3e\n\n ', params.grid_step(1),params.grid_step(2),params.grid_step(3));
        fprintf(file,'Parameter sets with SSD < minSSD*%.3f\n--------------\n%10s\t%10s\t%10s\t%10s\t%10s\n',...
            SSD_level_toshow,'Set #',p1name,p2name,p3name,'SSD');
        
        %============= plot of all the functions satisfying minSSD and save the parameters =========
 %       subplot(3,2,2);
 %       title(sprintf('SSD < SSD(%.2f%% noise RMSD )',noisepercent*SSD_level_toshow));
        hold on; 
    
        %---- collect parameter sets within SSD_levels from test set ---
        %--- do vectorized search ------
 
        % Storage for parameter sets:
        SSD_N=cell(length(SSD_levels),1);
        [I,J,K]=size(SSD_matrix);
        
        %-- find linear indices of SSDs to be shown and saved as equally
        %    good datasets 
        SSD_toshow_N=find(SSD_matrix<minSSD*SSD_level_toshow); 

        %--- find indices for different (next higher) level of SSD and put into separate categories ---
        SSD_N{1}=find(SSD_matrix<minSSD*SSD_levels(1));
        for level=2:length(SSD_levels)
                SSD_N{level}=find(SSD_matrix<minSSD*SSD_levels(level) & SSD_matrix > minSSD*SSD_levels(level-1));
        end
  
        %--- plot and save parameters that are considered to be equivalently good ---
        current_pars.param=zeros(3,1);
        for count=1:length(SSD_toshow_N)
            %-- find parameter indices --
            [i,j,k] = linear2triple_index(SSD_toshow_N(count),I,J,K);
            %-- get the parameters
            current_pars.param(1)= feval(params.funct{1},param_grid(1).array(i)) ;
            current_pars.param(2)= feval(params.funct{2},param_grid(2).array(j) ) ;
            current_pars.param(3)=  feval(params.funct{3},param_grid(3).array(k));
            
            
            %--- plot curves corresponding to these parameters
            
            start_index=1;
            for dataset=1:length(res_data.indices)

                x=res_data.Xs{dataset};
                y=res_data.Ynormalized{dataset};
                %-- extract composite X-data for line shape simulation
                next_start_index=start_index+res_data.lengths(dataset);
                x_composite=res_data.Xcomposite(start_index : next_start_index-1,:);

                y_fit=model(current_pars.param,x_composite);


                 plot(x,y_fit,'r','LineWidth',0.5,'DisplayName',sprintf('%d',count)); hold on;

                %    plot(x,y,'bo',x,y,'b:','LineWidth',1);

                start_index=next_start_index;
                


            end
           

          
            %Ynew_smooth=feval(funct, current_pars, x_smooth);
            
            %plot(x_smooth(:,X_column),Ynew_smooth,'k.', 'MarkerSize',1); hold on;
            
            %================  Save good parameter sets ==============
            %- this was a bug in original DPM
            fprintf(file,'%10d\t%10.3e\t%10.3e\t%10.3e\t%10.2e\n', count,...
                current_pars.param(1),current_pars.param(2),current_pars.param(3), SSD_matrix(SSD_toshow_N(count)) );
        end
        fclose(file);
        
        %--- plot test data on top
        %plot(X_taucp,Y,'o'); xlabel('X'); ylabel('Y');
        %hold on;

       %======== plot good parameter sets on maps ====
       if length(SSD_levels_colors)~= length(SSD_levels)
           error('SSD_levels_colors should have same number of elements as SSD_levels!');
       end

            %--- set limits ---
           % automatic
           p1max=max(param_grid(1).array);
           p1min=min(param_grid(1).array);
           p2max=max(param_grid(2).array);
           p2min=min(param_grid(2).array);
           p3max=max(param_grid(3).array);
           p3min=min(param_grid(3).array); 

%            % manual
%            p1max=3000;
%            p1min=0;
%            p2max=1;
%            p2min=0.75;
%            p3max=3000;
%            p3min=0; 

           HF3=figure; hold on 
           %============ plot p2 vs p1 ==============
           subplot(2,2,1);
           title('p2 vs p1');
           hold on;
           for level=length(SSD_levels):-1:1
               if ~isempty(SSD_N{level})>0
                   %--- extract values of parameters
                   [i,j,k]= linear2triple_index(SSD_N{level},I,J,K);
                   
                   p1=param_grid(1).array(i) ; 
                   p2=param_grid(2).array(j) ; 
                   %p3=param_grid(3).array(k);
                   
                  plot(p1,p2,SSD_levels_colors{level}, 'MarkerSize',10); hold on;
 %                  plot(p1,p2,SSD_levels_colors{level}); hold on;
                   
               end
           end
          
           set(gca,'XLim',[p1min p1max]);
           set(gca,'YLim',[p2min p2max]);
 
           %--- form a cross on the test parameter set ----
           h_x=[p1min p1max]; 
           point=feval(params.inverse_funct{2},params.param(2));
           h_y=[  point point ];
           plot(h_x,h_y,'r');
           
           v_y=[p2min p2max]; 

            point=feval(params.inverse_funct{1},params.param(1));
            v_x=[ point point ];
           plot(v_x,v_y,'r');
          xlabel(p1name); ylabel(p2name); %clabel(c,h);
 
  
          
                 %============ plot p3 vs p1 ==============
           subplot(2,2,2);
           title('p3 vs p1');
           hold on;
           for level=length(SSD_levels):-1:1
               if ~isempty(SSD_N{level})>0
                   %--- extract values of parameters
                   [i,j,k]= linear2triple_index(SSD_N{level},I,J,K);
                   
                   p1=param_grid(1).array(i) ; 
                   %p2=param_grid(2).array(j) ; 
                   p3=param_grid(3).array(k);
                   
                   plot(p1,p3,SSD_levels_colors{level}, 'MarkerSize',10); hold on;
                   
               end
           end
          
           set(gca,'XLim',[p1min p1max]);
           set(gca,'YLim',[p3min p3max]);
 
           %--- form a cross on the test parameter set ----
           h_x=[p1min p1max] ;
           point=feval(params.inverse_funct{3}, params.param(3));
           h_y=[  point point ];
           
           plot(h_x,h_y,'r');
           v_y=[p3min p3max] ;
           point=feval(params.inverse_funct{1}, params.param(1));
           v_x=[ point point];
           plot(v_x,v_y,'r');
          xlabel(p1name); ylabel(p3name);      
          
          
          
           %============ plot p3 vs p2 ==============
           subplot(2,2,3);
           title('p2 vs p3');
           hold on;
           for level=length(SSD_levels):-1:1
               if ~isempty(SSD_N{level})>0
                   %--- extract values of parameters
                   [i,j,k]= linear2triple_index(SSD_N{level},I,J,K);
                   
                   %p1=param_grid(1).array(i) ; 
                   p2=param_grid(2).array(j) ; 
                   p3=param_grid(3).array(k);
                   
                   plot(p3,p2,SSD_levels_colors{level}, 'MarkerSize',10); hold on;
                   
               end
           end
          
           set(gca,'XLim',[p3min p3max]);
           set(gca,'YLim',[p2min p2max]);
 
           %--- form a cross on the test parameter set ----
           h_x=[p2min p2max]; 
            point=feval(params.inverse_funct{3}, params.param(3));
          h_y=[ point point ];
           plot(h_y,h_x,'r');
           v_y=[p3min p3max]; 
            point=feval(params.inverse_funct{2}, params.param(2));
           v_x=[ point point ];
           plot(v_y,v_x,'r');
          xlabel(p3name); ylabel(p2name);        
          
          
          
          
          

          	
            %--- just the text data ----
            subplot(2,2,4);
            text(0.05,0.95,'Test parameters / grid step:');
            text(0.05,0.6,sprintf('%s= %.3g / %s %.3g\n%s= %.3g / %s %.3g\n%s= %.3g / %s %.3g',...
                p1name, params.param(1),  params.funct{1}, params.grid_step(1),...
                p2name, params.param(2),  params.funct{2}, params.grid_step(2),...
                p3name, params.param(3),  params.funct{3}, params.grid_step(3)));
            text(0.05, 0.2, sprintf('Grid sizes: '));
             text(0.05, 0.1, sprintf('%d x ', size(SSD_matrix)));
         
           

       
       
       %================= Save figures ===========================
          %saveas(HF1,sprintf('%s.png',project_name));
          saveas(HF1,fullfile(DPM_path,sprintf('%s_minSSD_sets.fig',project_name)));
          saveas(HF3,fullfile(DPM_path,sprintf('%s_DPM.fig',project_name)));
          saveas(HF1,fullfile(DPM_path,sprintf('%s_minSSD_sets.png',project_name)));
          saveas(HF3,fullfile(DPM_path,sprintf('%s_DPM.png',project_name)));

    otherwise
        error(sprintf('I cannot do %d dimensions!',length(params.param)));
end
        




