Test-Example~9.10: 3 Line Segments

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

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