Update src/seismic_hazard_forecasting.py

revert Update src/seismic_hazard_forecasting.py

gsim/lasocki_2013.py

@@ -0,0 +1,142 @@
+# -*- coding: utf-8 -*-
+# vim: tabstop=4 shiftwidth=4 softtabstop=4
+#
+# Copyright (C) 2013-2023 GEM Foundation
+#
+# OpenQuake is free software: you can redistribute it and/or modify it
+# under the terms of the GNU Affero General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# OpenQuake is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Affero General Public License for more details.
+#
+# You should have received a copy of the GNU Affero General Public License
+# along with OpenQuake. If not, see <http://www.gnu.org/licenses/>.
+
+"""
+Module exports :class:`Lasocki2013`.
+"""
+import numpy as np
+from scipy.constants import g
+
+from openquake.hazardlib.gsim.base import GMPE, CoeffsTable
+from openquake.hazardlib import const
+from openquake.hazardlib.imt import PGA, PGV, SA
+
+def get_magnitude_energy(mag):
+    """
+    Converts magnitude to energy-based magnitude term.
+    """
+    return 1.15 + 1.96 * mag
+
+def get_distance_term(coeffs, repi):
+    """
+    Computes the distance term using the given GMPE coefficients.
+    """
+    R_h = np.sqrt(repi ** 2 + coeffs["c7"] ** 2)
+    return np.log10(R_h)
+
+def get_standard_deviation(coeffs, coeffs_cov, magE, repi):
+    """
+    Computes the standard deviation term.
+    """
+    Cb = np.array(list(coeffs_cov)).reshape(3,3)
+    R_h = np.sqrt(repi ** 2 + coeffs["c7"] ** 2)
+    X0 = np.array([1, magE[0], np.log10(R_h[0]**2)])
+    variance_term = np.sqrt(X0 @ Cb @ X0 + coeffs["sigma"]**2)
+    return variance_term
+
+class Lasocki2013(GMPE):
+    """
+    Implement equation developed by Lasocki in "REPORT ON THE ATTENUATION 
+    RELATIONS OF PEAK GROUND ACCELERATION AND SPECTRAL ORDINATES OF GROUND 
+    MOTION FOR MINING-INDUCED SEISMIC EVENTS IN THE REGION OF THE ŻELAZNY 
+    MOST REPOSITORY", 2009.
+    Equation coefficients provided for the random horizontal component
+    """
+    #: Supported tectonic region type is induced, given
+    #: that the equations have been derived for the LGCD mining area
+    DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.INDUCED
+
+    #: Supported intensity measure types are spectral acceleration,
+    #: and peak ground acceleration
+    DEFINED_FOR_INTENSITY_MEASURE_TYPES = {PGA, SA}
+
+    #: Supported intensity measure component is random horizontal
+    #: :attr:`~openquake.hazardlib.const.IMC.RANDOM_HORIZONTAL`,
+    DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.RANDOM_HORIZONTAL
+    # DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.VERTICAL
+
+    #: Supported standard deviation type is total
+    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {const.StdDev.TOTAL}
+
+    #: site params are not required
+    REQUIRES_SITES_PARAMETERS = set()
+
+    #: Required rupture parameter is magnitude
+    REQUIRES_RUPTURE_PARAMETERS = {'mag'}
+
+    #: Required distance measure is epicentral distance
+    #: see paragraph 'Predictor Variables', page 6.
+    REQUIRES_DISTANCES = {'repi'}
+
+    def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
+        """
+        Computes mean ground motion values and standard deviations for Lasocki (2013) GMPE.
+        
+        Parameters:
+            repi (float): Epicentral distance (m)
+            mag (float): Earthquake magnitude
+            imts (list of str): List of intensity measure types (e.g., 'PHA', 'PVA')
+            mean (np.array): Array to store computed mean values
+            sig (np.array): Array to store computed standard deviations
+        """
+        # Loop through each IMT and compute values
+        for i, imt in enumerate(imts):            
+            C = self.COEFFS[imt]
+            C_cov = self.COEFFS_COV[imt]
+            mag = ctx.mag
+            repi = ctx.repi*1000.0 # Convert distance from km to m
+            # Compute magnitude energy term
+            magE = get_magnitude_energy(mag)
+
+            # Compute GMPE terms
+            mag_term = C['c1'] + C['c2'] * magE
+            dist_term = C['c5'] * get_distance_term(C, repi)
+            
+            # Compute mean ground motion
+            imean = mag_term + dist_term
+            mean_value = np.log((10 ** imean) / g)  # Convert to natural log scale and divide by g
+            mean[i] = mean_value
+
+            # Compute standard deviation
+            sigma = get_standard_deviation(C, C_cov, magE, repi)
+            sig[i] = sigma
+
+    #: coefficient table provided by report
+    COEFFS = CoeffsTable(sa_damping=5, table="""\
+    IMT  c1     c2      c5      c7    sigma     N        
+    pga  1.25   0.31    -1.34   558   0.196     1196
+    0.6  -3.86  0.55    -0.65   183   0.242     1194
+    1.0  -3.94  0.62    -0.66   308   0.262     1197
+    2.0  -1.14  0.47    -1.02   741   0.235     1195
+    5.0  0.99   0.31    -1.15   690   0.234     1206
+    10.0 3.06   0.27    -1.65   906   0.203     1192
+    20.0 2.62   0.27    -1.60   435   0.196     1191
+    50.0 2.09   0.27    -1.48   375   0.204     1191
+    """)
+    
+    COEFFS_COV = CoeffsTable(sa_damping=5, table="""\
+    IMT  Cb00      Cb01      Cb02      Cb10      Cb11      Cb12      Cb20      Cb21      Cb22        
+    pga  0.005586  -0.000376 -0.000752 -0.000376 0.000103  -0.000111 -0.000752 -0.000111 0.000440
+    0.6  0.007509  -0.000662 -0.000688 -0.000662 0.000161  -0.000154 -0.000688 -0.000154 0.000516
+    1.0  0.009119  -0.00075  -0.000948 -0.00075  0.000189  -0.000187 -0.000948 -0.000187 0.000657                 
+    2.0  0.008563  -0.000514 -0.001282 -0.000514 0.000147  -0.000164 -0.001282 -0.000164 0.000696         
+    5.0  0.008283  -0.000516 -0.001202 -0.000516 0.000145  -0.000161 -0.001202 -0.000161 0.000668
+   10.0  0.006904  -0.00036  -0.001145 -0.00036  0.000108  -0.000126 -0.001145 -0.000126 0.000578
+   20.0  0.005389  -0.000396 -0.000658 -0.000396 0.000104  -0.000107 -0.000658 -0.000107 0.000408
+   50.0  0.005874  -0.000449 -0.000678 -0.000449 0.000114  -0.000114 -0.000678 -0.000114 0.000428
+    """)

