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

SymDIRECT

Modul for solving the global optimization problem for a simmetric 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]} ]]
];

Weighted k -  Means

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

Sphere Representative

(* LAD(approach=1),  TLS(approach=2); Algebraic (approach=3); *)
SphRep[A_,r0_,approach_]:=Module[{C0={0,0,0},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
]

Data

(* MAPart4D is 4-partitiona of the data with known depths *)
SetDirectory[NotebookDirectory[]]; 
asd=mapa4=ReadList["C:/Rudolf/AB-STRUKA/1_RADOVI/B1_CLUSTERS/B1_PACO/Paco-2023/nb-IP3/MAPart4D.txt"];
IP0=ReadList["C:/Rudolf/AB-STRUKA/1_RADOVI/B1_CLUSTERS/B1_PACO/Paco-2019/DATA/IP0.txt"];
Print["Total data: ",Table[Length[asd[[1,j]]],{j,4}]," = ",Total[Table[Length[asd[[1,j]]],{j,4}]]]

Cluster selection

jj=3;data=asd[[1,jj,All,3]]; 
Print[Style["Cluster ",18,Blue], Style[jj,18,Blue]," - Length[data] = ",m=Length[data]];
data=SortBy[data,First];
Print["First: : ",First[data]]
Print["Last: ",Last[data]]
Print["Magnitudes: ",{Min[data[[All,10]]],Max[data[[All,10]]]}]

(* New dataset of earthquakes with known depths (without depth=0) *)
pod3=pod4=pod5=pod6=0;
Do[
If[data[[i,10]]<  4,pod3=pod3+1,
If[data[[i,10]]<5,pod4=pod4+1,pod5=pod5+1]]
,{i,Length[data]}];
Print["pod3=",pod3,"; pod4=",pod4,"; pod5=",pod5,"; podAll= ",pod3+pod4+pod5]
(* Kontrola *)
Nr=0;
Do[
If[6>data[[i,10]]≥5,Nr=Nr+1]
,{i,Length[data]}]
Print["Broj potresa 6 > Mg >= 5: ",Nr,"\nBroj potresa Mg>=6: ",pod5-Nr]

(* Transformation to[0,1] *)
C0={0,0,0}; Print["C0 = ",C0];
phi=data[[All,7]]; phi=Pi phi/180.;
th=data[[All,8]]; th=Pi th/180.;
depth=data[[All,9]]; 
 min=Min[depth]; (* the deepest *)
 max=Max[depth]; (* the shallowest *)
Print["{the deepest, the shallowest}: ",{min, max}," -> [0,1]"]
lin[t_]:=-t/max+1.
lin1[t_]=InverseFunction[lin][t];
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]];

Spheres

(* One sphere *)
Print[Style["One sphere:",18,Blue]]
Print[Style["Approach 1:  r = ",Blue,14],r1=Median[B],"; Actual depth = ", lin1[r1],"km"]
Print[Style["Approach 2:  r = ",Blue,14],r2=Mean[B],"; Actual depth = ", lin1[r2]," km"]
Print[Style["Approach A:  r = ",Blue,14],rA=Sqrt[Mean[Table[B[[s]]^2,{s,m}]]],"; Actual depth = ",lin1[rA]," km"]
								(* Figure *)
r0=r1; Print["r0 = ", r0]; 
C0={0,0,0};
sph=SphRep[A,r0,1];  (* Birao sam LAD distance "1" *)
Print[Style["Approach: ",Blue,14],sph[[1]],"\nr = ",sph[[2]],";  F=",sph[[3]],"; CPU=",sph[[4]]]
AG=Table[{-1,1},3];
slcSym=ListPlot[{{r0,0}},PlotMarkers->{"✶",15},PlotStyle->{Red,PointSize[.02],Opacity[.8]}];
Show[slB,slcSym,ImageSize->300]
Print["Stvarna dubina = ",lin1[sph[[2]]]];

							(* Three sphere *)
Print[Style["Three sphere:",18,Blue]]
				(* Compiled objective function *)
k=3; n=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}]]; 
				(* EXECUTION SimDIRECT and WKMeans *)
