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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Problem 5:  Solve a nonlinear least squares problems
%             with noisy data.
% The model is
%    y(t) = x(1) exp(alpha(1) t) + x(2) exp(alpha(2) t)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Modified 12-2008 to change "run" to "myrun" and
% eliminate unused variables.

global y t

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Read the data and initialize parameters for the problems.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fid = fopen('reaction','r');
A = fscanf(fid,'%g');
fclose(fid);
npts = length(A)/2;
A = reshape(A,2,npts)';
t = A(:,1);
y = A(:,2);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve the problem with 2 parameters and 5 
% different starting guesses using 
% Matlab's lsqnonlin.
% Plot the residuals and compare with residual using
% true parameters.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
myrun.init{1} = [-1,-2]';
myrun.init{2} = [-5,-6]';
myrun.init{3} = [-2,-6]';
myrun.init{4} = [ 0,-6]';
myrun.init{5} = [-1,-3]';

figure(1)
subplotno = 0;

for kk = 1:5,
   zz0 = myrun.init{kk};
   disp('')
   disp(sprintf('Initial alpha = %f, %f',zz0'))
   zz(:,kk) = lsqnonlin(@expeval2,zz0);
   f(:,kk) = expeval2(zz(:,kk));
   subplotno = subplotno + 1;
   subplot(3,2,subplotno)
   plot(t,f(:,kk))
   disp(sprintf('  Computed alpha = %f %f', zz(:,kk)'))
   disp(sprintf('  residual norm = %e', norm(f(:,kk)) ))
end

subplotno = subplotno + 1;
subplot(3,2,subplotno)
plot(t,y)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Determine the sensitivity of the solution by looking
% for nearby parameter values that produce almost the
% same residual norm. 
% 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

zcent = zz(:,1);
increment = .04;
step = .001;

beta1 = zcent(1)-increment:step:zcent(1)+increment;
beta2 = zcent(2)-increment:step:zcent(2)+increment;
m = length(beta1);

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

logresid = zeros(m,m);
   for i=1:m
      for j=1:m
         alpha = [beta1(i);beta2(j)];
         [x,r] = expfit(alpha,t,y);
         logresid(i,j) = log10(norm(r));
      end
   end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make a contour plot of the log of the residual.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   figure(2)
   [c,h] = contour(beta1,beta2,logresid,[-2.6:.05:-2.5,-2.36]);
   clabel(c,h)
   xlabel('\alpha_1')
   ylabel('\alpha_2')