seismo-relay/scratch/offset_hist_scan.py

#!/usr/bin/env python3
"""Offset detector — HISTOGRAM corpus (the other 90% of the archive).

`offset_scan3.py` measures the pre-trigger floor in *waveform* samples.  That
covers 6,577 of the archive's 70,112 unique series-3 files; the remaining
63,535 are **histograms**, which carry no samples — only a per-interval,
per-channel peak + half-period.  So the pre-trigger method cannot run on them.

The histogram analogue of "the resting floor" is the **low percentile of the
per-interval peaks**.  A histogram file is typically hours of continuous
monitoring, so the great majority of its intervals are definitionally quiet;
the bottom of that distribution is what the channel reads when nothing is
happening.  A healthy channel bottoms out at 0.000-0.005 in/s.  A channel
parked off zero cannot report a peak below its own displacement, so its floor
is pinned up.

⚠ The DC leakage into the histogram peak is PARTIAL.  Measured within-unit
against episodes already established from the waveform scan:

    BE18438  Vert  in-episode 0.0350 vs 0.0050 outside  (waveform pre = +0.18..+0.37)
    BE12599  Tran  in-episode 0.0250 vs 0.0050 outside  (waveform pre = +0.03..+0.49)

so the device's per-interval peak is evidently measured against a running /
AC-coupled baseline that removes most, but not all, of the DC.  The residual
is real and channel-specific, but the margin is ~5 quantisation counts rather
than the ~70 the waveform detector enjoys.  Do not carry the waveform
detector's 0.025 in/s floor across unexamined — calibrate on the CSV.

Because the absolute floor also moves with site noise (traffic, wind, a
generator), the statistic that matters most is the **cross-channel
differential**: a channel's floor minus the quietest of the other two geo
channels in the same file.  Site noise lifts all three together and cancels;
a DC offset lifts one.

This script does not decide anything.  It emits every candidate statistic per
(file, channel) so thresholds can be calibrated against the waveform-derived
ground truth in `offset_v3.csv` rather than guessed.

Usage:
    python scratch/offset_hist_scan.py --dir /home/serversdown/dl2-archive/files \
        --out /home/serversdown/dl2-archive/offset_hist.csv --jobs 4
"""
from __future__ import annotations

import argparse
import csv
import datetime
import logging
import re
import statistics
import sys
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  # noqa: E402

GEO = ("Tran", "Vert", "Long")
K = 10.0 / 32000.0          # ADC count -> in/s (see CLAUDE.md: full scale 32000)
_HIST = re.compile(r"\.[A-Za-z0-9]{2}0[Hh]$")
_STEM = re.compile(r"^([B-Z])(\d{3})")
_B36 = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"


_SERIAL_RE = re.compile(rb"\b([A-Z]{2}\d{3,6})\b")


def serial_of(name: str, path=None) -> str:
    """Real serial for a BW file.

    The filename encodes only the NUMBER: `<letter><3 digits>` where
    letter = chr(ord('B') + serial // 1000).  The two-letter family prefix
    ("BE", "BA", ...) is **not** in the filename, so it must be read out of
    the file body.  Four units in the DL2 archive are BA, not BE — assuming
    "BE" mislabels BA9229, BA10060, BA10895 and BA15957.
    """
    m = _STEM.match(name)
    if not m:
        return "?"
    num = (ord(m.group(1)) - ord("B")) * 1000 + int(m.group(2))
    if path is not None:
        try:
            for s in _SERIAL_RE.findall(Path(path).read_bytes()):
                s = s.decode()
                if s[2:].lstrip("0") == str(num):
                    return s
        except Exception:
            pass
    return f"BE{num}"          # last-resort fallback; prefix unverified


def stem_time(name: str):
    """Decode the filename's base-36 timestamp.  Epoch 1985-01-01, 1296 s/tick.

    Preferred over the file's own footer timestamp only because it costs
    nothing; the caller falls back to the decoded event when this fails.
    """
    try:
        base, ext = name.rsplit(".", 1)
        n = 0
        for c in base[4:8].upper():
            n = n * 36 + _B36.index(c)
        ab = _B36.index(ext[0].upper()) * 36 + _B36.index(ext[1].upper())
        return datetime.datetime(1985, 1, 1) + datetime.timedelta(seconds=n * 1296 + ab)
    except Exception:
        return None


