seismo-relay/scratch/verify_against_ascii.py

#!/usr/bin/env python3
"""Verify the series-3 decoder against preserved Blastware ASCII exports.

Pairs each `<stem>_<ext>_ASCII.TXT` with its binary `<stem>.<ext>`, decodes the
binary with the production codec, and compares against BW's own export:

  waveform  — per-channel sample counts, then every sample value
  histogram — interval count, then every per-interval channel peak

ADC counts convert as  ips = counts * geo_range_ips / 32000  (1 decoder unit =
16 counts = 0.005 in/s at the 10 in/s range; see CLAUDE.md).

Usage:
  python scratch/verify_against_ascii.py --dir <path> [--limit N] [--jobs N]
                                         [--out results.json] [--kind w|h|all]
"""
from __future__ import annotations

import argparse, json, re, sys, traceback
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_file, parse_report

GEO = ("Tran", "Vert", "Long")
_GEO_FS_COUNTS = 32000.0
_ASCII_SUFFIX_RE = re.compile(r"_ASCII\.TXT$", re.IGNORECASE)


def binary_for(ascii_path: Path) -> Path:
    """H907KXOW_WC0H_ASCII.TXT -> H907KXOW.WC0H"""
    stem = _ASCII_SUFFIX_RE.sub("", ascii_path.name)
    if "_" not in stem:
        return ascii_path.with_name(stem)
    head, _, ext = stem.rpartition("_")
    return ascii_path.with_name(f"{head}.{ext}")


def counts_to_ips(counts, geo_range_ips):
    r = geo_range_ips if geo_range_ips else 10.0
    return counts * r / _GEO_FS_COUNTS


def agrees(got, exp, geo_range_ips, tol=0.0006):
    """True when decoded `got` matches BW's exported `exp`.

    Saturation carve-out: when an event clips, BW clamps its export to the
    channel's range maximum (and writes OORANGE for the summary PPV), while
    the decoder faithfully reproduces raw counts that can sit a decoder unit
    or two past nominal full scale (32016 counts observed = 10.005 in/s on
    the 10 in/s range).  Same sign and both at/above the ceiling is agreement,
    not a decode error.
    """
    if abs(got - exp) <= tol:
        return True
    r = geo_range_ips if geo_range_ips else 10.0
    if abs(exp) >= r - tol and abs(got) >= r - tol and (got >= 0) == (exp >= 0):
        return True
    return False


def parse_interval_table(text: str):
    """Histogram interval rows: time, Tpk, Tfq, Vpk, Vfq, Lpk, Lfq, PVS, ..., micdB, micfq"""
    rows = []
    seen_header = False
    for line in text.splitlines():
        if "\t" not in line:
            continue
        cols = [c.strip().strip('"') for c in line.split("\t")]
        cols = [c for c in cols if c != ""]
        if not seen_header:
            if any(c in ("Tran", "Vert", "Long") for c in cols):
                seen_header = True
            continue
        if len(cols) < 7:
            continue
        if not re.match(r"^\d{1,2}:\d{2}:\d{2}$", cols[0]):
            continue
        def num(s):
            try:
                return float(s)
            except ValueError:
                return None
        rows.append({"time": cols[0], "Tran": num(cols[1]),
                     "Vert": num(cols[3]), "Long": num(cols[5])})
    return rows


