Mahalanobis Weighted k - means with partitions saving
Mahalanobis Incremental algorithm

Mahalanobis Clustering is a complete Mathematica-program for Maha-lanobis clustering using the normalized Mahalanobis distance-like function.All modules needed are listed in the preamble.One has to load one of the
suggested examples first,which will automatically perform the MWKMeans[]
algorithm.Then one can call the incremental algorithm MInc[] which will
produce several partitions sorted in the file PART.txt.Calling these partitions
by MCH[],MDB[],MSSWC[],and MArea[] modules,one obtains appropriate
indexes,based on which MAPart is suggested.

Modules

Needs["ComputationalGeometry`"]

Squared Mahalanobis distance

dm[x_,y_,S_]:=(x-y).Inverse[S].(x-y);
dM[x_,y_,S_]:=(Det[S]^(1/n))(x-y).Inverse[S].(x-y);

Mahalanobis objective functions

MWF[PI_,c_,Scov_]:=Module[{m=Length[Flatten[PI,1]],n=Length[Scov]},
Return[Sum[PI[[j,s,1]]*dM[PI[[j,s,2]],c[[j]],Scov[[j]]],{j,Length[PI]} ,{s,Length[PI[[j]]]}] ]
];

MWFmin[PI_,c_,Scov_]:=Module[{A=Flatten[PI,1][[All,2]],m=Length[A],k=Length[c]},
Sum[Flatten[PI,1][[All,1]][[i]]Min[Table[dM[c[[j]],A[[i]],Scov[[j]]],{j,k}]],{i,m}]
]

Mahalanobis weighted k - means algorithm

