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

%   Copyright 2020-2023 The MathWorks, Inc.
nq = size(D.k,1);
nfd = numel(D.Delay.Internal);
hasM = ~isempty(D.m);

if nfd==0
   % No internal delays
   MC = spones(D.c);
   if hasM
      MC = MC + spones(D.m);
   end
   FG = spones(D.f)+spones(D.g);
   [D.k,D.b,~,~,qkeep,ekeep] = smreal(D.k,D.b,FG,MC);
   jx = find(qkeep);
   if hasM
      D.m = D.m(ekeep,jx);
   end
   D.c = D.c(ekeep,jx);
   D.f = D.f(:,jx);
   D.g = D.g(:,jx);
else
   % Compute structurally minimal realization wrt both states and delays
   [m,c,k,b1,b2,f1,f2,g1,g2,d11,d12,d21,d22] = getBlockData(D);
   MC = blkdiag(spones(m)+spones(c),speye(nfd));
   FG1 = spones(f1)+spones(g1);
   FG2 = spones(f2)+spones(g2);
   [~,~,~,~,qdkeep,ekeep] = smreal([k b2;FG2 d22],[b1;d21],[FG1,d12],MC);
   qkeep = qdkeep(1:nq);
   ekeep = ekeep(1:nq);
   dkeep = qdkeep(nq+1:nq+nfd);
   if hasM
      D.m = D.m(ekeep,qkeep);
   end
   D.c = D.c(ekeep,qkeep);
   D.k = D.k(ekeep,qkeep);
   D.b = [b1(ekeep,:) b2(ekeep,dkeep)];
   D.f = [f1(:,qkeep) ; f2(dkeep,qkeep)];
   D.g = [g1(:,qkeep) ; g2(dkeep,qkeep)];
   D.d = [d11 d12(:,dkeep) ; d21(dkeep,:) d22(dkeep,dkeep)];
   D.Delay.Internal = D.Delay.Internal(dkeep,:);
end

% Update StateInfo
D.StateInfo = ltipack.util.selectStateInfo(D.StateInfo,qkeep);

% Eliminate pole/zero cancellation in
%   m q'' + c q' = b u , y = g q' + d u
% by rewriting it as
%   m v' + c v = b u , y = g v + d u  (v = q')
D = downshift(D);
