%-----------------------------------------------------------------------
% Nonlinear regression of dynamic data governed by the following ODE
%   Diffusion = Reaction
%   (d^2y/dx^2) = f(y)    BC: at x=0, dy/dx=0
%                             at x=1, y=1
% where
%   f(y) = alpha*y/(beta+y)
% Fit noisy y to find the parameters: alpha and beta
%-----------------------------------------------------------------------
% Programming Note: call "fmins" (old) or "fminsearch"
% Instructor: Nam Sun Wang
%-----------------------------------------------------------------------

% Start fresh ----------------------------------------------------------
      clear all

% Extend xDataz and yDataz to functions --------------------------------
      global xDataz yDataz

% Output a program header ----------------------------------------------
% Program header
      disp('--------------------------------------------------')
      disp('Parameter estimation via nonlinear regression in a')
      disp('one-dimensional diffusion-reaction system.        ')
      disp('  diffusion =  reaction                           ')
      disp('  d2y/dx2 = f(y)                                  ')
      disp('where                                             ')
      disp('  f(y) = alpha*y/(beta+y)                         ')
      disp('--------------------------------------------------')
      disp(' ')

% Generate test data, which we obtain experimentally in practice.
% Find the true concentration profile (by integrating the ODE with finite element method)
      alpha=1.;
      beta =1.;
      [x, y] = finitedb(alpha, beta);

% Add random noise to y-------------------------------------------------
% When noise=0.02 or higher, we fail to converge within the default number of
% iterations in FMINS, and FMINS gives wrong estimates of the parameters.
% When noise=0.01 or lower, we converge to O(1) numbers.  However, the
% estimated parameters may sometimes be quite wrong, although they are sometimes
% close.  Although the parameters are wrong, the fitted y is good in all cases.
      noise = 0.01;    % Noise level
      yDataz = y.*(1+(2*rand(size(y))-1)*noise);

% Initial guess of the parameters (alpha, beta) ------------------------
      guess = [0 0];

% Call FMINS|FMINSEARCH to find the set of parameters that minimize sse, which is
% returned in finiteds.m
%     param = fmins('finiteds', guess); ... old
      param = fminsearch('finiteds', guess);
      alphar = param(1);
      betar  = param(2);

% Generate the fitted function -----------------------------------------
      [xfit, yfit] = finitedb(alphar, betar);

% Find sum of squared errer --------------------------------------------
      sse = finiteds(param);

% Print the coefficients -----------------------------------------------
      disp('The originally assumed parameters are: ')
      disp([ 'alpha = ', num2str(alpha) ])
      disp([ 'beta  = ', num2str(beta)  ])
      disp('The calculated parameters from regression are: ')
      disp([ 'alpha = ', num2str(alphar) ])
      disp([ 'beta  = ', num2str(betar)  ])

% Plot and compare results ---------------------------------------------
      plot(x, yDataz, '.', x, y, '-', x, yfit, '-')
      xlabel('Position')
      ylabel('Concentration')
      title('Concentration Profile in Diffusion-Reaction System')
      text(0.1, 0.980, 'Parameters (True)')
      text(0.1, 0.965, [ '   \alpha = ' num2str(alpha) ])
      text(0.1, 0.950, [ '   \beta  = ' num2str(beta)  ])
      text(0.1, 0.935, 'Parameters (Regression)')
      text(0.1, 0.920, [ '   \alpha = ' num2str(alphar)])
      text(0.1, 0.905, [ '   \beta  = ' num2str(betar) ])
      text(0.1, 0.890, [ 'Error: sse = ' num2str(sse)   ])
      legend('Noisy Data', 'Function (True)', 'Function (Regressed)')