Seismogenic zoning IP
29-5-2023

Modules

Needs["ComputationalGeometry`"]

Squared Mahalanobis distance

n=2;
p=.7;  kappa=Quantile[ChiSquareDistribution[n],p];
dM[xpod_,ypod_,S_]:=(Det[S]^(1/n))(xpod-ypod).Inverse[S].(xpod-ypod);

Mahalanobis objective function

MWF[PI_,c_,Scov_]:=Module[{},
Return[Sum[PI[[j,s,1]]*dM[PI[[j,s,2]],c[[j]],Scov[[j]]],{j,Length[PI]} ,{s,Length[PI[[j]]]}] ]
];

Fmin[POD_,c_,Scov_]:=Module[{},
Sum[POD[[i,1]]Min[Table[dM[c[[j]],POD[[i,2]],Scov[[j]]],{j,Length[c]}]],{i,Length[POD]}]
]

Area index (for the data without weights)

(* Correction on p-ellipsoids *)
Clear[MArea]
MArea[PI_,z_,Scov_]:=Module[{k=Length[z],PP},
Ell=Table[Ellipsoid[z[[j]],kappa  Scov[[j]]],{j,k}];
					(* Subpartition *)
PP=Table[{},{j,k}];
Do[
Do[
If[PI[[j,s,2]]∈Ell[[j]],PP[[j]]=Append[PP[[j]],PI[[j,s]] ]]
,{s,Length[PI[[j]]]}]
,{j,k}];
(*  Index on subpartition PP *)
Print["areas: ",Table[Pi Det[ Scov[[j]]]//N,{j,k}],"; W: ",Table[Total[PP[[j,All,1]]],{j,k}]];
Sum[Det[kappa Scov[[j]]]/Total[PP[[j,All,1]]],{j,k}]
]

(* Stari index *)
Clear[MArea]
MArea[PI_,z_,Scov_]:=Module[{},
Print["areas: ",Table[Pi Det[ Scov[[j]]]//N,{j,Length[PI]}],"; W: ",Table[p Length[PI[[j]]],{j,Length[PI]}]];
Sum[Det[  Scov[[j]]]/(Length[PI[[j]]]),{j,Length[PI]}]
]

MMD Index

(* MMD Index;
Ako je Ind2≠0, iscrtacvaju se slike parova elipsoida *)
Clear[MMD]
MMD[PI_,z_,Scov_,Ind2_:0]:=Module[{k=Length[PI],P1,P2},
DDD=Table[0,{r,k},{s,k}];
PI0=Table[PI[[j,All,2]],{j,k}];
Ell=Table[Ellipsoid[z[[j]],kappa  Scov[[j]]],{j,k}];
	 (* M-udaljenosti (najprije se testira je li presjek prazan) *)
Do[
Do[
If[IntersectTest[Ell[[r,1]],Ell[[r,2]],Ell[[s,1]],Ell[[s,2]]],
min=FindMinimum[{EuclideanDistance[x,y],{x∈Ell[[r]],y∈Ell[[s]]}}, {x,y}];If[Ind2≠ 0,Print[min];{P1,P2}={x/.min[[2]],y/.min[[2]]};
		];
DDD[[r,s]]=min[[1]],
{P1,P2}={Ell[[r,1]],Ell[[r,1]]};DDD[[r,s]]=0
	];
If[Ind2≠ 0, 
If[n==2,
Print[Show[
Graphics[{boja[[r]],EdgeForm[Directive[Thickness[.004],boja[[r]]]],Opacity[0.075],Ell[[r]]}],
Graphics[{boja[[s]],EdgeForm[Directive[Thickness[.004],boja[[s]]]],Opacity[0.075],Ell[[s]]}],
ListPlot[PI0,PlotStyle->Table[{PointSize[.008],boja[[j]]},{j,Length[cen]}]], 
ListPlot[{P1,P2}, PlotStyle->{Red,PointSize[.025]}], 
Graphics[{Red,Thin,Line[{P1,P2}]}],
Axes->True, AxesOrigin->{0,0},AxesStyle->{Thickness[.005],Thickness[.005]},GridLines->Automatic,PlotRange->AG,AspectRatio->Automatic,ImageSize->250]],
Print[Show[
Graphics3D[{RGBColor[z[[r]]],Opacity[0.5],Ell[[r]]}],Graphics3D[{RGBColor[z[[s]]],Opacity[0.5],Ell[[s]]}],
ListPointPlot3D[PI0,PlotStyle->Table[{PointSize[.001],Opacity[.5],RGBColor[z[[j]]]},{j,Length[cen]}]],
ListPointPlot3D[{P1,P2}, PlotStyle->{Red,PointSize[Large]}],
Graphics3D[{Red,Thin,Line[{P1,P2}]}],
Axes->True,AxesStyle->{Red,Green,Blue},AxesLabel->{"R","G","B"},ImageSize->250]]
		];
			];
,{s,r+1,k}]
,{r,k}];
(* Print["DDD=",MatrixForm[DDD]," - distance matrix between clusters"]; *)
m=Sum[Length[PI[[j]]],{j,Length[PI]}];Print["m=",m];
mf=MWF[PI,z,Scov];Print["mf = ",mf//N];
(m-k)MinR[DDD] /mf
]

Minimal and the Maximal element of the upper triangular matrix

MinR[R_]:=Module[{n=Length[R[[1]]],min=Infinity},
Do[
red=Drop[R[[i]],{1,i}];
minn=Min[red];
If[minn<min,min=minn];
,{i,n-1}];
min]

MaxR[R_]:=Module[{n=Length[R[[1]]],max=-Infinity},
Do[
red=Drop[R[[i]],{1,i}];
maxn=Max[red];
If[maxn>max,max=maxn];
,{i,n-1}];
max]

IntersectTest

IntersectTest[c1_,cov1_,c2_,cov2_]:=Module[{FF,max},
FF[lambda_]:=lambda (1-lambda) (c2-c1).Inverse[lambda cov2+(1-lambda) cov1].(c2-c1);
max=First[FindMaximum[{FF[x],0<x<1},x]];
max>1
]

Weighted Mahalanobis k - means

WKMeans[Pod_,z_]:=Module[{m=Length[Pod],k=Length[z],w=Pod[[All,1]],PI0,PIpriv,tab,min,Iter=20,imin,c,c0,n=Length[z[[1]]],Scov,pom,corr,It=0,epsilon=.000000001},
c0=z;
						(* Računanje particije *)
PI=Table[{},{j,k}];Scov=Table[0,{j,k}];pom=Table[0,{j,k}];mj=Table[0,{j,k}];
Do[
tab=Table[Norm[Pod[[i,2]]-z[[j]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI];
                                                           (* Računanje centara *)
c=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,k}];
							(* LOOP *) 
While[
Chop[Norm[c0-c]]≠0&& It<Iter,c0=c;It=It+1;
				(* Racunanje matrica Scov *) 
Do[
mj[[j]]=Sum[PI[[j,s,1]],{s,Length[PI[[j]]]}];
Scov[[j]]=Sum[PI[[j,s,1]]*KroneckerProduct[(PI[[j,All,2]][[s]]-c[[j]]),(PI[[j,All,2]][[s]]-c[[j]])],{s,Length[PI[[j]]]}]/mj[[j]]
,{j,k}];
Do[
If[Min[Eigenvalues[Scov[[j]]]]<epsilon,Print["Cov.Mat. ",j,"-tog centra (",c[[j]],"): ",Scov[[j]]," je singularna"];Print["|PI[j]|: ",Length[PI[[j]]],", Eigen_min=",Min[Eigenvalues[Scov[[j]]]]];Scov[[j]]=IdentityMatrix[n] ]
,{j,k}];
PIpriv=PI;
PI=Table[{},{i,k}];
corr=Table[Log[m^2 Det[Scov[[jj]]]/mj[[jj]]^2],{jj,k}];
Do[
Do[
tab=Table[dM[PIpriv[[jj,s,2]],c[[j]],Scov[[j]]],{j,k}];
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],PIpriv[[jj,s]] ],
{s,Length[PIpriv[[jj]]]}],
{jj,k}];
PI=DeleteCases[PI,{}];
k=Length[PI];
c=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,k}];
]; 
				(* Slucaj smanjenja broja klastera u particiji *)
If[Length[c]<Length[c0],
F2[u_,v_]:=Sum[w[[i]]Min[Append[Table[Norm[Pod[[i,2]]-c[[s]]]^2,{s,Length[c]}],Norm[Pod[[i,2]]-{u,v}]^2]],{i,m}];
dom= TT[AG];
								(* Novi DIRECT centar*)
minp=DIRECT[F2,dom,DMin->10.^(-3),MaxIter->50, Eps->10^(-4)]//Timing; min=minp[[2]];Print["DIRECT: ",min];
c=Append[c,min[[1]]];
Print["Novi optimirani centar: ",c];
k=Length[c]; Print["Korigirani k: ",k];
								(* Nova particija *)
PI=Table[{},{j,k}];
Do[
tab=Table[Norm[Pod[[i,2]]-c[[j]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
];
{PI,c,Scov,It}
]

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];
(* skaliranje funkcije na hiperkocku duljine stranice 1 *)
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]]};
(* GLAVNA PETLJA *)
Monitor[
While[(dmaxrazni[[1]]>dtest)&&(iter<OptionValue[MaxIter]),
iter++;
(* KORAK 1: traženje potencijalno optimalnih hiperkvadara *)
dmaxrazni=Union[dmax];
(* svi podaci grupirani po d *)
pot0=Table[Select[Range[Length[fc]],dmax[[#]]==dmaxrazni[[i]]&],{i,Length[dmaxrazni]}];
(* podaci s najmanjim f grupirani po d *)
pot1=Table[pot0[[i,Position[fc[[pot0[[i]]]],Min[fc[[pot0[[i]]]]]][[All,1]]]],{i,Length[dmaxrazni]}];
(* provjera ostalih uvjeta *)
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];
(* KORAK 2: dijeljenje potencijalno optimalnih *)
Do[
(* smjerovi po kojima se dijeli *)
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]];];
(* slika za 2D slučaj *)
(*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]}]
(* vraćanje u početnu domenu *)
TMin[[1]]=TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]];
(* ispis slika za 2D slučaj *)
(*If[n==2,Print[ListAnimate[anim,AnimationRunning->False]]];*)
Return[TMin];
]

