classdef recursiveARMAX < idRecursivePolynomialModelEstimator
    %recursiveARMAX  Estimate ARMAX models using recursive estimation algorithms
    %
    %   Estimate the coefficients of the ARMAX models of the
    %   following form:
    %     A(q) y(t) = B(q) u(t-nk) +  C(q) e(t)
    %   where y(t) is the output measurement and u(t) is the input
    %   measurement. nk is the fixed time delay. A(q), B(q), and C(q) are
    %   polynomials of the time-delay operator 1/q, whose coefficients are
    %   estimated. e(t) is an unknown scalar that represents white noise.
    %
    %   recursiveARMAX can only estimate single-input single-output (SISO)
    %   models.
    %
    %   X = recursiveARMAX returns the recursiveARMAX estimator system
    %   object, X, to estimate 1 parameter in  A(q), B(q) and C(q) each.
    %   The initial guess for all parameters are 'eps'.
    %
    %   X = recursiveARMAX([na nb nc nk]) returns the recursiveARMAX
    %   estimator system object, X, to estimate na parameters in A(q), nb
    %   parameters in B(q) and nc parameters in C(q). na, nc and nk are
    %   non-negative scalars. nb must be positive.
    %
    %   X = recursiveARMAX([na nb nc nk],A0,B0,C0) returns the
    %   recursiveARMAX estimator system object, X, with initial guesses
    %   for the polynomial coefficients in A0, B0, and C0. A0 is a row
    %   vector of length na+1 and its first element must be 1. B0 is a row
    %   vector of length nb+nk with nk leading zeros. C0 must be a row
    %   vector of length nc+1. C0 must also be a stable discrete-time
    %   polynomial (all roots in the unit disk, i.e.
    %   all(abs(roots(C0))<1)==true).
    %
    %   X = recursiveARMAX([na nb nc nk],'PropertyName',PropertyValue)
    %   returns a recursiveARMAX estimator object, X, with each specified
    %   property set to the specified value.
    %
    %   step method syntax:
    %
    %   [A,B,C,yHat] = step(X, y, u) carries out the recursive algorithm
    %   over the measurements y(t) and u(t). A, B and C are the current
    %   parameter estimates for A(q), B(q) and  C(q) respectively. yHat is
    %   the estimate of y(t) based on A, B, C as well as the current and
    %   past values of y, u.
    %
    %   recursiveARMAX methods:
    %
    %   step     - See above description for use of this method
    %   release  - Allow property value and input characteristics changes
    %   clone    - Create another recursiveARMAX 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
    %
    %   recursiveARMAX properties:
    %
    %   InitialA                   - The initial polynomial coefficients of A(q)
    %   InitialB                   - The initial polynomial coefficients of B(q)
    %   InitialC                   - The initial polynomial coefficients of C(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 ARMAX model
    %      % Construct the estimator.
    %      X = recursiveARMAX([2 2 1 1]); % [na nb nc nk]
    %      % Provide initial parameter guesses. Optional, but recommended.
    %      X.InitialA = [1 -1.1 0.5];
    %      X.InitialB = [0 0.45 0.2];
    %      X.InitialC = [1 -0.7];
    %      % Provide variances of initial parameter guesses. Optional, but recommended.
    %      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,C,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 armax, recursiveARX, recursiveOE, recursiveBJ,
    % recursiveAR, recursiveARMA, recursiveLS    
    
    %   Copyright 2015-2018 The MathWorks, Inc.
    
    %#codegen
   
    methods
        function this = recursiveARMAX(varargin)
            this@idRecursivePolynomialModelEstimator(varargin{:});
        end
    end    
    
    methods(Access = protected)
        function setEstimator(this)
            coder.internal.errorIf(this.IsHistoryFinite_,...
                'Ident:simulink:recursiveEstimationFiniteHistoryNotSupported',...
                class(this));            
            
            this.Estimator_ = ...
                recursiveEstimation.internal.rpoly.rarmaxInfiniteHistory(...
                this.IsUsingFrames_,...
                this.getNumberOfInputs(this.nb),...
                this.getNumberOfOutputs(this.nb));
        end
        
        function setInitialParameters(this,tunablePropNVPairs,firstNVIndex,varargin)
            setInitialPolynomialCoefficients(this,'ABC',[2 3 4],tunablePropNVPairs,firstNVIndex,varargin{:});
        end
        
        function setInitialIOBuffers(~, ~)
            % Set initial input-output buffers during object construction
            
            % No op. Will be utilized with History='Finite' support
        end
        
        function validateOrdersPreSetImpl(~,val)
            % must have 4 non-negative integer elements (only SISO) in a
            % row vector. Nonnegative, integer, row checks are done in the
            % superclass
            validateattributes(val,{'numeric'},{'numel',4},'','Orders',1);
            % nb cannot contain 0s
            coder.internal.assert(all(val(2:numel(val)/2)>0),...
                'Ident:simulink:recursiveEstimationVectorNBCannotContain0');            
        end        
        
        function setOrdersImpl(this,val)            
            nInputs = (numel(val)-2)/2;
            this.na = val(1);
            this.nb = val(2:nInputs+1);
            this.nc = val(nInputs+2);
            this.nk = val(nInputs+2+(1:nInputs));
            this.nf = zeros(1,nInputs,'uint16');
        end        
        
        function [A,B,C,yHat] = stepImpl(this,y,u)       
            % Update estimates based on new y,u          
            
            [this.PastMeasurements, ... 
             this.PastJacobianEstimates, ...
             ~, ... % BufferInputs_
             ~, ... % BufferOutputs_
             this.Parameters_, ...
             this.ParameterCovariance_, ...
             A, B, C, e] = ...
                this.Estimator_.estimate(...
                    this.PastMeasurements, ...
                    this.PastJacobianEstimates, ...
                    cast(0,this.DataType), ... % BufferInputs_
                    cast(0,this.DataType), ... % BufferOutputs_
                    this.Parameters_, this.ParameterCovariance_, ...
                    y, u,... % new IO data
                    this.EnableAdaptation,...
                    getAdaptationParameter(this), ...
                    getAdaptationParameter2(this), ...
                    this.AlgorithmParameters_); % helper variables
                
            if nargout>=4
                yHat = y-e;
            end
        end
        
        function resetImpl(this)
            % Reset the states of the estimator
            resetImpl@idRecursivePolynomialModelEstimator(this);
            dataType = this.DataType;
            na_ = this.na;
            nb_ = this.nb;
            nc_ = this.nc;
            nk_ = this.nk;
            
            this.Parameters_ = this.InitialParameters_;
            if ~this.isInactivePropertyImpl('InitialParameterCovariance')
                this.ParameterCovariance_ = this.InitialParameterCovariance_;
            else
                this.ParameterCovariance_ = cast([],dataType);
            end
            this.PastMeasurements = ...
                zeros(1,recursiveARMAX.numberOfRequiredPastMeasurements(na_,nb_,nc_,nk_),dataType);
            this.PastJacobianEstimates = ...
                zeros(1,recursiveARMAX.numberOfRequiredPastJacobianEstimates(na_,nb_,nc_,nk_),dataType);
        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_ = recursiveARMAX.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;
            s.C = this.C;
        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 {'InitialD','InitialF','D','F','History','WindowLength'}
                    isInactive = true();
            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
            
            % Initial derivatives of the estimation error w.r.t the
            % estimated parameters
            [~,nam,nbm] = ...
                recursiveARMAX.numberOfRequiredPastJacobianEstimates(p.na,p.nb,p.nc,p.nk);
            % Data memory: phiMemory See localInitializeRarx() for detailed
            % explanation of phi. The elements in Phi that are required in
            % each time step
            necessaryDataForPhi = [(1:p.na) ...
                (p.na+p.nk+(1:p.nb)) ...
                (p.na+p.nk+p.nb+(1:p.nc))];
            % Derivatives memory: phiMemory The elements in Psi that are
            % required in each time step
            necessaryDataForPsi = [(1:p.na) ...
                (nam+p.nk+(1:p.nb)) ...
                (nam+nbm+(1:p.nc))];
            
            % Fcn outputs: Some of these are not needed for ARMA
            % estimation. However, the mex version of idrarmaxstep() requires
            % 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)+p.nc+p.nd+sum(p.nf,2));
            algorithmParams.necessaryDataForPhi = uint16(necessaryDataForPhi);
            algorithmParams.necessaryDataForPsi = uint16(necessaryDataForPsi);
            % B(q) parameters
            algorithmParams.b0PosInPhi = uint16(p.na+1);
            algorithmParams.b0PosInPsi = uint16(nam+1);
            % C(q) parameters
            algorithmParams.c1PosInTheta = uint16(p.na+sum(p.nb,2)+1);
            algorithmParams.c1PosInPhi = uint16(p.na+p.nb+p.nk+1);
            algorithmParams.c1PosInPsi = uint16(nam+nbm+1);
        end
        
        function [num,nam,nbm] = numberOfRequiredPastJacobianEstimates(na,nb,nc,nk)
            if na>0
                nam = max(na,nc);
            else
                nam = 0;
            end
            if nb>0
                nbm = max(nb+nk,nc+1);
            else
                nbm = 0;
            end
            num = nam+nbm+nc;
        end
        
        function num = numberOfRequiredPastMeasurements(na,nb,nc,nk)
            num = na+sum(nb+nk,2)+nc;
        end
    end
    
    methods(Static,Access=protected)
        function val = defaultOrders()
            val = [1 1 1 1];
        end
        
        function val = latestPossibleNVPairPosition()
            val = 5;
        end         
    end      
end

% LocalWords:  nk na nb nc nk na nb nc online ARMA
