2026-09-07 15:57:38 -04:00
2 changed files with 387 additions and 0 deletions
@@ -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. |
@@ -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()