Input

SetDirectory[NotebookDirectory[]]; 
asd=IP3=Import["C:/Rudolf/AB-STRUKA/1_RADOVI/B1_CLUSTERS/B1_PACO/Paco-2023/DATA/IP3.txt","Table"];
IP0=ReadList["C:/Rudolf/AB-STRUKA/1_RADOVI/B1_CLUSTERS/B1_PACO/Paco-2019/DATA/IP0.txt"];
isp=0;
Print["m0 = ",m0=Length[asd]]
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["Korištene bije redom:"];
boje={1,3.,8.,6.,5,9,2.1}/10;
Table[Graphics[{Hue[boje[[j]]],Disk[]},ImageSize->20],{j,Length[boje]}]

Earthquakes of magnitude M >= 3.5

MM=3.5;priv={};
Do[
If[asd[[i,10]]≥MM ,priv=Append[priv,asd[[i]]]]
,{i,m0}]
asd=priv;
Print["m0 = ",m0=Length[asd]]
Print["First: : ",First[asd]]
Print["Last: ",Last[asd]]
Print["Magnitudes: ",{Min[asd[[All,10]]],Max[asd[[All,10]]]}]

Detection data with unknown depths

nr=0; 
Nr=Table[0,{i,1998,2014}];
Mr=Table[0,{i,3}];
Godina=Table[1997+i,{i,Length[Nr]}];
Do[
If[asd[[i,9]]==0,nr=nr+1;(* Print[{asd[[i,1]],asd[[i,10]]}]*)];
Do[
Which[asd[[i,1]]==Godina[[j]],Nr[[j]]=Nr[[j]]+1];
,{j,Length[Nr]}];
Which[asd[[i,9]]==0&&asd[[i,10]]<4.5,Mr[[1]]=Mr[[1]]+1,
asd[[i,9]]==0&&(4.5≤ asd[[i,10]]<5),Mr[[2]]=Mr[[2]]+1,
asd[[i,9]]==0&&5≤ asd[[i,10]],Mr[[3]]=Mr[[3]]+1]
,{i,m0}]
tab=Table[{Godina[[j]],Nr[[j]]},{j,Length[Nr]}];
tabx=Table[1997+i,{i,Length[Nr]}];
taby=Table[10+50i,{i,0,10}];
ListPlot[tab, PlotStyle->{PointSize[.03]},GridLines->{tabx,taby}, ImageSize->Small]
Print["Nule: ",tab=Table[{Godina[[j]],Nr[[j]]},{j,Length[Nr]}]]
Print["{asd[[i,10]]<4.5, 4.5≤asd[[i,10]]<5.5, 5.5≤ asd[[i,10]]}: ",Mr]
Print["Total o dubina: ",{nr,N[nr/m0]}]

