Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Address pending review comments from #8318 #8369

Merged
merged 4 commits into from Jan 27, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
10 changes: 6 additions & 4 deletions openquake/calculators/event_based.py
Expand Up @@ -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
Expand Down
57 changes: 40 additions & 17 deletions openquake/hazardlib/shakemap/conditioned_gmfs.py
Expand Up @@ -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):
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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]},
Expand Down