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]]]}]//N]
];

FF[PI_]:=Module[{k=Length[PI]},
cen=Table[0,{j,k}];
Do[cen[[j]]=Mean[PI[[j]]],{j,k}]; 
Sum[Sum[Norm[PI[[j,s]]-cen[[j]]]^2,{s,Length[PI[[j]]]}]//N,{j,k}]
]

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


Fmin[A_,c_]:=Module[{m=Length[A],k=Length[c]},
Sum[Min[Table[Norm[A[[i]]-c[[j]]]^2,{j,k}]],{i,m}]
]

Distance between the point T and the circle Circ (c, r)

(* LAD distance *)
Dist1[T_,c_,r_]:=Module[{},Return[Abs[Norm[T-c]-r] ]
                                             ];       
				(* TLS quasi-distance *)
Dist2[T_,c_,r_]:=Module[{},Return[(Norm[T-c]-r)^2 ]
                                             ]; 
			(* Algebraic fit criterion  *)
DistA[T_,c_,r_]:=Module[{},Return[(Norm[T-c]^2-r^2)^2]
                                             ];

Hausdorff distance between two sets

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

Hausdorff distance between two sets of k circles

Input: Circles CC={(Subscript[C, j])((Subscript[c, j]),(Subscript[r, j])), j=1,⃛,k}; C1={(Subscript[(C^1), j])((Subscript[(c^1), j]),(Subscript[(r^1), j])), j=1,⃛,k};  (Subscript[C, j]), (Subscript[(C^1), j]) ((ϵR)^3)
Output: Hausdorff distance;
Submodules: none needed

HdistCircles[CC_,C1_]:=Module[{m=Length[CC]},
Dist[A_,B_]:=Norm[A[[1]]-B[[1]]]+Abs[A[[2]]-B[[2]]]; (* Udaljenost kružnica A i B*)
Max[Max[Table[Min[Table[Dist[CC[[i]],C1[[j]]],{j,m}]],{i,m}]],Max[Table[Min[Table[Dist[CC[[i]],C1[[j]]],{i,m}]],{j,m}]]]]

Sphere data generating

(* Following Marsaglia1972: A - dataset around the sphere; sph - dataset on the sphere *)
DataSph[c_,r_,sigma_,m_]:=Module[{sph,A},
ptsM=RandomVariate[UniformDistribution[{{-r,r},{-r,r}}],m];
Data={};
Table[If[Norm[ptsM[[i]]]^2<r^2,Data=Append[Data,ptsM[[i]]]],
{i,Length[ptsM]}];
sph=Table[{c[[1]]+(2/r) Data[[i,1]]Sqrt[r^2-Norm[Data[[i]]]^2], c[[2]]+(2/r) Data[[i,2]]Sqrt[r^2-Norm[Data[[i]]]^2],
c[[3]]+(1/r)(r ^2-2Norm[Data[[i]]]^2)},{i,Length[Data]}];
A=Table[RandomReal[MultinormalDistribution[sph[[i]],sigma IdentityMatrix[3]],{1}][[1]],{i,Length[sph]}];
{A,sph,Length[sph]}
]

(* Data generating with 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}
]

Sphere representative

