%DATDEMO Nonlinear data fitting demonstration.

%   Copyright 1990-2000 The MathWorks, Inc.
%   $Revision: 1.21 $  $Date: 2000/06/16 22:23:50 $

clf
fig=colordef(gcf, 'white');

Data = ...
   [ 0.0000    5.8955
   0.1000    3.5639
   0.2000    2.5173
   0.3000    1.9790
   0.4000    1.8990
   0.5000    1.3938
   0.6000    1.1359
   0.7000    1.0096
   0.8000    1.0343
   0.9000    0.8435
   1.0000    0.6856
   1.1000    0.6100
   1.2000    0.5392
   1.3000    0.3946
   1.4000    0.3903
   1.5000    0.5474
   1.6000    0.3459
   1.7000    0.1370
   1.8000    0.2211
   1.9000    0.1704
   2.0000    0.2636];

echo on;  clc
%   This example demonstrates fitting a nonlinear function to a
%   set of data.  It is used to demonstrate many of the different
%   methods which can be used in the OPTIMIZATION TOOLBOX.

pause   % Strike any key to continue.

%       Consider the following data:

Data
pause   % Strike any key to continue.

%   Let's plot this data.

pause   % Strike any key for plot.

hold off
t = Data(:,1);
y = Data(:,2);
axis([0 2 -0.5 6])
figure(fig)
hold on
plot(t,y,'ro','EraseMode','none'), 
title('Input data')

pause   % Strike any key to continue.
clc
%   We would like to fit the function
%
%     y =  c(1)*exp(-lam(1)*t) + c(2)*exp(-lam(2)*t)
%
%   to the data.  This function has 2 linear parameters c
%   and 2 nonlinear parameters lam.

%       Since the function has a combination of linear and nonlinear
%       parameters, we will separate the solving into two steps.
%       We will use one of the optimization routines such as LSQNONLIN to
%       solve for the nonlinear parameters, and inside our function we will
%       use "\" to solve for the linear parameters.
%
%       We write a function, called FITFUN2, that, given the nonlinear  
%       parameters lam and the data, solves for the current estimate 
%       of the linear parameters and then returns the error in the fit.  
%
%       Next, we provide a guess for initial estimates of the nonlinear 
%       parameters and invoke an optimization routine.

pause % Strike any key to see the M-file function FITFUN2:

type fitfun2

pause   % Strike any key to continue.

echo off
disp('Please Wait - Compiling Optimization Routines')

% test_long is a variable used for auto testing of this routine
if ~exist('test_long')  test_long = 0; end
if exist('method')~=1 method = 8; end
if ~length(method) method = 8; end
l = 2; 