def _pct(sorted_vals, q):
    """Nearest-rank percentile on an already-sorted list."""
    if not sorted_vals:
        return None
    i = min(len(sorted_vals) - 1, max(0, int(len(sorted_vals) * q / 100.0)))
    return sorted_vals[i]


def scan(path_str: str):
    logging.disable(logging.WARNING)   # per-worker: the codec warns on undecodables
    p = Path(path_str)
    try:
        ev = read_blastware_file(p)
    except Exception:
        return None
    s = ev.raw_samples or {}
    if not any(s.get(c) for c in GEO):
        return None

    ts = stem_time(p.name) or ev.timestamp
    stamp = ""
    if ts is not None:
        stamp = (f"{ts.year:04d}-{ts.month:02d}-{ts.day:02d}T"
                 f"{ts.hour:02d}:{ts.minute:02d}:{ts.second:02d}")

    # Per-channel floor candidates, in in/s.
    stats = {}
    for ch in GEO:
        v = sorted(s.get(ch) or [])
        if not v:
            continue
        stats[ch] = {
            "n": len(v),
            "min": v[0] * K,
            "p1": _pct(v, 1) * K,
            "p5": _pct(v, 5) * K,
            "p10": _pct(v, 10) * K,
            "p25": _pct(v, 25) * K,
            "med": statistics.median(v) * K,
            "peak": v[-1] * K,
            "zeros": sum(1 for x in v if x == 0) / len(v),
        }
    if len(stats) < 2:      # need at least one sibling channel for the differential
        return None

    # Mic floor as a site-noise proxy (raw counts; the dB conversion is not
    # needed — only its relative movement matters here).
    mic = sorted(s.get("MicL") or [])
    mic_p5 = _pct(mic, 5) if mic else ""

    rows = []
    for ch, st in stats.items():
        others = [stats[o]["p5"] for o in stats if o != ch]
        rows.append({
            "serial": serial_of(p.name, p),
            "timestamp": stamp,
            "filename": p.name,
            "channel": ch,
            "n_intervals": st["n"],
            "min": round(st["min"], 4),
            "p1": round(st["p1"], 4),
            "p5": round(st["p5"], 4),
            "p10": round(st["p10"], 4),
            "p25": round(st["p25"], 4),
            "median": round(st["med"], 4),
            "peak": round(st["peak"], 4),
            "frac_zero": round(st["zeros"], 4),
            # the site-noise-cancelling statistic: this channel's floor above
            # the quietest sibling geo channel in the same file
            "diff_p5": round(st["p5"] - min(others), 4),
            "mic_p5": mic_p5,
        })
    return rows


COLS = ["serial", "timestamp", "filename", "channel", "n_intervals",
        "min", "p1", "p5", "p10", "p25", "median", "peak", "frac_zero",
        "diff_p5", "mic_p5"]


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--dir", required=True)
    ap.add_argument("--out", required=True)
    ap.add_argument("--jobs", type=int, default=4)
    ap.add_argument("--limit", type=int, default=0, help="stop after N files (smoke test)")
    a = ap.parse_args()

    # Dedupe by basename — the DL2 export keeps a byte-identical `Sent/`
    # mirror of its root, which doubled two figures before it was caught.
    seen, files = set(), []
    for q in sorted(Path(a.dir).rglob("*")):
        if q.is_file() and _HIST.search(q.name) and q.name not in seen:
            seen.add(q.name)
            files.append(str(q))
    if a.limit:
        files = files[:a.limit]
    print(f"unique histogram binaries: {len(files)}", flush=True)

    rows, undecodable = [], 0
    with ProcessPoolExecutor(max_workers=a.jobs) as ex:
        futs = [ex.submit(scan, f) for f in files]
        for i, fut in enumerate(as_completed(futs), 1):
            r = fut.result()
            if r:
                rows.extend(r)
            else:
                undecodable += 1
            if i % 5000 == 0:
                print(f"  {i}/{len(files)}", flush=True)

    with open(a.out, "w", newline="") as fh:
        w = csv.DictWriter(fh, fieldnames=COLS)
        w.writeheader()
        w.writerows(rows)

    files_ok = len({r["filename"] for r in rows})
    units = len({r["serial"] for r in rows})
    ivals = sum(r["n_intervals"] for r in rows) // 3
    print(f"\ndecoded {files_ok}/{len(files)} files "
          f"({undecodable} undecodable), {units} units, ~{ivals/1e6:.1f}M intervals")
    print(f"wrote {a.out}  ({len(rows)} channel-rows)")


if __name__ == "__main__":
    main()