Modules

SymDIRECT

Modul for solving the global optimization problem for a symmetric Lipschitz-continuous function. 
Mathematica codes 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

(* For a given center and length of the side of the hyperrectangle the function test returns:
  1 If the whole hyperrectangle is located in the region \Delta
0 If a hyperrectangle lies in the region \Delta only partially
-1 If a hyperrectangle appears outside the region \Delta (or if a hyperrectangle has only one point in \Delta)*)
test[c_,dd_]:=Module[{n=Length[c],S,d=dd/2,unutra=1},
S=Table[1,{n}];
Do[
If[c[[i]]-d[[i]]<c[[i+1]]+d[[i+1]],
unutra=0;Break[];
]
,{i,1,n-1}];
If[unutra==0,
Do[
If[c[[i]]+S[[i]]*d[[i]]<c[[i+1]]+d[[i+1]],
If[c[[i]]+S[[i]]*d[[i]]>c[[i+1]]-d[[i+1]],
S[[i+1]]=-1,
unutra=-1;Break[];
]
]
,{i,1,n-1}];
];
Return[unutra];
]

Options[SymDIRECT]={Eps->0,MaxIter->20,DMin->10.^(-2)};
SymDIRECT[f_,dom_,OptionsPattern[]]:=Module[{n,ff,c,d,dmax,fc,dmaxrazni,pot0,pot1,pot2,temp,temp2,ind,novic,novid,novifc,noviunutra,red,TMin,dtest,anim={},unutra,brisi},
dtest=N[OptionValue[DMin]/Max[(dom[[All,2]]-dom[[All,1]])]];
SymDIRECT::maxiter="The maximal number `1` of iterations has been reached.";
iter=0; Nf=0;
n=Length[dom];
(* scaling the function on the hypercube with length of the side 1 *)
ff=Function[N[f@@({##}(dom[[All,2]]-dom[[All,1]])+dom[[All,1]]),50]];
(* initial center and diameters *)
c={Table[1/2,{i,n}]};
d={Table[1,{i,n}]};
unutra={test[c[[1]],d[[1]]]};
dmax={Max[d[[1]]]};
dmaxrazni=Union[dmax];
fc={ff@@c[[1]]};   Nf=Nf+1;
TMin={c[[1]],fc[[1]]};
(* MAIN LOOP *)
Monitor[
While[(dmaxrazni[[1]]>dtest)&&(iter<OptionValue[MaxIter]),
iter++;
(* STEP 1: searching for potentially optimal hyperrectangles *)
dmaxrazni=Union[dmax];
(* All data grouped according to d *)
pot0=Table[Select[Range[Length[fc]],dmax[[#]]==dmaxrazni[[i]]&],{i,Length[dmaxrazni]}];
(* Data with smallest f grouped according to d *)
pot1=Table[pot0[[i,Position[fc[[pot0[[i]]]],Min[fc[[pot0[[i]]]]]][[All,1]]]],{i,Length[dmaxrazni]}];
(* Checking other conditions *)
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 potentially optimal hyperrectangles *)
brisi={};
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]}];
Nf=Nf+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]}];
(* Eliminating those that are not inside *)
If[unutra[[pot]]==0,
(* Central hyperrectangle can go out or stay completely inside *)
If[test[c[[pot]],d[[pot]]]==-1,AppendTo[brisi,{pot}],unutra[[pot]]=test[c[[pot]],d[[pot]]]];
temp={};noviunutra={};
Do[
temp2=test[novic[[i]],novid[[i]]];
If[temp2==-1,AppendTo[temp,{i}],AppendTo[noviunutra,temp2]];
,{i,Length[novic]}];
novic=Delete[novic,temp];
novid=Delete[novid,temp];
novifc=Delete[novifc,temp];
,
(* If the divided hyperrectangle lies completely inside, all subhyperrectangles remain inside too *)
noviunutra=Table[1,{Length[novic]}];
];
c=Join[c,novic];
d=Join[d,novid];
dmax=Join[dmax,Map[Max,novid]];
fc=Join[fc,novifc];
unutra=Join[unutra,noviunutra];
(* Update of actual minimum *)
temp=Position[novifc,Min[novifc]][[1,1]];
If[novifc[[temp]]<TMin[[2]],TMin=N[{novic[[temp]],novifc[[temp]]},10]];
,{pot,pot2}];
(* Deletion of those hyperrectangles which went outside *)
c=Delete[c,brisi];
d=Delete[d,brisi];
dmax=Delete[dmax,brisi];
fc=Delete[fc,brisi];
unutra=Delete[unutra,brisi];
(*Print[iter,": ",{TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]],TMin[[2]]}];*)
If[iter==OptionValue[MaxIter],Message[SymDIRECT::maxiter,OptionValue[MaxIter]];];
(* Figure in 2D case *)
(*If[n==2 && iter<100,
AppendTo[anim,
Graphics[{{Green,Rectangle[]},Table[{Blue,Tooltip[Rectangle[c[[i]]-d[[i]]/2+d[[i]]/10,c[[i]]+d[[i]]/2-d[[i]]/10],unutra[[i]]]},{i,Length[d]}],Point[c]}]
];
];*)
];
,{iter,Length[pot2]}];
(* Back to the initial domain *)
TMin[[1]]=TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]];
(* Images printing in 2D case *)
(*If[n==2,Print[ListAnimate[anim,AnimationRunning->False]]];*)
Return[TMin];
]

