function [x,Resnorm,FVAL,EXITFLAG,OUTPUT,LAMBDA,JACOB] = lsqcurvefit(FUN,x,XDATA,YDATA,LB,UB,options,varargin)
%LSQCURVEFIT Solves non-linear least squares problems.
%   LSQCURVEFIT solves problems of the form:
%   min  sum {(FUN(X,XDATA)-YDATA).^2}  where X, XDATA, YDATA and the values 
%    X                                  returned by FUN can be vectors or 
%                                       matrices.
%
%   X=LSQCURVEFIT(FUN,X0,XDATA,YDATA) starts at X0 and finds
%   coefficients X to best fit the nonlinear function FUN(X,XDATA) 
%   to the data YDATA (in the least-squares sense). FUN is an M-file 
%   that computes a function of X and XDATA and returns a matrix of 
%   the objective function values: F=FUN(X,XDATA). 
%   NOTE: YDATA must be the same size as the matrix F returned by FUN. 
%
%   X=LSQCURVEFIT(FUN,X0,XDATA,YDATA,LB,UB) defines a set of lower and upper
%   bounds on the design variables, X, so that the solution is in
%   the range LB <= X <= UB.  Use empty matrices for LB and UB
%   if no bounds exist. Set LB(i) = -Inf if X(i) is unbounded below; 
%   set UB(i) = Inf if X(i) is unbounded above.
%
%   X=LSQCURVEFIT(FUN,X0,XDATA,YDATA,LB,UB,OPTIONS) minimizes with the default
%   parameters replaced by values in the structure OPTIONS, an argument
%   created with the OPTIMSET function.  See OPTIMSET for details.  Used
%   options are Display, TolX, TolFun, DerivativeCheck, Diagnostics, Jacobian,
%   JacobPattern, LineSearchType, LevenbergMarquardt, MaxFunEvals, MaxIter, 
%   DiffMinChange and DiffMaxChange, LargeScale, MaxPCGIter, PrecondBandWidth,
%   TolPCG, TypicalX. Use the Jacobian option to specify that FUN may be 
%   called with two output arguments where the second, J, is the Jacobian 
%   matrix: [F,J] = feval(FUN,X,XDATA). If FUN returns a vector (matrix) of m 
%   components when X has length n, then J is an m-by-n matrix where J(i,j) 
%   is the partial derivative of F(i) with respect to x(j). (Note that the 
%   Jacobian J is the transpose of the gradient of F.)
%
%   X=LSQCURVEFIT(FUN,X0,XDATA,YDATA,LB,UB,OPTIONS,P1,P2,..) passes the 
%   problem-dependent parameters P1,P2,... directly to the function FUN: 
%   FUN(X,XDATA,P1,P2,...).  Pass empty matrices for OPTIONS to use the 
%   default values.
%
%   [X,RESNORM]=LSQCURVEFIT(FUN,X0,XDATA,YDATA,...) returns the value of the
%   squared 2-norm of the residual at X: sum {(FUN(X,XDATA)-YDATA).^2}.
%
%   [X,RESNORM,RESIDUAL]=LSQCURVEFIT(FUN,X0,...) returns the value of residual,
%   FUN(X,XDATA)-YDATA, at the solution X. 
%
%   [X,RESNORM,RESIDUAL,EXITFLAG]=LSQCURVEFIT(FUN,X0,XDATA,YDATA,...) returns 
%   a string EXITFLAG that describes the exit condition of LSQCURVEFIT.  
%   If EXITFLAG is:
%      > 0 then LSQCURVEFIT converged to a solution X.
%      0   then the maximum number of function evaluations was reached.
%      < 0 then LSQCURVEFIT did not converge to a solution.
%
%   [X,RESNORM,RESIDUAL,EXITFLAG,OUTPUT]=LSQCURVEFIT(FUN,X0,XDATA,YDATA,...) 
%   returns a structure OUTPUT with the number of iterations taken in 
%   OUTPUT.iterations, the number of function evaluations in OUTPUT.funcCount, 
%   the algorithm used in OUTPUT.algorithm, the number of CG iterations (if used) in 
%   OUTPUT.cgiterations, and the first-order optimality (if used) in 
%   OUTPUT.firstorderopt.
%
%   [X,RESNORM,RESIDUAL,EXITFLAG,OUTPUT,LAMBDA]=LSQCURVEFIT(FUN,X0,XDATA,YDATA,...) 
%   returns the set of Lagrangian multipliers, LAMBDA, at the solution: 
%   LAMBDA.lower for LB and LAMBDA.upper for UB.
%
%   [X,RESNORM,RESIDUAL,EXITFLAG,OUTPUT,LAMBDA,JACOBIAN]=LSQCURVEFIT(FUN,X0,XDATA,YDATA,...)
%   returns the Jacobian of FUN at X.
%

