function [a,b,c,e,s,p] = aebalance(a,b,c,e,varargin)
%AEBALANCE  Balances (A,E) pair.
%
%   [A,B,C,E,S,P] = AEBALANCE(A,B,C,E) balances the (A,E) pair
%   and returns the balanced system (T\A*T,T\E*T,T\B,C*T) together 
%   with the scaling S and permutation P such that T(:,P) = diag(S).  
%  
%   [A,B,C,E,S,P] = AEBALANCE(A,B,C,E,Option1,Option2,...) specifies 
%   additional options as strings:
%     'perm'      Enables state permutation (default)
%     'noperm'    Prevents state permutation during balancing
%     'safebal'   Performs regularized row/column balancing
%                 (default)
%     'fullbal'   Performs full row/column balancing (equivalent 
%                 to BALANCE)
%
%   LOW-LEVEL UTILITY.

%   Authors: P. Gahinet
%   Copyright 1986-2020 The MathWorks, Inc.

% RE: Expects A,B,C,E to be matrices 2D arrays

% Get dimensions
nx = size(a,1);
ne = size(e,1);

% Quick exit when no state
if nx==0
   s = ones(0,1); 
   p = zeros(0,1);
   return
end

% Form matrix M = |A|+|E| to be scaled
% Note: Always best in DT, corresponds to emphasis on w=1 in CT
if ne==0
   mae = abs(a)+eye(nx);
else
   mae = abs(a)+abs(e);
end

% Perform scaling
[s,p] = mscale(mae,varargin{:}); 

% Apply scaling and permutation to (A,E)
is = 1./s;
a(p,p) = a .* (is * s');
if ne
   e(p,p) = e .* (is * s');
end

% Scale B,C if nonempty
if ~(isempty(b) && isempty(c))
   % B or C is nonempty
   b(p,:) = is .* b;
   c(:,p) = c .* s';
   bnorm = norm(b,1);
   cnorm = norm(c,inf);
   if bnorm>0 && cnorm>0
      sbc = pow2(round(log2(bnorm/cnorm)/2));
      b = b/sbc;  c = c*sbc;  s = s*sbc;
   end
end