src/seismic_hazard_forecasting.py

@@ -1,5 +1,205 @@
 # -*- coding: utf-8 -*-
+import numpy as np
+from scipy.stats import t, norm
+from scipy.optimize import root_scalar
+from timeit import default_timer as timer
+import logging
+logger = logging.getLogger(__name__)
+import mpmath
+from openquake.hazardlib.geo import Point #This class represents a geographical point in terms of longitude, latitude, and depth (with respect to the Earth surface).
+from openquake.hazardlib.geo.surface.planar import PlanarSurface
+from openquake.hazardlib.source.characteristic import CharacteristicFaultSource
+from openquake.hazardlib.mfd import ArbitraryMFD
+from openquake.hazardlib.tom import PoissonTOM
+from openquake.hazardlib.scalerel import WC1994 #Wells and Coppersmith magnitude – rupture area relationships
+from openquake.hazardlib.site import Site, SiteCollection
+from openquake.hazardlib.contexts import ContextMaker
+from openquake.hazardlib.valid import gsim
+from openquake.hazardlib.imt import PGA
+import sys
 
+def my_trivial_task(x):
+    """
+    A function that does a trivial task
+    """
+    return x * x
+
+def compute_IMT_exceedance(rx_lat, rx_lon, r, fr, p, lambdas, D, percentages_D, magnitudes, magnitude_pdf, magnitude_cdf, model, imt='PGA', IMT_min=0.01, IMT_max=2.0, rx_label=None, rtol=0.1, use_cython=False, **kwargs):
+    
+    n_events = len(r)
+        
+    try:
+        gmpes = [gsim(model)]
+    except:
+        msg = f"{model} was not found in the openquake gsim directory" 
+        logger.error(msg)
+        raise Exception(msg) 
+    
+    if model == 'Lasocki2013': #this model requires the number of earthquake records
+    
+        if imt=='PGA': #extract number of records for PGA
+            num_ground_motion_records = gmpes[0].COEFFS.non_sa_coeffs[PGA()]['N']
+        else: #extract number of records for SA()
+            freq = float(imt[imt.find('(')+1:imt.find(')')]) # get the desired frequency of SA
+            first_index = np.where(gmpes[0].COEFFS.get_coeffs('N')[0]==freq)[0][0]
+            num_ground_motion_records = gmpes[0].COEFFS.get_coeffs('N')[1][first_index][0]
+    else:
+        num_ground_motion_records = 0
+                               
+    #placeholder values that do not have any effect
+    Mag = 5.0 #placeholder mag, must be valid for that context; will be overwritten in loop
+    rupture_aratio = 1.5
+    Strike = 0
+    Dip = 90
+    Rake = 0
+    
+    Hypocenter = Point(rx_lon, rx_lat, 0.0) #does not matter in our case; just set eq location to be same as receiver
+    #according to the magnitude and MSR calculate planar surface
+    planar_surface = PlanarSurface.from_hypocenter(
+            hypoc=Hypocenter,
+            msr=WC1994(),
+            mag=Mag,
+            aratio=rupture_aratio,
+            strike=Strike,
+            dip=Dip,
+            rake=Rake,
+            )
+
+    # site for which we compute (receiver location)    
+    site_collection = SiteCollection([Site(location=Point(rx_lon, rx_lat, 0))])
+
+    imtls = {s: [0] for s in [imt]} #required for context maker, M = 2 IMTs
+
+    context_maker = ContextMaker('Induced', gmpes, {'imtls': imtls, 'mags': [Mag]}) #necessary contexts builder
+
+    src = CharacteristicFaultSource(source_id = 1,
+                                    name = 'rup',
+                                    tectonic_region_type = 'Induced',
+                                    mfd = ArbitraryMFD([Mag], [0.01]), #this does not have any effect
+                                    temporal_occurrence_model = PoissonTOM(50.), #this is also not really used
+                                    surface = planar_surface,
+                                    rake = Rake)
+
+    ctx = context_maker.from_srcs([src], site_collection)[0] #returns one context from the source for one rupture   
+
+    if use_cython:
+
+        from cython_exceedance import exceedance_core
+    
+        def exceedance_root_function(a):
+            return exceedance_core(a, r, fr, lambdas, D, percentages_D, magnitudes,
+                                   magnitude_pdf, magnitude_cdf, context_maker, ctx,
+                                   model, num_ground_motion_records) - p
+
+    else:
+    
+        def exceedance_root_function(a):
+            exceedance_prob_sum = 0
+            log_a = np.log(a) # Precompute log(a)
+        
+            for j in range(len(lambdas)): #loop through all lambdas
+                lambda_j = lambdas[j]
+                D_j_val = percentages_D[j] * D # Use a different name to avoid shadowing D
+                lambda_D_j = lambda_j * D_j_val
+                denom_j = (1 - np.exp(-lambda_D_j))
+                if denom_j == 0: # Avoid division by zero if lambda_D_j is very small or zero
+                   continue
+        
+                for i in range(n_events): #loop through all events
+                    ri = r[i] # Epicentral distance
+                    fr_i = fr[i] # Location probability f(r)
+                    ctx.repi = ri
+        
+                    # Precompute terms only dependent on lambda_j, D_j, m
+                    lambda_D_j_f_m = lambda_D_j * magnitude_pdf
+                    exp_term_m = np.exp(-lambda_D_j * (1 - magnitude_cdf))
+                    f_conditional_base_m = (lambda_D_j_f_m * exp_term_m) / denom_j
+        
+                    for k in range(len(magnitudes)):  #loop through all values of magnitude pdf and cdf
+                        m = magnitudes[k]
+                        ctx.mag = m # update context magnitude 
+        
+                        # Calculate f_conditional (simpler now)
+                        f_conditional = f_conditional_base_m[k]
+        
+                        mean, sig, _, _ = context_maker.get_mean_stds(ctx)
+                        log_gm_predicted = mean[0][0][0]
+                        variance_term = sig[0][0][0]
+                        residual = log_a - log_gm_predicted # Use precomputed log_a
+        
+                        if residual <= 0:
+                             exceedance_probability = 1.0
+                        else:
+                             # Avoid division by zero or very small numbers if variance_term is ~0
+                             if variance_term < 1e-15: # Adjust threshold as needed
+                                 exceedance_probability = 0.0
+                             else:
+                                t_value = residual / variance_term
+                                
+                                if model == 'Lasocki2013':
+                                    exceedance_probability = t.sf(t_value, num_ground_motion_records - 3) # student t distribution, degrees of freedom: n-3; sf = 1 - cdf
+                                else:
+                                    exceedance_probability = norm.sf(t_value) # equivalent to 1.0 - norm.cdf(t_value)
+        
+                        location_exceedance_prob = exceedance_probability * f_conditional * fr_i
+                        exceedance_prob_sum += location_exceedance_prob
+        
+            return exceedance_prob_sum - p
+    
+
+    
+    # Check function values at different test points
+    IMT_mid = (IMT_max-IMT_min)/2
+
+
+
+    lower_bound_value = exceedance_root_function(IMT_min)
+
+
+    mid_point_value = exceedance_root_function(IMT_mid)
+
+
+
+    upper_bound_value = exceedance_root_function(IMT_max)
+    
+
+
+
+    logger.info(f"Receiver: {str(rx_label)}")
+    logger.info(f"Function value at {imt} = {str(IMT_min)} : {lower_bound_value}")
+    logger.info(f"Function value at {imt} = {str(IMT_mid)} : {mid_point_value}")
+    logger.info(f"Function value at {imt} = {str(IMT_max)} : {upper_bound_value}")
+
+
+
+    if np.sign(lower_bound_value) == np.sign(upper_bound_value):
+        msg = "Function values at the interval endpoints must differ in sign for fsolve to work. Expand the interval or use a different model."
+        logger.error(msg)
+        gm_est = np.nan
+        return gm_est
+        # raise ValueError(msg)
+    
+
+
+    # Find root of function
+    start = timer()       
+        
+    try: 
+        method='brenth'
+        logger.debug("Now trying Scipy " + method)
+        output = root_scalar(exceedance_root_function, bracket=[IMT_min, IMT_max], rtol=rtol, method=method)
+        #output = root_scalar(exceedance_root_function, method=method)
+        gm_est = output.root
+    except Exception as error:
+        logger.error(f"An exception occurred: {error}")
+        logger.info("Set ground motion value to nan")
+        gm_est = np.nan
+        
+    end = timer()
+    logger.info(f"Ground motion estimation computation time: {round(end - start,1)} seconds")
+    logger.info(f"Estimated {imt}: {gm_est}")  
+    
+    return gm_est
 
 def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max,
          m_kde_method, xy_select, grid_dim, xy_win_method, rate_select, time_win_duration,
