Modules

Distance between the point T∈(ℝ^n)and the circle Circ(c, r)in (ℝ^n)

Input: point T∈(ℝ^n), center c ∈ (ℝ^n)and radius r >0;
Output: Distance from the point T to the circle Circ(c, r);
Submodules: None

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

Sphere data generating

(* A - dataset around the sphere; sph - dataset on the sphere *)
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
]

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];
(* Scaling the function onto the unit hypercube *)
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]]};
(* 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];
]

Data sphere (according to Marsaglia1972)

SeedRandom[7]
m=1000;c={5,5,5};r=4;sigma=.05;
AG=Table[{0,10},3];
ptsM=RandomVariate[UniformDistribution[{{-r,r},{-r,r}}],m];
Data={};
Table[If[ptsM[[i,1]]^2+ptsM[[i,2]]^2<r^2,Data=Append[Data,ptsM[[i]]]],
{i,Length[ptsM]}];
sph=Table[{c[[1]]+(2/r) Data[[i,1]]Sqrt[r^2-Data[[i,1]]^2-Data[[i,2]]^2], c[[2]]+(2/r) Data[[i,2]]Sqrt[r^2-Data[[i,1]]^2-Data[[i,2]]^2],
c[[3]]+(1/r)(r ^2-2(Data[[i,1]]^2+Data[[i,2]]^2))},{i,Length[Data]}];
A=Table[RandomReal[MultinormalDistribution[sph[[i]],sigma IdentityMatrix[3]],{1}][[1]],{i,Length[sph]}];
slc=ListPointPlot3D[{c}, PlotStyle->{Black,PointSize[.05],Opacity[.5]}];
Show[ListPointPlot3D[sph,PlotStyle->{PointSize->.015} ],
Graphics3D[{Opacity[.2],Sphere[c,r]}],slc,
PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},ImageSize->Small]
Show[ListPointPlot3D[A,PlotStyle->{PointSize->.015} ],
Graphics3D[{Opacity[.2],Sphere[c,r]}],slc,
PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},ImageSize->Small]

Data sphere (Mathematica)

SeedRandom[7]
m=1000;c={5,5,5};r=4;sigma=.05;
AG=Table[{0,10},3];
sph=RandomPoint[Sphere[c,r],m];
A=Table[RandomReal[MultinormalDistribution[sph[[i]],sigma IdentityMatrix[3]],{1}][[1]],{i,Length[sph]}];
slc=ListPointPlot3D[{c}, PlotStyle->{Orange,PointSize[.05],Opacity[.8]}];
Show[ListPointPlot3D[sph,PlotStyle->{PointSize->.015} ],
Graphics3D[{Opacity[.2],Sphere[c,r]}],slc,
PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},ImageSize->Small]
Show[ListPointPlot3D[A,PlotStyle->{PointSize->.015} ],
Graphics3D[{Opacity[.2],Sphere[c,r]}],slc,
PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},ImageSize->Small]

Data with outliers

SeedRandom[7]
m=1000;c={5,5,5};r=4;sigma=.05;
AG=Table[{0,10},3];
slSph=Graphics3D[{Opacity[.2],Sphere[c,r]}];
		(* Generating data around the sphere *)
data=DataSph[c,r,sigma,m];
A=data[[1]];
slc=ListPointPlot3D[{c}, PlotStyle->{Orange,PointSize[.05],Opacity[.8]}];
slA=ListPointPlot3D[A,PlotStyle->{PointSize->.015} ];
			(* Outliers with one stain *)
mout=200;
outliers=RandomVariate[UniformDistribution[{{7.5,10},{7.5,10},{7.5,10}}],mout];
slout=ListPointPlot3D[outliers,PlotStyle->{Red,PointSize->.015} ];
			(* Outliers with two stains  
mout=200;
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} ];*)
			(* data *)
A1=Union[A,outliers];
 m=Length[A]+mout;
Print["|A| = ",{Length[A],mout}," = ",m]
Show[slA,slout,slSph,slc,
PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},ImageSize->Small]

LAD approach

