The QuakeScope picker benchmarks, consolidated

Five notebooks benchmark the pickers behind the 2026 catalogue, each against a different reference and each rendered as its own report. This notebook reads their exported result tables (docs/benchmark/results/, written by an export cell at the end of each notebook when it is executed) and draws the figures a paper needs from all of them at once. Nothing is re-picked here and no number is typed in: change a benchmark, re-execute it, re-execute this.

study	notebook	reference	result files
US sequences	phasenet_sequence_comparison.ipynb	reviewed ComCat events, arrivals from SCEDC + NCEDC, manual only	us_sequences/
Global sequences	phasenet_global_sequences.ipynb	GeoNet, INGV, NOA bulletins, manual only	global_sequences/
Ridgecrest aftershocks	phasenet_aftershock_benchmark.ipynb	SCEDC analyst picks, one dense 30 min window	ridgecrest_aftershocks/
Ocean bottom	phasenet_obs_offshore_benchmark.ipynb	iasp91 predictions; Barcheck AACSE analyst picks; UW Axial picks; LDEO ML-DD events	obs_offshore/
Western reproduction	western_pick_validation.ipynb	the campaign's own Parquet, re-picked via FDSN	western_reproduction/

The protocol, the caveats and the reading of these figures are in docs/benchmark/README.md.

import json
from pathlib import Path

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

RES = Path("../docs/benchmark/results")
WEIGHTS = ["quakescope2026", "jma_wc", "original", "instance"]
COLORS = dict(zip(WEIGHTS, ["#2a78d6", "#eb6834", "#1baf7a", "#eda100"]))
C_P, C_S = "#2a78d6", "#eb6834"
plt.rcParams.update({"figure.dpi": 110, "axes.grid": True, "grid.alpha": 0.25, "grid.linewidth": 0.5,
                     "axes.axisbelow": True, "axes.spines.top": False, "axes.spines.right": False})

def load(study, name):
    p = RES / study / f"{name}.csv"
    return pd.read_csv(p) if p.exists() else None

meta = {s: json.loads((RES / s / "meta.json").read_text()) for s in
        ["us_sequences", "global_sequences", "ridgecrest_aftershocks", "obs_offshore", "western_reproduction"]
        if (RES / s / "meta.json").exists()}
print("result sets and when they were executed:")
for s, m in meta.items():
    print(f"  {s:24s} {m['executed']}  seisbench {m.get('seisbench')}")

result sets and when they were executed:
  us_sequences             2026-09-16T19:41:46+00:00  seisbench 0.12.3
  global_sequences         2026-09-16T19:31:15+00:00  seisbench 0.12.3
  ridgecrest_aftershocks   2026-09-16T19:21:05+00:00  seisbench 0.12.3
  obs_offshore             2026-09-16T19:34:01+00:00  seisbench 0.12.3
  western_reproduction     2026-09-16T19:46:12+00:00  seisbench 0.12.3

1. Recall against analyst picks at the shared 0.3 threshold

