function [res_data] = prepare_residue_data(Sparky_data_path, cur_resname, conditions, modes)
%
% 3/11/08
% Written by Evgenii Kovrigin, Medical College of Wisconsin
% by Evgenii Kovrigin 06/29/05
%
% - Structure 'conditions' passes concentrations, L/P ratio, temperature
% as well as starting parameter estimates for line widths
%       conditions.Ptotal, .PLratio, .Temperature, .LineWidth0
%
%
% - This function generates normalizes Y data for fitting and obtains 
% peak positions and start/end R2 if modes.estimate.peak_position=YES 
% and estimate.R2=YES set. 
% If 'NO' then it just fills the corresponding fields in the output data
% array with  numbers provided in modes.preset. w0,wN, R20, R2N
%
% - Normalization is done by either direct trace integration of  fitting 
% to lorentzian lineshapes (modes.normalization='integration'/'fit'). 
%
% Constant baseline offset
% correction is optional (modes.baseline_correction=YES/NO)
% 
% Plotting of raw and prepared data: modes.plot_raw_data=YES/NO
%
% Returns structure 'res_data' for easy access to raw data 
%    res_data fields
%
%           name - residue name
%           indices - array of indices of titration points corresponding 
%                       to ligand concentrations in main array
%                       (we may have some titration points dropped and then using index
 %                          we can extract corresponding concentrations)
%           Xs       -  cell array of X arrays, each one is a column (sequential order as read from the data folder.
%                          NOT the sequential titration order! That one is set by indices.)
%           Ys        - cell array of Y arrays 
%           
%           lengths - lengths of individual X arrays
%
 %         Ynormalized - cell array of arrays of normalizes intensity data
 %                                (in the order of reading off the disk)
 %         Xcomposite - X array containing frequency as well as relaxation
 %         rates, endpoint frequencies, concentrations etc
 %
 %         Ytotal - Y array for all the slices
%      
% ----------------------------
%  3/29/08 - index is added as an additional column to distinguish titration
%                   points
%               - trace number added to plot output


global  YES NO  X_columns ;
R2str='';
w0str='';

%--- messages ---
disp(sprintf('Setting up datasets for %s', cur_resname));
if modes.baseline_correction==YES
    disp('Doing constant baseline offset correction');
end
switch modes.normalization
    case 'fitting'
                disp('Doing normalization by LORENTZIAN FITTING (fast exchange limit)');
    case 'integration'
            disp('Doing normalization by INTEGRATION (isolated peak condition)');
     otherwise
        error('Wrong type of normalization: %s !!!', modes.normalization);
end
   
% Load the data 
[res_data presets ]= load_raw_residue_data(Sparky_data_path,...
    cur_resname, conditions);
% add storage for normalized data
res_data.Ynormalized=cell(1,length(res_data.indices));


%-- prepare figure 
if modes.plot_raw_data==YES
   Hf=figure;
   graph_title=sprintf('%s_normalization',cur_resname);
   set(Hf,'Filename',graph_title);
   hold on;

   if modes.baseline_correction==YES
       str=', baseline-offset corrected data (red)';
   else
       str='';
   end
   subplot(2,1,1);
   title(sprintf('%s : Raw data (blue) %s',cur_resname, str));
   subplot(2,1,2);
   title(sprintf('Normalized by %s; Parameter estimates: A (black), B (red)', modes.normalization));
   
end

%-- prepare room for the datasets
Ptotal_column=[];
LPratio_column=[];
R2A_column=[];
R2B_column=[];
wA_column=[];
wB_column=[];
Temp_column=[];
DataID_column=[];
Ytotal=[];
Xtotal=[];


