function [A, Y, H, cycles, HBindex] = fuzzy_c_mean (X, c, f, errstop, lim) 
% [A, Y, H, cycles, HBindex] = fuzzy_c_mean (X, c, f, errstop, lim). Esegue l' algoritmo di fuzzy c means secondo Bezdek
% Input:
% X      : matrice dei dati di input.
%          Formato: [n,m]=size(X) con n numero dei pattern ed m dimensione dei pattern.
% c      : numero delle classi (def. c = 2)
% f      : indice di fuzziness dell' algoritmo, 1<f<inf (def. f = 2)
% errstop: condizione di arresto  (def. = 0.001)
% lim    : limite valore di membership per assegnazione di un pattern ad un cluster
%          se lim = 0, assegnazione crisp. (def. lim = 0)
% Output:
% A      : Matrice delle membership function dei cluster.
%          Formato: [c,n]=size(A)
%          c = numero dei cluster
%          n = numero dei pattern di input
% Y      : Matrice dei centri dei cluster
%          Formato: [c,m]=size(Y)
%          c = numero dei cluster
%          m = dimensione dei pattern
% H      : Vettore di assegnazione dei pattern ai cluster
%          Formato: n=length(H)
%          n = numero dei pattern di input
%          H(i) memoriza il cluster  di appartenenza dell' iesimo pattern. 
%          L' etichetta del cluster corrisponde all' indice della riga
%          della matrice A. H(i) = 0 significa che l' iesimo pattern non è
%          stato assegnato ad alcun cluster
% cycles : numero dei cicli effettuati dall' algoritmo
% HBindex  : Xie-Beni index (misura la qualità della partizione)
%            minore è il valore di HBindex, migliore è la partizione

% Controlla Input e setta i default
if nargin == 1
    c = 2;
    f = 2;
    errstop = 0.001;
    lim = 0;
elseif nargin == 2
    f = 2;
    errstop = 0.001;
    lim = 0;
elseif nargin == 3
    errstop = 0.001;
    lim = 0;
elseif nargin == 4
    lim = 0;
end
[n,m]=size(X);
if n < c
    error('Number of patterns must be at least equal to the number of clusters');
end

% Init variables
err = Inf;
cycles = 0;

% Step1: Selezione di una pseudopartizione iniziale
A = DoInitPseudoPartition(c,n);

while err > errstop
	%Step2: Calcolo del centro dei cluster
	Y = ClusterCenters(A, X, f);

	% Salvataggio di A
	oldA = A;

	% Step3: Aggiorna la matrice di memberhip dei cluster
	A = UpdateMembership(X, Y, c, f);

	% Step4: Calcola l'errore fra la A e la membership oldA del passo precedente
	err = Error(A, oldA);
	cycles = cycles + 1;
end
H = AssignPatterns(A, lim);
HBindex = ComputeXieBeniIndex(X,Y,A,f);
% Disegna il grafico dei cluster
% drawfuzzycluster(X, A, lim);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Funzioni utilizzate da fuzzy_c_mean %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function U = DoInitPseudoPartition(c,n)
% Calcola una pseudopartizione iniziale
% Le membership sono assegnate secondo una distribuzione uniforme
% mantenendo il vincolo probabilistico di Ruspini
% 
% Input:
% c      : numero delle classi
% n      : numero dei pattern di input
% Output:
% U      : Matrice delle membership func  tion dei cluster.
%          Formato: [c,n]=size(U)
%          c = numero dei cluster
%          n = numero dei pattern di input
U = rand(c,n);
sumcol = sum(U);
for t=1:n
    U(:,t)=U(:,t)/sumcol(t);
end



function Y = ClusterCenters(A, X, f)
% Calcola il centro dei cluster
% Input:
% A      : Matrice delle membership funtion dei cluster.
%          Formato: [c,n]=size(A)
%          c = numero dei cluster
%          n = numero dei pattern di input 
% X      : matrice dei dati di input. 
%          Formato: [n,m]=size(X) con n numero dei pattern ed m dimensione dei pattern.
% f      : indice di fuzziness dell' algoritmo, 1<f<inf
% Output:
% Y      : Matrice dei centri dei cluster
%          Formato: [c,m]=size(Y)
%          c = numero dei cluster
%          m = dimensione dei pattern
[n,m]=size(X);
[c,n]=size(A);
Y = zeros(c,m);
% Si puo' migliorare con notazione matriciale *********
for i=1:c
  J = A(i,:).^f;
  %for t=1:n
  %  Y(i,:) = Y(i,:) + X(t,:)*J(t);
  %end
  Y(i,:) = Y(i,:) + J*X;
  Y(i,:) = Y(i,:)/sum(J);
end





function A = UpdateMembership(X, Y, c, f)
% Aggiorna la matrice di memberhip dei cluster
% Input:
% X      : matrice dei dati di input. 
%          Formato: [n,m]=size(X) con n numero dei pattern ed m dimensione dei pattern.
% Y      : Matrice dei centri dei cluster
%          Formato: [c,m]=size(Y)
%          c = numero dei cluster
%          m = dimensione dei pattern
% c      : numero dei cluster
% f      : indice di fuz  ziness dell' algoritmo, 1<f<inf
% Output:
% A      : Matrice delle membership funtion dei cluster.
%          Formato: [c,n]=size(A)
%          c = numero dei cluster
%          n = numero dei pattern di input 

