Trying to fix the HM2018 model

debian/changelog

@@ -1,4 +1,6 @@
   [Michele Simionato]
+  * Fixed a bug in calculations with a filtered site collection using the
+    HM2018CorrelationModel
   * Internal: raised a clear error message when get_composite_source_model is
    called without passing a datastore in presence of multifault sources
 

openquake/calculators/tests/event_based_test.py

@@ -40,7 +40,7 @@
     blocksize, case_1, case_2, case_3, case_4, case_5, case_6, case_7,
     case_8, case_9, case_10, case_12, case_13, case_14, case_15, case_16,
     case_17,  case_18, case_19, case_20, case_21, case_22, case_23, case_24,
-    case_25, case_26, case_27, case_28, case_29, case_30, src_mutex)
+    case_25, case_26, case_27, case_28, case_29, case_30, case_31, src_mutex)
 from openquake.qa_tests_data.event_based.spatial_correlation import (
     case_1 as sc1, case_2 as sc2, case_3 as sc3)
 
@@ -597,3 +597,9 @@ def test_30(self):
         out = self.run_calc(case_30.__file__, 'job.ini', exports='csv')
         [fname] = out['ruptures', 'csv']
         self.assertEqualFiles('expected/ruptures.csv', fname, delta=1E-6)
+
+    def test_31(self):
+        # HM2018CorrelationModel with filtered site collection
+        self.run_calc(case_31.__file__, 'job.ini', exports='csv')
+        [f] = export(('avg_gmf', 'csv'), self.calc.datastore)
+        self.assertEqualFiles('expected/avg_gmf.csv', f)

openquake/hazardlib/correlation.py

@@ -40,7 +40,7 @@ def apply_correlation(self, sites, imt, residuals, stddev_intra=0):
             represents sites (the length as ``sites`` parameter) and
             second one represents different realizations (samples).
         :param stddev_intra:
-            Intra-event standard deviation array. Note that different sites do
+            Intra-event standard deviation array (phi). Different sites do
             not necessarily have the same intra-event standard deviation.
         :returns:
             Array of the same structure and semantics as ``residuals``
@@ -171,17 +171,14 @@ def _get_correlation_matrix(self, sites, imt):
 
     def apply_correlation(self, sites, imt, residuals, stddev_intra):
         """
-        Apply correlation to randomly sampled residuals.
-
-        See Parent function
+        Apply correlation to randomly sampled residuals
         """
-        # stddev_intra is repeated if there is only one value
-        if len(stddev_intra) == 1:
-            stddev_intra = numpy.full(len(sites.complete), stddev_intra)
-        # Reshape 'stddev_intra' if needed
-        stddev_intra = stddev_intra.squeeze()
-        if not stddev_intra.shape:
-            stddev_intra = stddev_intra[None]
+        # TODO: the case of filtered sites is probably managed incorrectly
+        # NB: this is SLOW and we cannot use the cache as in JB2009 because
+        # we are not using the complete site collection
+        nsites = len(sites)
+        assert len(residuals) == len(stddev_intra) == nsites
+        D = numpy.diag(stddev_intra)  # phi as a diagonal matrix
 
         if self.uncertainty_multiplier == 0:   # No uncertainty
 
@@ -190,37 +187,30 @@ def apply_correlation(self, sites, imt, residuals, stddev_intra):
             # normalized, sampled from a standard normal distribution.
             # For this, every row of 'residuals' (every site) is divided by its
             # corresponding standard deviation element.
-            residuals_norm = residuals / stddev_intra[sites.sids, None]
-
-            # Lower diagonal of the Cholesky decomposition from/to cache
-            try:
-                cormaLow = self.cache[imt]
-            except KeyError:
-                # Note that instead of computing the whole correlation matrix
-                # corresponding to sites.complete, here we compute only the
-                # correlation matrix corresponding to sites.
-                cormaLow = numpy.linalg.cholesky(
-                       numpy.diag(stddev_intra[sites.sids]) @
-                       self._get_correlation_matrix(sites, imt) @
-                       numpy.diag(stddev_intra[sites.sids]))
-                self.cache[imt] = cormaLow
+            residuals_norm = residuals / stddev_intra[:, None]
+
+            # Lower diagonal of the Cholesky decomposition
+            # Note that instead of computing the whole correlation matrix
+            # corresponding to sites.complete, here we compute only the
+            # correlation matrix corresponding to sites
+            cormaLow = numpy.linalg.cholesky(
+                D @ self._get_correlation_matrix(sites, imt) @ D)
 
             # Apply correlation