574/602.

New database with earthquakes of known depths (without depth=0)
Evaluate this cell to consider only earthquakes of known depths
Before evaluating this cell, always first evaluate the preceding cell

asd1={};
Do[
If[asd[[i,9]]≠ 0,asd1=Append[asd1,asd[[i]]](*,Print[i," ",asd[[i]]]*)]
,{i,m0}]
isp=m0-Length[asd1];
asd=asd1;
m=Length[asd];
Print["Ispušteno je ",isp," podataka i sada raspolažemo s ",m, " podataka"]
Print["First: : ",First[asd]]
Print["Last: ",Last[asd]]
Print["Magnitudes: ",{Min[asd[[All,10]]],Max[asd[[All,10]]]}]
Print["Depths: ",{Min[asd[[All,9]]],Max[asd[[All,9]]]}]

Data

If[isp>0,Print[Style["Podaci s poznatim dubinama",16,Bold,Blue]],Print[Style["Svi podaci",16,Bold,Blue]]]
pod3=pod4=pod5=pod6=0;
Print["m = ",m=Length[asd]]
Do[
If[asd[[i,10]]<  4,pod3=pod3+1,
If[asd[[i,10]]<5,pod4=pod4+1,pod5=pod5+1]]
,{i,Length[asd]}];
Print["pod3=",pod3,"; pod4=",pod4,"; pod5=",pod5,"; podAll= ",pod3+pod4+pod5]
(* Kontrola *)
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]