Print[Style["SymDIRECT:",16,Blue]];
dom=Table[II,{k}];
tSym=Timing[
minSym=SymDIRECT[FFCompile,dom,DMin->10^(-4),MaxIter->20,Eps->10^(-4)]
];
Print["SymDIRECT Centers: ",Sort[minSym[[1]]]];
Print["F= ",minSym[[2]],"; IT = ",iter,"; CPU = ",tSym[[1]],"; Nf = ",Nf];
						(* WKMeans *)
tWKM=Timing[sol=WKMeans[B,w,minSym[[1]]]];
Print["\nk-means:  F= ",sol[[3]],"; IT = ",sol[[4]],"; CPU = ",tWKM[[1]]];	

slcSym=ListPlot[Table[{sol[[2,j]],0},{j, k}],PlotRange->{II,{-.05,.05}},Axes->{True,False},PlotMarkers->{"✶",18},PlotStyle->{Red,PointSize[.02],Opacity[.5]},Ticks->{{{0,"0"},{sol[[2,1]],"c3"},{sol[[2,2]],"c2"},{sol[[2,3]],"c1"},{1,"1"}},None},AspectRatio->Automatic];
Show[slcSym,slB,ImageSize->500]
Print["k-means Centers: ",cenB=Sort[minSym[[1]]]];

Analysis

(* Analiysis of results inside Unit Sphere *)
Print[Style["Analiysis of results inside Unit Sphere: ",18,Red]]
lin1[x_]:=InverseFunction[lin][x]
cenB=Sort[sol[[2]]];
Print[Style["Radia of spheres inside Unit Sphere: ",Blue,14],cenB]
ed={lin[min],Mean[{cenB[[2]],cenB[[3]]}],Mean[{cenB[[1]],cenB[[2]]}],lin[max]};
Print[Style["Radia of boundary spheres",Blue,14],Style["(they delimit clusters of transformed depths within the Unit Sphere):\n",Gray,12],ed]
Print[Style["Areas inside Unit Sphere (sphere garland) with earthquake events:\nNajpliće: ",Blue,14],{ed[[1]],ed[[2]]},"|",Length[sol[[1,1]]],Style["\nSredina: ",Blue,14],{ed[[2]],ed[[3]]},"|",Length[sol[[1,2]]],Style["\nNajdublje: ",Blue,14],{ed[[3]],ed[[4]]},"|",Length[sol[[1,3]]]]

					(* Analiysis of real results *)
Print[Style["Analiysis of real results: ",18,Red]]
Print[Style["Radia of spheres inside Real Sphere: ",Blue,14],lin1[cenB]]
ed1=lin1[ed];
Print[Style["Radia of boundary spheres",Blue,14],Style["(they delimit clusters of transformed depths within the Unit Sphere):\n",Gray,12],lin1[ed]]
Print[Style["Areas inside Real Sphere (sphere garland) with earthquake events:\nShallow: ",Blue,14],{ed1[[1]],ed1[[2]]},"|",Length[sol[[1,1]]],Style["\nIntermediate: ",Blue,14],{ed1[[2]],ed1[[3]]},"|",Length[sol[[1,2]]],Style["\nDeep: ",Blue,14],{ed1[[3]],ed1[[4]]},"|",Length[sol[[1,3]]]]

							(* Figure *)
cenDepth=lin1[cenB]; Print["Centri = ",cenDepth];slDepth=ListPlot[Table[{depth[[j]],0},{j,m}],PlotRange->{{min,max},{-.05,.05}},Axes->{True,False},Ticks->{{{min,"min"},{cenDepth[[3]],"c3"},{cenDepth[[2]],"c2"},{cenDepth[[1]],"c1"},{max,"Max"}},None},PlotStyle->{Red,PointSize[.002]},ImageSize->500];
slcDepth=ListPlot[{{cenDepth[[1]],0},{cenDepth[[2]],0},{cenDepth[[3]],0}},PlotRange->{{min,max},{-.05,.05}},Axes->{True,False},PlotMarkers->{"✶",20},PlotStyle->{Blue,PointSize[.02],Opacity[.5]},AspectRatio->Automatic];
Show[slDepth,slcDepth,ImageSize->700]