@@ -60,7 +260,7 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
     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
+    #from igfash.gm import compute_IMT_exceedance
     from igfash.compute import get_cdf, hellinger_dist, cols_to_rows
     from igfash.rate import lambda_probs, calc_bins, bootstrap_forecast_rolling
     from igfash.mc import estimate_mc
@@ -110,6 +310,7 @@ verbose: {verbose}")
     logger.debug(f"Igfash version {version('igfash')}")
     logger.debug(f"Rbeast version {version('rbeast')}")
     logger.debug(f"Dask version {version('dask')}")
+    logger.debug(f"Logging version {logging.__version__}")
     
     dask.config.set(scheduler='processes')
 
@@ -366,6 +567,8 @@ verbose: {verbose}")
 
         # Assign probabilities 
         lambdas, lambdas_perc = lambda_probs(act_rate, dates_calc, bin_edges)
+        lambdas = np.array(lambdas, dtype='d')
+        lambdas_perc = np.array(lambdas_perc, dtype='d')
         
         # print("Forecasted activity rates: ", lambdas, "events per", time_unit[:-1])
         logger.info(f"Forecasted activity rates: {lambdas} events per {time_unit} with percentages {lambdas_perc}")
@@ -452,7 +655,7 @@ verbose: {verbose}")
 
         start = timer()
 