Weighted objective functions

Input: partition Π={(Subscript[π, 1]),⃛,(Subscript[π, k])} of the set 𝒜 ⊂ (ℝ^n), (in case of WFmin[] only the set 𝒜), centers (Subscript[c, 1]),⃛,(Subscript[c, k]) ∈ (ℝ^n);
Output: weighted objective function values: ℱ(Π), \[FormalCapitalF](c), 𝒢(Π);
Submodules: none needed

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

WFmin[A_,w_,z_]:=Module[{},
Sum[w[[i]]Min[Table[Norm[z[[j]]-A[[i]]]^2,{j,Length[z]}]],{i,Length[A]}]
]

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

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

Weighted k -  Means

WKMeans[A_,w_,z_]:=Module[{m=Length[A],k=Length[z],n=Length[A[[1]]],c0,Pod,imin,PI,It=0,tab},
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]; 
						(* 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"]];
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];  
               ];
{PI,c0,WF[PI],It}]

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 Representative

(* LAD(approach=1),  TLS(approach=2); Algebraic (approach=3); eps = accuracy *)
SphRep[A_,C0_,r0_,eps_,approach_]:=Module[{p1,q1,t1,c1,r1,F1,T,min,out},
Which[approach==1,
F[rr_]:=Sum[Dist1[A[[i]],C0,rr],{i,Length[A]}];
T=Timing[min=FindMinimum[F[rr],{rr,r0},Method->"Newton"]];
r1=rr/.min[[2]]; F1=T[[2,1]];
out={approach,r1,F1,T[[1]]},
approach==2,
F[rr_]:=Sum[Dist2[A[[i]],C0,rr],{i,Length[A]}];
T=Timing[min=FindMinimum[F[rr],{rr,r0},Method->"Newton"]];
r1=rr/.min[[2]]; F1=T[[2,1]];
out={approach,r1,F1,T[[1]]},
approach==3,
F[rr_]:=Sum[DistA[A[[i]],C0,rr],{i,Length[A]}];
T=Timing[min=FindMinimum[F[rr],{rr,r0},Method->"Newton"]];
r1=rr/.min[[2]];F1=T[[2,1]];
out={approach,r1,F1,T[[1]]}
];
out
]

