function DPM_run(mode, res_data, best_fit_beta, minSSD, params, DPM_path)
% This function runs DPM calculations in LineShapeKin module
%
%  Input 
%   mode     -  computation or analysis - to allow for long (overnight)
%                     computations of SSD matrices
%   res_data - structure of the spectral raw and normalized data for the
%                     specific peak
%   best_fit_beta - the parameters from fitting
%   minSSD             - sum of squared residuals from fitting to judge
%                           similarity to other parameter sets
%   params     - parameter structure created in DPM_setup.m
%                   
%   figures_path - where to dump resulting figures

% ----- - old header ------
%
%--- Analysis of correlations between parameters of fitting functions ---
% by Evgenii Kovrigin   07/14/05
% 
% Brute force search for existence of parameter combinations that generate
% function within certain tolerance from data correspoinding to test 
% parameters.
%
%----- 
%  7/20/05      Working version. Reproduces Pat's results.

%  8/02/05      Grid scanning

%  8/05/05      Separated matrix computation and plotting

%  8/8/05        Vectorized matrix analysis for speed (more than 20x
%  acceleration)

% 8/9/05    Latest version: additional comments added 

% 8/23/05  The bug in R2 (/ for ./ ) fixed. Calculations for the paper
%                repeated

% 8/23/05  Manual plotting mode added: to plot interesting parameter sets
%
% 3/21/08 - Orthogonal X matrix now for compatibility with LineShapeKin

% =============== LINESHAPEKIN application ====================
%
% 3/21/08 Direct parameter mapping of line-shape analysis procedure
% (LineShapeKin ver. 3.0)
%
% 3/22/08 - working version after bug fixes
%
% 3/26/08  - now it computes SSD to EXPERIMENTAL DATA rather than to Y DATA
%                   produced with the model parameter set. This tells you
%                   what are other parameter sets that will describe you
%                   experimental data as good as the best-fit result.



global p1name  p2name  p3name X_column

global YES NO Model B_vector X_columns 


%========= INPUT YOUR PARAMETERS HERE ===============




%--- Note on meaning of 'test parameters': 
% We will analyze the parameter space for all the parameter sets
% that produce SSD close to the one produced by 'test parameters'.
% Results of this testing depend both on the values of 'test parameters'
% and on X variable as well

%-- 'Test parameters' as  A0, B0 etc  --

p1name=B_vector.parameter_names{1};  
p2name=B_vector.parameter_names{2};
p3name= B_vector.parameter_names{3};



%--- Enter here parameter values (a single value is OK)
%     NOTE: Put here real values of parameters regardless of the conversion
%     function you will put for grid search.
p1= [  best_fit_beta(1)     ];
p2= [  best_fit_beta(2)  ];
p3= [  best_fit_beta(3)  ]; 

SSD_levels= params.SSD_levels;  
SSD_levels_colors= params.SSD_levels_colors;  
SSD_level_toshow= params.SSD_level_toshow ;



%================   END OF USER INPUT ===================



%% ======= FIXED OPTIONS IN THIS PROJECT ======

%--- common prefix for output files ----
project_name=res_data.name;

%--  X data column-vector  ---
X=res_data.Xcomposite;  
%-- position of the column with actual X-values
X_column= X_columns.LPratio;


%--- This is an option to test a number of different combinations o
%       parameters. However, this does not work with LineShapeKin because
%       here we judge the quality by comparing to experimental data and
%       real SSD.
%
% p1= [    500 1000 2000 ];
% p2= [  0.75   0.95  0.99 ];
% p3= [  500 1000 2000  ]; 


%--- name of 'wrapper' function that calls 'model.m' and adapts parameters
%       for it.
funct='DPM_model';

%-- add field strength (may be just one field)---
fields=[500];

%--- number of intervals for smooth function plot presentation ---
smooth_x_intervals=100;

%--- number of synthetic sets for SSD determination
num_rnd_sets=100;

%========= END OF FIXED OPTIONS ===========








%-- calculate step of the grid for parameter space sampling ---                 
params.grid_step = [  ...
    (params.range(1,2)-params.range(1,1))/params.grid_points(1);...
    (params.range(2,2)-params.range(2,1))/params.grid_points(2);...
    (params.range(3,2)-params.range(3,1))/params.grid_points(3)
                                ]; 



%--- prepare X data like: 
%  taucp1   500
%  taucp2   500
%  taucp1   800
%  taucp2   800
   

%== show real parameter values to use
disp('Real parameter values to use')
format long e
disp(p1name);
param_num=1;
grid1=feval(params.funct{param_num}, ...
    [params.range(param_num,1):params.grid_step(param_num):params.range(param_num,2) ]);
disp(sprintf('%e\n',grid1 ));
disp(p2name);
param_num=2;
grid2=  feval(params.funct{param_num}, ...
    [params.range(param_num,1):params.grid_step(param_num):params.range(param_num,2) ]);
disp(sprintf('%e\n', grid2));
disp(p3name);
param_num=3;
grid3= feval(params.funct{param_num}, ...
    [params.range(param_num,1):params.grid_step(param_num):params.range(param_num,2) ]) ;
disp(sprintf('%e\n', grid3 ));
N1=length(grid1);
N2=length(grid2);
N3=length(grid3);

disp(sprintf('Grid size: %d x %d x %d = %d total points', N1,N2,N3, N1*N2*N3));


%--- add field dependence 

[N M]=size(X);
X_fields=[];
X_taucp=[];
for i=1:length(fields)
    X_field=ones(N,1)*fields(i)/fields(1);
    X_fields= [X_fields; X_field];
    X_taucp= [X_taucp; X];
end
X = [ X_taucp   X_fields];


%--- run analysis for every test point in the parameter space

disp(sprintf(' ========== Current run:   %s  ==============',mode));

for p1_counter=1:length(p1)
    for p2_counter=1:length(p2)
        for p3_counter=1:length(p3)
            %--- generate the name for the parameter set ----

            set_name=sprintf('%s.%s_%d.%s_%d.%s_%d',...
                            project_name,p1name,p1_counter,  p2name,p2_counter, p3name,p3_counter);
            disp(sprintf(' ------------- Test set %s ------------------',set_name));

            %-- make file names
            SSD_matrix_fname=fullfile(DPM_path,sprintf('%s_SSD_matrix.mat',set_name));
            

            %--- set current test parameters ----------
            params.param= [     p1(p1_counter)   ;     p2(p2_counter)    ;  p3(p3_counter)        ]; 
            
            switch mode
                %--------------------- SSD COMPUTATION -----------------
                case 'computation'
                        %-- extract experimental data
                        Y_experimental = res_data.Ytotal;
                        %-- compute SSD to experimental data for all
                        %    datasets on the grid
                        matrix_computation(funct, params, X, Y_experimental)

                        %--- rename saved data ---
                        movefile('SSD_matrix.mat',SSD_matrix_fname);


                 %----------------- ANALYSIS AND PLOTTING --------------
                case 'analysis'
                    %--- load saved data ----
                    try
                        load(SSD_matrix_fname);
                    catch
                        error('Error reading the data from %s!', SSD_matrix_fname);
                    end
                  
                   %--- data plotting and analysis ----
                    matrix_analysis_v3(res_data, minSSD, set_name, funct, params, X, fields, ...
                            smooth_x_intervals, num_rnd_sets,...
                            SSD_matrix, SSD_levels, SSD_level_toshow, SSD_levels_colors, DPM_path)
                            
                            

                 
            end
                    
            disp('Data files were written!');
        end
    end
end
             
             

 
 
 