Spherical Clustering: WKMeans and Incremental

(* 
Spherical Clustering is a complete Mathematica-program for spherical
clustering using the LS 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 WKMeans[] algorithm.Then one can call
the incremental algorithm Inc[] which will,among partitions obtained and
based on several indexes,propose the MAPart.If the original partition is
known,it will also calculate the Rand and Jaccard indexes,as well as the
Hausdorff distance between sets of original and calculated cluster centers. 
*)

Modules

Needs["ComputationalGeometry`"]

Objective functions

WF[PI_]:=Module[{m=Length[Flatten[PI,1]],c},
c=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}];
Return[Sum[PI[[j,s,1]]*Norm[PI[[j,s,2]]-c[[j]]]^2,{j,Length[PI]} ,{s,Length[PI[[j]]]}] ]
];

WG[PI_]:=Module[{m=Length[Flatten[PI,1]],c,set,cc},
c=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}];
set=Flatten[PI,1];
cc=Mean[WeightedData[set[[All,2]],set[[All,1]]]];
Return[Sum[Length[PI[[j]]]*Norm[cc-c[[j]]]^2,{j,Length[PI]} ]]
];

WFmin[PI_]:=Module[{A=Flatten[PI,1][[All,2]],k=Length[PI]},
c=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}];
Sum[Flatten[PI,1][[All,1]][[i]]Min[Table[Norm[c[[j]]-A[[i]]]^2,{j,k}]],{i,Length[A]}]
]

Weighted k - means algorithm

WKMeans[w_,A_,z_,Ind_:0]:=Module[{m=Length[A],k=Length[z],c0,Pod,imin,PI,It=0},
Pod=Table[{w[[i]],A[[i]]},{i,m}];
c0=z; 
PI=Table[{},{i,k}]; 
			(* Minimal distance principle *)
Do[
tab=Table[Norm[Pod[[i,2]]-c0[[j]]],{j,k}];
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI];
If[Ind≠ 0,Print[It,": F = ",WF[PI]]];
						(* New centers *)
While[
c=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}]; 
		(* The case of reducing the number of clusters in a partition
				in the case n=2 - can be extended to an arbitrary n *)
If[Length[c]<Length[c0],
Print["There was a case of reducing the number of clusters in the partition {",Length[c0]," -> ",Length[c],"}\nTry with a new initial centers"];
	];
Chop[Norm[c-c0]]≠0,c0=c;It=It+1;
				(* New partition *)
PI=Table[{},{i,k}]; 
Do[
tab=Table[Norm[Pod[[i,2]]-c0[[j]]],{j,k}];
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI]; 
If[Ind≠ 0,Print[It,": F = ",WF[PI]]];
               ];
{PI,c0,WF[PI],It}]

Incremental algorithm

Inc[w_,A_,cen_,K_,Ind_:0]:=Module[{k,centri,centriN,m=Length[A],n=Length[A[[1]]],dom},
AG=Table[{Min[A[[All,i]]]-.5,Max[A[[All,i]]]+.5},{i,n}];
					(* Initialization *)
centri=cen;
PART={};IndCH=IndDB=IndSSWC=IndDunn=IndH={};
				(* Iterative procedure *)
Do[
Print["\nTo existing ",Length[centri]," centers one new is added"];
	(* Compile in the case n=2 - can be extended to an arbitrary n *)
Clear[u,v];
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}];	
						(* Non Compile 
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}];
*)
					(* Applying of the DIRECT-algorithm *)
dom=AG;
minp=DIRECT[F2,dom,DMin->10.^(-3),MaxIter->50, Eps->10^(-4)]//Timing; min=minp[[2]];(* Print["DIRECT: ",min]; *)
centriN=Append[centri,min[[1]]];
If[Ind≠0,Print["Novi optimirani centar: ",centriN]];
If[Ind≠0,slcen=ListPlot[centriN,PlotMarkers->{"✶",20},PlotStyle->{Gray,PointSize[.025],Opacity[.7]}]];
					(* Call WKMeans *)
