Modules + Input

Mahalanobis distance - like functions

nn=2;
dm[x_,y_,S_]:=(x-y).Inverse[S].(x-y);
dM[x_,y_,S_]:=(Det[S]^(1/nn))(x-y).Inverse[S].(x-y);

Weight functions

wf[x_,y_,ksi_,eta_,d_]:=If[(((x-ksi))^2)+(((y-eta))^2)<(d^2),(((1-((Abs[((Sqrt[(((x-ksi))^2)+(((y-eta))^2)])/d)])^3)))^3),0];

LAD and OLS - Local Approximant

(*  
data: data of the form {xx,yy,zz};
nodes: points of the form {xi, eta};
sigma: parameters of the form Table[sig[[i]],{i,Length[data]}];
appr: 1 (LAD) or 2 (LS);
par: parameters of the form Table[{par0[[j]],par1[[j]],par2[[j]]},{j,Lqngrh[nodes]}];
If Ind≠0, printaju se vrijednosti parametara i funkcije G;
*)
LocalAppr[data_,nodes_,sigma_,appr_,Ind_]:=Module[{k=Length[nodes], m=Length[data],p0,p1,p2,v,min,par,time},
par={};
If[appr==1,Print[Style["LAD local approximants:",18,Blue]],Print[Style["OLS local approximants:",18,Blue]]];
time=Timing[
Do[
v=Table[wf[nodes[[j,1]],nodes[[j,2]],data[[ii,1]],data[[ii,2]],sigma[[ii]]]
,{ii,m}];
G[p0_,p1_,p2_]:=Sum[v[[i]]Abs[data[[i,3]]-p0-p1  data[[i,1]]-p2  data[[i,2]]]^appr,{i,m}];
min=NMinimize[G[p0,p1,p2],{p0,p1,p2}];
par=Append[par,{p0,p1,p2}/.min[[2]]];
If[Ind≠0,Print[j,". parameters: ",Last[par],";  G = ",min[[1]]]];
,{j,k}]];
Print["time = ",time[[1]]];
par
]

TTR

TTR[clu_]:=Module[{kappa=(365. 24 60 60),data=Table[0,{i,Length[clu]}],time=Table[0,{i,Length[clu]}],year=Table[0,{i,Length[clu]}]},
Do[
data[[i]]=Round[Take[clu[[i]],{1,6}]];
time[[i]]=AbsoluteTime[Round[Take[clu[[i]],{1,6}]]]/kappa;
year[[i]]=Floor[AbsoluteTime[Round[Take[clu[[i]],{1,6}]]]/kappa];
,{i,Length[clu]}];
SortBy[Table[{time[[i]],data[[i]],year[[i]],{clu[[i,8]],clu[[i,7]]},clu[[i,9]],clu[[i,10]]},{i,Length[clu]}],First]
]

Input

SetDirectory[NotebookDirectory[]]; 
asd=mapa4=ReadList["C:/Rudolf/AB-STRUKA/1_RADOVI/B1_CLUSTERS/B1_PACO/Paco-2023/nb-IP3/MAPart4D.txt"];
asd4=ReadList["C:/Rudolf/AB-STRUKA/1_RADOVI/B1_CLUSTERS/B1_PACO/Paco-2023/nb-IP3/MAPart4.txt"];
asD=ReadList["C:/Rudolf/AB-STRUKA/1_RADOVI/B1_CLUSTERS/B1_PACO/Paco-2023/nb-IP3/IP35D.txt"];
IP0=ReadList["C:/Rudolf/AB-STRUKA/1_RADOVI/B1_CLUSTERS/B1_PACO/Paco-2019/DATA/IP0.txt"];

Cluster selection

