%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Solution to CSE Your Homework Assignment Project 7    %
%  Fitting Exponentials: An Interest in Rates            %
%  problem4.m Dianne P. O'Leary   02/04                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Problem 4:  Solve some nonlinear least squares problems
%             and investigate the sensitivity to starting
%             guesses and number of parameters.
% The problems have data 
%    y(t) = x(1) exp(alpha(1) t) + x(2) exp(alpha(2) t)
%            + eta z(t)
% for t = 0:.01:6, 
% where z(t) is normally distributed with variance 1.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global y t

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize parameters for the problems.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

n = 2;           % number of exponential terms
terminate = 6;   % final value of t
npts = terminate*100 + 1;
t = linspace(0,terminate,npts)';
xtru = [.5;.5];
eta = 1.e-4      % noise level

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Loop over the alpha values.          
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

run.init{1} = [3,4,-1,-2]';
run.init{2} = [3,4,-5,-6]';

tru.alpha{1} = [-.3;-.4];
tru.alpha{2} = [-.3;-.31];

nparam(1) = 4;
nparam(2) = 2;

for kkk = 1:2,
   figure
   subplotno = 0;
   alphatru = tru.alpha{kkk};
   disp(sprintf('True alpha = %f, %f',alphatru))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the data ytru and the noisy data y.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   ytru = [exp(alphatru(1)*t),exp(alphatru(2)*t)]*xtru;
   y = ytru + eta*randn(size(ytru));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve the problem with 2 (expeval2) and 4 (expeval) 
% parameters and two different starting guesses (kk = 1:2) 
% using Matlab's lsqnonlin.
% Plot the residuals and compare with residual using
% true parameters.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   for kk = 1:2,
     zz0 = run.init{kk};
     disp('')
     disp(sprintf('Initial alpha = %f, %f',zz0(n+1:2*n)))
     zz(:,1) = lsqnonlin(@expeval,zz0);
     zz(n+1:2*n,2) = lsqnonlin(@expeval2,zz0(n+1:2*n));
     subplotno = subplotno + 1;
     f(:,1) = expeval (zz(:,1));
     subplot(3,2,subplotno)
     plot(t,f(:,1))
     subplotno = subplotno + 1;
     f(:,2) = expeval2(zz(n+1:2*n,2));
     subplot(3,2,subplotno)
     plot(t,f(:,2))

     for i=1:2,
       disp(sprintf('With %d parameters, computed alpha = %f %f', ...
           nparam(i),zz(n+1:2*n,i)'))
       disp(sprintf('  error in alpha = %e, residual norm = %e', ...
           norm(alphatru-sort(zz(n+1:2*n,i)) ), norm(f(:,i)) ))
     end
   end
   zz = [xtru;alphatru];
   f = expeval(zz);
   disp(sprintf('Norm of residual with true alpha = %e',norm(f)))
   subplotno = subplotno + 1;
   subplot(3,2,subplotno)
   plot(t,f)
   subplot(3,2,6)
   axis([0 5 0 5])
   axis off
   text(.1,4,'Top row: 1st guess')
   text(.1,3.2,'Middle row: 2nd guess')
   text(.1,2.4,'Left: 4 parameters')
   text(.1,1.6,'Right: 2 parameters')
   text(.1,.8,'Bottom row: True parameters')

end