%   Copyright (c) 1990-98 by The MathWorks, Inc.
%   $Revision: 1.18 $  $Date: 1998/10/22 20:11:10 $
%   Mary Ann Branch 8-22-96.

%   The default algorithm when OPTIONS.LargeScale = 'off' is the 
%   Levenberg-Marquardt method with a mixed quadratic and cubic line search procedure.  
%   A Gauss-Newton method is selected by setting OPTIONS.LargeScale='off' and 
%   OPTIONS.LevenbergMarquardt='off'. 
%

% ------------Initialization----------------
defaultopt = optimset('display','final','LargeScale','on', ...
   'TolX',1e-6,'TolFun',1e-6,'DerivativeCheck','off',...
   'Diagnostics','off',...
   'Jacobian','off','MaxFunEvals','100*numberOfVariables',...
   'DiffMaxChange',1e-1,'DiffMinChange',1e-8,...
   'PrecondBandWidth',0,'TypicalX','ones(numberOfVariables,1)','MaxPCGIter','max(1,floor(numberOfVariables/2))', ...
   'TolPCG',0.1,'MaxIter',400,'JacobPattern',[], ...
   'LineSearchType','quadcubic','LevenbergMarq','on'); 
% If just 'defaults' passed in, return the default options in X
if nargin==1 & nargout <= 1 & isequal(FUN,'defaults')
   x = defaultopt;
   return
end
xstart=x(:);
numberOfVariables=length(xstart);

large = 'large-scale';
medium = 'medium-scale';

if nargin < 7, options= [];
   if nargin < 6, UB = [];
      if nargin < 5, LB = [];
         if nargin < 4, 
            error('curvefit requires four input arguments');
         end, end, end, end
if nargout > 5
   computeLambda = 1;
else 
   computeLambda = 0;
end

options = optimset(defaultopt,options);
switch optimget(options,'display')
case {'off','none'}
   verbosity = 0;
case 'iter'
   verbosity = 2;
case 'final'
   verbosity = 1;
case 'testing'
   verbosity = Inf;
otherwise
   verbosity = 1;
end

[xstart,l,u,msg] = checkbounds(xstart,LB,UB,numberOfVariables);
if ~isempty(msg)
   EXITFLAG = -1;
   [Resnorm,FVAL,OUTPUT,LAMBDA,JACOB] = deal([]);
   x(:)=xstart;
   if verbosity > 0
      disp(msg)
   end
   return
end
lFinite = ~isinf(l);
uFinite = ~isinf(u);

if min(min(u-xstart),min(xstart-l)) < 0
   xstart = startx(u,l); 
end

% Need extra argname for XDATA
curvefitarg = 1;
lenVarIn = length(varargin);

diagnostics = isequal(optimget(options,'diagnostics','off'),'on');
gradflag =  strcmp(optimget(options,'Jacobian'),'on');
line_search = strcmp(optimget(options,'largescale','off'),'off'); % 0 means trust-region, 1 means line-search

% Convert to inline function as needed
if ~isempty(FUN)  % will detect empty string, empty matrix, empty cell array
   [funfcn, msg] = fprefcnchk(FUN,'lsqcurvefit',length(varargin),gradflag);
else
   errmsg = sprintf('%s\n%s', ...
      'FUN must be a function name, valid string expression, or inline object;', ...
      ' or, FUN may be a cell array that contains these type of objects.');
   error(errmsg)