switch modes.normalization

            
    case 'integration'
        %--- Normalization for WELL RESOLVED slow and intermediate exchangers
        %-- Assume that peak tails go to 0 on both sides and (optionally) subtract residual height as 
        %     a constant baseline offset.
        
       %--- what I can add here: calculate integrals for all of them,
       % AVERAGE and only then normalize. This way I will have fewer OUTLIER
       % 1D slices
            for i=1:length(res_data.indices)
                 
                x_current=res_data.Xs{i};
                y_current=res_data.Ys{i};
                y_len=res_data.lengths(i);

                %-- simple zero-order approximation ---
               baseline_0=min(y_current);
 
                if modes.baseline_correction==YES
                    %--- subtract baseline ---
                    y_corr=y_current - baseline_0;
                else
                    y_corr=y_current;
                end

                %-- Calculate integral of the line shape
                dx=x_current(2)-x_current(1);
                
                switch modes.integration_mode
                    
                    
                    
                    case 'rectangular'
                        %-- integrate by simple rectangular integration ---
                        Integral= sum(y_corr)*dx;
                        
                        
                        
                    case 'simpson'
                        %-- integration by Simpson's formula
                        
                        Simpsons_row=zeros(y_len,1);
                        for counter=1:y_len
                            if counter==1 || counter==y_len
                                Simpsons_multiplier=1;
                           elseif mod(counter,2)==0
                                Simpsons_multiplier=4;
                            else
                                Simpsons_multiplier=2;
                            end
                            
                            Simpsons_row(counter)=Simpsons_multiplier;
                        end
                        
                        Integral= sum(y_corr.*Simpsons_row)*dx/3;
                        
                        
                    otherwise
                        error('Wrong integration mode!')
                end
                
                %-- normalize to unity
                y_norm=(y_corr)/Integral;
                
                %-- scale intensities to compensate for incomplete peak
                %    envelope
                y_norm=y_norm/presets.scale_factor;
                

                %-- visualize raw and corrected data ---
                
                if modes.plot_raw_data==YES
                    subplot(2,1,1);
                    xlabel('Frequency, Hz');
                    hold on;
                    plotname=sprintf('Trace %d', res_data.indices(i));
                    plot(x_current,y_current,'bo','DisplayName',plotname);
                    plot(x_current,y_current,'b');
                    hold on;
                    if modes.baseline_correction==YES
                       plotname=sprintf('Trace %d - baseline', res_data.indices(i));
                       plot(x_current,y_corr,'r','DisplayName',plotname);
                    end
                    subplot(2,1,2);
                    hold on;
                    plotname=sprintf('Trace %d - normalized', res_data.indices(i));
                    plot(x_current,y_norm,'g','DisplayName',plotname);
                end

                %--- scale and save  --
                res_data.Ynormalized{i}=y_norm;

                %-- form data set for fitting
                Ytotal=[ Ytotal; y_norm ];
                Xtotal=[ Xtotal; x_current ];
                
                % -- generate columns of values that are specific to a 1D
                % slice (not sequential from conditions arrays but set by
                % indices - in case when some points are dropped from the
                % middle. 
                Ptotal_column    =[ Ptotal_column  ;  ones(y_len,1) *  conditions.Ptotal(res_data.indices(i))];
                LPratio_column  =[ LPratio_column ; ones(y_len,1) *  conditions.LPratio(res_data.indices(i))];
                Temp_column    =[ Temp_column    ; ones(y_len,1) *  conditions.Temperature(res_data.indices(i))];
                DataID_column    =[ DataID_column    ; ones(y_len,1) *  res_data.indices(i)];  %-- remembers which titration point it was
                 
               
            end

            

            %  prepare initial estimates
            dw0=conditions.LineWidth0;
            %--- take first peak ------
            x_A=res_data.Xs{1};
            y_A=res_data.Ynormalized{1};
            [maxval_A, I]=max(y_A);
            %-- fitting ----
            b_A=[maxval_A/dw0,dw0,x_A(I)];
            
            
            %--- take last peak ------
            x_B=res_data.Xs{length(res_data.indices)};
            y_B=res_data.Ynormalized{length(res_data.indices)};
            [maxval_B, I]=max(y_B);
            %-- fitting ----
            b_B=[maxval_B/dw0,dw0,x_B(I)];
            
            LARGE=1e16;  %- just a large number for maximum limits of estimated values

            
            
            %--- check if w0 and R2s were preset
            if     isempty(presets.R2s) || isempty(presets.w0s) %-- w0 or/and R2s are not preset. Go estimate them
               
                    % -- estimate R2s and v0s assuming that first spectrum  is
                    % one peak at a position of A form and the last
                    % spectrum is pure B form

                    disp('Estimate starting parameters using first titration point as pure state A and last point as pure state B.');
                   
                    
                     %--- set bounds and options
                    %lower_bounds=zeros(1,length(b_A));
                      lower_bounds=[0 0 min(x_A)];
                      upper_bounds=[LARGE LARGE max(x_A)];

                      options=optimset('Diagnostics','off', ...
                          'Display','off',...
                          'LargeScale','on',...
                          'LevenbergMarquardt', 'on',...
                          'TolFun',1e-12,...
                          'TolX',1e-12,...
                         'MaxFunEvals', 1e6);
                    
                   
                     %-- fit first peak
                     [b_A, resnorm, residual, exitflag, output,lambda, Jacobian]=...
                          lsqcurvefit(@classic_lorentzian, b_A, x_A, y_A,...
                          lower_bounds,upper_bounds,options);

                    %-- calculate parameters ----
                    R2_A=pi*b_A(2); 
                    w0_A=b_A(3);

                    
                    %-- fit last peak
                      lower_bounds=[0 0 min(x_B)];
                      upper_bounds=[LARGE LARGE max(x_B)];

                    
                     [b_B, resnorm, residual, exitflag, output,lambda, Jacobian]=...
                          lsqcurvefit(@classic_lorentzian, b_B, x_B, y_B,...
                          lower_bounds,upper_bounds,options);

                   %-- calculate parameters ----
                    R2_B=pi*b_B(2); 
                    w0_B=b_B(3);
                    
            
            end

            
            %- reset R2 if preset
            if ~isempty(presets.R2s)  % - R2 were supplied. Use them and adjust intensities to the same height
                disp('Resetting R2 to requested values and simply scale intensities')
                % scale intensities
                old_b_A2=b_A(2);
                old_b_B2=b_B(2);
                % 
                R2_A=presets.R2s(1);
                b_A(2)=R2_A/pi;
                R2_B=presets.R2s(2);
                b_B(2)=R2_B/pi;
                % scale them back now
                b_A(1)=(b_A(1)/old_b_A2)*b_A(2);
                b_B(1)=(b_B(1)/old_b_B2)*b_B(2);
              
                R2str='(fix)';
            end

            %-- reset w0 if preset
            if  ~isempty(presets.w0s)   % if w0 were preset - use them instead of optimized
                disp('Resetting w0 to requested values')
                w0_A=presets.w0s(1);
                b_A(3)=w0_A;
                w0_B=presets.w0s(2);
                b_B(3)=w0_B;
                
                w0str='(fix)';
            end
            
               
                
                   
    otherwise
        error('Wrong type of normalization: %s !!!', NORMALIZATION);
end






%-- visualize line shapes for starting and ending
%     parameters
if modes.plot_raw_data==YES
    dx_A=(max(x_A)-min(x_A))/100;
    x_A_smooth=[min(x_A) : dx_A : max(x_A)]';
    y_A_smooth=classic_lorentzian(b_A,x_A_smooth);
    
    dx_B=(max(x_B)-min(x_B))/100;
    x_B_smooth=[min(x_B) : dx_B : max(x_B)]';
    y_B_smooth=classic_lorentzian(b_B,x_B_smooth);

    
    
    %-- normalize line shapes if the R2 or w2 were preset
     if     ~isempty(presets.R2s) || ~isempty(presets.w0s)
         
            Integral_A = sum(y_A_smooth)*dx_A*presets.scale_factor;
            y_A_smooth=y_A_smooth/Integral_A;

            Integral_B = sum(y_B_smooth)*dx_B*presets.scale_factor;
            y_B_smooth=y_B_smooth/Integral_B;
     end
    
    
    
    subplot(2,1,2);
    plotname='A state estimated';
    plot1=plot(x_A_smooth,y_A_smooth,'k','DisplayName',plotname); 
    plotname='B state estimated';
    plot2=plot(x_B_smooth,y_B_smooth,'r','DisplayName',plotname); 

    annotation(Hf,'textbox','String',...
        sprintf('Estimates:  w0_A%s = %.3f Hz   R_{2,A}%s = %.3f s^{-1}    w0_B%s =%.3f Hz   R_{2,B}%s = %.3f s^{-1}', ...
        w0str,w0_A,R2str,R2_A,w0str,w0_B,R2str,R2_B),...
    'FitBoxToText','off','LineStyle','none',...
    'Position',[0.01855 0.02609 0.9499 0.0442]);

end
    



hold off;

if modes.plot_raw_data==YES
    %-- save figure
    fname=fullfile(conditions.figures_path,graph_title);
    saveas(Hf,  fname, 'fig');
    saveas(Hf,  fname, 'png');
end

%% Prepare main 'X' dataset for fitting


%-- generate columns of values that are constant for the spin system
%     througout the titration experiment
y_total_len=length(Ytotal);
R2A_column       = ones(y_total_len,1) * R2_A;
R2B_column        = ones(y_total_len,1) * R2_B;
wA_column         = ones(y_total_len,1) * w0_A;
wB_column         = ones(y_total_len,1) * w0_B;


%%- fill in main experimental data array

res_data.Xcomposite=ones(y_total_len, X_columns.number);

res_data.Xcomposite(:,X_columns.Ptotal)=Ptotal_column;
res_data.Xcomposite(:,X_columns.LPratio)=LPratio_column;
res_data.Xcomposite(:,X_columns.Xtotal )=Xtotal ;
res_data.Xcomposite(:,X_columns.R2A)=R2A_column;
res_data.Xcomposite(:,X_columns.R2B)=R2B_column;
res_data.Xcomposite(:,X_columns.wA)=wA_column;
res_data.Xcomposite(:,X_columns.wB)=wB_column;
res_data.Xcomposite(:,X_columns.Temp)=Temp_column;
res_data.Xcomposite(:,X_columns.DataID)=DataID_column;




res_data.Ytotal=Ytotal;
