classdef recursiveARX < idRecursivePolynomialModelEstimator
    %recursiveARX  Estimate ARX models using recursive estimation algorithms
    %
    %   Estimate the coefficients of the ARX models of the following form:
    %     A(q) y(t) = B(q) u(t-nk) + e(t)
    %   where y(t) is the output measurement and u(t) is the input
    %   measurement. nk is the fixed time delay. A(q) and B(q) are
    %   polynomials of the time-delay operator 1/q, whose coefficients are
    %   estimated. e(t) is an unknown scalar that represents white noise.
    %
    %   recursiveARX can estimate single-input single-output (SISO) and
    %   multiple-input single-output (MISO) models.
    %
    %   X = recursiveARX returns the recursiveARX estimator system
    %   object, X, to estimate 1 parameter in A(q) and B(q) each. The
    %   initial guess for all parameters are 'eps'.
    %
    %   X = recursiveARX([na nb nk]) returns the recursiveARX estimator
    %   system object, X, to estimate na parameters in A(q) and nb
    %   parameters in B(q). na is a scalar. nb and nk are scalars for SISO
    %   models. nb and nk are row vectors of the same length for MISO
    %   models. In this case, nb(i) and nk(i) correspond to the number of
    %   parameters in the B polynomial for the i-th input and the
    %   corresponding time-delay, respectively. nb cannot contain any
    %   zeros. na and all entries in nk must be non-negative.
    %
    %   X = recursiveARX([na nb nk],A0,B0) returns the recursiveARX
    %   estimator system object, X, with initial guesses for the polynomial
    %   coefficients in A0 and B0. A0 is a row vector of length
    %   na+1 and its first element must be 1. B0 is a matrix with nu
    %   (number of inputs) rows and max(nb+nk,[],2) columns. i-th row of B0
    %   must contain nk(i) zeros, followed by nb(i) parameter guesses.
    %   Entries beyond nk(i)+nb(i) are ignored.
    %
    %   X = recursiveARX([na nb nk],'PropertyName',PropertyValue) returns
    %   a recursiveARX estimator object, X, with each specified property
    %   set to the specified value.
    %
    %   step method syntax:
    %
    %   [A,B,yHat] = step(X, y, u) carries out the recursive algorithm over
    %   the measurements y(t) and u(t). A and B are the current parameter
    %   estimates for A(q) and B(q). yHat is the estimate of y(t) based on
    %   A, B as well as the current and past values of y, u.
    %
    %   recursiveARX methods:
    %
    %   step     - See above description for use of this method
    %   release  - Allow property value and input characteristics changes
    %   clone    - Create another recursiveARX System object with same property values
    %   isLocked - Locked status (logical)
    %   <a href="matlab:help matlab.System/reset   ">reset</a>    - Reset the internal states to initial conditions
    %
    %   recursiveARX properties:
    %
    %   InitialA                   - The initial polynomial coefficients of A(q)
    %   InitialB                   - The initial polynomial coefficients of B(q)
    %   EstimationMethod           - 'ForgettingFactor' (Default),
    %                                'KalmanFilter', 'Gradient' or
    %                                'NormalizedGradient'
    %   ForgettingFactor           - Forgetting factor, only used by the
    %                                ForgettingFactor algorithm
    %   ProcessNoiseCovariance     - The covariance of the process noise acting
    %                                on the parameter values. Only used by the
    %                                KalmanFilter algorithm
    %   AdaptationGain             - Adaptation gain, only used by the
    %                                Gradient and NormalizedGradient
    %                                algorithms
    %   NormalizationBias          - Bias on adaptation gain scaling. Only
    %                                used by the NormalizedGradient algorithm
    %   InitialParameterCovariance - Initial value for covariance of the
    %                                parameter estimation error
    %   EnableAdaptation           - Can be used to enable/disable the
    %                                parameter (theta) adaptation
    %
    %   EXAMPLE: Estimating a time-invariant ARX model
    %      % Construct the estimator.
    %      X = recursiveARX([3 3 1]); % [na nb nk]
    %      % Provide initial parameter guesses (optional).
    %      X.InitialA = [1 -0.8 0.2 0.1];
    %      X.InitialB = [0 0.5 0.4 0.2];
    %      % Provide variances of initial parameter guesses (optional).
    %      X.InitialParameterCovariance = 0.2;
    %      % Load the estimation data.
    %      load iddata1 z1;
    %      % Call the estimator's step command with one input-output data pair at a
    %      % time to update parameters. This illustrates online operation of the
    %      % estimator.
    %      for kk=1:numel(z1.u);
    %          [A,B,yHat] = step(X,z1.y(kk),z1.u(kk));
    %      end
    %      % Take a snapshot of the current model parameters.
    %      sysX = idpoly(X);
    %      % Set the sample time of the model snapshot from the data.
    %      sysX.Ts = z1.Ts;
    %      % Compare
    %      compare(z1,sysX);
    %
    %
    % See also arx, recursiveARMAX, recursiveOE, recursiveBJ,
    % recursiveAR, recursiveARMA, recursiveLS
    
    %   Copyright 2015-2018 The MathWorks, Inc.
    
    %#codegen
    
    properties(Dependent)
        %InitialOutputs The initial (past) values of outputs
        %   This property is only used when the History property is set to
        %   'Finite'.
        %
        %   This property stores the initial guesses for the output data.
        %   InitialOutputs is utilized at the at the very first estimation
        %   step and when the estimator is reset. It is used together with
        %   InitialInputs.
        %
        %   InitialOutputs must be a column vector with WindowLength+na-1
        %   elements. The default value is a vector of zeros.
        InitialOutputs
        
        %InitialInputs The initial (past) values of regressors
        %   This property is only used when the History property is set to
        %   'Finite'.
        %
        %   This property stores the initial guesses for the input data.
        %   InitialInputs is utilized at the at the very first estimation
        %   step and when the estimator is reset. It is used together with
        %   InitialOutputs.
        %
        %   InitialInputs must be a matrix, with WindowLength+max(nb(:))+max(nk(:)))
        %   rows, and as many columns as the number of inputs in the
        %   system. The default value is a matrix of zeros.
        InitialInputs
    end
    
    properties(Access=protected)
        InitialOutputs_;
        InitialInputs_;
        PastOutputs_;
        PastInputs_;      
    end
    
    methods
        function this = recursiveARX(varargin)
            this@idRecursivePolynomialModelEstimator(varargin{:});
        end
    end
    
    methods
        % get & set methods for dependent properties
        function set.InitialInputs(this,userValue)
            % This property can only be set if History='Finite'. Otherwise
            % its value is unused.
            %
            % An error is thrown because History option is nontunable,
            % hence any value set here is guaranteed to be not used.
            coder.internal.assert(this.IsHistoryFinite_,...
                'Ident:simulink:recursiveEstimationPropOnlyForFiniteHistory',...
                'InitialInputs');
            
            % Validate
            propertyName = 'InitialInputs';
            expectedSize = [...
                recursiveEstimation.internal.rpoly.rarxFiniteHistory.getInputBufferLength(...
                this.WindowLength,...
                this.nb,...
                this.nk) cast(this.getNumberOfInputs(this.nb), 'like', this.WindowLength)];            
            if isempty(coder.target)
                try 
                    validatedValue = controllib.internal.blocks.rls.validateInitialIOData(...
                        userValue,expectedSize,propertyName);
                catch ME
                    throw(ME)
                end
            else
                validatedValue = controllib.internal.blocks.rls.validateInitialIOData(...
                    userValue,expectedSize,propertyName);
            end
            
            % Assign
            this.InitialInputs_ = cast(validatedValue, this.DataType);
        end
        
        function val = get.InitialInputs(this)
            val = this.InitialInputs_;
        end
        
        function set.InitialOutputs(this,userValue)
            % This property can only be set if History='Finite'. Otherwise
            % its value is unused.
            %
            % An error is thrown because History option is nontunable,
            % hence any value set here is guaranteed to be not used.
            coder.internal.assert(this.IsHistoryFinite_,...
                'Ident:simulink:recursiveEstimationPropOnlyForFiniteHistory',...
                'InitialOutputs');
            % Validate
            propertyName = 'InitialOutputs';
            expectedSize = [...
                recursiveEstimation.internal.rpoly.rarxFiniteHistory.getOutputBufferLength(...
                this.WindowLength,...
                this.na) cast(this.getNumberOfOutputs(this.nb), 'like', this.WindowLength)];                
            if isempty(coder.target)
                try 
                    validatedValue = controllib.internal.blocks.rls.validateInitialIOData(...
                        userValue,expectedSize,propertyName);
                catch ME
                    throw(ME)
                end
            else
                validatedValue = controllib.internal.blocks.rls.validateInitialIOData(...
                    userValue,expectedSize,propertyName);
            end
            % Assign
            this.InitialOutputs_ = cast(validatedValue,this.DataType);
        end
        
        function val = get.InitialOutputs(this)
            val = this.InitialOutputs_;
        end
    end
    
    methods(Access = protected)
        function setEstimator(this)
            if this.IsHistoryFinite_
                this.Estimator_ = ...
                    recursiveEstimation.internal.rpoly.rarxFiniteHistory(...
                    this.IsUsingFrames_,...
                    this.WindowLength_,...
                    this.na,...
                    this.nb,...
                    this.nk,...
                    class(this)); % name of the estimator object, for warnings
            else
                this.Estimator_ = ...
                    recursiveEstimation.internal.rpoly.rarxInfiniteHistory(...
                    this.IsUsingFrames_,...
                    this.getNumberOfInputs(this.nb),...
                    this.getNumberOfOutputs(this.nb));
            end
        end
        
        
        function setInitialParameters(this,tunablePropNVPairs,firstNVIndex,varargin)
            setInitialPolynomialCoefficients(this,'AB',[2 3],tunablePropNVPairs,firstNVIndex,varargin{:});
        end
        
        function setInitialIOBuffers(this,tunablePropNVPairs)
            % Set initial input-output buffers during object construction
            
            % These properties are only used with History='Finite'
            if ~this.IsHistoryFinite_
                return;
            end
            
            % Assign InitialOutputs
            if ~isempty(tunablePropNVPairs.InitialOutputs)
                this.InitialOutputs = cast(tunablePropNVPairs.InitialOutputs, this.DataType);
            else
                this.InitialOutputs = cast(0, this.DataType); % scalar will be expanded in .set method
            end
            % Assign InitialInputs
            if ~isempty(tunablePropNVPairs.InitialInputs)
                this.InitialInputs = cast(tunablePropNVPairs.InitialInputs, this.DataType);
            else
                this.InitialInputs = cast(0, this.DataType); % scalar will be expanded in .set method
            end
        end
        
        function validateOrdersPreSetImpl(~,val)
            % must have 2*n+1 non-negative integer elements. Nonnegative
            % and integer checks are done in the superclass
            numelOrders = numel(val);
            coder.internal.errorIf(numelOrders<3 || mod(numelOrders,2)~=1,...
                'Ident:simulink:recursiveEstimationRARXInvalidOrdersDim');
            % nb cannot contain 0s
            coder.internal.assert(all(val(2:(numel(val)+1)/2)>0),...
                'Ident:simulink:recursiveEstimationVectorNBCannotContain0');
        end
        
        function setOrdersImpl(this,val)
            nInputs = (numel(val)-1)/2;
            this.na = val(1);
            this.nb = val(2:nInputs+1);
            this.nk = val(nInputs+1+(1:nInputs));
            this.nf = zeros(1,nInputs,'uint16');
        end
        
        function [A,B,yHat] = stepImpl(this,y,u)
            % Update estimates based on new y,u
            
            [this.PastMeasurements, ...
                this.PastInputs_, ...
                this.PastOutputs_, ...
                this.Parameters_, ...
                this.ParameterCovariance_, ...
                A, B, e] = ...
                this.Estimator_.estimate(...
                this.PastMeasurements, ...
                this.PastInputs_, this.PastOutputs_, ... % States
                this.Parameters_, this.ParameterCovariance_, ... % Outputs
                y, u,... % new IO data
                this.EnableAdaptation,...
                getAdaptationParameter(this), ...
                getAdaptationParameter2(this), ...
                this.AlgorithmParameters_); % helper variables

            if nargout>=3
                yHat = y-e;
            end
        end
        
        function resetImpl(this)
            % Reset the states of the estimator

            % Note: Initializing states, even if they are empty, is
            % necessary for code generation
            
            % PastJacobianEstimates: Not used by ARX
            this.PastJacobianEstimates = cast([],this.DataType);            
            
            if this.IsHistoryFinite_
                % States: 
                % * PastOutputs, PastInputs
                % * Parameters (only impacts initial residuals/output
                % predictions)
                %
                % Outputs, but not states:
                % * ParameterCovariance
                %
                % Unused:
                % * PastMeasurements
                
                % States
                %
                % PastOutputs
                this.PastOutputs_ = this.InitialOutputs;
                % PastInputs: Not used for time-series models
                if this.getNumberOfInputs(this.nb) > 0
                    this.PastInputs_ = this.InitialInputs;
                else
                    this.PastInputs_ = cast([],this.DataType);
                end
                               
                % Outputs that are not states
                this.Parameters_ = this.InitialParameters_;
                this.ParameterCovariance_ = this.InitialParameterCovariance_;
                      
                % Unused properties
                this.PastMeasurements = cast([],this.DataType);
            else               
                % States:
                % * Parameters
                % * PastMeasurements
                % * ParameterCovariance (FF, KF algorithms, otherwise unused)                
                % 
                % Outputs (and also states):
                % * Parameters
                % * ParameterCovariance (for FF, KF algorithms)
                %
                % Outputs, but not states:
                % * N/A
                %
                % Unused:
                % * PastOutputs, PastInputs

                % States (some of which are also outputs)
                %
                % Parameters
                this.Parameters_ = this.InitialParameters_;
                % PastMeasurements: Always initialized to zero
                this.PastMeasurements = ...
                    zeros(1,recursiveARX.numberOfRequiredPastMeasurements(this.na,this.nb,this.nk),this.DataType);                   
                % ParameterCovariance: Only used with FF, KF algorithms
                if isEstimationMethodForgettingFactor(this) || isEstimationMethodKalmanFilter(this)
                    this.ParameterCovariance_ = this.InitialParameterCovariance_;
                else
                    % Unused
                    this.ParameterCovariance_ = cast([],this.DataType);
                end

                % Unused
                this.PastOutputs_ = cast([],this.DataType);
                this.PastInputs_ = cast([],this.DataType);
            end

            % Reset the estimator internal states, update alg. outputs
            [this.Parameters_,this.ParameterCovariance_] = ...
                this.Estimator_.reset(...
                    this.PastMeasurements,...
                    this.PastInputs_,...
                    this.PastOutputs_,...
                    this.Parameters_,...
                    this.ParameterCovariance_);
        end
        
        function setupImpl(this)
            % One-time setup: Construct a structure containing all algorithm parameters
            p = struct(...
                'na', this.na, ...
                'nb', this.nb,...
                'nc', this.nc,...
                'nd', this.nd,...
                'nf', this.nf,...
                'nk', this.nk,...
                'estimationMethod',this.EstimationMethod_);
            this.AlgorithmParameters_ = recursiveARX.getAlgorithmParameters(p);
        end
        
        function s = getDiscreteStateImpl(this)
            % Present estimator states to the user
            s = getDiscreteStateImpl@idRecursivePolynomialModelEstimator(this);
            s.A = this.A;
            s.B = this.B;
        end
        
        function isInactive = isInactivePropertyImpl(this,propName)
            % Hide the unused properties of the superclass
            isInactive = isInactivePropertyImpl@idRecursivePolynomialModelEstimator(this,propName);
            if isInactive
                return;
            end
            switch propName
                case {'InitialC','InitialD','InitialF','C','D','F'}
                    isInactive = true();
                case {'InitialParameterCovariance'}
                    if this.IsHistoryFinite_
                        isInactive = true();
                    end
                case {'InitialOutputs','InitialInputs'}
                    if ~this.IsHistoryFinite_
                        isInactive = true();
                    end
            end
        end
    end
    
    methods(Static,Hidden)
        function algorithmParams = getAlgorithmParameters(p)
            % getAlgorithmParameters Pre-calculate the required parameters
            % for the estimation algorithm. This is used by this system
            % object and the recursive polynomial model estimator simulink
            % block.
            %
            %   algorithmParams = getAlgorithmParameters(p)
            %
            %   Inputs;
            %     p - A structure containing validated parameters
            
            % Data memory: phiMemory
            % At time k, the phiMemory vector is:
            % phiMemory[k] =  [-y(k-1)  -y(k-2) ... -y(k-na)
            %                  u_1(k) u_1(k-1) u_1(k-2) ... u_1(k-nb_1-nk_1+1)
            %                  u_2(k) u_2(k-1) u_2(k-2) ... u_2(k-nb_2-nk_2+1)
            %                  ...
            %                  u_N(k) u_N(k-1) u_N(k-2) ... u_N(k-nb_N-nk_N+1)];
            % where N is the total number of inputs. The subscript j=1,...,N denotes
            % the j-th input. nb_j and nk_j are the order of the denominator polynomial
            % and the delay associated with the j-th input.
            %
            % Only some data in phiMemory is required at time k for regression.
            % These are:
            % phi = [-y(k-1) -y(k-2) ... -y(k-na)
            %         u_j(k-nk_j) u_j(k-nk_j-1) u_j(k-nk_j-2) ... u_j(k-nb_j-nk_j+1) ...]
            % for _j=1,...,N. The remaining values in phi are required for future time
            % steps.
            %
            % We need to do two tasks at every time step:
            % a) Get phi from phiMemory. The variable necessaryDataForPhi captures the
            % indices of the necessary data in phi such that phi = phiMemory(necessaryDataForPhi)
            %
            % b) Update phiMemory with the data that becomes available at time k:
            % b1) Incorporate u(k) in phiMemory(finalInputPos)
            % b2) Execute the estimation step
            % b3) Shift the phiMemory vector to the right one step
            % b4) Incorporate -y(k) in phiMemory(1)
            %
            % b0PosInPhi holds the position of u_j(k), j=1,..,N in phiMemory. It is:
            % na+1+[nk1:nk1+nb1-1, nk1+nb1+nk2:nk1+nb1+nk2+nb2-1,...]
            nInputs = idRecursivePolynomialModelEstimator.getNumberOfInputs(p.nb);
            if nInputs>0
                b0PosInPhi = p.na+1+cumsum([0 p.nb(1:end-1)+p.nk(1:end-1)]);
                nbCumSum = cumsum([1 p.nb]);
                % Calculate:
                %   activeInputDataIndicesInPhi - Data necessary in phi, that are coming
                %                                 from the system inputs
                %   freeParameterIndicesInB     - The linear indices to the free parameters
                %                                 in the B polynomial matrix output This is
                %                                 only used when the initial parameters are
                %                                 supplied externally. Each column of the B
                %                                 polynomial contain nk_i leading zeros,
                %                                 nb_i free parameters, and some extra
                %                                 zeros for padding so that the all of the
                %                                 polynomial coefficients can be stored in
                %                                 a single matrix.
                activeInputDataIndicesInPhi = zeros(1, sum(p.nb,2)); % sum(nb,2)=# of estimated parameters in the numerators (assume single output model)
                for kk=1:nInputs
                    for kk2=0:p.nb(kk)-1
                        activeInputDataIndicesInPhi(nbCumSum(kk)+kk2) = b0PosInPhi(kk)+p.nk(kk)+kk2;
                    end
                end
            else
                b0PosInPhi = 0;
                activeInputDataIndicesInPhi = [];
            end
            necessaryDataForPhi = [(1:p.na) activeInputDataIndicesInPhi];
            
            % Fcn outputs: The data required inside the system identification block
            % Some of these are not needed for ar estimation. However, the mex version
            % of idrarxstep() require that algorithmParams have all these fields defined
            % at all times.
            algorithmParams.estimationMethod = uint8(p.estimationMethod);
            algorithmParams.na = uint16(p.na);
            algorithmParams.nb = uint16(p.nb);
            algorithmParams.nc = uint16(p.nc);
            algorithmParams.nd = uint16(p.nd);
            algorithmParams.nf = uint16(p.nf);
            algorithmParams.nk = uint16(p.nk);
            algorithmParams.nParameters = uint16(p.na+sum(p.nb,2));
            algorithmParams.necessaryDataForPhi = uint16(necessaryDataForPhi); % TBD - Lower mem req
            % B(q) parameters
            algorithmParams.b0PosInPhi = uint16(b0PosInPhi);
        end
        
        function num = numberOfRequiredPastMeasurements(na,nb,nk)
            num = na+sum(nb+nk,2);
        end
    end
    
    methods(Access=protected)
        function val = defaultOrders(~)
            val = [1 1 1];
        end
        
        function val = latestPossibleNVPairPosition(~)
            val = 4;
        end
    end
    
    methods(Hidden)
        % Methods for internal testing
        function s = getProtectedProperties(this)
            % Override base class method to have access to properties
            % defined in this method
            
            % Get properties that has protected or private GetAccess
            m = metaclass(this);
            s = struct();
            for kk=1:numel(m.PropertyList)
                if ~strcmp(m.PropertyList(kk).GetAccess,'public')
                    s.(m.PropertyList(kk).Name) = this.(m.PropertyList(kk).Name);
                end
            end
        end
    end    
end

% LocalWords:  nk MISO na nb nk MISO na nb nontunable online ARMA idrarxstep nd
