function D = utFoldDelay(D,InputDelay,OutputDelay)
% Folds specified subset of input and output delays into
% internal delays.

%   Author(s): P. Gahinet
%   Copyright 1986-2020 The MathWorks, Inc.
Delay = D.Delay;

% Process inputs
if isempty(InputDelay)
   InputDelay = zeros(size(Delay.Input));
end
infold = find(InputDelay>0); 
nfin = numel(infold);
if isempty(OutputDelay)
   OutputDelay = zeros(size(Delay.Output));
end
outfold = find(OutputDelay>0);
nfout = numel(outfold);
if nfin==0 && nfout==0
   return
end

% Sizes
[~,~,k,b1,b2,f1,f2,g1,g2,d11,d12,d21,d22] = getBlockData(D);
nfd = length(Delay.Internal);
[ny,nu] = size(d11);
nq = size(k,1);
m = nfd+nfin+nfout;

% Update delays
% RE: Supports one-shot specification and folding of input or output delays
Delay.Input = max(0,Delay.Input - InputDelay);
Delay.Output = max(0,Delay.Output - OutputDelay);
Delay.Internal = [...
      OutputDelay(outfold) ; ...
      Delay.Internal ; ...
      InputDelay(infold)];

% Build B,C,D
aux = [d11 ; d11(outfold,:) ; d21];
d = [ aux sparse(ny+nfout+nfd,nfout) [d12 ; d12(outfold,:) ; d22]  aux(:,infold) ; ...
   sparse(nfin,nu+m) ];
b2 = [sparse(nq,nfout) b2 b1(:,infold)];
f2 = [f1(outfold,:) ; f2 ; sparse(nfin,nq)];
g2 = [g1(outfold,:) ; g2 ; sparse(nfin,nq)];
b1(:,infold) = 0;   d(:,infold) = 0;
f1(outfold,:) = 0;  g1(outfold,:) = 0;   d(outfold,:) = 0;
Gy = speye(ny);
d(1:ny,nu+1:nu+nfout) = Gy(:,outfold);
Hu = speye(nu);  
d(ny+nfout+nfd+1:ny+m,1:nu) = Hu(infold,:);
D.b = [b1 , b2];
D.f = [f1 ; f2];
D.g = [g1 ; g2];
D.d = d;
D.Delay = Delay;