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.

Sign up for GitHub

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

Revert some variable names to match described notation #10430

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Revert some variable names to match described notation #10430

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
3070a50
Revert some variable names to match described notation
raoanirudh Mar 7, 2025
a818588
Merge branch 'master' into cond-refact
raoanirudh Mar 8, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
55 changes: 33 additions & 22 deletions openquake/hazardlib/calc/conditioned_gmfs.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
# -*- coding: utf-8 -*-
#
# Copyright (C) 2025 GEM Foundation
# Copyright (C) 2023-2025 GEM Foundation
#
# OpenQuake is free software: you can redistribute it and/or modify it
# under the terms of the GNU Affero General Public License as published
Expand Down Expand Up @@ -80,7 +80,7 @@
nominal_bias_stddev:
sqrt of the mean of cov_BD_BD_yD
mu_Y:
redicted mean of the intensity at the target sites
predicted mean of the intensity at the target sites
phi_Y:
predicted within-event standard deviation of the intensity at the target sites
tau_Y:
Expand Down Expand Up @@ -235,9 +235,9 @@ class TempResult:
native_data_available: bool
corr_HD_HD: numpy.ndarray = 0
cov_WD_WD_inv: numpy.ndarray = 0
D: numpy.ndarray = 0
phi_D_diag: numpy.ndarray = 0
T_D: numpy.ndarray = 0
zD: numpy.ndarray = 0
zeta_D: numpy.ndarray = 0


def _create_result(g, m, target_imt, observed_imts, station_data_filtered):
Expand Down Expand Up @@ -332,11 +332,11 @@ def create_result(g, m, target_imt, target_imts, observed_imts,
i * nss: (i + 1) * nss, 0]

# The raw residuals
t.zD = yD - mu_yD
t.D = numpy.diag(phi_D.flatten())
t.zeta_D = yD - mu_yD
t.phi_D_diag = numpy.diag(phi_D.flatten())

cov_WD_WD = compute_cov(station_sitecol, station_sitecol,
t.conditioning_imts, t.conditioning_imts, t.D, t.D)
t.conditioning_imts, t.conditioning_imts, t.phi_D_diag, t.phi_D_diag)

# Add on the additional variance of the residuals
# for the cases where the station data is uncertain
Expand Down Expand Up @@ -405,7 +405,7 @@ def get_mu_tau_phi(target_imt, gsim, mean_stds,
+ numpy.linalg.pinv(r.corr_HD_HD))

mu_HD_yD = numpy.linalg.multi_dot(
[cov_HD_HD_yD, r.T_D.T, r.cov_WD_WD_inv, r.zD])
[cov_HD_HD_yD, r.T_D.T, r.cov_WD_WD_inv, r.zeta_D])

# Compute the distribution of the conditional between-event
# residual B|Y2=y2
Expand All @@ -426,29 +426,40 @@ def get_mu_tau_phi(target_imt, gsim, mean_stds,

# Predicted uncertainty components at the target sites, from GSIM
tau_Y = mean_stds[2, 0][:, None]
Y = numpy.diag(mean_stds[3, 0])
phi_Y_diag = numpy.diag(mean_stds[3, 0])

# Compute the mean of the conditional between-event residual B|YD=yD
# for the target sites; the shapes are (nsites, nstations),
# (nstations, nsites), (nsites, nsites) respectively
# Compute the within-event covariance matrices for the
# target sites and observation sites; the shapes are
# (nsites, nstations) and (nstations, nsites) respectively
cov_WY_WD = compute_cov(target_sitecol, station_sitecol,
[target_imt], r.conditioning_imts, Y, r.D)
[target_imt], r.conditioning_imts, phi_Y_diag, r.phi_D_diag)
cov_WD_WY = compute_cov(station_sitecol, target_sitecol,
r.conditioning_imts, [target_imt], r.D, Y)
cov_WY_WY = compute_cov(target_sitecol, target_sitecol,
[target_imt], [target_imt], Y, Y)
r.conditioning_imts, [target_imt], r.phi_D_diag, phi_Y_diag)

# Compute the regression coefficient matrix [cov_WY_WD × cov_WD_WD_inv]
RC = cov_WY_WD @ r.cov_WD_WD_inv # shape (nsites, nstations)

# compute the mean, shape (nsites, 1)
mu = mu_Y + tau_Y @ mu_HD_yD[0, None] + RC @ (r.zD - mu_BD_yD)
# # Compute the mean of the conditional between-event residual B|YD=yD
# # for the target sites
# mu_HN_yD = mu_HD_yD[0, None]
# mu_BY_yD = tau_Y @ mu_HN_yD

# Compute the conditioned mean of the ground motion
# at the target sites; shape (nsites, 1)
mu_Y_yD = mu_Y + tau_Y @ mu_HD_yD[0, None] + RC @ (r.zeta_D - mu_BD_yD)

# Compute the within-event covariance matrix for the
# target sites (apriori) (nsites, nsites)
cov_WY_WY = compute_cov(target_sitecol, target_sitecol,
[target_imt], [target_imt], phi_Y_diag, phi_Y_diag)

# covariance matrices can contain extremely small negative values
# Both conditioned covariance matrices can contain extremely
# small negative values due to limitations of floating point
# operations (~ -10^-17 to -10^-15), these are clipped to zero

# Compute the conditioned within-event covariance matrix
# for the target sites clipped to zero, shape (nsites, nsites)
tau = (cov_WY_WY - RC @ cov_WD_WY).clip(min=0)
cov_WY_WY_wD = (cov_WY_WY - RC @ cov_WD_WY).clip(min=0)

# Compute the scaling matrix "C" for the conditioned between-event
# covariance matrix
Expand All @@ -460,8 +471,8 @@ def get_mu_tau_phi(target_imt, gsim, mean_stds,

# Compute the conditioned between-event covariance matrix
# for the target sites clipped to zero, shape (nsites, nsites)
phi = numpy.linalg.multi_dot([C, cov_HD_HD_yD, C.T]).clip(min=0)
return {(r.g, r.m): (mu, tau, phi, msg)}
cov_BY_BY_yD = numpy.linalg.multi_dot([C, cov_HD_HD_yD, C.T]).clip(min=0)
return {(r.g, r.m): (mu_Y_yD, cov_WY_WY_wD, cov_BY_BY_yD, msg)}


def get_me_ta_ph(cmaker, sdata, observed_imts, target_imts,
Expand Down