Export IP3  with known depths

(* Spremanje IP podataka magnitude M≥3.5 s poznatim dubinama. *)
SetDirectory[NotebookDirectory[]]; 
Export["IP35D.txt",asd]

Data set

pod3=pod4=pod5=pod6=0;
Print["m = ",m=Length[asd]]
Do[
If[asd[[i,10]]<  4,pod3=pod3+1,
If[asd[[i,10]]<5,pod4=pod4+1,pod5=pod5+1]]
,{i,Length[asd]}];
Print["pod3=",pod3,"; pod4=",pod4,"; pod5=",pod5,"; podAll= ",pod3+pod4+pod5]
(* Kontrola *)
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]

Construction Data set POD on : (-12, 5) x  (34, 44)

POD=Table[{asd[[i,10]],{asd[[i,8]],asd[[i,7]]},i},{i,Length[asd]}];
n=2;
w=Table[1,{i,m}]; W=Total[w];
w=Table[asd[[i,10]],{i,m}]; W=Total[w];
A=Table[{asd[[i,8]],asd[[i,7]]},{i,Length[asd]}];
AG=Table[{Min[POD[[All,2,i]]],Max[POD[[All,2,i]]]},{i,n}];
cA=Mean[WeightedData[POD[[All,2]],POD[[All,1]]]];
Print["Rectangle: ",AG,"; centroid(A) = ",cA; "; |POD|=",Length[POD]]
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,m}]
slpod=Show[slIP,ListPlot[pod4,PlotStyle->{Opacity[.6], PointSize[.005]},Axes->True,AspectRatio->Automatic],ListPlot[pod5,PlotStyle->{Opacity[.8],PointSize[.007]},Axes->True,AspectRatio->Automatic],ListPlot[pod6,PlotStyle->{Opacity[1],PointSize[.0095]},Axes->True,AspectRatio->Automatic],AxesOrigin->{-12,34},GridLines->Automatic,ImageSize->350]

