Test-Example~9.7:  5 points in R^2

Submodules: MWKMeans[], MWF[], MWFmin[], Rand[], Hdist[], dM[], DIRECT[];
Output: data set, optimal partition, Rand and Jaccard index, confusion matrix;

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]]}]