function[ xstar, info] = NLSQ_ERR( fcn, fo, fw, x0, uds, options)
%
%  NLSQ_ERR solves nonlinear least squares problems of the form
%
%                norm( ((fobs-fmodel(x)) ./ fweight) = MIN, 
%                     length(x) = n, length(f...) = m
%
%           with a global unconstained Gauss-Newton method with
%           error oriented convergence criterion and 
%           adaptive trust region strategy.
%
%           See Section 4.3.4 of the monograph
%           P. Deuflhard: Newton Methods for Nonlinear Problems.
%           Affine Invariance and Adaptive Algorithms.
%           Series Computational Mathematics 35 
%           Springer (2004)
%            
%
%  The input and output arguments of the M-function 
%
%    function[ xstar, info] = NLSQ_ERR( fcn, fo, fw, x0, uds, options)
%
%  are:
%
%  IN:
%
%    fcn     - name of a M-file, where the model function is defined
%    fo      - column array of length m holding the observed data
%    fw      - column array of length m holding the weights for the
%              mesurement data in fo
%    x0      - column array of length n holding the initial estimate 
%              for the unknown parameter vector
%    uds     - user defined structure, passed to the model function
%              use uds = [] as a place holder if not needed
%    options - structure with input options, see below
%              use options = [] to use all default options
%
%  OUT:
%
%    xstar   - last (accepted) iterate
%    info    - structure with output information, see below
%
%
%  The problem to be solved must be coded as a M-file with arbirary name.
%  The function must be of the form
%
%    function[ varargout] = problem_name( x, flag, uds)
%
%  where 
%
%    x    - the column vector of length n of parameters
%    flag - lower case string indicating the type of information 
%           that the function should return
%    uds  - a user defined structure
% 
%  According to the value of flag, the function must return
%    flag         varargout
%    ''         - [f]
%                 column vector of lentgth m with the values of the
%                 model function fmodel(x)
%    'init'     - [fobs, fwght, x0, options, uds] 
%                 column vector of lentgth m with observations
%                 column vector of lentgth m with weights
%                 column vector of lentgth n with an initial guess for
%                 the model parameters to be estimated
%                 structure with input options for NLSQ solver, may be []
%                 user defined structue
%    'jacobian' - J
%                 Jacobian (m,n)-matrix J(x) = dF/dx
%  
%
%  The fields of structure options used by NLSQ_ERR (missing fields
%  are created and set with default values) are:
%
%    options.xtol      - required accuracy for projected residual level
%           .kmax      - maximum iteration number
%           .lambdamin - minimal allowed damping factor
%           .lambda0   - initial damping factor
%           .numdif    - type of jacobian generation 
%                        0 - call of user function
%                        1 - internal numerical differentiation
%           .printmon  - level of monitor print information 
%                        0 - no print
%                        1 - Basic Info
%                        2 - Convergence monitor additionally
%                        3 - Further messages
%           .printsol  - level of solution print information 
%                        0 - no print
%                        1 - initial and final iterate
%                        2 - all iterates
%           .fid       - file identifier for print information
%           .iniscalx  - type of initial scaling for iterates 
%                        0 - xwght = max( abs(x0), xthrsh)
%                        1 - xwght = xthrsh
%                        2 - xwght = ones(n,1)
%           .rescalx   - type of rescaling for iterates
%                        0 - xwght = max( abs(xi), xthrsh);
%                        1 - xwght = xthrsh
%                        2 - xwght = ones(n,1)
%
%           .xthrsh    - threshold values for weighting in x
%
%
%  Information given
%
%    info.niter       - number of iterations
%        .convergence - convergence state ('false', 'true')
%        .ifail       - return code
%                       0 - NLSQ_ERR converged
%                      <0 - NLSQ_ERR failed to converge
%                      -1 - maximum iteration limit achieved
%                      -2 - maximum trial step limit achieved
%                     -11 - lambda < lambdamin in a posteriory estimate
%                     -12 - lambda < lambdamin in a priory estimate 
%        .nfcn        - number of function evaluation (numdif not counted)
%        .njac        - number of jacobian evaluations (user or numdif)
%        .ndec        - number of matrix decompositions (OR)
%        .nsol        - number of LSQ solves
%        .esterr      - estimated error of final parameter iterate
%        .xiter(:,1:niter+1)  - Gauss-Newton iterates
%        .flevel(1:niter+1)   - residual level functions F'F
%        .dxnorm(1:niter)     - weighted norm of corrections
%        .dxbarnorm(1:niter)  - weighted norm of simplified corrections
%        .theta(1:niter-1)    - contraction factors
%        .thetabar(1:niter  ) - modified contraction factors
%        .lambda(1:niter)     - damping factors
%        .options             - structure with used options
%
%  See also
%    other NLSQ solver: NLEQ_RES (M-file), NLSCON (FORTRAN)
%    NLSQ examples:     Feulgen, Enyzme  
%
% 
% 
%  Contact   - Bodo Erdmann, erdmann@zib.de
%  Version   - 0.95
%  Date      - 2004/04/16
%  Copyright - Zuse Institute Berlin 2003-2004
%
%

%  Set internal parameters
  imax = 5;
  delkcomp = 'off';
  format short e;
   
%  check input options
  [ options, iniinfo] = zgn_inicheck( fcn, x0, uds, options, 'NLSQ_ERR');

%  Set local copies of input parameters passed in structure options

  xtol = options.xtol;
  kmax = options.kmax;
  pmon = options.printmon;
  psol = options.printsol;
  luout = options.fid;
  numdif = options.numdif;
  lambda0 = options.lambda0;
  lambdamin = options.lambdamin;
  xthrsh = options.xthrsh;
  iniscalx = options.iniscalx;
  rescalx = options.rescalx;

% dimension arrays and matrices

  n = length( x0);
  m = length( fo);
  xk = zeros (n, 1);
  xkp1 = zeros (n, 1);
  dxk = zeros (n, 1);
  dxkm1 = zeros (n, 1);
  dxkp1 = zeros (n, 1);
  dxbarkp1 = zeros (n, 1);
  dxbark = zeros (n, 1);
  dxh = zeros (n, 1);
  xwght = zeros(n,1);
  fm = zeros (m, 1);
  Fk = zeros (m, 1);
  Fkp1 = zeros (m, 1);
  bk = zeros(m, 1);
  bkp1 = zeros(m, 1);
  v = zeros(m, 1);
%
  Jk = [];
  Jkm1 = [];
  Qk =[];
  Rk = [];
%
  info = [];
   
  flevel = zeros (1,kmax+1);
  dxnorm = zeros (1, kmax);
  dxbarnorm = zeros (1, kmax);
  theta = zeros (1, kmax-1);
  thetabar = zeros (1, kmax);
  lambda = zeros (1, kmax);

% initialize method

  sqrtn = sqrt(n);
  ifail = 0;
  convergence = 'false';
  esterr = 1;
  lambdak = lambda0;
  delk = 0;

  k = 0;
  xk = x0(:);
  xiter(:,k+1) = xk(:);
  xwght = zgn_iniscalx (x0, xthrsh, iniscalx);
  info.nfcn = 0;
  info.njac = 0;
  info.ndec = 0;
  info.nsol = 0;

if pmon > 0 
    fprintf(luout,' \n');
    fprintf(luout,' Global Gauss Newton Method, Error Oriented Version \n');
    fprintf(luout,'  - full rank assumed    \n');
    fprintf(luout,'  - unscaled linear systems    \n');
    fprintf(luout,'  - initial damping factor lambda0 = %8.2E \n', lambda0);
    if numdif == 0
      fprintf(luout,'  - analytically given Jacobian used \n');
    else
      fprintf(luout,'  - numerical differentiation for âJacobian \n');
    end
    fprintf(luout,' \n');
    fprintf(luout,' Number of observations  m = %4.0f \n', m);
    fprintf(luout,' Number of parameters    n = %4.0f \n', n);
    fprintf(luout,' \n');
  end
  if pmon > 1
    fprintf(luout, ...
     '  Iter    ResLvl    DxNorm   SDxNorm ThetaBar    Theta   lambda \n');
  end


% initial function evaluation

  [fm] = feval (fcn, xk, '', uds);
  info.nfcn = info.nfcn + 1;
   
  Fk = -(fo-fm) ./ fw;
  Fk_norm = norm(Fk);
  Fk_level = Fk_norm^2;
  flevel(k+1) = Fk_level;

% Top of Gauss Newton iteration loop
%==================================

  while k < kmax;
%
    if numdif == 0
      Jk = feval (fcn, xk, 'jacobian', uds);
    else
      Jk = zgn_numdif (fcn, xk, xwght, uds, fm);
    end
    info.njac = info.njac + 1;

    Jk = diag( fw.^(-1)) * Jk;

    [ Qk, Rk] = qr(Jk);
    info.ndec = info.ndec + 1;
%
    bk = Qk'*Fk;
    dxk = -Rk \ bk;
    info.nsol = info.nsol + 1;
    
    dxk_norm = norm(dxk ./ xwght)/sqrtn;
    dxnorm(k+1) = dxk_norm;
    
    if k > 0
      thetak = dxk_norm/dxkm1_norm;
      theta(k) = thetak;
      
%  compute a priory estimate for lambda

      v = Fk + Jkm1*dxbark;
      v = Qk'*v;
      dbark = -Rk \ v;
      info.nsol = info.nsol + 1;
      dxh = dxbark - dxk + dbark;
      omegabark = (norm(dxh ./ xwght)/sqrtn) / ...
                  (lambdakm1*dxkm1_norm*dxbark_norm);
      lambdak = min( 1, (1-delk)/(omegabark*dxk_norm));
      if lambdak < lambdamin 
        ifail = -12;
        break;
      end
      if strcmp( delkcomp, 'on') == 1
        if lambdakm1 == 1 & lambdak == 1
          delk = thetak;
        else
          v = Fkm1 + Jkm1*dxkm1;
          v = Qk'*v;
          v = -Rk \ v;
          info.nsol = info.nsol + 1;
          del1 = norm(v ./ xwght)/sqrtn;
          v = Qk'*Fkm1;
          v = -Rk \ v;
          info.nsol = info.nsol + 1;
          del2 = norm(v ./ xwght)/sqrtn;
          delk = del1/del2;
        end
      end
      
%      Jplusk = pinv(Jk);
%      Jpluskm1 = pinv(Jkm1);
%      v1 = -Jplusk*Fkm1+Jplusk*Jkm1*Jpluskm1*Fkm1;
%      v2 = -Jplusk*Fkm1;
%      delkh = norm(v1 ./ xwght)/norm(v2 ./ xwght)

    end  
      
    i = 0;
    
%  Top of trial loop
    while i < imax    

% Compute new iterate

      xkp1 = xk + lambdak*dxk;
      
      fm = feval (fcn, xkp1, '', uds);
      info.nfcn = info.nfcn + 1;

      Fkp1 = -(fo-fm) ./ fw;
      Fkp1_norm = norm(Fkp1);
      Fkp1_level = Fkp1_norm^2;

      bkp1 = Qk'*Fkp1;
      dxbarkp1 = -Rk \ bkp1;
      info.nsol = info.nsol + 1;

      dxbarkp1_norm = norm(dxbarkp1 ./ xwght)/sqrtn;
     
      thetabark = dxbarkp1_norm/dxk_norm;
      
      flevel(k+2) = Fkp1_level;
      dxbarnorm(k+1) = dxbarkp1_norm;
      thetabar(k+1) = thetabark;

% Convergence monitor print

      if pmon > 1
        if k == 0
          fprintf(luout,' %5.0f %9.2E %9.2E %9.2E %8.1E %17.1E \n',...
           k, Fk_level, dxk_norm, dxbarkp1_norm, thetabark, lambdak);
        else
          fprintf(luout,' %5.0f %9.2E %9.2E %9.2E %8.1E %8.1E %8.1E \n',...
           k, Fk_level, dxk_norm, dxbarkp1_norm, thetabark, thetak, lambdak);
        end
      end

% Convergence check

      if k > 10 & dxbarkp1_norm <= xtol
        convergence = 'true';
        k = k + 1;
        if pmon > 1
          fprintf(luout,' %5.0f %9.2E \n', ...
           k, Fkp1_level);
        end
        xiter(:,k+1) = xkp1(:);
        lambda(k) = lambdak;
        esterr = thetak*dxk_norm / (1-thetak);
        break;
      end

%  Natural Monotonicity test
 
      if thetabark > 1.0 - delk*lambdak
        if pmon > 2
          fprintf(luout,' natural monotonicity test failed \n');
        end
        dxh = dxbarkp1 - (1-lambdak)*dxk;
        hk = (2*norm(dxh ./ xwght)/sqrtn ) / (lambdak^2*dxk_norm);
        lambdak = min( lambdak/2, (1-delk)/hk);
        if lambdak < lambdamin 
          ifail = -11;
          break;
        end
        i = i + 1;
      else
        break
      end
%  end of trial loop
    end
 
    if i == imax
      ifail = -2;
    end
    if ifail < 0
      break
    end
    
    if strcmp (convergence, 'true') == 1
      break;
    end

    k = k + 1;
    xk = xkp1;
    Fkm1 = Fk;
    Fk = Fkp1;
    Fk_level = Fkp1_level;
    Jkm1 = Jk;
    dxkm1 = dxk;
    dxbark = dxbarkp1;
    lambdakm1 = lambdak;
    lambda(k) = lambdak;
    
% Update x-weights
    xwght = zgn_rescalx (xwght, xkp1, xthrsh, rescalx);
    dxkm1_norm = norm(dxkm1 ./ xwght)/sqrtn;
    dxbark_norm = norm(dxbark ./ xwght)/sqrtn;


    xiter(:,k+1) = xk(:);

    if psol > 1
      fprintf(luout,' \n');
      fprintf(luout,' Iterate is \n');
      for i = 1 : n
        fprintf(luout,' %5.0f  -  %20.13E \n', i, xk(i));
      end
    end

%  Check iteration count

    if k == kmax
      ifail = -1;
      break
    end

% Bottom of main Gauss Newton iteration loop
  end

  xstar = xkp1;
  info.convergence = convergence;
  info.niter = k;
  info.xiter = xiter;
  info.flevel = flevel(1:k+1);
  info.dxnorm = dxnorm(1:k);
  info.dxbarnorm = dxbarnorm(1:k);
  info.theta = theta(1:k-1);
  info.thetabar = thetabar(1:k);
  info.lambda = lambda(1:k);
  info.options = options;
  info.esterr = esterr;

  if strcmp (convergence, 'true') == 1
    if pmon > 0 
      fprintf(luout,' \n');
      fprintf(luout,' Gauss-Newton Method successfully completed \n');
    end
  if psol > 0
    fprintf(luout,' \n');
    fprintf(luout,' Solution is \n');
    for i = 1 : n
      fprintf(luout,' %5.0f  -  %20.13E \n', i, xkp1(i));
    end
  end
  else
    if pmon > 0 
      fprintf(luout,' \n');
      fprintf(luout,' Gauss-Newton Method failed to converge \n');
    end 
    if pmon > 1
      if ifail == -1
        fprintf(luout,'  - more than kmax = %5.0f iterations \n', ...
         kmax);
      elseif ifail == -2
        fprintf(luout,'  - more than imax = %5.0f trial steps \n', ...
         imax);
      elseif ifail == -11| ifail ==-12
        fprintf(luout,'  - lambda =%9.2E < lambdamin =%9.2E  \n', ...
         lambdak, lambdamin);
      end  
    end      
  end 
  if pmon > 0 
    fprintf(luout,'  \n');
    fprintf(luout,' Amount of work: \n');
    fprintf(luout,'  nIter =%5.0f \n  nJac =%5.0f,   nfcn =%5.0f,',... ;
            info.niter, info.njac, info.nfcn);
    fprintf(luout,'   nDec =%5.0f,   nSol =%5.0f \n', ... 
            info.ndec, info.nsol);
  end

% end of function NLSQ_ERR