(* 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[Hdist[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}
]

Expectation Maximization algorithm

(* Koristi se dm *)
EM[w_,A_,z_,Ind_:0]:=Module[{m=Length[A],n=Length[A[[1]]],k=Length[z],It=0,MaxIt=30,eps=.005,z0},
				(* Initialization *)
z0=z;Sigma0=Table[IdentityMatrix[n],{i,k}];P0=Table[1/k,{i,k}];lam=Table[0,{i,m}];
p[u_,mu_,Cov_]:=Exp[-(u-mu).Inverse[Cov].(u-mu)/2]/(Sqrt[2Pi Det[Cov]]);
				(* Iterative process *)
While[It=It+1;
PI=Table[{},{j,k}];
Do[
lam[[i]]=Table[(P0[[j]] p[A[[i]],z0[[j]],Sigma0[[j]]])/Sum[p[A[[i]],z0[[s]],Sigma0[[s]]]P0[[s]],{s,k}],{j,k}]
,{i,m}];
Do[
imin=Ordering[lam[[i]],-1][[1]];
PI[[imin]]=Append[PI[[imin]],A[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI];
If[Ind≠0,
slPI=ListPlot[Table[PI[[j]],{j,k}],PlotStyle->Table[Hue[j/k],{j,k}]];
sl0=ListPlot[{z0},PlotStyle->{Brown,Opacity[.8],PointSize[.025]}];
Do[
rad=Max[Table[dm[z0[[j]],PI[[j,s]],Sigma0[[j]]],{s,Length[PI[[j]]]}]];
slcon[j]=ContourPlot[dm[z0[[j]],{x,y},Sigma0[[j]]]==75rad/100,{x,xG[[1]],xG[[2]]},{y,yG[[1]],yG[[2]]}],{j,k}];
sl=Show[Table[slcon[j],{j,k}]];
Print[Show[slPI,sl,sl0,Ticks->{{-.5,.5},{-.5,.5}},AxesOrigin->{0,0},ImageSize->Small]];
];
z1=Table[Sum[lam[[i,j]]A[[i]],{i,m}]/Sum[lam[[i,j]],{i,m}],{j,k}];
Chop[Hdist[z0,z1]]>eps&&It≤ MaxIt,
Sigma1=Table[Sum[lam[[i,j]]KroneckerProduct[(A[[i]]-z0[[j]]),((A[[i]]-z0[[j]]))],{i,m}]/Sum[lam[[i,j]],{i,m}],{j,k}];
P1=Table[Sum[lam[[i,j]],{i,m}]/m,{j,k}];
If[Ind≠0,
Print[It,". ITERACIJA: \nCenters: ",z1,"\nCov: ",Table[Sigma1[[j]]//MatrixForm,{j,Length[z0]}],"\nP: ",P1,"\nHdist= ",Hdist[z0,z1]];
];
P0=P1;Sigma0=Sigma1;z0=z1;
];
{PI,z1,Sigma1,It,P1}
]

Mahalanobis weighted Incremental algorithm

(* Ako je poznata originalna particija PI0, možemo je usporediti s particijom dobivenom u brk-koraku *)
MInc[w_,A_,cen_,KK_,Ind_:0]:=Module[{centri,centriN,m=Length[A],n=Length[A[[1]]],k,dom,Scov,min,PI1},
Print["n=",n];
(* Initialization *)
xG={Min[A[[All,1]]]-.5,Max[A[[All,1]]]+.5};yG={Min[A[[All,2]]]-.5,Max[A[[All,2]]]+.5};
Scov=Table[IdentityMatrix[n],{i,Length[cen]}]; 
sol=MWKMeans[w,A,cen]; 
(*cen=sol[[2]];*)
centri=sol[[2]];
brk0=Length[centri];
PART={};
					(* Iterativni postupak *)
Do[
If[Ind≠0,Print["\nNa postojecih ",Length[centri]," centara dodaje se 1 novi"]];
k=Length[centri];
						(* Compile 
FCompile=Compile[{{u,_Real},{v,_Real}},
Evaluate[Sum[w[[i]]Min[Append[Table[Norm[A[[i]]-centri[[s]]]^2,{s,Length[centri]}],Norm[A[[i]]-{u,v}]^2]],{i,m}] ]];	
F2[u_,v_]:=Apply[FCompile,{u,v}];
*)
F2[u_,v_]:=Sum[w[[i]]Min[Append[Table[Norm[A[[i]]-centri[[s]]]^2,{s,Length[centri]}],Norm[A[[i]]-{u,v}]^2]],{i,m}];
dom={xG,yG};
	(* DIRECT *)
minp=DIRECT[F2,dom,DMin->10.^(-3),MaxIter->50, Eps->10^(-4)]//Timing; min=minp[[2]];If[Ind≠0,Print["DIRECT: ",min]];
centriN=Append[centri,min[[1]]];
If[Ind≠0,Print["Novi optimirani centar: ",centriN]];
(* Ovdje se vrijednost funkcije cilja ne moze racunati kao MWFmin jer jos nemam Scov[[k+1]] *)
k=Length[centriN];
					(* Poziv Mahalanobis k-means *)
t1=Timing[sol=MWKMeans[w,A,centriN]];
If[Ind≠0,Print["Novi MWKMeans profinjeni centar: ",sol[[2]]]];
Scov=sol[[3]];
If[Ind≠0,Print["MWKMeans iteracija = ",sol[[4]],"; MWF = ",MWF[sol[[1]],sol[[2]],sol[[3]]],";  MWFmin = ",MWFmin[sol[[1]],sol[[2]],sol[[3]]],"; mj: ",Table[Length[sol[[1,j]]],{j,Length[PI]}]]];
(* Ako je poznata originalna particija, ovdje ju možemo usporediti s particijom dobivenom u brk-koraku:
PI1=Table[sol[[1,j,All,2]],{j,Length[sol[[1]]]}];
If[brk≤Length[PI0],rand=Rand[PI1,PI0],rand=Rand[PI0,PI1]];
Print["Confusion Matrix = ",rand[[1]]//MatrixForm, "; Rand = ",rand[[2,1]],"; Jaccard = ",rand[[2,2]]];
*)
If[Ind≠0,Print["Time (DIRECT+WKM) = ",minp[[1]],"+",t1[[1]]," = ",minp[[1]]+t1[[1]]]];
Print[brk,"-particija"];
     (* FIGUREE *)
k=Length[sol[[1]]];
If[n==2,
slc=ListPlot[sol[[2]],PlotStyle->{Gray,PointSize[.05],Opacity[.7]}];
slPI=ListPlot[Table[sol[[1,j,All,2]],{j,k}],PlotStyle->{Table[PointSize[.015],{j,k}]},PlotStyle->{Table[Hue[j/k],{j,k}]}]
];
							(* Elipse *)
If[n==2,
Do[
mj=Length[sol[[1,j]]];
tab=Sort[Table[dM[sol[[2,j]],sol[[1,j,s,2]],sol[[3,j]]],{s,mj}]];
rad=tab[[Ceiling[95mj/100]]];
slcon[j]=ContourPlot[dM[sol[[2,j]],{x,y},sol[[3,j]]]==rad,{x,xG[[1]],xG[[2]]},{y,yG[[1]],yG[[2]]}]
,{j,k}];
sl=Show[Table[slcon[j],{j,k}],ImageSize->Small];
sl2=Print[Show[slPI,slc,sl,GridLines->Automatic,AspectRatio->Automatic,ImageSize->Small]];	
];		 
centri=sol[[2]];     
PART=Append[PART,sol[[1]]]; 
,{brk,KK}];
PART
]

Hausdorff distance between two sets

(* 
Po defaultu koristi se euklidska metrika. Ako se želi druga, treba to na dva mjesta upisati
Input: finite sets A,B ⊂ (ℝ^n); 
Output: Hausdorff distance;
Submodules: none needed 
*)
Hdist[a_,b_]:=Module[{m,n,a1,b1},
If[Length[a]>Length[b],a1=b;b1=a,a1=a;b1=b];
m=Length[a1];n=Length[b1];
Max[Max[Table[Min[Table[Norm[a1[[i]]-b1[[j]]],{j,n}]],{i,m}]],Max[Table[Min[Table[Norm[a1[[i]]-b1[[j]]],{i,m}]],{j,n}]]]]

Rand and Jaccard index and the confusion matrix

(* 
Input: two partitions (Π^((1))), (Π^((2))) of the set 𝒜 ⊂ (ℝ^n); 
Output: Rand index R, Jaccard index J, confusion matrix S, CPU-time;
Submodules: none needed 
*)
Rand[PI1_,PI2_]:=Module[{m=Length[Flatten[PI1]],S,S1=Table[0,{i,Length[PI1]}],S2=Table[0,{j,Length[PI2]}],a,b,c,d},
time=Timing[
(* Confusion matrix *)
S=Table[0,{r,Length[PI1]},{r,Length[PI2]}];
Do[
Do[
S[[r,s]]=Length[Intersection[PI1[[r]],PI2[[s]]]]
,{r,Length[PI1]}]
,{s,Length[PI2]}];
(* Numbers a, b, c, d *)
TT=Sum[S[[i,j]]^2,{i,Length[PI1]},{j,Length[PI2]}];
Do[
S1[[i]]=Sum[S[[i,j]],{j,Length[PI2]}]
,{i,Length[PI1]}];
Do[
S2[[j]]=Sum[S[[i,j]],{i,Length[PI1]}]
,{j,Length[PI2]}];
a=.5(TT-m);
b=.5(Sum[S1[[i]]^2,{i,Length[PI1]}]-TT);
c=.5(Sum[S2[[j]]^2,{j,Length[PI2]}]-TT);
d=.5(TT+m^2-Sum[S1[[i]]^2,{i,Length[PI1]}]-Sum[S2[[j]]^2,{j,Length[PI2]}]);
];
{(a+d)/Binomial[m,2],a/(a+b+c),S,time[[1]]}
]

DIRECT Modul

Options[DIRECT]={Eps->0,MaxIter->20,DMin->10.^(-2)};
DIRECT[f_,dom_,OptionsPattern[]]:=Module[{n,ff,c,d,dmax,fc,dmaxrazni,pot0,pot1,pot2,temp,ind,novic,novid,novifc,red,TMin,iter=0,dtest,anim={}},
dtest=N[OptionValue[DMin]/Max[(dom[[All,2]]-dom[[All,1]])]];
DIRECT::maxiter="Dosegnut maksimalan broj iteracija `1`.";
n=Length[dom];
(* skaliranje funkcije na hiperkocku duljine stranice 1 *)
ff=Function[N[f@@({##}(dom[[All,2]]-dom[[All,1]])+dom[[All,1]]),50]];
(* početni centar i promjeri *)
c={Table[1/2,{i,n}]};
d={Table[1,{i,n}]};
dmax={Max[d[[1]]]};
dmaxrazni=Union[dmax];
fc={ff@@c[[1]]};
TMin={c[[1]],fc[[1]]};
(* GLAVNA PETLJA *)
Monitor[
While[(dmaxrazni[[1]]>dtest)&&(iter<OptionValue[MaxIter]),
iter++;
(* KORAK 1: traženje potencijalno optimalnih hiperkvadara *)
dmaxrazni=Union[dmax];
(* svi podaci grupirani po d *)
pot0=Table[Select[Range[Length[fc]],dmax[[#]]==dmaxrazni[[i]]&],{i,Length[dmaxrazni]}];
(* podaci s najmanjim f grupirani po d *)
pot1=Table[pot0[[i,Position[fc[[pot0[[i]]]],Min[fc[[pot0[[i]]]]]][[All,1]]]],{i,Length[dmaxrazni]}];
(* provjera ostalih uvjeta *)
pot2={pot1[[-1]]};
Do[temp=Min[Table[(fc[[pot1[[i,1]]]]-fc[[pot2[[j,1]]]])/(dmax[[pot1[[i,1]]]]-dmax[[pot2[[j,1]]]]),{j,Length[pot2]}]];
If[temp>OptionValue[Eps]&&(i==1||Max[Table[(fc[[pot1[[i,1]]]]-fc[[pot1[[j,1]]]])/(dmax[[pot1[[i,1]]]]-dmax[[pot1[[j,1]]]]),{j,i-1}]]<temp)
,AppendTo[pot2,pot1[[i]]]];
,{i,Length[dmaxrazni]-1,1,-1}];
pot2=Flatten[pot2,1];
(* KORAK 2: dijeljenje potencijalno optimalnih *)
Do[
(* smjerovi po kojima se dijeli *)
ind=Flatten[Position[d[[pot]],Max[d[[pot]]]]];
novic=Table[c[[pot]],{2*Length[ind]}];
novid=Table[d[[pot]],{2*Length[ind]}];
d[[pot,ind]]/=3;
dmax[[pot]]/=3;
Do[
novic[[2i-1,ind[[i]]]]-=d[[pot,ind[[i]]]];
novic[[2i,ind[[i]]]]+=d[[pot,ind[[i]]]];
,{i,Length[ind]}];
novifc=Table[ff@@novic[[i]],{i,Length[novic]}];
red=Ordering[Table[Min[{novifc[[2i-1]],novifc[[2i]]}],{i,Length[ind]}]];
Do[novid[[2red[[j]]-1,ind[[red[[i]]]]]]=novid[[2red[[j]],ind[[red[[i]]]]]]/=3,{i,1,Length[red]},{j,i,Length[red]}];
c=Join[c,novic];
d=Join[d,novid];
dmax=Join[dmax,Map[Max,novid]];
fc=Join[fc,novifc];
If[novifc[[2red[[1]]-1]]<TMin[[2]],TMin=N[{novic[[2red[[1]]-1]],novifc[[2red[[1]]-1]]},10]];
If[novifc[[2red[[1]]]]<TMin[[2]],TMin=N[{novic[[2red[[1]]]],novifc[[2red[[1]]]]},10]];
,{pot,pot2}];
(*Print[iter,": ",{TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]],TMin[[2]]}];*)
If[iter==OptionValue[MaxIter],Message[DIRECT::maxiter,OptionValue[MaxIter]];];
(* slika za 2D slučaj *)
(*If[n==2,
AppendTo[anim,
Graphics[{{Green,Rectangle[]},Table[{Blue,Rectangle[c[[i]]-d[[i]]/2+d[[i]]/10,c[[i]]+d[[i]]/2-d[[i]]/10]},{i,Length[d]}],Point[c]}]
];
];*)
];
,{iter,Length[pot2]}]
(* vraćanje u početnu domenu *)
TMin[[1]]=TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]];
(* ispis slika za 2D slučaj *)
(*If[n==2,Print[ListAnimate[anim,AnimationRunning->False]]];*)
Return[TMin];
]

Test-Example~9.6: 3 Points in R^2

SeedRandom[33]; 
m0=100; k=3; m=k m0; n=2;
CC={{1,1},{3.5,3.5},{6,1}};
slC=ListPlot[CC,PlotStyle->{Gray,Opacity[.8],PointSize[.025]}];
cov={{{1,-.3},{-.3,1.}},{{1,.3},{.3,1}},{{1,.7},{.7,1}}};
z0={{2,0},{4,2},{7,2}};
slz0=ListPlot[z0,PlotStyle->{Gray,Opacity[.7],PointSize[.025]}];
PI0=Table[RandomReal[MultinormalDistribution[CC[[j]],cov[[j]]],m0],{j,k}];
							(* Execution *)
A=Flatten[PI0,1];
w=Table[1,{i,m}]; 
xG={Min[A[[All,1]]]-.5,Max[A[[All,1]]]+.5};yG={Min[A[[All,2]]]-.5,Max[A[[All,2]]]+.5};
cc=Mean[WeightedData[A,w]];
slPI0=Show[Table[ListPlot[PI0[[j]],PlotStyle->{Hue[j/k],PointSize[.015]}],{j,k}],PlotRange->{xG,yG}];
sldata=Show[slPI0,slC,AxesOrigin->{0,0},GridLines->Automatic,AspectRatio->Automatic,ImageSize->Small]
slA=ListPlot[A,PlotStyle->{PointSize[.015]},AxesOrigin->{0,0},GridLines->Automatic,AspectRatio->Automatic,ImageSize->Small]
sl2=Show[slz0,Table[ListPlot[PI0[[j]],PlotStyle->{PointSize[.01]}],{j,k}]];
(* k-menas *)
sol1=Timing[MWKMeans[w,A,z0]];
Print["\nTime = ",sol1[[1]]]; sol=sol1[[2]];
Print[sol[[4]],": MWF = ",MWF[sol[[1]],sol[[2]],sol[[3]]],";  MWFmin = ",MWFmin[sol[[1]],sol[[2]],sol[[3]]],"; mj: ",Table[Length[sol[[1,j]]],{j,Length[PI]}]];

								(* FIGURE *)	
Print["Mahalanobis k-means RECONSTRUCTION: "]
slc=ListPlot[sol[[2]],PlotStyle->{Gray,PointSize[.02]}];
slPI=ListPlot[Table[sol[[1,j,All,2]],{j,k}],PlotStyle->{PointSize[.01],Table[Hue[j/k],{j,k}]}];
PI1=Table[Union[sol[[1,j,All,2]]],{j,Length[sol[[1]]]}];						
eig=Table[Eigensystem[sol[[3,j]]],{j,k}];
slvec=Table[0,{j,k}];
Do[
slvec[[j]]=Graphics[{Gray,{Arrowheads[Small],Thickness[.008],Arrow[{sol[[2,j]],sol[[2,j]]+eig[[j,1,1]] eig[[j,2,1]]}], Arrowheads[Small],Arrow[{sol[[2,j]],sol[[2,j]]+eig[[j,1,2]]eig[[j,2,2]]}]}}]
,{j,k}];
Do[
mj=Length[sol[[1,j]]];
tab=Sort[Table[dM[sol[[2,j]],sol[[1,j,s,2]],sol[[3,j]]],{s,mj}]];
rad=tab[[Ceiling[95mj/100]]];
slcon[j]=ContourPlot[dM[sol[[2,j]],{x,y},sol[[3,j]]]==rad,{x,xG[[1]],xG[[2]]},{y,yG[[1]],yG[[2]]}],{j,k}];
sl=Show[Table[slcon[j],{j,k}]];
sl3=Show[slPI,sl,slc,slvec,AspectRatio->Automatic,AxesOrigin->{0,0}];
Print[GraphicsGrid[{{sldata,sl3}},ImageSize->Medium]]
					(* COMPARISON *)
Print["COMPARISON: "]
Print["Initial Centers: ",CC]
Print["Centers: ",sol[[2]]]
Print["\nInitial Covariance Matrices:\n"]
Table[MatrixForm[cov[[j]]],{j,k}]
Print["\nCovariance Matrices:"]
Table[MatrixForm[sol[[3,j]]],{j,k}]
rnd=Rand[PI0,PI1];
Print["Confusion Matrix = ",rnd[[3]]//MatrixForm,"; Rand-Jacard index: ",{rnd[[1]],rnd[[2]]}]

Test-Example~9.7: 5 Points

SeedRandom[33]; 
m0=100; k=5; m=k m0; n=2;
CC={{2,2},{4,6},{6,10},{8,7},{9,3}};
slC=ListPlot[CC,PlotStyle->{Gray,Opacity[.8],PointSize[.025]}];
cov={{{1,0},{0,1.2}},{{2,-1},{-1,.6}},{{5,0},{0,.5}},{{1.2,0},{0,1}},{{2,1},{1,1}}};
z0={{1,1},{3,6},{5,9},{9,6},{9,3}};
slz0=ListPlot[z0,PlotStyle->{Gray,Opacity[.7],PointSize[.025]}];
PI0=Table[RandomReal[MultinormalDistribution[CC[[j]],cov[[j]]],m0],{j,k}];
							(* Execution *)
A=Flatten[PI0,1];
w=Table[1,{i,m}]; 
xG={Min[A[[All,1]]]-.5,Max[A[[All,1]]]+.5};yG={Min[A[[All,2]]]-.5,Max[A[[All,2]]]+.5};
cc=Mean[WeightedData[A,w]];
slPI0=Show[Table[ListPlot[PI0[[j]],PlotStyle->{Hue[j/k],PointSize[.015]}],{j,k}],PlotRange->{xG,yG}];
sldata=Show[slPI0,slC,AxesOrigin->{0,0},GridLines->Automatic,AspectRatio->Automatic,ImageSize->Small]
slA=ListPlot[A,PlotStyle->{PointSize[.015]},AxesOrigin->{0,0},GridLines->Automatic,AspectRatio->Automatic,ImageSize->Small]
sl2=Show[slz0,Table[ListPlot[PI0[[j]],PlotStyle->{PointSize[.01]}],{j,k}]];
(* k-menas *)
sol1=Timing[MWKMeans[w,A,z0]];
Print["\nTime = ",sol1[[1]]]; sol=sol1[[2]];
Print[sol[[4]],": MWF = ",MWF[sol[[1]],sol[[2]],sol[[3]]],";  MWFmin = ",MWFmin[sol[[1]],sol[[2]],sol[[3]]],"; mj: ",Table[Length[sol[[1,j]]],{j,Length[PI]}]];

								(* FIGURE *)	
Print["Mahalanobis k-means RECONSTRUCTION: "]
slc=ListPlot[sol[[2]],PlotStyle->{Gray,PointSize[.02]}];
slPI=ListPlot[Table[sol[[1,j,All,2]],{j,k}],PlotStyle->{PointSize[.01],Table[Hue[j/k],{j,k}]}];
PI1=Table[Union[sol[[1,j,All,2]]],{j,Length[sol[[1]]]}];						
eig=Table[Eigensystem[sol[[3,j]]],{j,k}];
slvec=Table[0,{j,k}];
Do[
slvec[[j]]=Graphics[{Gray,{Arrowheads[Small],Thickness[.008],Arrow[{sol[[2,j]],sol[[2,j]]+eig[[j,1,1]] eig[[j,2,1]]}], Arrowheads[Small],Arrow[{sol[[2,j]],sol[[2,j]]+eig[[j,1,2]]eig[[j,2,2]]}]}}]
,{j,k}];
Do[
mj=Length[sol[[1,j]]];
tab=Sort[Table[dM[sol[[2,j]],sol[[1,j,s,2]],sol[[3,j]]],{s,mj}]];
rad=tab[[Ceiling[95mj/100]]];
slcon[j]=ContourPlot[dM[sol[[2,j]],{x,y},sol[[3,j]]]==rad,{x,xG[[1]],xG[[2]]},{y,yG[[1]],yG[[2]]}],{j,k}];
sl=Show[Table[slcon[j],{j,k}]];
sl3=Show[slPI,sl,slc,slvec,AspectRatio->Automatic,AxesOrigin->{0,0}];
Print[GraphicsGrid[{{sldata,sl3}},ImageSize->Medium]]
					(* COMPARISON *)
Print["COMPARISON: "]
Print["Initial Centers: ",CC]
Print["Centers: ",sol[[2]]]
Print["\nInitial Covariance Matrices:\n"]
Table[MatrixForm[cov[[j]]],{j,k}]
Print["\nCovariance Matrices:"]
Table[MatrixForm[sol[[3,j]]],{j,k}]
rnd=Rand[PI0,PI1];
Print["Confusion Matrix = ",rnd[[3]]//MatrixForm,"; Rand-Jacard index: ",{rnd[[1]],rnd[[2]]}]

Test-Example~9.8: 3 very elongated clusters

SeedRandom[33]; 
m0=100; k=3; m=k m0; n=2;
CC={{4,5},{8,4},{9,6}};
slC=ListPlot[CC,PlotStyle->{Gray,Opacity[.8],PointSize[.025]}];
cov={{{.4,0},{0,8.}},{{2,1},{1,.6}},{{10,0},{0,.05}}};
z0={{5,5},{7,4},{10,9}};
slz0=ListPlot[z0,PlotStyle->{Gray,Opacity[.7],PointSize[.025]}];
PI0=Table[RandomReal[MultinormalDistribution[CC[[j]],cov[[j]]],m0],{j,k}];
							(* Execution *)
A=Flatten[PI0,1];
w=Table[1,{i,m}]; 
xG={Min[A[[All,1]]]-.5,Max[A[[All,1]]]+.5};yG={Min[A[[All,2]]]-.5,Max[A[[All,2]]]+.5};
cc=Mean[WeightedData[A,w]];
slPI0=Show[Table[ListPlot[PI0[[j]],PlotStyle->{Hue[j/k],PointSize[.015]}],{j,k}],PlotRange->{xG,yG}];
sldata=Show[slPI0,slC,AxesOrigin->{0,0},GridLines->Automatic,AspectRatio->Automatic,ImageSize->Small]
slA=ListPlot[A,PlotStyle->{PointSize[.015]},AxesOrigin->{0,0},GridLines->Automatic,AspectRatio->Automatic,ImageSize->Small]
sl2=Show[slz0,Table[ListPlot[PI0[[j]],PlotStyle->{PointSize[.01]}],{j,k}]];
(* k-menas *)
sol1=Timing[MWKMeans[w,A,z0]];
Print["\nTime = ",sol1[[1]]]; sol=sol1[[2]];
Print[sol[[4]],": MWF = ",MWF[sol[[1]],sol[[2]],sol[[3]]],";  MWFmin = ",MWFmin[sol[[1]],sol[[2]],sol[[3]]],"; mj: ",Table[Length[sol[[1,j]]],{j,Length[PI]}]];

								(* FIGURE *)	
Print["Mahalanobis k-means RECONSTRUCTION: "]
slc=ListPlot[sol[[2]],PlotStyle->{Gray,PointSize[.02]}];
slPI=ListPlot[Table[sol[[1,j,All,2]],{j,k}],PlotStyle->{PointSize[.01],Table[Hue[j/k],{j,k}]}];
PI1=Table[Union[sol[[1,j,All,2]]],{j,Length[sol[[1]]]}];						
eig=Table[Eigensystem[sol[[3,j]]],{j,k}];
slvec=Table[0,{j,k}];
Do[
slvec[[j]]=Graphics[{Gray,{Arrowheads[Small],Thickness[.008],Arrow[{sol[[2,j]],sol[[2,j]]+eig[[j,1,1]] eig[[j,2,1]]}], Arrowheads[Small],Arrow[{sol[[2,j]],sol[[2,j]]+eig[[j,1,2]]eig[[j,2,2]]}]}}]
,{j,k}];
Do[
mj=Length[sol[[1,j]]];
tab=Sort[Table[dM[sol[[2,j]],sol[[1,j,s,2]],sol[[3,j]]],{s,mj}]];
rad=tab[[Ceiling[95mj/100]]];
slcon[j]=ContourPlot[dM[sol[[2,j]],{x,y},sol[[3,j]]]==rad,{x,xG[[1]],xG[[2]]},{y,yG[[1]],yG[[2]]}],{j,k}];
sl=Show[Table[slcon[j],{j,k}]];
sl3=Show[slPI,sl,slc,slvec,AspectRatio->Automatic,AxesOrigin->{0,0}];
Print[GraphicsGrid[{{sldata,sl3}},ImageSize->Medium]]
					(* COMPARISON *)
Print["COMPARISON: "]
Print["Initial Centers: ",CC]
Print["Centers: ",sol[[2]]]
Print["\nInitial Covariance Matrices:\n"]
Table[MatrixForm[cov[[j]]],{j,k}]
Print["\nCovariance Matrices:"]
Table[MatrixForm[sol[[3,j]]],{j,k}]
rnd=Rand[PI0,PI1];
Print["Confusion Matrix = ",rnd[[3]]//MatrixForm,"; Rand-Jacard index: ",{rnd[[1]],rnd[[2]]}]

Test-Example~9.9: 3 Line Segments

SeedRandom[33]; 
m0=100; n=2;k=3; m=k m0;var=.01; PI0=Table[{},{j,k}];
sq=2Sqrt[2]/2.;
L={Line[{-{sq,sq},{sq,sq}}],Line[{{.5,2},{2,-2}}],Line[{{0,-1.5},{3,-1.5}}]};
z0={{.5,.5},{1.5,0},{0,-1}};
slz0=ListPlot[{z0},PlotStyle->{Orange,Opacity[.8],PointSize[.03]}];
slLines=Show[Graphics[{L}],Axes->True,AxesOrigin->{0,0},ImageSize->200];
Do[
Do[
al=Random[];
PI0[[j]]=Append[PI0[[j]],RandomReal[MultinormalDistribution[al*L[[j,1,1]]+(1-al)*L[[j,1,2]],var{{1,0},{0,1}}]]],
{s,m0}]
,{j,k}]
slPI0=ListPlot[Table[PI0[[j]],{j,k}],PlotStyle->Table[{PointSize[.015],Hue[j/k]},{j,k}]];
sldata=Show[slPI0, slLines,AxesOrigin->{0,0},AspectRatio->Automatic,GridLines->Automatic,ImageSize->Small]
sl2=Show[slz0,Table[ListPlot[PI0[[j]],PlotStyle->{PointSize[.01]}],{j,k}]];

								(* EXECUTION *)
A=Flatten[PI0,1];
slA=ListPlot[A,PlotStyle->{PointSize[.015]},AxesOrigin->{0,0},GridLines->Automatic,AspectRatio->Automatic,ImageSize->Small]
w=Table[1,{i,m}]; 
xG={Min[A[[All,1]]]-.5,Max[A[[All,1]]]+.5};yG={Min[A[[All,2]]]-.5,Max[A[[All,2]]]+.5};
cc=Total[A]/Length[A];
sol1=Timing[MWKMeans[w,A,z0]];
Print["\nTime = ",sol1[[1]]]; sol=sol1[[2]];
Print[sol[[4]],": MWF = ",MWF[sol[[1]],sol[[2]],sol[[3]]],";  MWFmin = ",MWFmin[sol[[1]],sol[[2]],sol[[3]]],"; mj: ",Table[Length[sol[[1,j]]],{j,Length[PI]}]];

								(* FIGURE *)
Print["Mahalanobis k-means RECONSTRUCTION: "]
slc=ListPlot[sol[[2]],PlotStyle->{Black,PointSize[.02]}];
slPI=ListPlot[Table[sol[[1,j,All,2]],{j,k}],PlotStyle->Table[Hue[j/k],{j,k}]];
PI1=Table[Union[sol[[1,j,All,2]]],{j,Length[sol[[1]]]}];
eig=Table[Eigensystem[sol[[3,j]]],{j,k}];
slvec=Table[0,{j,k}];
Do[
slvec[[j]]=Graphics[{Gray,{Arrowheads[Small],Thickness[.008],Arrow[{sol[[2,j]],sol[[2,j]]+eig[[j,1,1]] eig[[j,2,1]]}], Arrowheads[Small],Arrow[{sol[[2,j]],sol[[2,j]]+eig[[j,1,2]]eig[[j,2,2]]}]}}]
,{j,k}];
Do[
mj=Length[sol[[1,j]]];
tab=Sort[Table[dM[sol[[2,j]],sol[[1,j,s,2]],sol[[3,j]]],{s,mj}]];
rad=tab[[Ceiling[95mj/100]]];
slcon[j]=ContourPlot[dM[sol[[2,j]],{x,y},sol[[3,j]]]==rad,{x,xG[[1]],xG[[2]]},{y,yG[[1]],yG[[2]]}],{j,k}];
sl=Show[Table[slcon[j],{j,k}]];
sl3=Show[slPI,sl,slc,slvec,AspectRatio->Automatic,AxesOrigin->{0,0},PlotRange->{xG,yG}];
Print[GraphicsGrid[{{sldata,sl3}},ImageSize->500]]

						(* COMPARISON *)
Print["COMPARISON: "]
vec=Table[0,{j,k}];
Do[
u=v=L[[j,1,1,2]]-L[[j,1,2,2]];v=L[[j,1,1,1]]-L[[j,1,2,1]];
vec[[j]]=eig[[j,2,Ordering[eig[[j,1]],-1][[1]]]];
If[v==0,Print["Slope of ",j,"-th line: Ifinity",";   Recognized: ", vec[[j,2]]/vec[[j,1]]], Print["Slope of ",j,"-th line: ",u/v,";   Recognized: ", vec[[j,2]]/vec[[j,1]]]]
,{j,k}]	
PI1=Table[sol[[1,j,All,2]],{j,k}];
ind=Rand[PI0,PI1];
Print["Konfusion Matrix: ",ind[[1]]]
rnd=Rand[PI0,PI1];
Print["Confusion Matrix = ",rnd[[3]]//MatrixForm,"; Rand-Jacard index: ",{rnd[[1]],rnd[[2]]}]

Test-Example~9.10: 3 Line Segments

SeedRandom[33]; 
m0=100; n=2;k=3; m=k m0;var=.0025; PI0=Table[{},{j,k}];
sq=2Sqrt[2]/2.;
L={Line[{-{sq,sq},{sq,sq}}],Line[{{.5,2},{1.5,-1.5}}],Line[{{-1,-.5},{2,-1.}}]};
z0={{.5,.5},{1.5,0},{.5,-.5}};
slz0=ListPlot[{z0},PlotStyle->{Orange,Opacity[.8],PointSize[.03]}];
slLines=Show[Graphics[{L}],Axes->True,AxesOrigin->{0,0},ImageSize->200];
Do[
Do[
al=Random[];
PI0[[j]]=Append[PI0[[j]],RandomReal[MultinormalDistribution[al*L[[j,1,1]]+(1-al)*L[[j,1,2]],var{{1,0},{0,1}}]]],
{s,m0}]
,{j,k}]
slPI0=ListPlot[Table[PI0[[j]],{j,k}],PlotStyle->Table[{PointSize[.015],Hue[j/k]},{j,k}]];
sldata=Show[slPI0, slLines,AxesOrigin->{0,0},AspectRatio->Automatic,GridLines->Automatic,ImageSize->Small]
sl2=Show[slz0,Table[ListPlot[PI0[[j]],PlotStyle->{PointSize[.01]}],{j,k}]];

								(* EXECUTION *)
A=Flatten[PI0,1];
slA=ListPlot[A,PlotStyle->{PointSize[.015]},AxesOrigin->{0,0},GridLines->Automatic,AspectRatio->Automatic,ImageSize->Small]
w=Table[1,{i,m}]; 
xG={Min[A[[All,1]]]-.5,Max[A[[All,1]]]+.5};yG={Min[A[[All,2]]]-.5,Max[A[[All,2]]]+.5};
cc=Total[A]/Length[A];
sol1=Timing[MWKMeans[w,A,z0]];
Print["\nTime = ",sol1[[1]]]; sol=sol1[[2]];
Print[sol[[4]],": MWF = ",MWF[sol[[1]],sol[[2]],sol[[3]]],";  MWFmin = ",MWFmin[sol[[1]],sol[[2]],sol[[3]]],"; mj: ",Table[Length[sol[[1,j]]],{j,Length[PI]}]];

								(* FIGURE *)
Print["Mahalanobis k-means RECONSTRUCTION: "]
slc=ListPlot[sol[[2]],PlotStyle->{Black,PointSize[.02]}];
slPI=ListPlot[Table[sol[[1,j,All,2]],{j,k}],PlotStyle->Table[Hue[j/k],{j,k}]];
PI1=Table[Union[sol[[1,j,All,2]]],{j,Length[sol[[1]]]}];
eig=Table[Eigensystem[sol[[3,j]]],{j,k}];
slvec=Table[0,{j,k}];
Do[
slvec[[j]]=Graphics[{Gray,{Arrowheads[Small],Thickness[.008],Arrow[{sol[[2,j]],sol[[2,j]]+eig[[j,1,1]] eig[[j,2,1]]}], Arrowheads[Small],Arrow[{sol[[2,j]],sol[[2,j]]+eig[[j,1,2]]eig[[j,2,2]]}]}}]
,{j,k}];
Do[
mj=Length[sol[[1,j]]];
tab=Sort[Table[dM[sol[[2,j]],sol[[1,j,s,2]],sol[[3,j]]],{s,mj}]];
rad=tab[[Ceiling[95mj/100]]];
slcon[j]=ContourPlot[dM[sol[[2,j]],{x,y},sol[[3,j]]]==rad,{x,xG[[1]],xG[[2]]},{y,yG[[1]],yG[[2]]}],{j,k}];
sl=Show[Table[slcon[j],{j,k}]];
sl3=Show[slPI,sl,slc,slvec,AspectRatio->Automatic,AxesOrigin->{0,0},PlotRange->{xG,yG}];
Print[GraphicsGrid[{{sldata,sl3}},ImageSize->Medium]]

						(* COMPARISON *)
Print["COMPARISON: "]
vec=Table[0,{j,k}];
Do[
u=v=L[[j,1,1,2]]-L[[j,1,2,2]];v=L[[j,1,1,1]]-L[[j,1,2,1]];
vec[[j]]=eig[[j,2,Ordering[eig[[j,1]],-1][[1]]]];
If[v==0,Print["Slope of ",j,"-th line: Ifinity",";   Recognized: ", vec[[j,2]]/vec[[j,1]]], Print["Slope of ",j,"-th line: ",u/v,";   Recognized: ", vec[[j,2]]/vec[[j,1]]]]
,{j,k}]	
PI1=Table[sol[[1,j,All,2]],{j,k}];
rnd=Rand[PI0,PI1];
Print["Confusion Matrix = ",rnd[[3]]//MatrixForm,"; Rand-Jacard index: ",{rnd[[1]],rnd[[2]]}]

Example~6.29: 2 Points + 3 Line Segments

SeedRandom[33]; 
							(* Tocke *)
n=2;m0=100; CC={{3,2},{8,6}}; kC=Length[CC];
slCC=ListPlot[CC,PlotStyle->{Red,Opacity[.7],PointSize[.025]}];
cov={.5{{2,-1},{-1,.9}},.3{{2,1},{1,1}}};
podCC=Table[RandomReal[MultinormalDistribution[CC[[j]],cov[[j]]],m0],{j,kC}];
							(* Segmenti *)
var=.1;
L={Line[{{1,9},{6,9}}],Line[{{3,5},{5,9}}],Line[{{3,6},{7,2}}]};
kL=Length[L];
slL=Show[Graphics[{L}],Axes->True,AxesOrigin->{0,0},ImageSize->200];
podL=Table[{},{j,kL}];
Do[
Do[
al=Random[];
podL[[j]]=Append[podL[[j]],RandomReal[MultinormalDistribution[al*L[[j,1,1]]+(1-al)*L[[j,1,2]],var{{1,0},{0,1}}]]],
{s,m0}]
,{j,kL}]
PI0={podCC[[1]],podCC[[2]],podL[[1]],podL[[2]],podL[[3]]};
k=k0=Length[PI0];
z0={{2,2},{9,5},{3,9},{4,7},{5,4}};
slz0=ListPlot[{z0},PlotStyle->{Orange,Opacity[.8],PointSize[.03]}];
slPI0=Show[Table[ListPlot[PI0[[j]],PlotStyle->{Hue[j/k],PointSize[.01]}],{j,k}]];
sldata=Show[slL,slCC,slPI0,AxesOrigin->{0,0},AspectRatio->Automatic];
sl2=Show[slz0,Table[ListPlot[PI0[[j]],PlotStyle->{PointSize[.01]}],{j,k}]];

								(* EXECUTION *)
A=Flatten[PI0,1];
m=Length[A];
w=Table[1,{i,m}];
xG={Min[A[[All,1]]]-.5,Max[A[[All,1]]]+.5};yG={Min[A[[All,2]]]-.5,Max[A[[All,2]]]+.5};
cc=Total[A]/Length[A];
sol1=Timing[MWKMeans[w,A,z0]];
Print["\nTime = ",sol1[[1]]]; sol=sol1[[2]];
Print[sol[[4]],": MWF = ",MWF[sol[[1]],sol[[2]],sol[[3]]],";  MWFmin = ",MWFmin[sol[[1]],sol[[2]],sol[[3]]],"; mj: ",Table[Length[sol[[1,j]]],{j,Length[PI]}]];

								(* FIGURE *)
Print["Mahalanobis k-means RECONSTRUCTION: "]
slc=ListPlot[sol[[2]],PlotStyle->{Black,PointSize[.02]}];
slPI=ListPlot[Table[sol[[1,j,All,2]],{j,k}],PlotStyle->Table[Hue[j/k],{j,k}]];
PI1=Table[Union[sol[[1,j,All,2]]],{j,Length[sol[[1]]]}];
(*
eig=Table[Eigensystem[sol[[3,j]]],{j,k}];
slvec=Table[0,{j,k}];
Do[
slvec[[j]]=Graphics[{Brown,{Arrowheads[Small],Arrow[{sol[[2,j]],sol[[2,j]]+eig[[j,1,1]] eig[[j,2,1]]}], Arrowheads[Small],Arrow[{sol[[2,j]],sol[[2,j]]+eig[[j,1,2]]eig[[j,2,2]]}]}}]
,{j,k}]; *)
Do[
mj=Length[sol[[1,j]]];
tab=Sort[Table[dM[sol[[2,j]],sol[[1,j,s,2]],sol[[3,j]]],{s,mj}]];
rad=tab[[Ceiling[95mj/100]]];
slcon[j]=ContourPlot[dM[sol[[2,j]],{x,y},sol[[3,j]]]==rad,{x,xG[[1]],xG[[2]]},{y,yG[[1]],yG[[2]]}],{j,k}];
sl=Show[Table[slcon[j],{j,k}]];
sl3=Show[slPI,sl,slc,AspectRatio->Automatic,AxesOrigin->{0,0}];
Print[GraphicsGrid[{{sldata,sl3}},ImageSize->500]]

						(* COMPARISON *)
Print["COMPARISON: "]
Print["Initial Centers: ",CC]
Print["Centers: ",Table[sol[[2,j]],{j,kC}]]
Print["\nInitial Covariance Matrices:\n"]
Table[MatrixForm[cov[[j]]],{j,kC}]
Print["\nCovariance Matrices:"]
Table[MatrixForm[sol[[3,j]]],{j,kC}]
Print["COMPARISON: "]
eig=Table[Eigensystem[sol[[3,j]]],{j,kC+1,k}];
vec=Table[0,{j,kL}];
Do[
u=v=L[[j,1,1,2]]-L[[j,1,2,2]];v=L[[j,1,1,1]]-L[[j,1,2,1]];
vec[[j]]=eig[[j,2,Ordering[eig[[j,1]],-1][[1]]]];
If[v==0,Print["Slope of ",j,"-th line: Ifinity",";   Recognized: ", vec[[j,2]]/vec[[j,1]]], Print["Slope of ",j,"-th line: ",u/v,";   Recognized: ", vec[[j,2]]/vec[[j,1]]]]
,{j,kL}]	
Print["Centers: ",sol[[2]]]
PIP=Table[sol[[1,j,All,2]],{j,k}];
rnd=Rand[PI0,PIP];
Print["Confusion Matrix = ",rnd[[3]]//MatrixForm,"; Rand-Jacard index: ",{rnd[[1]],rnd[[2]]}]

Example~6.29 - Figure~6.8: 2 Points + 3 Line Segments

z0={{3,3},{7,5}}; 
z0={{4,2},{4,8},{7,5}};
z0={{3,3},{3,7},{7,3},{7,5}};
z0={{2,2},{9,5},{3,9},{4,7},{5,4}};
z0={{2,2},{2,7},{5,3},{5,7},{7.5,3},{6.5,7.5}};

k=Length[z0];
slcen=ListPlot[z0,PlotStyle->{Orange,PointSize[.05],Opacity[.9]}];
slVoronoi=DiagramPlot[z0,LabelPoints->False];
slIni=Show[ListPlot[A,PlotStyle->{PointSize[.015]},AspectRatio->Automatic,ImageSize->Small],slVoronoi,slcen,GridLines->Automatic,ImageSize->Small]
sol=MWKMeans[w,A,z0];
slc=ListPlot[sol[[2]],PlotStyle->{Gray,PointSize[.05],Opacity[.7]}];
slPI=ListPlot[Table[sol[[1,j,All,2]],{j,k}],PlotStyle->{Table[PointSize[.015],{j,k}]},PlotStyle->{Table[Hue[j/k],{j,k}]}];
Do[
mj=Length[sol[[1,j]]];
tab=Sort[Table[dM[sol[[2,j]],sol[[1,j,s,2]],sol[[3,j]]],{s,mj}]];
rad=tab[[Ceiling[95mj/100]]];
slcon[j]=ContourPlot[dM[sol[[2,j]],{x,y},sol[[3,j]]]==rad,{x,xG[[1]],xG[[2]]},{y,yG[[1]],yG[[2]]}],{j,k}];
sl=Show[Table[slcon[j],{j,k}],ImageSize->Small];
slEnd=Show[slPI,slc,sl,AspectRatio->Automatic,GridLines->Automatic,ImageSize->Small];
Print[sol[[4]],": MWF = ",MWF[sol[[1]],sol[[2]],sol[[3]]],";  MWFmin = ",MWFmin[sol[[1]],sol[[2]],sol[[3]]],"; mj: ",Table[Length[sol[[1,j]]],{j,Length[PI]}]];
GraphicsGrid[{{slIni,slEnd}}]

Test-Example: Points in R^3

SeedRandom[13]; 
n=3;m0=100; k=5; m=k m0; PI0=Table[{},{j,k}];par=10;var=3;
(* par odredjuje raspon eigenvalues of covariance matrices *)
CC=RandomInteger[{0,10},{k,n}];
mpp=.2;While[mpp<.5,CC=RandomInteger[{0,10},{k,n}];pp=Table[{},{j,k}];
Do[
pp[[j]]=Min[Drop[Table[Norm[CC[[j]]-CC[[s]]],{s,k}],{j}]];
,{j,k}];
mpp=Min[pp]]
Print["Centri su OK!"]
u=Table[i,{i,n}];u=u/Norm[u];
U=IdentityMatrix[n]-2KroneckerProduct[u,u];
eig=Table[RandomReal[{0,par},n],{j,k}];
(*
Do[
eps=par/(10k);
seig=Sort[eig[[j]]];
While[Min[Table[seig[[s+1]]-seig[[s]],{s,k-1}]]<eps,eig[[j]]=RandomReal[{0,par},n];seig=Sort[eig[[j]]]  
        ]
,{j,k}]
*)
cov=var Round[Table[U.DiagonalMatrix[eig[[j]]].Transpose[U],{j,k}],.01];
Do[
PI0[[j]]=RandomReal[MultinormalDistribution[CC[[j]],cov[[j]]],m0];
,{j,k}];
w=Table[1,{i,m}]; W=Total[w];
z0=Table[CC[[j]]+RandomReal[{-1,1},n],{j,k}];
Print["Podaci su OK!"]
					(* Execution *)
A=Flatten[PI0,1];
w=Table[1,{i,m}];W=Total[w] ;
xG={Min[A[[All,1]]],Max[A[[All,1]]]};
yG={Min[A[[All,2]]],Max[A[[All,2]]]};
sol1=Timing[MWKMeans[w,A,z0]];
Print["\nTime = ",sol1[[1]]]; sol=sol1[[2]];
Print[sol[[4]],": MWF = ",MWF[sol[[1]],sol[[2]],sol[[3]]],";  MWFmin = ",MWFmin[sol[[1]],sol[[2]],sol[[3]]],"; mj: ",Table[Length[sol[[1,j]]],{j,Length[PI]}]];

					(* COMPARISON *)
Print["COMPARISON: "]
Print["Initial Centers: ",CC]
Print["Centers: ",CC1=sol[[2]]]
Print["Hausdorff= ",Hdist[CC,CC1]]
Print["\nEigenvalues of Initial Covariance Matrices:"]
Table[Eigenvalues[cov[[j]]],{j,k}]
Print["\nEigenvalues of Estimated Covariance Matrices:"]
Table[Eigenvalues[sol[[3,j]]],{j,k}]
PI1=Table[sol[[1,j,All,2]],{j,k}];
rnd=Rand[PI0,PI1];
Print["Confusion Matrix = ",rnd[[3]]//MatrixForm,"; Rand-Jacard index: ",{rnd[[1]],rnd[[2]]}]

Test primjer Earthquake in R^3 (HR)
Data: {Latitude, Longitude, Depth}, Weights: Magnitudes

SetDirectory[NotebookDirectory[]]; 
asd=ReadList["C:/Users/Rudolf Scitovski/Documents/AB-STRUKA/1_RADOVI/B1_CLUSTERS/B1_PACO/Paco-2019/DATA/HR35.txt"];
HR0=ReadList["C:/Users/Rudolf Scitovski/Documents/AB-STRUKA/1_RADOVI/B1_CLUSTERS/B1_PACO/Paco-2019/DATA/HR0.txt"];
m=Length[asd]
slHR=ListPlot[HR0,Joined->True,PlotStyle->{Orange}];
A=Table[{asd[[i,8]],asd[[i,7]],asd[[i,10]]},{i,m}];
w=Table[1,{i,m}]; W=Total[w];
w=Table[asd[[i,10]],{i,m}]; W=Total[w];
xG={Min[A[[All,1]]],Max[A[[All,1]]]};yG={Min[A[[All,2]]],Max[A[[All,2]]]};
zG={Min[A[[All,3]]],Max[A[[All,3]]]};
cc=Total[A]/Length[A];
Pod=Table[{asd[[i,8]],asd[[i,7]]},{i,m}];
pod4=pod5=pod6={};
Do[
If[w[[i]]≤4,pod4=Append[pod4,Pod[[i]]],
If[w[[i]]≤5,pod5=Append[pod5,Pod[[i]]],pod6=Append[pod6,Pod[[i]]]]]
,{i,m}]
slpod=Show[slHR,ListPlot[pod4,PlotStyle->{Black,Opacity[.4], PointSize[.003]},Axes->True,AspectRatio->1],ListPlot[pod5,PlotStyle->{Black, Opacity[.7],PointSize[.006]},Axes->True,AspectRatio->1],ListPlot[pod6,PlotStyle->{Black, Opacity[1],PointSize[.01]},Axes->True,AspectRatio->1],ImageSize->300];
sldepth=ListPlot3D[A];
GraphicsGrid[{{slpod,sldepth}},ImageSize->600]

Call Mahalanobis weighted k - means

(* Koristi se dM *)
sol1=Timing[MWKMeans[w,A,z0,1]];
Print[Style["\nTime = ",Blue,18],sol1[[1]]]; sol=sol1[[2]];
Print[sol[[4]],": MWF = ",MWF[sol[[1]],sol[[2]],sol[[3]]],";  MWFmin = ",MWFmin[sol[[1]],sol[[2]],sol[[3]]],"; mj: ",Table[Length[sol[[1,j]]],{j,Length[PI]}]];

								(* FIGURE *)	
Print["Mahalanobis k-means RECONSTRUCTION: "]
slc=ListPlot[sol[[2]],PlotStyle->{Gray,PointSize[.02]}];
slPI=ListPlot[Table[sol[[1,j,All,2]],{j,k}],PlotStyle->{PointSize[.01],Table[Hue[j/k],{j,k}]}];
PI1=Table[Union[sol[[1,j,All,2]]],{j,Length[sol[[1]]]}];						
eig=Table[Eigensystem[sol[[3,j]]],{j,k}];
slvec=Table[0,{j,k}];
Do[
slvec[[j]]=Graphics[{Gray,{Arrowheads[Small],Thickness[.008],Arrow[{sol[[2,j]],sol[[2,j]]+eig[[j,1,1]] eig[[j,2,1]]}], Arrowheads[Small],Arrow[{sol[[2,j]],sol[[2,j]]+eig[[j,1,2]]eig[[j,2,2]]}]}}]
,{j,k}];
Do[
mj=Length[sol[[1,j]]];
tab=Sort[Table[dM[sol[[2,j]],sol[[1,j,s,2]],sol[[3,j]]],{s,mj}]];
rad=tab[[Ceiling[95mj/100]]];
slcon[j]=ContourPlot[dM[sol[[2,j]],{x,y},sol[[3,j]]]==rad,{x,xG[[1]],xG[[2]]},{y,yG[[1]],yG[[2]]}],{j,k}];
sl=Show[Table[slcon[j],{j,k}]];
sl3=Show[slPI,sl,slc,AxesOrigin->{0,0},Ticks->{{-.5,.5},{-.5,.5}},ImageSize->Small]
Print[GraphicsGrid[{{sldata,sl3}},ImageSize->Medium]]
					(* COMPARISON *)
Print["COMPARISON: "]
Print["Initial Centers: ",CC]
Print["\nCovariance Matrices:"]
Table[MatrixForm[sol[[3,j]]],{j,k}]
rnd=Rand[PI0,PI1];
Print["Confusion Matrix = ",rnd[[3]]//MatrixForm,"; Rand-Jacard index: ",{rnd[[1]],rnd[[2]]}]

Call EM algorithm

(* Koristi se dm *)
time=Timing[em=EM[w,A,z0,0];];
Print[em[[4]],". ITERACIJA;    time = ",time[[1]],"\nCenters: ",em[[2]],"\nCov: ",Table[em[[3,j]]//MatrixForm,{j,Length[em[[2]]]}],"\nP: ",em[[5]],"\nHdist= ",Hdist[z0,em[[2]]]];
rand=Rand[PI0,em[[1]]];
Print["Confusion Matrix=",rand[[1]]//MatrixForm,"; Rand-Jacard index:",rand[[2]]]
(* Figures *)
k=Length[em[[2]]];
slPI=ListPlot[em[[1]],PlotStyle->Table[Hue[j/k],{j,k}]];
sl0=ListPlot[em[[2]],PlotStyle->{Brown,Opacity[.8],PointSize[.025]}];
Do[
rad=Max[Table[dm[em[[2,j]],em[[1,j,s]],em[[3,j]]],{s,Length[em[[1,j]]]}]];
slcon[j]=ContourPlot[dm[em[[2,j]],{x,y},em[[3,j]]]==75rad/100,{x,xG[[1]],xG[[2]]},{y,yG[[1]],yG[[2]]}],{j,k}];
sl=Show[Table[slcon[j],{j,k}]];
sl2=Show[slPI,sl,sl0,Ticks->{{-.5,.5},{-.5,.5}},AxesOrigin->{0,0},ImageSize->Small];
Print[sl2];
GraphicsGrid[{{sldata,sl2}},ImageSize->Medium]

Call Mahalanobis Incremental algorithm

(* Koristi se dM *)
cen={Mean[A]};KK=6;
inc=MInc[w,A,cen,KK,1];
(* Spremanje particija u.txt datoteku *)
SetDirectory[NotebookDirectory[]]; 
Export["PART.txt",PART]