diff --git a/debian/changelog b/debian/changelog
index 6868e97f592e..4bc123b65f3b 100644
--- a/debian/changelog
+++ b/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
diff --git a/openquake/calculators/tests/event_based_test.py b/openquake/calculators/tests/event_based_test.py
index 54e34c95e909..1f59d52bcde9 100644
--- a/openquake/calculators/tests/event_based_test.py
+++ b/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)
diff --git a/openquake/hazardlib/correlation.py b/openquake/hazardlib/correlation.py
index 0668ab4f08db..fc3a9b357389 100644
--- a/openquake/hazardlib/correlation.py
+++ b/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
diff --git a/openquake/qa_tests_data/event_based/case_31/__init__.py b/openquake/qa_tests_data/event_based/case_31/__init__.py
new file mode 100644
index 000000000000..e69de29bb2d1
diff --git a/openquake/qa_tests_data/event_based/case_31/expected/avg_gmf.csv b/openquake/qa_tests_data/event_based/case_31/expected/avg_gmf.csv
new file mode 100644
index 000000000000..27a5b9cb129f
--- /dev/null
+++ b/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
diff --git a/openquake/qa_tests_data/event_based/case_31/gmpe_logic_tree.xml b/openquake/qa_tests_data/event_based/case_31/gmpe_logic_tree.xml
new file mode 100644
index 000000000000..10bdddc339ff
--- /dev/null
+++ b/openquake/qa_tests_data/event_based/case_31/gmpe_logic_tree.xml
@@ -0,0 +1,37 @@
+
+
+
+
+
+
+
+
+ BooreEtAl2014
+ 0.6
+
+
+
+ CampbellBozorgnia2014
+ 0.4
+
+
+
+
+
+
+ BooreEtAl2014
+ 0.2
+
+
+
+ CampbellBozorgnia2014
+ 0.8
+
+
+
+
+
+
diff --git a/openquake/qa_tests_data/event_based/case_31/job.ini b/openquake/qa_tests_data/event_based/case_31/job.ini
new file mode 100644
index 000000000000..4ab4de3ae1dc
--- /dev/null
+++ b/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}
diff --git a/openquake/qa_tests_data/event_based/case_31/source_model.xml b/openquake/qa_tests_data/event_based/case_31/source_model.xml
new file mode 100644
index 000000000000..17faee8c0aaa
--- /dev/null
+++ b/openquake/qa_tests_data/event_based/case_31/source_model.xml
@@ -0,0 +1,114 @@
+
+
+
+
+
+
+
+
+
+
+ -5.0000000E-01 -5.0000000E-01 -3.0000000E-01 -1.0000000E-01 1.0000000E-01 2.0000000E-01 3.0000000E-01 -8.0000000E-01
+
+
+
+
+
+ 0.0000000E+00
+
+
+ 1.0000000E+01
+
+
+
+ WC1994
+
+
+ 1.0000000E+00
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ -6.0000000E-01 -6.0000000E-01 -4.0000000E-01 -2.0000000E-01 0.0000000E-01 1.0000000E-01 2.0000000E-01 -9.0000000E-01
+
+
+
+
+
+ 0.0000000E+00
+
+
+ 1.0000000E+01
+
+
+
+ WC1994
+
+
+ 1.0000000E+00
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/openquake/qa_tests_data/event_based/case_31/source_model_logic_tree.xml b/openquake/qa_tests_data/event_based/case_31/source_model_logic_tree.xml
new file mode 100644
index 000000000000..a489cd3c6d9d
--- /dev/null
+++ b/openquake/qa_tests_data/event_based/case_31/source_model_logic_tree.xml
@@ -0,0 +1,13 @@
+
+
+
+
+
+ source_model.xml
+ 1.0
+
+
+
+