pipekit-evaluate

Summary¶

Decision. Metrics live with their carriers: xrtoolz.metrics for xarray fields, statistics, events and physical checks; geotoolz for GeoTensor detection building blocks; pipekit-train for differentiable losses. pipekit-evaluate owns what no carrier library can: the taxonomy as tags on scorers, a scorer adapter over any metric operator, the evaluation report, and the benchmark ladder it already ships. There is no pipekit-array.metrics.classification and no third catalog type.

• A. The taxonomy and the report in pipekit-evaluate. Tag every Scorer with its unit, lens, stages and differentiability (A1); adapt any xrtoolz metric operator into a Scorer (A2); grow BenchmarkReport into an EvaluationReport that can be hashed, filtered, diffed and turned into a table (A3); fix the ladder’s rough edges (A4).

• B. The missing metrics, where their names are reserved. FSS and a shift-blind spectrum comparison (B1); POD, FAR, CSI and the contingency table, whose names xrtoolz.metrics.object already reserves (B2); NSE, index of agreement, Brier and the multivariate energy score (B3); Hausdorff distance and structure functions (B4).

• C. The missing wrappers. By season or month (C1), by event (C2) and by member (C3), beside the shipped lead-time, region and band wrappers.

• D. The trajectory unit. Particle advection, dispersion and FTLE in a new xrtoolz.lagrangian (D1), and trajectory scores in the reserved xrtoolz.metrics.lagrangian (D2).

• E. The event unit. Reference detectors (E1) and events as rows of geocatalog with an event profile, as Report 13 did for model states (E2).

• F. Integration. Differentiable scorers as pipekit-train losses (F1), per-step checks in cycles (F2), reports linked to models in the registry (F3), and presets from the four pilots (F4).

End state. A user picks scorers by unit, lens and stage, wraps them by lead, region, scale, season, event or member, and gets one report that records which model, which data and which configuration it scored, and can be compared with the last one. The metrics come from the carrier libraries; pipekit-evaluate composes and records them, and its ladder scores estimators rung by rung against an oracle.

The framework, as the scoping report put it¶

This section restates the scoping report’s framework, which this page adopts unchanged except where a box says otherwise.

Three failure modes. A point-wise loss penalises a correctly shaped but displaced feature twice, once where it is missing and once where it is spurious (the double penalty). Aggregate scores are dominated by large, high-variance scales, so a forecast can lose its small scales at little cost (scale-blind aggregation). And agreement on gridded fields says nothing about how parcels move through them (Eulerian blindness).

Three axes. An evaluation is a choice on each of three axes (Table 1):

A metric’s signature follows its unit: a field metric maps $(\hat{x}, x) \mapsto \mathbb{R}$; an event metric maps a predicted and a reference set of events to a contingency table; a trajectory metric maps a velocity field and drifters to dispersion statistics. Decomposition (by region, scale, lead time, season, event or member) is a wrapper that applies a metric on parts of the data and collects the results; it is not a lens.

Ownership rules¶

A metric lives with the carrier it reads; pipekit-evaluate composes metrics and records what they found (Figure 1).

pipekit-evaluate owns the taxonomy tags, the scorer adapter, the evaluation report and the benchmark ladder. It depends on pipekit and imports carrier libraries only through extras.

Neighbours follow the stack rules. For this page that means the metrics stay in xrtoolz (xarray carriers, with the decomposition wrappers) and geotoolz (GeoTensor detection blocks), neither of which imports pipekit-evaluate, and events become geocatalog rows (workstream E).

What the stack scores today¶

The scoping report’s matrix, filled with the names that ship at the pinned revisions (Table 2), and the decomposition wrappers beside it (Table 3). Every cell not listed below is empty.

The scoping report’s gap analysis, claim by claim, at the pinned revisions:

• “xr_toolz.metrics has good point-wise coverage.” Partly: it covers far more than point-wise. NSE and the index of agreement are missing.

• “Spectral comparison is missing.” Wrong. psd_score, PSDScore, WaveletPSDScore, resolved_scale and FrequencyBandSkill live in xrtoolz.metrics.spectral. What is missing is a comparison of two spectra that ignores phase (stage 1).

• “Structural (SSIM, FSS) is missing.” Partly: SSIM ships; FSS, Hausdorff and perceptual losses do not.

