function []=kmeans(a,b,K)
[x,n]=generare(a,b,K); // preluare date
C=init(x,n,K); // etapa de initializare
omega=[];   // structura pentru stocarea claselor
nomega=zeros(K);
d=[];
disp(K,"K=");
for t=1:10
    omega=[];   // structura pentru stocarea claselor
    nomega=zeros(1,K); // nr de elemente din fiecare clasa
    for el=1:n
//        disp('el'); disp(el);
        for k=1:K
           d(k)=dist(x(:,el),C(:,k));
        end   
//        disp('d='); disp(d);
        kmin=1;
        for k=2:K
            if d(k)<d(kmin) 
                kmin=k;
            end
        end    
//        disp(t,'t=',kmin,'kmin='); 
//        disp('nomega'); disp(nomega);
        nomega(kmin)=nomega(kmin)+1;
        omega(kmin,:,nomega(kmin))=x(:,el); 
    end
    for k=1:K
//        disp(k);
//        disp(omega(k));
        C(:,k)=medie(omega(k));
    end
end
// clasificare finala
omega=[];   // structura pentru stocarea claselor
nomega=zeros(1,K);
  for el=1:n
//        disp('el'); disp(el);
        for k=1:K
           d(k)=dist(x(:,el),C(:,k));
        end   
//        disp('d='); disp(d);
        kmin=1;
        for k=2:K
            if d(k)<d(kmin) 
                kmin=k;
            end
        end    
//        disp('kmin'); disp(kmin);
        nomega(kmin)=nomega(kmin)+1;
        omega(kmin,:,nomega(kmin))=x(:,el); 
    end
disp(omega,'omega=');    
culori=['ro','go','bo','ko','yo','co'];
culori2=['r*','g*','b*','k*','y*','c*'];
for k=1:K
    plot(omega(k,1,:)',omega(k,2,:)',culori(k));
    plot(C(1,k),C(2,k),culori2(k));
end 
endfunction

function med=medie(omega)
dim=size(omega);
med=[];
med=zeros(dim(1),1);
for i=1:dim(2)
    med=med+omega(:,i);
end
for i=1:dim(1)
    med(i,1)=med(i,1)/dim(2);
end
endfunction

// cautarea unui element
function rez=prezent(ind,m,i)
rez=0;
for j=1:m
    if ind(j)==ind(i) 
        rez=1;
    end
end    
endfunction        

// initializarea prototipurile prin selectie aleatoare din setul de date
function C=init(x,n,K)
ind=[];
for i=1:K
    ind(i)=grand(1,1,"uin",1,n);
    disp(ind(i));
    while prezent(ind,i-1,i)==1
        ind(i)=grand(1,1,"uin",1,n);
    end
end    
for i=1:K
   C(:,i)=x(:,ind(i));
end   
endfunction

// calcul functie distanta
function d=dist(x,y)
n=length(x);
d=0;
for i=1:n
    d=d+(x(i)-y(i))^2;
end    
endfunction

function [x,n]=generare(a,b,K)
x=[]; n=0; sigma=(b-a)/(4*K);
for k=1:K
  nr=grand(1,1,"uin",5,50);
  n=n+nr;
  xc=grand(1,1,"unf",a,b);
  yc=grand(1,1,"unf",a,b);
  for i=1:nr
    x=[x [grand(1,1,"nor",xc,sigma),grand(1,1,"nor",yc,sigma)]];
  end
end
endfunction

kmeans(0,1,4)
