function D = elimBranching(D,I,J)
% Eliminates input branching.
% Inputs U of D 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 2020 The MathWorks, Inc.

% Enforce input delay and offset consistency along branches
nuu = numel(I);  nu = numel(J);
idu = D.Delay.Input;
u0  = D.u0;
iduu = NaN(nuu,1);
uu0 = zeros(nuu,1);
m = zeros(nuu,1);
for ct=1:nu
   j = J(ct);
   iduu(j) = min(iduu(j),idu(ct));
   uu0(j) = uu0(j)+u0(ct);  m(j) = m(j)+1;
end
uu0 = uu0./m;  % take mean value of offsets for shared inputs
idfold = idu-iduu(J);
if any(idfold)
   D = utFoldDelay(D,idfold,[]);
end
u0fold = D.u0(1:nu,:)-uu0(J);
if any(u0fold)
   D.dx0 = D.dx0 - D.b(:,1:nu)*u0fold;
   D.y0 = D.y0 - D.d(:,1:nu)*u0fold;
end

% Eliminate branches
% UU = EYE * UU -> UU(J) = EYE(J,:) * UU (J has repetitions)
d = D.d;
M = eye(nuu,'like',d);  
M = M(J,:);
D.b = [D.b(:,1:nu) * M , D.b(:,nu+1:end)];
D.d = [d(:,1:nu) * M , d(:,nu+1:end)];
D.Delay.Input = iduu;
D.u0 = [uu0 ; D.u0(nu+1:end,:)];