Data Transformation

dd=Table[AG[[i,2]]-AG[[i,1]],{i,n}];
MD=DiagonalMatrix[1/dd]; 
MD1=DiagonalMatrix[dd];
vu=Table[Min[POD[[All,2,i]]],{i,n}];
TT[x_]:=MD.(x-vu); 
TT1[x_]:=MD1.x+vu; 
BOD=Table[{w[[i]],TT[POD[[i,2]]],asd[[i]]},{i,m}];
cB=Mean[WeightedData[BOD[[All,2]],BOD[[All,1]]]];
Print["Rectangle: ",TT[AG],"; centroid(B) = ",cB]
ListPlot[BOD[[All,2]],PlotRange->TT[AG],AspectRatio->Automatic, ImageSize->Small ]
Print["Example: \nai=POD[[3,2]] = ",ai=POD[[3,2]], "\nbi=TT[ai] = ",bi=TT[ai],"\nTT1[bi] = ",TT1[bi]  ]

Algorithm 2 (Incremental algorithm)

(* Inicial centers *)
n=2;cen=centri={TT[{-10,37}]};  KK=6;
IndMArea=IndMMD={};
Scov=Table[IdentityMatrix[n],{i,Length[cen]}];
brk0=Length[cen];
PART={};

					(* Iterative process *)
