function [A,B,C,E,StateNames] = smreal(A,B,C,E,StateNames)
%SMREAL  Compute a structurally minimal realization.
%
%   [A,B,C,E,NAMES] = SMREAL(A,B,C,E,NAMES) eliminates 
%   states that are not connected to any input or output.  
%   The resulting realization is structurally minimal, i.e., 
%   minimal when all nonzero entries of A,B,C,E are set to 1.
%   The argument NAMES contains the state names.

%   Author(s): P. Gahinet
%   Copyright (c) 1986-98 by The MathWorks, Inc.
%   $Revision: 1.1 $  $Date: 1998/05/05 13:56:17 $

% Uses characterization of the controllable subspace as smallest
% subspace S such that Im(B)<AS+ES, dim(S)=dim(AS+ES).
% This subspace is generated by
%     S0 = B;  Sk+1 = Sk + ASk + ESk
% To compute a structurally minimal realization, this recursion 
% is implemented on the boolean algebra (+ means OR, * means AND) 
% with A = A(A~=0) and b = any(B,2):
%     s(0) = b;  s(k+1) = A s(k) | E s(k) | s(k)
% The limit (stationary value) of the boolean vector s(k) indicates
% which states to keep.
%
% RE: The propagation equation is implemented as
%       ds(k+1) = (A+E) ds(k) & ~s(k)  
%       s(k+1) = s(k) + ds(k+1)
%     (ds(k+1) = newly activated states at iteration k)

% Dimensionality
sizes = size(A);
Colons = repmat({':'},[1 length(sizes)-2]);

% Derive structural connectivity information
Ab = any(A(:,:,:),3);
Bb = any(any(B(:,:,:),3),2);
Cb = any(any(C(:,:,:),3),1);
if ~isempty(E),
   Ab = Ab | any(E(:,:,:),3);
end

% Compute structurally minimal realization
% First identify structurally controllable states
xc = Bb;    % x(1)
dx = xc;    % dx(k) = x(k)-x(k-1)
while any(dx),
   Adx = any(Ab(:,dx),2);  % A dx(k) + E dx(k)
   dx = Adx & ~xc;         % dx(k+1)
   xc = xc | Adx;          % xc(k+1)
end

% Identify structurally observable states
xo = Cb;   
dx = xo;
while any(dx),
   Adx = any(Ab(dx,:),1);
   dx = Adx & ~xo;
   xo = xo | Adx;
end

% Minimal states are both controllable and observable
xco = xo & xc';

% Extract structurally minimal realization
if ~all(xco),
   A = A(xco,xco,Colons{:});
   B = B(xco,:,Colons{:});
   C = C(:,xco,Colons{:});
   if ~isempty(E),  
      E = E(xco,xco,Colons{:});
   end
   if nargin>4
      StateNames = StateNames(xco,1);
   end
end

