From 14e997b20cd96ab2b2e387cbb8c765370136856d Mon Sep 17 00:00:00 2001 From: serversdown Date: Fri, 28 Aug 2026 07:29:14 +0000 Subject: [PATCH] fix(histogram): partial final block no longer discards the correct stride MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit detect_multi_interval_stride() confirmed a candidate stride on a third block header whenever the body was long enough to contain one. But a body can exceed two strides and still hold only two real blocks: a partial final block leaves trailing padding. BE18193 T193L0XM.CI0H — 51 intervals at 2 s, i.e. one full 30-interval block plus a 21-interval remainder in a 2787-byte body — had every decisive check pass at stride 612 (header at 0, header at 612, block counter 256 -> 257) and was then rejected for the absent third header at 1224. It decoded to nothing. A missing third header now means end-of-stream rather than disqualification. The block-counter check is untouched — that is the test that prevents the false positives which once handed 9,082 standard-block files to the multi-interval walker. Found by running the full DL2 archive against its preserved Blastware ASCII exports (14,340 paired files, 11x the previous ground-truth corpus). Measured over 127,035 archive histogram binaries: recovered 8 files (strides 92, 252, 612; BE18193, BE18191, BE9557, BE9440) regressed 0 files Full-corpus verification: 14,337 -> 14,338 exact of 14,338 decodable pairs (the 2 excluded are series-4 IDF, a different codec). Also adds scratch/verify_against_ascii.py (per-sample decoder verification against BW exports, with a saturation carve-out — BW clamps clipped events to the range max while the decoder reports true counts) and scratch/offset_scan.py. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01HgTe8CamXAHcAmaQ6QNcog --- CHANGELOG.md | 21 +++ minimateplus/histogram_codec.py | 16 ++- scratch/offset_scan.py | 164 ++++++++++++++++++++++++ scratch/verify_against_ascii.py | 221 ++++++++++++++++++++++++++++++++ tests/test_histogram_codec.py | 35 +++++ 5 files changed, 453 insertions(+), 4 deletions(-) create mode 100644 scratch/offset_scan.py create mode 100644 scratch/verify_against_ascii.py diff --git a/CHANGELOG.md b/CHANGELOG.md index afa663f..1ab9bf9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,27 @@ All notable changes to seismo-relay are documented here. ## [Unreleased] ### Fixed +- **Sub-minute histograms with a partial final block decoded to nothing** + (`histogram_codec.detect_multi_interval_stride`). The stride search confirmed + itself on a third block header whenever the body was long enough to hold one — + but a body can exceed two strides and still contain only two real blocks, because + a *partial* final block leaves trailing padding. BE18193 `T193L0XM.CI0H` (51 + intervals at 2 s = one full 30-interval block plus a 21-interval remainder, in a + 2787-byte body) therefore had its correct stride of 612 discarded and produced an + empty decode. A missing third header now means end-of-stream rather than + disqualification; the block-counter check, which is what actually prevents the + false positives that once mis-dispatched 9,082 files, is unchanged. + + Found by decoding the full DL2 archive against its preserved Blastware ASCII + exports. Across 127,035 histogram binaries the fix recovers **8 files** (strides + 92, 252 and 612; units BE18193, BE18191, BE9557, BE9440) with **zero** files + regressed. Verification over all 14,340 archive pairs goes 14,337 → 14,338 exact, + the only remainder being two series-4 IDF files that belong to a different codec. + + ⚠ Prod stores hold `.h5` files generated before this fix. The 8 affected events + stay empty until `backfill_sidecars.py` is re-run — not urgent at 8 files, and + worth folding into the next backfill rather than doing one for this alone. + - **Histogram/waveform twin matching is now interval-based** (`find_twins`). A real trigger is recorded twice — as a triggered waveform (stamped at the trigger instant) and inside the scheduled histogram whose interval contains it (stamped at the 7am/7pm diff --git a/minimateplus/histogram_codec.py b/minimateplus/histogram_codec.py index d60f853..a24d281 100644 --- a/minimateplus/histogram_codec.py +++ b/minimateplus/histogram_codec.py @@ -413,10 +413,18 @@ def detect_multi_interval_stride(body: bytes) -> Optional[int]: if (_ctr(stride) - _ctr(0)) & 0xFFFF != 1: continue - # confirm on a third block when the body is long enough - if 2 * stride + _MULTI_HEADER_LEN <= len(body): - if not _is_multi_header(body, 2 * stride): - continue + # Confirm on a third block WHEN ONE IS ACTUALLY PRESENT. A body can + # be longer than two strides and still hold only two real blocks: a + # final *partial* block leaves trailing padding. E.g. 51 intervals at + # 2 s = one full 30-interval block + a 21-interval remainder, in a + # 2787-byte body — long enough to demand a third header at 1224 that + # does not exist. Requiring it unconditionally threw away the correct + # stride and the file decoded to nothing (BE18193 T193L0XM.CI0H). + # The block-counter check above is the decisive anti-false-positive + # test; this one is corroboration, so a missing third header means + # end-of-stream, not disqualification. + if (2 * stride + _MULTI_HEADER_LEN <= len(body) + and _is_multi_header(body, 2 * stride)): if (_ctr(2 * stride) - _ctr(stride)) & 0xFFFF != 1: continue return stride diff --git a/scratch/offset_scan.py b/scratch/offset_scan.py new file mode 100644 index 0000000..842bdcf --- /dev/null +++ b/scratch/offset_scan.py @@ -0,0 +1,164 @@ +#!/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 [--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() + + files = [p for p in Path(a.dir).rglob("*") if p.is_file() and _WAVE_RE.search(p.name)] + files.sort() + 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() diff --git a/scratch/verify_against_ascii.py b/scratch/verify_against_ascii.py new file mode 100644 index 0000000..da14211 --- /dev/null +++ b/scratch/verify_against_ascii.py @@ -0,0 +1,221 @@ +#!/usr/bin/env python3 +"""Verify the series-3 decoder against preserved Blastware ASCII exports. + +Pairs each `__ASCII.TXT` with its binary `.`, 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 [--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() diff --git a/tests/test_histogram_codec.py b/tests/test_histogram_codec.py index 558a0c9..1d79f0f 100644 --- a/tests/test_histogram_codec.py +++ b/tests/test_histogram_codec.py @@ -643,3 +643,38 @@ def test_multi_interval_matches_blastware_ascii_exactly(): assert hz is None elif not cell.startswith("<"): assert hz is not None and abs(hz - float(cell)) <= max(0.55, float(cell) * 0.02) + + +def test_partial_final_block_is_not_disqualified_by_missing_third_header(): + """A body can exceed two strides yet hold only two real blocks. + + Regression for BE18193 `T193L0XM.CI0H` — 51 intervals at 2 s = one full + 30-interval block plus a 21-interval remainder, in a body long enough to + demand a third block header at ``2 * stride`` that does not exist. The + third-block confirmation used to be mandatory whenever the body was long + enough, so the correct stride was discarded and the file decoded to + nothing. A missing third header means end-of-stream, not disqualification; + the block-counter check is the decisive anti-false-positive test. + """ + full = [(1, 1, 2, 2, 3, 3, 4, 4)] * 30 + partial = [(5, 5, 6, 6, 7, 7, 8, 8)] * 21 + body = (_mk_multi_block(full, ctr=256) + + _mk_multi_block(partial, ctr=257) + + b"\xff" * 700) # trailing padding past 2 * stride + stride = 12 + 20 * 30 + assert 2 * stride + 6 <= len(body), "padding must reach past two strides" + # the whole point: a third header is absent, and that must not disqualify + assert detect_multi_interval_stride(body) == stride + recs = walk_multi_interval_blocks(body) + assert len(recs) == 51 + assert recs[0]["t_peak"] == 1 + assert recs[-1]["t_peak"] == 5 + + +def test_third_block_still_rejects_a_mismatched_counter(): + """The corroboration must still bite when a third block IS present.""" + ivs = [(1, 1, 2, 2, 3, 3, 4, 4)] * 4 + body = (_mk_multi_block(ivs, ctr=256) + + _mk_multi_block(ivs, ctr=257) + + _mk_multi_block(ivs, ctr=999)) # counter jumps — not consecutive + assert detect_multi_interval_stride(body) != 12 + 20 * 4