function D = simplifyDelay(D)
% Replaces internal delays by input or output delays when possible.

%   Author(s): P. Gahinet
%   Copyright 1986-2020 The MathWorks, Inc.
Delay = D.Delay;
nfd = length(Delay.Internal);
if nfd>0
   % Data and size info
   K = D.k;  d = D.d;  g = D.g;
   MC = spones(D.c);
   if ~isempty(D.m)
      MC = MC + spones(D.m);
   end
   FG = spones(D.f)+spones(g);
   [nr,nc] = size(d);
   nu = nc-nfd;
   ny = nr-nfd;
   nuy = nu+ny;
   
   % Label transfers from [u;w] to [y;z] as zero (0), static (1), or
   % dynamic (2)
   [qkeep,ekeep] = iosmreal(K,D.b,FG,MC);
   IOMAP = zeros(nr,nc);
   for j=1:nc
      for i=1:nr
         qx = find(qkeep(:,i,j));
         if ~(d(i,j)==0 && isempty(qx))
            if nnz(MC(ekeep(:,i,j),qx))>0 || nnz(g(i,qx))>0
               IOMAP(i,j) = 2;
            else
               IOMAP(i,j) = 1;
            end
         end
      end
   end
   
   % Build directed graph with vertices {u,y,w,z} and edge values
   % 1 (static), 2 (dynamic), or 3 (delay)
   DG = eye(nu+ny+2*nfd);
   DG([nu+1:nuy,nuy+nfd+1:end],[1:nu,nuy+1:nuy+nfd]) = IOMAP;
   DG(nuy+1:nuy+nfd,nuy+nfd+1:end) = 3*eye(nfd);
   
   % Detect and delete edges involved in feedback loops (note: any cycle
   % must involve at least one delay)
   [p,q,r,c] = dmperm(DG);
   for ct=1:numel(c)-1
      bs = r(ct+1)-r(ct); % size of strongly connected component
      if bs>1
         DG(p(r(ct):r(ct+1)-1),q(c(ct):c(ct+1)-1)) = 0;
      end
   end
   DGO = digraph(logical(DG'));
   
   % A delay is "movable" to the input if no dynamic edge feeds into its
   % z input. Input delays are singleton movable successors of u in the
   % w subset.
   MOVABLE = ~any(DG(nuy+nfd+1:end,:)==2,2);
   MOVED = false(nfd,1);  % keeps track of internal delays moved to inputs
   for j=1:nu
      v = j;  % u(j)
      vnext = successors(DGO,v);  v = vnext(vnext~=v);
      while isscalar(v) && v>nuy
         k = v-nuy;  % index into w
         if k<=nfd && MOVABLE(k)
            Delay.Input(j) = Delay.Input(j) + Delay.Internal(k);
            MOVED(k) = true;
         end
         vnext = successors(DGO,v);  v = vnext(vnext~=v);
      end
   end
   Delay.Internal(MOVED) = 0;
   
   % A delay is "movable" to the output if its w output does not feed
   % into any dynamic edge. Output delays are singleton movable
   % predecessors of y in the z subset.
   MOVABLE = ~any(DG(:,nuy+1:nuy+nfd,:)==2,1);
   MOVED = false(nfd,1);  % keeps track of internal delays moved to inputs
   for i=1:ny
      v = nu+i;  % y(i)
      vprev = predecessors(DGO,v);  v = vprev(vprev~=v);
      while isscalar(v) && v>nuy
         k = v-(nuy+nfd);  % index into z
         if k>0 && MOVABLE(k)
            Delay.Output(i) = Delay.Output(i) + Delay.Internal(k);
            MOVED(k) = true;
         end
         vprev = predecessors(DGO,v);  v = vprev(vprev~=v);
      end
   end
   Delay.Internal(MOVED) = 0;
   
   % Eliminate zero internal delays
   if any(Delay.Internal==0)
      D.Delay = Delay;
      D = elimZeroDelay(D);
   end
end