Merge pull request #8369 from gem/cgmfs

Address pending review comments from #8318

openquake/calculators/event_based.py

@@ -443,10 +443,12 @@ def execute(self):
         proxies = [RuptureProxy(rec, scenario)
                    for rec in dstore['ruptures'][:]]
         if "station_data" in oq.inputs:
-            # this is meant to be used in scenario calculations with a single
-            # rupture; we are taking the first copy of the rupture (remember:
-            # _read_scenario_ruptures makes num_gmfs copies to parallelize)
-            # TODO: this is ugly andmust be improved upon!
+            # this is meant to be used in conditioned scenario calculations 
+            # with a single rupture; we are taking the first copy of the rupture 
+            # (remember: _read_scenario_ruptures makes num_gmfs copies in order to 
+            # parallelize, but the conditioning process is computationally quite 
+            # expensive, so we want to avoid repeating it num_gmfs times)
+            # TODO: this is ugly and must be improved upon!
             proxies = proxies[0:1]
         full_lt = self.datastore['full_lt']
         dstore.swmr_on()  # must come before the Starmap

openquake/hazardlib/shakemap/conditioned_gmfs.py

@@ -52,19 +52,39 @@
 # simulate the requested number of ground motion fields
 
 # Notation:
-# K: number of target sites at which ground motion is to be estimated
-# L: number of recording stations with ground motion observations
-# _D: subscript denoting data
-# _N: subscript denoting "native" data, for the target IMT
-# _NN: subscript denoting "nonnative" data, other IMTs at the stations
-# M: number of nonnative IMTs used in the conditioning process
-# Y_N: station data for the target IMT
-# Y_NN: station data for the 'M' nonnative IMTs
-# y2 = observations
-# φ = phi = within-event standard deviation
-# τ = tau = between-event standard deviation
-# Η (capital Greek η, not Latin H) = normalized between-event residual
-# δB = between-event residual, or bias
+# _D: subscript refers to the "Data" or observations
+# _Y: subscript refers to the target sites
+# yD: recorded values at the stations
+# var_addon_D: additional sigma for the observations that are uncertain, 
+#   which might arise if the values for this particular IMT were not directly
+#   recorded, but obtained by conversion equations or cross-correlation functions
+# mu_yD: predicted mean intensity at the observation points, from the specified GMM(s)
+# phi_D: predicted within-event uncertainty at the observation points, from the specified GMM(s)
+# tau_D: predicted between-event uncertainty at the observation points, from the specified GMM(s)
+# zeta_D: raw residuals at the observation points
+# cov_WD_WD: station data within-event covariance matrix, with the additional 
+#   variance of the residuals for the cases where the station data is uncertain
+# cov_WD_WD_inv: (pseudo)-inverse of the station data within-event covariance matrix
+# corr_HD_HD: cross-intensity measure correlations for the observed intensity measures 
+# mu_HD_yD: posterior mean of the (normalized) between-event residual
+# cov_HD_HD_yD: posterior covariance of the (normalized) between-event residual
+# mu_BD_yD: posterior mean of the between-event residual
+# cov_BD_BD_yD: posterior covariance of the conditional between-event residual
+# nominal_bias_mean: mean of mu_BD_yD, useful as a single value measure of the event bias,
+#   particularly in the heteroscedastic case
+# nominal_bias_stddev: sqrt of the mean of cov_BD_BD_yD
+# mu_Y: redicted mean of the intensity at the target sites
+# phi_Y: predicted within-event standard deviation of the intensity at the target sites
+# tau_Y: predicted between-event standard deviation of the intensity at the target sites
+# mu_BY_yD: mean of the conditional between-event residual for the target sites
+# cov_WY_WD and cov_WD_WY: within-event covariance matrices for the target sites and observation sites
+# cov_WY_WY: apriori within-event covariance matrix for the target sites
+# RC: regression coefficient matrix ("RC" = cov_WY_WD × cov_WD_WD_inv)
+# C: scaling matrix for the conditioned between-event covariance matrix
+# cov_WY_WY_wD: conditioned within-event covariance matrix for the target sites 
+# cov_BY_BY_yD: "conditioned between-event" covariance matrix for the target sites
+# mu_Y_yD: conditioned mean of the ground motion at the target sites
+# cov_Y_Y_yD: conditional covariance of the ground motion at the target sites
 
 
 class NoInterIntraStdDevs(Exception):
@@ -291,6 +311,9 @@ def get_conditioned_mean_and_covariance(
                       if imt_str not in ["MMI", "PGV"]}
     observed_imts = sorted([imt.from_string(imt_str)
                             for imt_str in observed_imtls])
+
+    # Generate the contexts and calculate the means and 
+    # standard deviations at the *station* sites ("_D")
     cmaker_D = ContextMaker(
         rupture.tectonic_region_type, [gmm],
         dict(truncation_level=0, imtls=observed_imtls,
@@ -466,13 +489,13 @@ def get_conditioned_mean_and_covariance(
         # Get the nominal bias and its variance as the means of the
         # conditional between-event residual mean and covariance
         nominal_bias_mean = numpy.mean(mu_BD_yD)
-        # ASK: why are we computing nominal_bias_stddev if it is not used??
         nominal_bias_stddev = numpy.sqrt(numpy.mean(cov_BD_BD_yD))
         logging.info(
-            "GMM: %s, IMT: %s, Nominal bias (mean): %.3f",
-            gmm_name, target_imt.string, nominal_bias_mean)
+            "GMM: %s, IMT: %s, Nominal bias mean: %.3f, Nominal bias stddev: %.3f",
+            gmm_name, target_imt.string, nominal_bias_mean, nominal_bias_stddev)
 
-        # From the GMMs, get the mean and stddevs at the target sites
+        # Generate the contexts and calculate the means and 
+        # standard deviations at the *target* sites ("_Y")
         cmaker_Y = ContextMaker(
             rupture.tectonic_region_type, [gmm], dict(
                 truncation_level=0, imtls={target_imt.string: [0]},