classdef (Abstract) OptimFunctions
    % Base class for objective and constraint functions object.
    % Nonlinear solvers are the clients of this object. The derived classes
    % should implement specific behavior for MATLAB functions, Simulation-based
    % functions, or external applications.

    %   Copyright 2023 The MathWorks, Inc.

    properties
        % Finite difference flags. Solver can override when options are not used to set some flags
        % e.g., 'chkFunEval' and 'chkComplexObj' (see barrier). lsq/fsolve 
        % will need to override 'isGrad' flag.
        finDiffFlags = struct('hasLBs',false,'hasUBs',false,'fwdFinDiff',true, ...
            'scaleObjConstr',false,'chkFunEval',false,'chkComplexObj',false, ...
            'isGrad',true)
    end
    
    properties(Access=private)
        % Set for finite difference
        ObjectiveFinDiff 
        ConstraintFinDiff
        % Flag to indicate if gradients are needed
        NeedGrads = false
        LB = []
        UB = []
        solver
        algorithm
        % Some basic size information
        sizes
        
        % Scaling for function outputs
        fscale = optimizationFunctionScaling;
        
        % Cache grads when objective/constraint is called (and fungrad is on)
        ObjEval = struct('X', [], 'Grad',[],'Hess',[])
        ConstrEval = struct('X',[], 'IneqGrad',[], 'EqGrad',[])
    end

    properties(SetAccess=private)
        
        % Classic optim-style (including cell-array) syntax for objecive
        ObjectiveFunctionCell
        % Classic optim-style (including cell-array) syntax for constraints
        ConstraintFunctionCell

        % Nonlinear objective function OR function_and_gradient
        ObjectiveFunAndGrad % fun and fungrad
        % Objective function gradient (separate function from ObjectiveFunAndGrad)
        ObjectiveGradient % _grad_
        % Hessian function (trust-region-reflective only)
        Hessian
        % Nonlinear constraint function OR function_and_gradient
        ConstraintFunAndGrad
        % Constraint function gradient (separate function from ConstraintFunAndGrad)
        ConstraintGradient % _grad_
        % Objective and Constraint function evaluated together
        %ObjectiveAndConstraint % one function for objfun and nonlcon        
    end

    properties(SetAccess=private,Hidden)
        % @todo in a future release.

        % These properties can be useful for a more robust handling of
        % function evaluations e.g., Simulink or external applications.
        
        % Call a function one time before optimization
        PreOptimizationFcns = {} 
        % @note: This is not possible to implement just yet but nice to
        % have in the future

        % Call a function one time after optimization
        PostOptimizationFcns = {} 

        % This function will get exactly the same input as the objective or
        % the constraints function. It can return one logical output to
        % indicate if optimization should halt. This function can be used to:
        % 1. Prepare model with new parameters (x) in external model,
        % writing to files, Simulink set_params, SimInputs, etc.
        % Additionally,
        % 2. If objective and constraints are part of the same evaluation,
        % then the fval, constr values can be stored for 'x'. Later, when
        % objective/constraints are called with the same 'x', return those
        % values.
        % 'x' will be a matrix if UseVectorized is true.
        
        % Function(s) to call before evaluating objective or constraint functions
        PreFunctionEvalFcns = {}        
        % Function(s) to call after evaluating objective or constraint functions
        PostFunctionEvalFcns = {}

    end

    properties(SetAccess=private,Hidden)

        FunArgs = struct('AdditionalParameters',{{}}, ... 
            'PreFunctionEvalFcns',{{}}, ... 
            'PostFunctionEvalFcns', {{}}, ...
            'OnErrorCallbackFcns', {{}})        
        % Classic data structure to keep function
        FunFcn = {}
        ConFcn = {}
        % will not need this after fmincon's full transition
        Flags = struct() % Used by fmincon but created here.

        % Options used for function evaluation and finite difference
        options = struct() % will be private after fmincon's full transition

        % Gradient estimation flags
        GradFlags = struct('ObjfunFinDiff', false, ...
            'ObjfunHaveGrad', false, ...
            'ObjfunHaveHessian',false, ...
            'ObjJacobian', false, ...
            'NonlconHaveGrad', false, ...
            'NonlconFinDiff', false)
    end

    properties(Constant,Hidden)
        % Must use fast because not all options related to function
        % evaluation will be part of an optimoptions for all the solvers.
        optimgetFlag = 'fast'
        % Default options.
        defaultopt = struct('DiffMaxChange',Inf, ...
            'DiffMinChange',0, ...
            'FinDiffRelStep', [], ...
            'FinDiffType','forward', ...
            'FunValCheck','off', ...
            'GradConstr','off', ...
            'GradObj','off', ...
            'DerivativeCheck', 'off', ...
            'TypicalX','ones(numberOfVariables,1)', ...
            'Hessian', [], ...
            'UseParallel',false, ...
            'Vectorized','off', ...
            'ProblemdefOptions', ... 
            optim.internal.constants.DefaultOptionValues.ProblemdefOptions)        
    end

    methods
        function self = OptimFunctions(PV)
            arguments
                %--- cell-style storage ---%
                PV.ObjectiveFunction = []
                PV.ConstraintFunction = []
                
                %--- user-facing new inputs ---%
                PV.Objective {functionHandleOrEmpty(PV.Objective)} = []
                PV.ObjectiveGradient {functionHandleOrEmpty(PV.ObjectiveGradient)} = []
                PV.ObjectiveAndGradient {functionHandleOrEmpty(PV.ObjectiveAndGradient)} = []

                PV.Constraint {functionHandleOrEmpty(PV.Constraint)} = []
                PV.ConstraintGradient {functionHandleOrEmpty(PV.ConstraintGradient)} = []
                PV.ConstraintAndGradient {functionHandleOrEmpty(PV.ConstraintAndGradient)} = []
                
                %PV.ObjectiveAndConstraint {functionHandleOrEmpty(PV.ObjectiveAndConstraint)} = []
                
                PV.PreOptimizationFcns {functionHandleOrEmpty(PV.PreOptimizationFcns)} = {}
                PV.PostOptimizationFcns {functionHandleOrEmpty(PV.PostOptimizationFcns)} = {}

                PV.AdditionalParameters cell = {}
                PV.PreFunctionEvalFcns {functionHandleOrEmpty(PV.PreFunctionEvalFcns)} = {}
                PV.PostFunctionEvalFcns {functionHandleOrEmpty(PV.PostFunctionEvalFcns)} = {}
                PV.OnErrorCallbackFcns {functionHandleOrEmpty(PV.OnErrorCallbackFcns)} = {}

            end

            self.FunArgs.AdditionalParameters = PV.AdditionalParameters;
            self.FunArgs.PreFunctionEvalFcns = PV.PreFunctionEvalFcns;
            self.FunArgs.PostFunctionEvalFcns = PV.PostFunctionEvalFcns;
            self.FunArgs.OnErrorCallbackFcns = PV.OnErrorCallbackFcns;

            self.PreOptimizationFcns = PV.PreOptimizationFcns;
            self.PostOptimizationFcns = PV.PostOptimizationFcns;
            self.ObjectiveFunctionCell = PV.ObjectiveFunction;
            self.ConstraintFunctionCell = PV.ConstraintFunction;

            % Store function handles in
            % ObjectiveFunctionCell/ConstraintFunctionCell only. More error 
            % checking could be done if this class is user facing.
            if isempty(self.ObjectiveFunctionCell)
                if ~isempty(PV.ObjectiveAndGradient)
                    self.ObjectiveFunctionCell = PV.ObjectiveAndGradient;
                elseif ~isempty(PV.Objective) && ~isempty(PV.ObjectiveGradient)
                    self.ObjectiveFunctionCell = {PV.Objective, PV.ObjectiveGradient};
                elseif ~isempty(PV.Objective)
                    self.ObjectiveFunctionCell = PV.Objective;
                end
            end
            if isempty(self.ConstraintFunctionCell)
                if ~isempty(PV.ConstraintAndGradient)
                    self.ConstraintFunctionCell = PV.ConstraintAndGradient;
                elseif ~isempty(PV.Constraint) && ~isempty(PV.ConstraintGradient)
                    self.ConstraintFunctionCell = {PV.Constraint, PV.ConstraintGradient};
                elseif ~isempty(PV.Constraint)
                    self.ConstraintFunctionCell = PV.Constraint;
                end
            end            

        end
        
        function self = setup(self, PV)
            % Set options related to how the functions and gradients are
            % evaluated. Some setup requires algorithm/solver names
            % (optimfcnchk) but it is very rare to pass these arguments. X0
            % and bounds are required for finite difference. Bounds are
            % assumed to be unbounded if not provided.
            arguments
                self 
                PV.NeedGrads logical = false
                PV.X0 double = []
                PV.LB double = []
                PV.UB double = []
                PV.Solver char = ''
                PV.Algorithm char = ''
                PV.Options struct = struct()
            end

           self.NeedGrads = PV.NeedGrads; 
           self = self.setOptions(PV.Solver, PV.Algorithm, PV.Options);
           self = self.setX0AndBounds(PV.X0, PV.LB, PV.UB);

           % Parse function inputs and set specific function evaluators
           % based on options and inputs.
           self = self.optimfcncheck();

           % Setup function handles on the workers if using PCT
           useParallel = optimget(self.options,'UseParallel',self.defaultopt,self.optimgetFlag);
           self.options.UseParallel = validateopts_UseParallel(useParallel,true,true);
           if self.options.UseParallel
               self = self.setupParallel();
           end

           if PV.NeedGrads && ~isempty(PV.X0)
               self = self.setFiniteDifferenceFlags();
           end

        end

        function self = setScale(self, fscale)
            % Provide a scaling object to be used for scaling function and
            % gradients after evaluation. This should be set by the
            % solvers.
            arguments
                self 
                fscale optimizationFunctionScaling
            end
            self.fscale = fscale;
            self.finDiffFlags.scaleObjConstr = ... 
                fscale.objIsScaled || fscale.constrIsScaled;
        end

        function [fval, self] = objective(self, X, fval, fgrad, hess)
            % Get the objective function value at X. If the function
            % returns higher order values (gradient and Hessian), they will
            % be stored so the caller must get the object back as the
            % second argument. The objgrad and objhessian function should be
            % called on the returned object.
            %
            % Input arguments fval, fgrad and hess can be passed as
            % initialized quantities (default is []) in which case the
            % return arguments have the same sizes.
            arguments
                self 
                X double
                fval double = []
                fgrad double = []
                hess double = []
            end

            % Insert PreFunctionEvalFcns calls

            [fval_, fgrad_, hess_] = self.ObjectiveFunAndGrad(self,self.FunFcn{3},...
                X,fgrad,hess);

            % Insert PostFunctionEvalFcns calls

            % Reshape output
            fval = reshapeOutputs(fval_,fval);

            if self.fscale.objIsScaled
                fval = self.fscale.scaleObjective(fval);
            end
            
            if self.GradFlags.ObjfunHaveGrad
                fgrad = reshapeOutputs(fgrad_,fgrad);
                % cache just the grad and Hess
                self.ObjEval.X = X;
                self.ObjEval.Grad = fgrad;
                if self.GradFlags.ObjfunHaveHessian
                    hess = reshapeOutputs(hess_,hess);
                    self.ObjEval.Hess = hess;
                end
            end
        end

        function [fgrad,numEvals,evalOK,self] = objgrad(self, X, fval, fgrad)
            % Get objective gradient. There are three possibilities: 
            % 1. If the gradient was returned when the objective() function
            % was called, it will be simply returned as output.
            % 2. If there is a separate function provided for objective
            % gradient calculation, it will be called and the value will be
            % returned.
            % 3. Lastly, finite difference will be used to estimate the
            % gradients. The Options argument provided in setup() function
            % call can be used to change the default parameters of finite
            % difference. 
            %
            % The output argument 'numEvals' indicates the number of new X
            % evaluated (for finite difference or grad). 
            % evalOK indicates if objective or constraint functions failed
            % during finite difference. If evalOK==false, then there was an
            % error (and solver may stop).
            %
            arguments
                self
                X double
                fval double = []
                fgrad double = []
            end

            if self.GradFlags.ObjfunFinDiff
                % Estimate using finite difference.
                if nargin ~= 4
                    % used in finite difference, fval is assumed to be
                    % non-empty
                    fgrad = self.initObjGrad(X,fval);
                end
                % (X,fval) are required inputs for finite difference but no
                % checking is done here.               
                [fgrad,~,~,numEvals,evalOK] = ...
                    computeFinDiffGradAndJac(X,self.ObjectiveFinDiff,[], ...
                    fval,[],[],fgrad,[],[],self.LB,self.UB,self.fscale, ...
                    self.options,self.finDiffFlags,self.sizes);
                % fgrad is scaled by fin-diff so we can return now.
                return
            end

            evalOK = true; % only relevant for finite difference
            
            if self.GradFlags.ObjfunHaveGrad
                % When objective is called before objgrad, this object uses
                % a cached copy for grad (ONLY when fungrad is true).
                useCache = isequal(X, self.ObjEval.X);
                if ~useCache
                    % Copy the cached values
                    oldObjEval = self.ObjEval;
                    [~, self] = self.objective(X, [], fgrad, []);
                    fgrad_ = self.ObjEval.Grad;
                    % Restore cache
                    self.ObjEval = oldObjEval;
                    numEvals = 1;
                else
                    fgrad_ = self.ObjEval.Grad;
                    numEvals = 0;
                end
                fgrad = reshapeOutputs(fgrad_,fgrad);
            elseif ~isempty(self.ObjectiveGradient)
                % Cache is not relevant
                % fun_then_grad or 'fun_then_grad_then_hess'
                 fgrad_ = self.ObjectiveGradient(self,self.FunFcn{4},X);
                 fgrad = reshapeOutputs(fgrad_,fgrad);
                 numEvals = 0; % assumption is that X is evaluated for fval
            end
            
            if self.fscale.objIsScaled            
                fgrad = self.fscale.scaleObjectiveGrad(fgrad);
            end

        end

        function [hess, self] = objhessian(self, X, hess)
            % Get Hessian of the objective. There are two possibilities: 
            % 1. If the Hessian was returned when the objective() function
            % was called, it will be simply returned.
            % 2. If there is a separate function provided for objective
            % Hessian calculation, it will be called and the value will be
            % returned.
            % If there is no Hessian from (1) or (2), [] will be returned.
            arguments
                self
                X double
                hess double = []
            end
            
            if self.GradFlags.ObjfunHaveHessian && isequal(X,self.ObjEval.X)
                % available from the cache 
                hess = self.ObjEval.Hess;
            elseif ~isempty(self.Hessian)
                hess_ = self.Hessian(self,self.FunFcn{5},X);
                hess = reshapeOutputs(hess_,hess);
            else
                hess = [];
            end
            
        end
        
        function [cin, ceq, self] = constraints(self, X, cin, ceq, cingrad, ceqgrad)
            % Get the constraint function value at X. If the function
            % returns gradient, they will be stored so the caller must get
            % the object back as the second argument. The constraintgrad()
            % function should be called on the returned object.
            %
            % Input arguments cin, ceq, cingrad, and ceqgrad can be passed
            % as initialized quantities (default is []) in which case the
            % return arguments have the same sizes.
            arguments
                self 
                X double 
                cin double = []
                ceq double = []
                cingrad double = []
                ceqgrad double = []
            end
            
            % Insert PreFunctionEvalFcns calls
            [cin_,ceq_,cingrad_,ceqgrad_] = self.ConstraintFunAndGrad(self,self.ConFcn{3},...
                X,cingrad,ceqgrad);
            % Insert PostFunctionEvalFcns calls

            cin = reshapeOutputs(cin_,cin);
            ceq = reshapeOutputs(ceq_,ceq);

            if self.fscale.constrIsScaled
                [cin, ceq] = self.fscale.scaleConstraints(cin,ceq);
            end

            if self.GradFlags.NonlconHaveGrad
                cingrad = reshapeOutputs(cingrad_,cingrad);
                ceqgrad = reshapeOutputs(ceqgrad_,ceqgrad);
                % Cache only the gradients and X
                self.ConstrEval.X = X;
                self.ConstrEval.IneqGrad = cingrad;
                self.ConstrEval.EqGrad = ceqgrad;
            end

        end

        function [JacCineqTrans,JacCeqTrans,numEvals,evalOK,self] = constraintgrad(... 
                self, X, cin, ceq, JacCineqTrans,JacCeqTrans)
            % Get constraint gradients. There are three possibilities: 
            % 1. If the gradient was returned when the constraints() function
            % was called, it will be simply returned as output.
            % 2. If there is a separate function provided for constraint
            % gradient calculation, it will be called and the value will be
            % returned.
            % 3. Lastly, finite difference will be used to estimate the
            % gradients if neither (1) or (2) applies. The Options argument
            % provided in setup() function call can be used to change the
            % default parameters of finite difference. 
            %
            % The output argument 'numEvals' indicates the number of new X
            % evaluated (for finite difference or grad). 
            % evalOK indicates if objective or constraint functions failed
            % during finite difference. If evalOK==false, then there was an
            % error (and solver may stop).
            %
            arguments
                self 
                X double
                cin double = []
                ceq double = []
                JacCineqTrans double = []
                JacCeqTrans double = []
            end
            
            if self.GradFlags.NonlconFinDiff
                if nargin ~= 6
                    [JacCineqTrans,JacCeqTrans] = self.initConstrGrad(X, cin, ceq);
                end
                % finite difference estimates
                [~,JacCineqTrans,JacCeqTrans,numEvals,evalOK] = ...
                    computeFinDiffGradAndJac(X,[],self.ConstraintFinDiff, ...
                    [],cin,ceq,[],JacCineqTrans,JacCeqTrans,self.LB,self.UB,...
                    self.fscale,self.options,self.finDiffFlags,self.sizes);
                % gradients are scaled by fin-diff so return now.
                return
            end
            
            evalOK = true; % only relevant for finite difference                        

            if self.GradFlags.NonlconHaveGrad
                % When constraints is called before constraintgrad, this
                % object uses a cached copy for grad (ONLY when fungrad is
                % true)
                useCache = isequal(X, self.ConstrEval.X);
                if ~useCache
                    % Copy the cached values
                    oldConstrEval = self.ConstrEval;
                    [~, ~, self] = self.constraints(X, [], [], JacCineqTrans, JacCeqTrans);
                    JacCineqTrans_ = self.ConstrEval.IneqGrad;
                    JacCeqTrans_ = self.ConstrEval.EqGrad;
                    % Restore cache
                    self.ConstrEval = oldConstrEval;
                    numEvals = 1;
                else
                    JacCineqTrans_ = self.ConstrEval.IneqGrad;
                    JacCeqTrans_ = self.ConstrEval.EqGrad;
                    numEvals = 0;
                end
                JacCineqTrans = reshapeOutputs(JacCineqTrans_,JacCineqTrans);
                JacCeqTrans = reshapeOutputs(JacCeqTrans_,JacCeqTrans);
            elseif ~isempty(self.ConstraintGradient)
                % Cache is not relevant
                % fun_then_grad or 'fun_then_grad_then_hess'
                [JacCineqTrans_,JacCeqTrans_] = ...
                    self.ConstraintGradient(self,self.ConFcn{4},X);
                JacCineqTrans = reshapeOutputs(JacCineqTrans_,JacCineqTrans);
                JacCeqTrans = reshapeOutputs(JacCeqTrans_,JacCeqTrans);
                numEvals = 0; % assumption is that cin,ceq was evaluated at X
            end

            if self.fscale.constrIsScaled
                [JacCineqTrans, JacCeqTrans] = ... 
                    self.fscale.scaleConstraintGrads(JacCineqTrans,JacCeqTrans);
            end

            
        end
        
        function [grad,JacCineqTrans,JacCeqTrans,numEvals,evalOK,self] = ... 
                computeGradAndJac(self,X,fval,cin,ceq,grad,... 
                JacCineqTrans,JacCeqTrans,findiffOnly)
            % Get objective gradient (grad) and constraint gradients
            % (JacCineqTrans,JacCeqTrans). There are three possibilities:
            %
            % 1. If the gradient was returned when the
            % objective/constraints() functions were called, they will be
            % simply returned as output. 
            % 2. If there is a separate function provided for gradient
            % calculation, they will be called and the value will be returned.
            % 3. Lastly, finite difference will be used to estimate the
            % gradients. The Options argument provided in setup() function
            % call can be used to change the default parameters of finite
            % difference. The output argument 'numEvals' indicates the
            % number of new X evaluated (for finite difference or grad).
            % 'evalOK' is relevant for finite difference, else it is always
            % 'true'.
            %
            % X is required.
            %
            % fval,grad are required if finite difference is used for
            % objective.
            %
            % cin,ceq,JacCineqTrans,JacCeqTrans are required for finite
            % difference is used for constraints.
            %
            % findiffOnly can be set to 'true' in which case, only finite
            % difference will be tried. This is useful when solver may call
            % this function at X0 when gradients are available but not
            % cached in this object (fun_then_grad). In that case, calling
            % this method with default inputs will trigger a call to the
            % gradient function.
            %
            % The output argument 'numEvals' indicates the number of new X
            % evaluated (for finite difference or grad). 
            % evalOK indicates if objective or constraint functions failed
            % during finite difference. If evalOK==false, then there was an
            % error (and solver may stop).
            %
            arguments
                self
                X 
                fval double = []
                cin double = []
                ceq double = []
                grad double = []
                JacCineqTrans double = []
                JacCeqTrans double = []
                findiffOnly logical = false
            end
            
            % Only finite-difference can set this to false later
            evalOK = true;
            numEvals = 0;

            if self.GradFlags.ObjfunFinDiff || self.GradFlags.NonlconFinDiff
                % finite difference estimation
                [grad,JacCineqTrans,JacCeqTrans,numFindiffEvals,evalOK] = ...
                    computeFinDiffGradAndJac(X,self.ObjectiveFinDiff,self.ConstraintFinDiff, ...
                    fval,cin,ceq,grad,JacCineqTrans,JacCeqTrans, ...
                    self.LB,self.UB,self.fscale,self.options,self.finDiffFlags, ...
                    self.sizes);
                
                numEvals = numEvals + numFindiffEvals;
            end

            if findiffOnly               
                return
            end

            if ~self.GradFlags.ObjfunFinDiff && ~isempty(fval)                
                % Get grad from fungrad (cached) or fun_then_grad                
                [grad,numFunEvals,evalOK,self] = ... 
                    self.objgrad(reshape(X,self.sizes.xShape),[],grad);
                numEvals = numEvals + numFunEvals;
            end
                
            if ~self.GradFlags.NonlconFinDiff && ... 
                    (~isempty(cin) || ~isempty(ceq))
                % Get grad from fungrad (cached) or fun_then_grad
                [JacCineqTrans,JacCeqTrans,numFunEvals,evalOK,self] = ... 
                    self.constraintgrad(reshape(X,self.sizes.xShape),[],[],... 
                    JacCineqTrans,JacCeqTrans);
                numEvals = numEvals + numFunEvals;
            end
        end

        function [fval,grad,hess,self] = objective_first_eval(self,X0)
            % Evaluate objective and analytical gradients.Evaluation is
            % done in a try-catch block and error will be thrown.

            try
                % This call will get objective but also store grad and Hess
                % if that is returned by the user function.
                [fval,self] = self.objective(X0);
            catch userFcn_ME
                optim_ME = MException('optimlib:fmincon:ObjectiveError', ...
                    getString(message('optimlib:fmincon:ObjectiveError')));
                userFcn_ME = addCause(userFcn_ME,optim_ME);
                rethrow(userFcn_ME)
            end

            if self.NeedGrads
                if self.GradFlags.ObjfunFinDiff
                    % We just need the correct size for grad
                    grad = self.initObjGrad(X0,fval);
                                        
                else % user supplied gradients
                    % If user grad is in a separate function it will be
                    % called, else we get grad from the last cache.
                    try
                        grad = self.objgrad(X0);
                    catch userFcn_ME
                        optim_ME = MException('optimlib:fmincon:GradientError', ...
                            getString(message('optimlib:fmincon:GradientError')));
                        userFcn_ME = addCause(userFcn_ME,optim_ME);
                        rethrow(userFcn_ME)
                    end
                end

                try
                    % hess will be returned as [] if there is no Hessian
                    hess = self.objhessian(X0);
                catch userFcn_ME
                    optim_ME = MException('optimlib:fmincon:HessianError', ...
                        getString(message('optimlib:fmincon:HessianError')));
                    userFcn_ME = addCause(userFcn_ME,optim_ME);
                    rethrow(userFcn_ME)
                end
            else
                grad = [];
                hess = [];
            end
        end

        function [cin,ceq,cingrad,ceqgrad,self] = constraint_first_eval(self, X0)
            % Evaluate both constraints and analytical gradients (if
            % available). Evaluation is done in a try-catch block.

            try
                % This call will get objective but also store grad and Hess
                % if that is returned by the user function
                [cin,ceq,self] = self.constraints(X0);
            catch userFcn_ME
                if strcmpi('MATLAB:maxlhs',userFcn_ME.identifier)
                    error(message('optimlib:fmincon:InvalidHandleNonlcon'))
                else
                    optim_ME = MException('optimlib:fmincon:NonlconError', ...
                        getString(message('optimlib:fmincon:NonlconError')));
                    userFcn_ME = addCause(userFcn_ME,optim_ME);
                    rethrow(userFcn_ME)
                end
            end

            noConstr = isempty(cin) && isempty(ceq);
            
            if self.NeedGrads
                nVar = numel(X0);
                % special case to avoid finite-difference calls for
                % nonlcon. This makes the code not do any unnecessary work.
                if noConstr
                    self.ConstraintFinDiff = [];
                end

                if self.GradFlags.NonlconFinDiff || noConstr
                    [cingrad,ceqgrad] = self.initConstrGrad(X0, cin, ceq);
                else % user supplied gradients
                    % If user grad is in a separate function it will be
                    % called, else we get grad from the cache.
                    try
                        [cingrad,ceqgrad,~,~,self] = self.constraintgrad(X0);
                    catch userFcn_ME
                        optim_ME = MException('optimlib:fmincon:NonlconFunOrGradError', ...
                            getString(message('optimlib:fmincon:NonlconFunOrGradError')));
                        userFcn_ME = addCause(userFcn_ME,optim_ME);
                        rethrow(userFcn_ME)
                    end
                    if isempty(cingrad)
                        cingrad = reshape(cingrad,nVar,0);
                    end
                    if isempty(ceqgrad)
                        ceqgrad = reshape(ceqgrad,nVar,0);
                    end                    
                end
                
                % Size checking                
                [cgrow,cgcol] = size(cingrad);
                [ceqgrow,ceqgcol] = size(ceqgrad);

                if cgrow ~= nVar || cgcol ~= numel(cin)
                    error(message('optimlib:fmincon:WrongSizeGradNonlinIneq', nVar, numel(cin)))
                end
                if ceqgrow ~= nVar || ceqgcol ~= numel(ceq)
                    error(message('optimlib:fmincon:WrongSizeGradNonlinEq', nVar, numel(ceq)))
                end
            else
                cingrad = [];
                ceqgrad=[];
            end
        end
        
    end

    methods (Access=private)

        function self = setFiniteDifferenceFlags(self)
            
            self.finDiffFlags.hasLBs = isfinite(self.LB);
            self.finDiffFlags.hasUBs = isfinite(self.UB);

            % Gather options needed for finitedifferences
            % Write checked DiffMaxChange, DiffMinChage, FinDiffType, FinDiffRelStep,
            % GradObj and GradConstr options back into struct for later use
            self.options.DiffMinChange = optimget(self.options,'DiffMinChange',self.defaultopt,self.optimgetFlag);
            self.options.DiffMaxChange = optimget(self.options,'DiffMaxChange',self.defaultopt,self.optimgetFlag);
            if self.options.DiffMinChange >= self.options.DiffMaxChange
                error(message('optimlib:fmincon:DiffChangesInconsistent', ...
                    sprintf( '%0.5g', self.options.DiffMinChange ), ...
                    sprintf( '%0.5g', self.options.DiffMaxChange )))
            end
            % Read in and error check option TypicalX
            [typicalx,ME] = getNumericOrStringFieldValue('TypicalX','ones(numberOfVariables,1)', ...
                ones(self.sizes.nVar,1),'a numeric value',self.options,self.defaultopt);
            if ~isempty(ME)
                throw(ME)
            end
            checkoptionsize('TypicalX', size(typicalx), self.sizes.nVar);
            self.options.TypicalX = typicalx;
            
            self.options.FinDiffType = optimget(self.options,'FinDiffType',self.defaultopt,self.optimgetFlag);
            self.options = validateFinDiffRelStep(self.sizes.nVar,self.options,self.defaultopt,self.optimgetFlag);

            % Create default structure of flags for finitedifferences:
            % This structure will (temporarily) ignore some of the features that are
            % algorithm-specific (e.g. scaling and fault-tolerance) and can be turned
            % on later for the main algorithm.
            self.finDiffFlags.fwdFinDiff = strcmpi(self.options.FinDiffType,'forward');
            self.finDiffFlags.chkFunEval = false;     
            self.finDiffFlags.chkComplexObj = false;

            % finDiffFlags.scaleObjConstr could be set via setScale before.
            self.finDiffFlags.scaleObjConstr = self.finDiffFlags.scaleObjConstr || false; 

            %@todo: For lsq/fsolve callers this have to be changed
            self.finDiffFlags.isGrad = true;
        end

        function self = setOptions(self, solver, algorithm, options)

            self.options = options;
            self.solver = solver;
            self.algorithm = algorithm;

            self.options.Vectorized = optimget(options,'Vectorized',self.defaultopt,self.optimgetFlag);
            self.options.UseParallel = optimget(options,'UseParallel',self.defaultopt,self.optimgetFlag);
            self.options.FunValCheck = optimget(options,'FunValCheck',self.defaultopt,self.optimgetFlag);
            self.options.DerivativeCheck = optimget(options,'DerivativeCheck',self.defaultopt,self.optimgetFlag);
            self.options.GradObj = optimget(self.options,'GradObj',self.defaultopt,self.optimgetFlag);
            self.options.GradConstr = optimget(self.options,'GradConstr',self.defaultopt,self.optimgetFlag);
            
        end

        function self = setX0AndBounds(self,X0,lb,ub)
            
            if ~isempty(X0)
                self.sizes.xShape = size(X0);
                self.sizes.nVar = numel(X0);
                
                if ~isempty(lb)
                    self.LB = lb(:);
                else
                    self.LB = -inf(self.sizes.nVar,1);
                end
                if ~isempty(ub)
                    self.UB = ub(:);
                else
                    self.UB = inf(self.sizes.nVar,1);
                end

                % if self.NeedGrads
                %     self.ObjEval.Grad = zeros(self.sizes.nVar,0);
                %     self.ConstrEval.IneqGrad = zeros(self.sizes.nVar,0);
                %     self.ConstrEval.EqGrad = zeros(self.sizes.nVar,0);
                % end
            end

        end
      
        function self = optimfcncheck(self)

            activeSet = 'active-set';
            sqp = 'sqp';
            trustRegionReflective = 'trust-region-reflective';
            %interiorPoint = 'interior-point';
            sqpLegacy = 'sqp-legacy';

            funValCheck = strcmpi(self.options.FunValCheck,'on');
            flags.grad = strcmpi(self.options.GradObj,'on');
            flags.constr = ~isempty(self.ConstraintFunctionCell);
            flags.gradconst = strcmpi(self.options.GradConstr,'on');

            % Notice that defaultopt.Hessian = [], so the variable "hessian" can be empty
            hessian = optimget(self.options,'Hessian',self.defaultopt,self.optimgetFlag);
            % If calling trust-region-reflective with an unavailable Hessian option value,
            % issue informative error message
            if strcmpi(self.algorithm,trustRegionReflective) && ...
                    ~( isempty(hessian) || strcmpi(hessian,'on') || strcmpi(hessian,'user-supplied') || ...
                    strcmpi(hessian,'off') || strcmpi(hessian,'fin-diff-grads')  )
                error(message('optimlib:fmincon:BadTRReflectHessianValue'))
            end

            flags.hess = ~iscell(hessian) && ...
                ( strcmpi(hessian,'user-supplied') || strcmpi(hessian,'on') );

            % Process objective function
            if ~isempty(self.ObjectiveFunctionCell)
                % will detect empty string, empty matrix, empty cell array
                % constrflag set to false because we're checking the
                % objective, not constraint.
                funfcn = optimfcnchk(self.ObjectiveFunctionCell,self.solver, ...
                    length(self.FunArgs.AdditionalParameters),funValCheck,flags.grad,flags.hess, ...
                    false,self.algorithm);                
            else
                funfcn = {'','','','',''};
            end

            % Process constraint function
            if flags.constr % Constraint is non-empty
                % hessflag set to false because hessian is never returned
                % by nonlinear constraint function. constrflag set to true
                % because we're checking the constraints.
                confcn = optimfcnchk(self.ConstraintFunctionCell,self.solver, ...
                    length(self.FunArgs.AdditionalParameters),funValCheck,flags.gradconst,false,true);                
            else
                flags.gradconst = false;
                confcn = {'','','','',''};
            end

            if strcmpi(self.solver,'fmincon') && ...
                    (strcmpi(self.algorithm,activeSet) || strcmpi(self.algorithm,sqp) || ...
                    strcmpi(self.algorithm,sqpLegacy))
                if flags.hess % conflicting options
                    flags.hess = false;
                    warning(message('optimlib:fmincon:HessianIgnoredForAlg', self.algorithm));
                    if strcmpi(funfcn{1},'fungradhess')
                        funfcn{1}='fungrad';
                    elseif  strcmpi(funfcn{1},'fun_then_grad_then_hess')
                        funfcn{1}='fun_then_grad';
                    end
                end
            end
            
            % Set the function evaluators for different funfcn{1}/confcn{1}
            % types
            self = self.setObjectiveEvaluator(funfcn);
            self = self.setConstraintEvaluator(confcn);
            
            % Store FunFcn for later use and solvers using this class just
            % for optimfcnchk.
            self.FunFcn = funfcn;
            self.ConFcn = confcn;
            self.Flags = flags;

        end

        function self = setObjectiveEvaluator(self, funfcn)
            switch funfcn{1}
                case 'fun'
                    self.ObjectiveFunAndGrad = @objfunEvaluator;
                    self.ObjectiveFinDiff = @(x) objfunFinDiffEvaluator(x,funfcn{3},self.FunArgs);
                    self.GradFlags.ObjfunFinDiff = true;

                case 'fungrad'
                    self.ObjectiveFunAndGrad = @objfungradEvaluator;
                    self.GradFlags.ObjfunHaveGrad = true;

                case 'fun_then_grad'
                    self.ObjectiveFunAndGrad = @objfunEvaluator;
                    self.ObjectiveGradient = @objgradEvaluator;

                case 'fungradhess'                    
                    self.ObjectiveFunAndGrad = @objfungradhessEvaluator;
                    self.GradFlags.ObjfunHaveGrad = true;
                    self.GradFlags.ObjfunHaveHessian = true;

                case 'fun_then_grad_then_hess'
                    self.ObjectiveFunAndGrad = @objfunEvaluator;
                    self.ObjectiveGradient = @objgradEvaluator;
                    self.Hessian = @objhessEvaluator;
                case ''
                    self.ObjectiveFunAndGrad = @(self,funfcn3,x,fgrad,hess) ... 
                        noObjective(self,funfcn3,x,zeros(0,1),fgrad,hess);
                otherwise
                    error(message('optimlib:fmincon:UndefinedCallType'));
            end
        end

        function self = setConstraintEvaluator(self, confcn)
            switch confcn{1}
                case 'fun'
                    self.ConstraintFunAndGrad = @constrfunEvaluator;
                    self.ConstraintFinDiff = @(x) constrfunFinDiffEvaluator(x,confcn{3},self.FunArgs);
                    self.GradFlags.NonlconFinDiff = true;

                case 'fungrad'
                    self.ConstraintFunAndGrad = @constrfungradEvaluator;
                    self.GradFlags.NonlconHaveGrad = true;

                case 'fun_then_grad'
                    self.ConstraintFunAndGrad = @constrfunEvaluator;
                    self.ConstraintGradient = @constrgradEvaluator;

                case ''
                    self.ConstraintFunAndGrad = @(self,confcn3,x,cingrad,ceqgrad) ... 
                        noConstraint(self,confcn3,x,zeros(0,1),zeros(0,1),cingrad,ceqgrad);
                otherwise
                    error(message('optimlib:fmincon:UndefinedCallType'));
            end
        end

        function [cingrad, ceqgrad ] = initConstrGrad(~, X, cin, ceq)
            cingrad = zeros(numel(X),numel(cin));
            ceqgrad = zeros(numel(X),numel(ceq));
        end

        function fgrad = initObjGrad(self, X, fval)
            if self.finDiffFlags.isGrad
                fgrad = zeros(numel(X),numel(fval));
            else
                fgrad = zeros(numel(fval),numel(X));
            end
        end        
        
    end

    methods(Abstract)
        self = setupParallel(self)       
    end

end

function [cin,ceq,cingrad,ceqgrad] = noConstraint(~,~,~,cin,ceq,cingrad,ceqgrad)
end

function [fval,fgrad,hess] = noObjective(~,~,~,fval,fgrad,hess)
end

function out = reshapeOutputs(out_,out)
    if isequal(out,[])
        out = out_;
    else
        out(:) = out_;
    end
end

function functionHandleOrEmpty(value)
if ~isempty(value) && ~isa(value,'function_handle')
    error(message('optimlib:fmincon:InvalidFUN'));
end
end