• “Detection is missing.” Partly: instance matching and an IoU-matched F1 ship; event-level POD, FAR and CSI are reserved stubs.

• “xr_toolz.calc has gradients but no balance checks.” Wrong: xrtoolz.metrics.physical has geostrophic-balance, divergence, static-stability and PV-conservation errors.

• “No decomposition wrappers.” Wrong: lead time, region and frequency band ship; season, event and member do not.

• “Trajectories and events are missing entirely.” Partly: trajectories are (xrtoolz.metrics.lagrangian is an empty module marked “V3”), but event building blocks exist across geotoolz, xrtoolz and geocatalog.

• “xr_toolz.budgets exists but is not evaluation-shaped” and “no stage annotation”: hold.

The running example: one methane plume, three ways to get it wrong¶

The stages score the stack’s recurring object, a methane plume, against its truth, using the shipped operators wherever they exist and marking hand-written code as user. The data are synthetic so that every score has a known right answer.

Stage 1 · The double penalty (Field × point-wise, structural, spectral)¶

TL;DR. Score a correctly shaped plume that is displaced, and a flat field that has no plume at all, and see which scores prefer which.

Problem. The truth $x(s)$ is a Gaussian plume of peak 50 ppb and width $\sigma = 3$ cells on a periodic transect of $N = 40$ cells (and a $64 \times 64$ grid in 2-D). Prediction A is the same plume shifted by $d = 6$ cells; prediction B is the flat field at the truth’s mean. For a shift $d$, the point-wise error of A is

which overtakes the flat field’s RMSE once $d \ge 5$. A score that sees the plume should prefer A.

1
2
3
4
5
6
7
8
9
10
11
12
13
# Shapes: 1-D N = 40 cells; 2-D 64 × 64; A = plume shifted 6 cells, B = flat mean
from xrtoolz.metrics import RMSE, SSIM, PhaseShiftError, psd_score

RMSE("ch4", dims="x")(A, truth)            # → 20.50   (B: 15.62)            today
PhaseShiftError("ch4", "x", periodic=True)(A, truth)
# → shift_x = 6, residual_rmse = 0.0       (B: shift 0, residual 15.62)     today
psd_score(A, truth, psd_dims=["x"])        # mean over k > 0: −1.105 (B: 0.0) today
fss(A, truth, threshold=10, window=9)      # → 0.625   (B: 0.0)              user
# relative spectrum difference Σ|P_A − P_x| / Σ P_x over k > 0
spectrum_gap(A, truth)                     # → 0.000   (B: 1.000)            user
# 2-D 64 × 64:
SSIM("ch4", ["y", "x"])(A2, truth2)        # → 0.885   (B: 0.336)            today
RMSE("ch4", dims=["y", "x"])(A2, truth2)   # → 4.671   (B: 4.096)            today

The double penalty is real and the stack can already see past it (Table 4): SSIM and PhaseShiftError both prefer the plume, and PhaseShiftError even reports the 6-cell displacement. Two lenses need code. FSS is not shipped (B1). And the shipped spectral score is not the lens the scoping report meant: psd_score compares the spectrum of the error with the spectrum of the truth, $1 - P_{\hat{x} - x} / P_x$, so it is phase-sensitive and suffers the double penalty too. A comparison of the two spectra, which ignores phase, prefers A perfectly; it is a few lines over the shipped power_spectrum and belongs beside psd_score (B1).

At the shift the scoping report’s examples suggest, $d = 4$, RMSE still prefers the plume (15.44 against 15.62): the flip needs $d \ge 5$ for this width, as (1) predicts.

Stage 2 · Scale-blind aggregation (Statistic × spectral, by scale)¶

TL;DR. Blur a methane enhancement map and see how much the RMSE notices, then look at the spectrum and at the error scale by scale.

Problem. The truth is a $128 \times 128$ field with an energy spectrum $E(k) \propto k^{-3}$, normalised to unit standard deviation. The prediction is the truth blurred by a Gaussian of $\sigma = 3$ pixels. A blur multiplies the spectrum by $e^{-(2 \pi k \sigma)^2}$, so the scale at which the error’s power reaches half the truth’s, the resolved scale, is known analytically: 12.03 pixels.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# Shapes: 128 × 128 pixels; k in cycles per pixel
from xrtoolz.metrics import FrequencyBandSkill, RMSE, SSIM, psd_score, resolved_scale
from xrtoolz.transforms import fit_spectral_slope, isotropic_power_spectrum

