clear all
close all
clc

%%% genera dati
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MU1=[0,0];
SIGMA1=[0.3 0;0 2];
Npoints=2000;
X1 = mvnrnd(MU1,SIGMA1,Npoints);
X1=X1';
%plot(X1(1,:),X1(2,:),'co')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%MU2=[7,7];
MU2=[5,5];
SIGMA2=[1.5 0.7;0.7 1.5];
X2 = mvnrnd(MU2,SIGMA2,Npoints);
X2=X2';
%hold on
%plot(X2(1,:),X2(2,:),'bo')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MU3=[-5,-10];
SIGMA3=[1.5 0.7;0.7 1.5];
X3 = mvnrnd(MU3,SIGMA3,Npoints);
X3=X3';
%hold on
%plot(X3(1,:),X3(2,:),'ro')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MU4=[10,-10];
SIGMA4=[1.5 0;0 1.5];
X4 = mvnrnd(MU4,SIGMA4,Npoints);
X4=X4';
%hold on
%plot(X4(1,:),X4(2,:),'mo')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%MU5=[0,-10];
%MU5=[-5,10];
MU5=[-5,5];
SIGMA5=[1 -0.8;-0.8 1];
X5 = mvnrnd(MU5,SIGMA5,Npoints);
X5=X5';
%hold on
%plot(X5(1,:),X5(2,:),'ko')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%X=[X1 X2 X3 X4 X5];
X=[X1 X2 X4 X5];
%X=[X1 X2];

[DIM,Ntot]=size(X);
plot(X(1,:),X(2,:),'bo')

%load  DatiGrazieLuigi
%[DIM,Ntot]=size(X);
%plot(X(1,:),X(2,:),'bo')

% figure
% plot(X(1,:))
hold on
% plot(X(2,:),'r')
%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Ktot=2; %%% numero clusters maximo 
%%%%%%%%%%%%%%%
%for k=1:Ktot
 k=4; %%% clusters ahora   

    %%% posiciones iniciales
C=-5+10*rand(2,k);
Ctot1=C(1,:);
Ctot2=C(2,:);

%plot(C(1,1),C(2,1),'sk','MarkerEdgeColor','r','MarkerFaceColor','k','MarkerSize',10)
%plot(C(1,2),C(2,2),'sb','MarkerEdgeColor','r','MarkerFaceColor','b','MarkerSize',10)
%plot(C(1,3),C(2,3),'sm','MarkerEdgeColor','r','MarkerFaceColor','m','MarkerSize',10)
 
%%% posiciones iniciales - plot
for j=1:k
  plot(C(1,j),C(2,j),'sk','MarkerEdgeColor','r','MarkerFaceColor','r','MarkerSize',10)
end
%%%%%%%
iter=10;
%%%%%%%%%
for r=1:iter
r

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=2;
dist=@(x,z) sum((x-z).^a);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:Ntot
    for j=1:k
       D(i,j)=dist(X(:,i),C(:,j));
    end
     clust_index(i)=find(D(i,:)==min(D(i,:)));
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%% for Ktot=3
%        if clust_index(i)==1
%         plot(X(1,i),X(2,i),'ok')
%        elseif clust_index(i)==2
%            plot(X(1,i),X(2,i),'ob')
%        else
%            plot(X(1,i),X(2,i),'om')
%        end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for j=1:k
  pos_j=find(clust_index==j);
    if isempty(pos_j)==0
      C(:,j)=mean(X(:,pos_j),2);
        E(r,j)=sum(var(X(:,pos_j)'));
    end
      
end
%plot(C(1,1),C(2,1),'sk','MarkerEdgeColor','r','MarkerFaceColor','k','MarkerSize',10)
%plot(C(1,2),C(2,2),'sb','MarkerEdgeColor','r','MarkerFaceColor','b','MarkerSize',10)
%plot(C(1,3),C(2,3),'sm','MarkerEdgeColor','r','MarkerFaceColor','m','MarkerSize',10)
 
Ctot1(r+1,:)=C(1,:);
Ctot2(r+1,:)=C(2,:);


for j=1:k
  plot(C(1,j),C(2,j),'dk','MarkerEdgeColor','k','MarkerFaceColor','c','MarkerSize',10)
end
   
pause

%cla

Energy_evolution=sum(E,2)

end %%% iter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%% posiciones finales
for j=1:k
  plot(C(1,j),C(2,j),'pk','MarkerEdgeColor','k','MarkerFaceColor','g','MarkerSize',15)
end

%%%%%%%%%%
%%%%%%%%%%
%  figure
%  plot(Ctot1(:,2),'r')
%  hold on
%  plot(Ctot1(:,1),'m')
%  plot(Ctot1(:,3),'b')
%  figure
%  plot(Ctot2(:,2),'r')
%  hold on
%  plot(Ctot2(:,1),'m')
%  plot(Ctot2(:,3),'b')

%Etot=sum(E,2);
%Etot(end)

%Elbow_AIC(k)=Etot(end)+k;
%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%figure 
%plot(exp(-Elbow_AIC),'o--')