% Levenberg-Marquardt example/test

%   Henri Gavin, Dept. Civil & Environ. Engineering, Duke Univ. 2 May 2016
%   modified from: ftp://fly.cnuce.cnr.it/pub/software/octave/leasqr/
%   Press, et al., Numerical Recipes, Cambridge Univ. Press, 1992, Chapter 15.

clear all
hold off
randn('seed',0);	% specify a particular random sequence for msmnt error

% *** For this demonstration example, simulate some artificial measurements by
% *** adding random errors to the curve-fit equation.  

example_number = 3                       % which example to run: 1 2 or 3  

consts = { example_number };             % a set of constants passed to lm_func

Npnt = 100;                              % number of data points

fn = sprintf('lm_examp_%d',example_number);

t = [1:Npnt]';                           % independent variable, column vector
% true value of parameters ...             column vectors
switch example_number
  case 1 
    a_true  = [ 20 ; -24 ; 30 ; -40 ]; 
  case 2 
    a_true  = [ 20 ;  10 ;  1 ;  50 ]; 
  case 3 
    a_true  = [  6 ;  20 ;  1 ;   5 ]; 
end

Npar = length(a_true);

msmnt_err = 0.5;
y_dat = lm_func(t,a_true,consts);
y_dat = y_dat + msmnt_err*randn(Npnt,1);  % add random measurement errors


% range of values for gridded paramter search
a1 = 0.1*a_true(1):0.1*a_true(1):2*a_true(1);
a2 = 0.1*a_true(2):0.1*a_true(2):2*a_true(2);
a3 = 0.1*a_true(3):0.1*a_true(3):2*a_true(3);
a4 = 0.1*a_true(4):0.1*a_true(4):2*a_true(4);

% gridded parameter search
for ia2 = 1:length(a2);
   for ia4 = 1:length(a4);
	at = [ a_true(1)  a2(ia2) a_true(3) a4(ia4) ];
	delta_y = ( y_dat - lm_func(t,at,consts) );
	X2(ia2,ia4) = (delta_y' * delta_y)/(Npnt-Npar);  % reduced chi-sq.
   end
end

figure(1); % ------------ plot shape of Chi-squared objective function
 clf
 mesh(a2,a4,log10(X2'), 'LineWidth', 2)
  xlabel('a_2')
  ylabel('a_4')
  zlabel('log_{10}(\chi^2_\nu)','rotation',90)
% print(sprintf('%sD.pdf',fn),'-dpdfcrop');

% *** Replace the lines above with a read-in of some
% *** experimentally-measured data.

% initial guess parameters  ...
switch example_number 
  case 1
     a_init  = [  4 ; -5 ; 6 ;    10 ];  
  case 2
     a_init  = [  5 ;  2 ;  0.2 ;  10 ]; 
  case 3
      a_init  = [ 10 ; 50 ;  5  ;  5.7 ]; 
end

% the proper value of the weight is 1/(squared measurement error)
 weight = 1/msmnt_err^2;

a_lb = -10*abs(a_init);
a_ub =  10*abs(a_init);

%{
figure(3)
 clf
 plot(t,y_dat,'o');
  xlabel('t')
  ylabel('y(t)')
%}

% Algorithmic Parameters
%         prnt MaxIter  eps1  eps2  eps3  eps4  lam0  lamUP lamDN UpdateType 
   opts = [  3,    100, 1e-3, 1e-3, 1e-1, 1e-1, 1e-2,    11,    9,        1 ];

[a_fit,Chi_sq,sigma_a,sigma_y,corr,R_sq,cvg_hst] = ...
	lm('lm_func',a_init,t,y_dat,weight,-0.01,a_lb,a_ub,consts,opts);

y_fit = lm_func(t,a_fit,consts);

disp('    initial    true       fit        sigma_a percent')
disp(' -------------------------------------------------------')
disp ([ a_init  a_true  a_fit sigma_a 100*abs(sigma_a./a_fit) ])

n = length(a_fit);

Chi_sq
R_sq
corr

lm_plots ( t, y_dat, y_fit, sigma_y, cvg_hst, fn );

% lm_examp.m  ---------------------------------------------------------------- 