-            return numpy.dot(cormaLow, residuals_norm)
+            return cormaLow @ residuals_norm
 
         else:   # Variability (uncertainty) is included
-            nsim = len(residuals[1])
-            nsites = len(residuals)
+            nsim = residuals.shape[1]
 
             # Re-sample all the residuals
             residuals_correlated = residuals * 0
             for isim in range(0, nsim):
+                # FIXME: the seed is not set!
                 corma = self._get_correlation_matrix(sites, imt)
-                cov = (numpy.diag(stddev_intra[sites.sids]) @ corma @
-                       numpy.diag(stddev_intra[sites.sids]))
+                # NB: corma is different at each loop since contains randomicity
                 residuals_correlated[0:, isim] = (
                     numpy.random.multivariate_normal(
-                        numpy.zeros(nsites), cov, 1))
+                        numpy.zeros(nsites), D @ corma @ D, 1))
 
             return residuals_correlated
 

openquake/qa_tests_data/event_based/case_31/expected/avg_gmf.csv

@@ -0,0 +1,102 @@
+#,,,,"generated_by='OpenQuake engine 3.18.0-git6ee84744a4', start_date='2023-10-18T08:46:31', checksum=193965386"
+site_id,lon,lat,gmv_SA(0.3),gsd_SA(0.3)
+0,-2.00061E+00,2.02920E-01,1.66161E-08,3.83684E+03
+1,-2.00061E+00,-2.46750E-01,7.08940E-08,6.33337E+03
+2,-2.00059E+00,6.52580E-01,1.04708E-10,2.66500E+00
+3,-2.00058E+00,-6.96410E-01,2.63923E-08,4.25542E+03
+4,-2.00053E+00,1.10224E+00,0.00000E+00,0.00000E+00
+5,-2.00053E+00,-1.14607E+00,1.47151E-09,6.14167E+02
+6,-2.00045E+00,1.55190E+00,0.00000E+00,0.00000E+00
+7,-2.00044E+00,-1.59573E+00,1.74891E-10,2.31377E+01
+8,-2.00034E+00,2.00156E+00,0.00000E+00,0.00000E+00
+9,-2.00033E+00,-2.04539E+00,0.00000E+00,0.00000E+00
+10,-1.55095E+00,2.02920E-01,2.08466E-03,1.91645E+02
+11,-1.55095E+00,-2.46750E-01,1.02306E-02,3.25967E+00
+12,-1.55090E+00,6.52580E-01,2.29330E-05,4.75726E+03
+13,-1.55089E+00,-6.96410E-01,9.26077E-03,6.89180E+00
+14,-1.55079E+00,1.10224E+00,1.96522E-10,3.25687E+01
+15,-1.55077E+00,-1.14607E+00,2.11505E-04,1.71172E+03
+16,-1.55062E+00,1.55190E+00,0.00000E+00,0.00000E+00
+17,-1.55061E+00,-1.59573E+00,1.23530E-07,7.75930E+03
+18,-1.55041E+00,2.00156E+00,0.00000E+00,0.00000E+00
+19,-1.55038E+00,-2.04539E+00,2.24175E-10,4.33195E+01
+20,-1.10129E+00,2.02920E-01,1.74861E-02,2.27762E+00
+21,-1.10128E+00,-2.46750E-01,2.21503E-02,2.58062E+00
+22,-1.10121E+00,6.52580E-01,5.44196E-03,3.57584E+01
+23,-1.10119E+00,-6.96410E-01,1.89384E-02,2.64911E+00
+24,-1.10104E+00,1.10224E+00,1.15034E-05,6.17678E+03
+25,-1.10102E+00,-1.14607E+00,1.30415E-02,2.48723E+00
+26,-1.10080E+00,1.55190E+00,1.06558E-10,3.01743E+00
+27,-1.10077E+00,-1.59573E+00,5.01971E-04,7.97006E+02
+28,-1.10047E+00,2.00156E+00,0.00000E+00,0.00000E+00
+29,-1.10043E+00,-2.04539E+00,6.65613E-08,6.19559E+03
+30,-6.51620E-01,2.02920E-01,2.94546E-02,2.31100E+00
+31,-6.51620E-01,-2.46750E-01,5.01483E-02,2.55007E+00
+32,-6.51520E-01,6.52580E-01,1.38517E-02,2.01028E+00
+33,-6.51500E-01,-6.96410E-01,3.85551E-02,2.56882E+00
+34,-6.51300E-01,1.10224E+00,1.14091E-03,3.25311E+02
+35,-6.51270E-01,-1.14607E+00,2.22350E-02,2.27361E+00
+36,-6.50970E-01,1.55190E+00,1.32426E-08,3.15846E+03
+37,-6.50940E-01,-1.59573E+00,1.02663E-02,5.49126E+00
+38,-6.50530E-01,2.00156E+00,0.00000E+00,0.00000E+00
+39,-6.50490E-01,-2.04539E+00,4.43623E-06,8.06900E+03
+40,-2.01960E-01,2.02920E-01,4.44669E-02,2.49045E+00
+41,-2.01950E-01,-2.46750E-01,1.36568E-01,2.77572E+00
+42,-2.01830E-01,6.52580E-01,1.96015E-02,2.24945E+00
+43,-2.01810E-01,-6.96410E-01,8.92375E-02,2.96026E+00
+44,-2.01560E-01,1.10224E+00,2.47168E-03,1.38391E+02
+45,-2.01520E-01,-1.14607E+00,2.99529E-02,2.31927E+00
+46,-2.01150E-01,1.55190E+00,6.95609E-08,5.94987E+03
+47,-2.01100E-01,-1.59573E+00,1.28681E-02,2.23196E+00
+48,-2.00600E-01,2.00156E+00,0.00000E+00,0.00000E+00
+49,-2.00540E-01,-2.04539E+00,4.95191E-05,3.73854E+03
+50,2.47710E-01,2.02920E-01,3.80965E-02,2.53644E+00
+51,2.47710E-01,-2.46750E-01,6.51902E-02,2.57147E+00
+52,2.47860E-01,6.52580E-01,1.66356E-02,2.44861E+00
+53,2.47890E-01,-6.96410E-01,5.34652E-02,2.59047E+00
+54,2.48190E-01,1.10224E+00,1.13762E-03,4.23517E+02
+55,2.48230E-01,-1.14607E+00,2.22328E-02,2.44921E+00
+56,2.48680E-01,1.55190E+00,7.28563E-08,6.31607E+03
+57,2.48740E-01,-1.59573E+00,1.19628E-02,2.29466E+00
+58,2.49340E-01,2.00156E+00,0.00000E+00,0.00000E+00
+59,2.49410E-01,-2.04539E+00,1.74536E-05,5.88702E+03
+60,6.97370E-01,2.02920E-01,1.95185E-02,2.46287E+00
+61,6.97380E-01,-2.46750E-01,2.32399E-02,2.24467E+00
+62,6.97550E-01,6.52580E-01,1.10836E-02,2.34874E+00
+63,6.97580E-01,-6.96410E-01,2.21145E-02,2.27606E+00
+64,6.97930E-01,1.10224E+00,1.31575E-05,6.37823E+03
+65,6.97980E-01,-1.14607E+00,1.40843E-02,2.33929E+00
+66,6.98500E-01,1.55190E+00,1.05567E-09,4.24703E+02
+67,6.98570E-01,-1.59573E+00,8.09910E-03,5.31183E+00
+68,6.99270E-01,2.00156E+00,0.00000E+00,0.00000E+00
+69,6.99360E-01,-2.04539E+00,3.81603E-08,4.76493E+03
+70,1.14703E+00,2.02920E-01,6.85712E-03,1.62976E+01
+71,1.14704E+00,-2.46750E-01,1.16328E-02,2.20476E+00
+72,1.14724E+00,6.52580E-01,2.20384E-04,1.60423E+03
+73,1.14728E+00,-6.96410E-01,1.17465E-02,2.27644E+00
+74,1.14768E+00,1.10224E+00,4.15107E-08,4.70646E+03
+75,1.14773E+00,-1.14607E+00,1.03014E-02,2.08519E+00
+76,1.14833E+00,1.55190E+00,1.12186E-10,4.10427E+00
+77,1.14841E+00,-1.59573E+00,7.78538E-07,8.60038E+03
+78,1.14921E+00,2.00156E+00,0.00000E+00,0.00000E+00
+79,1.14930E+00,-2.04539E+00,5.71386E-10,2.13829E+02
+80,1.59670E+00,2.02920E-01,1.52223E-06,8.73191E+03
+81,1.59671E+00,-2.46750E-01,2.63000E-05,4.79965E+03
+82,1.59693E+00,6.52580E-01,3.35592E-09,1.21846E+03
+83,1.59697E+00,-6.96410E-01,9.37763E-06,6.24420E+03
+84,1.59742E+00,1.10224E+00,0.00000E+00,0.00000E+00
+85,1.59748E+00,-1.14607E+00,1.27214E-08,3.10437E+03
+86,1.59816E+00,1.55190E+00,0.00000E+00,0.00000E+00
+87,1.59824E+00,-1.59573E+00,5.52216E-10,1.99022E+02
+88,1.59914E+00,2.00156E+00,0.00000E+00,0.00000E+00
+89,1.59925E+00,-2.04539E+00,0.00000E+00,0.00000E+00
+90,2.04636E+00,2.02920E-01,0.00000E+00,0.00000E+00
+91,2.04637E+00,-2.46750E-01,0.00000E+00,0.00000E+00
+92,2.04662E+00,6.52580E-01,0.00000E+00,0.00000E+00
+93,2.04666E+00,-6.96410E-01,1.11555E-10,3.83227E+00
+94,2.04716E+00,1.10224E+00,0.00000E+00,0.00000E+00
+95,2.04723E+00,-1.14607E+00,0.00000E+00,0.00000E+00
+96,2.04798E+00,1.55190E+00,0.00000E+00,0.00000E+00
+97,2.04808E+00,-1.59573E+00,0.00000E+00,0.00000E+00
+98,2.04908E+00,2.00156E+00,0.00000E+00,0.00000E+00
+99,2.04920E+00,-2.04539E+00,0.00000E+00,0.00000E+00