Do[
Print["\nNa postojecih ",Length[centri]," centara dodaje se 1 novi"];
k=Length[centri];
(*
F2[u_,v_]:=Sum[w[[i]]Min[Append[Table[Norm[BOD[[i,2]]-centri[[s]]]^2,{s,Length[centri]}],Norm[BOD[[i,2]]-{u,v}]^2]],{i,m}];
*)
FCompile=Compile[{{u,_Real},{v,_Real}},
Evaluate[Sum[BOD[[i,1]]Min[Append[Table[Norm[BOD[[i,2]]-centri[[s]]]^2,{s,Length[centri]}],Norm[BOD[[i,2]]-{u,v}]^2]],{i,m}]  ]];	
F2[u_,v_]:=Apply[FCompile,{u,v}];	
dom=TT[AG];
								(* DIRECT *)
minp=DIRECT[F2,dom,DMin->10.^(-3),MaxIter->50, Eps->10^(-4)]//Timing; min=minp[[2]];Print["DIRECT: ",min];
cen=Append[centri,min[[1]]];
Print["Novi optimirani centar: ",Table[TT1[cen[[j]]],{j,Length[cen]}]];
(* Ovdje se vrijednost funkcije cilja ne moze racunati kao Fmin jer jos nemam Scov[[k+1]] *)
k=Length[cen];
slcen=ListPlot[Table[TT1[cen[[j]]],{j,k}],PlotStyle->{Black,PointSize[.025],Opacity[.7]}];
sl1=Show[slcen,ListPlot[POD[[All,2]],PlotStyle->PointSize[.005]],slIP,AspectRatio->Automatic];
					(* Poziv Mahalanobis k-means *)
t1=Timing[sol=WKMeans[BOD,cen]];
Print["Novi WKMeans profinjeni centar: ",Table[TT1[sol[[2,j]]],{j,k}]];
Scov=sol[[3]];
Print["WKMeans iteracija = ",sol[[4]],"; MWF = ",MWF[sol[[1]],sol[[2]],sol[[3]]],";  Fmin = ",Fmin[BOD,sol[[2]],sol[[3]]],"; mj: ",Table[Length[sol[[1,j]]],{j,Length[PI]}]];
Print["Time (DIRECT+WKM) = ",minp[[1]],"+",t1[[1]]," = ",minp[[1]]+t1[[1]]];
marea=MArea[sol[[1]],sol[[2]],sol[[3]]]//N;
Print["MArea = ",marea];
IndMArea=Append[IndMArea,marea];
mmd=MMD[sol[[1]],sol[[2]],sol[[3]],0];
Print["MMD = ",mmd];
IndMMD=Append[IndMMD,mmd];
                                                                     (* FIGURES *)
k=Length[sol[[1]]];
slc=ListPlot[Table[TT1[sol[[2,j]]],{j,k}],PlotMarkers->{"✶",16},PlotStyle->{Blue,PointSize[.02],Opacity[.7]}];
slPI=ListPlot[Table[Table[TT1[sol[[1,j,All,2]][[s]]],{s,Length[sol[[1,j]]]}],{j,k}],PlotStyle->Table[{PointSize[.004],Hue[boje[[j]]]},{j,k}]];
							(* Elipses *)
slcon=Table[0,{j,k}];
Do[
mj=Length[sol[[1,j]]];
tab=Sort[Table[dM[TT1[sol[[2,j]]],TT1[sol[[1,j,s,2]]],MD1.sol[[3,j]].MD1],{s,mj}]];
rad=tab[[Ceiling[95mj/100]]];
slcon[[j]]=ContourPlot[dM[TT1[sol[[2,j]]],{x,y},MD1.sol[[3,j]].MD1]==rad,{x,AG[[1,1]],AG[[1,2]]},{y,AG[[2,1]],AG[[2,2]]}],{j,k}];
sl=Show[Table[slcon[[j]],{j,k}]];
sl2=Print[Show[slPI,slc,sl,slIP,AspectRatio->Automatic,GridLines->Automatic,AxesOrigin->{AG[[1,1]],AG[[2,1]]},PlotRange->AG,ImageSize->250]];
(*Print[GraphicsGrid[{{sl1,sl2}},ImageSize->600]];*)			 
centri=sol[[2]];     
PART=Append[PART,sol[[1]]]; 
,{brk,KK}]
							(* Indexes *)
br0=1;
Print["MArea = ",TabIndArea=Table[{j+br0,IndMArea[[j]]},{j,Length[IndMArea]}]];
Print["MMD = ",TabIndMMD=Table[{j+br0,IndMMD[[j]]},{j,Length[IndMMD]}]];
slMArea=ListPlot[TabIndArea,Joined->True,AxesOrigin->{1,.75Min[IndMArea]},PlotLabel->"Area",GridLines->Automatic,ImageSize->200];
slMMD=ListPlot[TabIndMMD,Joined->True,AxesOrigin->{1,.75Min[IndMMD]},PlotLabel->"MMD",GridLines->Automatic,ImageSize->200];
Print[Style[GraphicsGrid[{{slMArea,slMMD}}],Magnification->1.2]];
			(* Saving *)
SetDirectory[NotebookDirectory[]]; 
Export["PART.txt",PART]

MAPart (k=4)  Klusters: 1(1), 2(4), 3(2), 4(3)

LightBlue, EdgeForm@Thin privremeno promijenjen u Blue, EdgeForm@Thin

zN={{-9.510419284804572`,36.84883949346373`},{-3.61884990880011`,35.9116908458677`},{-1.7074441172943118`,42.80969343967378`},{2.368291266102675`,36.367830454994916`}};
zN=Table[TT[zN[[j]]],{j,Length[zN]}];
k=Length[zN];
sol1=Timing[WKMeans[BOD,zN]];
Print["Time = ",sol1[[1]]]; sol=sol1[[2]];
Print["MArea = ",MArea[sol[[1]],sol[[2]],sol[[3]]]];
Print["MMD = ",MMD[sol[[1]],sol[[2]],sol[[3]]]];
Print[sol[[4]],": MWF = ",MWF[sol[[1]],sol[[2]],sol[[3]]],";  Fmin = ",Fmin[BOD,sol[[2]],sol[[3]]],"; mj: ",Table[Length[sol[[1,j]]],{j,Length[sol[[1]]]}]];
Print["Centri normirani: ",sol[[2]]]
Print["Centri pravi: ",centroidi=Table[TT1[sol[[2,j]]],{j,k}]]
PI1=Table[Union[Table[TT1[sol[[1,j,All,2]][[s]]],{s,Length[sol[[1,j]]]}]],{j,Length[sol[[1]]]}];
							(* FIGURE *)
