function [D,xkeep] = sminreal(D)
% Eliminates states and delays that do not structurally contribute to 
% the I/O map. The boolean vector XKEEP indicates which states are kept 
% and which states are discarded.

%   Author(s): P. Gahinet
%   Copyright 1986-2023 The MathWorks, Inc.
nx = size(D.a,1);
[nyz,nuw] = size(D.d);
nfd = numel(D.Delay.Internal);
ny = nyz-nfd;
nu = nuw-nfd;
bStruct = spones(D.b);

if nfd==0
   % No internal delays
   % Note: Same as considering [b,dx0], but preserves nu and nu=0.
   bStruct = bStruct + spones(D.dx0);
   [D.a,~,D.c,D.e,xkeep,ekeep] = smreal(D.a,bStruct,D.c,D.e);
   ix = find(xkeep);
   ie = find(ekeep);
   D.b = D.b(ie,:);
else
   % Compute structurally minimal realization wrt both states and delays
   if isempty(D.e)
      E = [];
   else
      E = ltipack.util.blkdiag(D.e,eye(nfd));
   end
   bStruct(:,1:nu) = bStruct(:,1:nu) + spones(D.dx0);
   [~,~,~,~,xdkeep,edkeep] = smreal(...
      [D.a bStruct(:,nu+1:nuw);D.c(ny+1:nyz,:) D.d(ny+1:nyz,nu+1:nuw)],...
      [bStruct(:,1:nu);D.d(ny+1:nyz,1:nu)],[D.c(1:ny,:),D.d(1:ny,nu+1:nuw)],E);
   xkeep = xdkeep(1:nx);
   ix = find(xkeep);
   ie = find(edkeep(1:nx));
   id = find(xdkeep(nx+1:nx+nfd)');
   iuw = [1:nu nu+id];
   iyz = [1:ny ny+id];
   D.a = D.a(ie,ix);
   if ~isempty(D.e)
      D.e = D.e(ie,ix);
   end
   D.b = D.b(ie,iuw);
   D.c = D.c(iyz,ix);
   D.d = D.d(iyz,iuw);
   D.u0 = D.u0(iuw,:);
   D.y0 = D.y0(iyz,:);
   D.Delay.Internal = D.Delay.Internal(id,:);
end
D.x0 = D.x0(ix,:);
D.dx0 = D.dx0(ie,:);
