seismo-relay/scratch/verify_thor_against_csv.py

#!/usr/bin/env python3
"""Verify the Thor / Micromate (series-4) IDF decoder against Thor's own exports.

Sister harness to ``scratch/verify_against_ascii.py`` (series-3 / Blastware).

Ground truth is the ``.IDFW.csv`` / ``.IDFH.csv`` file Thor writes next to each
binary, under a sibling ``CSV/`` directory:

    <dir>/UM13981_20220207084555.IDFW
    <dir>/CSV/UM13981_20220207084555.IDFW.csv

For waveforms the CSV carries a per-sample block of four columns
(Tran, Vert, Long, Mic) in in/s and psi -- i.e. true per-sample ground truth,
exactly what the BW ASCII exports give us for series-3.  The leading 2-column
rows are the report header (PPV, sample rate, geo range, ...).

Usage:
    python scratch/verify_thor_against_csv.py [--root DIR] [--lsb FLOAT]
                                              [--limit N] [--kind idfw|idfh|both]
"""
from __future__ import annotations

import argparse
import csv
import os
import statistics
import sys
from collections import Counter, defaultdict

sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))

from micromate import idf_file as M

DEFAULT_ROOT = "/home/serversdown/thor-watcher/example-data"
GEO = ("Tran", "Vert", "Long")


def parse_export(path):
    """Return (header_dict, sample_rows) from a Thor CSV export."""
    hdr, rows = {}, []
    with open(path, newline="", encoding="utf-8", errors="replace") as fh:
        for rec in csv.reader(fh):
            if len(rec) == 2:
                hdr[rec[0].strip()] = rec[1].strip()
            elif len(rec) >= 3:
                try:
                    rows.append([float(x) for x in rec])
                except ValueError:
                    pass
    return hdr, rows


def index_corpus(root):
    """Map BASENAME.IDFW -> (binary_path, csv_path) for every paired file."""
    exports, binaries = {}, {}
    for dirpath, _dirs, files in os.walk(root):
        for name in files:
            up = name.upper()
            full = os.path.join(dirpath, name)
            if up.endswith(".IDFW.CSV") or up.endswith(".IDFH.CSV"):
                exports.setdefault(name[:-4].upper(), full)
            elif up.endswith(".IDFW") or up.endswith(".IDFH"):
                binaries.setdefault(up, full)
    return {k: (binaries[k], exports[k]) for k in binaries.keys() & exports.keys()}


def hdr_float(hdr, key):
    raw = hdr.get(key)
    if not raw:
        return None
    try:
        return float(raw.split()[0])
    except (ValueError, IndexError):
        return None


def verify_waveform(binpath, csvpath, lsb):
    """Compare one IDFW against its export.  Returns a result dict."""
    out = {"file": os.path.basename(binpath), "status": "ok"}
    try:
        res = M.read_idf_file(binpath)
    except NotImplementedError:
        out["status"] = "not-thor"
        return out
    except Exception as exc:  # noqa: BLE001 - harness reports, never raises
        out["status"] = "decode-error"
        out["detail"] = f"{type(exc).__name__}: {exc}"
        return out

    hdr, rows = parse_export(csvpath)
    if not rows:
        out["status"] = "no-gt-samples"
        return out

    gt = {ch: [r[i] for r in rows] for i, ch in enumerate(GEO)}
    out["gt_len"] = len(rows)
    out["geo_range"] = hdr.get("GeoRange")

    exact = total = 0
    lens, chan_status = {}, {}
    ppv_err = {}
    for ch in GEO:
        arr = res.samples.get(ch, [])
        ref = gt[ch]
        lens[ch] = len(arr)
        if len(arr) != len(ref):
            chan_status[ch] = "length"
            continue
        if not arr:
            chan_status[ch] = "empty"
            continue
        hits = sum(1 for c, v in zip(arr, ref) if abs(c * lsb - v) < 5e-5)
        exact += hits
        total += len(arr)
        chan_status[ch] = "exact" if hits == len(arr) else "value"
        gp = hdr_float(hdr, f"{ch}PPV")
        if gp:
            ppv_err[ch] = (max(abs(c) for c in arr) * lsb - gp) / gp

    out["lens"] = lens
    out["chan_status"] = chan_status
    out["exact"] = exact
    out["total"] = total
    out["ppv_err"] = ppv_err
    if all(v == "exact" for v in chan_status.values()):
        out["status"] = "exact"
    elif any(v == "length" for v in chan_status.values()):
        out["status"] = "length-mismatch"
    else:
        out["status"] = "value-mismatch"
    return out