end

JAC = [];
x(:) = xstart;
try
   switch funfcn{1}
   case 'fun'
      fuser = feval(funfcn{3},x,XDATA,varargin{:});
   case 'fungrad'
      [fuser,JAC] = feval(funfcn{3},x,XDATA,varargin{:});
   case 'fun_then_grad'
      fuser = feval(funfcn{3},x,XDATA,varargin{:}); 
      JAC = feval(funfcn{4},x,XDATA,varargin{:});
   otherwise
      error('Undefined calltype in LSQNONLIN');
   end
catch
   if isa(funfcn{3},'inline')
      
      errmsg = error(['Error in evaluating user supplied inline function: ',formula(funfcn{3})])
   else
      errmsg = error(['Error in evaluating user supplied function: ',funfcn{3}])
   end
   error(errmsg)
end

if ~isequal(size(fuser), size(YDATA))
   error('Function value and YDATA sizes are incommensurate.')
end
fuser = fuser - YDATA;  % preserve fuser shape until after subtracting YDATA 
f = fuser(:);
nfun=length(f);
if isempty(JAC)
   JAC = zeros(nfun, numberOfVariables);
end

% check size of JAC
[Jrows, Jcols]=size(JAC);
if Jrows~=nfun | Jcols ~=numberOfVariables
   errstr = sprintf('%s\n%s%d%s%d\n',...
      'User-defined gradient is not the correct size:',...
      '    the gradient matrix should be ',nfun,'-by-',numberOfVariables);
   error(errstr);
end

caller = 'lsqcurvefit';

% trustregion and enough equations (as many as variables) 
if ~line_search & nfun >= numberOfVariables 
   OUTPUT.algorithm = large;
   
   % trust region and not enough equations -- switch to line_search
elseif ~line_search & nfun < numberOfVariables 
   warnstr = sprintf('%s\n%s\n', ...
      'Large-scale method requires at least as many equations as variables; ',...
      '   switching to line-search method instead.  Upper and lower bounds will be ignored.');
   warning(warnstr);
   OUTPUT.algorithm = medium;
   
   % line search and no bounds  
elseif line_search & isempty(l(lFinite)) & isempty(u(uFinite))
   OUTPUT.algorithm = medium;
   
   % line search and  bounds  and enough equations, switch to trust region 
elseif line_search & (~isempty(l(lFinite)) | ~isempty(u(uFinite))) & nfun >= numberOfVariables
   warnstr = sprintf('%s\n%s\n', ...
      'Line-search method does not handle bound constraints; ',...
      '   switching to trust-region method instead.');
   warning(warnstr);
   OUTPUT.algorithm = large;
   
   % can't handle this one:   
elseif line_search & (~isempty(l(lFinite)) | ~isempty(u(uFinite)))  & nfun < numberOfVariables
   errstr = sprintf('%s\n%s\n%s\n', ...
      'Line-search method does not handle bound constraints ',...
      '   and large-scale method requires at least as many equations as variables; ',...
      '   aborting.');
   error(errstr);
end

if diagnostics > 0
   % Do diagnostics on information so far
   constflag = 0; gradconstflag = 0; non_eq=0;non_ineq=0;lin_eq=0;lin_ineq=0;
   confcn{1}=[];c=[];ceq=[];cGRAD=[];ceqGRAD=[];
   hessflag = 0; HESS=[];
   msg = diagnose('lsqcurvefit',OUTPUT,gradflag,hessflag,constflag,gradconstflag,...
      line_search,options,xstart,non_eq,...
      non_ineq,lin_eq,lin_ineq,l,u,funfcn,confcn,f,JAC,HESS,c,ceq,cGRAD,ceqGRAD);
end

% Execute algorithm
if isequal(OUTPUT.algorithm,large)
   if ~gradflag
      Jstr = optimget(options,'JacobPattern',[]);
      if isempty(Jstr)  
         % Put this code separate as it might generate OUT OF MEMORY error
         Jstr = sparse(ones(Jrows,Jcols));
      end
   else
      Jstr = [];
   end
   [x,FVAL,LAMBDA,JACOB,EXITFLAG,OUTPUT]=...
      snls(funfcn,x,l,u,verbosity,options,f,JAC,YDATA,caller,Jstr,computeLambda,XDATA,varargin{:});
