function [x,FVAL,EXITFLAG,OUTPUT,GRAD,HESSIAN] = fminunc1(FUN,x,options,varargin)
//FMINUNC  Finds the minimum of a function of several variables.
//   X=FMINUNC(FUN,X0) starts at the point X0 and finds a minimum X of the
//   function described in FUN. X0 can be a scalar, vector or matrix.
//   The function FUN (usually an M-file or inline object) should return a scalar 
//   function value F evaluated at X when called with feval: F=feval(FUN,X). 
//   See the examples below for more about FUN.
//
//   X=FMINUNC(FUN,X0,OPTIONS)  minimizes with the default optimization
//   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, GradObj,
//   HessPattern, LineSearchType, Hessian, HessUpdate, MaxFunEvals, MaxIter, DiffMinChange 
//   and DiffMaxChange, LargeScale, MaxPCGIter, PrecondBandWidth, TolPCG, TypicalX. 
//   Use the GradObj option to specify that FUN can be called with two output 
//   arguments where the second, G, is the partial derivatives of the
//   function df/dX, at the point X: [F,G] = feval(FUN,X).  Use Hessian
//   to specify that FUN can be called with three output arguments where the 
//   second, G, is the partial derivatives of the function df/dX, and the third H
//   is the 2nd partial derivatives of the function (the Hessian) at the 
//   point X: [F,G,H] = feval(FUN,X).  The Hessian is only used by the large-scale
//   method, not the line-search method. 
//
//   X=FMINUNC(FUN,X0,OPTIONS,P1,P2,...) passes the problem-dependent 
//   parameters P1,P2,... directly to the function FUN, e.g. FUN would be
//   called using feval as in: feval(FUN,X,P1,P2,...).  
//   Pass an empty matrix for OPTIONS to use the default values.
//
//   [X,FVAL]=FMINUNC(FUN,X0,...) returns the value of the objective 
//   function FUN at the solution X.
//
//   [X,FVAL,EXITFLAG]=FMINUNC(FUN,X0,...) returns a string EXITFLAG that 
//   describes the exit condition of FMINUNC.  
//   If EXITFLAG is:
//      > 0 then FMINUNC converged to a solution X.
//      0   then the maximum number of function evaluations was reached.
//      < 0 then FMINUNC did not converge to a solution.
//   
//   [X,FVAL,EXITFLAG,OUTPUT]=FMINUNC(FUN,X0,...) 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,FVAL,EXITFLAG,OUTPUT,GRAD]=FMINUNC(FUN,X0,...) returns the value 
//   of the gradient of FUN at the solution X.
//
//   [X,FVAL,EXITFLAG,OUTPUT,GRAD,HESSIAN]=FMINUNC(FUN,X0,...) returns the 
//   value of the Hessian of the objective function FUN at the solution X.
//
//   Examples 
//   Minimize the one dimensional function f(x) = sin(x) + 3:
//     To use an M-file, i.e. FUN = 'myfun', create a file myfun.m: 
//        function f = myfun(x)
//         f = sin(x)+3;
//     Then call FMINUNC to find a minimum of FUN near 2:
//        x = fminunc('myfun',2)
//     To minimize this function with the gradient provided, modify
//     the m-file myfun.m so the gradient is the second output argument:
//        function [f,g]= myfun(x)
//         f = sin(x) + 3;
//         g = cos(x);
//     and indicate the gradient value is available by creating an options
//     structure with OPTIONS.GradObj set to 'on' (using OPTIMSET):
//        options = optimset('GradObj','on');
//        x = fminunc('myfun',2,options);
//
//     To minimize the function f(x) = sin(x) + 3 using an inline object:
//        f = inline('sin(x)+3');
//        x = fminunc(f,2);
//     To use inline objects for the function and gradient, FUN is a cell array of two 
//     inline objects where the first is the objective and the second is the
//     gradient of the objective:
//        options = optimset('GradObj','on');
//        x = fminunc({ inline('sin(x)+3'), inline('cos(x)') } ,2,options);
//

//   When options.LargeScale=='on', the algorithm is a trust-region method.
//   When options.LargeScale=='off', the algorithm is the BFGS Quasi-Newton 
//   method with a mixed quadratic and cubic line search procedure. 
//

//   Copyright (c) 1990-98 by The MathWorks, Inc.
//   $Revision: 1.17 $  $Date: 1998/10/22 19:28:39 $
//   Andy Grace 7-9-90.


// ------------Initialization----------------
defaultopt = optimset('display','final','LargeScale','on', ...
   'TolX',1e-6,'TolFun',1e-6,'DerivativeCheck','off',...   
   'Diagnostics','off',...
   'GradObj','off','Hessian','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,'HessPattern','sparse(ones(numberofvariables))','HessUpdate','bfgs', ...
   'LineSearchType','quadcubic','Diagnostics','off'); 

