import io
from collections import Counter, defaultdict
from datetime import timezone
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import obspy
import pandas as pd
import seisbench.models as sbm
from obspy import UTCDateTime
from obspy.clients.fdsn import Client
from obspy.geodetics import gps2dist_azimuth
from s3fs import S3FileSystem

%matplotlib inline

# Five sequences, chosen to vary region, network, magnitude, and era.
# Stations were selected by testing that data actually returns, not from
# metadata alone - several catalogued stations serve nothing for these dates.
SEQUENCES = {
    "Ridgecrest": dict(
        time=UTCDateTime("2019-07-06T03:19:53"), lat=35.770, lon=-117.599, mag=7.1,
        picks_from=["SCEDC", "NCEDC"], min_mag=3.2, window_min=30,
        catalog_csv="comcat_reviewed_ridgecrest.csv",
        stations=[("CI", "CLC", "HH"), ("CI", "TOW2", "HH"), ("CI", "SRT", "HH"),
                  ("CI", "WRC2", "HH"), ("CI", "JRC2", "HH")],
        note="dense aftershock sequence, events seconds apart",
    ),
    "San Simeon": dict(
        time=UTCDateTime("2003-12-22T19:15:56"), lat=35.700, lon=-121.100, mag=6.5,
        picks_from=["SCEDC", "NCEDC"], min_mag=2.0, window_min=120,
        catalog_csv="comcat_reviewed_san_simeon.csv",
        stations=[("BK", "PKD", "HH"), ("CI", "PHL", "HH"), ("CI", "SMM", "HH"),
                  ("BK", "SAO", "HH"), ("CI", "LCP", "HH")],
        note="2003 - sparser network; SCEDC and NCEDC locate this one with "
             "different station sets, so both are queried",
    ),
    "Monte Cristo": dict(
        time=UTCDateTime("2020-05-15T11:03:27"), lat=38.169, lon=-117.850, mag=6.5,
        picks_from=["SCEDC", "NCEDC"], min_mag=2.0, window_min=120,
        stations=[("NN", "BRS2", "HH"), ("NN", "LHV", "HH"), ("NN", "ION4", "HH"),
                  ("NN", "DSP", "HH"), ("NN", "Q09A", "HH")],
        note="Nevada, Basin and Range - a different crust and a different network",
    ),
    "Mendocino 2024": dict(
        time=UTCDateTime("2024-12-05T18:44:21"), lat=40.374, lon=-125.022, mag=7.0,
        picks_from=["SCEDC", "NCEDC"], min_mag=2.0, window_min=120,
        catalog_csv="comcat_reviewed_mendocino.csv",
        stations=[("BK", "PETL", "HH"), ("NC", "KCT", "HH"), ("BK", "WLKR", "HH"),
                  ("NC", "KMPB", "HH"), ("BK", "RBOW", "HH")],
        note="offshore - every station is one-sided and 55 km or more away",
    ),
    "Monroe WA": dict(
        time=UTCDateTime("2019-07-12T09:51:38"), lat=47.873, lon=-122.016, mag=4.6,
        picks_from=None, min_mag=1.0, window_min=30,
        stations=[("UW", "SP2", "HH"), ("UW", "BERY", "HH"), ("UW", "RATT", "HH"),
                  ("UW", "BST16", "HH"), ("UW", "BST20", "HH")],
        note="Cascadia, moderate magnitude - no analyst arrivals published, so "
             "this one is scored on model agreement only",
    ),
}

# Each network is served by its home archive. Asking the wrong one silently
# returns nothing, which is how station lists end up looking empty.
#
# CI, BK and NC are read straight from the SCEDC and NCEDC public S3 buckets -
# the same path QuakeScope uses in production, and faster than FDSN. UW and NN
# have no equivalent open bucket, so they come over FDSN.
S3_BUCKET = {"CI": "scedc", "BK": "ncedc", "NC": "ncedc", "NP": "ncedc"}
ROUTE = {"CI": "SCEDC", "BK": "NCEDC", "NC": "NCEDC", "NP": "NCEDC",
         "UW": "EARTHSCOPE", "NN": "EARTHSCOPE", "LB": "EARTHSCOPE",
         "IM": "EARTHSCOPE"}

# Channel band, best first. 2003 is the reason this is a list: the SCEDC bucket
# carries only BH and LH for the CI stations that year, so a hard-coded HH
# would quietly drop them.
CHANNEL_PREFERENCE = ["HH", "BH"]

WEIGHTS = ["quakescope2026", "jma_wc", "original", "instance"]

# Aftershock window starts after the mainshock coda. Its length is set per
# sequence, because aftershock productivity and how much of it the analysts
# worked through vary enormously - 30 minutes at Ridgecrest yields a larger
# reference than two hours at San Simeon.
WINDOW_START = 600      # seconds after origin, all sequences

