Core Modules

Core Modules#

This page documents the main PyEEPAS workflow modules.

PPE Learning#

PPE Model Parameter Learning

Learn PPE (Proximity to Past Earthquakes) model parameters using Maximum Likelihood Estimation.

The PPE model decomposes earthquake occurrence rate into: - Influence of historical earthquakes (seismic kernels) - Uniform background seismicity

Parameters learned: a (intensity), d (distance decay, km), s (background rate)

ppe_learning.ppe_learning_fast(config_file='config.json', catalog_start_year=None, learn_start_year=None, learn_end_year=None, m0=None, mT=None, fit_mode='joint', spatial_samples=50, use_fast_mode=True, ppe_ref_mag='mT')[source]#

Learn PPE model parameters using Maximum Likelihood Estimation.

Parameters:

  • config_file – Configuration JSON file path

  • catalog_start_year – Catalog start year (optional, overrides config)

  • learn_start_year – Learning period start year (optional, overrides config)

  • learn_end_year – Learning period end year (optional, overrides config)

  • m0 – Completeness magnitude (optional, overrides config)

  • mT – Target magnitude threshold (optional, overrides config)

  • fit_mode – ‘joint’ (optimize a,d,s) or ‘decoupled_gr’ (fix a=1)

  • spatial_samples – Grid resolution for integration (default: 50)

  • use_fast_mode – Use fast trapezoidal integration (default: True)

  • ppe_ref_mag – Magnitude reference: ‘mT’ or ‘m0’ (default: ‘mT’)

Returns:

Fitted parameters containing {‘a’: float, ‘d’: float, ‘s’: float, ‘ln_L’: float}

Return type:

dict

Aftershock Parameter Fitting#

Aftershock Model Parameter Estimation

Estimate aftershock model parameters (ν, κ) to distinguish independent events from aftershocks using an ETAS-like model.

Model: λ’(t,m,x,y) = ν·λ₀ + κ·Σλᵢ’

Independence weights wᵢ = ν·λ₀ / λ’ are calculated for use in EEPAS learning.

fit_aftershock_params.fit_aftershock_params_fast(config_file='config.json', ppe_ref_mag=None, target_mag=None, integration_mode='fast', spatial_samples=50)[source]#

Main function for aftershock parameter estimation.

Workflow:

  1. Load parameters a, d, s learned from Step 1 (PPE)

  2. Prepare earthquake catalogs (precursors, targets, PPE sources)

  3. Optimize v, k parameters using maximum likelihood estimation (MLE)

  4. Save results for use in Step 3 (EEPAS)

Parameters:

  • config_file – Path to configuration file

  • ppe_ref_mag – PPE reference magnitude, ‘m0’ or ‘mT’ (default: ‘mT’)

  • integration_mode – Integration mode, ‘fast’ (trapezoidal) or ‘accurate’ (dblquad)

  • spatial_samples – Grid resolution for fast mode

  • target_mag – Target event magnitude threshold, ‘m0’ or ‘mT’ (default: ‘m0’), controls the lower bound of integration and magnitude range of summed events

Returns:

Dictionary containing optimized parameters {v, k, nll, success, …}

Return type:

dict

EEPAS Learning#

EEPAS Parameter Learning with Automatic Boundary Adjustment.

Learn EEPAS model parameters with automatic boundary adjustment to avoid local optima caused by parameters hitting bounds during optimization.

Defaults: Three-stage + Multistart (3 starts) + fminsearchcon optimizer

Optimization Modes (in priority order):

  1. Custom Stages (if customStages enabled in config)

  2. Single-Stage (if –single-stage specified)

  3. Three-Stage (default, or if –three-stage specified)

Examples

Default (three-stage):

python3 eepas_learning_auto_boundary.py --config config.json

Single-stage optimization:

python3 eepas_learning_auto_boundary.py --config config.json --single-stage

Custom stages (enable in config):

# Set "customStages": {"enable": true, "stages": [...]} in config.json
python3 eepas_learning_auto_boundary.py --config config_custom.json
eepas_learning_auto_boundary.parse_result_csv(csv_path)[source]#

Parse EEPAS result CSV file.

Parameters:

csv_path (str) – Path to CSV file

Returns:

Dictionary containing 9 EEPAS parameters (am, bm, Sm, at, bt, St, ba, Sa, u) and ln_likelihood

Return type:

dict

