diff --git a/docs/offset_investigation.md b/docs/offset_investigation.md index 20da227..537691c 100644 --- a/docs/offset_investigation.md +++ b/docs/offset_investigation.md @@ -58,6 +58,14 @@ Companion material: sensor-check failures. Do not try to use it as a screen. - **Cause is still unsettled.** Instantel's autozero fixes the minority of cases; the rest are hardware. We cannot yet tell which is which remotely. +- **The histogram corpus (63,535 files, 9.7x the waveforms) is now scanned too** — + see §8b. It independently confirms BE18438 and BE9558 with a clean 2.5x + separation, but detects only **2 of the 5** confirmed units, cannot attribute a + channel, and resolves time to ~a month. **A negative histogram result is not + evidence of health** — DC leakage into the interval peak varies 45x between units. +- **`offset_scan3.py` has a label defect** (§8b): its spread gate discards 18.8% of + high-|pre| rows onto units currently counted as clean. Re-cut before quoting any + precision number again. - **Best open lead:** `SUB 0x0E` (channel sensor data, 8 channels × 10 bytes, unimplemented) may carry the very numbers Instantel says to check against **2027–2069**. Untested. @@ -495,6 +503,165 @@ doubled two reported figures before it was caught. --- +## 8b. The histogram corpus — the other 90% of the archive (2026-09-04) + +Every result above §8 comes from **waveform** files. `offset_scan3.py` filters on +`\.[A-Za-z0-9]{2}0[Ww]$`, so the corpus it scanned is 6,577 unique binaries. The +archive also holds **63,535 unique histograms** — 9.7x more files — which the +pre-trigger method cannot touch, because a histogram carries no samples: only a +per-interval, per-channel peak and half-period. + +`scratch/offset_hist_scan.py` scans them. **63,505 of 63,535 decoded (99.95%), +43 units, 77.9M intervals.** Two of the 45 units have no histograms at all. +Output: `/home/serversdown/dl2-archive/offset_hist.csv` (190,515 channel-rows). + +### The premise, and how far it actually holds + +A histogram file is hours of continuous monitoring, so most of its intervals are +definitionally quiet, and a channel parked off zero cannot report a peak below +its own displacement. The signal is real — two within-unit contrasts, siblings +unmoved in both: + +| unit | channel | in-episode floor | outside | waveform \|pre\| same window | +|---|---|---|---|---| +| BE18438 | Vert | 0.0350 | 0.0050 | +0.18 .. +0.37 | +| BE12599 | Tran | 0.0250 | 0.0050 | +0.03 .. +0.49 | + +But the **leakage from a waveform pedestal into the histogram floor is bimodal, +not merely partial**: measured ratio ~0.9 on BE18438 Vert, ~0.7 on BE9558, +**~0.02 on BE12599** — two orders of magnitude on one instrument. The device +evidently measures each interval peak against a running baseline, and how much +DC survives that varies per unit. **Consequence: a negative histogram result +carries almost no information.** Do not read "clean in the histograms" as clean. + +### The detector that survived + + dmin(file, ch) = min[ch] - min over the other two geo channels, SAME file + gates (both hard): n_intervals >= 60 AND mic_p5 <= 5 raw counts + day statistic: median of dmin over that day's qualifying files + flag day at dmin >= 0.020 in/s (4 A/D counts) + episode at >= 3 CONSECUTIVE observed days + +**Result: BE18438|Vert, BE9558|Tran, BE9558|Long.** Threshold-insensitive — +the journal's own test for a real signal against a tuned one — and this is the +first operating point in the investigation that passes it cleanly. The identical +answer holds across: statistic `min` or `p5`; length gate 10/30/60/120/300; mic +gate 3/5/8/10; threshold 0.015–0.035 (a 2.3x span); persistence K = 2,3,4,5,7. + +Separation, ranked by highest floor sustained over 3 consecutive gated days +across all 135 unit-channels: + +| unit-channel | best3 | +|---|---| +| BE18438 Vert | 0.1650 | +| BE9558 Long | 0.0350 | +| BE9558 Tran | 0.0250 | +| *(2.5x gap)* | | +| BE7145 Tran | 0.0100 | +| entire rest of fleet | <= 0.0050 (one quantisation count) | + +Day-level false alarm: **37 of 99,432 gated unit-channel-days = 0.037%.** + +### What it does NOT do — read this before trusting it + +- **It finds 2 of the 5 confirmed units, not 5.** The site-quiet gate is what + makes it work and it is also what costs BE11529 and BE12599. BE11529's + four-day single-axis ramp (Tran 0.025 -> 0.055, both siblings pinned at 0.005) + is the most offset-shaped thing in the corpus outside the two detections, and + the gate discards it. +- **The positive class is two units.** Every threshold here is fitted to + BE18438 and BE9558, which contribute 22 of the 37 flagged days in the entire + corpus. No cross-validation is possible at n=2. +- **Per-channel attribution is NOT established.** Rotating the three geo channel + labels within each file — preserving every value, file and day, destroying + only channel identity — reproduces the episode *count* with p = 0.769 and the + label agreement at p = 0.038–0.077. Report a **unit and a window**; do not + name a geophone axis on the strength of this detector alone. +- **Timing resolution is ~1 month, not ~1 day.** A 30-day label shift still + scores 2 of 9 episode hits; the signal dies only past ~60 days. The day-level + series look far crisper than they are. +- **Ground truth here is a sibling detector, not a service record.** Agreement + between the two corpora is corroboration of a shared method. Nothing in this + section has been checked against an actual repair, calibration or RMA. + +### Dead ends — keep these dead + +- **Absolute floor (min / p1 / p5 / p10 / p25, thresholded alone) — RETIRED.** + Not fleet-comparable and mostly not about the channel. Scoring each cell using + *only the other two channels* — a statistic containing zero information about + the suspect channel — reaches AUC 0.746 against the same labels, versus 0.872 + for the absolute floor itself. **66% of its apparent discrimination is "that + day was noisy at that site."** Interval size alone moves its p99 7x (0.0350 at + 1 min vs 0.0050 at 2 s). And of all files with any channel above 0.025, 56.5% + have **all three** channels above it — common-mode, i.e. the wrong physics. +- **Zero-fraction — STRUCTURALLY IMPOSSIBLE, not merely weak.** The device never + reports a zero histogram interval peak. The value is a max over hundreds of + samples of a channel that always carries at least 1 count of noise, so it is + clamped at 1 A/D count (0.005 in/s). There is no zero to count. +- **Interval size, sample rate, geo range, firmware — refuted as confounds for + the differential.** All four are *file-level scalars*: they move all three geo + channels together, so they cannot produce a single-channel lift and the + within-file differential is immune to them by construction. Geo range is + identical across the three geo channels in **63,535 of 63,535** binaries. + (Interval size remains fatal to the *absolute*-floor version, above.) + +### Two findings that are independent of the histogram detector + +**1. `offset_scan3.py`'s `spread <= 0.02` gate is discarding real signal.** +It rejects **113 of the 600 channel-rows with |pre| >= 0.025 (18.8%)**, and the +rejections are not random — 92 of them fall across 41 unit-channels currently +labelled NEGATIVE. Four would become sustained positives under an +amplitude-only >=3-consecutive rule: **BE12599|Long (run of 8), BE18003|Vert +(4), BE10895|Vert (3), BE12844|Tran (3).** Until this is re-cut, the fleet label +is **three-state — POSITIVE / NEGATIVE / SPREAD-REJECTED(unknown)** — and the +third state should be excluded from both TP and FP counts rather than silently +scored as healthy. Every precision figure computed against the two-state label, +in this section and in §3, is affected. + +**2. The waveform corpus sees ~7% of the days a unit was deployed.** 2,627 +(unit, day) observations against the histogram corpus's 35,105 — 13.4x — with a +per-unit median ratio of 0.070. BE12599, a confirmed unit, is waveform-observed +on 39 of its 1,666 histogram-observed days (**2.3%**). Any statement of the form +"the fault was absent before date X" that rests on waveform coverage alone is +much weaker than its event count suggests. + +### BE10895 — reclassified (see also §4) + +Previously dismissed as a transient. The histogram record shows its **Vert** +quiet-minute floor at 0.005 on 62/62 qualifying files from 2023-07-07, then +0.010–0.015 on 48/58 files from 2023-08-03 to 08-27, while Tran moves on 2/58 +and Long on 9/58 and the site mic floor never leaves 1–3 counts. Independently, +**42 of its 85 waveform events (49.4%) are single-axis-dominant** — one geo peak +>= 10x both siblings and >= 0.05 in/s — the **highest rate in the 45-unit +fleet** (BE13117 36.1%, BE18438 29.4%), and **100% of it on Vert**. Vert +excursions of 0.1–1.5 in/s with Tran/Long at 0.005–0.035 are not ground motion. + +This is a genuine Vert-channel hardware fault, but **not the classic pedestal** — +the differential is only one A/D count. Caveat: its entire histogram record is a +single 52-day deployment ending 2023-08-27, so nothing says whether it +persisted, was serviced, or resolved. + +The other six marginal units — BE11007, BE17354, BE18004, BE18104, BE9557, +BE18003 — are **clean**. All seven cap at +0.005 to +0.007 (one A/D count) +lifetime under the quiet-site gate, against +0.175 for BE18438 Vert and +0.062 +for BE9558 Long. Three individual waveform flags fall in windows with **zero** +histogram coverage and are NO-DATA, not clean: BE18004|Tran 2024-10-16, +BE9557|Tran 2021-06-28, BE9557|Vert 2025-06-12. + +### Still open in this section + +- **The 11 thin-coverage units were not screened** (BE10202, BE11462, BE13779, + BE15760, BE15957, BE16754, BE16758, BE8081, BE8626, BE9229, BE9887 — each + under 20 waveform events, several with hundreds of histograms). This is the + population most likely to hold a previously unknown offset, and it is the one + slice of the plan that did not run. BE11462 was incidentally scored clean by + the full-archive pass; BE10202 has no histogram files at all. +- **No completeness audit was run** over the above. +- Re-cutting the ground truth three-state (finding 1) and re-scoring everything + against it. + +--- + ## 9. Chronology | date | event | @@ -512,3 +679,8 @@ doubled two reported figures before it was caught. | 2026-08-28 | Bimodality established; sensor check proven **blind** to offsets; `SUB 0x0E` identified as the best open lead. | | 2026-08-28 | **v1 detector retracted.** Brian challenged the "come and go" finding against field experience. Two flaws found: dominant-axis-only scoring and mean-instead-of-median. Corrected detector shows persistent pedestals on **8 of 45 units**, and the gaps are service windows. | | 2026-08-28 | **Detector v3 (Brian's method):** pre-trigger floor + pre/mid/end consistency. Healthy channels proven to sit at 0.000 +/-1 unit (94.5%), confirming no decoder zero-point bias. Final: **5 of 45 units (11%)**, threshold-insensitive. | +| 2026-09-04 | **Histogram corpus scanned** — 63,505 of 63,535 files, 43 units, 77.9M intervals (9.7x the waveform corpus). `scratch/offset_hist_scan.py`. | +| 2026-09-04 | Absolute-floor statistic **retired**: 66% of its discrimination is a day/site confound (other-channels-only AUC 0.746 vs 0.872). Zero-fraction shown **structurally impossible** — the device clamps every interval peak at >= 1 count. | +| 2026-09-04 | Site-quiet-gated cross-channel differential established: **BE18438 Vert, BE9558 Tran+Long**, threshold-insensitive over a 2.3x span. Finds only **2 of the 5** confirmed units — leakage into the histogram floor is bimodal (0.9 to 0.02), so a negative result carries almost no information. Per-channel attribution **not** established (channel-scramble p = 0.769). | +| 2026-09-04 | **BE10895 reclassified** from transient to a genuine Vert fault of a different subtype — 49.4% single-axis-dominant events, the highest in the fleet, 100% on Vert. The other six marginal units are clean. | +| 2026-09-04 | **Defect found in `offset_scan3.py`**: its `spread <= 0.02` gate discards 18.8% of rows with \|pre\| >= 0.025, concentrated on 41 negative unit-channels; 4 would be sustained positives without it. The fleet label is three-state, not two. | diff --git a/scratch/offset_hist_scan.py b/scratch/offset_hist_scan.py new file mode 100644 index 0000000..1a2b9a1 --- /dev/null +++ b/scratch/offset_hist_scan.py @@ -0,0 +1,215 @@ +#!/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" + + +def serial_of(name: str) -> str: + """`P036L318.C80H` -> `BE14036`. See CLAUDE.md, serial encoding.""" + m = _STEM.match(name) + if not m: + return "?" + return f"BE{(ord(m.group(1)) - ord('B')) * 1000 + int(m.group(2))}" + + +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), + "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()