# Inference runs once at a low floor, keeping each pick's peak probability, so
# any threshold can be applied afterwards without re-running the models. The
# headline tables use REPORT_THRESHOLD; section 8 sweeps the rest.
DETECT_FLOOR = 0.02
REPORT_THRESHOLD = 0.3
THRESHOLD_SWEEP = [0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7]
MATCH_TOL = 0.5         # seconds, for counting an analyst pick as recovered

# CVD-validated categorical slots; every chart also carries direct value
# labels, so identity never rests on colour alone.
COLORS = ["#2a78d6", "#eb6834", "#1baf7a", "#eda100"]
C_P, C_S = "#2a78d6", "#eb6834"

_fs = S3FileSystem(anon=True)
_clients = {}


def _scedc_key(net, sta, cha, comp, year, doy, loc=""):
    """SCEDC: station padded to five, location padded to three, one file per channel-day."""
    base = f"{net}{sta.ljust(5, '_')}{cha}{comp}{loc.ljust(3, '_')}{year}{doy:03d}.ms"
    return f"scedc-pds/continuous_waveforms/{year}/{year}_{doy:03d}/{base}"


def _ncedc_key(net, sta, cha, comp, year, doy, loc=""):
    """NCEDC: network is a directory level and the day folder uses a dot."""
    return (f"ncedc-pds/continuous_waveforms/{net}/{year}/{year}.{doy:03d}/"
            f"{sta}.{net}.{cha}{comp}.{loc}.D.{year}.{doy:03d}")


_KEY_FN = {"scedc": _scedc_key, "ncedc": _ncedc_key}


_listings = {}


def _day_listing(bucket, net, year, doy):
    """Object names for one network-day, listed once and reused.

    Listing the whole day prefix rather than globbing per station does three
    things: it finds the location code instead of guessing it, it lets a
    missing channel be detected without a failed request, and it avoids an
    s3fs trap - a glob populates the directory cache with only the matching
    entries, after which opening a sibling object raises FileNotFoundError
    even though it exists.
    """
    if bucket == "scedc":
        prefix = f"scedc-pds/continuous_waveforms/{year}/{year}_{doy:03d}/"
    else:
        prefix = f"ncedc-pds/continuous_waveforms/{net}/{year}/{year}.{doy:03d}/"
    if prefix not in _listings:
        try:
            _listings[prefix] = {k.split("/")[-1] for k in _fs.ls(prefix)}
        except Exception:
            _listings[prefix] = set()
    return prefix, _listings[prefix]


def _object_name(bucket, net, sta, cha, comp, year, doy, names):
    """The object holding this channel-day, whatever its location code."""
    if bucket == "scedc":
        head = f"{net}{sta.ljust(5, '_')}{cha}{comp}"
        tail = f"{year}{doy:03d}.ms"
        hits = [n for n in names if n.startswith(head) and n.endswith(tail)]
    else:
        head = f"{sta}.{net}.{cha}{comp}."
        tail = f".D.{year}.{doy:03d}"
        hits = [n for n in names if n.startswith(head) and n.endswith(tail)]
    return sorted(hits)[0] if hits else None


def read_s3(net, sta, cha, t0, t1):
    """Whole-day objects from the public bucket, trimmed to the window."""
    bucket = S3_BUCKET.get(net)
    if bucket is None:
        return None
    year, doy = t0.year, int(t0.strftime("%j"))
    prefix, names = _day_listing(bucket, net, year, doy)
    if not names:
        return None
    st = obspy.Stream()
    for comp in "ZNE":
        name = _object_name(bucket, net, sta, cha, comp, year, doy, names)
        if name is None:
            return None                      # incomplete three-component set
        try:
            with _fs.open(prefix + name) as fh:
                st += obspy.read(io.BytesIO(fh.read()))
        except Exception:
            return None
    st.merge(fill_value=0)
    st.trim(t0, t1)
    return st if len(st) >= 3 else None


def client_for(net):
    """FDSN client for whichever archive holds this network."""
    provider = ROUTE.get(net, "EARTHSCOPE")
    if provider not in _clients:
        _clients[provider] = Client(provider, timeout=300)
    return _clients[provider]


def fetch_sequence(seq):
    """Waveforms for every station of one sequence, over the aftershock window."""
    t0 = seq["time"] + WINDOW_START
    t1 = t0 + seq["window_min"] * 60
    streams, rates = {}, {}
    for net, sta, _ in seq["stations"]:
        got = None
        for cha in CHANNEL_PREFERENCE:
            got = read_s3(net, sta, cha, t0, t1)          # bucket first
            if got is None and net not in S3_BUCKET:       # no bucket: FDSN
                try:
                    got = client_for(net).get_waveforms(net, sta, "*", cha + "?", t0, t1)
                    got.merge(fill_value=0)
                    if len(got) < 3:
                        got = None
                except Exception:
                    got = None
            if got is not None:
                break
        if got is None:
            print(f"    {net}.{sta}: no data")
            continue
        expected = seq["window_min"] * 60 * got[0].stats.sampling_rate
        if got[0].stats.npts < 0.5 * expected:
            print(f"    {net}.{sta}: incomplete")
            continue
        key = f"{net}.{sta}"
        streams[key] = got
        rates[key] = got[0].stats.sampling_rate
    return streams, t0, t1, rates