(* LAD(approach=1),  TLS(approach=2); Algebraic (approach=3)*)
SphRep[Data_,approach_]:=Module[{C0,r0,p1,q1,t1,c1,r1,F1,T,min,out},
Which[approach==1,
F[p_,q_,t_,rr_]:=Sum[Dist1[Data[[i]],{p,q,t},rr],{i,Length[Data]}];
C0=Median[Data]; r0=Median[Table[Norm[C0-Data[[i]]],{i,Length[Data]}]];
T=Timing[min=FindMinimum[F[p,q,t,rr],{p,C0[[1]]},{q,C0[[2]]},{t,C0[[3]]},{rr,r0},Method->"Newton"]];
p1=p/.min[[2]]; q1=q/.min[[2]]; t1=t/.min[[2]];
c1={p1,q1,t1}; r1=rr/.min[[2]];F1=T[[2,1]];
out={{c1,r1,F1,T[[1]]}},
approach==2,
F[p_,q_,t_,rr_]:=Sum[Dist2[Data[[i]],{p,q,t},rr],{i,Length[Data]}];
C0=Mean[Data]; r0=Mean[Table[Norm[C0-Data[[i]]],{i,Length[Data]}]];
T=Timing[min=FindMinimum[F[p,q,t,rr],{p,C0[[1]]},{q,C0[[2]]},{t,C0[[3]]},{rr,r0},Method->"Newton"]];
p1=p/.min[[2]]; q1=q/.min[[2]]; t1=t/.min[[2]];
c1={p1,q1,t1}; r1=rr/.min[[2]]; F1=T[[2,1]];
out={{c1,r1,F1,T[[1]]}},
approach==3,
F[p_,q_,t_,rr_]:=Sum[DistA[Data[[i]],{p,q,t},rr],{i,Length[Data]}];
C0=Mean[Data]; r0=Sqrt[Mean[Table[Norm[C0-Data[[i]]]^2,{i,Length[Data]}]]];
T=Timing[min=FindMinimum[F[p,q,t,rr],{p,C0[[1]]},{q,C0[[2]]},{t,C0[[3]]},{rr,r0}]];
p1=p/.min[[2]]; q1=q/.min[[2]]; t1=t/.min[[2]];
c1={p1,q1,t1}; r1=rr/.min[[2]];F1=T[[2,1]];
out={{c1,r1,F1,T[[1]]}}
];
out
]

KCS (k - closest sphere algorithm) - general

(* If Ind≠0, shows intermediate steps with pictures *)
KCS[A_,Circ0_,Iter_,Ind_]:=Module[{it=0,eps=0.05,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];
				 (* Calculate circle-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[
(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[.5],Sphere[Sph[[1,1]],Sph[[1,2]]]},{Opacity[.5],Sphere[Sph[[3,1]],Sph[[3,2]]]},{Opacity[.5],Sphere[Sph[[2,1]],Sph[[2,2]]]}}];
Print[Show[
ListPointPlot3D[PI,PlotStyle->{PointSize->.008}],sl5,slout,PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},ImageSize->Small]];];
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}
]

KCS (k - closest sphere algorithm) - including comparison with the original

KCS[A_,Circ0_,Iter_,Ind_]:=Module[{it=0,PI,Sph,tab,fc,min,imin,F0,F1,eps=0.0005},
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,Length[A]}];
PI=DeleteCases[PI,{}];
k=Length[PI]; 
F0=F00=Fcir[PI,Circ0]; 
(* This delhat F is only for text *)
Print["hat F=",F0,"; delhat F = ",100(F0-Forig)/Forig];
				 (* Calculate circle-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[
(*Chop[F0-F1]≠0&&it<Iter,*)
(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,"%\nHd(S0,Sph): ",HdistCircles[S0,Sph]];
sl5= Graphics3D[{{Opacity[.3],Sphere[Sph[[1,1]],Sph[[1,2]]]},{Opacity[.3],Sphere[Sph[[2,1]],Sph[[2,2]]]},{Opacity[.3],Sphere[Sph[[3,1]],Sph[[3,2]]]}}];
Print[Show[
ListPointPlot3D[{PI[[1]],PI[[3]],PI[[2]]},PlotStyle->{PointSize->.008}],sl5,slout,PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},ImageSize->Small]];];
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,Length[A]}];
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,(F0-F1)/F00}
]

DIRECT

Modul for solving the global optimization problem for a Lipschitz-continuous function. 
Mathematica code written by I.Vazler;
Input: ordinary or compiled function, Domain , {Eps->10^(-4),MaxIter->20,DMin->10.^(-2)},
Output: optimal variables values, optimal objective function value , CPU-time, ;
Submodules: none needed

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

Three well separated spheres

n=3; k=3;
C0={{6,6,6},{2,3,3},{2,5,8.5}}; R0={3,2,1.5};
A=Table[0,{j,k}]; AG=Table[{0,10},3];
sl3= Graphics3D[{{Opacity[.3],Sphere[C0[[1]],R0[[1]]]},{Opacity[.3],Sphere[C0[[2]],R0[[2]]]},{Opacity[.3],Sphere[C0[[3]],R0[[3]]]}},Axes->True,AxesLabel->{x,y,z},PlotRange->AG,ImageSize->Small];
slc=ListPointPlot3D[{C0}, PlotStyle->{Gray,PointSize[.04],Opacity[.5]}];
Show[sl3,slc]