t1=Timing[sol=WKMeans[w,A,centriN]];
k=Length[sol[[1]]]; 
If[Ind≠0,Print["Novi WKMeans profinjeni centar: ",sol[[2]]]];
If[Ind≠0,Print["WKMeans: {WF,WG}= ",{WF[sol[[1]]],WG[sol[[1]]]},"; It=",sol[[4]],"; mj: ",Table[Length[sol[[1,j]]],{j,Length[sol[[1]]]}]]];
If[Ind≠0,Print["Time (DIRECT+WKM) = ",minp[[1]],"+",t1[[1]]," = ",minp[[1]]+t1[[1]]]];
vch=VCH[sol[[1]]];IndCH=Append[IndCH,vch];
vdb=VDB[sol[[1]]];IndDB=Append[IndDB,vdb];
svc=SSWC[sol[[1]]];IndSSWC=Append[IndSSWC,svc];
vd=VD[sol[[1]]];IndDunn=Append[IndDunn,vd];
If[Ind≠0,Print["k=",Length[centriN],"; CH=",vch,"; DB=",vdb,"; SSWC=",svc,"; Dunn=",vd]];
If[Ind≠0,
tab=Table[sol[[1,j,All,2]],{j,Length[sol[[1]]]}];
ind=Rand[podaci,tab];
Print["Hausdorff dist. = ",Hdist[sol[[2]],c]];
Print["{Rand, Jaccard} index = ",{ind[[1]],ind[[2]]}];
Print["Confusion matrix S= ",ind[[3]]//MatrixForm];
];
                                                                     (* FIGURES *)
If[n==2,
k=Length[sol[[1]]];
slc=ListPlot[sol[[2]],PlotMarkers->{"✶",20},PlotStyle->{Gray,PointSize[.025],Opacity[.7]}];
slPI=ListPlot[Table[sol[[1,j,All,2]],{j,k}],PlotStyle->Table[{PointSize[.01],Hue[j/k]},{j,k}]]; 
	(* This is just for drawing the next center in the case n=2
 - can be extended to an arbitrary n *)
						(* Compile *)
Clear[u,v,F2]; 
FCompile=Compile[{{u,_Real},{v,_Real}},
Evaluate[Sum[Min[Append[Table[Norm[A[[i]]-centriN[[s]]]^2,{s,Length[centriN]}],Norm[A[[i]]-{u,v}]^2]],{i,m}] ]];	
F2[u_,v_]:=Apply[FCompile,{u,v}];
mini=DIRECT[F2,dom,DMin->10.^(-3),MaxIter->20, Eps->10^(-4)];
nextc=ListPlot[{mini[[1]]},PlotStyle->{Red,PointSize[.04],Opacity[.8]}];
								(*  *)
sl2=Print[Show[slPI,slc,nextc,AspectRatio->Automatic,ImageSize->Small]]
];			      
centri=sol[[2]];     
PART=Append[PART,sol[[1]]]; 
,{brk,K}];
If[Ind≠0,Print["IndCH: ",IndCH,"\nIndDB: ",IndDB,"\nIndSSWC: ",IndSSWC,"\nIndDunn: ",IndDunn]];
If[Ind≠0,Print[GraphicsGrid[{{ListPlot[Table[{i+1,IndCH[[i]]},{i,Length[IndCH]}],Joined->True,AxesOrigin->{0,1000}],ListPlot[Table[{i+1,IndDB[[i]]},{i,Length[IndDB]}],Joined->True,AxesOrigin->{0,.5}],ListPlot[Table[{i+1,IndSSWC[[i]]},{i,Length[IndSSWC]}],Joined->True,AxesOrigin->{0,.75}],ListPlot[Table[{i+1,IndDunn[[i]]},{i,Length[IndDunn]}],Joined->True,AxesOrigin->{0,0}]}},ImageSize->Large]]];
{PART,IndCH,IndDB,IndSSWC,IndDunn}
]

Hausdorff distance between two sets of centroids

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

Calinski - Harabasz index

Clear[VCH]
VCH[PI_]:=Module[{k=Length[PI],m,FF,mj,F,G},
FF=Table[0,{j,Length[PI]}];
m=Sum[Length[PI[[j]]],{j,Length[PI]}];
z=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}];
set=Flatten[PI,1];
cc=Mean[WeightedData[set[[All,2]],set[[All,1]]]];
Do[
mj=Length[PI[[j]]];
FF[[j]]=Sum[PI[[j,s,1]]*Norm[PI[[j,s,2]]-z[[j]] ]^2,{s,mj}];
   ,{j,Length[PI]}];