RMSE("z", dims=["y", "x"])(blur, truth)     # → 0.098 of the truth's std     today
SSIM("z", ["y", "x"])(blur, truth)          # → 0.885                         today
# slope of E(k) over k ∈ [0.05, 0.15] cycles / px
fit_spectral_slope(isotropic_power_spectrum(truth.z, ["y", "x"]),
                   k_min=0.05, k_max=0.15)  # → −3.10  (blur: −9.39)          today
score = psd_score(blur.z, truth.z, psd_dims=["y", "x"], isotropic=True)
resolved_scale(score, "freq_r", level=0.5)  # frequency coordinate, 50 % level
# → 11.98 px  (analytic 12.03 px)                                            today
bands = {"large": (0, 0.05), "meso": (0.05, 0.125), "small": (0.125, 0.5)}
FrequencyBandSkill("z", ["y", "x"], bands)(blur, truth)
# → RMSE 0.053, 0.067, 0.049; divided by each band's truth std: 0.05, 0.60, 0.99

The RMSE is under 10 % of the signal and SSIM is 0.885, yet the prediction has lost almost all of its variance below 8 pixels and its spectral slope is −9.4 where the truth’s is −3.1. The shipped tools see this: resolved_scale recovers the analytic 12-pixel scale to 0.4 %, and FrequencyBandSkill shows the small-scale error at 99 % of that band’s signal. The scoping report’s spectral lens and “by scale” wrapper exist; the one missing piece is a score that compares the two slopes directly (B1).

Stage 3 · Did it find the plumes? (Event × detection)¶

TL;DR. Detect plumes in a predicted and a reference map, match them, and count hits, misses and false alarms.

Problem. The reference map has two plumes. The prediction has the first one shifted by 3 pixels, misses the second, and adds one false alarm. With $a$ hits, $b$ false alarms and $c$ misses, the detection scores are

1
2
3
4
5
6
7
8
9
10
11
12
13
# Shapes: two label images; masks (n_objects, H, W)
from xrtoolz.metrics import instance_f1_at_iou, mask_iou_matrix, match_instances
labels_ref, labels_pred = label(ref > 0.3), label(pred > 0.3)    # scipy   user
mask_iou_matrix(pred_masks, ref_masks)       # → [[0.438, 0], [0, 0]]       today
match_instances(pred_masks, ref_masks, iou_threshold=0.3)
# → tp = 1, fp = 1, fn = 1                                                  today
instance_f1_at_iou(pred_masks, ref_masks, iou_threshold=0.3)  # → F1 = 0.5   today
# POD = a / (a + c) = 0.5;  FAR = b / (a + b) = 0.5               user
# CSI = a / (a + b + c) = 0.333                                  user
CriticalSuccessIndex()
# → NotImplementedError: … a name reservation; implementation lands with V5
centroid_displacement(pred_objects, ref_objects)
# → pairs object 2 with object 2: 30.6 px between unrelated plumes

Matching works: match_instances counts the hit, the miss and the false alarm correctly, and the scores in (2) follow from its counts in one line each. They are not shipped, although their names are reserved in xrtoolz.metrics.object (B2). CentroidDisplacement pairs objects by label id rather than by matching, so it compares unrelated plumes when the labels are not aligned; it should take match_instances’ pairing (B2). Detection itself is assembled from geotoolz (segment.Threshold, measure.LabelConnectedComponents, plume.PlumeMask) and xrtoolz (remove_small_objects_2d); what does not exist is an event as an object with an id, a time span and attributes (E).

Stage 4 · The benchmark ladder (pipekit-evaluate)¶

TL;DR. Score three estimators of the same linear-Gaussian retrieval with the ladder pipekit-evaluate ships: an exact oracle, a prior-only baseline and an over-confident shortcut.

Problem. The state $x \in \mathbb{R}^{40}$ is drawn from $\mathcal{N}(0, B)$ with a squared-exponential covariance of length 3 cells; every fourth cell is observed, $y = H x + \varepsilon$ with $\varepsilon \sim \mathcal{N}(0, 0.3^2 I)$. The oracle is the exact posterior mean $m = B H^\top (H B H^\top + R)^{-1} y$ and covariance. The ladder places each estimator in a cell (level, rung, complexity) and scores it against the reference its rung calls for: the known truth for model-based inference, the oracle’s posterior for an amortized predictor.

