Mahalanobis weighted 𝓀 - means algorithm

(* 
Input: data set 𝒜 ⊂ (ℝ^n), data weights w ∈ (ℝ^n), initial centers (Subscript[z, 1]),⃛,(Subscript[z, k]) ∈ (ℝ^n), ;
Output: partition Π={(Subscript[π, 1]),⃛,(Subscript[π, k])} (each partition consists of clus-ters with corresponding weights associated 
to their elements), centers (Subscript[c, 1]),⃛,(Subscript[c, k]) ∈ (ℝ^n), covariance matrices (∑_j)∈(ℝ^(n\[Cross]n)), j=1⃛k, number of iterations It;
Submodules: dM[], DIRECT[] 
*)

(* Koristi se dM *)
MWKMeans[w_,A_,z_,Ind_:0]:=Module[{m=Length[w],k=Length[z],W=Total[w],PI0,PIpriv,tab,min,imin,Pod,c,c0,n=Length[z[[1]]],Scov,pom,corr,It=0,epsilon=.000000001},
Pod=Table[{w[[i]],A[[i]]},{i,m}];
c0=z;
						(* Računanje particije *)
PI=Table[{},{j,k}];Scov=Table[0,{j,k}];pom=Table[0,{j,k}];mj=Table[0,{j,k}];
Do[
tab=Table[Norm[Pod[[i]][[2]]-z[[j]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI];
                                                           (* Računanje centara *)
c=Table[Sum[PI[[j,s,1]]*PI[[j,s,2]],{s,Length[PI[[j]]]}]/Sum[PI[[j,s,1]],{s,Length[PI[[j]]]}],{j,k}];
							(* LOOP *)
While[
Chop[Norm[c0-c]]≠0,c0=c;It=It+1;
				(* Racunanje matrica Scov *)
Do[
mj[[j]]=Sum[PI[[j,s,1]],{s,Length[PI[[j]]]}];
Scov[[j]]=Sum[PI[[j,s,1]]*KroneckerProduct[(PI[[j,All,2]][[s]]-c[[j]]),(PI[[j,All,2]][[s]]-c[[j]])],{s,Length[PI[[j]]]}]/mj[[j]]
,{j,k}];
Do[
If[Min[Eigenvalues[Scov[[j]]]]<epsilon,Print["Cov.Mat. ",j,"-tog centra (",c[[j]],"): ",Scov[[j]]," je singularna"];Print["|PI[j]|: ",Length[PI[[j]]],", Eigen_min=",Min[Eigenvalues[Scov[[j]]]]];Scov[[j]]=IdentityMatrix[n] ]
,{j,k}];
PIpriv=PI;
PI=Table[{},{i,k}];
corr=Table[Log[m^2 Det[Scov[[jj]]]/mj[[jj]]^2],{jj,k}];
Do[
Do[
tab=Table[dM[PIpriv[[jj,s,2]],c[[j]],Scov[[j]]],{j,k}];
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],PIpriv[[jj,s]] ],
{s,Length[PIpriv[[jj]]]}],
{jj,k}];
PI=DeleteCases[PI,{}];
k=Length[PI];
c=Table[Sum[PI[[j,s,1]]*PI[[j,s,2]],{s,Length[PI[[j]]]}]/Sum[PI[[j,s,1]],{s,Length[PI[[j]]]}],{j,k}];
]; 
				(* Slucaj smanjenja broja klastera u particiji *)
If[Length[c]<Length[c0],
F2[u_,v_]:=Sum[w[[i]]Min[Append[Table[Norm[A[[i]]-c[[s]]]^2,{s,Length[c]}],Norm[A[[i]]-{u,v}]^2]],{i,m}];
dom={xG,yG};
								(* Novi DIRECT centar*)
minp=DIRECT[F2,dom,DMin->10.^(-3),MaxIter->50, Eps->10^(-4)]//Timing; min=minp[[2]];Print["DIRECT: ",min];
c=Append[c,min[[1]]];
Print["Novi optimirani centar: ",c];
k=Length[c]; Print["Korigirani k: ",k];
								(* Nova particija *)
PI=Table[{},{j,k}];
Do[
tab=Table[Norm[Pod[[i]][[2]]-c[[j]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
];
{PI,c,Scov,It}
]