Mahalanobis Weighted k-2 means  for the data with additional informations

Input: data set B⊂ (ℝ^3)x (Z^2)centers (Subscript[z, 1]),(Subscript[z, 2]),(Subscript[z, 3]) ∈ (([0,1])^3), weights (Subscript[w, 1]),⃛,(Subscript[w, m]) >0,number of iterations KK. Each additional center is obtained solving a GOP
using Compile;
Output: Optimal partition,  centers, covariance matrices number of iteration, Mahalanobis objective function value;
Submodules: dM[],WF[],[MWF[], DIRECT[]

MWKMeans2[w_,B_,z_,Ind_:0]:=Module[{m=Length[B],k=Length[z],W=Total[w],tab,min,imin,PI,c0,c,Scov,PIpriv,Pod,n=Length[z[[1]]],pom,It=0,epsilon=.000000001,eps=.0025},
Pod=Table[{w[[i]],B[[i,1]],B[[i,2]],i},{i,Length[B]}]; 
AG=Table[{0,1},{i,3}];
c0=z;
PI=Table[{},{j,k}];Scov=Table[0,{j,k}];mj=Table[0,{j,k}]; 
							(* 2-separable nethod *)
rest=Pod;
pol=Mean[c0];
zz=c0[[2]]-c0[[1]];
PI[[1]]=Select[rest,zz.(#[[2]]-pol)<=0&];
PI[[2]]=Select[rest,zz.(#[[2]]-pol)>0&];
					 (* Centers determining *)
c=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}];  
F0=WF[PI];  
If[Ind≠ 0,Print[It," |rest| = ",Length[rest],": F0(WF) = ",F0]];
If[Length[c]<Length[c0],
Print["There was a case of reducing the number of clusters in the partition"]];
							(* LOOP *)
While[
c0=c;It=It+1;
				(* Covariance matrices determining *)
Do[
mj[[j]]=Sum[PI[[j,s,1]],{s,Length[PI[[j]]]}];
Scov[[j]]=Covariance[PI[[j,All,2]]]
,{j,k}];
Do[
If[Min[Eigenvalues[Scov[[j]]]]<epsilon,Print["Cov.Mat. ",j,"-th center (",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}];
Do[
Do[
tab=Table[Sqrt[dM[PIpriv[[jj,s,2]],c[[j]],Scov[[j]]]],{j,k}];
imin=PositionSmallest[tab][[1]];
PI[[imin]]=Append[PI[[imin]],PIpriv[[jj,s]] ],
{s,Length[PIpriv[[jj]]]}],
{jj,k}];
k=Length[PI];
c=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}];
F1= MWF[PI,c,Scov];delF=Abs[F0-F1]/F0;absc=ChessboardDistance[c,c0];
If[Ind≠ 0,Print[It,": F1(MWF) = ",F1,"; |c-c0| = ",absc,"; delF = ",delF]];
absc>eps&&It<IT, F0=F1; 
]; 
{PI,c,Scov,It,F1}
]