Seven sequences with a published analyst reference: four in the western United States (SCEDC and NCEDC, manually reviewed picks on curated event lists) and three abroad (the operators' bulletins). Six stations per sequence abroad, five in the US; windows of 30 min to 3 h starting 10 min after the mainshock; a pick counts as recovered within 0.5 s. The number under each group is the analyst pick count, the sample size of that recall.

us = load("us_sequences", "recall_at_threshold")
gl = load("global_sequences", "recall_at_threshold")
us["region"], gl["region"] = "US", "abroad"
bench = pd.concat([us, gl], ignore_index=True)
ORDER = [s for s in ["Ridgecrest", "San Simeon", "Monte Cristo", "Mendocino 2024", "Kaikoura 2016", "Norcia 2016", "Thessaly 2021"]
         if s in set(bench.sequence)]

fig, axes = plt.subplots(2, 1, figsize=(12.5, 7.2), sharex=True)
for ax, phase in zip(axes, ("P", "S")):
    sub = bench[bench.phase == phase]
    w = 0.8 / len(WEIGHTS)
    for i, name in enumerate(WEIGHTS):
        xs, ys = [], []
        for j, s in enumerate(ORDER):
            r = sub[(sub.sequence == s) & (sub.weights == name)]
            if len(r):
                xs.append(j + (i - 1.5) * w); ys.append(float(r.recall.iloc[0]))
        ax.bar(xs, ys, width=w * 0.92, color=COLORS[name], label=name)
        for x, y in zip(xs, ys):
            ax.text(x, y + 0.012, f"{y:.2f}", ha="center", fontsize=6.6, color="#3d3d3d")
    ax.set_ylim(0, 1.08); ax.set_ylabel(f"{phase} recall at conf >= 0.3"); ax.grid(axis="x")
axes[0].legend(frameon=False, ncol=4, fontsize=9, loc="upper right")
n_of = lambda s, ph: int(bench[(bench.sequence == s) & (bench.phase == ph)].analyst.iloc[0])
axes[1].set_xticks(range(len(ORDER)))
axes[1].set_xticklabels([f"{s}\nn = {n_of(s, 'P')} P, {n_of(s, 'S')} S" for s in ORDER], fontsize=8.5)
axes[1].axvline(3.5, color="#8a8a8a", lw=0.8, ls=":"); axes[0].axvline(3.5, color="#8a8a8a", lw=0.8, ls=":")
axes[0].set_title("Recall against analyst picks at a shared threshold (US left, abroad right)", loc="left", fontsize=11)
fig.tight_layout()

tab = (bench.pivot_table(index=["phase", "sequence"], columns="weights", values="recall")
            .reindex(columns=WEIGHTS).reindex([(p, s) for p in ("P", "S") for s in ORDER]))
tab["analyst"] = [int(bench[(bench.phase == p) & (bench.sequence == s)].analyst.iloc[0]) for p, s in tab.index]
tab.round(3)

	weights	quakescope2026	jma_wc	original	instance	analyst
phase	sequence
P	Ridgecrest	0.605	0.582	0.648	0.216	347
	San Simeon	0.823	0.855	0.855	0.887	62
	Monte Cristo	0.750	0.812	0.812	0.875	16
	Mendocino 2024	0.830	0.900	0.680	0.840	200
	Kaikoura 2016	0.661	0.703	0.658	0.644	357
	Norcia 2016	0.763	0.795	0.736	0.618	701
	Thessaly 2021	0.719	0.747	0.724	0.724	395
S	Ridgecrest	0.517	0.487	0.725	0.181	298
	San Simeon	0.750	0.750	0.938	0.750	16
	Monte Cristo	0.417	0.500	0.667	0.500	12
	Mendocino 2024	0.619	0.681	0.735	0.558	113
	Kaikoura 2016	0.502	0.520	0.705	0.533	458
	Norcia 2016	0.744	0.777	0.907	0.669	656
	Thessaly 2021	0.390	0.402	0.429	0.393	331

2. A threshold is not an operating point

The models' probabilities sit on different scales, so at 0.3 they emit different numbers of picks and the ranking above partly measures liberality. Recall against picks emitted takes the threshold out: each curve is the same model swept from 0.02 to 0.7, and the comparison is made at equal pick counts.

sweep = pd.concat([load("us_sequences", "threshold_sweep").assign(region="US"),
                   load("global_sequences", "threshold_sweep").assign(region="abroad")], ignore_index=True)

seqs = [s for s in ORDER if s in set(sweep.sequence)]
fig, axes = plt.subplots(2, len(seqs), figsize=(3.3 * len(seqs), 6.4), squeeze=False)
for row, phase in enumerate(("P", "S")):
    for ax, s in zip(axes[row], seqs):
        sub = sweep[(sweep.phase == phase) & (sweep.sequence == s)]
        for name in WEIGHTS:
            d = sub[sub.weights == name].sort_values("emitted")
            if not len(d):
                continue
            ax.plot(d.emitted, d.recall, marker="o", ms=3, lw=1.6, color=COLORS[name], label=name)
            star = d[np.isclose(d.thr, 0.3)]
            if len(star):
                ax.plot(star.emitted, star.recall, marker="*", ms=11, color=COLORS[name], mec="#16150f", mew=0.5, zorder=5)
        ax.set_title(f"{s}  {phase}", fontsize=9.5, loc="left"); ax.set_ylim(0, 1)
        if not len(sub):
            n = int(bench[(bench.sequence == s) & (bench.phase == phase)].analyst.iloc[0])
            ax.text(0.5, 0.5, f"reference too small\nto sweep (n = {n})", ha="center", va="center",
                    transform=ax.transAxes, fontsize=8.5, color="#7a7973")
            ax.set_xticks([])
        if row == 1:
            ax.set_xlabel("picks emitted", fontsize=9)
    axes[row][0].set_ylabel(f"{phase} recall")
axes[0][0].legend(frameon=False, fontsize=7.5)
fig.suptitle("Recall against picks emitted; the star is the shared 0.3 threshold", x=0.01, ha="left", fontsize=11)
fig.tight_layout()

Recall at a matched pick budget

For each sequence and phase, the budget is the midpoint of the range every model can reach, and each model's recall is read off its curve there. This is the table the choice of weights should be made from.

def matched(sweep, phase, sequence, names=WEIGHTS):
    sub = sweep[(sweep.phase == phase) & (sweep.sequence == sequence)]
    present = [n for n in names if n in set(sub.weights)]
    if len(present) < 2:                    # the sweep skips references under 20 picks
        return None
    lo = max(sub[sub.weights == n].emitted.min() for n in present)
    hi = min(sub[sub.weights == n].emitted.max() for n in present)
    if not np.isfinite([lo, hi]).all() or hi <= lo:
        return None
    mid = 0.5 * (lo + hi)
    row = {"budget": int(round(mid))}
    for n in present:
        d = sub[sub.weights == n].sort_values("emitted")
        row[n] = float(np.interp(mid, d.emitted, d.recall))
    return row

rows = []
for phase in ("P", "S"):
    for s in seqs:
        r = matched(sweep, phase, s)
        if r is None:                       # a model tops out below the others' floor
            r = matched(sweep, phase, s, [n for n in WEIGHTS if n != "instance"])
            if r is None:
                continue
            r["instance"] = np.nan
        rows.append(dict(phase=phase, sequence=s, **r))
budget = pd.DataFrame(rows).set_index(["phase", "sequence"])
budget.round(3)

		budget	quakescope2026	jma_wc	original	instance
phase	sequence
P	Ridgecrest	462	0.689	0.671	0.631	NaN
	San Simeon	395	0.830	0.817	0.857	0.903
	Mendocino 2024	490	0.832	0.814	0.690	0.888
	Kaikoura 2016	1214	0.560	0.587	0.588	0.698
	Norcia 2016	910	0.660	0.680	0.621	0.773
	Thessaly 2021	1333	0.662	0.688	0.653	0.756
S	Ridgecrest	411	0.654	0.635	0.624	NaN
	Mendocino 2024	320	0.580	0.579	0.616	0.650
	Kaikoura 2016	902	0.452	0.452	0.423	0.667
	Norcia 2016	936	0.651	0.669	0.612	0.849
	Thessaly 2021	1063	0.401	0.403	0.359	0.445

Blank instance cells are sequences where it cannot reach the others' pick counts even with its threshold on the floor: a ceiling, not a calibration offset, and the reason it is not the default weight despite winning most of the matched comparisons where it does reach.

pair = budget.dropna(subset=["quakescope2026", "jma_wc"]).copy()
pair["v7 minus jma_wc, matched budget"] = pair.quakescope2026 - pair.jma_wc
at_thr = bench.pivot_table(index=["phase", "sequence"], columns="weights", values="recall")
pair["v7 minus jma_wc, at 0.3"] = (at_thr.quakescope2026 - at_thr.jma_wc).reindex(pair.index)

fig, ax = plt.subplots(figsize=(8.5, 4.2))
y = np.arange(len(pair))
ax.axvline(0, color="#8a8a8a", lw=1)
ax.plot(pair["v7 minus jma_wc, at 0.3"], y, marker="o", ls="none", ms=6, color="#b8b8b3", label="at the shared 0.3")
ax.plot(pair["v7 minus jma_wc, matched budget"], y, marker="o", ls="none", ms=7, color="#2a78d6", label="at a matched pick budget")
ax.set_yticks(y); ax.set_yticklabels([f"{s}  {p}" for p, s in pair.index], fontsize=8.5)
ax.set_xlabel("recall of quakescope2026 minus recall of jma_wc"); ax.grid(axis="y")
ax.legend(frameon=False, fontsize=9, loc="lower left")
ax.set_title("The fine-tune against its parent, every sequence and phase", loc="left", fontsize=11)
fig.tight_layout()
print(f"mean difference at matched budget: P {pair.xs('P')['v7 minus jma_wc, matched budget'].mean():+.3f}, "
      f"S {pair.xs('S')['v7 minus jma_wc, matched budget'].mean():+.3f}; "
      f"positive in {(pair['v7 minus jma_wc, matched budget'] > 0).sum()} of {len(pair)} cases")

mean difference at matched budget: P -0.004, S -0.000; positive in 5 of 11 cases

3. Timing

Mean absolute residual of matched picks against the analyst time at the shared threshold. The fine-tune was selected on timing; this is where an edge would show.

mae = bench.pivot_table(index=["phase", "sequence"], columns="weights", values="MAE").reindex(columns=WEIGHTS)
mae = mae.reindex([(p, s) for p in ("P", "S") for s in ORDER])
(mae * 1000).round(0).astype("Int64").rename(columns=lambda c: f"{c} (ms)")

	weights	quakescope2026 (ms)	jma_wc (ms)	original (ms)	instance (ms)
phase	sequence
P	Ridgecrest	40	33	58	48
	San Simeon	67	66	58	83
	Monte Cristo	47	51	49	89
	Mendocino 2024	71	74	119	89
	Kaikoura 2016	139	139	154	141
	Norcia 2016	49	48	57	56
	Thessaly 2021	129	135	139	139
S	Ridgecrest	53	64	66	79
	San Simeon	20	35	103	94
	Monte Cristo	219	236	146	174
	Mendocino 2024	86	101	124	120
	Kaikoura 2016	162	163	177	166
	Norcia 2016	64	69	85	76
	Thessaly 2021	129	132	146	148

4. Dense aftershocks, one station cluster

A separate, smaller benchmark on a 30 min Ridgecrest window with overlapping events, scored the same way; its tolerance sweep shows the recall numbers are not sensitive to the matching window between 0.25 and 2 s.

ra = load("ridgecrest_aftershocks", "recall_at_threshold")
sens = load("ridgecrest_aftershocks", "tolerance_sensitivity")
display(ra.pivot(index="weights", columns="phase", values=["analyst", "recall", "MAE"]).reindex(WEIGHTS).round(3))
sens.set_index("tolerance_s").round(3)

	analyst		recall		MAE
phase	P	S	P	S	P	S
weights
quakescope2026	125.0	144.0	0.928	0.792	0.034	0.045
jma_wc	125.0	144.0	0.936	0.806	0.030	0.053
original	125.0	144.0	0.872	0.868	0.037	0.056
instance	125.0	144.0	0.784	0.729	0.036	0.080

	quakescope2026	jma_wc	original	instance
tolerance_s
0.25	0.792	0.799	0.868	0.715
0.50	0.792	0.806	0.868	0.729
1.00	0.792	0.806	0.868	0.729
2.00	0.792	0.806	0.868	0.729

5. Ocean-bottom stations

Three questions, three references. Do OBS-trained weights beat land weights on OBS data (iasp91-predicted P within 10 s, three deployments)? Do the obs campaign's stored picks agree with published analyst picks (Barcheck's AACSE tables, 1 s)? And what happens on caldera microseismicity (Axial Seamount, against the UW real-time picks and the LDEO ML-DD events)?

