official-apps/SeismicHazardForecasting
11
0
Fork 2
Code Issues Pull Requests Activity

Update src/seismic_hazard_forecasting.py

Browse Source 
update to activity rate calculation based on June 25 email
This commit is contained in:
ftong
2026-06-30 17:47:42 +02:00
parent 27396a7d82
commit 75a488a54e
1 changed files with 16 additions and 13 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/seismic_hazard_forecasting.py
+16 -13
View File
@@ -152,19 +152,19 @@ def apply_beast(act_rate):
prob : A list of probabilities and change points extracted from BEAST results. prob : A list of probabilities and change points extracted from BEAST results.
""" """
mirror_len = int(np.ceil(0.20 * len(act_rate))) #mirror_len = int(np.ceil(0.20 * len(act_rate)))
left_mirror = act_rate[:mirror_len][::-1] #left_mirror = act_rate[:mirror_len][::-1]
right_mirror = act_rate[-mirror_len:][::-1] #right_mirror = act_rate[-mirror_len:][::-1]
act_rate_mirrored = np.concatenate([left_mirror, act_rate, right_mirror]) #act_rate_mirrored = np.concatenate([left_mirror, act_rate, right_mirror])
mcmc_th = int(np.clip(np.ceil(len(act_rate) / 100), 2, 15)) mcmc_th = int(np.clip(np.ceil(len(act_rate) / 100), 2, 15))
beast_result = rb.beast(act_rate_mirrored, period=0, beast_result = rb.beast(act_rate, period=0,
tcp_minmax=[0, tcp_max], tcp_minmax=[0, tcp_max],
torder_minmax=[torder_min, torder_max], torder_minmax=[torder_min, torder_max],
tseg_minlength=2, mcmc_chains=10, tseg_minlength=2, mcmc_chains=10,
mcmc_thin=mcmc_th, mcmc_seed=10) mcmc_thin=mcmc_th, mcmc_seed=10)
# User-driven ncp selection # User-driven ncp selection
if ncp_choice == 'median': if ncp_choice == 'median':
ncp = beast_result.trend.ncp_median ncp = beast_result.trend.ncp_median
if np.isnan(ncp) or ncp == 0: if np.isnan(ncp) or ncp == 0:
@@ -187,17 +187,20 @@ def apply_beast(act_rate):
return beast_result, np.array([]) return beast_result, np.array([])
ncp = int(ncp) ncp = int(ncp)
# Filter NaNs first — BEAST fills unused cp slots with nan
valid_cps = beast_result.trend.cp[~np.isnan(beast_result.trend.cp)]
cps = beast_result.trend.cp[:ncp] cps = beast_result.trend.cp[:ncp]
# Discard mirrored zone changepoints and correct indices # Discard mirrored zone changepoints and correct indices
valid_mask = (cps > mirror_len) & (cps <= mirror_len + len(act_rate)) #valid_mask = (cps > mirror_len) & (cps <= mirror_len + len(act_rate))
cps = cps[valid_mask] - mirror_len #cps = cps[valid_mask] - mirror_len
# Discard changepoints too close to the end (artifacts of end-mirroring). # Discard changepoints too close to the start or end (artifacts of mirroring).
# bins_after_cp sets the minimum number of bins required after a changepoint. # bins_after_cp / bins_before_cp set the minimum buffer bins required at each end.
bins_after_cp = 1 #bins_before_cp = 2
if len(cps) > 0: #bins_after_cp = 2
cps = cps[cps <= len(act_rate) - bins_after_cp] #if len(cps) > 0:
# cps = cps[(cps >= bins_before_cp) & (cps <= len(act_rate) - bins_after_cp)]
return beast_result, np.sort(cps) return beast_result, np.sort(cps)