else
   [x,FVAL,JACOB,EXITFLAG,OUTPUT] = ...
      nlsq(funfcn,x,verbosity,options,f,JAC,YDATA,caller,XDATA,varargin{:});
   LAMBDA.ub=[]; LAMBDA.lb=[];   
end
Resnorm = FVAL'*FVAL;

% Reset FVAL to shape of the user-function, fuser
FVAL = reshape(FVAL,size(fuser));

%--end of lsqcurvefit--

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [allfcns,msg] = fprefcnchk(funstr,caller,lenVarIn,gradflag)
%PREFCNCHK Pre- and post-process function expression for FUNCHK.
%   [ALLFCNS,MSG] = PREFUNCHK(FUNSTR,CALLER,lenVarIn,GRADFLAG) takes
%   the (nonempty) expression FUNSTR from CALLER with LenVarIn extra arguments,
%   parses it according to what CALLER is, then returns a string or inline
%   object in ALLFCNS.  If an error occurs, this message is put in MSG.
%
%   ALLFCNS is a cell array: 
%    ALLFCNS{1} contains a flag 
%    that says if the objective and gradients are together in one function 
%    (calltype=='fungrad') or in two functions (calltype='fun_then_grad')
%    or there is no gradient (calltype=='fun'), etc.
%    ALLFCNS{2} contains the string CALLER.
%    ALLFCNS{3}  contains the objective function
%    ALLFCNS{4}  contains the gradient function (transpose of Jacobian).
%  
%    NOTE: we assume FUNSTR is nonempty.
% Initialize
msg='';
allfcns = {};
funfcn = [];
gradfcn = [];

if gradflag
   calltype = 'fungrad';
else
   calltype = 'fun';
end

% {fun}
if isa(funstr, 'cell') & length(funstr)==1
   % take the cellarray apart: we know it is nonempty
   if gradflag
      calltype = 'fungrad';0
   end
   [funfcn, msg] = fcnchk(funstr{1},lenVarIn);
   if ~isempty(msg)
      error(msg);
   end
   
   % {fun,[]}      
elseif isa(funstr, 'cell') & length(funstr)==2 & isempty(funstr{2})
   if gradflag
      calltype = 'fungrad';
   end
   [funfcn, msg] = fcnchk(funstr{1},lenVarIn);
   if ~isempty(msg)
      error(msg);
   end  
   
   % {fun, grad}   
elseif isa(funstr, 'cell') & length(funstr)==2 % and ~isempty(funstr{2})
   
   [funfcn, msg] = fcnchk(funstr{1},lenVarIn);
   if ~isempty(msg)
      error(msg);
   end  
   [gradfcn, msg] = fcnchk(funstr{2},lenVarIn);
   if ~isempty(msg)
      error(msg);
   end
   calltype = 'fun_then_grad';
   if ~gradflag
      warnstr = ...
         sprintf('%s\n%s\n%s\n','Jacobian function provided but OPTIONS.Jacobian=''off'';', ...
         '  ignoring Jacobian function and using finite-differencing.', ...
         '  Rerun with OPTIONS.Jacobian=''on'' to use Jacobian function.');
      warning(warnstr);
      calltype = 'fun';
   end   
   
elseif ~isa(funstr, 'cell')  %Not a cell; is a string expression, function name string or inline object
   [funfcn, msg] = fcnchk(funstr,lenVarIn);
   if ~isempty(msg)
      error(msg);
   end   
   if gradflag % gradient and function in one function/M-file
      gradfcn = funfcn; % Do this so graderr will print the correct name
   end  
else
   errmsg = sprintf('%s\n%s', ...
      'FUN must be a function name, valid string expression, or inline object;', ...
      ' or, FUN may be a cell array that contains these type of objects.');
   error(errmsg)
end

allfcns{1} = calltype;
allfcns{2} = caller;
allfcns{3} = funfcn;
allfcns{4} = gradfcn;
allfcns{5}=[];