1
2
3
4
5
6
7
8
9
10
11
12
13
# Shapes: x (40,); y (10,)
import pipekit_evaluate as pev

report = pev.Benchmark(
    task=LGTask(),                                   # 40 cells, every 4th  user
    estimators=[ExactOI(), ZeroPrior(), WrongRidge()],                    # user
    scorers=[RMSE(), GaussCoverage()],               # pev.Scorer protocol  user
    level="L4",
).run()                                                                 # today
# L4/r2/c1  model_based_inference  ref=known_theta       rmse 0.406  cov90 0.825
# L4/r2/c0  model_based_inference  ref=known_theta       rmse 0.569  cov90 0.975
# L4/r5/c0  amortized_predictor    ref=oracle_posterior  rmse 0.597
report.to_dict()           # {cell label: {metric: value}}; drops reference and notes

The ladder runs and resolves references as its design says: the exact oracle beats the prior-only baseline against the truth, and the over-confident shortcut is scored against the oracle’s posterior, not the truth. It took about 80 lines of glue, and the run exposed six frictions, all in Section ??:

• No carriers ship. The Gaussian prior, the estimate, the observation operator and a no-op forward model for a static problem are all user code, and Benchmark never calls the task’s prior(), simulator() or observation_operator().

• One reference per cell. The amortized rung cannot be scored against both the oracle and the truth in one run.

• Estimate has no covariance. A coverage scorer must read a covariance by duck typing or sample.

• The complexity axis scores the oracle against a worse baseline (RMSE 0.735 of the exact posterior against the prior), and its note for a step-0 cell says “generative rung”, which a step-0 cell is not.

• The docstring and the code disagree on which estimator becomes the oracle (run.py:77 says the lowest complexity; the code takes the highest).

• to_dict drops the reference, the note and the axis, and nothing hashes, tables or logs the report.

Report 17 runs the same ladder on rungs 2 to 5 of one retrieval, including the emulator rungs 3 and 4 this stage leaves empty.

Stage 5 · What cannot be scored yet: trajectories¶

TL;DR. The Lagrangian unit has no producer in the stack.

Problem. For a velocity field $u(s, t)$, a particle path solves $\dot{s} = u(s, t)$; dispersion is the mean squared separation of initially close pairs, and the finite-time Lyapunov exponent is the growth rate of that separation. Scoring a forecast’s transport skill needs particles advected through both the forecast and the reference velocity.

Nothing in the pinned stack advects a particle. xrtoolz.metrics.lagrangian is an empty module reserved for “V3”; PVConservationError takes trajectories as input but nothing produces them; and viz.validation.EulerianLagrangianPanel plots them. Workstream D adds the producer in xrtoolz and the scores beside the reserved module. No run here.

The scoping report’s names, checked¶

Table 5 checks the stack names the scoping report’s examples use.

Current-state audit¶

Of fourteen pieces checked, eight ship and work, three ship with defects or gaps, and three are missing (Table 6); Figure 2 counts the page’s calls by marker. The audit read the pinned revisions on 13 October 2026.

Smaller findings:

• No metric is differentiable or stage-tagged. xrtoolz’s kernels are NumPy (_pixel_kernels.py: “Backend: numpy. JAX / CuPy variants are out of scope”), and none carries a flag; differentiable losses live in pipekit_train.loss.

• pot_exceedances returns exceedances, not events: no ids or durations, so marine-heatwave-style events need a detector (E1).

• The pilots name libraries this audit did not check (somax, vardax, plumax, xtremax, filterax, gaussx); their scorer lists (SBC, coverage, twCRPS, C2ST, Okubo–Weiss detection, triple collocation) are the best catalogue of what the ladder will need.

• pipekit-evaluate’s extras are cycle, xskill and proper; the xskill and proper extras have no code behind them yet.

Workstream A · The taxonomy and the report in pipekit-evaluate¶

Decision. pipekit-evaluate’s Scorer carries the taxonomy as data, any metric operator becomes a scorer through one adapter, and the ladder’s report becomes the scoping report’s EvaluationReport (Table 7).

