From 792bef067ea421ba0ca94788285fcb066fe946e4 Mon Sep 17 00:00:00 2001
From: Marco Pagani <marco.pagani@globalquakemodel.org>
Date: Tue, 13 Aug 2024 18:37:19 +0200
Subject: [PATCH] Add example implementation of basin correction

---
 .../hazardlib/gsim/mgmpe/modifiable_gmpe.py   | 18 +++++++++++
 .../tests/gsim/mgmpe/modifiable_gmpe_test.py  | 31 +++++++++++++++++++
 2 files changed, 49 insertions(+)

diff --git a/openquake/hazardlib/gsim/mgmpe/modifiable_gmpe.py b/openquake/hazardlib/gsim/mgmpe/modifiable_gmpe.py
index c3f431ce3b55..72366c82365c 100644
--- a/openquake/hazardlib/gsim/mgmpe/modifiable_gmpe.py
+++ b/openquake/hazardlib/gsim/mgmpe/modifiable_gmpe.py
@@ -41,6 +41,17 @@
 
 # ################ BEGIN FUNCTIONS MODIFYING mean_stds ################## #
 
+def m9_basin_term(ctx, imt, me, si, ta, phi):
+    """
+    Implements the NSHM "M9" basin term to be applied to other GMMs for
+    Cascadia.
+    """
+    if imt.period > 1.9:
+        fb_m9 = np.log(2.0)
+        idx = ctx.z2pt5 > 6.0
+        me[idx] += fb_m9
+
+
 def sigma_model_alatik2015(ctx, imt, me, si, ta, ph,
                            ergodic, tau_model, phi_ss_coetab, tau_coetab):
     """
@@ -223,6 +234,7 @@ class ModifiableGMPE(GMPE):
     DEFINED_FOR_REFERENCE_VELOCITY = None
 
     def __init__(self, **kwargs):
+
         # Create the original GMPE
         [(gmpe_name, kw)] = kwargs.pop('gmpe').items()
         self.params = kwargs  # non-gmpe parameters
@@ -248,6 +260,12 @@ def __init__(self, **kwargs):
             setattr(self, 'DEFINED_FOR_STANDARD_DEVIATION_TYPES',
                     {StdDev.TOTAL, StdDev.INTRA_EVENT, StdDev.INTER_EVENT})
 
+        if ('m9_basin_term' in self.params) and (
+                'z2pt5' not in self.gmpe.REQUIRES_SITES_PARAMETERS):
+            tmp = list(self.gmpe.REQUIRES_SITES_PARAMETERS)
+            tmp.append('z2pt5')
+            self.gmpe.REQUIRES_SITES_PARAMETERS = frozenset(tmp)
+
         # This is required by the `sigma_model_alatik2015` function
         key = 'sigma_model_alatik2015'
         if key in self.params:
diff --git a/openquake/hazardlib/tests/gsim/mgmpe/modifiable_gmpe_test.py b/openquake/hazardlib/tests/gsim/mgmpe/modifiable_gmpe_test.py
index cc353f8b2e98..035dec23c41a 100644
--- a/openquake/hazardlib/tests/gsim/mgmpe/modifiable_gmpe_test.py
+++ b/openquake/hazardlib/tests/gsim/mgmpe/modifiable_gmpe_test.py
@@ -213,3 +213,34 @@ def test(self):
 
             # Check the computed intra-event stdev
             aae(phi[0], [exp_stdev])
+
+
+class ModifiableGMPEBasinTerm(unittest.TestCase):
+
+    def get_ctx(self, cmaker):
+        ctx = cmaker.new_ctx(4)
+        ctx.mag = 6.0
+        ctx.hypo_depth = 10.
+        ctx.rhypo = ctx.rrup = ctx.repi = np.array([1., 10., 30., 70.])
+        ctx.z2pt5 = np.array([1.0, 2.0, 6.0, 7.0])
+        ctx.vs30 = np.array([760, 760, 760, 760])
+        return ctx
+
+    def test_basin(self):
+
+        # Create the Modifiable GMPE and original GMM
+        gmm_name = 'ParkerEtAl2020SInter'
+        gmm = ModifiableGMPE(gmpe={gmm_name: {}},
+                             m9_basin_term={})
+        gmpe = valid.gsim(gmm_name)
+
+        # Context maker
+        cmaker = simple_cmaker([gmm, gmpe], ['SA(2.0)'])
+        ctx = self.get_ctx(cmaker)
+
+        # Compute ground motion
+        mea, _, _, _ = cmaker.get_mean_stds([ctx])
+
+        # Check the computed mean
+        aae(mea[0, 0, 0:3], mea[1, 0, 0:3])
+        aae(mea[0, 0, 3:], mea[1, 0, 3:] + 0.6931471806)