seismo-relay/scratch/offset_scan.py

#!/usr/bin/env python3
"""Scan series-3 waveform binaries for the 'offset' hardware fault.

A healthy geophone trace is centred on zero.  An offset unit sits displaced,
so the channel mean approaches its own peak.  Detector (unchanged from the
2026-08-25 run, see memory note `offset-archive-analysis-backlog`):

    dominant-axis  |mean| / peak > 0.7
    AND            |mean| >= 0.9 * the unit's geo trigger level

Trigger level is read from a paired _ASCII.TXT where one exists, otherwise
from a per-serial median learned across that unit's ASCII files, otherwise
--default-trigger.

Serial is decoded from the BW filename: prefix letter encodes thousands
(chr(ord('B') + n)), next 3 digits the remainder -- T193 -> BE18193.

Usage:
  python scratch/offset_scan.py --dir <path> [--jobs N] --out offsets.csv
"""
from __future__ import annotations

import argparse, csv, json, re, sys
from collections import defaultdict
from concurrent.futures import ProcessPoolExecutor, as_completed
from pathlib import Path

sys.path.insert(0, str(Path(__file__).resolve().parent.parent))

from minimateplus.event_file_io import read_blastware_file
from minimateplus.bw_ascii_report import parse_report

GEO = ("Tran", "Vert", "Long")
_GEO_FS_COUNTS = 32000.0
_WAVE_RE = re.compile(r"\.[A-Za-z0-9]{2}0[Ww]$")
_STEM_RE = re.compile(r"^([B-Z])(\d{3})")

MEAN_OVER_PEAK_MIN = 0.7
TRIGGER_FRACTION   = 0.9


def serial_from_name(name: str):
    m = _STEM_RE.match(name)
    if not m:
        return None
    letter, digits = m.group(1), m.group(2)
    return f"BE{(ord(letter) - ord('B')) * 1000 + int(digits)}"


def counts_to_ips(c, gr):
    return c * (gr or 10.0) / _GEO_FS_COUNTS


def scan_one(path_str: str, default_trigger: float) -> dict | None:
    p = Path(path_str)
    try:
        gr, trig = 10.0, None
        ap = p.with_name(p.name.replace(".", "_", 1) + "_ASCII.TXT") \
            if False else p.parent / (p.stem + "_" + p.suffix.lstrip(".") + "_ASCII.TXT")
        if ap.exists():
            rep = parse_report(ap.read_text(errors="replace"))
            gr = rep.geo_range_ips or 10.0
            trig = rep.geo_trigger_level_ips
        ev = read_blastware_file(p)
        s = ev.raw_samples or {}
        if not all(s.get(c) for c in GEO):
            return None
        best = None
        for ch in GEO:
            arr = s[ch]
            n = len(arr)
            if n == 0:
                continue
            mean = sum(arr) / n
            peak = max(abs(v) for v in arr)
            if peak == 0:
                continue
            ratio = abs(mean) / peak
            if best is None or peak > best["peak_counts"]:
                best = {"channel": ch, "mean_counts": mean,
                        "peak_counts": peak, "ratio": ratio}
        if best is None:
            return None
        ts = ev.timestamp
        return {
            "serial": serial_from_name(p.name) or "?",
            "timestamp": (f"{ts.year:04d}-{ts.month:02d}-{ts.day:02d}T"
                          f"{ts.hour:02d}:{ts.minute:02d}:{ts.second:02d}") if ts else "",
            "filename": p.name,
            "channel": best["channel"],
            "offset_ips": round(counts_to_ips(best["mean_counts"], gr), 4),
            "peak_ips": round(counts_to_ips(best["peak_counts"], gr), 4),
            "mean_over_peak": round(best["ratio"], 3),
            "trigger_level_ips": trig if trig is not None else "",
            "geo_range_ips": gr,
        }
    except Exception:
        return None


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--dir", required=True)
    ap.add_argument("--jobs", type=int, default=4)
    ap.add_argument("--limit", type=int, default=0)
    ap.add_argument("--default-trigger", type=float, default=0.2)
    ap.add_argument("--out", required=True)
    a = ap.parse_args()

    # The DL2 export keeps a byte-identical `Sent/` mirror of the root, so
    # enumerate paths but keep only the first occurrence of each basename —
    # otherwise every event is counted twice.
    seen = set()
    files = []
    for q in sorted(Path(a.dir).rglob("*")):
        if q.is_file() and _WAVE_RE.search(q.name) and q.name not in seen:
            seen.add(q.name)
            files.append(q)
    if a.limit:
        files = files[: a.limit]
    print(f"waveform binaries to scan: {len(files)}", flush=True)

    rows = []
    with ProcessPoolExecutor(max_workers=a.jobs) as ex:
        futs = [ex.submit(scan_one, str(p), a.default_trigger) for p in files]
        for n, f in enumerate(as_completed(futs), 1):
            r = f.result()
            if r:
                rows.append(r)
            if n % 2000 == 0:
                print(f"  {n}/{len(files)}", flush=True)

    # learn per-serial trigger levels from the rows that had an ASCII
    by_serial = defaultdict(list)
    for r in rows:
        if r["trigger_level_ips"] != "":
            by_serial[r["serial"]].append(float(r["trigger_level_ips"]))
    med = {}
    for k, v in by_serial.items():
        v.sort()
        med[k] = v[len(v) // 2]

    for r in rows:
        if r["trigger_level_ips"] == "":
            r["trigger_level_ips"] = med.get(r["serial"], a.default_trigger)
        r["suspect"] = int(
            r["mean_over_peak"] > MEAN_OVER_PEAK_MIN
            and abs(r["offset_ips"]) >= TRIGGER_FRACTION * float(r["trigger_level_ips"])
        )

    cols = ["serial", "timestamp", "filename", "channel", "offset_ips", "peak_ips",
            "mean_over_peak", "trigger_level_ips", "geo_range_ips", "suspect"]
    with open(a.out, "w", newline="") as fh:
        w = csv.DictWriter(fh, fieldnames=cols)
        w.writeheader()
        w.writerows(rows)

    sus = [r for r in rows if r["suspect"]]
    print(f"\nscanned {len(rows)} decodable waveforms")
    print(f"suspect events: {len(sus)}")
    per = defaultdict(int)
    for r in sus:
        per[r["serial"]] += 1
    print(f"units with >=1 suspect event: {len(per)} of {len({r['serial'] for r in rows})}")
    for s, n in sorted(per.items(), key=lambda x: -x[1])[:20]:
        print(f"  {s:10} {n}")
    print(f"\nwrote {a.out}")


if __name__ == "__main__":
    main()