classdef (Hidden = true) AbstractSS < ltipack.ltidata
   % Class definition for @AbstractSS (dense or sparse state-space data)
   
   %   Author(s): P. Gahinet
   %   Copyright 2020 The MathWorks, Inc.

   % E dx = dx0 + A (x-x0) + B (u-u0)
   %    y = y0 + C (x-x0) + D (u-u0)
   properties
      a
      b
      c
      d
      e
      dx0
      x0
      u0
      y0
      Scaled = false;  % true if model has been scaled
   end
   
   methods
      
      function boo = isfinite(D)
         % Returns TRUE if the model data does not contain NaN's or Inf's.
         boo = allfinite(D.d);
      end
      
      function boo = isreal(D)
         % Returns TRUE if model has real data.
         boo = isreal(D.a) && isreal(D.b) && isreal(D.c) && ...
            isreal(D.d) && isreal(D.e) && isreal(D.dx0) && ...
            isreal(D.x0) && isreal(D.u0) && isreal(D.y0);
      end
      
      function boo = isstatic(D)
         % Returns TRUE if model is a pure gain.
         if isempty(D.e)
            boo = isempty(D.a) && ~hasdelay(D);
         else
            boo = nnz(D.e)==0 && ~hasdelay(D); % no differential variables
         end
      end

      function boo = hasOffset(D)
         % Returns TRUE if model has offsets.
         boo = any(D.dx0) || any(D.x0) || any(D.u0) || any(D.y0);
      end
      
      function nx = order(D)
         % Computes order of SS models
         nx = size(D.a,1);
      end
      
      function [ny,nu] = iosize(D)
         % Returns I/O size.
         %   [NY,NU] = IOSIZE(SYS)
         %   S = IOSIZE(SYS) returns S = [NY NU].
         [ny,nu] = size(D.d);  % NOTE: Data defines I/O size
         nfd = length(D.Delay.Internal);
         ny = ny-nfd;   nu = nu-nfd;
         if nargout<2
            ny = [ny nu];
         end
      end


      function D = fixOffsets(D)
         % REVISIT: Only there to mitigate incorrect size
         [ny,nu] = size(D.d);
         nx = size(D.a,1);
         if ~any(D.dx0) && numel(D.dx0)~=nx
            %fprintf('Fixing DX0 offset length\n')
            D.dx0 = zeros(nx,1);
         end
         if ~any(D.x0) && numel(D.x0)~=nx
            %fprintf('Fixing X0 offset length\n')
            D.x0 = zeros(nx,1);
         end
         if ~any(D.u0) && numel(D.u0)~=nu
            %fprintf('Fixing U0 offset length\n')
            D.u0 = zeros(nu,1);
         end
         if ~any(D.y0) && numel(D.y0)~=ny
            %fprintf('Fixing Y0 offset length\n')
            D.y0 = zeros(ny,1);
         end
      end

      function [f,g] = getAffineOffsets(D)
         % Gets affine offsets f,g such that
         % E dx = A x + B u + f
         %    y = C x + D u + g
         f = D.dx0-D.a*D.x0-D.b*D.u0;
         g = D.y0-D.c*D.x0-D.d*D.u0;
      end
            
      function D = xperm(D,perm)
         % Applies state permutation
         nx = size(D.a,1);
         if numel(perm)~=nx
            error(message('Control:transformation:xperm',nx))
         end
         D.a = D.a(perm,perm);
         D.b = D.b(perm,:);
         D.c = D.c(:,perm);
         if ~isempty(D.e)
            D.e = D.e(perm,perm);
         end
         D.dx0 = D.dx0(perm,:);
         D.x0 = D.x0(perm,:);
      end
      
      function D = ioperm(D,yperm,uperm)
         % Applies I/O permutations
         D.Delay.Input = D.Delay.Input(uperm,:);
         D.Delay.Output = D.Delay.Output(yperm,:);
         nfd = length(D.Delay.Internal);
         if nfd>0
            [ny,nu] = size(D.d);
            yperm = [yperm(:) ; (length(yperm)+1:ny).'];
            uperm = [uperm(:) ; (length(uperm)+1:nu).'];
         end
         D.b = D.b(:,uperm);
         D.c = D.c(yperm,:);
         D.d = D.d(yperm,uperm);
         D.u0 = D.u0(uperm,:);
         D.y0 = D.y0(yperm,:);
      end
      
      function D = scaleTime(D,sf)
         % Time vector scaling tnew = sf * told
         D = scaleTime@ltipack.ltidata(D,sf);
         if D.Ts==0
            % Reflect change in dx/dt
            if isempty(D.e)
               D.a = D.a/sf;   D.b = D.b/sf;   D.dx0 = D.dx0/sf;
            else
               D.e = sf*D.e;
            end
            D.Scaled = false;
         end
      end
      
      function D = uminus(D)
         % Computes -D (model whose output is -y
         nfd = length(D.Delay.Internal);
         if nfd==0
            D.c = -D.c;
            D.d = -D.d;
            D.y0 = -D.y0;
         else
            ny = size(D.d,1)-nfd;
            D.c(1:ny,:) = -D.c(1:ny,:);
            D.d(1:ny,:) = -D.d(1:ny,:);
            D.y0(1:ny,:) = -D.y0(1:ny,:);
         end
      end
      
      function D = conj(D)
         % Forms model with complex conjugate coefficients.
         D.a = conj(D.a);
         D.b = conj(D.b);
         D.c = conj(D.c);
         D.d = conj(D.d);
         D.e = conj(D.e);
         D.dx0 = conj(D.dx0);
         D.x0 = conj(D.x0);
         D.u0 = conj(D.u0);
         D.y0 = conj(D.y0);
      end
      
      function D = transpose(D)
         % Transposition of state-space models.
         % Offset interpretation: transpose of linear model with state
         % x-x0, input u-u0, output y-y0
         Ts = D.Ts;
         if (Ts==0 && any(D.dx0)) || (Ts~=0 && any(D.dx0-D.x0))
            error(message('Control:transformation:Transpose6'))
         end
         D.a = D.a.';
         aux = D.b.';  D.b = D.c.';  D.c = aux;
         aux = D.u0;  D.u0 = D.y0;  D.y0 = aux;
         D.d = D.d.';
         D.e = D.e.';
         id = D.Delay.Input;
         D.Delay.Input = D.Delay.Output;
         D.Delay.Output = id;
      end

      % Delay management -----------------------------------------
      function boo = hasdelay(D)
         % Returns TRUE if model has delays.
         Delay = D.Delay;
         boo = any(Delay.Input) || any(Delay.Output) || any(Delay.Internal);
      end
      
      function boo = hasInternalDelay(D)
         % Returns T if model has internal delays.
         boo = any(D.Delay.Internal);
      end
      
      function D = checkDelay(D)
         % Checks validity and consistency of delay data.
         [Ny,Nu] = iosize(D);
         Ts = D.Ts;
         Delay = D.Delay;
         Delay.Input = ltipack.util.checkInputDelay(Delay.Input,Nu,Ts);
         Delay.Output = ltipack.util.checkOutputDelay(Delay.Output,Ny,Ts);
         Delay.Internal = ltipack.util.checkInternalDelay(Delay.Internal,Ts);
         D.Delay = Delay;
      end
            
      function [D,tau] = getInternalDelayModel(D)
         tau = D.Delay.Internal;
         nfd = numel(tau);
         D.Delay.Internal = zeros(0,1);
         D.Delay.Input = [D.Delay.Input ; zeros(nfd,1)];
         D.Delay.Output = [D.Delay.Output ; zeros(nfd,1)];
      end
      
      function D = setInternalDelayModel(D,tau)
         % Special LFT with pure delays (used by setDelayModel)
         nfd = numel(tau);
         [rs,cs] = iosize(D);
         nu = cs-nfd;  ny = rs-nfd;
         D.Delay.Internal = [tau + D.Delay.Input(nu+1:cs,:) + ...
            D.Delay.Output(ny+1:rs,:) ; D.Delay.Internal];
         D.Delay.Input = D.Delay.Input(1:nu,:);
         D.Delay.Output = D.Delay.Output(1:ny,:);
      end
      
      function D = utFoldDelay(D,InputDelay,OutputDelay)
         % Folds specified subset of input and output delays into internal delays.
         [ny,nu] = iosize(D);
         if isempty(InputDelay)
            InputDelay = zeros(nu,1);
         end
         if isempty(OutputDelay)
            OutputDelay = zeros(ny,1);
         end
         [D.b,D.c,D.d,D.u0,D.y0,D.Delay] = ltipack.foldDelay('lti',...
            D.b,D.c,D.d,D.u0,D.y0,D.Delay,InputDelay,OutputDelay);
      end
   
   end
   
   methods (Access=protected)
      
      function D = elimZeroDelay(D)
         % Folds zero internal delays into H(s) or H(z) (rational part)
         isZero = (D.Delay.Internal==0);
         idx = find(isZero);
         ndx = length(idx);
         if ndx>0
            % Close all LFT loops associated with zero internal delays
            % NOTE: This may increase the order or produce a NaN model.
            % Push non-appx delays to last row/columns of D.d
            [ny,nu] = iosize(D);
            dperm = [idx ; find(~isZero)];
            yperm = [(1:ny).' ; ny+dperm];
            uperm = [(1:nu).' ; nu+dperm];
            D.b = D.b(:,uperm);
            D.c = D.c(yperm,:);
            D.d = D.d(yperm,uperm);
            D.u0 = D.u0(uperm,:);
            D.y0 = D.y0(yperm,:);
            D.Delay.Internal(idx,:) = [];
            D.Delay.Input(nu+ndx,1) = 0;
            D.Delay.Output(ny+ndx,1) = 0;
            D = lft(D,createGain(D,eye(ndx)),nu+1:nu+ndx,ny+1:ny+ndx,1:ndx,1:ndx);
         end
      end

      function MaxDelay = getMaxDelay(D)
         % Estimates max I/O delay.
         MaxDelay = 0;
         for ct=1:length(D)
            Delay = D(ct).Delay;
            % Note: crude estimate only, does not take into
            % account structure of internal delay model.
            MaxDelay = max(MaxDelay,...
               max([0;Delay.Input]) + max([0;Delay.Output]) + sum(Delay.Internal));
         end
      end
      
      function boo = utIsIOScaling(D)
         % True if D is of the form diag(dj*exp(-s*tauj))
         [ny,nu] = size(D.d);
         if ny==nu && isempty(D.a) && isempty(D.Delay.Internal) && ~hasOffset(D)
            [i,j] = find(D.d);
            boo = isequal(i,j);
         else
            boo = false;
         end
      end

   end

   
   methods (Static, Access=protected)
      
      function [Delay,fid,fod,ffd] = discretizeDelay(Delay,Ts)
         % Extracts discrete delays for discretization purposes.
         %
         %   [DELAYD,FID,FOD] = DISCRETIZEDELAY(DELAYC,TS) extracts the
         %   discrete input and output delays when discretizing D
         %   with sampling interval Ts.  The structure DELAYD contains
         %   the integer-valued discrete delays and the vectors
         %   FID and FOD contain the residual normalized fractional
         %   delays at the inputs and outputs (FID and FOD take value
         %   in [0,1)). The internal delays are set to zero.
         %
         %   [DDELAY,FID,FOD,FFD] = DISCRETIZEDELAY(DELAYC,TS) also
         %   discretizes the internal delays.
         
         % Add static method to be included for compiler
         %#function ltipack.util.splitDelay
         
         % Compute discrete input and output delays
         [Delay.Input,fid] = ltipack.util.splitDelay(Delay.Input,Ts);
         [Delay.Output,fod] = ltipack.util.splitDelay(Delay.Output,Ts);
         if nargout<4
            Delay.Internal(:) = 0;
         else
            [Delay.Internal,ffd] = ltipack.util.splitDelay(Delay.Internal,Ts);
         end
      end
            
   end
   
   
end