KCS (k - closest spheres ) - 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,tab,p0,q0,t0,r0,fc,min,imin,F00,F0,F1},
Sph=Table[0,{j,k}];
PI=Table[{},{i,k}];
Do[
tab=Table[Dist2[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, "; k= ",k, Table[Length[PI[[j]]],{j,k}]];
				 (* Sphere centers - LAD distance *)
Do[
r0=Median[Table[Norm[ C0-PI[[j,s]]],{s,Length[PI[[j]]]}]];
fc[rr_]:=Sum[Dist2[PI[[j,s]],C0,rr],{s,Length[PI[[j]]]}];
min=FindMinimum[fc[rr],{rr,r0}];
Sph[[j]]={C0,rr/.min[[2]]};
,{j,k}];
F1=Fcir[PI,Sph] ;  
Print["Sph: ",Sph,"; F1 = ",F1];
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[.3],Sphere[Sph[[2,1]],Sph[[2,2]]]},{Opacity[.8],Sphere[Sph[[3,1]],Sph[[3,2]]]}}];
Print[Show[
ListPointPlot3D[PI,PlotStyle->{{Blue,PointSize->.013},{Red,PointSize->.013},{Green,PointSize->.013}}],sl5,slSph1,PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},ImageSize->Small]];];
F0=F1;
PI=Table[{},{j,k}];
Do[
tab=Table[Dist2[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[
r0=Median[Table[Norm[ C0-PI[[j,s]]],{s,Length[PI[[j]]]}]];
fc[rr_]:=Sum[Dist2[PI[[j,s]],C0,rr],{s,Length[PI[[j]]]}];
min=FindMinimum[fc[rr],{rr,r0}];
Sph[[j]]={C0,rr/.min[[2]]};
,{j,k}];
F1=Fcir[PI,Sph]; 
];
{PI,Sph,F1,it}
]

Data + Initial approximation

SetDirectory[NotebookDirectory[]]; 
asd=Import["7-.txt","Table"];
asd[[1]]
m=Length[asd];Print["m = ",m]

Transformation to [0, 1]

C0={0,0,0}; Print["C0 = ",C0];
phi=asd[[All,3]]; phi=Pi phi/180.;
th=asd[[All,4]]; th=Pi th/180.;
depth=-asd[[All,5]]; 
 min=Min[depth]; (* the deepest *)
 max=Max[depth]; (* the shallowest *)
Print["{the deepest, the shallowest}: ",{min, max}," -> [0,1]"]
lin[t_]:=-t/min+1.
B=lin[depth];
II={0,1};
w=Table[1,{i,Length[depth]}];
n=1;k=3; 
line=Graphics[{Thickness[0.004], Line[{{0,0},{1,0}}]}];
slB=ListPlot[Table[{B[[j]],0},{j,m}],PlotStyle->{PointSize[.01]}];
slB=Show[line,slB,PlotRange->{II,{-.05,.05}},Axes->{True,False},Ticks->{None,None},AspectRatio->Automatic,ImageSize->Medium]
(* The set A *)
A=Table[B[[i]]{Cos[phi[[i]]] Cos[th[[i]]],Cos[phi[[i]]] Sin[th[[i]]],Sin[phi[[i]]]},{i,m}];
slA=ListPointPlot3D[A,PlotStyle->{Red,PointSize->.015} ];
Print[Style["{Norm[A[[j]]] = B[[j]]} for all j !!!",Blue,16]];

Initial approximation for k=1

k=1;
cFunk0[var_]:=Compile[var,Evaluate[Sum[Min[Table[Norm[Table[(Subscript[x, j]),{j,1+(jk-1)n,1+jk n-1}]-B[[i]]]^2,{jk, 1,k}]],{i,m}]]];
FFCompile=cFunk0[Table[{(Subscript[x, i]),_Real},{i,n k}]];

(* SymDIRECT *)
Print[Style["SymDIRECT:",16,Blue]];
dom=Table[II,{k}];
tSym=Timing[
minSym=SymDIRECT[FFCompile,dom,DMin->10^(-4),MaxIter->20,Eps->10^(-4)]
];
cenSymDir=minSym[[1]];
Print[Style["SymDIRECT Centers: ",Green,16],cenSymDir];
Print["F= ",minSym[[2]],"; IT = ",iter,"; CPU = ",tSym[[1]],"; Nf = ",Nf];
slcSym=ListPlot[Table[{cenSymDir[[j]],0},{j, k}],(* PlotRange->{II,{-.05,.05}},Axes->{True,False},*)PlotMarkers->{"✶",18},PlotStyle->{Red,PointSize[.02],Opacity[.8]}(*,Ticks->{{{0,"0"},{cenSymDir[[1]],"c1"},{cenSymDir[[2]],"c2"},{cenSymDir[[3]],"c3"},{1,"1"}},None},
Ticks->{{{0," "},{1," "}},None},AspectRatio->Automatic *)];
Show[slB,slcSym,Ticks->{{{0," "},{1," "}},None},ImageSize->300]

Initial approximation for k=3

k=3;
cFunk0[var_]:=Compile[var,Evaluate[Sum[Min[Table[Norm[Table[(Subscript[x, j]),{j,1+(jk-1)n,1+jk n-1}]-B[[i]]]^2,{jk, 1,k}]],{i,m}]]];
FFCompile=cFunk0[Table[{(Subscript[x, i]),_Real},{i,n k}]];

(* SymDIRECT *)
Print[Style["SymDIRECT:",16,Blue]];
dom=Table[II,{k}];
tSym=Timing[
minSym=SymDIRECT[FFCompile,dom,DMin->10^(-4),MaxIter->20,Eps->10^(-4)]
];
cenSymDir=minSym[[1]];
Print[Style["SymDIRECT Centers: ",Green,16],cenSymDir];
Print["F= ",minSym[[2]],"; IT = ",iter,"; CPU = ",tSym[[1]],"; Nf = ",Nf];
slcSym=ListPlot[Table[{cenSymDir[[j]],0},{j, k}],PlotMarkers->{"✶",15},PlotStyle->{Red,PointSize[.02],Opacity[.8]}(*,PlotRange->{II,{-.05,.05}},Axes->{True,False},Ticks->{{{0,"0"},{cenSymDir[[1]],"c1"},{cenSymDir[[2]],"c2"},{cenSymDir[[3]],"c3"},{1,"1"}},None},AspectRatio->Automatic*)];
Show[slB,slcSym,ImageSize->300]

One sphere

r0=cenSymDir[[1]]; Print["r0 = ", r0]; (* The bigest radius - the shallowest *)
lin1[t_]:=InverseFunction[lin][t];

Do[
sph=SphRep[A,C0,r0,-.0005,aa];
Print[Style["Approach: ",Blue,14],sph[[1]],"\nr = ",sph[[2]],";  F=",sph[[3]],"; CPU=",sph[[4]]]
,{aa,3}]

Figure

r0=cenSymDir[[1]]; Print["r0 = ", r0]; 
sph=SphRep[A,C0,r0,-.0005,1];  (* LAD distance "1" is choosen *)
Print[Style["Approach: ",Blue,14],sph[[1]],"\nr = ",sph[[2]],";  F=",sph[[3]],"; CPU=",sph[[4]]]
AG=Table[{-1,1},3];
slSph1=Graphics3D[{LightBlue,Opacity[.1],Sphere[C0,1]}];
slSph=Graphics3D[{LightBlue,Opacity[.3],Sphere[C0,sph[[2]]]}];
slc=ListPointPlot3D[{C0}, PlotStyle->{Brown,PointSize[.06],Opacity[.8]}];
Show[slA,slSph,slSph1,
PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},ImageSize->Small]
Print["Real depth = ",Round[-lin1[sph[[2]]]]"km"];

