--- title: "Spectral matching workflow: selection, shaping, and suite scaling" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Spectral matching workflow: selection, shaping, and suite scaling} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", eval = FALSE, purl = FALSE) ``` This vignette documents the complete spectral-matching workflow: how a suite of recorded ground motions is selected from a processed database, adjusted so that its response spectra match a target uniform-hazard spectrum (UHS), and delivered as traceable file products. Two runners drive the workflow from one JSON file per site: ```{r overview} library(gmsp) runStage0(file = "gmsp/match.json", root = "/path/to/project") runMatch(file = "gmsp/match.json", root = "/path/to/project") ``` ## Method overview The workflow has three stages. **Spectral preselection (`runStage0()`).** Every candidate record in the processed database pool is scored by the compatibility of its PSA spectrum (H1 component) against the target band defined by the UHS mean and its lower and upper envelope quantiles. The `N` most compatible records form the suite. Optional seismological screens (magnitude, distance, provider, one record per event) restrict the pool before ranking. **Modal shaping (per record).** Each selected record is decomposed into intrinsic mode functions (see `TS2IMF()`), the retained modes are re-weighted by positive modal factors (`shape.modalFactor`) to improve spectral compatibility, and the reconstruction is normalized so the quadratic acceleration norm of the original record is preserved. The in-memory algorithm is exported as `fitModalFactor()`. **Record scaling (suite).** One positive coefficient $c_i$ per record (`shape.recordFactor`) scales the normalized reconstructed history. The in-memory algorithm is exported as `fitRecordFactor()`. The coefficients solve $$ \min_{c_i > 0} \; \frac{1}{J}\sum_{j=1}^{J} w_j\, r_j^{2} \;+\; \frac{\lambda_B}{N J}\sum_{i=1}^{N}\sum_{j=1}^{J} v_{ij}^{2}, $$ where $r_j$ is the target-normalized residual of the suite-mean spectrum at period $T_j$, $v_{ij}$ is the excursion of the scaled record spectrum outside the compatibility envelope, $\lambda_B$ (`shape.bandWeight`) weights envelope containment, and $w_j$ applies `shape.deficitWeight` when the suite mean falls below the target. The containment term conditions every coefficient individually, so no artificial coefficient bounds are needed; an admissible range, when configured, acts only as a backstop. ## Prerequisites - A processed record database with a per-process index, e.g. built with `runProcess()`; the match config points at it through `path.database` and `processID`. - A target table (`path.target`, CSV or Rds) with columns `ID`, `TR`, `Vs30`, `p`, `Tn`, and the spectral value named by `target.value`. It must provide the mean spectrum and the envelope quantiles (for example `p05` and `p95`) for the configured `target.ID` and `target.envelope.ID`. ## The site JSON One JSON file per site drives both runners: ```json { "id": "gac30", "path": { "database": "~/kashimaDB/gmdb.v2", "selection": "gmsp/selection/gac30/selection.csv", "target": "oq/data/UHSTable.Rds", "out": "gmsp/match/gac30" }, "processID": "F25", "stage0": { "N": 35 }, "target": { "ID": "GAC30", "TR": 10000, "Vs30": 1500, "p": "mean", "siteID": "gac30", "value": "SaF", "units": "g", "envelope": { "ID": "GAC30.max", "lower": "p05", "upper": "p95" } }, "ocid": "H1", "xi": 0.05, "Tn": { "min": 0.04, "max": 5 }, "shape": { "method": "a3", "bandWeight": 2, "deficitWeight": 4, "modalFactor": { "min": 0.001 } }, "imf": { "includeResidue": false }, "units": { "source": "mm", "target": "mm", "psa": "mm/s2" }, "plot": { "audit": true, "standard": true }, "parallel": true, "workers": 14, "override": true } ``` `path.selection` is written by `runStage0()` and read by `runMatch()`: the selection is a product of the workflow, not a hand-maintained list. ## Step 1: spectral preselection ```{r stage0} runStage0(file = "gmsp/match.json", root = "/path/to/project") ``` Outputs, written next to `path.selection`: - `selection.csv` — the top-`stage0.N` records by band compatibility; - `stage0.spectral.csv` — the full scored pool (band metrics, `Mw`, `Repi`, `PGA`, `AI` per record) for auditing the selection; - `stage0.psa.csv` — the pool spectra cache. The first run sweeps the database (roughly an hour for several thousand records, checkpointed in batches); any later run — a changed screen, a different `N` — reuses the cache and completes in seconds. Sites of the same project share the target period grid, so the cache can be copied between site selection folders to skip the sweep entirely. ### Seismological screening The `candidate` block restricts the pool before ranking: ```json "candidate": { "ownerID": ["NGAW", "CESMD"], "pool": "gmsp/selection/broad/selection.csv", "filter": { "EventMagnitude": { "min": 6, "max": 8 }, "Repi": { "max": 200 } }, "exclude": ["95ef1184f1bea4fa"], "oneRecordPerEvent": true } ``` - `ownerID` — provider whitelist; - `pool` — a broad `runSelect()` output used as the candidate pool, for criteria beyond simple ranges (for example unions of magnitude-distance windows); - `filter` — numeric `min`/`max` ranges over any index column; - `exclude` — record blacklist; - `oneRecordPerEvent` — keeps the best-ranked record of each event. Before fixing `N`, check that enough genuinely compatible records exist: ```{r gate} library(data.table) DT <- fread("gmsp/selection/gac30/stage0.spectral.csv") DT[, .(le05 = sum(bandRMSE <= 0.05, na.rm = TRUE), le10 = sum(bandRMSE <= 0.10, na.rm = TRUE), le20 = sum(bandRMSE <= 0.20, na.rm = TRUE))] ``` If the compatible tail is thinner than the desired suite size, reduce `N` rather than padding the suite with incompatible records. ## Step 2: match ```{r match} runMatch(file = "gmsp/match.json", root = "/path/to/project") ``` Products under `path.out`: ```text data/TSW.csv, PSW.csv, IMW.csv matched histories, spectra, intensities data/ITS.csv, IPSA.csv per-mode traceability products metadata/ selection, indexes, the JSON used match/target.csv the target band actually used match/scale.csv c_i per record match/factors.csv modal shaping and energy factors match/stage.csv record-scaling summary diagnostics match/shaping.csv modal shaping QA (see below) match/metrics.csv, eval.csv per-record optimization diagnostics match/mean.csv suite-mean spectrum match/recordCompatibility.csv per-record band diagnostics match/compatibility.csv suite-level band summary ``` ## Step 3: acceptance diagnostics `match/stage.csv` summarizes the suite solve: - `RMSE` — target-normalized residual of the suite mean; - `pgaRatio` — suite-mean PGA over target PGA; - `insideFraction`, `meanViolation`, `maxViolation` — envelope containment of the scaled record spectra (`maxViolation` should stay below about 1); - `scaleMin/Median/Max`, `nAtMin`, `nAtMax` — the coefficient distribution. With band conditioning the coefficients cluster near their natural scale and `nAtMin`/`nAtMax` are zero unless a backstop range was configured and binds — which warrants investigation, not acceptance. `match/shaping.csv` guards the modal shaping: per record, the minimum and median ratio of the shaped to the source spectrum across the target grid. Ratios far below one at some periods mean the modal shaping removed real spectral content there (for example by excluding the lowest-frequency mode); the selection compatibility paid for by the preselection must survive the shaping. ## The in-memory API The algorithm core of the match runner is exported for direct use on R objects — no files or JSON involved. Three functions cover the two fitting problems and their acceptance diagnostics: ```{r api} library(gmsp) library(data.table) # One record's IMF rows (from TS2IMF), its source TSW (from AT2TS), # and the target band: Tn (including 0 = PGA), Sa.low, Sa.mean, Sa.high # in the PSA units of the signals. Band <- data.table(Tn = c(0, 0.05, 0.10, 0.20), Sa.low = 0.6 * c(100, 50, 35, 20), Sa.mean = c(100, 50, 35, 20), Sa.high = 1.8 * c(100, 50, 35, 20)) # Modal shaping of one record: fits the shape.modalFactor coefficients. Fit <- fitModalFactor(.x = IMF, target = Band, source = TSW, key = Selection[1L], factorMin = 0.001) Fit$coeff # per-mode factors: b, energyFactor, stageFactor Fit$metrics # one-row fit diagnostics (RMSE, band containment, ...) Fit$PSA # source / start / final spectra, long format # Record scaling of the suite: fits one shape.recordFactor per record. PSA <- rbindlist(lapply(Fits, `[[`, "PSA")) # all fitted records Suite <- fitRecordFactor(.x = PSA, target = Band, records = Selection$RecordID, factorMin = 0.5, factorMax = 1.6) Suite$scale # RecordID, scaleFactor (c_i) Suite$stage # one-row solver diagnostics # Compatibility diagnostics of the scaled suite against the band. QA <- matchCompatibility(.x = PSW, target = Band, scale = Suite$scale, factorMin = 0.5, factorMax = 1.6) QA$record # per-record PGA position, containment, RMSE QA$summary # suite aggregate ``` `runMatch()` composes exactly these functions over a project record set, adding file orchestration and the product tree below. Historical note: the two fitting problems appear as "Stage A" and "Stage B" in older logs, and `match/stage.csv` keeps the frozen product label `stage = "B"`. ## Composition with the other runners ```text runProcess() builds the processed database the pool comes from runSelect() optional broad metadata screen feeding candidate.pool runStage0() spectral preselection into path.selection runMatch() modal shaping + record scaling into path.out runPlot() optional raw QA widgets from the products runExport() packages products for delivery ``` Report-grade figures and tables are rendered outside `gmsp` by the reporting layer, from the materialized CSV products listed above.