Merge branch 'master' into ftong-patch-cython

Reviewed-on: #18
Reviewed-by: asia <asia@noreply.example.org>

LICENCE.txt

@@ -0,0 +1,13 @@
+Copyright 2025 Institute of Geophysics, Polish Academy of Sciences
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.

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
@@ -45,22 +45,19 @@ 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.io import read_mat_cat, read_mat_m, read_mat_mc, read_mat_pdf, read_csv
     from igfash.window import win_CTL, win_CNE
     import igfash.kde as kde
     from igfash.gm import compute_IMT_exceedance
@@ -84,6 +81,7 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
         m_range = [None]
     else:
         m_range = read_mat_m(m_file)
+        m_max = m_range[-1] # take m_max from the m_file
 
     if verbose:
         logger.setLevel(logging.DEBUG)
@@ -105,25 +103,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')
 
@@ -145,12 +131,18 @@ 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
         elif mc_file != None:
             logger.info("Mc estimation output file provided; selecting largest Mc from the list")
-            mc = read_mc(mc_file)
+            mc = read_mat_mc(mc_file)
             trim_to_mc = True
         else:
             logger.info("No Mc provided; using all magnitudes from the catalog")
@@ -166,9 +158,10 @@ 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
+            logger.info(f"No m_max was given. Therefore m_max is automatically set to: {m_max}")
 
         start = timer()
 
@@ -233,8 +226,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)
@@ -330,8 +321,8 @@ verbose: {verbose}")
     lambdas = [None]
     if custom_rate != None and forecast_select:
         logger.info(f"Using activity rate specified by user: {custom_rate} per {time_unit}")
-        lambdas = [custom_rate]
-        lambdas_perc = [1]
+        lambdas = np.array([custom_rate], dtype='d')
+        lambdas_perc = np.array([1], dtype='d')
 
     elif rate_select:
         logger.info(f"Activity rate modeling selected")
@@ -349,6 +340,12 @@ verbose: {verbose}")
         elif time_unit == 'years':
             multiplicator = 1 / 365
 
+        # 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)
+
         # 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]]
@@ -365,7 +362,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)
@@ -377,18 +374,29 @@ 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:
+            if m_file is None: 
+                msg = f"The selected ground motion model {model} is only valid for magnitudes up to 4.5. Please select a lower maximum magnitude."
+            else:
+                msg = f"The selected ground motion model {model} is only valid for magnitudes up to 4.5, but the provided magnitude file includes values up to {m_max}. Please adjust the magnitude range in the file accordingly."
+            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)
@@ -426,7 +434,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:
@@ -455,7 +466,7 @@ verbose: {verbose}")
                 imls = [dask.delayed(compute_IMT_exceedance)(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, p, lambdas,
                                                             forecast_len, lambdas_perc, m_range, m_pdf, m_cdf, model,
                                                             log_level=logging.DEBUG, imt=imt, IMT_min=0.0, IMT_max=2.0, rx_label=i,
-                                                            rtol=0.1, use_cython=False) for i in iter]
+                                                            rtol=0.1, use_cython=True) for i in iter]
                 iml = dask.compute(*imls)
                 iml_grid_raw.append(list(iml))
 
@@ -467,7 +478,7 @@ verbose: {verbose}")
                 for i in iter:
                     iml_i = compute_IMT_exceedance(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, p, lambdas, forecast_len, 
                                                     lambdas_perc, m_range, m_pdf, m_cdf, model, imt=imt, IMT_min = 0.0,
-                                                    IMT_max = 2.0, rx_label = i, rtol = 0.1, use_cython=False)
+                                                    IMT_max = 2.0, rx_label = i, rtol = 0.1, use_cython=True)
                     iml.append(iml_i)
                     logger.info(f"Estimated {imt} at rx {i} is {iml_i}")
                 iml_grid_raw.append(iml)
@@ -475,6 +486,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