F=Total[FF];
G=Sum[Sum[PI[[j,s,1]],{s,Length[PI[[j]]]}]*Norm[z[[j]]-cc]^2,{j,Length[PI]}];
(m-Length[PI])*G/((Length[PI]-1)F)
                                     ]

Davies - Bouldin

Clear[VDB]
VDB[PI_]:=Module[{FF,VV,mj,sum=0},
FF=VV=Table[0,{j,Length[PI]}];
z=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}];
Do[
mj=Length[PI[[j]]];
FF[[j]]=Sum[PI[[j,i,1]]*Norm[PI[[j,i,2]]-z[[j]]]^2,{i,mj}];
VV[[j]]=Sqrt[FF[[j]]/Sum[PI[[j,i,1]],{i,mj}]]
   ,{j,Length[PI]}];
Do[
max={};
Do[
If[i≠j,max=Append[max,(VV[[i]]+VV[[j]])/Norm[z[[i]]-z[[j]]]]]
,{i,Length[PI]}];
sum=sum+Max[max]
,{j,Length[PI]}];
sum/Length[PI]
            ]

Simplified Silhouette With Criterion

Clear[SSWC]
SSWC[PI_]:=Module[{k=Length[PI],a,b,dd,w},
set=Flatten[PI,1];m=Length[set];
a=b=Table[0,{i,m}];
w=set[[All,1]];A=set[[All,2]];
z=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}];
Do[
dd=Sort[Table[w[[i]] Norm[A[[i]]-z[[j]]]^2,{j,k}]];
a[[i]]=dd[[1]];b[[i]]=dd[[2]];
,{i,m}];
Sum[(b[[i]]-a[[i]])/Max[{a[[i]],b[[i]]}],{i,m}]/m
]

Dunn Index

Clear[VD]
VD[PI_]:=Module[{z,DD,diam},
z=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}];
DD=Table[If[i≠j,Min[Table[Norm[PI[[i,p,2]]-PI[[j,q,2]]],{p,Length[PI[[i]]]},{q,Length[PI[[j]]]}]],Infinity],{i,Length[PI]},{j,Length[PI]}];
diam=Table[Max[Table[Norm[PI[[i,p,2]]-PI[[i,q,2]]],{p,Length[PI[[i]]]},{q,Length[PI[[i]]]}]],{i,Length[PI]}];
Min[DD/Max[diam]]
]

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

Example-2.25:  5 clusters in R^2

SeedRandom[1213];
c={{2,2},{3,5},{6,7},{7,3},{8,8}};
var=.5; broj=100; AG={{0,10},{0,10}};
podaci=Table[RandomVariate[MultinormalDistribution[ci,var IdentityMatrix[2]],broj],{ci,c}];
pod=Flatten[podaci,1]; 
m=Length[pod];
w=Table[1,{i,m}];
						(* PODACI *)
A=Flatten[podaci,1];
Print["Dimensions: ",Dimensions[pod], "Mean[A] = ",Mean[A]];	
xG={Min[A[[All,1]]]-.5,Max[A[[All,1]]]+.5};yG={Min[A[[All,2]]]-.5,Max[A[[All,2]]]+.5};
AG={xG,yG};
m=Length[A];
w=Table[1,{i,m}]; 
slA=ListPlot[A,PlotStyle->{PointSize[.01]},AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG,ImageSize->Small];
slc=ListPlot[c,PlotStyle->{Gray,PointSize[.03],Opacity[.6]}];
slpod1=Show[ListPlot[podaci,PlotStyle->{PointSize[.01]}],slc,AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG,ImageSize->Small]
z0={{3,3},{4,5},{6,6},{8,4},{8,6}};k=Length[z0];
slz0=ListPlot[z0,PlotStyle->{Gray,PointSize[.04],Opacity[.8]}];
slpod2=Show[slA,slz0,AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG]
slpod3=Show[slA,ListPlot[{Mean[A]},PlotMarkers->{"✶",20},PlotStyle->{Gray,PointSize[.025],Opacity[.7]}],AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG]
GraphicsGrid[{{slpod1,slpod2}},ImageSize->Medium]
							(* k-means *)