(* Nakon pokretanja "Izabrani klaster", ne daje MAPart4D s contours - ovo treba provjeriti ! *)
jj=1;
data=asd[[jj,All,3]]; 
Print[Style["CLUSTER ",18,Hue[1]], Style[jj,18,Red]," - Length[data] = ",m=Length[data]];
data=SortBy[data,First];
Print["First: : ",First[data]]
Print["Last: ",Last[data]]
Print["Range of Magnitudes: ",{Min[data[[All,10]]],Max[data[[All,10]]]}]
slIP=ListPlot[IP0,Joined->True,AxesOrigin->{-12,34},PlotStyle->{Thickness[.004],Orange},ImageSize->Small];
					(* Structure according magnitudes *)
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]
							(* Control *)
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]

				(* Data construction of the selected cluster;
				    Depths are defined as negative numbers;
					POD: 3D data set;
					pod: 2D subset of POD *)
n=3; h=1;
POD=Table[{data[[i,8]],data[[i,7]], -data[[i,9]]},{i,Length[data]}];
w=Table[data[[i,10]],{i,Length[data]}]; W=Total[w];
AG=Table[{Min[POD[[All,i]]],Max[POD[[All,i]]]},{i,n}];
Print["RECTANGLE: ",AG,"; |POD|=",Length[POD]]
pod=Table[{data[[i,8]],data[[i,7]]},{i,Length[data]}];
slpod=ListPlot[pod,PlotStyle->{PointSize[.01]}];
ag=Table[{Min[POD[[All,i]]]-h,Max[POD[[All,i]]]+h},{i,n-1}];
cA=Mean[WeightedData[pod,w]];
slcA=ListPlot[{cA},PlotMarkers->{"✶",18},PlotStyle->{Orange,PointSize[.05]}];
Print["rectangle: ",ag,"; |pod|=",Length[pod],";  cA = ",cA]
								(* FIGURE *)
Scov=Sum[w[[s]]*KroneckerProduct[(pod[[s]]-cA),(pod[[s]]-cA)],{s,Length[pod]}]/W;
tab=Table[dM[cA,pod[[s]],Scov],{s,Length[pod]}];
rad=Quantile[tab,95/100];
slcon=ContourPlot[dM[cA,{x,y},Scov]==rad,{x,ag[[1,1]],ag[[1,2]]},{y,ag[[2,1]],ag[[2,2]]}];
sl1a=Show[slcon/.Line[x_]->{LightBlue,EdgeForm@Thin,FaceForm[Opacity@.5],Polygon[x]},PlotRange->ag];
sljj=Show[slIP,slcon,sl1a,slpod,slcA, AspectRatio->Automatic, PlotRange->ag,AxesOrigin->{-12,34},GridLines->Automatic,ImageSize->300 ];

									(* Nodes *)
nodes={cA};
Do[
If[rr==1,ka={0,1,-1},ka={0,.4,1,-.4,-1,100}];
Do[
nsol=NSolve[{dM[cA,{x,y},Scov]==rr, y-cA[[2]]==kk(x-cA[[1]])},{x,y}];
P1={x,y}/.nsol[[1]];
P2={x,y}/.nsol[[2]];
nodes=Append[nodes,P1];nodes=Append[nodes,P2]
,{kk,ka}]
,{rr,{0.2 rad,.7rad}}]
							(* Figures *)
slData=ListPointPlot3D[POD,ViewPoint->{1.2,2.4,1.8},  PlotStyle->{{PointSize[.015],GrayLevel[.2]}}];
slnodes=ListPlot[nodes,PlotStyle->{Orange, PointSize[.02]},AspectRatio->Automatic];
sl2=Show[sljj, slnodes];
GraphicsGrid[{{slnodes,sl2, slData}},ImageSize->700]
Print["Površina = ",Pi Det[Scov],";    |pod| = ",Length[pod],";    |nodes| = ",Length[nodes],";    Gustoća = ",Length[pod]/(Pi Det[Scov])]

OLS or LAD Local approximants

(* za appr=1 dobiva se LAD-Loc.Appr., a za appr=2 OLS-Loc.Appr. *)
appr=2;
sigma=Table[5,{i,Length[POD]}];
loc=LocalAppr[POD,nodes,sigma,appr,1];

