Correlation modules refactor #11 – conditioned GMFs calculator - #11724
Merged
Conversation
Compute GSIM mean, between-event standard deviation, and within-event standard deviation for every conditioned target IMT. Index those arrays with the target IMT being processed instead of reusing the first IMT's values throughout a multi-IMT calculation. Add a focused regression test with IMT-dependent statistics and zero station residuals so that both the computed arrays and the conditioned means expose incorrect first-IMT reuse.
Remove the hardcoded PGV exclusions from both target and observed intensity measures. Ground-motion conditioning now excludes only MMI and lets the selected GSIM and correlation models validate the remaining IMTs. Use one shared filter for both paths and add a focused test showing that PGA and PGV are retained while MMI is excluded.
Flatten the station observation-error variances before adding them to the within-event covariance diagonal. The previous column-vector shape broadcast against the diagonal and could cause fill_diagonal to reuse the wrong variance for later observations. Add a two-station regression test with unequal uncertainties and compare the resulting inverse covariance with the explicitly constructed matrix.
Assemble one all-IMT covariance system for station observations. The matrix combines the selected within-event model, observation-error variances, and between-event residuals represented in the small unique-IMT latent space. Missing station observations are excluded consistently, and the station pseudoinverse is computed once for repeated solves. Add a deterministic Du-Ning test covering two IMTs, heterogeneous phi and tau values, unequal observation errors, residual ordering, and the independent between-event contribution.
Build the dense all-IMT target prior and target-station covariance blocks needed as a correctness reference for genuine joint correlation models. Expose posterior mean and covariance evaluation together with Matheron substitution against the precomputed station system. Verify the affine transformation deterministically: applying it to every column of a small joint Cholesky factor reproduces the direct Schur-complement covariance without Monte Carlo sampling error.
Route genuine spatial-cross-IMT models through one all-IMT Gaussian posterior. Draw a single unconditional target-and-station prior, apply Matheron substitution once, and preserve model-defined cross-IMT dependence instead of resetting the random stream for each target IMT. Keep the historical spatial-only path and finite truncated-MVN behavior unchanged. Update the memory preflight to account for the full joint covariance and factor, so oversized dense cases fail clearly until the structured sampler is available. Add a focused integration test that prevents fallback to the legacy per-IMT method.
Compute the exact all-IMT conditional mean in bounded target-site chunks. Joint mean-only calculations no longer allocate or share a complete target-station distance matrix; each chunk builds and applies only its own analytic covariance block. Choose the chunk size from the number of target IMTs, observed IMTs, and stations, and verify deterministically that multiple chunk sizes reproduce the fully materialized posterior mean.
Update the four conditioned QA cases affected by the correctness fixes. Multi-IMT calculations now use the GSIM mean, tau, and phi associated with each target IMT rather than the first target IMT's arrays. Cases with unequal station uncertainties also receive the variance belonging to each observation instead of a broadcast value from an earlier row. Keep all generated reference changes isolated so reviewers can distinguish intentional numerical updates from the conditioning implementation.
Use the joint conditioning path for every genuine spatial-cross-IMT model, including deterministic calculations when truncated_mvn is enabled. Reject finite truncated joint sampling explicitly until that algorithm is implemented instead of silently reverting to the legacy per-IMT approximation. Sample the dense posterior directly and fall back to a tolerance-checked eigendecomposition when it is positive semidefinite but singular. This supports collocated, zero-error station systems that are valid under pseudoinverse conditioning and avoids constructing a larger singular target-and-station prior.
Stop building and sharing the station-to-target distance matrix when a joint model uses the target-to-station block only. Keep deterministic joint conditioning bounded by the configured chunk size. Estimate dense covariance temporaries, target-station blocks, random samples, and float32 results before launching workers. This replaces the earlier two-matrix estimate, which could admit dense calculations that would exhaust memory, while avoiding false rejection of chunked mean-only calculations.
Cover incomplete station records, actual PGA-PGV conditioning with Du and Ning (2021), finite-truncation rejection, and singular station systems with deterministic matrix tests. Add an independently assembled two-IMT posterior with unequal phi and tau values and non-diagonal between-event correlation. Its expected station, target-station, target, posterior-mean, and posterior-covariance matrices do not reuse the implementation mapping arrays.
Create the conditioned computer event and realization arrays before estimating the random-sample workspace, since the estimate depends on the resulting event count. Preserve the existing initialization before worker submission by moving it rather than duplicating it. Validate the ordering with the South Napa DuNing2021 benchmark reduced to one deterministic field; it completed while reporting a 183.38 MB chunked joint-conditioning workspace.
Keep the historical float32 covariance blocks on the legacy per-IMT path, but allow the shared covariance builder to emit float64 arrays for dense joint station, target-station, and target systems. This avoids discarding model coefficient precision before the pseudoinverse and Schur-complement operations that define the joint reference posterior.
Build each deterministic joint-conditioning chunk from local positional indices rather than global site IDs. SiteCollection.filtered expects positions, so using IDs selected the wrong rows or failed after earlier site filtering produced non-contiguous IDs. Exercise multiple chunks on a target collection with site IDs 0, 2, 4, 6, and 8 and compare the result with the fully materialized posterior mean.
When the station covariance is singular, require the observed residual vector to lie in its column space before applying the pseudoinverse. Incompatible zero-error observations represent a probability-zero event and do not define a valid Gaussian conditional distribution. Also verify that each target-station covariance block lies in the station covariance row space. Test both rejection of contradictory collocated observations and acceptance of a compatible singular system.
Include four dense float64 station-scale matrices in the joint conditioning preflight for the within-event covariance, total station covariance, pseudoinverse, and matrix-product or factorization workspace. Apply the term to both random-field and chunked mean-only estimates, and add a station-rich case where the quadratic station system dominates the target arrays.
Rewrite the module overview around the residual model, calculation flow, and notation used by Engler et al. (2022), especially Appendix B. Use the paper's capital-tau and covariance names in the joint implementation, document the all-IMT generalization of equations B16 and B17, and map Engler's grid dimensions to OpenQuake's established array axes. This is a naming-only change with no numerical effect.
Give every conditioned-GMF notation entry its own definition and map both calculation paths to the relevant equations in Engler et al. (2022), Appendix B. Document OpenQuake's station-error, pseudoinverse, separable-covariance, clipping, and truncation choices so that the boundaries of the implementation are explicit. Verify that the all-IMT conditional-MVN implementation matches an independently assembled evaluation of equations B8-B9 and B16-B17. This commit does not change production numerical behavior.
Replace sentence-length conditioned-GMF test names with concise identifiers. Retain the less obvious intent in nearby comments so that test reports remain readable without losing the rationale for each regression check.
Remove the fixed eight-million-element conditioning block. Derive a whole-target-site chunk size from memory.conditioned_gmf_gb after reserving the estimated station-system and result workspaces and dividing the remaining budget across GSIM tasks. Pass the calculated block size from the admission check to the worker-side mean calculation so estimation and execution use the same plan. Report the selected site count, clarify the configuration comment and rejection message, and verify both adaptive and minimum one-site blocks.
Update scenario case 26 after conditioned targets began using the GSIM statistics for each requested IMT. PGA remains unchanged, while the SA fields now use their period-specific mean, tau, and phi values instead of reusing the PGA statistics. This expected output was missed by the earlier QA refresh because its numerical comparison is skipped on macOS.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Correct joint conditioning of ground-motion fields across sites and IMTs
Part of #11230. This PR establishes correct joint conditioning of ground-motion fields across sites and intensity measures. It allows the newly added spatial-cross-IMT within-event correlation models to be used without resorting to the historical per-IMT, separable approximation. It also fixes several correctness issues uncovered while checking the conditioned-GMF workflow and substantially expands the test coverage. It also begins to lay the groundwork for efficient memory management and performance improvements for calculations involving conditioning on observations.
Motivation
The conditioned-GMF calculator historically processed each target IMT independently. Cross-IMT dependence was approximated by multiplying a spatial correlation with a separate cross-IMT coefficient, and realizations for different target IMTs were not drawn from one joint posterior. The reasons for this were both performance considerations, as well as the absence of any joint spatial-cross-IMT correlation models in the engine.
With the recent addition of joint spatial-cross-IMT models such as Loth and Baker (2013), Markhvida et al. (2018), or Du and Ning (2021), the historical approximations are no longer necessary, and the calculator is rewired in this PR to be able to use the joint spatial-cross-IMT models that are now available in the engine. The existing approach also made it difficult to introduce scalable sampling algorithms because the conditioning logic and historical sampling path are tightly coupled.
During this work, several other pre-existing correctness issues were identified:
These issues are corrected here and covered by focused regression tests.
Joint spatial-cross-IMT conditioning architecture
For a genuine spatial-cross-IMT within-event model, the calculator now constructs one station system containing all usable observed IMTs and station sites. This system combines:
The target IMTs and sites are then treated as a single IMT-major Gaussian vector. The calculator constructs the total station, target-to-station, and target covariance blocks and applies the standard conditional multivariate-normal equations.
This is algebraically equivalent to the partitioned formulation in equations B8–B9 and B16–B17 of Engler et al. (2022), generalized to multiple jointly correlated target IMTs. A deterministic matrix test independently evaluates the Engler partition and confirms that it produces the same posterior mean and covariance as the new joint formulation.
The historical spatial-only path remains available and retains its existing per-IMT and finite truncated-normal behavior, but should be removed in a dedicated follow-up PR. This was not done in the current PR since the numbers for several tests would change.
Bugfixes and addition of numerical safeguards
The PR additionally:
Since this module has undergone several changes since its original implementation, the definitions of the notations used in the module have been updated. The module documentation now also provides an overview of the basic calculation flow.
Memory and performance groundwork
Deterministic joint calculations no longer construct the complete target-to-station or target-to-target matrices. The exact posterior mean is evaluated in target-site chunks.
The conditioned-calculation now distinguishes between:
It also initializes the event and realization arrays before estimating memory, so the estimate uses the actual number of fields instead of the rougher approximation used previously.
Oversized calculations now fail fast before allocating the full posterior rather than exhausting system memory.
For deterministic joint conditioning, the chunk size is derived from the configured
memory.conditioned_gmf_gbworkspace budget. The planner first reserves memory for the station covariance system and final result arrays, divides the remaining budget across GSIM tasks, and then determines how many complete target sites can fit while accounting for three temporary float64 covariance arrays. The calculated block size is passed to the worker, ensuring that memory estimation and execution use exactly the same plan.If even a single target-site chunk cannot fit, the calculation is rejected with an error reporting the required and configured memory. The user can then reduce the number of target or station sites or adjust
memory.conditioned_gmf_gbfor the machine running the calculation.Note that the configuration
memory.conditioned_gmf_gbrepresents the workspace available to conditioned-GMF calculations, not the machine’s total physical memory.QA changes
Expected results are updated for four existing conditioned QA cases. These changes are intentional consequences of the bug fixes described above:
Current limitations
What should come next
This work also prepares for scalable conditional-simulation methods. The immediate next step should be to investigate conditional simulation approaches for replacing construction of the dense posterior covariance with an unconditional joint field followed by a station-based correction. We'd need to evaluate scalable methods for generating that unconditional field, including circulant embedding as proposed by Bailey et al. (2021), alongside related memory-efficient approaches such as the successive conditional simulation algorithm of Verros et al. (2017).
Performance could be evaluated using an event such as 2014 M6.0 South Napa with a smaller footprint, and a stretch case such as the 2023 M7.8 Kahramanmaraş earthquake with a much larger footprint in
oq-risk-tests, targeting approximately:Once the joint sampling path is scalable and stable, the historical per-IMT conditioning implementation and its separable approximation can be removed.