216 lines
8.5 KiB
Python
216 lines
8.5 KiB
Python
|
# -----------------
|
||
|
|
# Copyright © 2024 ACK Cyfronet AGH, Poland.
|
||
|
|
# -----------------
|
||
|
|
|
||
|
|
import numpy as np
|
||
|
|
|
||
|
|
class M_max_models:
|
||
|
|
def __init__(self, data = None, f_name = None, time_win = None, space_win = None,
|
||
|
|
Mc = None, b_method = None, num_bootstraps = None,
|
||
|
|
G = None, Mu = None,
|
||
|
|
dv = None, Mo = None, SER = None,
|
||
|
|
cl = None,
|
||
|
|
ssd = None, C = None,
|
||
|
|
):
|
||
|
|
|
||
|
|
self.data = data # Candidate data table: 2darray n x m, for n events and m clonums: x, y, z, t, mag
|
||
|
|
self.f_name = f_name # Feature's name to be calculated, check: def ComputeFeaure(self)
|
||
|
|
self.time_win = time_win # Time window whihc a feature is computed in
|
||
|
|
self.space_win = space_win # Space window ...
|
||
|
|
|
||
|
|
self.Mc = Mc # Magnitude of completeness for computing b-positive
|
||
|
|
self.b_method = b_method # list of b_methods
|
||
|
|
self.num_bootstraps = num_bootstraps # Num of bootstraps for standard error estimation of b-value
|
||
|
|
|
||
|
|
self.G = G # Shear modulus
|
||
|
|
self.Mu = Mu # Friction coefficient
|
||
|
|
self.dv = dv # Injected fluid
|
||
|
|
self.SER = SER
|
||
|
|
|
||
|
|
self.Mo = Mo # Cumulative moment magnitude
|
||
|
|
self.cl = cl # Confidence level
|
||
|
|
|
||
|
|
self.ssd = ssd # Static stress drop (Shapiro et al. 2013)
|
||
|
|
self.C = C # Geometrical constant (Shapiro et al. 2013)
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
def b_value(self, b_flag):
|
||
|
|
if b_flag == '1':
|
||
|
|
return 1, None
|
||
|
|
|
||
|
|
# maximum-likelihood estimate (MLE) of b (Deemer & Votaw 1955; Aki 1965; Kagan 2002):
|
||
|
|
elif b_flag == 'b':
|
||
|
|
X = self.data[np.where(self.data[:,-1]>self.Mc)[0],:]
|
||
|
|
if X.shape[0] > 0:
|
||
|
|
b = 1/((np.mean(X[:,-1] - self.Mc))*np.log(10))
|
||
|
|
std_error = b/np.sqrt(X.shape[0])
|
||
|
|
else:
|
||
|
|
raise ValueError("All events in the current time window have a magnitude less than 'completeness magnitude'. Use another value either for 'time window', 'minimum number of events' or 'completeness magnitude'.")
|
||
|
|
return b, std_error
|
||
|
|
|
||
|
|
# B-positive (van der Elst 2021)
|
||
|
|
elif b_flag == 'bp':
|
||
|
|
|
||
|
|
# Function to perform bootstrap estimation
|
||
|
|
def bootstrap_estimate(data, num_bootstraps):
|
||
|
|
estimates = []
|
||
|
|
for _ in range(num_bootstraps):
|
||
|
|
# Generate bootstrap sample
|
||
|
|
bootstrap_sample = np.random.choice(data, size=len(data), replace=True)
|
||
|
|
# Perform maximum likelihood estimation on bootstrap sample
|
||
|
|
diff_mat = np.diff(bootstrap_sample)
|
||
|
|
diff_mat = diff_mat[np.where(diff_mat>0)[0]]
|
||
|
|
estimate = 1/((np.mean(diff_mat - np.min(diff_mat)))*np.log(10))
|
||
|
|
estimates.append(estimate)
|
||
|
|
return np.array(estimates)
|
||
|
|
|
||
|
|
diff_mat = np.diff(self.data[:,-1])
|
||
|
|
diff_mat = diff_mat[np.where(diff_mat>0)[0]]
|
||
|
|
bp = 1/((np.mean(diff_mat - np.min(diff_mat)))*np.log(10))
|
||
|
|
bootstrap_estimates = bootstrap_estimate(diff_mat, self.num_bootstraps)
|
||
|
|
std_error = np.std(bootstrap_estimates, axis=0)
|
||
|
|
return bp, std_error
|
||
|
|
|
||
|
|
# Tapered Gutenberg_Richter (TGR) distribution (Kagan 2002)
|
||
|
|
elif b_flag == 'TGR':
|
||
|
|
from scipy.optimize import minimize
|
||
|
|
|
||
|
|
# The logarithm of the likelihood function for the TGR distribution (Kagan 2002)
|
||
|
|
def log_likelihood(params, data):
|
||
|
|
beta, Mcm = params
|
||
|
|
n = len(data)
|
||
|
|
Mt = np.min(data)
|
||
|
|
l = n*beta*np.log(Mt)+1/Mcm*(n*Mt-np.sum(data))-beta*np.sum(np.log(data))+np.sum(np.log([(beta/data[i]+1/Mcm) for i in range(len(data))]))
|
||
|
|
return -l
|
||
|
|
X = self.data[np.where(self.data[:,-1]>self.Mc)[0],:]
|
||
|
|
M = 10**(1.5*X[:,-1]+9.1)
|
||
|
|
initial_guess = [0.5, np.max(M)]
|
||
|
|
bounds = [(0.0, None), (np.max(M), None)]
|
||
|
|
|
||
|
|
# Minimize the negative likelihood function for beta and maximum moment
|
||
|
|
result = minimize(log_likelihood, initial_guess, args=(M,), bounds=bounds, method='L-BFGS-B',
|
||
|
|
options={'gtol': 1e-12, 'disp': False})
|
||
|
|
beta_opt, Mcm_opt = result.x
|
||
|
|
eta = 1/Mcm_opt
|
||
|
|
S = M/np.min(M)
|
||
|
|
dldb2 = -np.sum([1/(beta_opt-eta*S[i])**2 for i in range(len(S))])
|
||
|
|
dldbde = -np.sum([S[i]/(beta_opt-eta*S[i])**2 for i in range(len(S))])
|
||
|
|
dlde2 = -np.sum([S[i]**2/(beta_opt-eta*S[i])**2 for i in range(len(S))])
|
||
|
|
std_error_beta = 1/np.sqrt(dldb2*dlde2-dldbde**2)*np.sqrt(-dlde2)
|
||
|
|
return beta_opt*1.5, std_error_beta*1.5
|
||
|
|
|
||
|
|
def McGarr(self):
|
||
|
|
b_value, b_stderr = self.b_value(self.b_method)
|
||
|
|
B = 2/3*b_value
|
||
|
|
if B < 1:
|
||
|
|
sigma_m = ((1-B)/B)*(2*self.Mu)*(5*self.G)/3*self.dv
|
||
|
|
Mmax = (np.log10(sigma_m)-9.1)/1.5
|
||
|
|
if b_stderr:
|
||
|
|
Mmax_stderr = b_stderr/np.abs(np.log(10)*(1.5*b_value-b_value**2))
|
||
|
|
else:
|
||
|
|
Mmax_stderr = None
|
||
|
|
else:
|
||
|
|
Mmax = None
|
||
|
|
Mmax_stderr = None
|
||
|
|
|
||
|
|
return b_value, b_stderr, Mmax, Mmax_stderr
|
||
|
|
|
||
|
|
def Hallo(self):
|
||
|
|
b_value, b_stderr = self.b_value(self.b_method)
|
||
|
|
B = 2/3*b_value
|
||
|
|
if b_value < 1.5:
|
||
|
|
sigma_m = self.SER*((1-B)/B)*(2*self.Mu)*(5*self.G)/3*self.dv
|
||
|
|
Mmax = (np.log10(sigma_m)-9.1)/1.5
|
||
|
|
if b_stderr:
|
||
|
|
Mmax_stderr = self.SER*b_stderr/np.abs(np.log(10)*(1.5*b_value-b_value**2))
|
||
|
|
else:
|
||
|
|
Mmax_stderr = None
|
||
|
|
else:
|
||
|
|
Mmax = None
|
||
|
|
Mmax_stderr = None
|
||
|
|
|
||
|
|
return b_value, b_stderr, Mmax, Mmax_stderr
|
||
|
|
|
||
|
|
def Li(self):
|
||
|
|
sigma_m = self.SER*2*self.G*self.dv - self.Mo
|
||
|
|
Mmax = (np.log10(sigma_m)-9.1)/1.5
|
||
|
|
if Mmax < 0:
|
||
|
|
return None
|
||
|
|
else:
|
||
|
|
return Mmax
|
||
|
|
|
||
|
|
def van_der_Elst(self):
|
||
|
|
b_value, b_stderr = self.b_value(self.b_method)
|
||
|
|
# Seismogenic_Index
|
||
|
|
X = self.data
|
||
|
|
si = np.log10(X.shape[0]) - np.log10(self.dv) + b_value*self.Mc
|
||
|
|
if b_stderr:
|
||
|
|
si_stderr = self.Mc*b_stderr
|
||
|
|
else:
|
||
|
|
si_stderr = None
|
||
|
|
|
||
|
|
Mmax = (si + np.log10(self.dv))/b_value - np.log10(X.shape[0]*(1-self.cl**(1/X.shape[0])))/b_value
|
||
|
|
if b_stderr:
|
||
|
|
Mmax_stderr = (np.log10(X.shape[0]) + np.log10(X.shape[0]*(1-self.cl**(1/X.shape[0]))))*b_stderr
|
||
|
|
else:
|
||
|
|
Mmax_stderr = None
|
||
|
|
|
||
|
|
return b_value, b_stderr, si, si_stderr, Mmax, Mmax_stderr
|
||
|
|
|
||
|
|
def L_Shapiro(self):
|
||
|
|
from scipy.stats import chi2
|
||
|
|
|
||
|
|
X = self.data[np.isfinite(self.data[:,1]),1:4]
|
||
|
|
# Parameters
|
||
|
|
STD = 2.0 # 2 standard deviations
|
||
|
|
conf = 2 * chi2.cdf(STD, 2) - 1 # covers around 95% of population
|
||
|
|
scalee = chi2.ppf(conf, 2) # inverse chi-squared with dof=#dimensions
|
||
|
|
|
||
|
|
# Center the data
|
||
|
|
Mu = np.mean(X, axis=0)
|
||
|
|
X0 = X - Mu
|
||
|
|
|
||
|
|
# Covariance matrix
|
||
|
|
Cov = np.cov(X0, rowvar=False) * scalee
|
||
|
|
|
||
|
|
# Eigen decomposition
|
||
|
|
D, V = np.linalg.eigh(Cov)
|
||
|
|
order = np.argsort(D)[::-1]
|
||
|
|
D = D[order]
|
||
|
|
V = V[:, order]
|
||
|
|
|
||
|
|
# Compute radii
|
||
|
|
VV = V * np.sqrt(D)
|
||
|
|
R1 = np.sqrt(VV[0, 0]**2 + VV[1, 0]**2 + VV[2, 0]**2)
|
||
|
|
R2 = np.sqrt(VV[0, 1]**2 + VV[1, 1]**2 + VV[2, 1]**2)
|
||
|
|
R3 = np.sqrt(VV[0, 2]**2 + VV[1, 2]**2 + VV[2, 2]**2)
|
||
|
|
|
||
|
|
L = (1/3*(1/R1**3+1/R2**3+1/R3**3))**(-1/3)
|
||
|
|
|
||
|
|
return R1, R2, R3, L
|
||
|
|
|
||
|
|
def Shapiro(self):
|
||
|
|
R1, R2, R3, L = self.L_Shapiro()
|
||
|
|
|
||
|
|
return R1, R2, R3, L, np.log10(R3**2)+(np.log10(self.ssd)-np.log10(self.C)-9.1)/1.5
|
||
|
|
|
||
|
|
def All_models(self):
|
||
|
|
|
||
|
|
return self.McGarr()[2], self.Hallo()[2], self.Li(), self.van_der_Elst()[-2], self.Shapiro()[-1]
|
||
|
|
|
||
|
|
|
||
|
|
def ComputeModel(self):
|
||
|
|
if self.f_name == 'max_mcg':
|
||
|
|
return self.McGarr()
|
||
|
|
if self.f_name == 'max_hlo':
|
||
|
|
return self.Hallo()
|
||
|
|
if self.f_name == 'max_li':
|
||
|
|
return self.Li()
|
||
|
|
if self.f_name == 'max_vde':
|
||
|
|
return self.van_der_Elst()
|
||
|
|
if self.f_name == 'max_shp':
|
||
|
|
return self.Shapiro()
|
||
|
|
if self.f_name == 'max_all':
|
||
|
|
return self.All_models()
|