(* Ovdje su sve točke iste veličine
slPI=ListPlot[Table[Table[TT1[sol[[1,j,All,2]][[s]]],{s,Length[sol[[1,j]]]}],{j,k}],PlotStyle->Table[{PointSize[.006],Hue[boje[[j]]]},{j,k}]]; 
*)
(* Ovdje veličine točaka ovise o težini, tj. magnitudi *)
clu=Table[{},{j,k}];
Do[
pod4=pod5=pod6={};
Do[
If[sol[[1,j,i,1]]≤4,pod4=Append[pod4,TT1[sol[[1,j,i,2]]]],
If[sol[[1,j,i,1]]≤5,pod5=Append[pod5,TT1[sol[[1,j,i,2]]]],pod6=Append[pod6,TT1[sol[[1,j,i,2]]]]]]
,{i,Length[sol[[1,j]]]}];
clu[[j]]=Show[slIP,ListPlot[pod4,PlotStyle->{Opacity[.6], PointSize[.005],Hue[boje[[j]]]},Axes->True,AspectRatio->Automatic],ListPlot[pod5,PlotStyle->{Opacity[.8],PointSize[.007],Hue[boje[[j]]]},Axes->True,AspectRatio->Automatic],ListPlot[pod6,PlotStyle->{Opacity[1],PointSize[.0095],Hue[boje[[j]]]},Axes->True,AspectRatio->Automatic],AxesOrigin->{-12,34},GridLines->Automatic,ImageSize->250]
,{j,Length[sol[[1]]]}]
slPI=Show[Table[clu[[j]],{j,k}]];
							(* 95% ellipses *)
(* Treba izabrati ellipses koje će biti filled, i.e. {2,6,7,8,9} 
    Kod crtanja elipsi treba malo proširiti AG: 
ne smije se dogoditi da elipsa ima više od jednog prekida *)
slcon=Table[0,{j,k}];
Do[
mj=Length[sol[[1,j]]];
tab=Table[dM[TT1[sol[[2,j]]],TT1[sol[[1,j,s,2]]],MD1.sol[[3,j]].MD1],{s,mj}];
rad=Quantile[tab,95/100];
slcon[[j]]=ContourPlot[dM[TT1[sol[[2,j]]],{x,y},MD1.sol[[3,j]].MD1]==rad,{x, AG[[1,1]]-1, AG[[1,2]]+1},{y,AG[[2,1]]-1,AG[[2,2]]+1}],{j,k}];
sl=Show[Table[slcon[[j]],{j,k}]]; 
(*                    (* Bojanje klastera *)
sl1a=Show[Table[slcon[j],{j,{6}}]/.Line[x_]->{LightBlue,EdgeForm@Thin,FaceForm[Opacity@.5],Polygon[x]},PlotRange->AG];
sl1b=Show[Table[slcon[j],{j,{2}}]/.Line[x_]->{LightBrown,EdgeForm@Thin,FaceForm[Opacity@.5],Polygon[x]},PlotRange->AG];
*)
slcentroidi=ListPlot[centroidi,PlotMarkers->{"✶",16},PlotStyle->{Blue,Opacity[.7]}];
sl2=Print[Show[slPI,slIP,sl,slcentroidi,AspectRatio->Automatic,PlotRange->AG,GridLines->Automatic,AxesOrigin->{AG[[1,1]],AG[[2,1]]},ImageSize->350]];

					(* Glavne Elipse *)
