classdef recursiveBJ < idRecursivePolynomialModelEstimator
    %recursiveBJ  Estimate Box-Jenkins models using recursive estimation algorithms
    %
    %   Estimate the coefficients of the Box-Jenkins (BJ) models of the
    %   following form:
    %     y(t) = [B(q)/F(q)] u(t-nk) +  [C(q)/D(q)]e(t)
    %   where y(t) is the output measurement and u(t) is the input
    %   measurement. nk is the fixed time delay. B(q), C(q), D(q) and F(q)
    %   are polynomials of the time-delay operator 1/q, whose coefficients
    %   are estimated. e(t) is an unknown scalar that represents white
    %   noise.
    %
    %   recursiveBJ can only estimate single-input single-output (SISO)
    %   models.
    %
    %   X = recursiveBJ returns the recursiveBJ estimator system
    %   object, X, to estimate 1 parameter in  B(q), C(q), D(q) and F(q)
    %   each. The initial guess for all parameters are 'eps'.
    %
    %   X = recursiveBJ([nb nc nd nf nk]) returns the recursiveBJ estimator
    %   system object, X, to estimate nb parameters in B(q), nc parameters
    %   in C(q), nd parameters in D(q) and nf parameters in F(q). nc, nd
    %   and nf are non-negative scalars. nb and nk are positive scalars.
    %
    %   X = recursiveBJ([nb nc nd nf nk],B0,C0,D0,F0) returns the
    %   recursiveBJ estimator system object, X, with initial guesses for
    %   the polynomial coefficients in B0, C0, D0, F0. B0 is a row vector
    %   of length nb+nk with nk leading zeros. C0, D0 and F0 must all be
    %   row vectors of length nc+1, nd+1 and nf+1 respectively. C0, D0 and
    %   F0 must also be stable discrete-time polynomials (all roots in the
    %   unit disk, i.e. all(abs(roots(C0))<1)==true).
    %
    %   X = recursiveBJ([nb nc nd nf nk],'PropertyName',PropertyValue)
    %   returns a recursiveBJ estimator object, X, with each specified
    %   property set to the specified value.
    %
    %   step method syntax:
    %
    %   [B,C,D,F,yHat] = step(X, y, u) carries out the recursive algorithm
    %   over the measurements y(t) and u(t). B, C, D and F are the current
    %   parameter estimates for B(q), C(q), D(q) and F(q) respectively.
    %   yHat is the estimate of y(t) based on B, C, D, F as well as the
    %   current and past values of y, u.
    %
    %   recursiveBJ methods:
    %
    %   step     - See above description for use of this method
    %   release  - Allow property value and input characteristics changes
    %   clone    - Create another recursiveBJ 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
    %
    %   recursiveBJ properties:
    %
    %   InitialB                   - The initial polynomial coefficients of B(q)
    %   InitialC                   - The initial polynomial coefficients of C(q)
    %   InitialD                   - The initial polynomial coefficients of D(q)
    %   InitialF                   - The initial polynomial coefficients of F(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 Box-Jenkins model
    %      % Construct the estimator.
    %      X = recursiveBJ([2 1 1 2 1]); % [nb nc nd nf nk]
    %      % Provide initial parameter guesses. Optional, but recommended.
    %      X.InitialB = [0 0.8 0.2];
    %      X.InitialC = [1 0.1];
    %      X.InitialD = [1 -0.3];
    %      X.InitialF = [1 -1.1 0.6];
    %      % Provide variances of initial parameter guesses. Optional, but recommended.
    %      X.InitialParameterCovariance = 0.1;
    %      % 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);
    %          [B,C,D,F,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 bj, recursiveARX, recursiveARMAX, recursiveOE,
    % recursiveAR, recursiveARMA, recursiveLS
    
    %   Copyright 2015-2018 The MathWorks, Inc.
    
    %#codegen
    
    methods
        function this = recursiveBJ(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.rbjInfiniteHistory(...
                this.IsUsingFrames_,...
                this.getNumberOfInputs(this.nb),...
                this.getNumberOfOutputs(this.nb));
        end
        
        function setInitialParameters(this,tunablePropNVPairs,firstNVIndex,varargin)
            setInitialPolynomialCoefficients(this,'BCDF',[2 3 4 5],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 5 non-negative integer elements (only SISO) in a
            % row vector. Nonnegative, integer, row checks are done in the
            % superclass
            validateattributes(val,{'numeric'},{'numel',5},'','Orders',1);
            % nk must be at least 1 for BJ
            nInputs = (numel(val)-2)/3;
            coder.internal.assert(all(val(2*nInputs+2+(1:nInputs))>0),...
                'Ident:simulink:recursiveEstimationRBJNkLowerBound');
            % nb cannot contain 0s
            coder.internal.assert(all(val(1:nInputs)>0),...
                'Ident:simulink:recursiveEstimationVectorNBCannotContain0');
        end
        
        function setOrdersImpl(this,val)
            nInputs = (numel(val)-2)/3;
            this.nb = val(1:nInputs);
            this.nc = val(nInputs+1);
            this.nd = val(nInputs+2);
            this.nf = val(nInputs+2+(1:nInputs));
            this.nk = val(2*nInputs+2+(1:nInputs));
        end
        
        function [B,C,D,F,yHat] = stepImpl(this,y,u)
            % Update estimates based on new y,u

            [this.PastMeasurements, ... 
             this.PastJacobianEstimates, ...
             ~, ... % BufferInputs_
             ~, ... % BufferOutputs_
             this.Parameters_, ...
             this.ParameterCovariance_, ...
             B, C, D, F, 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>=5
                yHat = y-e;
            end
        end
        
        function resetImpl(this)
            % Reset the states of the estimator
            resetImpl@idRecursivePolynomialModelEstimator(this);
            dataType = this.DataType;
            nb_ = this.nb;
            nc_ = this.nc;
            nd_ = this.nd;
            nf_ = this.nf;
            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,recursiveBJ.numberOfRequiredPastMeasurements(nb_,nc_,nd_,nf_,nk_),dataType);
            this.PastJacobianEstimates = ...
                zeros(1,recursiveBJ.numberOfRequiredPastJacobianEstimates(nb_,nc_,nd_,nf_,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_ = recursiveBJ.getAlgorithmParameters(p);
        end
        
        function s = getDiscreteStateImpl(this)
            % Present estimator states to the user
            s = getDiscreteStateImpl@idRecursivePolynomialModelEstimator(this);
            s.B = this.B;
            s.C = this.C;
            s.D = this.D;
            s.F = this.F;
        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 {'InitialA','A','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
            [~,nbm,ndm] = ...
                recursiveBJ.numberOfRequiredPastJacobianEstimates(p.nb,p.nc,p.nd,p.nf,p.nk);
            % Data memory: phiMemory
            % See localInitializeRarx() for detailed explanation of phi.
            % The elements in Phi that is required in the current time step
            necessaryDataForPhi = [(p.nk:p.nb+p.nk-1) ...
                (nbm+1:nbm+p.nc) ...
                (nbm+p.nc+1:nbm+p.nc+p.nd) ...
                (nbm+p.nc+ndm+1:nbm+p.nc+ndm+p.nf)];
            % Derivatives memory: phiMemory
            % The elements in Psi that is required in the current time step
            necessaryDataForPsi = necessaryDataForPhi;
            
            % Fcn outputs: The data required inside the system identification block
            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(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(1);
            algorithmParams.b0PosInPsi = uint16(1);
            % C(q) parameters
            algorithmParams.c1PosInTheta = uint16(p.nb+1);
            algorithmParams.c1PosInPhi = uint16(nbm+1);
            algorithmParams.c1PosInPsi = uint16(nbm+1);
            % D(q) parameters
            algorithmParams.d1PosInTheta = uint16(p.nb+p.nc+1);
            algorithmParams.d1PosInPhi = uint16(nbm+p.nc+1);
            algorithmParams.d1PosInPsi = uint16(nbm+p.nc+1);
            % F(q) parameters
            algorithmParams.f1PosInTheta = uint16(p.nb+p.nc+p.nd+1);
            algorithmParams.fnfPosInTheta = uint16(p.nb+p.nc+p.nd+p.nf);
            algorithmParams.f1PosInPhi = uint16(nbm+p.nc+ndm+1);
            algorithmParams.f1PosInPsi = uint16(nbm+p.nc+ndm+1);
        end
        
        function [num,nbm,ndm] = numberOfRequiredPastJacobianEstimates(nb,nc,nd,nf,nk)
            nbm = max([nb+nk-1 nf+nc nd]);
            ndm = max(nc,nd);
            nfm = max(nd,nf+nc);
            num = nbm+ndm+nfm+nc;
        end
        
        function num = numberOfRequiredPastMeasurements(nb,nc,nd,nf,nk)
            nbm = max([nb+nk-1 nf+nc nd]);
            ndm = max(nc,nd);
            nfm = max(nd,nf+nc);
            num = nbm+ndm+nfm+nc;
        end
    end
    
    methods(Static,Access=protected)
        function val = defaultOrders()
            val = [1 1 1 1 1];
        end
        
        function val = latestPossibleNVPairPosition()
            val = 6;
        end
    end
end

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