classdef (Abstract) QPActiveSetWarmStart < optim.warmstart.base.LinConstrWarmStart
%QPACTIVESETWARMSTART Warm start base class for the Optimization Toolbox
%   active set algorithm.  Provides the backend for solving sequences of
%   quadratic programming problems.
%
%   See also OPTIMWARMSTART, OPTIMOPTIONS

%   Copyright 2020-2023 The MathWorks, Inc.

    properties (Constant, Access = protected)
        NumSolveOutputs = 4;
    end
    
    properties (Access = private)
        % Arrays containing 1 = active, 0 = inactive.  Stored as double
        % arrays for call to builtin function.
        Wineq_
        Wupper_
        Wlower_
                
        % Inequality constraints from the previous solve.
        Aold_
        bold_
        lbold_
        ubold_
    end
    
    properties (Access = protected)
        % MATLAB struct copy of options for call to builtin function.
        OptionsStruct_
    end
    
    methods
        
        % Extra inputs are for upper bounds on linear
        % equalities/inequalities if we need them.
        function obj = QPActiveSetWarmStart(x0, options, MaxLinearEqualities, MaxLinearInequalities)
            
            obj = obj@optim.warmstart.base.LinConstrWarmStart(x0, options, ...
                MaxLinearEqualities, MaxLinearInequalities);
            
            nVar = numel(x0);
                        
            % Denote that there is no prior history for this warm start
            % object to start from.
            obj.Wineq_ = false(0, 1);
            obj.Wupper_ = false(nVar, 1);
            obj.Wlower_ = false(nVar, 1);
            
            % Set initial arrays
            obj.Aold_ = [];
            obj.bold_ = [];
            obj.lbold_ = -Inf(nVar, 1);
            obj.ubold_ = Inf(nVar, 1);
            
            % Options format for the builtin driver.
            obj.OptionsStruct_.Display = options.Display;
            obj.OptionsStruct_.ConvexCheck = 'on';
            obj.OptionsStruct_.MaxIterations = options.MaxIterations;
            obj.OptionsStruct_.ConstraintTolerance = options.ConstraintTolerance;
            obj.OptionsStruct_.OptimalityTolerance = options.OptimalityTolerance;
            obj.OptionsStruct_.StepTolerance = options.StepTolerance;
            obj.OptionsStruct_.ObjectiveLimit = options.ObjectiveLimit;
            obj.OptionsStruct_.PricingTolerance = options.PricingTolerance;
        end
    end
    
    methods (Hidden)
        function out = PackUnconstrainedSolution(obj, X)
            out = obj;
            out.X = X;

            % Update old constraints and active set.
            out.Aold_ = [];
            out.bold_ = [];
            out.lbold_(:) = -Inf;
            out.ubold_(:) = Inf;
            
            out.Wineq_ = false(0, 1);
            out.Wlower_(:) = false;
            out.Wupper_(:) = false;
        end
    end
    
    methods (Access = protected)
               
        function [out, varargout] = performWarmStart(obj, varargin)
            out = obj;
            
            % Get the separate inputs
            [H, f, A, b, Aeq, beq, lb, ub] = deal(varargin{:});
            H = full(H);
            f = full(f(:));
            A = full(A);
            b = full(b(:));
            Aeq = full(Aeq);
            beq = full(beq(:));
            lb = full(lb(:));
            ub = full(ub(:));
            
            [out.Wineq_, out.Wlower_, out.Wupper_] = ...
                optim.warmstart.utils.qpactiveset.updateWorkingSet...
                (out.Wineq_, out.Wlower_, out.Wupper_, ...
                A, b, out.Aold_, out.bold_, lb, out.lbold_, ub, out.ubold_, out.OptionsStruct_.ConstraintTolerance);
            
            % Update constraints
            out.Aold_ = A;
            out.bold_ = b;
            out.lbold_ = lb;
            out.ubold_ = ub;

            % Update MaxIterations unless it is manually set.
            if strcmpi(out.OptionsStruct_.MaxIterations, '10*(numberOfVariables + numberOfConstraints)')
                nVar = numel(out.X);
                mConstr =  numel(b) + numel(beq) + sum(isfinite(lb)) + sum(isfinite(ub));
                out.OptionsStruct_.MaxIterations = 10*(nVar + mConstr);
            end
            
            % Convert for call into builtin.
            activeIneq = find(out.Wineq_);
            activeLB = find(out.Wlower_);
            activeUB = find(out.Wupper_);
            
            % Call warm start builtin
            [out.X, ~, activeIneq, activeLB, activeUB, fval, exitflag, lambda, mexoutput] = ...
                optim.warmstart.internal.activesetqpws(...
                H, f, A, b, Aeq, beq, lb, ub, out.X, ...
                out.OptionsStruct_, activeIneq, activeLB, activeUB);
            
            % Update active constraints
            out.Wineq_(:) = false;
            out.Wlower_(:) = false;
            out.Wupper_(:) = false;

            out.Wineq_(activeIneq) = true;
            out.Wlower_(activeLB) = true;
            out.Wupper_(activeUB) = true;
            
            % Needed for quadprog messaging
            mexoutput.MaxIterations = out.OptionsStruct_.MaxIterations;
            
            varargout = cell(out.NumSolveOutputs, 1);
            [varargout{:}] = deal(fval, exitflag, lambda, mexoutput);
        end
        
        function out = setConvex(obj)
            % we have a guarentee that these problems are convex
            out = obj;
            out.OptionsStruct_.ConvexCheck = 'off'; 
        end

    end
end

