function [N,M,K,r,m]=normform(A)

% [N,M,K,r,m]=normform(A)
%       N es la forma normal de la matriz reducible nxn no negativa A.
%                Los estados de las clases transitorias van al principio.
%       M es la matriz de permutacion tal que N=M.' *A *M.
%       K es el cell array de las clases ordenadas.
%       r es el numero total de las clases transitorias + los estados sin retorno.
%       m es el numero total de bloques en la diagonal de N, o de clases +
%           estados sin retorno (asociados a un bloque diagonal (0) 1x1).
%
%    Ultima fecha de depuracion de errores: 12 de mayo de 2002
%    Autor: Juan-Miguel Gracia.

clc

n=length(A);


if irredble(A)==1,
    disp(' ')
    disp(' ')
    disp('                           ¡La matriz es irreducible! Su forma normal es ella misma.')
    disp('                                                      Final de la ejecucion.')
    disp('')
disp('                        Introduzca como argumento de NORMFORM una matriz reducible. ')
disp(' ')
disp(' ')
disp(' ')
disp(' ')
disp(' ')
disp(' ')

    break
end
disp(' ')
disp(' ')
disp(' ')
disp(' ')
disp(' ')
disp(' ')
disp(' ')
disp(' ')
disp(' ')
disp(' ')
disp('                                    ---------------------------------------------------------------------  ')
disp('                                                               ¡espere por favor! ')
disp('                                    ---------------------------------------------------------------------  ')
pause(2)


C=clases(A);
[O,F,S]=hasse(C,A);
T=clausura(O);

ini=transito(T);
%% Clases ordenadas segun el orden de los bloques en la diagonal de N.
z(1)=length(ini{1});
for k=2:length(ini)
      z(k)=z(k-1)+length(ini{k});
end

for i=1:length(ini)
    if i==1,
        for  j=1:z(1)
            K{j}=C{ini{1}(j)};
        end
    else
        if length(z)>=2,
            for  j=z(i-1)+1:z(i)
                K{j}=C{ini{i}(j-z(i-1))};
            end
        end
    end
end

if  isempty(S),
    M=permuta(K);
    N=M.'*A*M;
    r=length(S);
    m=length(K); break
else
    for j=1:length(F)
        K{z(length(ini))+j}=C{F(j)};
    end
    M=permuta(K);
    N=M.'*A*M;
    r=length(S);
    m=length(K);
end
%%===========================================================
%% Subfunciones requeridas por 'normform'
%%===========================================================
function C=clases(A)

% C=clases(A)
% da las clases de comunicacion de la matriz
% no negativa nxn A. Esto es equivalente a dar
% las componentes fuertemente conexas del
% digrafo cuya matriz de adyacencia es A.
%


clc

n=length(A);

R=1:n;%Conjunto a clasificar.
j=1; % Contador de clases.

for k=length(R):-1:1
      if isempty(R), break, end
      disp(' ')
      r=rem(k,3);
      switch r
      case 0
          disp('                                      ')
          disp('                                                                    ¡espere por favor! ')
          disp('                                      ')
      case 1
          disp('                                      ')
          disp('                                  ¡espere por favor! ')
          disp('                                      ')
      case 2
          disp('                                      ')
          disp('¡espere por favor! ')
          disp('                                      ')
      end
    c=nchoosek(R,k);
      s=size(c);
      i=1;
      while i<=s(1)
            if  length(c(i,:))>1,
                  if  irredble(A(c(i,:),c(i,:)))==1 
                        C{j}=c(i,:);  j=j+1;
                        R=setdiff(R,c(i,:)); 
                        if length(R)<s(2), % Las combinaciones restantes no interesan.
                              i=s(1)+1;
                        end
                  end
            else
                  C{j}=c(i,:); j=j+1;
                  R=setdiff(R,c(i,:));
                  if length(R)<s(2), % Las combinaciones restantes no interesan.
                        i=s(1)+1;
                  end
                  
              end
              i=i+1;         
      end
      
end



%%===========================================================
function [O,F,S]=hasse(C,A)
% [O,F,S]=hasse(C) da una matriz booleana O que describe la relacion binaria
%       correspondiente a un diagrama de Hasse de la relacion de orden
%       entre las clases C{1},C{2},...,C{m} de la matriz A.
%       Tambien da los indices de las clases finales F y los de las clases 
%       transitorias (y estados sin retorno) S.
%       Despues habra que hallar la clausura reflexiva y transitoria  de O.

m=length(C);
O=zeros(m);
for I=1:m
      for J=1:m
            if I==J,
                  O(I,J)=0;   
            elseif  any(any(A(C{I},C{J}))),
                  O(I,J)=1;
            else
                  O(I,J)=0;
            end
      end
end

F=[];S=[];
for I=1:m
      if sum(O(I,:))==0,
            F=[F,I];
      else
            S=[S,I];
      end
end
%%========================================================
function T=clausura(O)
% T=clasura(O) toma la matriz booleana O (una relacion binaria) y
%                    devuelve su clausura transitiva y reflexiva.

t=length(O);

T=eye(t);

for k=1:t
      P=potboole(O,k);
      for i=1:t
            for j=1:t
                  T(i,j)=T(i,j) | P(i,j);
            end
      end
end

%%===========================================================

function ini=transito(T)
% ini=transito(T) da la cadena inicial bajo el orden parcial T de clases
%   transitorias y estados sin retorno.

v=sum(T.');
t=length(T);
for i=1:t
      k=find(v==max(v));
      ini{i}=k;
      v(k)=0;
      if all(v<=1), break, end
end

%%============================================================
function M=permuta(K)
sigma=[];
for j=1:length(K)
      sigma=[sigma,K{j}];
end
p=length(sigma);

for i=1:p
      for j=1:p
            if i==sigma(j),
                  M(i,j)=1;
            else
                  M(i,j)=0;
            end
      end
end