Convex combination of local approximants (takes approx 0.5 min)

(*
jj=1; {x,-13,-5},{y,34.5,39};sigma=5;
jj=2: {x,-8,3},{y,34.5,41.5}; sigma=7;
jj=3: {x,-6,-1.3},{y,33.7,37.3};sigma=4;
jj=4: {x,-1.5,6.5},{y,35.,37.8}; sigma=2;
jj=5: {x,-11,-4},{y,38.5,45.5}; sigma=4;
jj=6: {x,-5.5,5.5},{y,41.2,44.6}; sigma=4;
*)
k=Length[nodes];
sigma=Table[5,{i,Length[POD]}];
lokalni[ii_,xi_,yi_]:=loc[[ii,1]]+loc[[ii,2]] xi+loc[[ii,3]] yi;
ufunk[ii_,xi_,yi_]:=wf[xi,yi,nodes[[ii]][[1]],nodes[[ii]][[2]],sigma[[ii]]];
omega[ii_,xi_,yi_]:=((ufunk[ii,xi,yi])/((∑_(jj=1)^k)ufunk[jj,xi,yi]));
Sploha[xi_,yi_]:=(∑_(jj=1)^k)lokalni[jj,xi,yi] omega[jj,xi,yi];
					(* FIGURES *)
If[appr==1,Print[Style["LAD local approximants ",18,Blue]," for jj = ",jj],
Print[Style["OLS local approximants ",18,Blue]," for jj = ",jj]];
slfi=Plot3D[Sploha[x,y],{x,-12,5},{y,34.,44}];
slCon=ContourPlot[Sploha[x,y],{x,-13,-5},{y,34.5,39}];
slCon2=Show[slCon,sljj];
Print[GraphicsGrid[{{slfi,slCon2}},ImageSize->500]]

The whole area

Print["m0 = ",m0=Length[asD]]
Print["First: : ",First[asD]]
Print["Last: ",Last[asD]]
Print["Magnitudes: ",{Min[asD[[All,10]]],Max[asD[[All,10]]]}]
boje={1,1.5,2.,3.,4,5,6.,7,8.5,9,10}/10;
boje={7,8,10,2,3,5}/10//N
Table[Graphics[{Hue[boje[[j]]],Disk[]},ImageSize->20],{j,Length[boje]}]
slIP=ListPlot[IP0,Joined->True,AxesOrigin->{-12,34},PlotStyle->{Thickness[.008],Red,Opacity[.9]},ImageSize->Small]

Data construction

n=3; h=0; k=4;
POD=Table[{asD[[i,8]],asD[[i,7]], -asD[[i,9]]},{i,Length[asD]}];
w=Table[asD[[i,10]],{i,Length[asD]}]; W=Total[w];
A=Table[{asD[[i,8]],asD[[i,7]]},{i,Length[asD]}];
AG=Table[{Min[POD[[All,i]]],Max[POD[[All,i]]]},{i,n}];
ag=Table[{Min[POD[[All,i]]]-h,Max[POD[[All,i]]]+h},{i,n-1}];
Print["RECTANGLE: ",AG,"; |POD|=",Length[POD]]
pod=Table[{asD[[i,8]],asD[[i,7]]},{i,Length[asD]}];
cA=Mean[WeightedData[pod,w]];
slcA=ListPlot[{cA},PlotMarkers->{"✶",18},PlotStyle->{Orange,PointSize[.05]}];
Print["rectangle: ",ag,"; |pod|=",Length[pod],";  cA = ",cA]
							(* Transformations *)