// If just 'defaults' passed in, return the default options in X
if nargin==1 && nargout <= 1 && isequal(FUN,'defaults')
   x = defaultopt;
   return
end

if nargin < 2, error('FMINUNC requires two input arguments');end
if nargin < 3, options=[]; end 
XOUT=x(:);
numberofvariables=length(XOUT);
medium = 'medium-scale: Quasi-Newton line search'; 
large = 'large-scale: trust-region Newton'; 

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
diagnostics = isequal(optimget(options,'diagnostics','off'),'on');
gradflag =  strcmp(optimget(options,'GradObj'),'on');
hessflag = strcmp(optimget(options,'Hessian'),'on');
line_search = strcmp(optimget(options,'largescale','off'),'off'); // 0 means trust-region, 1 means line-search
computeLambda = 0;

// Convert to inline function as needed
if ~isempty(FUN)  // will detect empty string, empty matrix, empty cell array
   [funfcn, msg] = fprefcnchk(FUN,'fminunc',length(varargin),gradflag,hessflag);
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

GRAD=zeros(numberofvariables,1);
HESS = [];

switch funfcn{1}
case 'fun'
   f = feval(funfcn{3},x,varargin{:});
case 'fungrad'
   [f,GRAD(:)] = feval(funfcn{3},x,varargin{:});
case 'fungradhess'
   [f,GRAD(:),HESS] = feval(funfcn{3},x,varargin{:});
case 'fun_then_grad'
   f = feval(funfcn{3},x,varargin{:}); 
   GRAD(:) = feval(funfcn{4},x,varargin{:});
case 'fun_then_grad_then_hess'
   f = feval(funfcn{3},x,varargin{:}); 
   GRAD(:) = feval(funfcn{4},x,varargin{:});
   HESS = feval(funfcn{5},x,varargin{:});
otherwise
   error('Undefined calltype in FMINUNC');
end

// Determine algorithm
// If line-search and no hessian,  then call line-search algorithm
if line_search  && ...
      (~isequal(funfcn{1}, 'fun_then_grad_then_hess') && ~isequal(funfcn{1}, 'fungradhess'))
  output.algorithm = medium; 
    
// Line-search and Hessian -- no can do, so do line-search after warning: ignoring hessian.   
elseif line_search && ...
      (isequal(funfcn{1}, 'fun_then_grad_then_hess') || isequal(funfcn{1}, 'fungradhess'))
   warning('Line-search algorithm requested (OPTIONS.LargeScale = 0) so user-provided Hessian will be ignored.')
   if isequal(funfcn{1}, 'fun_then_grad_then_hess')
      funfcn{1} = 'fun_then_grad';
   elseif isequal(funfcn{1}, 'fungradhess')
      funfcn{1} = 'fungrad';
   end
     output.algorithm = medium;
// If not line-search (trust-region) and Hessian, call trust-region   
elseif ~line_search && ...
      (isequal(funfcn{1}, 'fun_then_grad_then_hess') || isequal(funfcn{1}, 'fungradhess'))
   l=[]; u=[]; Hstr=[];
   output.algorithm = large; 
// If not line search (trust-region) and no Hessian but grad, use sparse finite-differencing.
elseif ~line_search && ...
      (isequal(funfcn{1}, 'fun_then_grad') || isequal(funfcn{1}, 'fungrad'))
   n = length(XOUT); 
   Hstr = optimget(options,'HessPattern',[]);
   if isempty(Hstr)
      // Put this code separate as it might generate OUT OF MEMORY error
      Hstr = sparse(ones(n));
   end
   if ischar(Hstr)  // should be a matrix
      Hstr = eval(Hstr);
   end
   l=[]; u=[];
   output.algorithm = large;
   
// Trust region but no grad, no can do; warn and use line-search    
elseif ~line_search
   warnstr = sprintf('//s\n//s\n', ...
      'Gradient must be provided for trust-region method; ',...
      '   using line-search method instead.');
   warning(warnstr);
   output.algorithm = medium;
else
   error('Problem not handled by fminunc.')   
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;
   LB=[]; UB =[];confcn{1}=[];c=[];ceq=[];cGRAD=[];ceqGRAD=[];
   msg = diagnose('fminunc',output,gradflag,hessflag,constflag,gradconstflag,...
      line_search,options,XOUT,non_eq,...
      non_ineq,lin_eq,lin_ineq,LB,UB,funfcn,confcn,f,GRAD,HESS,c,ceq,cGRAD,ceqGRAD);
   
end