def verify_histogram(binpath, csvpath, lsb):
    out = {"file": os.path.basename(binpath), "status": "ok"}
    try:
        res = M.read_idf_file(binpath)
    except NotImplementedError:
        out["status"] = "not-thor"
        return out
    except Exception as exc:  # noqa: BLE001
        out["status"] = "decode-error"
        out["detail"] = f"{type(exc).__name__}: {exc}"
        return out
    hdr, _rows = parse_export(csvpath)
    out["n_intervals"] = len(res.intervals or [])
    errs = {}
    for ch, attr in (("Tran", "transverse_ips"), ("Vert", "vertical_ips"),
                     ("Long", "longitudinal_ips")):
        gp = hdr_float(hdr, f"{ch}PPV")
        dv = getattr(res.event.peaks, attr, None)
        if gp and dv:
            errs[ch] = (dv - gp) / gp
    out["ppv_err"] = errs
    out["status"] = "peaks" if errs else "no-gt-peaks"
    return out


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--root", default=DEFAULT_ROOT)
    ap.add_argument("--lsb", type=float, default=M._GEO_LSB_IPS)
    ap.add_argument("--limit", type=int, default=0)
    ap.add_argument("--kind", choices=("idfw", "idfh", "both"), default="both")
    ap.add_argument("--show", type=int, default=15, help="worst-N detail rows")
    args = ap.parse_args()

    pairs = index_corpus(args.root)
    keys = sorted(pairs)
    if args.kind != "both":
        keys = [k for k in keys if k.endswith(args.kind.upper())]
    if args.limit:
        keys = keys[: args.limit]

    print(f"root: {args.root}")
    print(f"geo LSB under test: {args.lsb!r} in/s per count")
    print(f"paired files: {len(keys)}\n")

    wf, hg = [], []
    for k in keys:
        binpath, csvpath = pairs[k]
        if k.endswith(".IDFW"):
            wf.append(verify_waveform(binpath, csvpath, args.lsb))
        else:
            hg.append(verify_histogram(binpath, csvpath, args.lsb))

    if wf:
        st = Counter(r["status"] for r in wf)
        ex = sum(r.get("exact", 0) for r in wf)
        tot = sum(r.get("total", 0) for r in wf)
        print("=" * 68)
        print(f"WAVEFORM (IDFW): {len(wf)} files")
        for s, n in st.most_common():
            print(f"    {s:16} {n:5d}  ({100*n/len(wf):5.1f}%)")
        if tot:
            print(f"  per-sample exact: {ex}/{tot} = {100*ex/tot:.3f}%")
        errs = [e for r in wf for e in r.get("ppv_err", {}).values()]
        if errs:
            print(f"  PPV rel-error: median {statistics.median(errs):+.4%}  "
                  f"mean {statistics.mean(errs):+.4%}  "
                  f"max|.| {max(abs(e) for e in errs):.4%}")
        bad = [r for r in wf if r["status"] not in ("exact",)]
        if bad:
            print(f"\n  worst {min(args.show, len(bad))} of {len(bad)} non-exact:")
            for r in bad[: args.show]:
                print(f"    {r['file']:42} {r['status']:16} "
                      f"lens={r.get('lens')} gt={r.get('gt_len')} "
                      f"{r.get('detail','')}")

    if hg:
        st = Counter(r["status"] for r in hg)
        print("=" * 68)
        print(f"HISTOGRAM (IDFH): {len(hg)} files")
        for s, n in st.most_common():
            print(f"    {s:16} {n:5d}  ({100*n/len(hg):5.1f}%)")
        errs = [e for r in hg for e in r.get("ppv_err", {}).values()]
        if errs:
            print(f"  PPV rel-error: median {statistics.median(errs):+.4%}  "
                  f"mean {statistics.mean(errs):+.4%}  "
                  f"max|.| {max(abs(e) for e in errs):.4%}")
            within = lambda t: 100*sum(1 for e in errs if abs(e) <= t)/len(errs)
            print(f"  within 0.5%: {within(0.005):.1f}%   "
                  f"within 2%: {within(0.02):.1f}%   within 5%: {within(0.05):.1f}%")
    return 0


if __name__ == "__main__":
    raise SystemExit(main())