def check_one(ascii_path_str: str) -> dict:
    ap = Path(ascii_path_str)
    bp = binary_for(ap)
    res = {"ascii": ap.name, "binary": bp.name, "status": "?",
           "kind": None, "detail": ""}
    try:
        if not bp.exists():
            res["status"] = "no_binary"
            return res
        text = ap.read_text(errors="replace")
        rep = parse_report(text, parse_samples=True)
        ev = read_blastware_file(bp)
        gr = rep.geo_range_ips
        res["kind"] = kind = ("histogram"
                              if (rep.event_type or "").lower().startswith(("full histogram", "histogram"))
                              else "waveform")
        samples = ev.raw_samples or {}
        dec_n = {c: len(samples.get(c) or []) for c in GEO}

        if kind == "histogram":
            rows = parse_interval_table(text)
            res["n_ascii"] = len(rows)
            res["n_decoded"] = dec_n["Tran"]
            if not rows:
                res["status"] = "no_ascii_table"
                return res
            if dec_n["Tran"] == 0:
                res["status"] = "decode_empty"
                return res
            if dec_n["Tran"] != len(rows):
                res["status"] = "count_mismatch"
                res["detail"] = f"decoded {dec_n['Tran']} vs ascii {len(rows)}"
                return res
            bad = 0
            worst = 0.0
            for i, row in enumerate(rows):
                for ch in GEO:
                    exp = row[ch]
                    if exp is None:
                        continue
                    got = counts_to_ips(samples[ch][i], gr)
                    if not agrees(got, exp, gr):
                        bad += 1
                        worst = max(worst, abs(got - exp))
            res["worst_abs"] = round(worst, 6)
            res["status"] = "exact" if bad == 0 else "value_mismatch"
            if bad:
                res["detail"] = f"{bad} interval-channel values off"
            return res

        # waveform
        asc = rep.samples or []
        res["n_ascii"] = len(asc)
        res["n_decoded"] = dec_n["Tran"]
        if not asc:
            res["status"] = "no_ascii_table"
            return res
        if dec_n["Tran"] == 0:
            res["status"] = "decode_empty"
            return res
        if len({dec_n[c] for c in GEO}) != 1:
            res["status"] = "channel_len_mismatch"
            res["detail"] = str(dec_n)
            return res
        if dec_n["Tran"] != len(asc):
            res["status"] = "count_mismatch"
            res["detail"] = f"decoded {dec_n['Tran']} vs ascii {len(asc)}"
            return res
        bad = 0
        worst = 0.0
        for i, quad in enumerate(asc):
            for j, ch in enumerate(GEO):
                exp = quad[j]
                got = counts_to_ips(samples[ch][i], gr)
                if not agrees(got, exp, gr):
                    bad += 1
                    worst = max(worst, abs(got - exp))
        res["worst_abs"] = round(worst, 6)
        res["status"] = "exact" if bad == 0 else "value_mismatch"
        if bad:
            res["detail"] = f"{bad} sample values off"
        return res
    except Exception as e:
        res["status"] = "error"
        res["detail"] = f"{type(e).__name__}: {e}"
        return res


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--dir", required=True)
    ap.add_argument("--limit", type=int, default=0)
    ap.add_argument("--jobs", type=int, default=8)
    ap.add_argument("--kind", choices=["w", "h", "all"], default="all")
    ap.add_argument("--out", default=None)
    a = ap.parse_args()

    root = Path(a.dir)
    files = sorted(p for p in root.rglob("*")
                   if p.is_file() and p.name.upper().endswith("_ASCII.TXT"))
    if a.kind != "all":
        want = "0W" if a.kind == "w" else "0H"
        files = [p for p in files
                 if _ASCII_SUFFIX_RE.sub("", p.name).upper().endswith(want)]
    if a.limit:
        files = files[: a.limit]
    print(f"pairs to check: {len(files)}", flush=True)

    out = []
    from collections import Counter
    tally = Counter()
    with ProcessPoolExecutor(max_workers=a.jobs) as ex:
        futs = {ex.submit(check_one, str(p)): p for p in files}
        for n, f in enumerate(as_completed(futs), 1):
            r = f.result()
            out.append(r)
            tally[(r["kind"], r["status"])] += 1
            if n % 500 == 0:
                print(f"  {n}/{len(files)}", flush=True)

    print("\n=== results ===")
    for (kind, status), n in sorted(tally.items(), key=lambda x: -x[1]):
        print(f"  {str(kind):10} {status:22} {n}")
    if a.out:
        Path(a.out).write_text(json.dumps(out, indent=1))
        print(f"\nwrote {a.out}")


if __name__ == "__main__":
    main()