eepas_learning_auto_boundary.eepas_with_auto_boundary(config_path, m0=None, max_adjustment_rounds=3, disable_boundary_adjustment=False, tolerance=0.01, expansion_factor=2.0, nll_convergence_threshold=0.1, single_stage=True, multi_start=True, n_starts=3, use_basinhopping=False, basinhopping_niter=20, optimizer='fminsearchcon', ppe_ref_mag=None, use_fast_mode=False, magnitude_samples=50, lead_time_days=None)[source]#

EEPAS parameter learning with automatic boundary adjustment.

Workflow:

  1. Execute EEPAS optimization with current boundaries (three-stage or single-stage)

  2. Check if any result parameters hit boundaries

  3. If hit, relax boundaries and save new configuration

  4. Check if NLL has converged

  5. Repeat until convergence or maximum rounds reached

Parameters:

  • config_path (str) – Configuration file path

  • m0 (float) – Completeness magnitude (if not specified, uses value from config file)

  • max_adjustment_rounds (int) – Maximum boundary adjustment rounds (default: 3), usually converges in 2-3 rounds

  • tolerance (float) – Boundary hit detection tolerance (default: 0.01 = 1%), parameter distance to boundary < 1% is considered hitting

  • expansion_factor (float) – Boundary expansion factor (default: 2.0), lower bound: new = old / 2.0, upper bound: new = old * 2.0

  • nll_convergence_threshold (float) – NLL convergence threshold (default: 0.1), two consecutive rounds with improvement < 0.1 → automatic stop

  • single_stage (bool) – Whether to use single-stage full parameter optimization (default: True, uses single-stage)

  • multi_start (bool) – Whether to use multi-start optimization (default: True)

  • n_starts (int) – Number of multi-start points (default: 3)

  • use_basinhopping (bool) – Whether to use basinhopping global optimization (default: False)

  • basinhopping_niter (int) – Number of basinhopping iterations (default: 20)

  • optimizer (str) – Optimizer to use: ‘fminsearchcon’, ‘SLSQP’, ‘L-BFGS-B’ (default: ‘fminsearchcon’)

  • ppe_ref_mag (str) – PPE reference magnitude: ‘mT’ or ‘m0’ (default: None, uses config)

  • use_fast_mode (bool) – Whether to use fast magnitude integration (default: False)

  • magnitude_samples (int) – Number of magnitude samples for fast mode (default: 50)

  • lead_time_days (float) – Fixed lead time L in days for FLEEPAS (default: None) - None: Read from config’s timeComp.lead_time_days, or use all historical precursors if not set - float: Only use earthquakes within [t-L, t-delay] as precursors (FLEEPAS mode)

  • disable_boundary_adjustment (bool)

Returns:

Final parameter dictionary containing optimized EEPAS parameters

Return type:

dict

eepas_learning_auto_boundary.main()[source]#

PPE Forecast#

PPE Earthquake Forecasting

Generate PPE model seismicity rate forecasts for testing region grid cells.

Uses learned PPE parameters (a, d, s) to calculate background earthquake occurrence rate λ₀(t,x,y) based on proximity to past earthquakes.

ppe_make_forecast.compute_grid_integral_fast(x_bounds, y_bounds, event_x, event_y, event_m, a, d, s, ppe_ref_mag, coeff, spatial_samples=20)[source]#

Fast computation of PPE integral using trapezoidal integration (replaces dblquad).

Parameters:

  • x_bounds – Grid x range [x_min, x_max]

  • y_bounds – Grid y range [y_min, y_max]

  • event_x – Event coordinates and magnitudes

  • event_y – Event coordinates and magnitudes

  • event_m – Event coordinates and magnitudes

  • a – PPE parameters

  • d – PPE parameters

  • s – PPE parameters

  • ppe_ref_mag – PPE parameters

  • coeff – Temporal and magnitude integration coefficient

  • spatial_samples – Number of sampling points per dimension (default 20, adjustable for precision)

Returns:

Integral value

Return type:

float

Note

Refactored to use fast_kernel_sum_2d and trapezoidal_2d from unified module

ppe_make_forecast.ppe_make_forecast(config_file='config.json', catalog_start_year=None, forecast_start_year=None, forecast_end_year=None, celln=None, m0=None, mT=None, delay=None, use_rolling_update=None, forecast_period_days=None, ppe_ref_mag='mT', use_fast_mode=True, spatial_samples=50)[source]#

Main function for PPE model earthquake forecasting.

Workflow:

  1. Load PPE parameters (a, d, s) from the learning phase

  2. Load historical earthquake catalog

  3. For each forecast time point:

    • Calculate seismic kernel contribution for each grid cell

    • Add background seismicity rate

    • Generate seismicity rate map for that time point

  4. Save results in MATLAB format for subsequent analysis

