% script for testing levenberg-marquardt method using the resid function.
% Steps:
%	-	read in data from daily reference and absolute reference files
%	-	shift/stretch - ie use lsqnonlin with marq method and resid function

global resid_plot_results	%temp

global lowest
global ss_count 

% Globals from INIT.M:
global n_abs n_pix abs_gas filenames dir

% Read in absolute reference spectrum:
disp(strcat('Absolute reference spectrum:', std_name))
% Initialize array:
ref = zeros(1,n_pix);
% Read in the absolute reference spectrum:
ref = read_spectrum(strcat(dir.ip, std_name));
% Verify that the correct number of points were read in:
if length(ref) ~= n_pix
   disp('Hey! something wrong with # of pixels in absolute ref spectrum')
end

% Read in daily reference spectrum:
disp(strcat('Daily reference spectrum:', file_nms(lowest,1:12)))
% Initialize array:
sp = zeros(1,n_pix);
% Read in the daily reference spectrum:
sp = read_spectrum(strcat(dir.data, file_nms(lowest,1:12)));
% Verify that the correct number of points were read in:
if length(sp) ~= n_pix
   disp('Hey! something wrong with # of pixels in daily ref spectrum')
end

% Now we try using the marq method to optimize stretch and shift:
shift = 0.5;
stretch = 4.0;
ss(1) = shift;
ss(2) = stretch;
cba_sp = [0 b_sp a_sp];			% a_sp, b_sp read in from A_B.REF
cba_ref = [0 b_ref a_ref];		% a_ref, b_ref read in from A_B.REF

global wv_o3 wv_no2 wv_oclo wv_o4 wv_h2o wv_bro wv_ring wv_no3
global sm_o3 sm_no2 sm_oclo sm_o4 sm_h2o sm_bro sm_ring sm_no3
% Give index numbers of gases to include in the fit:
abs_fit = [1 2 3 4 5 6 7];
load(strcat(dir.cs_op, 'cs_o3.mat'));
load(strcat(dir.cs_op, 'cs_no2.mat'));
load(strcat(dir.cs_op, 'cs_oclo.mat'));
load(strcat(dir.cs_op, 'cs_o4.mat'));
load(strcat(dir.cs_op, 'cs_bro.mat'));
load(strcat(dir.cs_op, 'cs_h2o.mat'));
load(strcat(dir.cs_op, 'cs_ring.mat'));

return	% debugging, Aug 21

% test resid:
%res = resid(ss, sp, ref, cba_sp, cba_ref, abs_fit)
%return

% Apply Marq method:
ss_UB = [1 2];			% upper bound on ss
ss_LB = [-1 0];		% lower bound on ss
% note, if we are using fsolve, upper and lower bounds are not used
% Create global counter to keep track of shift/stretch modifications:
ss_count = 0;
clear global ss_record
global ss_record

%option_Marq = optimset('LevenbergMarquardt', 'on');
option_Marq = optimset('display','final','LargeScale','on', ...
   'TolX',1e-6,'TolFun',0.05,'DerivativeCheck','off',...
   'Jacobian','off','MaxFunEvals','100*numberOfVariables',...
   'Diagnostics','off',...
   'DiffMaxChange',1e-1,'DiffMinChange',1e-8,...
   'PrecondBandWidth',0,'TypicalX','ones(numberOfVariables,1)',...
   'MaxPCGIter','max(1,floor(numberOfVariables/2))', ...
   'TolPCG',0.1,'MaxIter',400,'JacobPattern',[], ...
   'LineSearchType','quadcubic','LevenbergMarq','on'); 

%[ss_new, chisq] = lsqnonlin(resid(ss, sp, ref, cba_sp, cba_ref, abs_fit),...
%   ss, ss_LB, ss_UB, option_Marq);
%[ss_new, chisq] = fsolve(resid(ss, sp, ref, cba_sp, cba_ref, abs_fit),...
%   ss, ss_LB, ss_UB, option_Marq);

%[ss_new, chisq] = fsolve(resid(ss, sp, ref, cba_sp, cba_ref, abs_fit),...
%   ss, option_Marq);
%[ss_new, chisq] = fsolve(resid(ss, sp, ref, cba_sp, cba_ref, abs_fit),...
%   ss, option_Marq, sp, ref, cba_sp, cba_ref, abs_fit);
%[ss_new, chisq, exitflag] = fsolve(resid(ss, sp, ref, cba_sp, cba_ref, abs_fit),...
%   ss, option_Marq, sp, ref, cba_sp, cba_ref, abs_fit);

% HEY!!! This code actually works: (!!)
resid_plot_results = 0;
f = inline('resid(ss, sp, ref, cba_sp, cba_ref, abs_fit)', 'ss',...
   'sp', 'ref', 'cba_sp', 'cba_ref', 'abs_fit');
[ss_new, chisq, exitflag] = fsolve(f,...
   ss, option_Marq, sp, ref, cba_sp, cba_ref, abs_fit);
% Try lsqnonlin instead of fsolve:
%[ss_new, chisq, exitflag] = lsqnonlin(f,...
%   ss, option_Marq, sp, ref, cba_sp, cba_ref, abs_fit);
% Although, matrix becomes singular. So might not be doing what I think it is.
% But this is progress! (July 19)

% note, matrix is badly scaled (july 19)

% Show user the results:
disp('New shift and stretch:')
disp(strcat('	shift = ', num2str(ss(1)), ', stretch = ', num2str(ss(2))))
disp(strcat('Chi squared = ', num2str(chisq)))
if exitflag > 0
   disp('Fit was convergent')
elseif exitflag == 0
   disp('Max # of function evaluations reached - fit did not converge')
else
   disp('Fit was divergent')
end

% plot results:
resid_plot_results = 1;
res = resid([0 1], sp, ref, cba_sp, cba_ref, abs_fit);