function [A,B,C,D,E,dx0,x0,uu0,y0,Delays] = elimBranchingDF(sys,I,J,varargin)
% Eliminates input branching.
% Inputs U of SYS are assumed to be branches U = UU(J) of the unique
% set of input signals UU = U(I). Returns an equivalent model
% whose inputs are the unique signals UU.

%   Copyright 2022-2023 The MathWorks, Inc.
[A,B,C,D,E,dx0,x0,u0,y0,Delays] = sys.DataFunction(varargin{:});
nu = numel(J);   % number of branched inputs

% UU = EYE * UU -> U = UU(J) = EYE(J,:) * UU = M * UU (J has repetitions).
M = eye(numel(I));
M = M(J,:);

% Enforce input delay consistency along branches.
if ~isempty(Delays)
   % When a branched input has a varying delay on one branch, the varying
   % delays on all branches must be folded.
   VARYING = ~isnan(Delays.Input);
   nr = sum(M,1)';  % branching count
   FOLD = VARYING & nr(J)>1;
   if any(FOLD)
      idf = NaN(nu,1);
      idf(FOLD) = Delays.Input(FOLD);
      [B,C,D,u0,y0,Delays] = ltipack.foldDelay('ltv',B,C,D,u0,y0,Delays,idf,NaN(size(Delays.Output)));
      Delays.Input(FOLD) = NaN;
   end
   Delays.Input = Delays.Input(I,:);
end

% Equalize input offsets along branches.
% Note: Use offset averaging for consistency with PARALLEL
if isempty(u0)
   uu0 = [];
else
   % Using M'*M = diag(sum(M,1)), we get M'*U0 = diag(sum(M,1)) * UU0
   nr = sum(M,1)';  % counts repetitions
   uu0 = [nr.\(M'*u0(1:nu,:)) ; u0(nu+1:end,:)];
   delta = uu0(J,:)-u0(1:nu,:); % to be folded
   if isempty(dx0)
      dx0 = B(:,1:nu)*delta;
   else
      dx0 = dx0+B(:,1:nu)*delta;
   end
   if isempty(y0)
      y0 = D(:,1:nu)*delta;
   else
      y0 = y0+D(:,1:nu)*delta;
   end
end

% Consolidate B,D
B = [B(:,1:nu) * M , B(:,nu+1:end)];
D = [D(:,1:nu) * M , D(:,nu+1:end)];