// If line-search and no hessian,  then call line-search algorithm
if isequal(output.algorithm, medium)
   [x,FVAL,GRAD,HESSIAN,EXITFLAG,OUTPUT] = fminusub(funfcn,x,verbosity,options,f,GRAD,HESS,varargin{:});
   
elseif isequal(output.algorithm, large)
   [x,FVAL,LAMBDA,EXITFLAG,OUTPUT,GRAD,HESSIAN] = sfminbx(funfcn,x,l,u,verbosity,options,defaultopt,computeLambda,f,GRAD,HESS,Hstr,varargin{:});
   OUTPUT.algorithm = large; // override sfminbx output: not using the reflective part of the method   
end

endfunction

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

function [allfcns,msg] = fprefcnchk(funstr,caller,lenVarIn,gradflag,hessflag)
//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
//    ALLFCNS{5}  contains the hessian function.
//  
//    NOTE: we assume FUNSTR is nonempty.
// Initialize
msg='';
allfcns = {};
funfcn = [];
gradfcn = [];
hessfcn = [];
if gradflag && hessflag 
   calltype = 'fungradhess';
elseif 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
      error('Gradient function expected but not found.')
   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
      error('Gradient function expected but not found.')
   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','Gradient function provided but OPTIONS.GradObj==''off'';', ...
         '  ignoring gradient function and using finite-differencing.', ...
         '  Rerun with OPTIONS.GradObj=''on'' to use gradient function.');
      warning(warnstr);
      calltype = 'fun';
   end

   
// {fun, [], []}   
elseif isa(funstr, 'cell') && length(funstr)==3 ...
      && ~isempty(funstr{1}) && isempty(funstr{2}) && isempty(funstr{3})
   if gradflag
      error('Gradient function expected but not found.')
   end
   if hessflag
      error('Hessian function expected but not found.')
   end
   
   [funfcn, msg] = fcnchk(funstr{1},lenVarIn);
   if ~isempty(msg)
      error(msg);
   end  
   
// {fun, grad, hess}   
elseif isa(funstr, 'cell') && length(funstr)==3 ...
      && ~isempty(funstr{2}) && ~isempty(funstr{3})
   [funfcn, msg] = fcnchk(funstr{1},lenVarIn);
   if ~isempty(msg)
      error(msg);
   end  
   [gradfcn, msg] = fcnchk(funstr{2},lenVarIn);
   if ~isempty(msg)
      error(msg);
   end
   [hessfcn, msg] = fcnchk(funstr{3},lenVarIn);
   if ~isempty(msg)
      error(msg);
   end
   calltype = 'fun_then_grad_then_hess';
   if ~hessflag && ~gradflag
      warnstr = ...
         sprintf('//s\n//s\n//s\n','Hessian and gradient functions provided ', ...
         '  but OPTIONS.Hessian and OPTIONS.GradObj==''off''; ignoring Hessian and gradient functions.', ...
         '  Rerun with OPTIONS.Hessian=''on'' and OPTIONS.GradObj=''on'' to use derivative functions.');
      warning(warnstr)
      calltype = 'fun';
   elseif hessflag && ~gradflag
      warnstr = ...
         sprintf('//s\n//s\n//s\n','Hessian and gradient functions provided ', ...
         '  but OPTIONS.GradObj==''off''; ignoring Hessian and gradient functions.', ...
         '  Rerun with OPTIONS.Hessian=''on'' and OPTIONS.GradObj=''on'' to use derivative functions.');
      warning(warnstr)
      calltype = 'fun';
   elseif ~hessflag && gradflag
      warnstr = ...
         sprintf('//s\n//s\n//s\n','Hessian function provided but OPTIONS.Hessian=''off'';', ...
         '  ignoring Hessian function,', ...
         '  Rerun with OPTIONS.Hessian=''on'' to use Hessian function.');
      warning(warnstr);
      calltype = 'fun_then_grad';
   end
   
// {fun, grad, []}   
elseif isa(funstr, 'cell') && length(funstr)==3 ...
      && ~isempty(funstr{2}) && isempty(funstr{3})
   if hessflag
      error('Hessian function expected but not found.')
   end
   [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','Gradient function provided but OPTIONS.GradObj==''off'';', ...
         '  ignoring gradient function and using finite-differencing.', ...
         '  Rerun with OPTIONS.GradObj=''on'' to use gradient function.');
      warning(warnstr);
      calltype = 'fun';
   end

// {fun, [], hess}   
elseif isa(funstr, 'cell') && length(funstr)==3 ...
      && isempty(funstr{2}) && ~isempty(funstr{3})
   error('Hessian function given without gradient function.')

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} = hessfcn;

endfunction
