function [x,FVAL,EXITFLAG,OUTPUT,GRAD,HESSIAN] = fminunc(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,computeLambda,f,GRAD,HESS,Hstr,varargin{:});
   OUTPUT.algorithm = large; % override sfminbx output: not using the reflective part of the method   
end


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

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;