det = load("obs_offshore", "detection_vs_iasp91")
aacse = load("obs_offshore", "aacse_campaign_vs_analyst")
resc = pd.read_csv(RES / "obs_offshore" / "aacse_rescoring.csv", header=[0, 1])
resc.columns = [a if b.startswith("Unnamed") or b == "P_median_dt" and a == "median_P_dt_s" else f"{a} {b}" for a, b in resc.columns]
resc = resc.rename(columns={resc.columns[0]: "weights"}).set_index("weights")
axial = load("obs_offshore", "axial_event_detection_by_magnitude")

fig, axes = plt.subplots(1, 3, figsize=(15, 4), gridspec_kw=dict(width_ratios=[1.2, 1, 1]))

ax = axes[0]
piv = det.pivot(index="weights", columns="experiment", values="rate")
piv = piv.reindex(["pickblue_phasenet", "pickblue_eqt", "obstransformer", "quakescope2026", "original"])
x = np.arange(len(piv)); w = 0.26
for k, (exp, c) in enumerate(zip(piv.columns, ["#2a78d6", "#eb6834", "#1baf7a"])):
    ax.bar(x + (k - 1) * w, piv[exp], width=w * 0.92, color=c, label=exp)
ax.set_xticks(x); ax.set_xticklabels(piv.index, rotation=20, ha="right", fontsize=8)
ax.set_ylim(0, 1); ax.set_ylabel("P within 10 s of iasp91, conf >= 0.3"); ax.legend(frameon=False, fontsize=8)
ax.set_title("Five models, three deployments", loc="left", fontsize=10.5); ax.grid(axis="x")

