function [x,FVAL,LAMBDA,EXITFLAG,OUTPUT,GRAD,HESSIAN] = sfminbx(funfcn,x,l,u,verb, ...
    options,defaultopt,computeLambda,initialf,initialGRAD,initialHESS,Hstr, ...
    detailedExitMsg,computeConstrViolForOutput,makeExitMsg,varargin)
%

%SFMINBX Nonlinear minimization with box constraints.
%
% Locate a local minimizer to
%
%               min { f(x) :  l <= x <= u }.
%
%	where f(x) maps n-vectors to scalars.

%   Copyright 1990-2019 The MathWorks, Inc.

caller = funfcn{2};
% Set algorithm in order to populate output.algorithm
algorithmName = 'trust-region-reflective'; % when caller is fmincon
if strcmpi(caller,'fminunc')
    % When caller is fminunc
    algorithmName = 'trust-region';
end

% Check to see if the function value at the initial point is well-defined.
% If not, then terminate immediately.
if ~isfinite(initialf) || ~isreal(initialf)
    error(message('optimlib:sfminbx:UsrObjUndefAtX0', caller));
end

%   Initialization
xcurr = x(:);  % x has "the" shape; xcurr is a vector
n = length(xcurr);
iter = 0;
numFunEvals = 1;  % feval of user function in fmincon or fminunc

if n == 0
    error(message('optimlib:sfminbx:InvalidN'))
end

if isempty(l)
   l = -inf*ones(n,1);
else
    % Some variables are possibly bounded (i.e. lb ~= [])
    arg2 = (l <= -1e10);
    l(arg2) = -inf;
end
if isempty(u)
   u = inf*ones(n,1);
else
    % Some variables are possibly bounded (i.e. ub ~= [])
   arg = (u >= 1e10);
   u(arg) = inf;
end

if any(u == l)
    error(message('optimlib:sfminbx:InvalidBounds'))
end

% Add ActiveConstrTol and Preconditioner to defaultopt. These are not
% documented options, and are not added to defaultopt at the user-facing
% function level.
defaultopt.ActiveConstrTol = sqrt(eps);
defaultopt.Preconditioner = @hprecon;

% Read in user options
numberOfVariables = n;

pcmtx = optimget(options,'Preconditioner',defaultopt,'fast');

mtxmpy = optimget(options,'HessMult',defaultopt,'fast');
% Check if name clash
functionNameClashCheck('HessMult',mtxmpy,'hessMult_optimInternal', ...
    'optimlib:sfminbx:HessMultNameClash');

% Use internal Hessian-multiply function if user does not provide HessMult function
% or options.Hessian is off
if isempty(mtxmpy) || (~strcmpi(funfcn{1},'fungradhess') && ~strcmpi(funfcn{1},'fun_then_grad_then_hess'))
    mtxmpy = @hessMult_optimInternal;
end

% Handle the output
outputfcn = optimget(options,'OutputFcn',defaultopt,'fast');
if isempty(outputfcn)
    haveoutputfcn = false;
else
    haveoutputfcn = true;
    xOutputfcn = x; % Last x passed to outputfcn; has the input x's shape
    % Parse OutputFcn which is needed to support cell array syntax for OutputFcn.
    outputfcn = createCellArrayOfFunctions(outputfcn,'OutputFcn');
end

% Handle the plot functions
plotfcns = optimget(options,'PlotFcns',defaultopt,'fast');
if isempty(plotfcns)
    haveplotfcn = false;
else
    haveplotfcn = true;
    xOutputfcn = x; % Last x passed to outputfcn; has the input x's shape
    % Parse PlotFcn which is needed to support cell array syntax for PlotFcn.
    plotfcns = createCellArrayOfFunctions(plotfcns,'PlotFcns');
end