m=4500.;
p={5,10,20}; mout=m p/100
m+mout

Three spheres : Dataset

SeedRandom[17]
		(* Sphere 1 *)
sigma=0.05;m1=2000;
data=DataSph[C0[[1]],R0[[1]],sigma,m1];
A1=data[[1]];
		(* Sphere 2 *)
sigma=0.05;m2=1500;
data=DataSph[C0[[2]],R0[[2]],sigma,m2];
A2=data[[1]];
		(* Sphere 3 *)
sigma=0.05;m3=1000;
data=DataSph[C0[[3]],R0[[3]],sigma,m3];
A3=data[[1]];
				(* outliers: uniformly distributed *)
mout=250; mout=450; mout=900;   
mout=250;
outliers=RandomVariate[UniformDistribution[{{0,10},{0,10},{0,10}}],mout];
slout=ListPointPlot3D[outliers,PlotStyle->{Red,PointSize->.012} ];
				(* outliers: non-uniformly distributed 
mout=170;
outliers1=RandomVariate[UniformDistribution[{{7.5,10},{7.5,10},{7.5,10}}],mout];
outliers2=RandomVariate[UniformDistribution[{{7.5,10},{0,3},{0,3}}],mout];
outliers=Join[outliers1,outliers2];
slout=ListPointPlot3D[outliers ,PlotStyle->{Red,PointSize->.015} ];*)
						(* Dataset *)
A=Join[A1,A2,A3,outliers];
Print["|Sphere1|= ",Length[A1],"; |Sphere2|= ",Length[A2],"; |Sphere3|= ",Length[A3],"; Out=",mout, ";   |A|= ",Length[A]," with ", Round[100. mout/Length[A],-.01],"% outliers"]
slData=Show[sl3,
ListPointPlot3D[{A1,A2,A3},PlotStyle->{PointSize->.008}],slout,PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},ImageSize->Small]
Print[Style["Data set A:",16,Blue]];
slA=Show[ListPointPlot3D[{A1,A3,A2},PlotStyle->{PointSize->.008}],slout,PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},ImageSize->Small]

Determining of Forig

(* Concatenation of original spheres  *)
S0=Table[{C0[[j]],R0[[j]]},{j,k}];  
(* Determining corresponding clusters and the function value F0 *)
PI=Table[{},{i,k}];
Do[
tab=Table[Dist1[A[[i]],S0[[j,1]],S0[[j,2]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],A[[i]]],{i,Length[A]}];
Forig=Fcir[PI,S0]; 
Print["Forig = ",Forig];

Initial approximation and Solution 
(including relative deviation from the original)

n=3; k=3;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->10, 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["|A|=",Length[A],"; cenDIR: ",Round[cenDIR,-.001]]
				(* Initial spheres of radii 1 and Hausdorff distance Hdist[Cp,cenDIR]*)
CircDIR=Table[{cenDIR[[j]],1},{j,k}];
Print["Hdist[C0,cenDIR]=",Round[Hdist[C0,cenDIR],-.001],"; F0 = ",min[[2]],"; T(DIR)=",T[[1]]]
sl4= Graphics3D[{{Opacity[.4],Sphere[CircDIR[[1,1]],CircDIR[[1,2]]]},{Opacity[.4],Sphere[CircDIR[[2,1]],CircDIR[[2,2]]]},{Opacity[.4],Sphere[CircDIR[[3,1]],CircDIR[[3,2]]]}}];
Print[Show[
ListPointPlot3D[PI,PlotStyle->{PointSize->.008}],sl4,slout,PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},ImageSize->Small]];

					(* Call of KCS and Solution *)
Iter=20;
TTotal=Timing[kcs=KCS[A,CircDIR,Iter,1]];
If[kcs[[4]]<Iter,
Print["Spheres: ",C1=Round[kcs[[2]],-.001],"\nF1 = ",Round[kcs[[3]],-.001],"; del F=",kcs[[5]],"; HdistCircles[S0,C1]=",Round[HdistCircles[S0,kcs[[2]]],-.001],"\n It= ",kcs[[4]],"; Time(DIR)=",T[[1]],"; Time(KCS)=",TTotal[[1]]],
Print["In ",Iter," iterations no solution found"]
]