openquake/qa_tests_data/event_based/case_31/gmpe_logic_tree.xml

@@ -0,0 +1,37 @@
+<?xml version="1.0" encoding="UTF-8"?>
+
+<nrml xmlns:gml="http://www.opengis.net/gml"
+      xmlns="http://openquake.org/xmlns/nrml/0.4">
+    <logicTree logicTreeID='lt1'>
+
+            <logicTreeBranchSet uncertaintyType="gmpeModel" branchSetID="bs1"
+                    applyToTectonicRegionType="Active Shallow Crust">
+
+                <logicTreeBranch branchID="BA">
+                    <uncertaintyModel>BooreEtAl2014</uncertaintyModel>
+                    <uncertaintyWeight>0.6</uncertaintyWeight>
+                </logicTreeBranch>
+
+                <logicTreeBranch branchID="TA">
+                    <uncertaintyModel>CampbellBozorgnia2014</uncertaintyModel>
+                    <uncertaintyWeight>0.4</uncertaintyWeight>
+                </logicTreeBranch>
+
+            </logicTreeBranchSet>
+            <logicTreeBranchSet uncertaintyType="gmpeModel" branchSetID="bs2"
+                    applyToTectonicRegionType="Subduction Interface">
+
+                <logicTreeBranch branchID="BA">
+                    <uncertaintyModel>BooreEtAl2014</uncertaintyModel>
+                    <uncertaintyWeight>0.2</uncertaintyWeight>
+                </logicTreeBranch>
+
+                <logicTreeBranch branchID="TA">
+                    <uncertaintyModel>CampbellBozorgnia2014</uncertaintyModel>
+                    <uncertaintyWeight>0.8</uncertaintyWeight>
+                </logicTreeBranch>
+
+            </logicTreeBranchSet>
+
+    </logicTree>
+</nrml>

