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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Problem 2:  Investigate sensitivity of x parameters
%             to errors in the data.
% Generate 100 problems with 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 uniformly distributed in [-1,1].
% Study the behavior of the solution as a function of
%    -- error level eta
%    -- closeness of the rate constants in alpha
%    -- coordinate system
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Modification 12-2008:  changed variable "cd" to "pcd".

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

n = 2;           % number of exponential terms
nsamp = 100;     % number of trials
terminate = 6;   % final value of t
npts = terminate*100 + 1;
t = linspace(0,terminate,npts)';
xtru = [.5;.5];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize parameters for the plots.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pcd.s{1} = 'r*';
pcd.s{2} = 'b+';
pcd.s{3} = 'mo';
plotno = 0;

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

for kkk = 1:2,

   if (kkk == 1)
       alphatru = [-.3;-.4]
   else
       alphatru = [-.3, -.31]
   end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the data ytru, the data matrix A, and its SVD.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   ytru = [exp(alphatru(1)*t),exp(alphatru(2)*t)]*xtru;

   for j=1:length(alphatru)
      A(:,j) = exp(alphatru(j)*t);
   end
   [u,s,v] = svd(A);
   s = diag(s);
   condA = s(1)/s(n);
   
   disp(sprintf('Condition number of A = %e',condA))
   disp(sprintf('Right singular vectors of A:'))
   disp(v)
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For each noise level eta, solve nsamp problems.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   plotno = plotno + 1;
   figure(plotno)
   subplotno = 0;
   for eta = [1.e-6, 1.e-4, 1.e-2]
      subplotno = subplotno + 1;
      x = zeros(2,nsamp);
      for i=1:nsamp,
         y = ytru + eta*2*(rand(size(ytru))-.5);
         beta = u'*y;
         w(:,i) = beta(1:n) ./ s;
         x(:,i) = v*w(:,i);
         delxy(i) = norm(x(:,i)-xtru)/norm(y-ytru);
      end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the results, as an offset from the true values.  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
      wtru = v'*xtru;
      subplot(3,2,subplotno*2-1)
      plot(x(1,:)-xtru(1),x(2,:)-xtru(2),pcd.s{subplotno})
      xlabel('change in x_1')
      ylabel('change in x_2')

      subplot(3,2,subplotno*2)
      plot(w(1,:)-wtru(1),w(2,:)-wtru(2),pcd.s{subplotno})
      xlabel('change in w_1')
      ylabel('change in w_2')
      title(sprintf('\eta = %e',eta))

      disp(sprintf('Maximum rel_delx/rel_dely = %f', ...
           max(delxy)/norm(xtru)*norm(ytru)))
   end
   
end