CATALOG_DIR = Path("../docs/rerun_2026")   # team-curated reviewed event lists
ORIGIN_LEAD = 180     # s: an origin this far before t0 can still put arrivals inside it
# Tolerance for calling two origins the same earthquake. Ridgecrest aftershocks
# arrive 0.35 s apart at the tightest, so time alone is not enough to identify an
# event: at 5 s, 17 returned origins matched more than one curated event and the
# coverage count silently undercounted. Requiring proximity too brings that to 2.
ORIGIN_TOL = 2.0      # s
ORIGIN_DIST_KM = 30.0
CATALOG_RADIUS = 1.5  # deg: generous, because the curated list does the filtering


def _curated_origins(seq, t0, t1):
    """Origin times of the team's reviewed events for this sequence."""
    df = pd.read_csv(CATALOG_DIR / seq["catalog_csv"])
    df = df[df["status"] == "reviewed"].copy()
    df["t"] = pd.to_datetime(df["time"], format="mixed", utc=True)
    lo = (t0 - ORIGIN_LEAD).datetime.replace(tzinfo=timezone.utc)
    hi = t1.datetime.replace(tzinfo=timezone.utc)
    sel = df[(df["t"] >= lo) & (df["t"] <= hi)]
    return [(UTCDateTime(r["t"].isoformat()), r["latitude"], r["longitude"])
            for _, r in sel.iterrows()]


def _match_curated(origin, curated):
    """Index of the curated event this origin is, or None. Nearest in time among
    those within both tolerances, so neighbouring aftershocks are not conflated."""
    best, best_dt = None, None
    for j, (ct, cla, clo) in enumerate(curated):
        dt = abs(origin.time - ct)
        if dt > ORIGIN_TOL:
            continue
        if gps2dist_azimuth(origin.latitude, origin.longitude,
                            cla, clo)[0] / 1000.0 > ORIGIN_DIST_KM:
            continue
        if best_dt is None or dt < best_dt:
            best, best_dt = j, dt
    return best


def analyst_picks(seq, t0, t1):
    """Reviewed analyst P and S picks in the window, keyed by (station, phase).

    The curated list in `docs/rerun_2026/comcat_reviewed_*.csv` decides **which
    earthquakes count** - every row is a `status == reviewed` ComCat origin the
    team selected. Arrivals are then harvested from **every** archive in
    `picks_from` that independently located those earthquakes, matched to the
    list by origin time within ORIGIN_TOL.

    Both halves matter. Anchoring the event set on a committed file is what
    makes the reference reproducible: an identical live radius/magnitude query
    returned 110 events one day and 123 another, because the event services
    serve a revisable catalog. Harvesting from every archive is what makes it
    complete: SCEDC and NCEDC locate the same California earthquake with
    different station sets, and ComCat carries only the authoritative solution.
    Taking the authoritative origin alone costs San Simeon 18 of its 20 S picks.

    Only `evaluation_mode == "manual"` picks are kept, so the reference is
    analyst work rather than the preliminary automatic picks the catalogs also
    carry. Duplicates are collapsed at MATCH_TOL, the tolerance used to score
    the models.
    """
    archives = seq["picks_from"]
    if not archives:
        return None
    if isinstance(archives, str):
        archives = [archives]
    curated = _curated_origins(seq, t0, t1) if seq.get("catalog_csv") else None
    keep = {f"{n}.{s}" for n, s, _ in seq["stations"]}

    out = defaultdict(list)
    stats = Counter()
    covered = set()
    for archive in archives:
        try:
            kwargs = dict(starttime=t0 - ORIGIN_LEAD, endtime=t1,
                          latitude=seq["lat"], longitude=seq["lon"],
                          includearrivals=True)
            if curated is None:      # no list: fall back to a magnitude floor
                kwargs.update(maxradius=1.0, minmagnitude=seq["min_mag"])
            else:
                kwargs.update(maxradius=CATALOG_RADIUS)
            cat = Client(archive, timeout=300).get_events(**kwargs)
        except Exception as exc:
            print(f"    {archive}: no catalog ({type(exc).__name__})")
            continue

        matched = added = 0
        for ev in cat:
            origin = ev.preferred_origin() or (ev.origins[0] if ev.origins else None)
            if origin is None:
                continue
            if curated is not None:
                hit = _match_curated(origin, curated)
                if hit is None:
                    stats["origin not in the curated list"] += 1
                    continue
                covered.add(hit)
            matched += 1
            by_id = {p.resource_id.id: p for p in ev.picks}
            for arr in origin.arrivals:
                pick = by_id.get(arr.pick_id.id)
                if pick is None or not arr.phase:
                    continue
                phase = arr.phase[0].upper()
                if phase not in ("P", "S"):
                    continue
                wid = pick.waveform_id
                key = f"{wid.network_code}.{wid.station_code}"
                if key not in keep or not (t0 <= pick.time <= t1):
                    continue
                if str(pick.evaluation_mode) == "automatic":
                    stats["automatic, dropped"] += 1
                    continue
                if any(abs(pick.time - u) <= MATCH_TOL for u in out[(key, phase)]):
                    stats["duplicate across archives"] += 1
                    continue
                out[(key, phase)].append(pick.time)
                added += 1
        note = f", {matched} matched the curated list" if curated is not None else ""
        print(f"    {archive}: {len(cat)} events{note}, +{added} picks on these stations")

    if curated is not None:
        print(f"    curated list: {len(covered)}/{len(curated)} reviewed events "
              f"located by at least one archive")
    if stats:
        print("    (" + ", ".join(f"{v} {k}" for k, v in stats.items()) + ")")
    return {k: sorted(v) for k, v in out.items()} if out else None