while 1
   lam = [1 0]'; 
   if ~test_long
      clc
      disp(' ')
      disp('   Choose any of the following methods to perform the data fit')
      disp('')        
      disp('        UNCONSTRAINED:    1) Broyden-Fletcher-Golfarb-Shanno')
      disp('                          2) Davidon-Fletcher-Powell')
      disp('                          3) Steepest Descent')
      disp('                          4) Simplex Search')
      disp('         LEAST SQUARES:   5) Gauss-Newton  ')
      disp('                          6) Levenberg-Marquardt ')
      disp('         MINIMAX          7) Seq. Quadratic Progr.')
      disp('')
      disp('                          0) Quit')
      disp('')
      disp('Note: Options 1:6 perform a least squares fit')
      disp('      Option 7 minimizes the worst case error')
      disp('      Gauss-Newton is the fastest method')
   end
   
   
   if test_long 
      if l>=2 
         method=method-1;
         l = 0;
      end
   else 
      method=-1; 
   end
   while (method<0 | method>7)
      method = [];
      while ~length(method)
         method = input('Select a method number: ');
      end
   end
   if (method == 0) 
      disp('End of demo')
      break
   end
    OPTIONS=optimset('LargeScale','off','Display','off');
   if method==2, OPTIONS=optimset(OPTIONS,'HessUpdate','dfp');
   elseif method==3, OPTIONS=optimset(OPTIONS,'HessUpdate','steepdesc');
   elseif method==5, OPTIONS=optimset(OPTIONS,'LevenbergMarquardt','off');
   end
   if test_long
      l = l + 1;  
   else 
      l = [];
   end
   if method~=4&method~=7
      disp('')
      disp('    Choose any of the following line search methods')
      disp('')
      disp('           1) Mixed Polynomial Interpolation')
      disp('           2) Cubic Interpolation')
      disp('')
      
      while ~length(l)
         l = input('Select a line search number: ');
      end
      if l==2, OPTIONS=optimset(OPTIONS,'LineSearch','cubicpoly');  end
   end
   
   disp('')
   OPTIONS = optimset(OPTIONS,'TolX',1e-3);
   t0=clock;
   hold off
   h = plot(t,y,'b-','EraseMode','xor');
   axis([0 2 -0.5 6])
   hold on
   plot(t,y,'ro','EraseMode','none'), 
   title('Input data')
   disp(' ')
   if method==5|method==6
      disp('[lam,resnorm,residual,exitflag,output]= lsqnonlin(''fitfun2'',lam,[],[],OPTIONS,Data,h);');
      [lam,resnorm,residual,exitflag,output]= lsqnonlin('fitfun2',lam,[],[],OPTIONS,Data,h);
      fval = resnorm;
   elseif method==4
      disp('f = inline(''norm(fitfun2(x,Data,h))'',''x'',''Data'',''h'')');
      disp('[lam,fval,exitflag,output] = fminsearch(f,lam,OPTIONS,Data,h)');
      f = inline('norm(fitfun2(x,Data,h))','x','Data','h');
      [lam,fval,exitflag,output] = fminsearch(f,lam,OPTIONS,Data,h);
      fval = fval.* fval;
   elseif method~=7 % 1,2,3
      disp('f = inline(''norm(fitfun2(x,Data,h))'',''x'',''Data'',''h'')');
      disp('[lam,fval,exitflag,output] = fminunc(f,lam,OPTIONS,Data,h)');
      f = inline('norm(fitfun2(x,Data,h))','x','Data','h');
      [lam,fval,exitflag,output] = fminunc(f,lam,OPTIONS,Data,h);
      fval = fval.*fval;
   else % method = 7
      disp('OPTIONS = optimset(OPTIONS,''MinAbsMax'',length(t));');
      disp('f = inline(''fitfun2(x,Data,h)'',''x'',''Data'',''h'');');
      disp('[lam,allfvals,maxfval,exitflag,output]=fminimax(f,lam,[],[],[],[],[],[],[],OPTIONS,Data,h);')
      OPTIONS = optimset(OPTIONS,'MinAbsMax',length(t));
      f = inline('fitfun2(x,Data,h)','x','Data','h');
      [lam,allfvals,maxfval,exitflag,output]=fminimax(f,lam,[],[],[],[],[],[],[],OPTIONS,Data,h);
      fval = allfvals'*allfvals;   
   end
   if test_long
      if method==7
         if maxfval-0.03*(method==7)-(method==3) > 0.15, 
            error(sprintf('Optimization Toolbox in datdemo: maxfval is %g, method is %g',maxfval,method)), 
         end
         
      else
         
         if fval-0.03*(method==7)-(method==3) > 0.15, 
            error(sprintf('Optimization Toolbox in datdemo: fval is %g, method is %g',fval,method)), 
         end
      end
   
   end
   funEvals=sprintf('\n Number of iterations: %g.  Number of function evaluations: %g.', output.iterations, output.funcCount);
   disp(funEvals)
   disp(sprintf(' Norm of the residual at solution ( norm(FITFUN2(x,Data,h))^2 ): %g ',fval));
   execution_time=etime(clock, t0);
   disp(sprintf(' Execution time: %g \n',execution_time));
%   disp(sprintf(' Nonlinear lam: %g %g, linear c:%g %g \n',lam(1),lam(2),c(1),c(2)));

   disp('Strike any key for menu')
   pause
end
close(fig)