ax = axes[1]
tols = [0.25, 0.5, 1.0, 2.0]
for _, r in aacse.iterrows():
    ax.plot(tols, [r[f"recall_{t:g}s"] for t in tols], marker="o", ms=4, lw=1.6,
            color=C_P if r.pha == "P" else C_S, ls="-" if r.stations == "OBS" else "--",
            label=f"{r.stations} {r.pha} (n={int(r.n):,})")
ax.set_xscale("log"); ax.set_xticks(tols); ax.set_xticklabels([f"{t:g}" for t in tols])
ax.set_xlabel("tolerance (s)"); ax.set_ylabel("recall of manual picks, covered station-days"); ax.set_ylim(0, 1)
ax.legend(frameon=False, fontsize=8); ax.set_title("obs campaign vs AACSE analysts", loc="left", fontsize=10.5)

ax = axes[2]
ax.plot(range(len(axial)), axial.detected, marker="o", ms=5, lw=1.8, color="#2a78d6", label="P on >= 3 stations")
ax.plot(range(len(axial)), axial.detected_conf05, marker="o", ms=5, lw=1.8, color="#eb6834", label="same, conf >= 0.5")
for i, n in enumerate(axial.n):
    ax.text(i, 0.02, f"{int(n):,}", ha="center", fontsize=6.5, color="#555", rotation=90, va="bottom")