def at_threshold(store, thr):
    """Pick times from a (time, probability) store, at a given threshold."""
    return {k: sorted(t for t, v in vals if v >= thr) for k, vals in store.items()}


def match(reference, candidate, tol=MATCH_TOL):
    """Greedy nearest match; each candidate pick is consumed at most once."""
    used, residuals = set(), []
    for a in reference:
        best_i = best_d = None
        for i, m in enumerate(candidate):
            if i in used:
                continue
            d = m - a
            if abs(d) <= tol and (best_d is None or abs(d) < abs(best_d)):
                best_i, best_d = i, d
        if best_i is not None:
            used.add(best_i)
            residuals.append(best_d)
    return residuals, len(candidate) - len(used)

available = sbm.PhaseNet.list_pretrained()
models = {}
for name in WEIGHTS:
    if name not in available:
        print(f"{name:<16} not installed - skipping")
        continue
    try:
        models[name] = sbm.PhaseNet.from_pretrained(name)
        print(f"{name:<16} loaded")
    except Exception as exc:
        print(f"{name:<16} could not load ({type(exc).__name__}) - skipping")

if not models:
    raise RuntimeError("no weight sets available")
names = list(models)

2026-08-18 20:30:07,837 | seisbench | WARNING | The download precheck failed with a timeout. This is not an error itself, but might indicate a subsequent error. If you encounter an error, this might be caused by the firewall setup of your network. Please check https://github.com/seisbench/seisbench#known-issues for details.

As a mitigation, you might be able to switch to our backup repository at the cost of lower download speeds. To switch, run 'seisbench.use_backup_repository'.

quakescope2026   loaded
jma_wc           loaded
original         loaded
instance         loaded

results = {}
for label, seq in SEQUENCES.items():
    print(f"{label} (M{seq['mag']}, {seq['note']})")
    streams, t0, t1, rates = fetch_sequence(seq)
    print(f"    {len(streams)}/{len(seq['stations'])} stations: "
          + ", ".join(f"{k.split(chr(46))[1]}@{v:g}Hz" for k, v in rates.items()))
    if not streams:
        continue
    reference = analyst_picks(seq, t0, t1)
    if reference is not None:
        n_p = sum(len(v) for (_, ph), v in reference.items() if ph == 'P')
        n_s = sum(len(v) for (_, ph), v in reference.items() if ph == 'S')
        print(f"    analyst picks: {n_p} P, {n_s} S")
    else:
        print("    no analyst arrivals - model agreement only")

    picks = {}
    for name, model in models.items():
        per_station = defaultdict(list)
        for sta, st in streams.items():
            try:
                out = model.classify(st, P_threshold=DETECT_FLOOR,
                                     S_threshold=DETECT_FLOOR)
            except Exception as exc:
                print(f"    {name} {sta}: {type(exc).__name__}")
                continue
            for p in out.picks:
                per_station[(sta, p.phase)].append((p.peak_time, float(p.peak_value)))
        picks[name] = dict(per_station)
        tot = sum(1 for v in per_station.values() for t, q in v if q >= REPORT_THRESHOLD)
        print(f"    {name:<16} {tot} picks at {REPORT_THRESHOLD}")

    results[label] = dict(streams=streams, picks=picks, reference=reference,
                          t0=t0, t1=t1, rates=rates)
    print()

