Mahalanobis Weighted k - 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[]

MWKMeansAdd[w_,B_,z_,Ind_:0]:=Module[{m=Length[B],k=Length[z],W=Total[w],PIpriv,tab,min,imin,PI,Scov,Pod,c,c0,n=Length[z[[1]]],pom,corr,It=0,epsilon=.000000001,eps=.0025},
Pod=Table[{w[[i]],B[[i,1]],B[[i,2]]},{i,Length[B]}]; 
AG=Table[{0,1},{i,3}];
c0=z; 
						(* Partition determining *)
PI=Table[{},{j,k}];Scov=Table[0,{j,k}];mj=Table[0,{j,k}];
Do[
tab=Table[EuclideanDistance[Pod[[i,2]],c0[[j]]],{j,k}]//N;
imin=PositionSmallest[tab][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI];API=PI;  
                                                           (* 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,": 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,"-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[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}];
PI=DeleteCases[PI,{}];
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; 
]; 
		(* The case of reducing the number of clusters in a partition *)
If[Length[c]<Length[c0],
F2[u_,v_,t_]:=Sum[w[[i]]Min[Append[Table[SquaredEuclideanDistance[B[[i,1]],c[[s]]],{s,Length[c]}],Norm[B[[i,1]]-{u,v,t}]^2]],{i,m}];
dom=AG;
								(* New DIRECT centers*)
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["New optimizer centars: ",c];
k=Length[c]; Print["Korigirani k: ",k];
								(* New partition *)
PI=Table[{},{j,k}];
Do[
tab=Table[EuclideanDistance[Pod[[i]][[2]],c[[j]]],{j,k}]//N;
imin=PositionSmallest[tab][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
];
{PI,c,Scov,It,F1}
]