classdef recursiveOE < idRecursivePolynomialModelEstimator
    %recursiveOE  Estimate Output-Error models using recursive estimation algorithms
    %
    %   Estimate the coefficients of the Output-Error (OE) models of the
    %   following form:
    %     y(t) = [B(q)/F(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. B(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.
    %
    %   recursiveOE can only estimate single-input single-output (SISO)
    %   models.
    %
    %   X = recursiveOE returns the recursiveOE estimator system
    %   object, X, to estimate 1 parameter in  B(q) and F(q) each. The
    %   initial guess for all parameters are 'eps'.
    %
    %   X = recursiveOE([nb nf nk]) returns the recursiveOE estimator
    %   system object, X, to estimate nb parameters in B(q) and nf
    %   parameters in F(q). nb and nk are positive scalars. nf is a
    %   non-negative scalar.
    %
    %   X = recursiveOE([nb nf nk],B0,F0) returns the recursiveOE
    %   estimator system object, X, with initial guesses for the polynomial
    %   coefficients in B0 and F0. B0 is a row vector of length nb+nk with
    %   nk leading zeros. F0 must be a row vectors of length nf+1
    %   respectively. F0 must also be stable discrete-time polynomials (all
    %   roots in the unit disk, i.e. all(abs(roots(F0))<1)==true).
    %
    %   X = recursiveOE([nb nf nk],'PropertyName',PropertyValue)
    %   returns a recursiveOE estimator object, X, with each specified
    %   property set to the specified value.
    %
    %   step method syntax:
    %
    %   [B,F,yHat] = step(X, y, u) carries out the recursive algorithm over
    %   the measurements y(t) and u(t). B and F are the current parameter
    %   estimates for B(q) and F(q) respectively. yHat is the estimate of
    %   y(t) based on B, F as well as the current and past values of y, u.
    %
    %   recursiveOE methods:
    %
    %   step     - See above description for use of this method
    %   release  - Allow property value and input characteristics changes
    %   clone    - Create another recursiveOE 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
    %
    %   recursiveOE properties:
    %
    %   InitialB                   - The initial polynomial coefficients of B(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 Output-Error model
    %      % Construct the estimator.
    %      X = recursiveOE([2 2 1]); % [nb nf nk]
    %      % Provide initial parameter guesses. Optional, but recommended.
    %      X.InitialB = [0 -0.5 0.9];
    %      X.InitialF = [1 -0.7 0.1];
    %      % Provide variances of initial parameter guesses. Optional, but recommended.
    %      X.InitialParameterCovariance = 0.2;
    %      % Load the estimation data.
    %      load iddata3 z3;
    %      % 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(z3.u);
    %          [B,F,yHat] = step(X,z3.y(kk),z3.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 = z3.Ts;
    %      % Compare
    %      compare(z3,sysX);
    %
    %
    % See also oe, recursiveARX, recursiveARMAX, 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+nf-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+nk,[],'all')
        %   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 = recursiveOE(varargin)
            this@idRecursivePolynomialModelEstimator(varargin{:});
        end
        
        % 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.roeFiniteHistory.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.roeFiniteHistory.getOutputBufferLength(...
                this.WindowLength,...
                this.nf) 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.roeFiniteHistory(...
                    this.IsUsingFrames_,...
                    this.WindowLength_,...
                    this.nb,...
                    this.nf,...
                    this.nk,...
                    class(this)); % name of the estimator object, for warnings
            else
                this.Estimator_ = ...
                    recursiveEstimation.internal.rpoly.roeInfiniteHistory(...
                    this.IsUsingFrames_,...
                    this.getNumberOfInputs(this.nb),...
                    this.getNumberOfOutputs(this.nb));
            end
        end        
        
        function setInitialParameters(this,tunablePropNVPairs,firstNVIndex,varargin)
            setInitialPolynomialCoefficients(this,'BF',[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(this,val)
            % must have 3 non-negative integer elements (only SISO) in a
            % row vector. Nonnegative, integer, row checks are done in the
            % superclass
            validateattributes(val,{'numeric'},{'numel',3},'','Orders',1);
            nInputs = numel(val)/3;
            % nk must be at least 1 for OE, for infinite history estimator
            if ~this.IsHistoryFinite_        
                coder.internal.assert(all(val(2*nInputs+(1:nInputs))>0),...
                'Ident:simulink:recursiveEstimationROENkLowerBound');
            end
            % nb cannot contain 0s
            coder.internal.assert(all(val(1:nInputs)>0),...
                'Ident:simulink:recursiveEstimationVectorNBCannotContain0');
        end
        
        function setOrdersImpl(this,val)
            nInputs = numel(val)/3;
            this.nb = val(1:nInputs);
            this.nf = val(nInputs+(1:nInputs));
            this.nk = val(2*nInputs+(1:nInputs));
        end
        
        function [B,F,yHat] = stepImpl(this,y,u)
            % Update estimates based on new y,u
            
            [this.PastMeasurements, ...
                this.PastJacobianEstimates, ...
                this.PastInputs_, ...
                this.PastOutputs_, ...
                this.Parameters_, ...
                this.ParameterCovariance_, ...
                B, F, e] = ...
                this.Estimator_.estimate(...
                this.PastMeasurements, ...
                this.PastJacobianEstimates, ...
                this.PastInputs_, this.PastOutputs_, ... % States
                this.Parameters_, this.ParameterCovariance_, ...
                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
            resetImpl@idRecursivePolynomialModelEstimator(this);
            dataType = this.DataType;
            
            if this.IsHistoryFinite_
                % States: 
                % * PastOutputs, PastInputs
                %
                % Outputs, but not states:
                % * Parameters
                % * ParameterCovariance
                %
                % Unused:
                % * PastMeasurements
                % * PastJacobianEstimates
                
                % States
                %
                % PastOutputs, PastInputs
                this.PastOutputs_ = this.InitialOutputs;
                this.PastInputs_ = this.InitialInputs;
                               
                % Outputs that are not states
                this.Parameters_ = this.InitialParameters_;
                this.ParameterCovariance_ = this.InitialParameterCovariance_;
                      
                % Unused properties
                this.PastMeasurements = cast([], dataType);
                this.PastJacobianEstimates = cast([], dataType);
            else               
                % States:
                % * Parameters
                % * PastMeasurements
                % * PastJacobianEstimates
                % * 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_;
                % ParameterCovariance: Only used with FF, KF algorithms
                if ~this.isInactivePropertyImpl('InitialParameterCovariance')
                    this.ParameterCovariance_ = this.InitialParameterCovariance_;
                else
                    % Unused
                    this.ParameterCovariance_ = cast([],dataType);
                end         
                
                this.PastMeasurements = ...
                    zeros(1,recursiveOE.numberOfRequiredPastMeasurements(this.nb,this.nf,this.nk),dataType);
                this.PastJacobianEstimates = ...
                    zeros(1,recursiveOE.numberOfRequiredPastJacobianEstimates(this.nb,this.nf,this.nk),dataType);                
                
                % Unused
                this.PastOutputs_ = cast([],dataType);
                this.PastInputs_ = cast([],dataType);
            end

            % Reset the estimator internal states, update alg. outputs
            [this.Parameters_,this.ParameterCovariance_] = ...
                this.Estimator_.reset(...
                    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_ = recursiveOE.getAlgorithmParameters(p);
        end
        
        function s = getDiscreteStateImpl(this)
            % Present estimator states to the user
            s = getDiscreteStateImpl@idRecursivePolynomialModelEstimator(this);
            s.B = this.B;
            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','InitialC','InitialD','A','C','D'}
                    isInactive = true();
                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
            
            [~,nbm] = recursiveOE.numberOfRequiredPastJacobianEstimates(p.nb,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.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)+sum(p.nf,2));
            algorithmParams.necessaryDataForPhi = uint16(necessaryDataForPhi);
            algorithmParams.necessaryDataForPsi = uint16(necessaryDataForPsi);
            % B(q) parameters
            algorithmParams.b0PosInPhi = uint16(1);
            algorithmParams.b0PosInPsi = uint16(1);
            % F(q) parameters
            algorithmParams.f1PosInTheta = uint16(p.nb+1);
            algorithmParams.fnfPosInTheta = uint16(p.nb+p.nf);
            algorithmParams.f1PosInPhi = uint16(nbm+1);
            algorithmParams.f1PosInPsi = uint16(nbm+1);
        end
        
        function [num,nbm] = numberOfRequiredPastJacobianEstimates(nb,nf,nk)
            if nb>0
                nbm = max(nb+nk-1,nf);
            else
                nbm = 0;
            end
            num = nbm+nf;
        end
        
        function num = numberOfRequiredPastMeasurements(nb,nf,nk)
            if nb>0
                nbm = max(nb+nk-1,nf);
            else
                nbm = 0;
            end
            num = nbm+nf;
        end
    end
    
    methods(Static,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 nb nk nb nd na nc