Ridgecrest (M7.1, dense aftershock sequence, events seconds apart)

    5/5 stations: CLC@100Hz, TOW2@100Hz, SRT@100Hz, WRC2@100Hz, JRC2@100Hz

    SCEDC: 141 events, 138 matched the curated list, +645 picks on these stations

    NCEDC: 2 events, 2 matched the curated list, +0 picks on these stations
    curated list: 138/138 reviewed events located by at least one archive
    (3 origin not in the curated list, 1 duplicate across archives)
    analyst picks: 347 P, 298 S

    quakescope2026   643 picks at 0.3

    jma_wc           617 picks at 0.3

    original         1010 picks at 0.3

    instance         198 picks at 0.3

San Simeon (M6.5, 2003 - sparser network; SCEDC and NCEDC locate this one with different station sets, so both are queried)

    5/5 stations: PKD@80Hz, PHL@20Hz, SMM@20Hz, SAO@80Hz, LCP@20Hz

    SCEDC: 19 events, 14 matched the curated list, +57 picks on these stations

rows = []
for label, res in results.items():
    if res['reference'] is None:
        continue
    for name in names:
        if name not in res['picks']:
            continue
        thresholded = {name: at_threshold(res['picks'][name], REPORT_THRESHOLD)}
        for phase in ('P', 'S'):
            n_ref = n_hit = n_extra = 0
            residuals = []
            for sta in res['streams']:
                ref = res['reference'].get((sta, phase), [])
                got = thresholded[name].get((sta, phase), [])
                r, extra = match(ref, got)
                n_ref += len(ref); n_hit += len(r); n_extra += extra
                residuals += r
            if n_ref == 0:
                continue
            rows.append(dict(
                sequence=label, weights=name, phase=phase,
                analyst=n_ref, matched=n_hit,
                recall=round(n_hit / n_ref, 3),
                MAE=round(float(np.mean(np.abs(residuals))), 3) if residuals else np.nan,
                extra=n_extra,
            ))

bench = pd.DataFrame(rows)
if len(bench):
    print(bench.to_string(index=False))
else:
    print('no analyst-referenced results')

      sequence        weights phase  analyst  matched  recall   MAE  extra
    Ridgecrest quakescope2026     P      347      210   0.605 0.040    142
    Ridgecrest quakescope2026     S      298      154   0.517 0.053    137
    Ridgecrest         jma_wc     P      347      202   0.582 0.033    136
    Ridgecrest         jma_wc     S      298      145   0.487 0.064    134
    Ridgecrest       original     P      347      225   0.648 0.058    265
    Ridgecrest       original     S      298      216   0.725 0.066    304
    Ridgecrest       instance     P      347       75   0.216 0.048     25
    Ridgecrest       instance     S      298       54   0.181 0.079     44
    San Simeon quakescope2026     P       62       51   0.823 0.067    322
    San Simeon quakescope2026     S       16       12   0.750 0.020    179
    San Simeon         jma_wc     P       62       53   0.855 0.066    556
    San Simeon         jma_wc     S       16       12   0.750 0.035    225
    San Simeon       original     P       62       53   0.855 0.058    331
    San Simeon       original     S       16       15   0.938 0.103    389
    San Simeon       instance     P       62       55   0.887 0.083    243
    San Simeon       instance     S       16       12   0.750 0.094    104
  Monte Cristo quakescope2026     P       16       12   0.750 0.047    714
  Monte Cristo quakescope2026     S       12        5   0.417 0.219    457
  Monte Cristo         jma_wc     P       16       13   0.812 0.051    856
  Monte Cristo         jma_wc     S       12        6   0.500 0.236    504
  Monte Cristo       original     P       16       13   0.812 0.049    734
  Monte Cristo       original     S       12        8   0.667 0.146    933
  Monte Cristo       instance     P       16       14   0.875 0.089    453
  Monte Cristo       instance     S       12        6   0.500 0.174    270
Mendocino 2024 quakescope2026     P      200      166   0.830 0.071    318
Mendocino 2024 quakescope2026     S      113       70   0.619 0.086    297
Mendocino 2024         jma_wc     P      200      180   0.900 0.074    544
Mendocino 2024         jma_wc     S      113       77   0.681 0.101    353
Mendocino 2024       original     P      200      136   0.680 0.119    283
Mendocino 2024       original     S      113       83   0.735 0.124    433
Mendocino 2024       instance     P      200      168   0.840 0.089    199
Mendocino 2024       instance     S      113       63   0.558 0.120    149

if len(bench):
    for phase in ('S', 'P'):
        sub = bench[bench.phase == phase]
        if not len(sub):
            continue
        piv = sub.pivot(index='sequence', columns='weights', values='recall')
        piv = piv.reindex(columns=[n for n in names if n in piv.columns])
        print(f'\n{phase} recall by sequence')
        print(piv.to_string())

S recall by sequence
weights         quakescope2026  jma_wc  original  instance
sequence                                                  
Mendocino 2024           0.619   0.681     0.735     0.558
Monte Cristo             0.417   0.500     0.667     0.500
Ridgecrest               0.517   0.487     0.725     0.181
San Simeon               0.750   0.750     0.938     0.750