esp = 1/(f-1);
[n,m]=size(X);
A=zeros(c,n);
distYX=zeros(c,n);

% Calcola le distanze fra i centri dei cluster ed ogni pattern
for k=1:n
  for i=1:c 
    distYX(i,k) = Distance(Y(i,:),X(k,:));
  end
end
% Cerca le distanze nulle su ogni colonna di distYK.
% Una distanza e' considerata nulla se inferiore a eps
% Se le trova mette nelle A(i,k) corrispondenti un valore val = 1/numerozeri. Pone quindi
% gli altri A(i,k) a zero.
for k=1:n
  posiz0 = find(abs(distYX(:,k))<eps);
  numzeros = length(posiz0);
  if numzeros ~= 0  % esistono distanze nulle per il kesimo pattern
     val = 1/numzeros;
     for t=1:c
       if distYX(t,k) < eps
         A(t,k) = val;
       else
         A(t,k) = 0;
       end
     end     
  else % non esistono distanze nulle per il kesimo pattern
  % Calcola gli A(i,k) per cui distYK<>0   
    dden = 0;
    % usare calcolo vettoriale per efficienza ******
	for j=1:c
	    dden = dden + (1/distYX(j,k))^esp;		  
    end
    for i=1:c 
    	dnum = (1/distYX(i,k))^esp;        
		A(i,k) = dnum/dden;
	end 
  end % end else
end

   
function d = Distance(x,y)
% Calcola il quadrato della distanza Euclidea fra i vettori x ed y
d = (x-y).^2;
d = sum(d);
 


function err = Error(A, oldA)
% Calcola l' errore, in termini di massimo delle distanze fra le membership a passi successivi
%Input:
% A      : Matrice delle membership funtion dei cluster.
%          Formato: [c,n]=size(A)
%          c = numero dei cluster
%          n = numero dei pattern di input 
% oldA   : Matrice delle membership funtion dei cluster al passo precedente
%Output:
% err    : errore calcolato come distanza massima fra le membership a passi successivi
oldA = abs(A-oldA);
err = max(max(oldA));


function H = AssignPatterns(A, lim)
% Assegna i pattern alle classi in modo crisp
%Input:
% A      : Matrice delle membership funtion dei cluster.
%          Formato: [c,n]=size(A)
%          c = numero dei cluster
%          n = numero dei pattern di input 
% lim   :  limite valore di membership per assegnazione di un pattern ad un cluster
%          se lim = 0, assegnazione al cluster con membership massima. 
%Output:
% H      : Vettore di assegnazione dei pattern ai cluster
%          Formato: n=length(H)
%          n = numero dei pattern di input
%          H(i) memoriza il cluster  di appartenenza dell' iesimo pattern. 
%          L' etichetta del cluster corrisponde all' indice della riga
%          della matrice A. H(i) = 0 significa che l' iesimo pattern non è
%          stato assegnato ad alcun cluster

% 
[m,n]=size(A);
for k=1:n
    maxvalue = -1;
    maxcluster = 0;
    for j=1:m
        if A(j,k)>maxvalue
            maxvalue = A(j,k);
            maxcluster = j;
        end
    end
    if maxvalue < lim
        maxcluster = 0;
    end
    MaxVal(k).value = maxvalue;
    MaxVal(k).cluster = maxcluster;
end

H = [MaxVal.cluster];
    

function HBindex = ComputeXieBeniIndex(X,Y,A,f) 
% HBindex = ComputeXieBeniIndex(X,Y,A). Calcola l' indice di Xie-Beni
% L' indice di Xie-Beni misura la qualità di una partizione secondo il
% criterio della massima uniformità all' interno di un cluster e della
% massima diversità fra i cluster. In pratica è il rapporto fra la somma
% pesata delle distanza dei pattern dai centri dei cluster sommata su tutti
% i cluster e la distanza minima fra i centri dei cluster.
%Input:
% X      : matrice dei dati di input.
%          Formato: [n,m]=size(X) con n numero dei pattern ed m dimensione dei pattern.
% Y      : Matrice dei centri dei cluster
%          Formato: [c,m]=size(Y)
%          c = numero dei cluster
%          m = dimensione dei pattern
% A      : Matrice delle membership funtion dei cluster.
%          Formato: [c,n]=size(A)
%          c = numero dei cluster
%          n = numero dei pattern di input 
% f      : indice di fuzziness dell' algoritmo, 1<f<inf
%Output:
% HBindex  : Xie-Beni index
%            minore è il valore di HBindex, migliore è la partizione

% Calcola le distanze fra i centri dei cluster ed ogni pattern
[c,n] = size(A);
WithinClusterDistance = 0;
for i=1:c 
    for k=1:n
    WithinClusterDistance = WithinClusterDistance + A(i,k)^f * Distance(Y(i,:),X(k,:));
  end
end

DistCenters = zeros(1,(c*(c-1))/2);
k=0;
for i=1:c-1 
    for j=i+1:c
        k=k+1;
        DistCenters(k) = Distance(Y(i,:),Y(j,:));
    end
end
 
minDistCenters = min(DistCenters);
HBindex = WithinClusterDistance/(c*minDistCenters);


                