function BlocksMarked = minjacobian(J)
% MINJACOBIAN Using the information from the block reduction from Simulink,
% continue to reduce the size of the model by using the actual
% linearization to deal with blocked paths.

%  Author(s): John Glass
%  Copyright 1986-2004 The MathWorks, Inc.
% $Revision: 1.1.6.1 $ $Date: 2004/08/01 00:12:20 $

%% Get the block dimensions
[ny,nu] = size(J.D);
nxz = size(J.A,1);

%% Get the block handles
BlockHandles = J.Mi.BlockHandles;

%% Get the blocks that have been marked as in the path by the internal
%% algorithm
BlocksInPath_Pass1 = J.Mi.ForwardMark & J.Mi.BackwardMark;
% BlocksInPath_Pass1 = ones(size(BlocksInPath_Pass1));
%% Extract index elements, tack on the number of states, inputs, and
%% outputs to account for the last block in the list.  Convert from
%% C-indexing to Matlab-indexing.
InputIdx = [J.Mi.InputIdx+1;nu+1];
OutputIdx = [J.Mi.OutputIdx+1;ny+1];

%% Construct the input to block list
InputList = zeros(size(J.Mi.E,1),1);
ctr = 1;
for ct = 1:(length(InputIdx)-1)
    if (InputIdx(ct) ~= InputIdx(ct+1)) 
        lgnth = InputIdx(ct+1)-InputIdx(ct);
        InputList(ctr:ctr+lgnth-1) = ...
            BlockHandles(ct)*ones(lgnth,1)*BlocksInPath_Pass1(ct);
        ctr = ctr + lgnth;
    end
end

%% Construct the output to block list
OutputList = zeros(size(J.Mi.E,2),1);
ctr = 1;
for ct = 1:(length(OutputIdx)-1)
    if (OutputIdx(ct) ~= OutputIdx(ct+1))
        lgnth = OutputIdx(ct+1)-OutputIdx(ct);
        OutputList(ctr:ctr+lgnth-1) = ...
            BlockHandles(ct)*ones(lgnth,1)*BlocksInPath_Pass1(ct);
        ctr = ctr + lgnth;
    end
end

%% Eliminate blocks that have been removed
BlocksRemovedu = find(InputList);
BlocksRemovedy = find(OutputList);
InputList = InputList(BlocksRemovedu);
OutputList = OutputList(BlocksRemovedy);

B = J.B;
D = J.D;
E = any(J.Mi.E,3);
F = J.Mi.F;
G = J.Mi.G;
B = B(:,BlocksRemovedu);
D = D(:,BlocksRemovedu);
D = D(BlocksRemovedy,:);
E = E(BlocksRemovedu,:);
E = E(:,BlocksRemovedy);
F = F(BlocksRemovedu,:);
G = G(:,BlocksRemovedy);

%% Create the matrix that maps blocks with states and nonzero inputs to
%% outputs
[ny,nu] = size(D);
ThetaCB = sparse(ny,nu);

[row,indB] = find(B);
for ct = 1:length(indB);
    Block = InputList(indB(ct));
    if Block
        OutInd = find(Block == OutputList);
        ThetaCB(OutInd,indB(ct)) = 1;
    end
end    
Theta = any(D,3) + ThetaCB;

%% Forward loop
%% Absorb Theta into the F and G matrices
Bbc = any(F,2);
Bbo = any(Theta*F,2);
Cbc = any(G*Theta,1);
Cbo = any(G,1);

%% Compute structurally minimal realization
%% First identify structurally input signals
xcc = Bbc;    % x(1)
dx = xcc;    % dx(k) = x(k)-x(k-1)
while any(dx),
   Thetadx = any(Theta(:,dx),2);  
   Edx = any(E(:,Thetadx),2);
   dx = Edx & ~xcc;         % dx(k+1)
   xcc = xcc | Edx;          % xc(k+1)
end

%% Identify structurally observable signals
xoc = Cbc;   
dx = xoc;
while any(dx),
   Adx = any(E(dx,:),1);
   Thetadx = any(Theta(Adx,:),1);
   dx = Thetadx & ~xoc;
   xoc = xoc | Thetadx;
end

%% Minimal states are both controllable and observable
InputMask = xoc & xcc';
InputBlocksHit = unique(InputList(find(InputMask)));

%% Compute structurally minimal realization
%% Second identify structurally output signals
xco = Bbo;    % x(1)
dx = xco;    % dx(k) = x(k)-x(k-1)
while any(dx),
    Edx = any(E(:,dx),2);
    Thetadx = any(Theta(:,Edx),2);
    dx = Thetadx & ~xco;         % dx(k+1)
    xco = xco | Thetadx;          % xc(k+1)
end
  
%% Identify structurally observable signals
xoo = Cbo;   
dx = xoo;
while any(dx),
    Thetadx = any(Theta(dx,:),1);
    Edx = any(E(Thetadx,:),1);
    dx = Edx & ~xoo;
    xoo = xoo | Edx;
end

%% Minimal states are both controllable and observable
OutputMask = xoo & xco';
OutputBlocksHit = unique(OutputList(find(OutputMask)));

%% Loop over each block to determine if it is in the list
BlocksMarked = zeros(length(BlockHandles),1);

%% Get the indices of blocks in the first pass
ind_Pass1 = find(BlocksInPath_Pass1);

%% Identify the blocks that are in the linearization path
for ct = 1:length(BlockHandles(ind_Pass1))
    if (any(find(BlockHandles(ind_Pass1(ct)) == InputBlocksHit))) && ... 
            (any(find(BlockHandles(ind_Pass1(ct)) == OutputBlocksHit)));
        BlocksMarked(ind_Pass1(ct)) = 1;
    end
end