P recall by sequence
weights         quakescope2026  jma_wc  original  instance
sequence                                                  
Mendocino 2024           0.830   0.900     0.680     0.840
Monte Cristo             0.750   0.812     0.812     0.875
Ridgecrest               0.605   0.582     0.648     0.216
San Simeon               0.823   0.855     0.855     0.887

if len(bench):
    sub = bench[bench.phase == 'S']
    seqs = [s for s in SEQUENCES if s in set(sub.sequence)]
    fig, ax = plt.subplots(figsize=(11, 4.4))
    width = 0.8 / max(len(names), 1)
    for i, name in enumerate(names):
        vals, xs = [], []
        for j, s in enumerate(seqs):
            row = sub[(sub.sequence == s) & (sub.weights == name)]
            if len(row):
                vals.append(float(row['recall'].iloc[0]))
                xs.append(j + (i - (len(names) - 1) / 2) * width)
        bars = ax.bar(xs, vals, width=width * 0.92, color=COLORS[i], label=name)
        for x, v in zip(xs, vals):
            ax.text(x, v + 0.015, f'{v:.2f}', ha='center', fontsize=8)
    ax.set_xticks(range(len(seqs)))
    ax.set_xticklabels(seqs, fontsize=9)
    ax.set_ylabel('S recall against analyst picks')
    ax.set_ylim(0, 1.05)
    ax.grid(alpha=0.25, lw=0.5, axis='y')
    ax.set_axisbelow(True)
    ax.legend(frameon=False, fontsize=9, ncol=len(names))
    ax.set_title('S recall by sequence', fontsize=11, loc='left')
    fig.tight_layout()
    plt.show()

no_ref = [l for l, r in results.items() if r['reference'] is None]
for label in no_ref:
    res = results[label]
    print(f'{label}')
    for name in names:
        if name not in res['picks']:
            continue
        thr_v = at_threshold(res['picks'][name], REPORT_THRESHOLD)
        p = sum(len(v) for (_, ph), v in thr_v.items() if ph == 'P')
        s = sum(len(v) for (_, ph), v in thr_v.items() if ph == 'S')
        print(f'    {name:<16} {p:4d} P, {s:4d} S')
    print()
    pairs = [(a, b) for i, a in enumerate(names) for b in names[i + 1:]]
    rows = []
    for a, b in pairs:
        if a not in res['picks'] or b not in res['picks']:
            continue
        for phase in ('P', 'S'):
            shared = tot_a = 0
            for sta in res['streams']:
                pa = at_threshold(res['picks'][a], REPORT_THRESHOLD).get((sta, phase), [])
                pb = at_threshold(res['picks'][b], REPORT_THRESHOLD).get((sta, phase), [])
                r, _ = match(pa, pb)
                shared += len(r); tot_a += len(pa)
            rows.append(dict(pair=f'{a} vs {b}', phase=phase,
                             shared=shared, of_first=tot_a,
                             frac=round(shared / tot_a, 3) if tot_a else np.nan))
    print(pd.DataFrame(rows).to_string(index=False))
    print('\nfrac = fraction of the first model\'s picks that the second also found,'
          f' within {MATCH_TOL} s')

Monroe WA
    quakescope2026     39 P,   25 S
    jma_wc             60 P,   37 S
    original           24 P,   30 S
    instance           25 P,   12 S

                      pair phase  shared  of_first  frac
  quakescope2026 vs jma_wc     P      34        39 0.872
  quakescope2026 vs jma_wc     S      24        25 0.960
quakescope2026 vs original     P      21        39 0.538
quakescope2026 vs original     S      17        25 0.680
quakescope2026 vs instance     P      19        39 0.487
quakescope2026 vs instance     S       9        25 0.360
        jma_wc vs original     P      22        60 0.367
        jma_wc vs original     S      23        37 0.622
        jma_wc vs instance     P      20        60 0.333
        jma_wc vs instance     S       9        37 0.243
      original vs instance     P      17        24 0.708
      original vs instance     S      11        30 0.367

frac = fraction of the first model's picks that the second also found, within 0.5 s

