classdef polynomial < pmodel.Generic
   % Linear polynomial model structure: Ay = B/F u + C/D e
   % A, C, D and F polynomials are required to be monic.
   
   %   Author(s): Rajiv Singh
   %   Copyright 2009-2024 The MathWorks, Inc.
      
   properties (Dependent)
      % A polynomial (array of row-vector parameters)
      %
      % Use this property to read the current values of the polynomials, to
      % initialize them and/or to fix specific entries. It can be set to []
      % to denote defaults: diagonal polynomials are 1 and off-diagonal
      % polynomials are 0. Diagonal polynomials (A(i,i).Value) must be
      % monic. For discrete-time systems the off-diagonal entries must
      % start with zeros. For continuous-time systems, the off-diagonal
      % polynomials must be zero.
      A
      % B polynomial (Ny-by-Nu array of row-vector parameters)
      %
      % Use this property to read the current values of the polynomials, to
      % initialize them and/or to fix specific entries. It can be set to []
      % for systems with no inputs. If B = [], F must also be [].
      B
      % C polynomial (Ny-by-1 array of row-vector parameters)
      %
      % Use this property to read the current values of the polynomials, to
      % initialize them and/or to fix specific entries. C = [] denotes the
      % default polynomial (C(s) or C(q^-1) = 1). The leading coefficient
      % of each C polynomial is fixed to 1.
      C
      % D polynomial (Ny-by-1 array of row-vector parameters)
      %
      % Use this property to read the current values of the polynomials, to
      % initialize them and/or to fix specific entries. D = [] denotes the
      % default polynomial (D(s) or D(q^-1) = 1). The leading coefficient
      % of each D polynomial is fixed to 1.
      D
      % F polynomial (Ny-by-Nu array of row-vector parameters)
      %
      % Use this property to read the current values of the polynomials, to
      % initialize them and/or to fix specific entries. It can be set to []
      % for default values: if B = [], which denotes Nu=0, F = [] denotes
      % absence of F polynomial from model; if Nu>0, f = [] denotes that
      % the default value is in use (F(s) = F(q^-1) = 1). The leading
      % coefficient of each F polynomial is fixed to 1.
      F
      % Input-output delay (Ny-by-Nu double matrix). Use IODelay to specify
      % delay in applicable time units (continuous-time) or number of
      % samples (discrete-time).
      IODelay
   end
   
   properties
      % Flag denoting presence/absence of integration on noise channels
      % (Ny-by-1 logical vector). IntegrateNoise(i) = TRUE implies that the
      % noise channel for i:th output contains an integrator. This means
      % that the D polynomial contains an additional term (which is not
      % contained in the "d" property value) equal to [1 0] for
      % continuous-time system, and equal to [1 -1] for discrete-time
      % systems.
      IntegrateNoise
   end
   
    properties (Hidden)
      % A polynomial. Alias for "A" property.
      a
      % B polynomial. Alias for "B" property.
      b
      % C polynomial. Alias for "C" property.
      c
      % D polynomial. Alias for "D" property.
      d
      % F polynomial. Alias for "F" property.
      f
      % Input-output delay. Alias for "IODelay" property.
      ioDelay      
   end
   
   methods (Static, Hidden)
      function D = array(size)
         % Create a ssdata array of a given size
         D = pmodel.polynomial.newarray(size);
      end

      function M = fixNkInB(M,Ts,nk)
         % Fix leading zeros entries of B polynomials (required by
         % IDPOLY in discrete-time case). If a polynomial is all zero, fix
         % the values to maintain backward compatibility.
         
         % If nk is provided as an input argument, only nk leading zeros
         % are fixed. Then, M.b must contain nk or more leading zeros.
         
         CustomNk = nargin>2;
         b = M.b;
         if ~isempty(b) && Ts~=0
            % fix leading zeros of B polynomial
            for ct = 1:numel(b)
               if CustomNk
                  b(ct).Free(1:nk(ct)) = false;
               else
                  Ind = find(b(ct).Value,1,'first');
                  if ~isempty(Ind)
                     b(ct).Free(1:Ind-1) = false;
                  else
                     b(ct).Free = false; % all zeros are fixed
                  end
               end
            end
            M.b = b;
         end
      end % fixNkInB
      
      function P = value2poly(Name, Value)
         % Create an array of param.Continuous objects to represent the
         % polynomial of given name and value. Value is assumed to be a
         % cell array, even in SISO case. Each element is assumed to be a
         % non-Inf double full row vector.
         Name = upper(Name);
         if isempty(Value)
            P = []; return
         elseif isnumeric(Value)
            Value = {Value};
         end
         
         [Ny, Nc] = size(Value);
         IsB = strcmp(Name,'B');
         
         P(Ny, Nc) = param.Continuous;
         for ct = 1:Ny*Nc
            % No need for checking value format/attribute since polynomial
            % constructor is invoked only by IDPOLY.
            P(ct) = param.Continuous(Name,Value{ct});
            if ~IsB
               P(ct).Free(1) = false;
            end
         end
      end
      
      function P = checkPolyParameter(P, Name, Ts)
         % Validate the value of given polynomial (array).
         
         if isequal(P,[]), return, end
         
         if ~isa(P,'param.Continuous')
            error(message('Control:pmodel:badp',Name))
         elseif any(~strcmp({P.Name},Name))
            error(message('Control:pmodel:rename',Name))
         end
         
         ny = size(P,1);
         IsA = strcmp(Name,'A');
         IsB = strcmp(Name,'B');
         
         for ct = 1:numel(P)
            Val = P(ct).Value;
            if isempty(Val)
               error(message('Control:pmodel:emptypolyval',Name))
            elseif ~isvector(Val) || any(isinf(Val))
               error(message('Control:pmodel:finitevecrequired',Name))
            else
               Val2 = double(full(Val(:).'));
               if isequal(size(Val),size(Val2))
                  P(ct).Value = Val2;
               else
                  P2 = param.Continuous(Name,Val2);
                  P2.Free = P(ct).Free.';
                  P2.Minimum = P(ct).Minimum.';
                  P2.Maximum = P(ct).Maximum.';
                  P2.Scale = P(ct).Scale.';
                  P2.Info = P(ct).Info';
                  P(ct) = P2;
               end
            end
            
            if ~IsB
               if ~IsA && Val(1)~=1
                  error(message('Control:pmodel:MonicPoly',Name))
               elseif IsA && rem(ct,ny+1)==1 && Val(1)~=1
                  % Diagonal A polynomial must be monic
                  if ny==1
                     error(message('Control:pmodel:MonicPoly',Name))
                  else
                     error(message('Control:pmodel:MonicAPoly'))
                  end
               end
               P(ct).Free(1) = false;
            end
         end
         
         if ~IsA || ny<=1, return, end
         
         % Check off-diagonal A polynomials if ny>1.
         if Ts~=0
            % Discrete-time: off-diagonal polynomials must start with 0.
            for i = 1:ny
               for j = 1:ny
                  if i~=j && P(i,j).Value(1)~=0
                     error(message('Control:pmodel:OffdiagAPolyDT'))
                  end
               end
            end
         else
            % Continuous-time: non-diagonal array not supported
            I = 1:ny^2; I(1:ny+1:end) = [];            
            for i = I
               if norm(P(i).Value)~=0
                  error(message('Control:pmodel:CTIdpolyNondiagA'))
               end
            end
            
            %{
            % Continuous-time: polyorder(P(i,j)) < polyorder(P(i,i)), j~=i
            % Future (row-reduced) form
            for i = 1:ny
               Li = length(P(i,i).Value);
               for j = 1:ny
                  if i~=j && length(P(i,j).Value)>=Li
                     error(message('Control:pmodel:OffdiagAPolyCT')
                  end
               end
            end
            %}
         end
      end % function
      
   end % static methods
   
   methods (Access = protected)
      
      % PARAMETER SERIALIZATION INTERFACE (PMODEL.GENERIC)
      function ps = getParamSet(M)
         % Get vector of parameters in scalar polynomial object M.
         
         % Keep parameters for each output are together.
         ny = iosize(M,1);
         ap = M.a; bp = M.b; cp = M.c; dp = M.d; fp = M.f;
         ps = [];
         for ct = 1:ny
            if ~isempty(ap), ps = [ps, ap(ct,:)]; end
            if ~isempty(bp), ps = [ps, bp(ct,:)]; end
            if ~isempty(cp), ps = [ps, cp(ct)];   end
            if ~isempty(dp), ps = [ps, dp(ct)];   end
            if ~isempty(fp), ps = [ps, fp(ct,:)]; end
         end
         ps = ps';
      end
      
      function M = setParamSet(M,ps)
         % Set parameter values in scalar polynomial object M.
         IOSize = size(M.ioDelay);
         Ny = IOSize(1); Nu = IOSize(2);
         Offset = 0;
         for ky = 1:Ny
            if ~isempty(M.a)
               M.a(ky,:) = ps(Offset+1:Offset+Ny)';
               Offset = Offset + Ny;
            end
            
            if ~isempty(M.b)
               M.b(ky,:) = ps(Offset+1:Offset+Nu)';
               Offset = Offset + Nu;
            end
            
            if ~isempty(M.c)
               M.c(ky) = ps(Offset+1);
               Offset = Offset + 1;
            end
            
            if ~isempty(M.d)
               M.d(ky) = ps(Offset+1);
               Offset = Offset + 1;
            end
            
            if ~isempty(M.f)
               M.f(ky,:) = ps(Offset+1:Offset+Nu)';
               Offset = Offset + Nu;
            end
         end
      end
      
      function dispScalarModel(M)
         % Display of a scalar ss object.
         ff = 'abcdf';
         for ct = 1:length(ff)
            val = M.(ff(ct));
            if ~isequal(val,[])
               S.(upper(ff(ct))) = val;
            end
         end
         S.IODelay = M.ioDelay;
         S.IntegrateNoise = M.IntegrateNoise;
         disp(S)
      end
      
      function dispModelType(M,Size)
         % Display model type information.
         if prod(Size)==1
            Type = getPolynomialType(M);
            if Type(1)=='G'
               Type = 'Polynomial';
            end
            
            disp(ctrlMsgUtils.message('Control:pmodel:msgPolyDisp',Type))
         else
            disp(ctrlMsgUtils.message('Control:pmodel:msgPolyDispArray',mat2str(Size)))
         end
      end
      
   end % protected methods
   
   methods
      % Constructor
      function M = polynomial(A,B,C,D,F,ioDelay)
         % Construct polynomial object based on specified polynomials.
         % Use [] to denote a non-existent polynomial in the model
         % structure. For active polynomials, use cell arrays of row vectors.
         %
         % Polynomials are retained as-is; no automatic truncation of
         % trailing zeros is performed. Polynomials C, D and F must be
         % monic. Diagonal polynomials in A must be monic while the
         % off-diagonal ones must start with zero. No dimension consistency
         % checks are performed. The input arguments must use prequalified,
         % expected values.
         
         if nargin==0
            return
         end
         
         if ~isempty(A), M.a = pmodel.polynomial.value2poly('A',A); end
         if ~isempty(B), M.b = pmodel.polynomial.value2poly('B',B); end
         if ~isempty(C), M.c = pmodel.polynomial.value2poly('C',C); end
         if ~isempty(D), M.d = pmodel.polynomial.value2poly('D',D); end
         if ~isempty(F), M.f = pmodel.polynomial.value2poly('F',F); end
         
         M.ioDelay = ioDelay;
         M.IntegrateNoise = false(size(ioDelay,1),1);
      end
      
      function Value = get.A(M)
         %GET method for "A" property
         Value = M.a;
      end
      
      function M = set.A(M, Value)
         %SET method for "A" property
         M.a = Value;
      end
      
      function Value = get.B(M)
         %GET method for "B" property
         Value = M.b;
      end
      
      function M = set.B(M, Value)
         %SET method for "B" property
         M.b = Value;
      end
      
      function Value = get.C(M)
         %GET method for "C" property
         Value = M.c;
      end
      
      function M = set.C(M, Value)
         %SET method for "C" property
         M.c = Value;
      end
      
      function Value = get.D(M)
         %GET method for "D" property
         Value = M.d;
      end
      
      function M = set.D(M, Value)
         %SET method for "D" property
         M.d = Value;
      end
      
      function Value = get.F(M)
         %GET method for "F" property
         Value = M.f;
      end
      
      function M = set.F(M, Value)
         %SET method for "F" property
         M.f = Value;
      end
      
      function Value = get.IODelay(M)
         %GET method for "IODelay" property
         Value = M.ioDelay;
      end
      
      function M = set.IODelay(M, Value)
         %SET method for "IODelay" property
         M.ioDelay = Value;
      end
      
      function [M, NVScale, K] = scaleTime(M,sf,Ts)
         % Time vector scaling tnew = sf * told
         K = [];
         if Ts==0
            iod = M.ioDelay;
            [ny, nu] = size(iod); NVScale = zeros(1,ny);
            M.ioDelay =  M.ioDelay * sf;
            ap = M.a; bp = M.b; cp = M.c; dp = M.d; fp = M.f;
            
            for ky = 1:ny
               na = 0;
               if ~isempty(ap)
                  na = length(ap(ky,ky).Value)-1;
                  powA = cumprod([1 sf(:,ones(na,1))]);
                  % RE: A is diagonal for Ts==0
                  sc = 1./powA;
                  ap(ky,ky) = pmodel.scaleParameter(ap(ky,ky),sc);
                  K = [K, sc(2:end)];
               end
               
               Kf = [];
               for ku = 1:nu
                  nf = 0;
                  nb = length(bp(ky,ku).Value); 
                  if ~isempty(fp), nf = length(fp(ky,ku).Value)-1; end
                  powBF = cumprod([1 sf(:,ones(max(nb,nf),1))]);
                  if ~isempty(fp)
                     sc = 1./powBF(1:nf+1);
                     fp(ky,ku) = pmodel.scaleParameter(fp(ky,ku),sc);
                     Kf = [Kf, sc(2:end)];
                  end
                  Gain = sf^(nb-nf-na-1);
                  sc = Gain./powBF(1:nb);
                  bp(ky,ku) = pmodel.scaleParameter(bp(ky,ku),sc);
                  K = [K, sc];
               end
               
               nc = 0; nd = 0;
               if ~isempty(cp)
                  nc = length(cp(ky).Value)-1;
                  sc = 1./cumprod([1 sf(:,ones(nc,1))]);
                  cp(ky) = pmodel.scaleParameter(cp(ky),sc);
                  K = [K, sc(2:end)];
               end
               if ~isempty(dp)
                  nd = length(dp(ky).Value)-1;
                  sc = 1./cumprod([1 sf(:,ones(nd,1))]);
                  dp(ky) = pmodel.scaleParameter(dp(ky),sc);
                  K = [K, sc(2:end)];
               end
               
               if ~isempty(Kf)
                  K = [K, Kf];
               end
               
               NVScale(ky) = sf^(nc-nd-na);               
            end
            NVScale = diag(NVScale.^2);
            M.a = ap; M.b = bp; M.c = cp; M.d = dp; M.f = fp;
         else
            NVScale = 1/sf;
         end
      end
      
   end
   
   methods (Hidden)
      function Type = getPolynomialType(M)
         % Determine a type and name for the polynomial structure. Types
         % are one of: ar, ari, arma, arima, arx, arma, arix, arimax,
         % bj, oe, ararx, ararmax, fir, generic.
         [na, nb, nc, nd, nf] = polyorder(M);
         
         IsActive = ~[all(na(:)<=1), all(nb(:)==0), all(nc<=1),...
            all(nd(:)<=1), all(nf(:)<=1)];
         IsActive(6) = any(M.IntegrateNoise);
         
         if isequal(IsActive,[0 1 0 0 1 0]) || isequal(IsActive,[0 0 0 0 1 0]) 
            % y = B/F u + e
            Type = 'OE';
         elseif isequal(IsActive,[0 1 1 1 1 0])
            % y = B/F u + C/D e
            Type = 'BJ';
         elseif isequal(IsActive,[1 1 0 0 0 0])
            % Ay = Bu + e
            Type = 'ARX';
         elseif isequal(IsActive,[1 1 1 0 0 0])
            % Ay = Bu + Ce;
            Type = 'ARMAX';
         elseif isequal(IsActive,[0 1 0 0 0 0])
            % y = Bu + e
            Type = 'FIR';
         elseif isequal(IsActive,[1 1 0 1 0 0])
            % Ay = Bu + 1/D e
            Type = 'ARARX';
         elseif isequal(IsActive,[1 1 1 1 0 0])
            % Ay = Bu + C/D e
            Type = 'ARARMAX';
         elseif isequal(IsActive,[1 0 0 0 0 0]) || isequal(IsActive,[0 0 0 0 0 0])
            % Ay = e
            Type = 'AR';
         elseif isequal(IsActive,[1 0 1 0 0 0]) 
            % Ay = Ce
            Type = 'ARMA';
         elseif isequal(IsActive, [0 0 1 0 0 0])
            % y = Ce
            Type = 'MA';
         elseif isequal(IsActive,[1 0 0 0 0 1])
            % A(1-z^-1)y = e
            Type = 'ARI';
         elseif isequal(IsActive,[1 0 1 0 0 1])
            % A(1-z^-1)y = Ce
            Type = 'ARIMA';
         elseif isequal(IsActive,[1 1 0 0 0 1])
            % A(1-z^-1)y = B(1-z^-1)u + e
            Type = 'ARIX';
         elseif isequal(IsActive,[1 1 1 0 0 1])
            % A(1-z^-1)y = B(1-z^-1)u + Ce
            Type = 'ARIMAX';
         elseif ~IsActive(6)
            % Ay = B/F u + C/D e
            Type = 'Polynomial';
         else % isActive(6)
            % A(1-z^-1)y = B/F u + C/D e
            Type = 'Polynomial with Noise Integration';
         end
      end % getPolynomialType
      
      %------------------------------------------------------------------
      function varargout = polyorder(M, PolyName)
         % Determine polynomial sizes. These are raw sizes that do not
         % compensate for fixed leading entries.
         %  [na, nb, nc, nd, nf] = polyorder(M)
         %  nx = polyorder(M, PolyX)
         
         if nargin==1
            varargout{1} = localGetPolySize(M.a);
            if nargout<2, return, end
            
            varargout{2} = localGetPolySize(M.b);
            if nargout<3, return, end
            
            varargout{3} = localGetPolySize(M.c);
            if nargout<4, return, end
            
            varargout{4} = localGetPolySize(M.d);
            if nargout<5, return, end
            
            varargout{5} = localGetPolySize(M.f);
         else
            varargout{1} = localGetPolySize(M.(PolyName));
         end
         
      end % order
      
      %------------------------------------------------------------------
      function s = iosize(M,dim)
         % I/O sizes
         s = size(M.ioDelay);
         if nargin>1
            if dim<3
               s = s(dim);
            else
               s = 1;
            end
         end
      end % iosize
      
      %--------------------------------------------------------------------
      function M = iocat(dim, M1, M2)
         % Concatenate models along input (2) or output (1) dimension.
         
         if dim==1
            % Output concatenation.
            % Input sizes of M1 and M2 are assumed to match.
            M = M1;
            IOSize1 = iosize(M1);
            IOSize2 = iosize(M2);
            Ny = IOSize1(1)+IOSize2(1);
            Nu = IOSize1(2);
            
            
            % block diagonal concatenation
            if ~isempty(M1.a) || ~isempty(M2.a)
               ap = pmodel.polynomial.value2poly('A',num2cell(eye(Ny)));
               if ~isempty(M1.a)
                  ap(1:IOSize1(1),1:IOSize1(1)) = M1.a;
               end
               
               if ~isempty(M2.a)
                  ap(IOSize1(1)+1:Ny,IOSize1(1)+1:Ny) = M2.a;
               end
               M.a = ap;
            end
            
            if ~isempty(M1.b)
               % Both models have active b polynomials (b is [] only when
               % Nu = 0).
               M.b = [M1.b; M2.b];
            end
            
            if ~isempty(M1.c) || ~isempty(M2.c)
               cp = pmodel.polynomial.value2poly('C',num2cell(ones(Ny,1)));
               if ~isempty(M1.c)
                  cp(1:IOSize1(1)) = M1.c;
               end
               
               if ~isempty(M2.c)
                  cp(IOSize1(1)+1:Ny) = M2.c;
               end
               M.c = cp;
            end
            
            if ~isempty(M1.d) || ~isempty(M2.d)
               dp = pmodel.polynomial.value2poly('D',num2cell(ones(Ny,1)));
               if ~isempty(M1.d)
                  dp(1:IOSize1(1)) = M1.d;
               end
               
               if ~isempty(M2.d)
                  dp(IOSize1(1)+1:Ny) = M2.d;
               end
               M.d = dp;
            end
            
            if ~isempty(M1.f) || ~isempty(M2.f)
               fp = pmodel.polynomial.value2poly('F',num2cell(ones(Ny,Nu)));
               if ~isempty(M1.f)
                  fp(1:IOSize1(1),:) = M1.f;
               end
               
               if ~isempty(M2.f)
                  fp(IOSize1(1)+1:Ny,:) = M2.f;
               end
               M.f = fp;
            end
            
            M.ioDelay = [M1.ioDelay; M2.ioDelay];
            M.IntegrateNoise = [M1.IntegrateNoise; M2.IntegrateNoise];
            
         else
            % Input Concatenation.
            % Requires models with no noise dynamics (A=C=D=1);
            % Output sizes of M1 and M2 are assumed to match.
            
            [na, nb, nc, nd, nf] = polyorder(M1);
            I1 = ~[all(na(:)==0), all(nb(:)==0), all(nc==0),...
               all(nd(:)==0), all(nf(:)==0)];
            I1(6) = any(M1.IntegrateNoise);
            
            [na, nb, nc, nd, nf] = polyorder(M2);
            I2 = ~[all(na(:)==0), all(nb(:)==0), all(nc==0),...
               all(nd(:)==0), all(nf(:)==0)];
            I2(6) = any(M2.IntegrateNoise);
            
            Valid = isequal(I1,[0 1 0 0 1 0]) || isequal(I1,[0 0 0 0 1 0])...
               || isequal(I1,[0 1 0 0 0 0]) || isequal(I1,[0 0 0 0 0 0])...
               || isequal(I2,[0 1 0 0 1 0]) || isequal(I2,[0 0 0 0 1 0])...
               || isequal(I2,[0 1 0 0 0 0]) || isequal(I2,[0 0 0 0 0 0]);
            if ~Valid
               % Valid for OE, FIR and trivial AR forms
               error(message('Control:pmodel:polynomialHorzcat'))
            end            
            M = M1;
            M.b = [M.b, M2.b];
            f1 = M.f; f2 = M2.f; 
            I = [isequal(f1,[]),isequal(f2,[])];
            if xor(I(1),I(2))
               if I(1)
                  ios = iosize(M);
                  f0 = pmodel.polynomial.value2poly('f',num2cell(ones(ios)));
                  fp = [f0, f2];
               else
                  ios = iosize(M2);
                  f0 = pmodel.polynomial.value2poly('f',num2cell(ones(ios)));
                  fp = [f1, f0];
               end
            else
               fp = [f1, f2];
            end
            M.f = fp;
            M.ioDelay = [M.ioDelay, M2.ioDelay];
         end
      end
      
      %--------------------------------------------------------------------
      function M = setIODelayValue(M, Value)
         % Update IO delay value
         M.ioDelay = Value;
      end
      
      %--------------------------------------------------------------------
      function M = getsubsys(M,rowIndex,colIndex)
         % Extracts polynomial model for given I/O pair.
         
         if isempty(rowIndex)
            M.a = []; M.b = []; M.c = []; M.d = []; M.f = [];
         else
            if ~isequal(M.a,[])
               Ny = size(M.a,1);
               if Ny>1 && ~strcmp(rowIndex,':') && ~isequal(rowIndex,1:Ny)
                  av = reshape({M.a.Value},[Ny Ny]);
                  Separable = all(all(cellfun(@(x)isempty(x) || isequal(x,0),...
                     av(rowIndex,setdiff(1:Ny,rowIndex)))));
                  if ~Separable
                     error(message('Control:pmodel:polynomialGetsubsys1'))
                  end
               end
               
               ru = unique(rowIndex);
               if length(ru)<length(rowIndex)   
                  [r, I] = sort(rowIndex);
                  Ny2 = numel(r);
                  ap(Ny2,Ny2) = param.Continuous('A',0);  
                  ap = pmodel.changeParName(ap,'A');
                  L = 0;
                  for ct = 1:numel(ru)
                     rct = ru(ct);
                     N = sum(r==rct);
                     thisa = M.a(rct,:); % row of A to be copied
                     L2 = L*(Ny2+1);
                     ap(L2+1:Ny2+1:L2+N+(Ny2*(N-1))) = thisa(rct);
                     if rct>1
                        ap(L+1:L+N,1:rct-1) = repmat(thisa(1:rct-1),[N 1]);
                     end
                     if rct<Ny
                        ap(L+1:L+N,rct+N-1+(1:Ny-rct)) = repmat(thisa(rct+1:end),[N 1]);
                     end
                     L = L+N;
                  end
                  
                  if ~isequal(I(:)', 1:Ny2)
                     % Restore requested order of output channels 
                     ap = ap(I,I);
                  end
                                    
                  M.a = ap;
               else
                  % no repetitions
                  M.a = M.a(rowIndex, rowIndex);
               end
            end
            
            if ~isequal(M.c,[]), M.c = M.c(rowIndex); end
            if ~isequal(M.d,[]), M.d = M.d(rowIndex); end
            if isempty(colIndex)
               M.b = []; M.f = [];
            else
               if ~isequal(M.b,[]), M.b = M.b(rowIndex,colIndex); end
               if ~isequal(M.f,[]), M.f = M.f(rowIndex,colIndex); end
            end
         end
         
         M.ioDelay = M.ioDelay(rowIndex,colIndex);
         M.IntegrateNoise = M.IntegrateNoise(rowIndex,1);
      end
      
      %--------------------------------------------------------------------
      function M = checkDataType(M, Ts)
         % Validates data type of A,B,C,D,F parameters and ioDelay,
         % IntegrateNoise values.
         
         M.a = pmodel.polynomial.checkPolyParameter(M.a, 'A', Ts);
         M.b = pmodel.polynomial.checkPolyParameter(M.b, 'B', Ts);
         M.c = pmodel.polynomial.checkPolyParameter(M.c, 'C', Ts);
         M.d = pmodel.polynomial.checkPolyParameter(M.d, 'D', Ts);
         M.f = pmodel.polynomial.checkPolyParameter(M.f, 'F', Ts);
         
         NI = M.IntegrateNoise;
         if ~islogical(NI) && ~isequal(NI, logical(NI))
            error(message('Control:pmodel:badnoiseintg'))
         else
            M.IntegrateNoise = logical(NI);
         end
         
         M.ioDelay = pmodel.checkIODelay(M.ioDelay,Ts);
      end

      %--------------------------------------------------------------------
      function M = utRemoveZeroPad(M, Ts)
         % Remove extra zeros from numerator and denominator.
         % If Ts==0, leading zeros from numerator are removed.
         % If Ts~=0, trailing zeros from numerator and denominator are
         % removed. Note that in the discrete-time case, the numerator and
         % the denominator polynomials are expressed in z^-1, so that the
         % trailing zeros are superfluous.
         Ap = M.a; Bp = M.b; Cp = M.c; Dp = M.d; Fp = M.f;
         if Ts==0
            % Remove leading zeros from the B polynomial (others are monic).
            for ct = 1:numel(Bp)
                v = Bp(ct).Value;
                len = length(v);
               Ind = find(v,1,'first');
               if isempty(Ind)
                  % If numerator is all zero, leave only one.
                  Ind = len;
               end
               if Ind>1
                  Bp(ct) = subsrefParameter(Bp(ct),Ind:len);
               end
            end
            M.b = Bp;
         else
            % Remove trailing zeros from all polynomials
            M.a = localRemoveTrailingZeros(Ap);
            M.b = localRemoveTrailingZeros(Bp);
            M.c = localRemoveTrailingZeros(Cp);
            M.d = localRemoveTrailingZeros(Dp);
            M.f = localRemoveTrailingZeros(Fp);
         end
      end
      
      %--------------------------------------------------------------------
      function M = checkDataSize(M)
         % Check data consistency. Note that size inconsistency with IO
         % delay causes IO delay to be updated.
         Ap = M.a; Bp = M.b; Cp = M.c; Dp = M.d; Fp = M.f;
         
         Ny = max([size(Ap,1), size(Bp,1), size(Cp,1), size(Dp,1), size(Fp,1)]);
         Nu = max([size(Bp,2), size(Fp,2)]);
         
         if ~((isequal(Ap,[]) || isequal(size(Ap),[Ny Ny])) && ...
               (isequal(Bp,[]) || isequal(size(Bp),[Ny Nu])) &&...
               (isequal(Cp,[]) || isequal(size(Cp),[Ny 1])) && ...
               (isequal(Dp,[]) || isequal(size(Dp),[Ny 1])) && ...
               (isequal(Fp,[]) || isequal(size(Fp),[Ny Nu])))
            error(message('Control:pmodel:polybadsize'))
         elseif ~isempty(Fp) && isempty(Bp)
            % A time series model being turned into an IO model via
            % modification of F (if via modification of B, no special
            % handling required since F = [] implies a reinterpreted
            % default value)
            M.b = pmodel.polynomial.value2poly('b',ones([Ny, Nu]));
         end
         
         % Update I/O delay
         OldSize = iosize(M);
         if ~isequal([Ny Nu],OldSize)
            if ~any(M.ioDelay(:))
               M.ioDelay = zeros(Ny,Nu);
            elseif (Ny~=1 || Nu~=1) && isscalar(M.ioDelay)
               M.ioDelay = M.ioDelay(ones(Ny,1),ones(Nu,1)); % silent expansion
            else
               error(message('Control:pmodel:badiodelaysize'))
            end
         end
         
         % Reset IntegrateNoise if at default value
         NI = M.IntegrateNoise;
         if ~isvector(NI)
            error(message('Control:pmodel:badnoiseintg'))
         elseif length(NI)~=Ny
            if ~any(NI)
               % silent update
               M.IntegrateNoise = false(Ny,1);
            else
               error(message('Control:pmodel:badnoiseintg'))
            end
         end
      end
   end % methods
   
end %class

%--------------------------------------------------------------------------
%                   Local Functions
%--------------------------------------------------------------------------
function N = localGetPolySize(P)
% Compute polynomial lengths contained in polynomial array P
Sz = size(P);
N = zeros(Sz);
for ct = 1:numel(N)
   N(ct) = prod(getSize(P(ct))); %numel(P(ct).Value);
end
end

%---------------------------------------------------------------------------
function P = localRemoveTrailingZeros(P)
% Remove trailing zeros from a parameter array.

for ct = 1:numel(P)
   v = P(ct).Value;
   len = length(v);
   Ind = find(v,1,'last');
   if isempty(Ind)
      % If polynomial is all zero (can happen for "B"), leave only one.
      Ind = min(1,len);
   end
   if Ind<len
      P(ct) = subsrefParameter(P(ct),1:Ind);
   end
end
end