Three spheres

AG=Table[{-1.,1.},3];
Circ0=Table[{C0,cenSymDir[[j]]},{j,3}]
kcs=Timing[KCS[A,Circ0,3,1]];
Print["CPU(KCS) = ",kcs[[1]]," seconds"]

Analysis of results in the Unit Sphere model

cenB=Sort[kcs[[2,2,All,2]]];
Print[Style["Radii of spheres inside Unit Sphere: ",Blue,14],cenB]
ed={lin[min],Mean[{cenB[[1]],cenB[[2]]}],Mean[{cenB[[2]],cenB[[3]]}],lin[max]};
Print[Style["Radii of boundary spheres",Blue,14],Style["(they delimit clusters of transformed depths within the Unit Sphere):\n",Gray,12],ed]
Print[Style["Spherical regions inside Unit Sphere with the number of earthquake events:\n Deep: ",Blue,14],{ed[[1]],ed[[2]]}," | ",Length[kcs[[2,1,3]]],Style["\n Intermediate: ",Blue,14],{ed[[2]],ed[[3]]}," | ",Length[kcs[[2,1,2]]],Style["\n Shallow: ",Blue,14],{ed[[3]],ed[[4]]}," | ",Length[kcs[[2,1,1]]]]

Analysis of real results

Print[Style["Depths (in kilometers) of cluster centers in the Geoid: ",Blue,14],Sort[-lin1[cenB]]]
ed1=lin1[ed];
 (* Print[Style["Deepest earthquake and depths separating Deep from Intermediate and Intermediate from Shallow earthquakes: ",Blue,14](* Style["(they delimit clusters of transformed depths within the Unit Sphere):\n" Gray,12]*),-lin1[ed]] *) 
Print[Style["Depth ranges (in kilometers) with numbers of earthquake events:\n Shallow: ",Blue,14],{Round[-ed1[[4]]],Round[-ed1[[3]]]}," | ",Length[kcs[[2,1,1]]],Style["\n Intermediate: ",Blue,14],{Round[-ed1[[3]]],Round[-ed1[[2]]]}," | ",Length[kcs[[2,1,2]]],Style["\n Deep: ",Blue,14],{Round[-ed1[[2]]],Round[-ed1[[1]]]}," | ",Length[kcs[[2,1,3]]]]