dd=Table[ag[[i,2]]-ag[[i,1]],{i,n-1}];
MD=DiagonalMatrix[1/dd];
MD1=DiagonalMatrix[dd];
vu=Table[Min[pod[[All,i]]],{i,n-1}];
TT[x_]:=MD.(x-vu); 
TT1[x_]:=MD1.x+vu; 
Print["Example: \nai=pod[[300]] = ",ai=pod[[3]], "\nbi=TT[ai] = ",bi=TT[ai],"\nTT1[bi] = ",TT1[bi]  ]
								(* FIGURE *)
pod4=pod5=pod6={};
Do[
If[w[[i]]≤4,pod4=Append[pod4,A[[i]]],
If[w[[i]]≤5,pod5=Append[pod5,A[[i]]],pod6=Append[pod6,A[[i]]]]]
,{i,Length[asD]}]
slpod=Show[slIP,ListPlot[pod4,PlotStyle->{Black,Opacity[.8], PointSize[.005]},Axes->True,AspectRatio->Automatic],ListPlot[pod5,PlotStyle->{Black,Opacity[.5],PointSize[.01]},Axes->True,AspectRatio->Automatic],ListPlot[pod6,PlotStyle->{Black,Opacity[.3],PointSize[.015]},Axes->True,AspectRatio->Automatic],PlotRange->ag,GridLines->Automatic,AxesOrigin->{ag[[1,1]],ag[[2,1]]},ImageSize->350]

Knotes

nodes={};
Do[
Do[
nodes=Append[nodes,{i,j}]
,{i,-11,4,1}]
,{j,35,43,1}]
Length[nodes]
					(* Figures *)
slData=ListPointPlot3D[POD, PlotStyle->{{PointSize[.01],GrayLevel[.2]}}];
slnodes=ListPlot[nodes,PlotStyle->{Orange, PointSize[.015]},AxesOrigin->{ag[[1,1]],ag[[2,1]]},GridLines->Automatic,AspectRatio->Automatic];
GraphicsGrid[{{slnodes,Show[slIP, slnodes,GridLines->Automatic,Frame->True], slData}},ImageSize->700]
Print["|pod| = ",Length[pod],";    |nodes| = ",Length[nodes]]

MAPart4D

j=1;asd[[1,j,All,2]]

PIP=cc=Scov=Table[0,{j,k}];
Do[
PIP[[j]]=asd[[1,j,All,2]];
cc[[j]]=Mean[WeightedData[asd[[1,j,All,2]],asd[[1,j,All,1]]]];
Scov[[j]]=Sum[asd[[1,j,s,1]]*KroneckerProduct[(asd[[1,j,s,2]]-cc[[j]]),(asd[[1,j,s,2]]-cc[[j]])],{s,Length[asd[[1,j]]]}]/Total[asd[[1,j,All,1]]];
,{j,k}]
slPI=ListPlot[Table[Table[TT1[PIP[[j,s]]],{s,Length[PIP[[j]]]}],{j,k}],PlotStyle->Table[{PointSize[.006],Hue[boje[[j]]]},{j,k}]];
Do[
mj=Length[PIP[[j]]];
tab=Table[dM[TT1[cc[[j]]],TT1[PIP[[j,s]]],MD1.Scov[[j]].MD1],{s,mj}];
rad=Quantile[tab,95/100];
slcon[j]=ContourPlot[dM[TT1[cc[[j]]],{x,y},MD1.Scov[[j]].MD1]==rad,{x,ag[[1,1]]-1, ag[[1,2]]+1},{y,ag[[2,1]]-1,ag[[2,2]]+1},ContourStyle->{Thickness[.008],Gray}]
	,{j,k}];
sl=Show[Table[slcon[j],{j,k}]];
sl2=Show[slPI,sl,slIP,AspectRatio->Automatic,PlotRange->ag,GridLines->Automatic,AxesOrigin->{ag[[1,1]],ag[[2,1]]},ImageSize->350]

OLS or LAD Local approximants

(* for appr=1 is obtained LAD-Loc.Appr., and for appr=2 OLS-Loc.Appr. *)
appr=1;
sigma=Table[5,{i,Length[POD]}];
loc=LocalAppr[POD,nodes,sigma,appr,1];

