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 X0 and finds a minimum X of the function
%   FUN. FUN accepts input X and returns a scalar function value F evaluated 
%   at X. X0 can be a scalar, vector or matrix. 
%
%   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, HessMult, HessUpdate, MaxFunEvals, 
%   MaxIter, DiffMinChange and DiffMaxChange, LargeScale, MaxPCGIter, 
%   PrecondBandWidth, TolPCG, TypicalX. Use the GradObj option to specify that 
%   FUN also returns a second output argument G that is the partial 
%   derivatives of the function df/dX, at the point X. Use the Hessian option
%   to specify that FUN also returns a third output argument H that
%   is the 2nd partial derivatives of the function (the Hessian) at the 
%   point 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
%     FUN can be specified using @:
%        X = fminunc(@myfun,2)
%
%   where MYFUN is a MATLAB function such as:
%
%       function F = myfun(x)
%       F = sin(x) + 3;
%
%     To minimize this function with the gradient provided, modify
%     the MYFUN 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);
%
%     FUN can also be an inline object:
%        x = fminunc(inline('sin(x)+3'),2);
%
%   See also OPTIMSET, FMINSEARCH, FMINBND, FMINCON, @, INLINE.

%   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 1990-2000 The MathWorks, Inc. 
%   $Revision: 1.21 $  $Date: 2000/06/16 01:52:56 $
%   Andy Grace 7-9-90.


% ------------Initialization----------------
defaultopt = struct('Display','final','LargeScale','on', ...
   'TolX',1e-6,'TolFun',1e-6,'DerivativeCheck','off',...   
   'Diagnostics','off',...
   'GradObj','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,...
   'Hessian','off','HessMult',[],...
   'HessPattern','sparse(ones(numberOfVariables))',...
   'HessUpdate','bfgs', ...
   'LineSearchType','quadcubic'); 

% 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'; 

switch optimget(options,'Display',defaultopt,'fast')
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',defaultopt,'fast'),'on');
mtxmpy = optimget(options,'HessMult',defaultopt,'fast');
if isequal(mtxmpy,'hmult')
   warnstr = sprintf('%s\n%s\n%s\n', ...
            'Potential function name clash with a Toolbox helper function:',...
            'Use a name besides ''hmult'' for your HessMult function to',...
            'avoid errors or unexpected results.');
   warning(warnstr)
end
gradflag =  strcmp(optimget(options,'GradObj',defaultopt,'fast'),'on');
hessflag = strcmp(optimget(options,'Hessian',defaultopt,'fast'),'on');
line_search = strcmp(optimget(options,'LargeScale',defaultopt,'fast'),'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] = optimfcnchk(FUN,'fminunc',length(varargin),gradflag,hessflag);
else
   errmsg = sprintf('%s\n%s', ...
      'FUN must be a function 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'))
    warnstr = sprintf('%s\n%s\n', ...
        'Medium-scale method is a Quasi-Newton method and does not use analytic Hessian.',...
        'Hessian flag in options will be ignored (user-supplied Hessian will not be used).');
    warning(warnstr)
    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',defaultopt,'fast');
   if ischar(Hstr) 
      if isequal(lower(Hstr),'sparse(ones(numberofvariables))')
      % Put this code separate as it might generate OUT OF MEMORY error
         Hstr = sparse(ones(n));
      else
         error('Option ''HessPattern'' must be a matrix if not the default.')
      end
   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,defaultopt,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,defaultopt,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


