KCS (k-means algorithm for multiple spheres)

Input: Data set 𝒜 ⊂ (ℝ^3), concatenation of sphere-centers, Nr.Iter., Ind (If Ind≠0 then the intermediate results and figures are shown );
Output: Optimal partition of the set  𝒜 ,spheres,the  objective function value, nr.Iter.;
Submodules: Dist1[], Fcir[]

(* If Ind≠0 then the intermediate results and the related figures are shown  *)
KCS[A_,Circ0_,Iter_,Ind_]:=Module[{it=0,eps=0.000005,m=Length[A],k=Length[Circ0],Sph,PI,tab,p0,q0,t0,r0,fc,min,imin,F00,F0,F1},
Sph=Table[0,{j,k}];
PI=Table[{},{i,k}];
Do[
tab=Table[Dist1[A[[i]],Circ0[[j,1]],Circ0[[j,2]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],A[[i]]],{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI]; 
F0=F00=Fcir[PI,Circ0]; 
Print["hat F=",F0];
				 (* Calculation of the sphere-centers *)
Do[
 {p0,q0,t0}=Median[PI[[j]]];
r0=Median[Table[Norm[ {p0,q0,t0}-PI[[j,s]]],{s,Length[PI[[j]]]}]];
fc[p_,q_,t_,rr_]:=Sum[Dist1[PI[[j,s]],{p,q,t},rr],{s,Length[PI[[j]]]}];
min=FindMinimum[fc[p,q,t,rr],{p,p0},{q,q0},{t,t0},{rr,r0}];
Sph[[j]]={{p/.min[[2]],q/.min[[2]],t/.min[[2]]},rr/.min[[2]]};
,{j,k}];
F1=Fcir[PI,Sph] ;  
						(* While loop *)
While[
(F0-F1)/F00>eps&&it<Iter,
it=it+1;
If[Ind≠0,
Print["It ",it,": Sphere: ", Sph,";  F = ",{F0,F1},"; del F = ",100(F0-F1)/F00];
sl5= Graphics3D[{{Opacity[.2],Sphere[Sph[[1,1]],Sph[[1,2]]]},{Opacity[.2],Sphere[Sph[[2,1]],Sph[[2,2]]]},{Opacity[.2],Sphere[Sph[[3,1]],Sph[[3,2]]]}}];
Print[Show[
ListPointPlot3D[PI,PlotStyle->{PointSize->.008}],sl5,slout,PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},ImageSize->300]];];
F0=F1;
PI=Table[{},{j,k}];
Do[
tab=Table[Dist1[A[[i]],Sph[[j,1]],Sph[[j,2]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],A[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI];
Sph=Table[0,{j,k}];
Do[
{p0,q0,t0}=Median[PI[[j]]];
r0=Median[Table[Norm[ {p0,q0,t0}-PI[[j,s]]],{s,Length[PI[[j]]]}]];
fc[p_,q_,t_,rr_]:=Sum[Dist1[PI[[j,s]],{p,q,t},rr],{s,Length[PI[[j]]]}];
min=FindMinimum[fc[p,q,t,rr],{p,p0},{q,q0},{t,t0},{rr,r0}];
Sph[[j]]={{p/.min[[2]],q/.min[[2]],t/.min[[2]]},rr/.min[[2]]};
,{j,k}];
F1=Fcir[PI,Sph]; 
];
{PI,Sph,F1,it}
]

Example

Modules

Needs["ComputationalGeometry`"]

Fcir[PI_,Circ_]:=Module[{},
Return[
Sum[Sum[Dist1[PI[[j,s]],Circ[[j,1]],Circ[[j,2]]],{s,Length[PI[[j]]]}]
,{j,Length[PI]}]
]]

	(* LAD distance *)
Dist1[T_,c_,r_]:=Module[{},Return[Abs[Norm[T-c]-r] ]
                                             ];       

(* Data generating by Mathematica *)
DataSph[c_,r_,sigma_,m_]:=Module[{sph,A},
sph=RandomPoint[Sphere[c,r],m];
A=Table[RandomReal[MultinormalDistribution[sph[[i]],sigma IdentityMatrix[3]],{1}][[1]],{i,Length[sph]}];
{A,sph}
]

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

(* KCS *)
(* If Ind≠0 then the intermediate results and the related figures are shown *)
KCS[A_,Circ0_,Iter_,Ind_]:=Module[{it=0,eps=0.000005,m=Length[A],k=Length[Circ0],Sph,PI,tab,p0,q0,t0,r0,fc,min,imin,F00,F0,F1},
Sph=Table[0,{j,k}];
PI=Table[{},{i,k}];
Do[
tab=Table[Dist1[A[[i]],Circ0[[j,1]],Circ0[[j,2]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],A[[i]]],{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI]; 
F0=F00=Fcir[PI,Circ0]; 
Print["F0=",F0];
				 (* Calculation of the sphere-centers *)
Do[
 {p0,q0,t0}=Median[PI[[j]]];
r0=Median[Table[Norm[ {p0,q0,t0}-PI[[j,s]]],{s,Length[PI[[j]]]}]];
fc[p_,q_,t_,rr_]:=Sum[Dist1[PI[[j,s]],{p,q,t},rr],{s,Length[PI[[j]]]}];
min=FindMinimum[fc[p,q,t,rr],{p,p0},{q,q0},{t,t0},{rr,r0}];
Sph[[j]]={{p/.min[[2]],q/.min[[2]],t/.min[[2]]},rr/.min[[2]]};
,{j,k}];
F1=Fcir[PI,Sph] ;  
						(* While loop *)
While[
(F0-F1)/F00>eps&&it<Iter,
it=it+1;
If[Ind≠0,
Print["It ",it,": Sphere: ", Sph,";  F = ",{F0,F1},"; del F = ",100(F0-F1)/F00];
sl5= Graphics3D[{{Opacity[.2],Sphere[Sph[[1,1]],Sph[[1,2]]]},{Opacity[.2],Sphere[Sph[[2,1]],Sph[[2,2]]]},{Opacity[.2],Sphere[Sph[[3,1]],Sph[[3,2]]]}}];
Print[Show[
ListPointPlot3D[PI,PlotStyle->{PointSize->.008}],sl5,slout,PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},ImageSize->300]];];
F0=F1;
PI=Table[{},{j,k}];
Do[
tab=Table[Dist1[A[[i]],Sph[[j,1]],Sph[[j,2]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],A[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI];
Sph=Table[0,{j,k}];
Do[
{p0,q0,t0}=Median[PI[[j]]];
r0=Median[Table[Norm[ {p0,q0,t0}-PI[[j,s]]],{s,Length[PI[[j]]]}]];
fc[p_,q_,t_,rr_]:=Sum[Dist1[PI[[j,s]],{p,q,t},rr],{s,Length[PI[[j]]]}];
min=FindMinimum[fc[p,q,t,rr],{p,p0},{q,q0},{t,t0},{rr,r0}];
Sph[[j]]={{p/.min[[2]],q/.min[[2]],t/.min[[2]]},rr/.min[[2]]};
,{j,k}];
F1=Fcir[PI,Sph]; 
];
{PI,Sph,F1,it}
]

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="The maximal number `1` of iterations has been reached.";
n=Length[dom];
(* Scaling the function onto the unit hypercube *)
ff=Function[N[f@@({##}(dom[[All,2]]-dom[[All,1]])+dom[[All,1]]),50]];
(* Initial centers and diameters *)
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]]};
(* MAIN LOOP *)
Monitor[
While[(dmaxrazni[[1]]>dtest)&&(iter<OptionValue[MaxIter]),
iter++;
(* STEP 1: Searching for the potentially optimal hyperrectangle *)
dmaxrazni=Union[dmax];
(* all data grouped according to d *)
pot0=Table[Select[Range[Length[fc]],dmax[[#]]==dmaxrazni[[i]]&],{i,Length[dmaxrazni]}];
(* data with the lawest f-value grouped according to d *)
pot1=Table[pot0[[i,Position[fc[[pot0[[i]]]],Min[fc[[pot0[[i]]]]]][[All,1]]]],{i,Length[dmaxrazni]}];
(* checking the rest of the requirements *)
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];
(* STEP 2: dividing the potentially optimal directions *)
Do[
(* dividing directions *)
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]];];
(* figure for the 2D case *)
(*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]}];
(* back to the original domain *)
TMin[[1]]=TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]];
(* printing the figure in the 2D case *)
(*If[n==2,Print[ListAnimate[anim,AnimationRunning->False]]];*)
Return[TMin];
]

Data

n=3; k=3;
C0={{6,6,6},{3,4,5},{3,5,7.5}};   R0={3,2,1.5};
A=Table[0,{j,k}]; AG=Table[{0,10},3];
sl3= Graphics3D[{{Opacity[.2],Sphere[C0[[1]],R0[[1]]]},{Opacity[.2],Sphere[C0[[2]],R0[[2]]]},{Opacity[.2],Sphere[C0[[3]],R0[[3]]]}},Axes->True,AxesLabel->{y,x,z},PlotRange->AG,ImageSize->Small];
(* Data with outliers *)
SeedRandom[13]
	(* Sphere 1 *)
sigma=0.05;m1=1000;
data=DataSph[C0[[1]],R0[[1]],sigma,m1];
A1=data[[1]];
		(* Sphere 2 *)
sigma=0.05;m2=750;
data=DataSph[C0[[2]],R0[[2]],sigma,m2];
A2=data[[1]];
		(* Sphere 3 *)
sigma=0.05;m3=500;
data=DataSph[C0[[3]],R0[[3]],sigma,m3];
A3=data[[1]];
(* outliers *)
mout=90;
outliers=RandomVariate[UniformDistribution[{{0,10},{0,10},{0,10}}],mout];
slout=ListPointPlot3D[outliers,PlotStyle->{Red,PointSize->.012} ];
A=Union[A1,A2,A3,outliers];
Print["|Sphere1|= ",Length[A1],"; |Sphere2|= ",Length[A2],"; |Sphere3|= ",Length[A],"; |Outliers| = ",mout, ";   |A|= ",Length[A]]
slData=Show[
ListPointPlot3D[{A1,A2,A3},PlotStyle->{PointSize->.008}],sl3,slout,PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},PlotLabel->"Data"];
slA=Show[ListPointPlot3D[A,PlotStyle->{PointSize->.008}],slout,PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},PlotLabel->"Data set A"];
GraphicsGrid[{{slData,slA}},ImageSize->500]

Execution

(* Initial approximation construction *)
dom=Table[{0,10},{j,k n}];
FCompile=Compile[{{z1,_Real},{z2,_Real},{z3,_Real},{z4,_Real},{z5,_Real},{z6,_Real},{z7,_Real},{z8,_Real},{z9,_Real}},
Evaluate[
Sum[Min[Norm[{z1,z2,z3}-A[[i]]],Norm[{z4,z5,z6}-A[[i]]],Norm[{z7,z8,z9}-A[[i]]]],{i,Length[A]}]]];
f[z1_,z2_,z3_,z4_,z5_,z6_,z7_,z8_,z9_]:=Apply[FCompile,{z1,z2,z3,z4,z5,z6,z7,z8,z9}];
T=Timing[min=DIRECT[f,dom,DMin->10.^(-3),MaxIter->20, Eps->10^(-4)]];
zz=min[[1]];
cenDIR=Table[0,{j,k}];
Do[
jk=n j-n+1;
cenDIR[[j]]=Table[zz[[jk+(s-1)]],{s,n}];
,{j,k}];
Print["cenDIR: ",Round[cenDIR,-.001]];
	(* Initial sphere-radus 1 and Hausdorff distance Hdist[Cp,cenDIR]*)
CircDIR=Table[{cenDIR[[j]],1},{j,k}];
Print["Hdist[C0,cenDIR]=",Round[Hdist[C0,cenDIR],-.001]];				
                   (* call the KCS and the Solution *)
		(* If Ind≠0 then the intermediate results and the related figures are shown *)
kcs=Timing[KCS[A,CircDIR,20,0]];
Print["Spheres: ",C1=Round[kcs[[2,2]],-.001],"\nF1 = ",Round[kcs[[2,3]],-.001],"; Hdist[C0,C1]=",Round[Hdist[C0,C1[[All,1]]],-.001],"; It= ",kcs[[2,4]],"; Time(DIR)=",T[[1]],"; Time(KCS)=",kcs[[1]]]