1
2
3
4
5
6
7
8
9
10
11
12
# Shapes: as stage 1; xm is xrtoolz.metrics
scorers = [
    pev.MetricScorer(xm.RMSE("ch4", dims="x"), unit="field", lens="pointwise",
                     stages={"validation", "final"}),          # proposed (A2)
    pev.MetricScorer(xm.SSIM("ch4", ["y", "x"]), unit="field", lens="structural",
                     stages={"final"}),                         # proposed (A2)
    pev.MetricScorer(xm.FSS("ch4", threshold=10, window=9), unit="field",
                     lens="structural", stages={"final"}),      # proposed (A2, B1)
]
report = pev.evaluate(prediction, truth, scorers, model_ref=h)  # proposed (A3)
report.filter(lens="structural").to_pandas()                   # proposed (A3)
report.diff(last_week)                  # per-entry changes      proposed (A3)

Workstream B · The missing metrics¶

Decision. Each missing metric goes into the xrtoolz module that already holds its lens, under the name it reserves where one is reserved (Table 8).

Workstream C · The missing wrappers¶

Decision. Three more wrappers join SkillByLeadTime, EvaluateByRegion and FrequencyBandSkill, with the same shape: a metric in, a dataset with one more dimension out (Table 9).

Workstream D · The trajectory unit¶

Decision. Particle advection is a transform over velocity fields and lives in xrtoolz; trajectory scores live in the module xrtoolz reserves for them (Table 10).

Workstream E · The event unit¶

Decision. Events are detected with the stack’s segmentation tools and stored as catalog rows with an event profile, not in a third catalog type (Table 11).

Workstream F · Integration¶

Decision. Evaluation connects to training, cycles and the registry through the seams those packages already have or that their own reports propose (Table 12).

Testing gates¶

Each workstream merges only when its gates pass; the existing gates are checks the code passes today (Table 13).

Sequencing and risks¶

The taxonomy and report come first, because they need no new metrics and make every later metric land in a structured place; metrics fill the matrix next, then the two new units, then the integrations (Figure 3).

• Two homes for metrics. Tags in pipekit-evaluate and code in xrtoolz could drift. Fallback: the tags live on the MetricScorer wrapper, not on xrtoolz operators, so xrtoolz never depends on the taxonomy.

• Differentiable lenses do not exist yet. Every xrtoolz metric is NumPy, so the Stage axis’s “training” column is empty except for pipekit-train’s losses. Fallback: F1 accepts only pipekit_train.loss objects and JAX scorers until xrtoolz grows JAX kernels.

• Spectral scores are sensitive to windows and detrending. Stage 1 and 2 passed window=None and detrend=None; defaults change the numbers. Fallback: every spectral scorer records its spectral options in its report entry.

• Event detection is domain-specific. Plumes, heatwaves and eddies need different detectors. Fallback: ship two reference detectors (threshold plus connected components, and the heatwave rule) and leave the rest to domain libraries.

• Lagrangian evaluation is expensive. Advecting many particles over long windows is final-evaluation-only. Fallback: tag trajectory scorers stages={"final"} and provide subsampled variants.

Open questions¶

1. Does the taxonomy belong on the metric or on the scorer? This page puts it on the scorer (A1). If xrtoolz operators carried it too, a catalogue of the whole matrix could be generated from the code.

2. One report type or two? BenchmarkReport (rung by rung) and the scoping report’s EvaluationReport (lens by lens) could be one class with two views, or two classes.

3. Should events live in geocatalog? The event profile follows Report 13’s decision for states. Events also carry lineage to the field they were detected in; is a source_uri column enough?

4. Where do differentiable metrics come from? JAX versions of xrtoolz kernels, a pipekit-jax metric set, or pipekit-train losses only?

5. What should the capability index cover? Indexing every xrtoolz subpackage adds hundreds of names. All of them, or the public namespaces only (metrics, transforms, geo, ocn, budgets)?

Upcoming work¶

Table 14 lists the follow-ups.

References¶

Sources¶

• The scoping report, “Report 14 — pipekit-evaluate: multidimensional evaluation as first-class infrastructure” (supplied by the author).

• pipekit-evaluate 0.0.2 in references/geostack/pipekit/packages/pipekit-evaluate: its source, tests, README, CHANGELOG, benchmark_ladder.md and the four pilot designs.

• xrtoolz 0.0.5 (xrtoolz.metrics, .transforms, .geo, .budgets), geotoolz and geocatalog at the pinned revisions.

• The running example’s scripts and outputs, described in Section ??.