ax.set_xticks(range(len(axial))); ax.set_xticklabels(axial.magnitude, fontsize=8); ax.set_xlabel("LDEO magnitude")
ax.set_ylim(0, 1); ax.set_ylabel("fraction of LDEO events detected"); ax.legend(frameon=False, fontsize=8)
ax.set_title("Axial Seamount, event level", loc="left", fontsize=10.5)
fig.tight_layout()

print("AACSE 2018: the five models re-scored on the same analyst windows, plus the campaign's stored pick")
resc.round(3)

AACSE 2018: the five models re-scored on the same analyst windows, plus the campaign's stored pick

	conf>=0.3 P	conf>=0.3 S	conf>=0.2 P	conf>=0.2 S	median_P_dt_s
weights
pickblue_phasenet	0.695	0.651	0.770	0.737	0.03
pickblue_eqt	0.652	0.569	0.735	0.675	0.04
obstransformer	0.752	0.842	0.805	0.866	0.07
quakescope2026	0.638	0.292	0.742	0.459	0.06
original	0.448	0.478	0.488	0.560	0.10
campaign (stored)	0.802	0.703	0.825	0.809	0.01

6. Does the stored catalogue reproduce?

Station-days of the western campaign re-picked through the production code path with ObsPy/FDSN as the data source on a different CPU architecture, and matched to the stored Parquet on (station, band, phase, peak time) exactly.

rep = load("western_reproduction", "per_station_day")
if rep is None:
    print("western_reproduction/ not present: execute tutorials/western_pick_validation.ipynb first")