pcflags = optimget(options,'PrecondBandWidth',defaultopt,'fast') ;
tol2 = optimget(options,'TolX',defaultopt,'fast') ;
tol1 = optimget(options,'TolFun',defaultopt,'fast') ;
tolFunValue = optimget(options,'TolFunValue',defaultopt,'fast');
maxiter = optimget(options,'MaxIter',defaultopt,'fast') ;
maxfunevals = optimget(options,'MaxFunEvals',defaultopt,'fast') ;
pcgtol = optimget(options,'TolPCG',defaultopt,'fast') ;  % pcgtol = .1;
kmax = optimget(options,'MaxPCGIter', defaultopt,'fast') ;
if ischar(kmax)
    if isequal(lower(kmax),'max(1,floor(numberofvariables/2))')
        kmax = max(1,floor(numberOfVariables/2));
    else
        error(message('optimlib:sfminbx:InvalidMaxPCGIter'))
    end
end

if ischar(maxfunevals)
    if isequal(lower(maxfunevals),'100*numberofvariables')
        maxfunevals = 100*numberOfVariables;
    else
        error(message('optimlib:sfminbx:InvalidMaxFunEvals'))
    end
end

dnewt = [];
done = false;
posdef = 1; npcg = 0;

pcgit = 0; delta = 10;nrmsx = 1; ratio = 0;
oval = inf;  g = zeros(n,1); newgrad = g; Z = [];

if ~isempty(Hstr)  % use sparse finite differencing
    switch funfcn{1}
        case 'fungrad'
            val = initialf; g(:) = initialGRAD;
        case 'fun_then_grad'
            val = initialf; g(:) = initialGRAD;
        otherwise
            error(message('optimlib:sfminbx:UndefinedCalltype'))
    end
    % Determine coloring/grouping for sparse finite-differencing
    p = colamd(Hstr)';
    p = (n+1) - p;
    group = color(Hstr,p);
    % pass in the user shaped x
    H = sfd(x,g,Hstr,group,[],options.DiffMinChange,options.DiffMaxChange,funfcn,varargin{:});
else % user-supplied computation of H or dnewt
    switch funfcn{1}
        case 'fungradhess'
            val = initialf; g(:) = initialGRAD; H = initialHESS;
        case 'fun_then_grad_then_hess'
            val = initialf; g(:) = initialGRAD; H = initialHESS;
        otherwise
            error(message('optimlib:sfminbx:UndefinedCalltype'))
    end
end

pp = size(g,2);

%   Extract the Newton direction?
if pp == 2 
    dnewt = g(1:n,2); 
end

% Initialize the output function.
if haveoutputfcn || haveplotfcn
   [xOutputfcn, optimValues, stop] = callOutputAndPlotFcns(outputfcn,plotfcns,xcurr,xOutputfcn,'init',iter,numFunEvals, ...
        val,[],[],[],pcgit,[],[],[],delta,l,u,varargin{:});
    if stop
        [x,FVAL,LAMBDA,EXITFLAG,OUTPUT,GRAD,HESSIAN] = ...
            cleanUpInterrupt(xOutputfcn,optimValues,npcg,verb,algorithmName,caller,makeExitMsg);
        return;
    end
end

if verb > 2
    fprintf(['\n                                Norm of      First-order \n',...
        ' Iteration        f(x)          step          optimality   CG-iterations\n']);
    formatstrFirstIter = ' %5.0f      %13.6g                  %12.3g                ';
    formatstr = ' %5.0f      %13.6g  %13.6g   %12.3g     %7.0f';
end

% Create the fault tolerance structure
faultTolStruct = createFaultTolStruct(true);
newFaultTolStruct = faultTolStruct;