time=Timing[sol=WKMeans[w,A,z0]];
Print["F= ",WF[sol[[1]]],   ";  G= ",WG[sol[[1]]],"; CPU=",time[[1]], "; It = ",sol[[4]]]
							(* SLIKE *)
Print["k-means reconstruction for selected centers z0: ",z0]
slc=ListPlot[sol[[2]],PlotMarkers->{"✶",20},PlotStyle->{Gray,PointSize[.025],Opacity[.7]}];
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]]]}];
sl3=Show[slA,slc,AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG];
Print[GraphicsGrid[{{slpod1,sl3}},ImageSize->Medium]];
tab=Table[sol[[1,j,All,2]],{j,Length[sol[[1]]]}];
ind=Rand[podaci,tab];
Print["Hausdorff dist. = ",Hdist[sol[[2]],c]]
Print["{Rand, Jaccard} index = ",{ind[[1]],ind[[2]]}]
Print["Confusion matrix S= ",ind[[3]]//MatrixForm]

Test-example~9.2:  5 clusters in R^2

SeedRandom[1213];
c={{2.5,3},{2,8.5},{5,6},{8,5},{8.5,1.5}};
var=1; broj=200; AG={{0,10},{0,10}};
podaci=Table[RandomVariate[MultinormalDistribution[ci,var IdentityMatrix[2]],broj],{ci,c}];
pod=Flatten[podaci,1]; 
m=Length[pod];
w=Table[1,{i,m}];
Dimensions[pod]	
						(* PODACI *)
A=Flatten[podaci,1];
xG={Min[A[[All,1]]]-.5,Max[A[[All,1]]]+.5};yG={Min[A[[All,2]]]-.5,Max[A[[All,2]]]+.5};
AG={xG,yG};
m=Length[A];
w=Table[1,{i,m}]; 
slA=ListPlot[A,PlotStyle->{PointSize[.008]},AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG,ImageSize->Small];
slc=ListPlot[c,PlotStyle->{Gray,PointSize[.03],Opacity[.6]}];
slpod1=Show[ListPlot[podaci,PlotStyle->{PointSize[.008]}],slc,AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG,ImageSize->Small]
z0={{3,3},{4,9},{4,5},{8,4},{8,1}};k=Length[z0];
slz0=ListPlot[z0,PlotStyle->{Gray,PointSize[.04],Opacity[.9]}];
slpod2=Show[slA,AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG,ImageSize->Small]
GraphicsGrid[{{slpod1,slpod2}},ImageSize->Medium]
							(* k-means *)
time=Timing[sol=WKMeans[w,A,z0]];
Print["F= ",WF[sol[[1]]],   ";  G= ",WG[sol[[1]]],"; CPU=",time[[1]], "; It = ",sol[[4]]]
							(* SLIKE *)
Print["k-means reconstruction for selected centers z0: ",z0]
slc=ListPlot[sol[[2]],PlotMarkers->{"✶",20},PlotStyle->{Gray,PointSize[.025],Opacity[.7]}];
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]]]}];
sl3=Show[slA,slc,AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG];
Print[GraphicsGrid[{{slpod1,sl3}},ImageSize->Medium]];
tab=Table[sol[[1,j,All,2]],{j,Length[sol[[1]]]}];
ind=Rand[podaci,tab];
Print["Hausdorff dist. = ",Hdist[sol[[2]],c]]
Print["{Rand, Jaccard} index = ",{ind[[1]],ind[[2]]}]
Print["Confusion matrix S= ",ind[[3]]//MatrixForm]

Test-Example~9.3:  6 clusters in R^2

SeedRandom[1213];
c={{-8,-6},{-3,-8},{-1,-1},{8,4},{6,9},{-1,8}};
var=1.5; broj=200; AG={{-15,15},{-15,15}};
podaci=Table[RandomVariate[MultinormalDistribution[ci,var IdentityMatrix[2]],broj],{ci,c}];
pod=Flatten[podaci,1]; 
m=Length[pod];
w=Table[1,{i,m}];
Dimensions[pod]	
						(* PODACI *)
A=Flatten[podaci,1];
xG={Min[A[[All,1]]]-.5,Max[A[[All,1]]]+.5};yG={Min[A[[All,2]]]-.5,Max[A[[All,2]]]+.5};
AG={xG,yG};
m=Length[A];
w=Table[1,{i,m}]; 
slA=ListPlot[A,PlotStyle->{PointSize[.01]},AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG,ImageSize->Small]
slc=ListPlot[c,PlotStyle->{Gray,PointSize[.03],Opacity[.6]}];
slpod1=Show[ListPlot[podaci,PlotStyle->{PointSize[.01]}],slc,AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG,ImageSize->Small]
z0={{-10,-8},{0,-8},{0,0},{0,8},{8,2},{8,10}};k=Length[z0];
slz0=ListPlot[z0,PlotStyle->{Red,PointSize[.04],Opacity[.8]}];
slpod2=Show[slA,slz0,AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG];
GraphicsGrid[{{slpod1,slpod2}},ImageSize->Medium]
							(* k-means *)
time=Timing[sol=WKMeans[w,A,z0]];
Print["F= ",WF[sol[[1]]],   ";  G= ",WG[sol[[1]]],"; CPU=",time[[1]],, "; It = ",sol[[4]]]
							(* SLIKE *)
Print["k-means reconstruction for selected centers z0: ",z0]
slc=ListPlot[sol[[2]],PlotMarkers->{"✶",20},PlotStyle->{Red,PointSize[.025],Opacity[.7]}];
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]]]}];
sl3=Show[slA,slc,AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG];
Print[GraphicsGrid[{{slpod1,sl3}},ImageSize->Medium]];
tab=Table[sol[[1,j,All,2]],{j,Length[sol[[1]]]}];
ind=Rand[podaci,tab];
Print["Hausdorff dist. = ",Hdist[sol[[2]],c]]
Print["{Rand, Jaccard} index = ",{ind[[1]],ind[[2]]}]
Print["Confusion matrix S= ",ind[[3]]//MatrixForm]

Test-Example~9.4:  4 clusters in R^3

(* U ovom slučaju za Incewmwntal algorithm potrebno je u Compile prilagoditi funkciju F 3 varijable *)
SeedRandom[1213];
c={{2.5,3,2},{2,8.5,4},{5,6,4},{8,5,6}};
var=1.; broj=100;  n=Length[c[[1]]];
podaci=Table[RandomVariate[MultinormalDistribution[ci,var IdentityMatrix[3]],broj],{ci,c}];
pod=Flatten[podaci,1]; 
m=Length[pod];
w=Table[1,{i,m}];
Dimensions[pod]
						(* DATA *)
A=Flatten[podaci,1];
AG=Table[{Min[A[[All,i]]]-.5,Max[A[[All,i]]]+.5},{i,n}];
m=Length[A];
w=Table[1,{i,m}]; 
slc=ListPointPlot3D[c,PlotStyle->{Black,PointSize[.04],Opacity[.8]},PlotRange->AG,ViewPoint->{1.3,-2.4,2.},ImageSize->Small];
slpod1=Show[ListPointPlot3D[podaci,PlotStyle->{PointSize[.012]}],slc,AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG,ImageSize->Small]
slA=ListPointPlot3D[A,PlotRange->AG,PlotStyle->{PointSize[.012]},ViewPoint->{1.3,-2.4,2.},ImageSize->Small]
slA=ListPointPlot3D[A,PlotRange->AG,PlotStyle->{PointSize[.012]},ViewPoint->{-2.3,1.4,1},ImageSize->Small]
z0={{3,3,3},{3,5,3},{5,5,5},{9,3,5}};k=Length[z0];
(* k-means *)
time=Timing[sol=WKMeans[w,A,z0]];
Print["F= ",WF[sol[[1]]],   ";  G= ",WG[sol[[1]]],"; CPU=",time[[1]], "; It = ",sol[[4]]]
							(* SLIKE *)
Print["k-means reconstruction for selected centers z0: ",z0]
slc=ListPointPlot3D[sol[[2]], ViewPoint->{1.3,-2.4,2.},PlotStyle->{Gray,PointSize[.04],Opacity[.7]}];
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]]]}];
sl3=Show[slc,slA,AspectRatio->Automatic,GridLines->Automatic,PlotRange->AG];
Print[GraphicsGrid[{{slpod1,sl3}},ImageSize->Medium]];
tab=Table[sol[[1,j,All,2]],{j,Length[sol[[1]]]}];
ind=Rand[podaci,tab];
Print["Hausdorff dist. = ",Hdist[sol[[2]],c]]
Print["{Rand, Jaccard} index = ",{ind[[1]],ind[[2]]}]
Print["Confusion matrix S= ",ind[[3]]//MatrixForm]

Test-example-5:  small set A={...} in R^2

A={{1,5},{2,4},{3,6},{2,7},{2,1},{3,2},{4,1},{3,1},{6,5},{8,5},{7,4},{7,6}}//N; w=Table[1,{i,Length[A]}];   z0={{3,2},{2,5},{6,3}}; 
slpod2=ListPlot[{A,z0},PlotStyle->{PointSize[.05]},AspectRatio->Automatic,GridLines->Automatic];
time=Timing[sol=WKMeans[w,A,z0]];
Print["F= ",WF[sol[[1]]]//N,   ";  G= ",WG[sol[[1]]]//N,"; CPU=",time[[1]], "; It = ",sol[[4]]];
							(* Figures *)
Print["k-means reconstruction for selected centers z0: ",z0]
slc=ListPlot[sol[[2]],PlotMarkers->{"✶",20},PlotStyle->{Gray,PointSize[.025],Opacity[.7]}];
slPI=ListPlot[Table[sol[[1,j,All,2]],{j,Length[sol[[2]]]}],PlotStyle->Table[{PointSize[.05],Hue[j/k]},{j,Length[sol[[2]]]}]];
sl3=Show[slPI,slc,AspectRatio->Automatic,AxesOrigin->{0,0}];
Print[GraphicsGrid[{{slpod2,sl3}},ImageSize->Medium]]

Call Incremental algorithm

inc=Inc[w,A,{Mean[A]},7,1];
IndCH=Round[inc[[2]],-.1];IndDB=Round[inc[[3]],-.01];IndSSWC=Round[inc[[4]],-.01];IndDunn=Round[inc[[5]],-.001]; 
ListPlot[Table[{i+1,IndCH[[i]]},{i,Length[IndCH]}],Joined->True,AxesOrigin->{0,1000},GridLines->Automatic,ImageSize->Small]
ListPlot[Table[{i+1,IndDB[[i]]},{i,Length[IndCH]}],Joined->True,AxesOrigin->{0,.5},GridLines->Automatic,ImageSize->Small]
ListPlot[Table[{i+1,IndSSWC[[i]]},{i,Length[IndCH]}],Joined->True,AxesOrigin->{0,.75},GridLines->{Automatic,{.8,.85,.9}},ImageSize->Small]
ListPlot[Table[{i+1,IndDunn[[i]]},{i,Length[IndDunn]}],Joined->True,AxesOrigin->{0,0},GridLines->Automatic,ImageSize->Small]
Print["IndCH: ",IndCH,"\nIndDB: ",IndDB,"\nIndSSWC: ",IndSSWC,"\nIndDunn: ",IndDunn];

Optimal partition

(* for z0 take optimal centers from Incremental algorithm which corresponds the MAPart *)
z0={{7.832862177883973`,7.895305506607625`},{1.9116806835626434`,2.210985686097703`},{7.010453103277853`,2.881342410340607`},{3.0043499006777736`,4.960522416761811`},{5.832883341424164`,6.870614316633205`}};
k=Length[z0];
(* k-means *)
sol=WKMeans[w,A,z0];
Print["F= ",WF[sol[[1]]],   ";  G= ",WG[sol[[1]]]]		
							(* SLIKE *)
Print["k-means reconstruction for selected centers z0: ",z0]
slc=ListPlot[sol[[2]],PlotMarkers->{"✶",20},PlotStyle->{Gray,PointSize[.025],Opacity[.7]}];
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]]]}];
sl3=Show[slPI,slc,AspectRatio->Automatic,AxesOrigin->{0,0}];
Print[GraphicsGrid[{{slpod1,sl3}},ImageSize->500]];
tab=Table[sol[[1,j,All,2]],{j,Length[sol[[1]]]}];
ind=Rand[podaci,tab];
Print["Hausdorff dist. = ",Hdist[sol[[2]],c]]
Print["{Rand, Jaccard} index = ",{ind[[1]],ind[[2]]}]
Print["Confusion matrix S= ",ind[[3]]//MatrixForm]