else:
    m = meta["western_reproduction"]["totals"]
    print(f"{m['station_days_with_picks']} of {m['station_days_targeted']} targeted station-days hold campaign picks; on those:")
    print(f"  campaign picks {m['campaign_picks']:,}, re-picked {m['repicked']:,}, matched exactly {m['matched_exact']:,} "
          f"({m['matched_exact'] / m['campaign_picks']:.4%} of the campaign's), campaign-only {m['campaign_only']}, re-pick-only {m['repick_only']:,}")
    r = rep[rep.n_prod > 0].sort_values(["sequence", "kind", "tid"]).reset_index(drop=True)
    fig, ax = plt.subplots(figsize=(9, 0.22 * len(r) + 1.2))
    y = np.arange(len(r))
    ax.barh(y, r.recall, color=[C_P if k == "mainshock" else C_S for k in r.kind], height=0.72)
    for i, row in r.iterrows():
        ax.text(min(row.recall, 1.0) + 0.0004, i, f"{int(row.n_prod):,}", va="center", fontsize=7, color="#3d3d3d")
    ax.set_yticks(y); ax.set_yticklabels([f"{a}  {b}  {c}" for a, b, c in zip(r.sequence, r.tid, r.day)], fontsize=7)
    ax.set_xlim(0.95, 1.004); ax.set_xlabel("fraction of the campaign's picks reproduced exactly (axis starts at 0.95)")
    ax.legend(handles=[plt.Rectangle((0, 0), 1, 1, color=C_P, label="mainshock day"),
                       plt.Rectangle((0, 0), 1, 1, color=C_S, label="quiet day")], frameon=False, fontsize=8, loc="lower left")
    ax.set_title("Exact reproduction by station-day; the number is the campaign's pick count", loc="left", fontsize=10.5)
    ax.grid(axis="y"); fig.tight_layout()

33 of 52 targeted station-days hold campaign picks; on those:
  campaign picks 59,315, re-picked 69,063, matched exactly 59,298 (99.9713% of the campaign's), campaign-only 17, re-pick-only 9,765

7. Everything in one table

Recall at the matched budget where one exists, else at the shared threshold, for the three land weights and the incumbent; the sample size; and the study each row comes from. This is the table to cite.

summary = budget.copy()
summary["n_analyst"] = [int(bench[(bench.phase == p) & (bench.sequence == s)].analyst.iloc[0]) for p, s in summary.index]
summary["study"] = ["global_sequences" if bench[bench.sequence == s].region.iloc[0] == "abroad" else "us_sequences"
                    for p, s in summary.index]
cols = ["budget", "n_analyst"] + WEIGHTS + ["study"]
summary = summary[cols]
summary.to_csv(RES / "summary_matched_budget.csv")
summary.round(3)

		budget	n_analyst	quakescope2026	jma_wc	original	instance	study
phase	sequence
P	Ridgecrest	462	347	0.689	0.671	0.631	NaN	us_sequences
	San Simeon	395	62	0.830	0.817	0.857	0.903	us_sequences
	Mendocino 2024	490	200	0.832	0.814	0.690	0.888	us_sequences
	Kaikoura 2016	1214	357	0.560	0.587	0.588	0.698	global_sequences
	Norcia 2016	910	701	0.660	0.680	0.621	0.773	global_sequences
	Thessaly 2021	1333	395	0.662	0.688	0.653	0.756	global_sequences
S	Ridgecrest	411	298	0.654	0.635	0.624	NaN	us_sequences
	Mendocino 2024	320	113	0.580	0.579	0.616	0.650	us_sequences
	Kaikoura 2016	902	458	0.452	0.452	0.423	0.667	global_sequences
	Norcia 2016	936	656	0.651	0.669	0.612	0.849	global_sequences
	Thessaly 2021	1063	331	0.401	0.403	0.359	0.445	global_sequences

---

Regenerate: execute the five benchmark notebooks (their runtimes are in reports/README.md), each of which writes its tables to docs/benchmark/results/<study>/, then execute this notebook. meta.json in each folder records when and with which SeisBench version.