%   MAIN LOOP: GENERATE FEAS. SEQ.  xcurr(iter) S.T. f(x(iter)) IS DECREASING.
while ~done
    if any(~isfinite(g)) || ~isreal(g)
        error('optimlib:sfminbx:InvalidUserFunction', ...
            getString(message('optimlib:commonMsgs:GradInfNaNComplex',caller)));
    end

    %     Update
    [v,dv] = definev(g(:,1),xcurr,l,u);
    gopt = v.*g(:,1); 
    gnrm = norm(gopt,inf);
    r = abs(min(u-xcurr,xcurr-l)); 
    degen = min(r + abs(g(:,1)));
    % If no upper and lower bounds (e.g., fminunc), then set degen flag to -1.
    if all( (l == -inf) & (u == inf) )
        degen = -1;
    end

    % Display
    if verb > 2
        if iter > 0
            if faultTolStruct.undefObj
                fprintf(getString(message('optimlib:commonMsgs:ObjInfNaNComplex', ...
                    faultTolStruct.undefValue)));
            end
            currOutput = sprintf(formatstr,iter,val,nrmsx,gnrm,pcgit);
        else
            currOutput = sprintf(formatstrFirstIter,iter,val,gnrm);
        end
        disp(currOutput);
    end    

    % OutputFcn call
    if haveoutputfcn || haveplotfcn
       [xOutputfcn, optimValues, stop] = callOutputAndPlotFcns(outputfcn,plotfcns,xcurr,xOutputfcn,'iter',iter,numFunEvals, ...
            val,nrmsx,g,gnrm,pcgit,posdef,ratio,degen,delta,l,u,varargin{:});
        if stop  % Stop per user request.
            [x,FVAL,LAMBDA,EXITFLAG,OUTPUT,GRAD,HESSIAN] = ...
                cleanUpInterrupt(xOutputfcn,optimValues,npcg,verb,algorithmName,caller,makeExitMsg);
            return;
        end
    end
    
    %     TEST FOR CONVERGENCE
    diff = abs(oval-val);
    oval = val;
    if ((gnrm < tol1) && (posdef == 1) )
        done = true;
        
        EXITFLAG = 1;
        if iter == 0
            msgID = 'optimlib:sfminbx:Exit100detailed';
        else
            msgID = 'optimlib:sfminbx:Exit1detailed';
        end
        % Call createExitMsg with createExitMsgExitflag = 100 if x0 is
        % optimal, otherwise createExitMsgExitflag = 1
        msgData = {{replace(msgID,'detailed','basic')}, ...
            {msgID,gnrm,tol1},verb >= 2,detailedExitMsg};
    elseif (nrmsx < .9*delta) && (ratio > .25) && (diff < tolFunValue*(1+abs(oval)))
        done = true;
        EXITFLAG = 3;
        msgData = {{'optimlib:sfminbx:Exit3basic',caller}, ...
            {'optimlib:sfminbx:Exit3detailed',tolFunValue}, ...
            verb > 1,detailedExitMsg};
    elseif (iter > 1) && (nrmsx < tol2)
        done = true;
        EXITFLAG = 2;
        msgData = {{'optimlib:sfminbx:Exit2basic',caller}, ...
            {'optimlib:sfminbx:Exit2detailed',nrmsx,tol2}, ...
            verb > 1,detailedExitMsg};
    elseif numFunEvals > maxfunevals
        done = true;
        EXITFLAG = 0;
        msgData = {{'optimlib:commonMsgs:Exit0basic',caller,maxfunevals}, ...
                    {},verb > 0,false};
    elseif iter > maxiter
        done = true;
        EXITFLAG = 0;
        % Call createExitMsg with createExitMsgExitflag = 10 for MaxIter exceeded
        msgData = {{'optimlib:commonMsgs:Exit10basic',caller,maxiter}, ...
                    {},verb > 0,false};
    end

    % Reset fault tolerance structure
    faultTolStruct = newFaultTolStruct;

    %     Step computation
    if ~done
        if haveoutputfcn % Call output functions (we don't call plot functions with 'interrupt' flag)
            [~, ~, stop] = callOutputAndPlotFcns(outputfcn,{},xcurr,xOutputfcn,'interrupt',iter,numFunEvals, ...
                val,nrmsx,g,gnrm,pcgit,posdef,ratio,degen,delta,l,u,varargin{:});
            if stop  % Stop per user request.
                [x,FVAL,LAMBDA,EXITFLAG,OUTPUT,GRAD,HESSIAN] = ...
                    cleanUpInterrupt(xOutputfcn,optimValues,npcg,verb,algorithmName,caller,makeExitMsg);
                return;
            end
        end
        %       Determine trust region correction
        dd = abs(v); 
        D = sparse(1:n,1:n,full(sqrt(dd)));
        theta = max(.95,1-gnrm);
        oposdef = posdef;
        [sx,snod,qp,posdef,pcgit,Z] = trdog(xcurr, g(:,1),H,D,delta,dv,...
            mtxmpy,pcmtx,pcflags,pcgtol,kmax,theta,l,u,Z,dnewt,'hessprecon',varargin{:});
        if isempty(posdef) 
            posdef = oposdef; 
        end
        nrmsx = norm(snod); 
        npcg  = npcg + pcgit;
        newx  = xcurr + sx;

        % Perturb?
        [~,newx] = perturbTrustRegionReflective(newx,l,u);

        % Make newx conform to user's input x
        x(:) = newx;
        % Evaluate f, g, and H
        if ~isempty(Hstr) % use sparse finite differencing
            switch funfcn{1}
                case 'fungrad'
                    [newval,newgrad(:)] = feval(funfcn{3},x,varargin{:});
                case 'fun_then_grad'
                    newval = feval(funfcn{3},x,varargin{:});
                    newgrad(:) = feval(funfcn{4},x,varargin{:});
                otherwise
                    error(message('optimlib:sfminbx:UndefinedCalltype'))
            end
            newH = sfd(x,newgrad,Hstr,group,[],options.DiffMinChange,options.DiffMaxChange, ...
                       funfcn,varargin{:});

        else % user-supplied computation of H or dnewt
            switch funfcn{1}
                case 'fungradhess'
                    [newval,newgrad(:),newH] = feval(funfcn{3},x,varargin{:});
                case 'fun_then_grad_then_hess'
                    newval = feval(funfcn{3},x,varargin{:});
                    newgrad(:) = feval(funfcn{4},x,varargin{:});
                    newH = feval(funfcn{5},x,varargin{:});
                otherwise
                    error(message('optimlib:sfminbx:UndefinedCalltype'))
            end

        end
        numFunEvals = numFunEvals + 1;

        % Update fault tolerance structure
        faultTolStruct = updateFaultTolStruct(faultTolStruct, ...
            newval, verb > 2);

        % Update trust region radius
        if ~faultTolStruct.currTrialWellDefined
            % Shrink the trust region if the function value at the new step
            % is not defined (i.e. it's inf/NaN/complex)
            delta = min(nrmsx/20,delta/20);
        else
            % Update trust region using change in objective and in
            % quadratic model
            [~,pp] = size(newgrad);
            aug = .5*snod'*((dv.*abs(newgrad(:,1))).*snod);
            ratio = (newval + aug -val)/qp;

            if (ratio >= .75) && (nrmsx >= .9*delta)
                delta = 2*delta;
            elseif ratio <= .25
                delta = min(nrmsx/4,delta/4);
            end
        end

        % Accept the step if the function value is defined at the new step
        % and the function value is lower than at the current step
        if faultTolStruct.currTrialWellDefined && newval < val
            xcurr = newx;
            val = newval;
            g = newgrad;
            H = newH;
            Z = [];

            %          Extract the Newton direction?
            if pp == 2, dnewt = newgrad(1:n,2); end
        end
        iter = iter + 1;
    end % if ~done
end % while

if haveoutputfcn || haveplotfcn
   callOutputAndPlotFcns(outputfcn,plotfcns,xcurr,xOutputfcn,'done',iter,numFunEvals, ...
        val,nrmsx,g,gnrm,pcgit,posdef,ratio,degen,delta,l,u,varargin{:});
    % Do not check value of 'stop' as we are done with the optimization
    % already.
end
HESSIAN = H;
GRAD = g;
FVAL = val;

OUTPUT.iterations = iter;
OUTPUT.funcCount = numFunEvals;
OUTPUT.stepsize = nrmsx;
OUTPUT.cgiterations = npcg;
OUTPUT.firstorderopt = gnrm;
OUTPUT.algorithm = algorithmName;
% Look up exit message, if requested
if makeExitMsg
    OUTPUT.message = createExitMsg(msgData{:});
end
if computeConstrViolForOutput
    OUTPUT.constrviolation = max([0; (l-xcurr); (xcurr-u) ]);
else
    OUTPUT.constrviolation = [];
end

x(:) = xcurr;

if computeLambda    
    LAMBDA = formBoxLambda(g,l,u);
    LAMBDA.ineqlin = [];
    LAMBDA.eqlin = [];
    LAMBDA.ineqnonlin = [];
    LAMBDA.eqnonlin = [];
else
    LAMBDA = [];
end

%--------------------------------------------------------------------------
function [xOutputfcn, optimValues, stop] = callOutputAndPlotFcns(outputfcn,plotfcns,xvec,xOutputfcn,state, ...
    iter,numFunEvals, ...
    val,nrmsx,g,gnrm,pcgit,posdef,ratio,degen,delta,l,u,varargin)
% CALLOUTPUTANDPLOTFCNS assigns values to the struct OptimValues and then calls the
% outputfcn/plotfcns.
%
% The input STATE can have the values 'init','iter','interrupt', or 'done'.
%
% For the 'done' state we do not check the value of 'stop' because the
% optimization is already done.

optimValues.iteration = iter;
optimValues.funccount = numFunEvals;
optimValues.fval = val;
optimValues.stepsize = nrmsx;
optimValues.gradient = g;
optimValues.firstorderopt = gnrm;
optimValues.cgiterations = pcgit;
optimValues.positivedefinite = posdef;
optimValues.ratio = ratio;
optimValues.degenerate = min(degen,1);
optimValues.trustregionradius = delta;
optimValues.procedure = '';
optimValues.constrviolation = max([0; (l-xvec); (xvec-u) ]);

xOutputfcn(:) = xvec;  % Set x to have user expected size
stop = false;
if ~isempty(outputfcn)
    switch state
        case {'iter','init','interrupt'}
            stop = callAllOptimOutputFcns(outputfcn,xOutputfcn,optimValues,state,varargin{:}) || stop;
        case 'done'
            callAllOptimOutputFcns(outputfcn,xOutputfcn,optimValues,state,varargin{:});
    end
end
% Call plot functions
if ~isempty(plotfcns)
    switch state
        case {'iter','init'}
            stop = callAllOptimPlotFcns(plotfcns,xOutputfcn,optimValues,state,varargin{:}) || stop;
        case 'done'
            callAllOptimPlotFcns(plotfcns,xOutputfcn,optimValues,state,varargin{:});
    end
end
%--------------------------------------------------------------------------
function [x,FVAL,LAMBDA,EXITFLAG,OUTPUT,GRAD,HESSIAN] = cleanUpInterrupt(xOutputfcn,optimValues,npcg, ...
                                                        verb,algorithmName,caller,makeExitMsg)
% CLEANUPINTERRUPT updates or sets all the output arguments of SFMINBX when the optimization
% is interrupted.  The HESSIAN and LAMBDA are set to [] as they may be in a
% state that is inconsistent with the other values since we are
% interrupting mid-iteration.

% Call plot function driver to finalize the plot function figure window. If
% no plot functions have been specified or the plot function figure no
% longer exists, this call just returns.
callAllOptimPlotFcns('cleanuponstopsignal');

x = xOutputfcn;
FVAL = optimValues.fval;
EXITFLAG = -1;
OUTPUT.iterations = optimValues.iteration;
OUTPUT.funcCount = optimValues.funccount;
OUTPUT.stepsize = optimValues.stepsize;
OUTPUT.algorithm = algorithmName;
OUTPUT.firstorderopt = optimValues.firstorderopt;
OUTPUT.cgiterations = npcg; % total number of CG iterations
if makeExitMsg
    OUTPUT.message = createExitMsg({'optimlib:commonMsgs:ExitNeg1basic',caller}, ...
                                    {},verb > 0,false);
end
OUTPUT.constrviolation = optimValues.constrviolation;
GRAD = optimValues.gradient;
HESSIAN = []; % May be in an inconsistent state
LAMBDA = []; % May be in an inconsistent state