-        use_pp = False
+        use_pp = True
 
         if use_pp:         # use dask parallel computing
             pbar = ProgressBar()
@@ -463,13 +666,16 @@ verbose: {verbose}")
             for imt in products:
                 logger.info(f"Estimating {imt}")
 
-                imls = [dask.delayed(compute_IMT_exceedance)(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, p, lambdas,
+                imls = [dask.delayed(compute_IMT_exceedance)(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, p, lambdas, forecast_len,
-                                                            forecast_len, lambdas_perc, m_range, m_pdf, m_cdf, model,
+                                                             lambdas_perc, m_range, m_pdf, m_cdf, model, imt=imt, IMT_min=0.0,
-                                                            log_level=logging.DEBUG, imt=imt, IMT_min=0.0, IMT_max=2.0, rx_label=i,
+                                                             IMT_max=2.0, rx_label=i, rtol=0.1, use_cython=False) for i in iter[:2]]
-                                                            rtol=0.1, use_cython=True) for i in iter]
+                #imls = [dask.delayed(my_trivial_task)(i) for i in iter]
                 iml = dask.compute(*imls)
                 iml_grid_raw.append(list(iml))
 
+
+
+
         else:
             iml_grid_raw = []
             iter = indices
@@ -478,7 +684,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=True)
+                                                    IMT_max = 2.0, rx_label = i, rtol = 0.1, use_cython=False)
                     iml.append(iml_i)
                     logger.info(f"Estimated {imt} at rx {i} is {iml_i}")
                 iml_grid_raw.append(iml)
@@ -486,11 +692,15 @@ verbose: {verbose}")
         end = timer()
         logger.info(f"Ground motion exceedance computation time: {round(end - start, 1)} seconds")
 
+        logger.info("Test3 run finished")
+
         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

src/shf_wrapper.py

@@ -58,6 +58,7 @@ def main(argv):
     parser.add_argument("--verbose", type=str2bool)
 
     args = parser.parse_args()
+    print(args)
     shf(**vars(args))
     return