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, and pipekit-train for differentiable losses. pipekit-evaluate owns what no carrier library can, which is 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. We 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 which can be hashed, filtered, diffed and turned into a table (A3), and 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, the index of agreement, Brier and the multivariate energy score (B3), and the 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 which records the model, the data and the configuration it scored, and which 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 we adopt unchanged on this page 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 very little cost (scale-blind aggregation). Lastly, 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 a prediction and a reference, $(\hat{\vect{x}}, \vect{x}) \mapsto \mathbb{R}$ with both in $\mathbb{R}^{D_x}$, whereas an event metric maps a predicted and a reference set of events to a contingency table, and 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 which applies a metric on parts of the data and collects the results, so it is not a lens.

Ownership rules¶

A metric lives with the carrier it reads, and 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 only imports carrier libraries through extras.

The 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¶

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

We go through 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, although 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 which ignores phase (stage 1).

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

• “Detection is missing.” Partly: instance matching and an IoU-matched F1 ship, but 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 missing (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. We use the shipped operators wherever they exist and mark hand-written code as user. The data are synthetic, so 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) \in \mathbb{R}$ is a Gaussian plume of peak 50 ppb and width $\sigma_\text{plume} = 3$ cells on a periodic transect of $N_{\Omega_x} = 40$ cells, with $s \in \Omega$ the position along the transect (and a $64 \times 64$ grid in 2-D). Prediction A, $\hat{x}_A$, is the same plume shifted by $d = 6$ cells, and prediction B, $\hat{x}_B$, is the flat field at the truth’s mean. For a shift $d$, the point-wise error of A is

where $\langle \cdot \rangle$ is the mean over the transect, and it overtakes the flat field’s RMSE once $d \ge 5$. So 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 κ > 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 κ > 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. However, two lenses need code. FSS is not shipped (B1). In addition, the shipped spectral score is not the lens which the scoping report meant. psd_score compares the spectrum of the error with the spectrum of the truth, $1 - P_{\hat{x} - x}(\kappa) / P_x(\kappa)$ with $P_x(\kappa) \in \mathbb{R}_{\ge 0}$ the power spectral density at wavenumber $\kappa$, so it is phase-sensitive and suffers the double penalty too. A comparison of the two spectra, which ignores phase, prefers A perfectly, and it is a few lines over the shipped power_spectrum which belong beside psd_score (B1).

At the shift which the scoping report’s examples suggest, $d = 4$, RMSE still prefers the plume (15.44 against 15.62), because 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(\kappa) \propto \kappa^{-3}$, normalised to unit standard deviation, where $\kappa \in \mathbb{R}_{\ge 0}$ is the radial wavenumber in cycles per pixel. The prediction is the truth blurred by a Gaussian of width $\sigma_\text{blur} = 3$ pixels. A blur multiplies the spectrum by $e^{-(2 \pi \kappa \sigma_\text{blur})^2}$, so the scale at which the error’s power reaches half the truth’s (the resolved scale) is known analytically, at 12.03 pixels.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# Shapes: 128 × 128 pixels; κ 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(κ) over κ ∈ [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. So the scoping report’s spectral lens and “by scale” wrapper exist, and 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 $N_\text{hit}$ hits, $N_\text{false}$ false alarms and $N_\text{miss}$ misses (all counts), the detection scores are (2), (3) and (4),

where each score is a scalar in $[0, 1]$.

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 = N_hit / (N_hit + N_miss) = 0.5;  FAR = N_false / (N_hit + N_false) = 0.5   user
# CSI = N_hit / (N_hit + N_false + N_miss) = 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, since match_instances counts the hit, the miss and the false alarm correctly, and the scores in (2) to (4) follow from its counts in one line each. However, 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, and it should take the pairing of match_instances instead (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 $\vect{x} \in \mathbb{R}^{D_x}$, with $D_x = 40$ cells, is drawn from $\mathcal{N}(\vect{0}, \gmat{\Sigma}_x)$ with a squared-exponential covariance $\gmat{\Sigma}_x \in \mathbb{R}^{D_x \times D_x}$ of length 3 cells. Every fourth cell is observed, so $\vect{y} = \mat{H} \vect{x} + \gvect{\varepsilon}$ with $\vect{y} \in \mathbb{R}^{D_y}$, $D_y = 10$, the selection matrix $\mat{H} \in \mathbb{R}^{D_y \times D_x}$ and $\gvect{\varepsilon} \sim \mathcal{N}(\vect{0}, \gmat{\Sigma}_\varepsilon)$ with $\gmat{\Sigma}_\varepsilon = 0.3^2 \mat{I}_{D_y}$. The oracle is the exact posterior mean $\gvect{\mu}_{x \mid y} = \gmat{\Sigma}_x \mat{H}^\top (\mat{H} \gmat{\Sigma}_x \mat{H}^\top + \gmat{\Sigma}_\varepsilon)^{-1} \vect{y} \in \mathbb{R}^{D_x}$ and its covariance $\gmat{\Sigma}_{x \mid y} \in \mathbb{R}^{D_x \times D_x}$. The ladder places each estimator in a cell (level, rung, complexity) and scores it against the reference which its rung calls for, namely the known truth for model-based inference and 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 and not the truth. However, it took about 80 lines of glue, and the run exposed six frictions, which are 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. So 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 which 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 $\vect{u}(\vect{s}, t) \in \mathbb{R}^{D_s}$ at location $\vect{s} \in \Omega \subset \mathbb{R}^{D_s}$ and time $t \in \mathcal{T}$, a particle path $\vect{s}(t)$ solves $\mathrm{d}\vect{s}(t) / \mathrm{d}t = \vect{u}(\vect{s}(t), t)$. Dispersion is the mean squared separation of initially close pairs, and the finite-time Lyapunov exponent is the growth rate of that separation. So scoring a forecast’s transport skill needs particles which are 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. There is 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 the fourteen pieces we checked, eight ship and work, three ship with defects or gaps, and three are missing (Table 6), and 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”), none of them carries a flag, and the differentiable losses live in pipekit_train.loss.

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

• The pilots name libraries which 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, and 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 which 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, and they have the same shape: a metric goes in, and a dataset with one more dimension comes out (Table 9).

Workstream D · The trajectory unit¶

Decision. Particle advection is a transform over velocity fields, so it lives in xrtoolz, and trajectory scores live in the module which 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, and not in a third catalog type (Table 11).

Workstream F · Integration¶

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

Testing gates¶

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

Sequencing and risks¶

The taxonomy and report come first, because they need no new metrics and they make every later metric land in a structured place. The metrics fill the matrix next, then the two new units, and 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 and 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 only accepts pipekit_train.loss objects and JAX scorers until xrtoolz grows JAX kernels.

• Spectral scores are sensitive to windows and detrending. Stages 1 and 2 passed window=None and detrend=None, and the 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: we 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 for final evaluation only. Fallback: we 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). However, 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, so 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. Should it cover all of them, or only the public namespaces (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 ??.