Parameters:

  • config_file – Configuration file path

  • catalog_start_year – Catalog start year

  • forecast_start_year – Forecast start year

  • forecast_end_year – Forecast end year

  • celln – Number of spatial grid cells

  • m0 – Completeness magnitude

  • mT – Target magnitude (for calculating seismic kernel)

  • delay – Delay in days (how long after an earthquake before it is counted)

  • use_rolling_update – Whether to use rolling update True: Each time point uses the latest data (default) False: All time points use the same historical data

  • forecast_period_days – Forecast period length (days, default 91.31 ≈ 3 months)

  • use_fast_mode – Whether to use fast mode (default True) True: Use Numba JIT accelerated grid integration (60-70x faster) False: Use accurate dblquad numerical integration (slow but precise)

  • spatial_samples – Number of sampling points for fast mode (default 20) Increasing improves precision but reduces speed

Returns:

Forecast results saved as .mat file

EEPAS Forecast#

EEPAS Earthquake Forecasting

Generate EEPAS model forecasts combining PPE background rate with precursory signals from historical earthquakes.

Final forecast: λ = μ·λ₀ + Σ wᵢ·ηᵢ·λᵢ

Each historical earthquake contributes a precursory kernel with temporal (lognormal), magnitude (Gaussian), and spatial (exponential) distributions.

eepas_make_forecast.compute_spatial_contribution_fast(X1, X2, Y1, Y2, xee, yee, sigma_spatial, qq_other, W_use)[source]#

Compute spatial contributions for all grid cells using Numba acceleration.

Parameters:

  • X1 – Grid boundary arrays (celln,)

  • X2 – Grid boundary arrays (celln,)

  • Y1 – Grid boundary arrays (celln,)

  • Y2 – Grid boundary arrays (celln,)

  • xee – Event coordinate arrays (n_events,)

  • yee – Event coordinate arrays (n_events,)

  • sigma_spatial – Spatial standard deviation array (n_events,)

  • qq_other – Product of other contribution terms (scale×magnitude×time×weight) (n_events,)

  • W_use – Weight array (n_events,)

Returns:

Expected number of events in each grid cell (celln,)

Return type:

numpy.ndarray

eepas_make_forecast.apply_time_compensation(ppe_rate, eepas_rate_with_eta, u, time_comp_config, forecast_start_time, catalog_end_time, N_tre_M, eepas_params, forecast_period_days=91.31, catalog_start_time=None)[source]#

Apply time compensation to EEPAS forecast.

Implements LEEPAS (Lead-time compensated EEPAS) method from Rhoades et al. (2019).

Parameters:

  • ppe_rate – PPE background rate λ₀ (N_rows, N_cells)

  • eepas_rate_with_eta – EEPAS time-varying rate Σ η(mᵢ)λᵢ (N_rows, N_cells) where η(mᵢ) already includes (1-μ) for magnitude completeness

  • u – Failure-to-predict rate (μ parameter)

  • time_comp_config – Time compensation configuration dict

  • forecast_start_time – Forecast start time (decimal years)

  • catalog_end_time – Catalog end time (decimal years)

  • N_tre_M – Number of time windows

  • eepas_params – EEPAS parameters dict

  • forecast_period_days – Length of each forecast window (days, default 91.31)

  • catalog_start_time – Catalog start time (decimal years). If provided, Lead Time (L) is automatically calculated as the full catalog history length (forecast_start_time - catalog_start_time). This is the RECOMMENDED way to use LEEPAS, as EEPAS precursors can occur years to decades before target earthquakes. If not provided, falls back to config’s lead_time_days.

Returns:

Compensated forecast rate (N_rows, N_cells)

Return type:

compensated_rate

Formulas (from Rhoades et al. 2019, Equations 12-14):
Mode A (background compensated):

λ_A = [μ + (1-μ)(1-p)] λ₀ + Σ η(mᵢ)λᵢ

Mode B (time-varying compensated):

λ_B = μ λ₀ + (1/p) Σ η(mᵢ)λᵢ

Mode C (optimal mixture):

λ_C = ω λ_A + (1-ω) λ_B

where:

p(T,L,m) = completeness function (0-1) η(mᵢ) = magnitude scaling function, includes (1-μ) for magnitude completeness

Note

The (1-μ) in η is for ensuring Gutenberg-Richter distribution (magnitude completeness), NOT the same as the mixing weight in standard EEPAS formula.

