From eb14b7d2353c689894f00f9134840d83e5874c3d Mon Sep 17 00:00:00 2001
From: ftong <ftong@noreply.example.org>
Date: Tue, 1 Jul 2025 16:06:11 +0200
Subject: [PATCH 01/14] update package version logging, remove cpu tracking

---
 src/seismic_hazard_forecasting.py | 32 ++++++++------------------------
 1 file changed, 8 insertions(+), 24 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index e1ad7ac..c5468fd 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -45,20 +45,16 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
     """
 
     import sys
+    from importlib.metadata import version
     import logging
     from base_logger import getDefaultLogger
     from timeit import default_timer as timer
     from math import ceil, floor, isnan
     import numpy as np
-    import scipy
-    import obspy
-    import dask
     from dask.diagnostics import ProgressBar  # use Dask progress bar
     import kalepy as kale
     import utm
     from skimage.transform import resize
-    import psutil
-    import openquake.engine
     import igfash
     from igfash.io import read_mat_cat, read_mat_m, read_mat_pdf, read_csv
     from igfash.window import win_CTL, win_CNE
@@ -105,25 +101,13 @@ verbose: {verbose}")
 
     # print key package version numbers
     logger.debug(f"Python version {sys.version}")
-    logger.debug(f"Numpy version {np.__version__}")
-    logger.debug(f"Scipy version {scipy.__version__}")
-    logger.debug(f"Obspy version {obspy.__version__}")
-    logger.debug(f"Openquake version {openquake.engine.__version__}")
-    logger.debug(f"Igfash version {igfash.__version__}")
-
-    # print number of cpu cores available
-    ncpu = psutil.cpu_count(logical=False)
-    logger.debug(f"Number of cpu cores available: {ncpu}")
-    for process in psutil.process_iter():
-        with process.oneshot():
-
-            # cpu = process.cpu_percent()
-            cpu = process.cpu_percent() / ncpu
-
-            if cpu > 1:
-                logger.debug(f"{process.name()}, {cpu}")
-
-    logger.debug(f"BASELINE CPU LOAD% {psutil.cpu_percent(interval=None, percpu=True)}")
+    logger.debug(f"Numpy version {version('numpy')}")
+    logger.debug(f"Scipy version {version('scipy')}")
+    logger.debug(f"Obspy version {version('obspy')}")
+    logger.debug(f"Openquake version {version('openquake.engine')}")
+    logger.debug(f"Igfash version {version('igfash')}")
+    logger.debug(f"Rbeast version {version('rbeast')}")
+    logger.debug(f"Dask version {version('dask')}")
 
     dask.config.set(scheduler='processes')
 

From b287715f44d1cfa2776a6ecb35854bd223ade37b Mon Sep 17 00:00:00 2001
From: ftong <ftong@noreply.example.org>
Date: Tue, 1 Jul 2025 16:09:24 +0200
Subject: [PATCH 02/14] use distances in km instead of m

---
 src/seismic_hazard_forecasting.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index c5468fd..31cdd32 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -410,7 +410,10 @@ verbose: {verbose}")
             rx_lat[i], rx_lon[i] = utm.to_latlon(x_rx[i], y_rx[i], utm_zone_number,
                                                  utm_zone_letter)  # get receiver location as lat,lon
 
-        # experimental - compute ground motion only at grid points that have minimum probability density of thresh_fxy
+        # convert distances from m to km because openquake ground motion models take input distances in kilometres
+        distances = distances/1000.0
+
+        # compute ground motion only at grid points that have minimum probability density of thresh_fxy
         if exclude_low_fxy:
             indices = list(np.where(fxy.flatten() > thresh_fxy)[0])
         else:

From 921de3f4230c25ebd7a95edd8add5c9e561d5e79 Mon Sep 17 00:00:00 2001
From: ftong <ftong@noreply.example.org>
Date: Tue, 1 Jul 2025 16:24:59 +0200
Subject: [PATCH 03/14] fix - activity rate figure no longer cut off at the
 bottom

---
 src/seismic_hazard_forecasting.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 31cdd32..30d0132 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -349,7 +349,7 @@ verbose: {verbose}")
         act_rate, bin_counts, bin_edges, out, pprs, rt, idx, u_e = calc_bins(np.array(datenum_data), time_unit,
                                                                              time_win_duration, dates_calc,
                                                                              rate_forecast, rate_unc_high, rate_unc_low,
-                                                                             multiplicator, quiet=True)
+                                                                             multiplicator, quiet=True, figsize=(14,9))
 
         # Assign probabilities 
         lambdas, lambdas_perc = lambda_probs(act_rate, dates_calc, bin_edges)

From bc4d57a5ab1e115dc2c2f7d937e21ff30bc635ec Mon Sep 17 00:00:00 2001
From: ftong <ftong@noreply.example.org>
Date: Wed, 2 Jul 2025 14:27:33 +0200
Subject: [PATCH 04/14] catch null result from ground motion forecasting

---
 src/seismic_hazard_forecasting.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 30d0132..8fae701 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -462,6 +462,11 @@ verbose: {verbose}")
         end = timer()
         logger.info(f"Ground motion exceedance computation time: {round(end - start, 1)} seconds")
 
+        if np.isnan(iml_grid_raw).all():
+            msg = "No valid ground motion intensity measures were forecasted. Try a different ground motion model."
+            logger.error(msg)
+            raise Exception(msg)
+
         # create list of one empty list for each imt
         iml_grid = [[] for _ in range(len(products))]  # final ground motion grids
         iml_grid_prep = iml_grid.copy()  # temp ground motion grids

From bd1ad26693761f3a80098c4d0a4b957e4a1904ee Mon Sep 17 00:00:00 2001
From: ftong <ftong@noreply.example.org>
Date: Wed, 2 Jul 2025 15:33:42 +0200
Subject: [PATCH 05/14] change default m_max from 3 to 1 mag unit above max mag
 in catalog

---
 src/seismic_hazard_forecasting.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 8fae701..d0b0605 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -150,9 +150,9 @@ verbose: {verbose}")
             lat = np.delete(lat, indices)
             lon = np.delete(lon, indices)
 
-        # if user does not provide a m_max, set m_max to 3 magnitude units above max magnitude in catalog
+        # if user does not provide a m_max, set m_max to 1 magnitude unit above max magnitude in catalog
         if m_max == None:
-            m_max = mag.max() + 3.0
+            m_max = mag.max() + 1.0
 
         start = timer()
 

From f4d2cfc3cd31e7efb499a121f4a3ec7a89d48adb Mon Sep 17 00:00:00 2001
From: ftong <ftong@noreply.example.org>
Date: Wed, 2 Jul 2025 15:43:05 +0200
Subject: [PATCH 06/14] add check that m_max is not too large for the ground
 motion model

---
 src/seismic_hazard_forecasting.py | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index d0b0605..a7fad20 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -111,6 +111,14 @@ verbose: {verbose}")
 
     dask.config.set(scheduler='processes')
 
+    # validate m_max against the grond motion model if forecasting was selected
+    models_anthro_limited = ['Lasocki2013', 'Atkinson2015', 'ConvertitoEtAl2012Geysers'] # these models require that m_max<=4.5
+    if forecast_select:
+        if m_max > 4.5 and model in models_anthro_limited: 
+            msg = "Selected ground motion model only valid up to a maximum magnitude of 4.5. Please try again with a lower maximum magnitude."
+            logger.error(msg)
+            raise Exception(msg)
+
     # run magnitude distribution modeling if selected by user and no magnitude pdf file provided
     if m_select and m_range[0] == None and m_pdf[0] == None:
         logger.info("Magnitude distribution modeling selected")

From a94193949339cdaff40899e4d26c3b74adbbc6e1 Mon Sep 17 00:00:00 2001
From: ftong <ftong@noreply.example.org>
Date: Mon, 7 Jul 2025 13:25:05 +0200
Subject: [PATCH 07/14] remove processing of depth as it is not used

---
 src/seismic_hazard_forecasting.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index a7fad20..2fc2db7 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -225,8 +225,6 @@ verbose: {verbose}")
         y_min = y.min()
         x_max = x.max()
         y_max = y.max()
-        z_min = depth.min()
-        z_max = depth.max()
 
         grid_x_max = int(ceil(x_max / grid_dim) * grid_dim)
         grid_x_min = int(floor(x_min / grid_dim) * grid_dim)

From 84440152eb5527a3de1b00a8036a0cae71e71847 Mon Sep 17 00:00:00 2001
From: ftong <ftong@noreply.example.org>
Date: Mon, 7 Jul 2025 13:40:02 +0200
Subject: [PATCH 08/14] Update src/seismic_hazard_forecasting.py

---
 src/seismic_hazard_forecasting.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 2fc2db7..3c0bb9d 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -51,6 +51,7 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
     from timeit import default_timer as timer
     from math import ceil, floor, isnan
     import numpy as np
+    import dask
     from dask.diagnostics import ProgressBar  # use Dask progress bar
     import kalepy as kale
     import utm

From 1ea6c85ab205bbaa82c26463b5a7d7534cb46864 Mon Sep 17 00:00:00 2001
From: ftong <ftong@noreply.example.org>
Date: Mon, 7 Jul 2025 13:45:17 +0200
Subject: [PATCH 09/14] Update src/seismic_hazard_forecasting.py

---
 src/seismic_hazard_forecasting.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 3c0bb9d..3417616 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -17,7 +17,7 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
                     then the program looks for a label of 'Mw' for magnitude in the catalog
         mc: The magnitude of completeness (Mc) of the catalog
         m_max:M_max. The magnitude distribution is estimated for the range from Mc to M_max. If no value is provided,
-                    then the program sets M_max to be 3 magnitude units above the maximum magnitude value in the catalog.
+                    then the program sets M_max to be 1 magnitude units above the maximum magnitude value in the catalog.
         m_kde_method: The kernel density estimator to use.
         xy_select: If True, perform an estimation of the magnitude distribution using KDE with the chosen KDE method
         grid_dim: The grid cell size (in metres) of the final ground motion product map. A smaller cell size will
@@ -160,9 +160,10 @@ verbose: {verbose}")
             lon = np.delete(lon, indices)
 
         # if user does not provide a m_max, set m_max to 1 magnitude unit above max magnitude in catalog
+        logger.info(f"mag max of catalog is: {mag.max()}")
         if m_max == None:
             m_max = mag.max() + 1.0
-
+        logger.info(f"m_max is: {m_max}")
         start = timer()
 
         t_windowed, r_windowed = win_CNE(time, [lon, lat, mag], win_size=len(mag), win_overlap=0, min_events=1)

From 60ae1c96cdef4191068181e8835bb525e42ec48c Mon Sep 17 00:00:00 2001
From: ftong <ftong@noreply.example.org>
Date: Mon, 7 Jul 2025 14:43:42 +0200
Subject: [PATCH 10/14] raise exception for null magnitude values

---
 src/seismic_hazard_forecasting.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 3417616..15826b3 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -138,6 +138,12 @@ verbose: {verbose}")
 
         time, mag, lat, lon, depth = read_mat_cat(catalog_file, mag_label=mag_label, catalog_label='Catalog')
 
+        # check for null magnitude values
+        m_null_idx = np.where(np.isnan(mag))[0]
+        if len(m_null_idx) > 0:
+            msg = f"There are null values in the magnitude column of the catalog at indices {m_null_idx}"
+            logger.error(msg)
+            raise Exception(msg)
         if mc != None:
             logger.info("Mc value provided by user")
             trim_to_mc = True

From 70c08f47ab2f263c0b29729f58c397bb3c509682 Mon Sep 17 00:00:00 2001
From: ftong <ftong@noreply.example.org>
Date: Mon, 7 Jul 2025 14:49:35 +0200
Subject: [PATCH 11/14] add info msg about m_max

---
 src/seismic_hazard_forecasting.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 15826b3..77ef2bd 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -166,10 +166,9 @@ verbose: {verbose}")
             lon = np.delete(lon, indices)
 
         # if user does not provide a m_max, set m_max to 1 magnitude unit above max magnitude in catalog
-        logger.info(f"mag max of catalog is: {mag.max()}")
         if m_max == None:
             m_max = mag.max() + 1.0
-        logger.info(f"m_max is: {m_max}")
+            logger.info(f"No m_max was given. Therefore m_max is automatically set to: {m_max}")
         start = timer()
 
         t_windowed, r_windowed = win_CNE(time, [lon, lat, mag], win_size=len(mag), win_overlap=0, min_events=1)

From 1244655a687e6fcfa775744591e0b84750c5bd85 Mon Sep 17 00:00:00 2001
From: ftong <ftong@noreply.example.org>
Date: Mon, 7 Jul 2025 15:07:33 +0200
Subject: [PATCH 12/14] Update src/seismic_hazard_forecasting.py

---
 src/seismic_hazard_forecasting.py | 25 +++++++++++++------------
 1 file changed, 13 insertions(+), 12 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 77ef2bd..806ca6c 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -112,14 +112,6 @@ verbose: {verbose}")
 
     dask.config.set(scheduler='processes')
 
-    # validate m_max against the grond motion model if forecasting was selected
-    models_anthro_limited = ['Lasocki2013', 'Atkinson2015', 'ConvertitoEtAl2012Geysers'] # these models require that m_max<=4.5
-    if forecast_select:
-        if m_max > 4.5 and model in models_anthro_limited: 
-            msg = "Selected ground motion model only valid up to a maximum magnitude of 4.5. Please try again with a lower maximum magnitude."
-            logger.error(msg)
-            raise Exception(msg)
-
     # run magnitude distribution modeling if selected by user and no magnitude pdf file provided
     if m_select and m_range[0] == None and m_pdf[0] == None:
         logger.info("Magnitude distribution modeling selected")
@@ -169,6 +161,7 @@ verbose: {verbose}")
         if m_max == None:
             m_max = mag.max() + 1.0
             logger.info(f"No m_max was given. Therefore m_max is automatically set to: {m_max}")
+
         start = timer()
 
         t_windowed, r_windowed = win_CNE(time, [lon, lat, mag], win_size=len(mag), win_overlap=0, min_events=1)
@@ -374,18 +367,26 @@ verbose: {verbose}")
 
     if forecast_select:
         products = products_string.split()
-        logger.info(
-            f"Ground motion forecasting selected with ground motion model {model} and IMT products {products_string}")
+        logger.info(f"Ground motion forecasting selected with ground motion model {model} and IMT products {products_string}")
+
+        # validate m_max against the grond motion model
+        models_anthro_limited = ['Lasocki2013', 'Atkinson2015', 'ConvertitoEtAl2012Geysers'] # these models require that m_max<=4.5
+        if m_max > 4.5 and model in models_anthro_limited: 
+            msg = f"Selected ground motion model {model} is only valid up to a maximum magnitude of 4.5. Please try again with a lower maximum magnitude."
+            logger.error(msg)
+            raise Exception(msg)
 
         if not xy_select:
             msg = "Event location distribution modeling was not selected; cannot continue..."
             logger.error(msg)
             raise Exception(msg)
-        elif m_pdf[0] == None:
+
+        if m_pdf[0] == None:
             msg = "Magnitude distribution modeling was not selected and magnitude PDF file was not provided; cannot continue..."
             logger.error(msg)
             raise Exception(msg)
-        elif lambdas[0] == None:
+        
+        if lambdas[0] == None:
             msg = "Activity rate modeling was not selected and custom activity rate was not provided; cannot continue..."
             logger.error(msg)
             raise Exception(msg)

From 36378f4d6ca404437ad9706619fccb2c6f6b1a12 Mon Sep 17 00:00:00 2001
From: ftong <ftong@noreply.example.org>
Date: Mon, 7 Jul 2025 16:36:40 +0200
Subject: [PATCH 13/14] Update src/seismic_hazard_forecasting.py

---
 src/seismic_hazard_forecasting.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 806ca6c..e41d787 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -339,6 +339,12 @@ verbose: {verbose}")
         elif time_unit == 'years':
             multiplicator = 1 / 365
 
+        # Raise an exception when time_range from the user is too large
+        if time_range/multiplicator > 0.5*(time[-1] - time[0]):
+            msg = "Activity rate estimation time window must be less than half the catalog length. Use a shorter time window."
+            logger.error(msg)
+            raise Exception(msg)
+
         # Selects dates in datenum format and procceeds to forecast value
         start_date = datenum_data[-1] - (2 * time_win_duration / multiplicator)
         dates_calc = [date for date in datenum_data if start_date <= date <= datenum_data[-1]]

From 664ab1025b5095ce3ddccc0d614129301d7bd88a Mon Sep 17 00:00:00 2001
From: ftong <ftong@noreply.example.org>
Date: Mon, 7 Jul 2025 16:42:57 +0200
Subject: [PATCH 14/14] catch time_win_duration values that are too large

---
 src/seismic_hazard_forecasting.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index e41d787..2213c35 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -339,8 +339,8 @@ verbose: {verbose}")
         elif time_unit == 'years':
             multiplicator = 1 / 365
 
-        # Raise an exception when time_range from the user is too large
-        if time_range/multiplicator > 0.5*(time[-1] - time[0]):
+        # Raise an exception when time_win_duration from the user is too large relative to the catalog
+        if time_win_duration/multiplicator > 0.5*(time[-1] - time[0]):
             msg = "Activity rate estimation time window must be less than half the catalog length. Use a shorter time window."
             logger.error(msg)
             raise Exception(msg)