def record_section(label, res, weight, span=300):
    seq = SEQUENCES[label]
    fig, ax = plt.subplots(figsize=(11, 5))
    order = sorted(res['streams'],
                   key=lambda s: [f'{n}.{t}' for n, t, _ in seq['stations']].index(s)
                   if s in [f'{n}.{t}' for n, t, _ in seq['stations']] else 99)
    for row, sta in enumerate(order):
        tr = res['streams'][sta].select(component='Z')
        if not tr:
            continue
        tr = tr[0].copy()
        tr.trim(res['t0'], res['t0'] + span)
        x = tr.data.astype(float)
        peak = np.abs(x).max()
        if peak > 0:
            x = x / peak * 0.42
        ax.plot(tr.times(), x + row, color='#3d3d3d', lw=0.45)
        ax.annotate(sta, (-span * 0.015, row + 0.16), fontsize=8,
                    color='#52514e', ha='right')
        for (s2, phase), times in at_threshold(res['picks'][weight], REPORT_THRESHOLD).items():
            if s2 != sta:
                continue
            for t in times:
                dt = t - res['t0']
                if 0 <= dt <= span:
                    ax.plot([dt, dt], [row - 0.34, row + 0.34],
                            color=C_P if phase == 'P' else C_S,
                            lw=1.1, ls='-' if phase == 'P' else '--', alpha=0.85)
    handles = [plt.Line2D([], [], color=C_P, lw=1.4, label='P'),
               plt.Line2D([], [], color=C_S, lw=1.4, ls='--', label='S')]
    ax.legend(handles=handles, frameon=False, fontsize=9, ncol=2, loc='upper right')
    ax.set_xlim(-span * 0.06, span)
    ax.set_ylim(-0.8, len(order) - 0.2)
    ax.set_yticks([])
    ax.set_xlabel(f'seconds into the aftershock window '
                  f'(starts {WINDOW_START // 60} min after the mainshock)')
    ax.set_title(f'{label} - M{seq["mag"]} - picks from {weight}',
                 fontsize=11, loc='left')
    ax.grid(alpha=0.2, lw=0.5, axis='x')
    fig.tight_layout()
    return fig


for label, res in results.items():
    best = max(res['picks'], key=lambda n: sum(
        1 for v in res['picks'][n].values() for t, q in v if q >= REPORT_THRESHOLD))
    record_section(label, res, best)
    plt.show()

sweep_rows = []
for label, res in results.items():
    if res['reference'] is None:
        continue
    for name in names:
        if name not in res['picks']:
            continue
        for thr in THRESHOLD_SWEEP:
            view = at_threshold(res['picks'][name], thr)
            for phase in ('P', 'S'):
                hit = tot = emitted = 0
                for sta in res['streams']:
                    ref = res['reference'].get((sta, phase), [])
                    got = view.get((sta, phase), [])
                    r, _ = match(ref, got)
                    hit += len(r); tot += len(ref); emitted += len(got)
                if tot < 20:          # too few analyst picks to be meaningful
                    continue
                sweep_rows.append(dict(sequence=label, weights=name, phase=phase,
                                       thr=thr, recall=hit / tot, emitted=emitted))
sweep = pd.DataFrame(sweep_rows)
print(f'{len(sweep)} rows; sequences with enough analyst picks to sweep: '
      f"{sorted(set(sweep.sequence))}")

180 rows; sequences with enough analyst picks to sweep: ['Mendocino 2024', 'Ridgecrest', 'San Simeon']

def matched_budget(sweep, phase, sequence, n_points=4):
    """Recall for each model at a common number of emitted picks."""
    sub = sweep[(sweep.phase == phase) & (sweep.sequence == sequence)]
    present = [n for n in names if n in set(sub.weights)]
    if len(present) < 2:
        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:
        ceiling = min(present, key=lambda n: sub[sub.weights == n].emitted.max())
        cap = int(sub[sub.weights == ceiling].emitted.max())
        print(f'{sequence}: no common budget - {ceiling} tops out at {cap} picks, '
              f'below where the others start. See the saturation table.')
        return None
    rows = []
    for target in np.linspace(lo, hi, n_points):
        row = {'picks_emitted': int(round(target))}
        for n in present:
            d = sub[sub.weights == n].sort_values('emitted')
            row[n] = round(float(np.interp(target, d.emitted, d.recall)), 3)
        rows.append(row)
    return pd.DataFrame(rows)


for sequence in sorted(set(sweep.sequence)):
    tab = matched_budget(sweep, 'S', sequence)
    if tab is None:
        continue
    print(f'{sequence} - S recall at matched pick budgets')
    print(tab.to_string(index=False))
    print()

Mendocino 2024 - S recall at matched pick budgets
 picks_emitted  quakescope2026  jma_wc  original  instance
           141           0.396   0.372     0.389     0.483
           261           0.528   0.519     0.553     0.606
           380           0.627   0.635     0.665     0.675
           500           0.692   0.709     0.727     0.761

Ridgecrest: no common budget - instance tops out at 193 picks, below where the others start. See the saturation table.

seqs = [s for s in sorted(set(sweep.sequence))]
fig, axes = plt.subplots(1, len(seqs), figsize=(5.3 * len(seqs), 4.2), squeeze=False)
for ax, sequence in zip(axes[0], seqs):
    sub = sweep[(sweep.phase == 'S') & (sweep.sequence == sequence)]
    for name, color in zip(names, COLORS):
        d = sub[sub.weights == name].sort_values('emitted')
        if not len(d):
            continue
        ax.plot(d.emitted, d.recall, marker='o', ms=4, lw=1.8,
                color=color, label=name)
        star = d[np.isclose(d.thr, REPORT_THRESHOLD)]
        if len(star):
            ax.plot(star.emitted, star.recall, marker='*', ms=15,
                    color=color, mec='#16150f', mew=0.6, zorder=5)
    ax.set_title(f'{sequence} - S', fontsize=11, loc='left')
    ax.set_xlabel('S picks emitted')
    ax.grid(alpha=0.25, lw=0.5)
