function [beta_opt,beta_error_opt, resnorm_opt]=fit_run(resStr, beta0, Xfull, Yfull, conditions, modes)
%------------------- LineShapeKin FIT_ALL_DATA -----------------------------
%-- version 3.0 
%
% 3/17/08
% Written by Evgenii Kovrigin, Medical College of Wisconsin

%-- Fitting a dataset to the model performing a 
%- simple grid-search and plotting  the result --
%  
%
%- Parameters:
%    resStr - formal name of the data set
%    Xfull  - concatenated dataset of experimental variables
%    Yfull  - concatenated dataset of peak intensities
%    conditions - structure of experimental conditions
%    modes       - structure of control parameters (fitting, display, model
%                        switching modes)
%
% 	estimate_error - YES/NO  -  a flag to skip error determination when doing a grid search 
%                             just for speed of calculations (if you want to)
%                             
% 	Return parameters: 
%                   (1) optimized vector of function parameters 
%                   (2) vector of their errors
%                   (3) norm - sum of squares from fitting
% 


global YES NO B_vector Model;


disp(sprintf('\nFitting  to <%s>', Model.name));
 
%--- create arrays ---
grid_N = length(B_vector.grid_values);
pars_N = length(beta0);


beta_error=zeros(pars_N,1); %-- error values

beta_grid = zeros(pars_N,grid_N);     %-- grid search results
beta_error_grid = zeros(pars_N,grid_N);     %-- grid search errors

resnorm_grid = zeros(1,grid_N); %-- array of resnorm values of grid search

%--- Define what mode is in use: grid search or final optimization?
if modes.grid_search==YES
    disp(sprintf('Grid search:'));
    passes_N=grid_N;
    %--- pretty display of Koffs in a grid ---
    disp(sprintf('Current values in a grid search:'));
    outstring='';
    for grid_element=1:length(B_vector.grid_values)
        outstring=sprintf('%s %.2f\t',outstring,B_vector.grid_values(grid_element));
    end
    disp(outstring);
    disp(sprintf('\nGrid search parameter: [%d]  %s', B_vector.grid_search_parameter, ...
        B_vector.parameter_names{B_vector.grid_search_parameter}));
else %-- final optimization mode - one fitting pass starting from optimized values
    disp(sprintf('\nFinal optimization:'));
    passes_N=1;
end


%% Main fitting loop    
for grid_element = 1:passes_N

  if modes.grid_search==YES 
              %--- assign starting values for variable parameters ---
              beta0(B_vector.grid_search_parameter)=B_vector.grid_values(grid_element);
              disp(sprintf('\nSearch %d',grid_element));
  end

  
  
  
  
  %------- use LSQCURVEFIT -------------
  lower_bounds=B_vector.min';
  upper_bounds=B_vector.max';
  
  options=optimset('Diagnostics','off', ...
      'LargeScale','on',...
      'LevenbergMarquardt', 'on',...
      'TolFun',1e-12,...
      'TolX',1e-12,...
     'MaxFunEvals', 1e6);
  [beta, resnorm, residual, exitflag, output,lambda, Jacobian]=...
      lsqcurvefit(@model, beta0, Xfull, Yfull,...
      lower_bounds,upper_bounds,options);

  
% %   %-- test on resnorm computation: test passed, norm computation is correct
%   
%   y_fitted= model(beta, Xfull);
%   new_norm=sum( (y_fitted-Yfull).^2 );
%   disp(sprintf('TESTING: my norm is %e', new_norm))
% keyboard;
  
  
  % --- Exitflag interpretation ----
  if exitflag > 0
      disp('LSQCURVFIT converged to a solution');
  elseif exitflag < 0
      disp('LSQCURVFIT DID NOT converge to a solution!');
  else
      disp('Maximum number of iterations reached!');
  end
  
  %--- Skip this section for speed of grid search ---
  if modes.estimate_error==YES
      disp('Determination of parameter uncertainties...');
      ci=nlparci(beta,residual,Jacobian);
      beta_error=ci(:,2)-beta;
  else 
      disp('Skipping error determination');
  end
  


  %--- save norm of residuals, parameters and their errors
  resnorm_grid(grid_element)=resnorm;
  beta_grid(:,grid_element)=beta;
  beta_error_grid(:,grid_element)=beta_error;

  % report results at current grid point.
  disp(sprintf('Norm (sum of squared residuals): %.2e',resnorm));
  for  beta_element=1:length(beta);
       if modes.estimate_error==YES
             disp(sprintf('[%d] %s:  %.5e  +/-  %.1e',  beta_element, B_vector.parameter_names{beta_element},...
                 beta( beta_element), beta_error( beta_element)));
       else
             disp(sprintf('[%d] %s:  %.5e ',  beta_element, B_vector.parameter_names{beta_element},...
                 beta( beta_element)));
       end
  end
  
%- end of a loop over a grid
end
%


    

 
%----- Export grid results to a plain text file -------
fName = sprintf('%s.fit.txt',resStr);

if modes.grid_search==YES   %--- report grid search results
 
    %--- Find minimal resnorm values and corresponding fit results
    [resnorm,imin] = min(resnorm_grid);
    resnorm_opt = resnorm_grid(imin);
    beta_opt = beta_grid(:,imin);  
    beta_error_opt = beta_error_grid(:,imin);

    disp('Minimal norm result:')
     for  beta_element=1:length(beta);
          disp(sprintf('[%d] %s:  %.3e  ',  beta_element, B_vector.parameter_names{beta_element},...
              beta_opt( beta_element)));
     end

    
    
    %---- First run: reporting  Grid search results ----
	fp = fopen(fullfile(conditions.fit_results_path,fName),'w');
	fprintf(fp,'Residue: \t%s\n', resStr);
    fprintf(fp,'Grid search parameter: [%d]  %s\n\n', B_vector.grid_search_parameter, ...
        B_vector.parameter_names{B_vector.grid_search_parameter});
    

	for grid_element=1:grid_N
           %--- Just pretty mark at minimum line ---
           if grid_element==imin
               special_mark='<----------------------- minimum';
           else
               special_mark='';
           end
        %-- print the current grid parameters to the file
        fprintf(fp,'----------------------\nNorm(%d): %.2e  %s \n', grid_element,  resnorm_grid(grid_element), special_mark);
        for beta_element=1:length(beta);
               fprintf(fp, '\t[%d] %s:  %.5e\n', beta_element, B_vector.parameter_names{beta_element},...
                  beta_grid(beta_element, grid_element));
        end

	end
	
	fclose(fp);
else
        %--- Second run (single fitting session): determination of errors of optimized values and appending results to the file---

        %-- take results
        resnorm_opt = resnorm;
        beta_opt = beta;  
        beta_error_opt = beta_error;

        
        
        
        fp = fopen(fullfile(conditions.fit_results_path,fName),'a');
        fprintf(fp,'\n\n\nOptimum norm: %.2e\n', resnorm_opt);
        for beta_element=1:length(beta);
           fprintf(fp, '\t[%d] %s:  %.5e  +/-  %.1e\n',  beta_element, B_vector.parameter_names{beta_element},...
               beta_opt(beta_element), beta_error_opt(beta_element));
        end
        fclose(fp);
        

        %--- save sum of squares of residual for model comparison as:
        %----    Number_of_experimental_titration_points   Sum_of_squares_of_residuals
       
        fName2=sprintf('%s.SS.txt', resStr);
        fp2 = fopen(fullfile(conditions.SS_path,fName2),'w');
        fprintf(fp2,'%d %e\n', length(conditions.LPratio), resnorm); 
        fclose(fp2);
        
        
end

