fix(offset): retract the v1 detector — per-channel median, not dominant-axis mean
Brian challenged the v1 finding that offsets "come and go", against field experience that a unit which develops one stays broken until the geophone is replaced. He was right; v1 had two flaws, both of which manufactured false recoveries: 1. It scored only the axis with the largest peak, so a real event on one axis hid a persistent pedestal on another. BE12599 on 2026-08-21 read "clean" because Long had a 1.065 in/s event, while Tran sat at +0.4732 in/s and was never examined. 2. It used the mean, which a real transient perturbs. The median is the resting baseline and a blast does not move it. Same event, Long channel: mean +0.0783 vs median -0.0050. offset_scan2.py flags a CHANNEL when |median| >= 0.025 in/s (5 A/D counts, Instantel's own criterion) and treats >=3 consecutive flagged events as the real signal. No m/p ratio guard is needed — that existed only to compensate for the mean. Corrected results: units with any flagged event 6 -> 19 of 45 units with a sustained pedestal 8 of 45 (18%) runs >=3 consecutive 29; 1-2 event runs (noise) 69 Also corrected: the affected channel is most often Vert, not Tran (v1 named whichever axis had the largest peak, so it was frequently wrong). BE10895 and BE18003 were invisible to v1. BE12599's fault began 2026-08-14, not 08-17. The decode itself was never in question and is confirmed against Blastware's own ASCII export: on K558LJN3.BK0W, BW shows Tran parked at +0.265..+0.375 in/s for the entire record while Vert and Long sit at ~0.005 — Instantel's "parallel lines above or below the zero line". Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HgTe8CamXAHcAmaQ6QNcog
This commit is contained in:
@@ -0,0 +1,108 @@
|
||||
#!/usr/bin/env python3
|
||||
"""Offset detector v2 — per-channel MEDIAN pedestal.
|
||||
|
||||
Supersedes the dominant-axis / mean detector in offset_scan.py, which had two
|
||||
flaws that manufactured false "recoveries":
|
||||
|
||||
1. It scored only the axis with the largest peak, so a real event on one axis
|
||||
hid a persistent pedestal on another. BE12599 2026-08-21 read "clean"
|
||||
because Long had a 1.065 in/s event, while Tran sat at +0.47 in/s.
|
||||
2. It used the MEAN, which a real transient perturbs. The median is the
|
||||
resting baseline: most samples sit at it, so a blast does not move it.
|
||||
Same event, Long: mean +0.0783 vs median -0.0050.
|
||||
|
||||
Flags a CHANNEL when |median| >= --floor in/s (default 0.025 = 5 A/D counts,
|
||||
Instantel's own criterion; 1 A/D count = 0.005 in/s).
|
||||
|
||||
Emits one row per (event, channel) so persistence can be tracked per channel.
|
||||
"""
|
||||
from __future__ import annotations
|
||||
import argparse, csv, re, statistics, 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
|
||||
|
||||
GEO = ("Tran", "Vert", "Long")
|
||||
K = 10.0 / 32000.0 # ADC counts -> in/s at the 10 in/s range
|
||||
_WAVE_RE = re.compile(r"\.[A-Za-z0-9]{2}0[Ww]$")
|
||||
_STEM_RE = re.compile(r"^([B-Z])(\d{3})")
|
||||
|
||||
|
||||
def serial_from_name(n):
|
||||
m = _STEM_RE.match(n)
|
||||
return f"BE{(ord(m.group(1))-ord('B'))*1000+int(m.group(2))}" if m else "?"
|
||||
|
||||
|
||||
def scan_one(ps):
|
||||
p = Path(ps)
|
||||
try:
|
||||
ev = read_blastware_file(p)
|
||||
s = ev.raw_samples or {}
|
||||
if not all(s.get(c) for c in GEO):
|
||||
return None
|
||||
ts = ev.timestamp
|
||||
stamp = (f"{ts.year:04d}-{ts.month:02d}-{ts.day:02d}T"
|
||||
f"{ts.hour:02d}:{ts.minute:02d}:{ts.second:02d}") if ts else ""
|
||||
out = []
|
||||
for ch in GEO:
|
||||
a = s[ch]
|
||||
out.append({
|
||||
"serial": serial_from_name(p.name), "timestamp": stamp,
|
||||
"filename": p.name, "channel": ch,
|
||||
"median_ips": round(statistics.median(a) * K, 4),
|
||||
"mean_ips": round(statistics.fmean(a) * K, 4),
|
||||
"peak_ips": round(max(abs(v) for v in a) * K, 4),
|
||||
})
|
||||
return out
|
||||
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("--floor", type=float, default=0.025)
|
||||
ap.add_argument("--out", required=True)
|
||||
a = ap.parse_args()
|
||||
|
||||
seen, files = set(), []
|
||||
for q in sorted(Path(a.dir).rglob("*")):
|
||||
if q.is_file() and _WAVE_RE.search(q.name) and q.name not in seen:
|
||||
seen.add(q.name); files.append(str(q))
|
||||
print(f"unique waveform binaries: {len(files)}", flush=True)
|
||||
|
||||
rows = []
|
||||
with ProcessPoolExecutor(max_workers=a.jobs) as ex:
|
||||
for n, f in enumerate(as_completed([ex.submit(scan_one, p) for p in files]), 1):
|
||||
r = f.result()
|
||||
if r: rows.extend(r)
|
||||
if n % 2000 == 0: print(f" {n}/{len(files)}", flush=True)
|
||||
|
||||
for r in rows:
|
||||
r["offset"] = int(abs(r["median_ips"]) >= a.floor)
|
||||
|
||||
cols = ["serial","timestamp","filename","channel","median_ips","mean_ips","peak_ips","offset"]
|
||||
with open(a.out, "w", newline="") as fh:
|
||||
w = csv.DictWriter(fh, fieldnames=cols); w.writeheader(); w.writerows(rows)
|
||||
|
||||
from collections import defaultdict
|
||||
ev_flagged = {(r["serial"], r["filename"]) for r in rows if r["offset"]}
|
||||
ev_all = {(r["serial"], r["filename"]) for r in rows}
|
||||
per = defaultdict(set)
|
||||
for r in rows:
|
||||
if r["offset"]: per[r["serial"]].add(r["filename"])
|
||||
tot = defaultdict(set)
|
||||
for r in rows: tot[r["serial"]].add(r["filename"])
|
||||
print(f"\nfloor = {a.floor} in/s ({a.floor/0.005:.0f} A/D counts)")
|
||||
print(f"events with >=1 offset channel: {len(ev_flagged)} of {len(ev_all)}")
|
||||
print(f"units affected: {len(per)} of {len(tot)}")
|
||||
for s in sorted(per, key=lambda s: -len(per[s])):
|
||||
print(f" {s:9} {len(per[s]):4} / {len(tot[s]):4} events")
|
||||
print(f"\nwrote {a.out}")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
Reference in New Issue
Block a user