Convex combination of local approximants (takes approx 0.5 min)

(*
sigma=5;
*)
k=Length[nodes];
sigma=Table[5,{i,Length[POD]}];
lokalni[ii_,xi_,yi_]:=loc[[ii,1]]+loc[[ii,2]] xi+loc[[ii,3]] yi;
ufunk[ii_,xi_,yi_]:=wf[xi,yi,nodes[[ii]][[1]],nodes[[ii]][[2]],sigma[[ii]]];
omega[ii_,xi_,yi_]:=((ufunk[ii,xi,yi])/((∑_(jj=1)^k)ufunk[jj,xi,yi]));
Sploha[xi_,yi_]:=(∑_(jj=1)^k)lokalni[jj,xi,yi] omega[jj,xi,yi];
					(* FIGURES *)
If[appr==1,Print[Style["LAD local approximants ",18,Blue]],
Print[Style["OLS local approximants ",18,Blue]]];
slfi=Plot3D[Sploha[x,y],{x,-12,5},{y,34.,44},ColorFunction->"RustTones",PlotStyle->{Opacity[0.7]},PlotRange->All]
slCon=ContourPlot[Sploha[x,y],{x,-12,5},{y,34.,44},PlotTheme->"Web",ContourStyle->None,PlotRange->All,PlotLegends->Automatic]
(*Print[GraphicsGrid[{{slfi,slCon}},ImageSize->500]]
Print["Time for slfi = ",time1[[1]], ";  Time for slCon = ",time2[[1]]]*)

Image construction

If[appr==1,Print[Style["LAD local approximants ",18,Blue]],
Print[Style["OLS local approximants ",18,Blue]]];Print[Show[slCon,slpod,sl,GridLines->Automatic,AspectRatio->.7]]
Print[Show[slfi,AspectRatio->.7]]

pod4=pod5=pod6={};
Do[
If[w[[i]]≤4,pod4=Append[pod4,POD[[i]]],
If[w[[i]]≤5,pod5=Append[pod5,POD[[i]]],pod6=Append[pod6,A[[i]]]]]
,{i,Length[asD]}]
slData=Show[ListPointPlot3D[pod4,PlotStyle->{Opacity[.8], PointSize[.005],Black}],ListPointPlot3D[pod5,PlotStyle->{Opacity[.7],PointSize[.01],Black}],ListPointPlot3D[pod6,PlotStyle->{Opacity[.5],PointSize[.015],Black}],PlotRange->AG];Show[slData,slfi, PlotRange->{{-11.91`,4.972`},{34.05`,43.999`},{-45,0}},ImageSize->350]

Strong earthquakes in (Subscript[π, 1])

jj=1;ASD=asd4[[1,jj,All,3]];
nn=2;
Print[Style["Dobar potres\n",16,Blue],ASD[[159]]]
Print["m0 = ",m0=Length[ASD],Style[" (broj svih potresa)",12,Blue], " u ",jj,". tom klasteru"]
Print["First: : ",First[ASD]]
Print["Last: ",Last[ASD]]
Print["Magnitudes: ",{Min[ASD[[All,10]]],Max[ASD[[All,10]]]}]
slIP=ListPlot[IP0,Joined->True,AxesOrigin->{-12,34},PlotStyle->{Thickness[.004],Orange},ImageSize->Small];
Print[Style["The image is made with all earthquakes",20,Blue]]

Podaci