F[p_,q_,t_,rr_]:=Sum[Dist1[A1[[i]],{p,q,t},rr],{i,Length[A1]}]
C0=Median[A1]; r0=Median[Table[Norm[C0-A1[[i]]],{i,m}]];Print["C0=",Round[C0,-.05],"; ||C-C1|| = ",Round[Norm[c-C0],-.05]"; r0=",Round[r0,-.05],"; del r0 = ",Round[100Abs[r-r0]/r,-.05],"%   F0=",Round[F[C0[[1]],C0[[2]],C0[[3]],r0],-.05]]
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]];
Print["C1=",Round[c1,-.05],"; ||C-C1|| = ",Round[Norm[c-c1],-.05],"; r1=",Round[r1,-.05],"; del r = ",Round[100Abs[r-r1]/r,-.1],"%  F1=",Round[T[[2,1]],-.05],"; CPU=",T[[1]]];
list=Table[Norm[A[[i]]-c1]-r1,{i,Length[A]}]; 
Print["Stdev. = ",StandardDeviation[list]];
slSph1=Graphics3D[{Green,Opacity[.2],Sphere[c1,r1]}];
slc1=ListPointPlot3D[{c1}, PlotStyle->{Red,PointSize[.03],Opacity[.8]}];
Show[slA,slout,slSph1,slc,slc1,
PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},ImageSize->Small]

TLS approach

F[p_,q_,t_,rr_]:=Sum[Dist2[A1[[i]],{p,q,t},rr],{i,Length[A1]}]
C0=Mean[A1]; r0=Mean[Table[Norm[C0-A1[[i]]],{i,m}]];Print["C0=",Round[C0,-.05],"; ||C-C0|| = ",Round[Norm[c-C0],-.05]"; r0=",Round[r0,-.05],"; del r0 = ",Round[100Abs[r-r0]/r,-.05],"%   F0=",Round[F[C0[[1]],C0[[2]],C0[[3]],r0],-.05]]
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]];
Print["C1=",Round[c1,-.05],"; ||C-C1|| = ",Round[Norm[c-c1],-.05],"; r1=",Round[r1,-.05],"; del r = ",Round[100Abs[r-r1]/r,-.1],"%  F1=",Round[T[[2,1]],-.05],"; CPU=",T[[1]]];
list=Table[Norm[A[[i]]-c1]-r1,{i,Length[A]}]; 
Print["Stdev. = ",StandardDeviation[list]];
slSph1=Graphics3D[{Green,Opacity[.2],Sphere[c1,r1]}];
slc1=ListPointPlot3D[{c1}, PlotStyle->{Red,PointSize[.03],Opacity[.8]}];
Show[slA,slout,slSph1,slc,slc1,
PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},ImageSize->Small]

Algebraic approach

F[p_,q_,t_,rr_]:=Sum[DistA[A1[[i]],{p,q,t},rr],{i,Length[A1]}]
C0=Mean[A1]; r0=Sqrt[Mean[Table[Norm[C0-A1[[i]]]^2,{i,m}]]];Print["C0=",Round[C0,-.05],"; ||C-C1|| = ",Round[Norm[c-C0],-.05]"; r0=",Round[r0,-.05],"; del r0 = ",Round[100Abs[r-r0]/r,-.1],"%   F0=",Round[F[C0[[1]],C0[[2]],C0[[3]],r0],-.05]]
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]];
slSph1=Graphics3D[{Green,Opacity[.2],Sphere[c1,r1]}];
Print["C1=",Round[c1,-.05],"; ||C-C1|| = ",Round[Norm[c-c1],-.05],"; r1=",Round[r1,-.05],"; del r = ",Round[100Abs[r-r1]/r,-.1],"%  F1=",Round[T[[2,1]],-.05],"; CPU=",T[[1]]];
list=Table[Norm[A[[i]]-c1]-r1,{i,Length[A]}]; 
Print["Stdev. = ",StandardDeviation[list]];
slc1=ListPointPlot3D[{c1}, PlotStyle->{Red,PointSize[.03],Opacity[.8]}];
Show[slA,slout,slSph1,slc,slc1,
PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},ImageSize->Small]

The best representative - call of modul

SphRep[A1,1]
SphRep[A1,2]
SphRep[A1,3]