Three intersecting spheres

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[.3],Sphere[C0[[1]],R0[[1]]]},{Opacity[.3],Sphere[C0[[2]],R0[[2]]]},{Opacity[.3],Sphere[C0[[3]],R0[[3]]]}},Axes->True,AxesLabel->{x,y,z},PlotRange->AG,ImageSize->Small];
slc=ListPointPlot3D[{C0}, PlotStyle->{Gray,PointSize[.04],Opacity[.5]}];
Show[sl3,slc]

Three spheres : Dataset

(* SeedRandom[13]	*)
SeedRandom[13]
	(* Sphere 1 *)
sigma=0.05;m1=2000;
data=DataSph[C0[[1]],R0[[1]],sigma,m1];
A1=data[[1]];
		(* Sphere 2 *)
sigma=0.05;m2=1500;
data=DataSph[C0[[2]],R0[[2]],sigma,m2];
A2=data[[1]];
		(* Sphere 3 *)
sigma=0.05;m3=1000;
data=DataSph[C0[[3]],R0[[3]],sigma,m3];
A3=data[[1]];
						(* outliers: uniformly distributed *)
 mout=250;    mout=450; 
mout=250;
outliers=RandomVariate[UniformDistribution[{{0,10},{0,10},{0,10}}],mout];
slout=ListPointPlot3D[outliers,PlotStyle->{Red,PointSize->.012} ];
				(* outliers: non-uniformly distributed 
mout=30;
outliers1=RandomVariate[UniformDistribution[{{7.5,10},{7.5,10},{7.5,10}}],mout];
outliers2=RandomVariate[UniformDistribution[{{7.5,10},{0,3},{0,3}}],mout];
outliers=Join[outliers1,outliers2];
slout=ListPointPlot3D[outliers ,PlotStyle->{Red,PointSize->.015} ];*)
						(* Dataset *)
A=Union[A1,A2,A3,outliers];
Print["|Sphere1|= ",Length[A1],"; |Sphere2|= ",Length[A2],"; |Sphere3|= ",Length[A3],"; mout = ",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},ImageSize->Small]
Print[Style["Data set A:",16,Blue]];
slA=Show[ListPointPlot3D[{A1,A2,A3},PlotStyle->{PointSize->.005}],slout,PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},ImageSize->Small]

Determining of Forig

(* Concatenation of original spheres  *)
S0=Table[{C0[[j]],R0[[j]]},{j,k}];  
(* Determining corresponding clusters and the function value F0 *)
PI=Table[{},{i,k}];
Do[
tab=Table[Dist1[A[[i]],S0[[j,1]],S0[[j,2]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],A[[i]]],{i,Length[A]}];
Forig=Fcir[PI,S0]; 
Print["Forig = ",Forig];

Initial approximation and Solution

n=3; k=3;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 spheres of radii 1 and Hausdorff distance Hdist[Cp,cenDIR]*)
CircDIR=Table[{cenDIR[[j]],1},{j,k}];
Print["Hdist[C0,cenDIR]=",Round[Hdist[C0,cenDIR],-.001],"; F0 = ",min[[2]],"; T(DIR)=",T[[1]]]
sl4= Graphics3D[{{Opacity[.4],Sphere[CircDIR[[1,1]],CircDIR[[1,2]]]},{Opacity[.4],Sphere[CircDIR[[2,1]],CircDIR[[2,2]]]},{Opacity[.4],Sphere[CircDIR[[3,1]],CircDIR[[3,2]]]}}];
Print[Show[
slA,sl4,slout,PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},ImageSize->Small]];

		(* Call of KCS and Solution *)			
 Iter=20;
TTotal=Timing[kcs=KCS[A,CircDIR,Iter,1]];
If[kcs[[4]]<Iter,
Print["Spheres: ",C1=Round[kcs[[2]],-.001],"\nF1 = ",Round[kcs[[3]],-.001],"; del F=",kcs[[5]],"; HdistCircles[S0,C1]=",Round[HdistCircles[S0,kcs[[2]]],-.001],"\n It= ",kcs[[4]],"; Time(DIR)=",T[[1]],"; Time(KCS)=",TTotal[[1]]],
Print["In ",Iter," iterations no solution found"]
]