axes[0][0].set_ylabel('recall against analyst picks')
axes[0][0].legend(frameon=False, fontsize=9)
fig.suptitle('Same curve, different operating points - stars mark the shared 0.3 threshold',
             fontsize=10, x=0.01, ha='left')
fig.tight_layout()
plt.show()

sat_rows = []
for label, res in results.items():
    for name in names:
        if name not in res['picks']:
            continue
        row = {'sequence': label, 'weights': name}
        for thr in (DETECT_FLOOR, 0.1, REPORT_THRESHOLD):
            view = at_threshold(res['picks'][name], thr)
            row[f'S@{thr}'] = sum(len(v) for (_, ph), v in view.items() if ph == 'S')
        sat_rows.append(row)
sat = pd.DataFrame(sat_rows)
print('S picks emitted, by threshold')
print(sat.to_string(index=False))
print(f'\nS@{DETECT_FLOOR} is the ceiling: what the model finds with the threshold'
      ' on the floor.')

S picks emitted, by threshold
      sequence        weights  S@0.02  S@0.1  S@0.3
    Ridgecrest quakescope2026     684    459    291
    Ridgecrest         jma_wc     638    447    279
    Ridgecrest       original     832    646    520
    Ridgecrest       instance     246    149     98
    San Simeon quakescope2026    1265    579    191
    San Simeon         jma_wc    1546    726    237
    San Simeon       original     986    640    404
    San Simeon       instance     571    277    116
  Monte Cristo quakescope2026    1482    922    462
  Monte Cristo         jma_wc    1594   1023    510
  Monte Cristo       original    1952   1352    941
  Monte Cristo       instance     823    558    276
Mendocino 2024 quakescope2026    1267    757    367
Mendocino 2024         jma_wc    1555    914    430
Mendocino 2024       original    1036    726    516
Mendocino 2024       instance     659    398    212
     Monroe WA quakescope2026     138     56     25
     Monroe WA         jma_wc     215     85     37
     Monroe WA       original     137     58     30
     Monroe WA       instance      37     20     12

S@0.02 is the ceiling: what the model finds with the threshold on the floor.


`quakescope2026`	the v7 fine-tune, this project's production candidate
`jma_wc`	the SeisBench Japanese model v7 was distilled from — the baseline that matters
`original`	Zhu & Beroza (2019), the published reference
`instance`	trained on the Italian INSTANCE dataset — the weight the 2025 campaign actually ran, so it is the incumbent rather than a control

Field	Ridgecrest / SCEDC	San Simeon / NCEDC
`evaluation_mode`	manual 3195, automatic 27	manual 3977, automatic 36
`evaluation_status`	reviewed 3195, preliminary 27	final 3679, reviewed 298, preliminary 36
`creation_info.agency_id`	`CI` (all)	`NC` (all)
`creation_info.author`	`None` (all)	`None` (all)
`method_id`	`None` (all)	`None` (all)
`onset`	emergent 2411, impulsive 799	emergent 1197, impulsive 763, unset 2053
`polarity`	undecidable 2783, positive 302, negative 137	negative 1953, positive 1147, undecidable 913
`time_errors.uncertainty`	set on all	set on all

Phase pickers across five earthquake sequences¶

Weight sets¶

How it is scored¶

1. Configuration¶

2. Data access¶

3. Load the weight sets¶

4. Run every model over every sequence¶

4a. Where the reference picks come from¶

What the QuakeML actually carries¶

Two consequences for the scoring¶

The event set comes from the team's curated lists¶

Arrivals are harvested from every archive, not just the authoritative one¶

5. Recall against analyst picks¶

S recall, side by side¶

6. Monroe WA — agreement without a reference¶

7. What the picks look like¶

8. Is one model better, or just more liberal?¶

Can every model reach the same budget?¶

What the curves say¶

Reading the result¶

Sequence	Date	M	Setting
Ridgecrest	2019-07-06	7.1	Eastern California, dense aftershocks seconds apart
San Simeon	2003-12-22	6.5	Central Coast, 2003 network and instrumentation
Monte Cristo	2020-05-15	6.5	Nevada, Basin and Range, different network
Mendocino	2024-12-05	7.0	Offshore, one-sided geometry, every station 55 km+
Monroe WA	2019-07-12	4.6	Cascadia, moderate magnitude

arrivals harvested from	P	S
the authoritative archive only	28	2
every archive that located the event	61	20