% donnees de base: R et D
% version basee sur calcul matriciel
% a fixer trunc = 0 ou 1  
% a fixer coup, niveau de coupure (pour utopia: 0.01)

[n,m] = size(R) ;
L1 = ones(n,1);
C1 = ones(m,1);
D  = D';

% somme des colonnes de R sans 0

RBAR = diag(1 ./ max(L1,R * C1))*R;
DTAR = diag(1 ./ max(C1,D * L1))*D;

GBAR = RBAR*DTAR ;

I = eye(n) ;

% les trois matrices de structure

S = (I-0.5*GBAR)^-1 ;
SR = S*RBAR ;
SD = S'*DTAR' ;

% S tilda 
% on met a 0 ce qui est inferieur a coup
 
ST = S;

if trunc
  for i=1:n,
   for j=1:n,
     x = ST(i,j);
     if x< coup 
       ST(i,j) = 0;
      else
       ST(i,j) = 1;
     end
   end 
   %ST(i,i)=0;
  end

end

% Calcul des coordonnees, affichage et tri

CX = sum(ST,2);
CY = sum(ST,1);

CX = n * L1 - CX;
CY = n * L1 - CY';

plot(CX,CY,'bo');

coord = sortrows([CX CY (1:n)']);

% donnees de base ST, n,m
% agragation des unites d'information liees (cos < coup1)
% (sorte d'analyse cluster)
% liW nombre de composantes
% W contient par ligne: 
% le nombre de composante et le no de chacune d'entre elle
% P contient par colonne les unites reellement liees
% Isol contient les unites isolees
% on pourrait calculer les longueurs une fois pour toute (corrcoef)
% ou toutes les correlations puis choisir les meilleures !

% symetrisation
V = ST + ST';

% mise a 1
V = V & 1;

% calcul des normes des lignes
poids = sqrt(sum(V,2));
 
coup1 = 0.3;

% unites triees par poids decroissant
% 1ere colonne les normes
% 2eme colonne no des UI

% version top-down
V = flipud(sortrows ([poids (1:n)' V])); 

% version bottom-up
%V = sortrows ([poids (1:n)' V]); 

% exclusion des isolees (norme = 1)
lmax = n;
%for i = 1:n;
%  if V(i,1) == 1
%    lmax = i-1;
%    break;
%  end
% end

%if lmax ~= n,
%  Isol = V(lmax+1:n,2) ;
%end

W = zeros(lmax,lmax+1);
P = zeros(lmax+1,lmax);

% demarrage de l'agregation
% unite W courante: 

noW = V(1,2);
%nbre de W definie 
liW = 1;
W(1,1) = 1;
W(1,2) = noW;

% P les composantes en vertical (1e ligne:les poids)
P(:,1) = [V(1,1) V(1,3:lmax+2)]';

for i=2:lmax,
  % produit de tous les W avec vi
  AUX = (((V(i,3:lmax+2) * P(2:lmax+1,1:liW) ) 
                            ./ P(1,1:liW) )./ V(i,1));
  LP = [AUX' (1:liW)'];
  LP = flipud(sortrows(LP));
  
  if(LP(1,1) > coup1)

  % agreger V(i,:) a W(LP(1,2)) et a P 
  % recalculer la norme de P(:,LP(1,2))
  
    ct = LP(1,2);
   
    W(ct,1) = W(ct,1)+1;
    W(ct,W(ct,1)+1) = V(i,2);
    P(2:lmax+1,ct) = (P(2:lmax+1,ct) + V(i,3:lmax+2)') & 1;
    P(1,ct) = sqrt(sum(P(2:lmax+1,ct)));
  else 
    noW = V(i,2);
    liW = liW+1;
    W(liW,1) = 1;
    W(liW,2) = noW;
    P(:,liW) = [V(i,1) V(i,3:lmax+2)]';
   end
end

W = W(1:liW,:);
W1 = flipud(sortrows (W));

P = P(:,1:liW);

Idées d'amélioration:

1. mettre une variable (booléenne qui commande la version top-down ou bottom-up.
2. introduire une répétition: fixer une valeur élevée de coup1, calculer P, recommencer avec une valeur inférieure de coup1 et V = P', etc. A voir pour le critère d'arrêt.