pod3=pod4=pod5=pod6=0;
Print["m = ",m=Length[ASD],Style[" (Broj izabranih potresa s poznatom ili i nepoznatom dubinom )",12,Blue]]
(* Struktura po magnitudi *)
Print[Style["Potresi magnitude Mg ≥ 5:",16,Blue]]
Do[
If[ASD[[i,10]]<  4,pod3=pod3+1,
If[ASD[[i,10]]<5,pod4=pod4+1,pod5=pod5+1;Print[i,". ",ASD[[i]]]]]
,{i,Length[ASD]}];
Print["pod3=",pod3,"; pod4=",pod4,"; pod5=",pod5,"; podAll= ",pod3+pod4+pod5]
Nr=0;
Do[
If[6>ASD[[i,10]]≥5,Nr=Nr+1]
,{i,Length[ASD]}]
Print["Broj potresa 6 > Mg >= 5: ",Nr,"\nBroj potresa Mg>=6: ",pod5-Nr]
						(* Data and Figure *)
POD=Table[{ASD[[i,10]],{ASD[[i,8]],ASD[[i,7]]},i},{i,Length[ASD]}];
nn=2;
w=POD[[All,1]]; W=Total[w];
AG=Table[{Min[POD[[All,2,i]]],Max[POD[[All,2,i]]]},{i,nn}];
cA=Mean[WeightedData[POD[[All,2]],POD[[All,1]]]];
Print["Rectangle: ",AG,"; centroid(A) = ",cA; "; |POD|=",Length[POD]]
Print["Partition ",jj,": {M3, M4, M5}: ",Style[{pod3, pod4, pod5},16,Blue]]
					(* Position of selected cluster *)
Show[slIP,ListPlot[POD[[All,2]],PlotRange->All],PlotStyle->{PointSize[.5]}, AspectRatio->Automatic, AxesOrigin->{-12,34},GridLines->Automatic,ImageSize->Small ]
(* Izabrani klsater *)
Scov=Sum[POD[[s,1]]*KroneckerProduct[(POD[[All,2]][[s]]-cA),(POD[[All,2]][[s]]-cA)],{s,Length[POD]}]/W;
tab=Table[dM[cA,POD[[s,2]],Scov],{s,Length[POD]}];
rad=Quantile[tab,90/100];
AG2=Table[{Floor[Min[POD[[All,2,i]]]],Ceiling[Max[POD[[All,2,i]]]]},{i,nn}];
slcon=ContourPlot[dM[cA,{x,y},Scov]==rad,{x,AG2[[1,1]],AG2[[1,2]]},{y,AG2[[2,1]],AG2[[2,2]]}];
sl1a=Show[slcon/.Line[x_]->{LightBlue,EdgeForm@Thin,FaceForm[Opacity@.5],Polygon[x]},PlotRange->AG];
sl159=Show[slIP,slCon,sl1a,ListPlot[POD[[All,2]],PlotStyle->{PointSize[.01]}], AspectRatio->Automatic, PlotRange->AG2,AxesOrigin->{-12,34},GridLines->Automatic,ImageSize->200 ]
(* Transformation to[0,1] *)
C0={0,0,0}; Print["C0 = ",C0];
slSph1=Graphics3D[{LightBlue,Opacity[.1],Sphere[C0,1]}]; 
phi=ASD[[All,7]]; phi=Pi phi/180.;
th=ASD[[All,8]]; th=Pi th/180.;
depth=-ASD[[All,9]]; 
 min=Min[depth]; (* the deepest *)
 max=Max[depth]; (* the shallowest *)
Print["{the deepest, the shallowest}: ",{min, max}," -> [0,1]"]
lin[t_]:=-t/min+1.
lin1[t_]=InverseFunction[lin][t];
B=lin[depth];
II={0,1};
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]];
						(* Are there in B outliers? *)
dubine=-depth;
med=Median[dubine]+.01;
mad=1.483Median[Abs[dubine-med]]; Print["MAD=",mad];
dubinetil=Abs[dubine-med]/mad;
Outliers={};
Do[
If[dubinetil[[i]]>2.5,Outliers=Append[Outliers,{i,dubine[[i]]}]]
,{i,Length[dubine]}]
Outliers
Print["Od ukupno ",m," podataka postoji ",Length[Outliers]," outliers, tj., ",100.Length[Outliers]/m,"% outliers"]