openquake/qa_tests_data/event_based/case_31/job.ini

@@ -0,0 +1,37 @@
+[general]
+
+description = HM2018 with a filtered site collection
+calculation_mode = event_based
+ses_seed = 24
+
+[geometry]
+region = -2.0 -2.0, -2.0 2.0, 2.0 2.0, 2.0 -2.0
+region_grid_spacing = 50.0
+
+[logic_tree]
+number_of_logic_tree_samples = 0
+
+[erf]
+rupture_mesh_spacing = 5
+width_of_mfd_bin = 1.0
+area_source_discretization = 10.0
+
+[site_params]
+reference_vs30_type = measured
+reference_vs30_value = 600.0
+reference_depth_to_2pt5km_per_sec = 5.0
+reference_depth_to_1pt0km_per_sec = 100.0
+
+[calculation]
+source_model_logic_tree_file = source_model_logic_tree.xml
+gsim_logic_tree_file = gmpe_logic_tree.xml
+investigation_time = 50.0
+intensity_measure_types = SA(0.3)
+truncation_level = 3
+maximum_distance = 200.0
+minimum_magnitude = 5.5
+
+[event_based_params]
+ses_per_logic_tree_path = 1
+ground_motion_correlation_model = HM2018
+ground_motion_correlation_params =  {'uncertainty_multiplier': 0}