Per-period time-lag: Each forecast period may have a different time-lag T in non-rolling update mode. This function correctly handles varying T across periods.

eepas_make_forecast.eepas_make_forecast(config_file='config.json', catalog_start_year=None, forecast_start_year=None, forecast_end_year=None, celln=None, m0=None, weight_flag=None, delay=None, use_causal_ew=None, use_rolling_update=None, forecast_period_days=None, use_fast_mode=True, magnitude_samples=50, ppe_ref_mag='mT', lead_time_days=None, override_at=None, override_Sa=None, output_suffix=None, _hybrid_submodel=False)[source]#

Main function for EEPAS model earthquake forecasting.

Workflow:

  1. Load all parameters from three learning stages:

    • PPE parameters (a,d,s): Step 1

    • Aftershock parameters (ν,κ): Step 2

    • EEPAS parameters (am,bm,Sm,at,bt,St,ba,Sa,u): Step 3

  2. Calculate weights wᵢ for each historical earthquake (aftershock down-weighting)

  3. For each forecast time window:

    • Calculate precursory signal contribution for each grid cell

    • Add PPE background rate

    • Generate complete seismicity rate map

  4. Save in MATLAB format

Parameters:

  • config_file – Path to configuration file

  • catalog_start_year – Starting year of earthquake catalog

  • forecast_start_year – Starting year of forecast period

  • forecast_end_year – Ending year of forecast period

  • celln – Number of spatial grid cells (uses CELLE count if not specified)

  • m0 – Completeness magnitude

  • weight_flag – Weight calculation mode 0 = Uniform weights (wᵢ=1, no aftershock consideration) 1 = Compute weights (aftershock down-weighting using ETAS model)

  • delay – Delay in days (how long after an earthquake before counting precursory signal)

  • use_causal_ew – EW weight calculation mode 1 (dynamic): Each period uses all events up to (t1-delay) to compute EW 0 (fixed): Each period uses all events up to (t2) to compute EW (fixed within period) Both modes are causal, source events restricted to before t1-delay

  • use_rolling_update – Whether to use rolling update for source event selection True (default): Each time window dynamically uses all events up to (t1-delay) False: All time windows use the same snapshot of events (up to T3dec-delay) Note: This is independent of use_causal_ew which controls EW weight calculation

  • forecast_period_days – Length of each forecast window (days, default 91.31≈3 months)

  • use_fast_mode – Whether to use fast mode (default True) True: Use Numba JIT-accelerated midpoint rectangle method (fast) False: Use exact quad_vec integration (slow but precise)

  • magnitude_samples – Number of samples for fast mode (default 20, increase for higher accuracy)

  • ppe_ref_mag – PPE reference magnitude selection (‘m0’ or ‘mT’, default uses mT)

  • lead_time_days – Fixed lead time L in days for FLEEPAS (default None) - None: Use all historical earthquakes as precursors (original EEPAS) - float: Only use earthquakes within [t-L, t-delay] as precursors (FLEEPAS)

  • override_at – Override fitted at (a_T) parameter for Hybrid EEPAS (default None) - None: Use fitted value from Fitted_par_EEPAS_*.csv - float: Override with specified value

  • override_Sa – Override fitted Sa (σ_A) parameter for Hybrid EEPAS (default None) - None: Use fitted value from Fitted_par_EEPAS_*.csv - float: Override with specified value

  • output_suffix – Suffix for output filename for Hybrid EEPAS (default None) - None: Use default output filename - str: Append suffix to output filename (e.g., ‘_model1’)

  • _hybrid_submodel – Internal flag to prevent recursive hybrid calls (default False) - False: Normal execution, will check config for hybrid settings - True: Running as submodel, skip hybrid processing

Hybrid EEPAS Config (in modelParams.hybrid):
{

“enable”: true, “delta_aT_values”: [-0.5, 0.0, 0.5], // Offsets along trade-off line “weights”: “equal” // Or array like [0.33, 0.34, 0.33]

}

Returns:

EEPAS forecast rate matrix with shape (N_rows, N_cols+1) where:

  • N_rows = num_time_windows × num_magnitude_bins (e.g., 40 periods × 25 mag bins = 1000)

  • N_cols = num_spatial_cells (e.g., 177 grid cells for Italy)

  • Column 0: Time window index (1, 2, 3, …)

  • Columns 1 to N_cols: Earthquake rate for each spatial cell

Each rate value represents expected number of earthquakes in that spatiotemporal-magnitude bin during the forecast period.

Return type:

np.ndarray

See Also#

Contents