El=Table[{},{j,k}];rad=Table[0,{j,k}]; eig=Table[0,{j,k}];
Do[
mj=Length[sol[[1,j]]];
rad[[j]]=Sum[dM[TT1[sol[[2,j]]],TT1[sol[[1,j,s,2]]],MD1.sol[[3,j]].MD1],{s,mj}]/mj;
Do[If[dM[TT1[sol[[2,j]]],TT1[sol[[1,j,s,2]]],MD1.sol[[3,j]].MD1]≤ rad[[j]],El[[j]]=AppendTo[El[[j]],TT1[sol[[1,j,s,2]]]]  ];
,{s,mj}];
,{j,k}];
	(* Struktura optimalne particije *)
(*
pi=Table[{},{j,k}];
Do[
Do[
Do[
If[A[[i]]==PI1[[j,s]],pi[[j]]=AppendTo[pi[[j]],i]]
,{s,Length[PI1[[j]]]}];
,{i,m}]; 
,{j,k}]
*)
Print["Centroidi: ",centroidi]
Print[k,": Br.el.po-clust: ",tt=Table[Length[sol[[1,j]]],{j,k}],";   Ukupno: ",Total[tt]]
Print[k,": Br.el.po HEllipse: ",ttE=Table[Length[El[[j]]],{j,k}],";   Ukupno: ",Total[ttE]]
Print["M-radijusi Glavnih Elipsi: ",Sqrt[rad]]
Do[
Print["Granice Glavnih Elipsi (po atributima} ",j,"-tog klastera:"];
Print[Table[{s,Round[{Min[El[[j,All,s]]],Max[El[[j,All,s]]]}]},{s,n}]];
,{j,k}];
sdw=Table[Sqrt[Sum[sol[[1,j,s,1]]dM[centroidi[[j]],TT1[sol[[1,j,s,2]]],MD1.sol[[3,j]].MD1],{s,Length[sol[[1,j]]]}]/Sum[sol[[1,j,s,1]],{s,Length[sol[[1,j]]]}]],{j,k}];
sd=Table[Sqrt[Sum[dM[centroidi[[j]],TT1[sol[[1,j,s,2]]],MD1.sol[[3,j]].MD1],{s,Length[sol[[1,j]]]}]/tt[[j]]],{j,k}];
Print["Prosj.Rasipanje po clusteru {sd,sdw}: ",{sd,sdw}, ";  Suma = ",Total[sd]]
Print["Eigenvalues i Eigenvectors normirane Scov po klasterima: "]
Do[
Print["EigenValues[Scov] ",j,": ",sv=Eigenvalues[sol[[3,j]]],"; Izduzenost=",Max[sv]/Min[sv]];
,{j,k}];
Do[
Print["EigenVectors[Scov] ",j,": ",Eigenvectors[sol[[3,j]]]];
,{j,k}];

Cluster renumbering

nsol={{sol[[1,1]],sol[[1,2]],sol[[1,4]],sol[[1,3]]},{sol[[2,1]],sol[[2,2]],sol[[2,4]],sol[[2,3]]},{sol[[3,1]],sol[[3,2]],sol[[3,4]],sol[[3,3]]}};

j=4;Show[ListPlot[Table[TT1[nsol[[1,j,s,2]]],{s,Length[nsol[[1,j]]]}],PlotRange->AG],slIP]

Export MAPart4  and   MAPart4D  with known depths

(* Spremanje MAPart4 is poznatim dubinama. Particija je zadana kao: (w_i, a_i, Data_i); Nakon toga Centroids and Covariance matrices *)
SetDirectory[NotebookDirectory[]]; 
Export["MAPart4D.txt",nsol]

Earthquakes in clusters

Print["Br.eq.po klasteru:",tt,"; Total=",Total[tt]]
ind=Table[{},{j,k}];
Do[
Do[
ind[[j]]=Append[ind[[j]],sol[[1,j,j1,3]]]
,{j1,Length[sol[[1,j]]]}]
,{j,k}]
j=4;
Length[ind[[j]]]
ind[[j]]