---
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.
