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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Problem 3:  Investigate sensitivity of x parameters by
% making a contour plot of residuals of the least squares 
% fit as a function of the two estimated rate constants
% alpha(1) and alpha(2), where 
%    y(t) = x(1) exp(alpha(1) t) + x(2) exp(alpha(2) t)
% Let the observation interval vary from 
% t_final = 1 to 6 sec, with 100 observations per second. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize parameters for the problem.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

noiselevel = 0;

xtru = [.5;.5];
alphatru = [-.3;-.7];

beta = -.1:-.05:-.8;
m = length(beta);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make contour plots of the norm of the
% residual as a function of various estimates of
% alpha, using the optimal values of x_1 and x_2.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Vary the length of time for the observations.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for terminate = 1:6,
   npts = terminate*100 + 1;
   t = linspace(0,terminate,npts)';
   ytru = [exp(alphatru(1)*t),exp(alphatru(2)*t)]*xtru;
   y = ytru + noiselevel*randn(size(ytru));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute solutions for various values of alpha.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   maxcondA = 1;
   for i=1:m
      [x,r,A] = expfit(beta(i),t,y);
      resid(i,i) = log10(norm(r));
      for j=i+1:m
         alpha = [beta(i);beta(j)];
         [x,r,A,condA] = expfit(alpha,t,y);
         logresid(i,j) = log10(norm(r));
         logresid(j,i) = logresid(i,j);
         maxcondA = max(maxcondA,condA);
      end
   end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make a contour plot of the log of the residual,
% and display the worst condition number encountered.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   subplot(3,2,terminate)
   [c,h]=contour(beta,beta,logresid,-10:4:-2);
%  clabel(c,h)
   xlabel('\alpha_1')
   ylabel('\alpha_2')
   text(-.4,-.2,(sprintf('t_{final} = %d',terminate)))
   disp(sprintf('For tfinal= %d, maximum condition number = %e', ...
           terminate,maxcondA))

end