diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index b74cf34..5e2b674 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -152,19 +152,19 @@ def apply_beast(act_rate): prob : A list of probabilities and change points extracted from BEAST results. """ - mirror_len = int(np.ceil(0.20 * len(act_rate))) - left_mirror = act_rate[:mirror_len][::-1] - right_mirror = act_rate[-mirror_len:][::-1] - act_rate_mirrored = np.concatenate([left_mirror, act_rate, right_mirror]) + #mirror_len = int(np.ceil(0.20 * len(act_rate))) + #left_mirror = act_rate[:mirror_len][::-1] + #right_mirror = act_rate[-mirror_len:][::-1] + #act_rate_mirrored = np.concatenate([left_mirror, act_rate, right_mirror]) 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], torder_minmax=[torder_min, torder_max], tseg_minlength=2, mcmc_chains=10, mcmc_thin=mcmc_th, mcmc_seed=10) - # User-driven ncp selection +# User-driven ncp selection if ncp_choice == 'median': ncp = beast_result.trend.ncp_median if np.isnan(ncp) or ncp == 0: @@ -187,17 +187,20 @@ def apply_beast(act_rate): return beast_result, np.array([]) 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] # Discard mirrored zone changepoints and correct indices - valid_mask = (cps > mirror_len) & (cps <= mirror_len + len(act_rate)) - cps = cps[valid_mask] - mirror_len + #valid_mask = (cps > mirror_len) & (cps <= mirror_len + len(act_rate)) + #cps = cps[valid_mask] - mirror_len - # Discard changepoints too close to the end (artifacts of end-mirroring). - # bins_after_cp sets the minimum number of bins required after a changepoint. - bins_after_cp = 1 - if len(cps) > 0: - cps = cps[cps <= len(act_rate) - bins_after_cp] + # Discard changepoints too close to the start or end (artifacts of mirroring). + # bins_after_cp / bins_before_cp set the minimum buffer bins required at each end. + #bins_before_cp = 2 + #bins_after_cp = 2 + #if len(cps) > 0: + # cps = cps[(cps >= bins_before_cp) & (cps <= len(act_rate) - bins_after_cp)] return beast_result, np.sort(cps)