7 Commits
26 changed files with 134 additions and 2888 deletions
+4 -229
View File
@@ -4,224 +4,7 @@ All notable changes to seismo-relay are documented here.
---
## v0.30.0 — 2026-09-12
**The series-4 correctness release** — the Thor / Micromate counterpart to
v0.26.0's series-3 work. The decoder is now verified per-sample against
Thor's own CSV exports: **459 waveform files, 3,807,158 / 3,807,165 samples
exact** across three independent ground-truth corpora, and production IDFW is
**575/575** with zero truncations and zero decode failures. Series-3
re-verified **unchanged at 14,338/14,338** after every shared-codec change.
⚠ **This release owes the prod store a Thor backfill.** Every stored
series-4 geophone value is **3.3% low**, and histogram peaks from monitoring
runs longer than ~4 hours can be far worse (the interval cap discarded the
tail, frequently the part holding the peak). Run
`scripts/backfill_thor_events.py` — `TOOL_VERSION` is bumped to `0.30.0`, so
regeneration is gated correctly and **no `--force` is needed**. DB backup
first. Series-3 events are untouched by this release and do not need
re-running.
⚠ **Terra-View displays these values.** Series-4 geophone readings will rise
~3.3% after the backfill, and some histogram PPVs will rise a great deal more.
That is a correction, not a regression.
### Fixed — event-report PDF used a per-trace geo Y scale
The waveform plot scaled each geo lane to its own peak, so a small channel
filled its lane and looked as large as a big one, and the `Geo: X in/s/div`
footer reflected only whichever channel was measured first — wrong for the
other two. All three geo lanes now share one symmetric scale (max |sample|
across them, padded, 0.05 in/s floor), matching the event modal and BW's
single amp/div; the footer reflects that shared scale. Mic keeps its own psi
scale. Large events are unchanged.
### Fixed — series-4 (Thor / Micromate) decoder is now per-sample exact
Verified against **Thor's own CSV exports**, which carry a per-sample
four-column block beside every binary (`CSV/<name>.IDFW.csv`) — 1,012 paired
files that had been sitting in the corpus unused. Previous notes asserted
"Thor has no ASCII ground truth", which is why the decoder stayed pinned to a
superseded walker with an unverifiable scale factor.
| metric | before | after |
|---|---|---|
| IDFW per-sample exact | 39.1% | **100.000%** (1,057,536/1,057,536) |
| IDFW files fully exact | 0/153 | **153/153** |
| IDFW PPV median error | −3.32% | **−0.002%** |
| IDFH within 2% of Thor PPV | 51.1% | **100.0%** (858/858) |
| prod IDFW PPV median error (8 units) | −3.3% | **−0.001%** |
| decode cost | — | 6 ms/file |
Four independent root causes:
- **Geo LSB was `0.0003`, should be `0.000310308`** — the old value was Thor's
4-decimal *display rounding* of the LSB mistaken for the LSB, so every
series-4 geophone sample read **3.3% low**. Pinned to ±6e-11 by
intersecting 991,415 rounding constraints; corroborated by the ±full-scale
seed (`±32226`) in unwritten IDFH slots. Applies to IDFH too, which had a
separate (also wrong) `10.0/32768`.
- **IDFH histograms were capped at 250 intervals** — the segment validator
required the interval counter's high byte to be zero, but the counter is a
uint16 cumulative index, so every segment past interval 255 was rejected.
Any run over ~4 hours lost its tail, often the part holding the peak.
540/858 corpus files affected.
- **Record mode `00 00` (raw int16, 10-byte header) was unhandled** — the
record fell through the dispatch, silently dropping each channel's first
512 samples. This produced the long-standing "loud events truncate"
symptom. `MODE_ABSOLUTE` is now also accepted as a segment-0 preamble.
- **Body-offset search matched `00 02 00` inside record headers** — picking a
candidate part-way down the chain, which decodes a rotation-shifted body
that drops each channel's segment 0. The search now anchors on record
headers and takes the chain head.
Also fixes the separately-tracked "UM-series decodes ~1000× low" bug
(`UM11402_20260406130113.IDFW` now matches its device report exactly).
Series-3 re-verified **unchanged at 14,338/14,338 exact** after the shared
`waveform_codec` change.
⚠ **This is a codec change: the Thor store owes a regeneration.** Run
`scripts/backfill_thor_events.py` (bump `TOOL_VERSION` first, or pass
`--force`), DB backup first. All stored series-4 `.h5`/sidecar peaks are
currently ~3.3% low, and histogram peaks for runs over ~4 hours may be
badly low.
⚠ **Thor's histogram PPV has a 0.0050 in/s display floor** — 41.4% of prod
IDFH sidecars report a component PPV larger than their own vector sum. On
quiet files the decoder is now *more* accurate than that reference.
New: `scratch/verify_thor_against_csv.py`, `tests/test_idf_binary_codec.py`
(10 tests, fixtures under `tests/fixtures/thor-idf/`).
### Fixed — mic-disabled (3-channel) units
Verified on a second corpus (`9-10-26-csv-req`: UM11402, UM12947, UM20147) —
**139/139 waveforms per-sample exact (1,273,380 samples), 877/877 histograms
within 2%** (was 66.9% and 56.6%).
- **Waveform body head sat below the scan floor.** A 3-channel unit's shorter
header puts the record chain head at `0x0dba`, under the old
`_BODY_SCAN_FLOOR` of `0x0E00`. The scan couldn't see it and fell through
to the Vert segment-0 record, decoding a body shifted one position around
the channel rotation — Vert came up exactly 512 samples short. Floor
lowered to `0x0C00`; body-offset scoring now accepts 3 channels as "equal"
instead of demanding 4.
- **Histogram interval record is 56 bytes, not 72.** It is
`16 × n_channels + 8`, so mic-disabled units pack 56. Assuming 72 read 7
intervals out of every 10-interval segment then walked off alignment into
garbage decoding as ~10 in/s peaks (errors up to +191,000%). The interval
count now comes from the segment's cumulative counter and the stride is
derived from it; also recovers 4 files that decoded no intervals at all.
Combined across both corpora: **292/292 waveform files, 2,330,916/2,330,916
samples exact.** Production IDFW truncations 41 → 22.
### Fixed — `40 NN` int16 blocks with NN > 8
`data_block_len()` rejected any `40 NN` block with `NN > 0x08`. The cap had
no evidence behind it: every corpus available when it was written used only
NN ∈ {1,2,3,4,8}, so it was never exercised. Loud UM12947 events use NN of
12, 16, 20 … up to 196, and because the block walker stops at the first
unrecognised tag rather than raising, rejecting them surfaced as **silently
short channels** (e.g. Tran 1812 / Vert 2132 / Long 2324 on a file whose
export has 2324 for all three). The bound is the buffer, not a constant.
Verified against Thor exports for UM12947 (2025-07-14 … 09-25, 167
waveforms): length mismatches **22 → 0**, **1,476,242/1,476,249** samples
exact. These are not truncated recordings — the exports carry full sample
counts.
`tests/test_waveform_codec.py` asserted the cap as intended behaviour; that
assertion was wrong and has been replaced with one pinning the opposite,
carrying the evidence.
### Result across all three ground-truth corpora
**459 waveform files, 3,807,158 / 3,807,165 samples exact.** Production
IDFW: **575/575**, zero truncations, zero decode failures, median PPV error
−0.0007% across 8 units. Series-3 re-verified **unchanged at 14,338/14,338**
after every shared-codec change.
The 7 residual samples each differ by one 4th-decimal tick and are **Thor's
own rounding**: intersecting the per-sample rounding constraints over that
corpus is infeasible (the binding pair contradict by 2.3e-11, 7e-5 relative),
so no single linear LSB reproduces every printed value. `_GEO_LSB_IPS` is
already pinned to ~1e-11 — do not retune it to chase these.
---
## v0.29.0 — 2026-09-04
First release to reach prod since **v0.27.0**, so it ships **both** the
`false_trigger_reason` column below *and* the v0.28.0 offset (DC-baseline)
detector: v0.28.0 was version-bumped in-tree (`TOOL_VERSION`, CHANGELOG) but
never tagged or deployed, so 0.29.0 is the first build to carry either to prod.
Pairs with Terra-View ≥ 0.24.0. The `false_trigger_reason` column auto-migrates
on startup; the offset detector still needs the shape backfill on the prod store
(`scripts/backfill_event_shape.py`) to populate `shape_offset*` on existing rows.
### Added
- **`events.false_trigger_reason` — optional FT cause.** A nullable `TEXT`
column recording *why* an event is a false trigger (e.g. `"offset"`), as a
subtype of the FT flag: setting a reason via the sidecar review PATCH implies
`false_trigger=1`, and the reason is cleared whenever FT ends up 0
(confirm-real, clear-FT, `set_false_trigger(false)`). `propagate_review_to_twins`
carries the reason to the histogram/waveform twin alongside the flag.
Auto-migrated (`_SCHEMA` + `_migrate` ADD COLUMN — not the Migration-1
rebuild); exposed via `/db/events`. Terra-View surfaces it as a manual
"Flag as offset" action + an `FT · offset` badge.
### Fixed
- **BlastMate serials — the family prefix is read from the file, not guessed.**
The Blastware filename encodes only the serial *number* (`L895…` → 10895);
the two-letter prefix is not in it. `waveform_store` synthesised `"BE"`, so
an imported **BlastMate** (serials `BA…`) was filed under a MiniMate Plus
serial that does not exist — silently, and Terra-View read it straight
through. `save_imported_bw` now resolves serial as hint → file body →
filename guess, via a new `_serial_from_bw_bytes` that accepts a candidate
only when its numeric part matches the filename. `client._decode_0a_partial_header`
likewise matched a literal `b"BE"` in monitor-log partial records; on a
BlastMate that returned −1 and skipped the whole block, losing the **geo
threshold** along with the serial. It now matches any two-letter prefix and
requires the NUL terminator — stricter than the search it replaces.
BlastMate is the MiniMate Plus's larger Series III sibling and its files are
byte-compatible: all 1,493 in the DL2 archive decode through the existing
codec at 100%, same four channels. **The serial string was the only thing
blocking BlastMate support in SFM.** Four archive units were affected —
BA9229, BA10060, BA10895, BA15957.
**No backfill and no `TOOL_VERSION` bump**: this changes which serial an
*import* is filed under, not any decoded value, so existing sidecars and
`.h5` files are untouched. **No migration either** — prod holds no BlastMate
events (the archive's BA units last recorded 2018-10 through 2023-11; the
prod backfill reaches back only to ~May 2025).
---
## v0.28.0 — 2026-09-02
**Offset (DC-baseline) false-trigger detector.** Productionizes the validated
pre-trigger detector: a geophone event whose baseline sits off zero and stays
flat across the record (sensor bumped / settled / drifted) is now flagged and
surfaced in Terra-View as an `offset` false-trigger reason — catching offsets the
crest/near-peak spike rule misses (an offset is low-crest and flat).
### Added
- `shape_metrics.offset_from_samples` / `offset_from_h5`: per geophone channel,
`|median(pre-trigger)| ≥ 0.025 in/s` AND `pre/mid/end spread ≤ 0.02` → offset;
the consistency test rejects transients (a real event moves one third). Reads
the `.h5` samples + the `pretrig_samples` attr, range-aware via the in/s float
samples. Constants `OFFSET_FLOOR` / `OFFSET_MAX_SPREAD` are tunable.
- `events.shape_offset` / `shape_offset_axis` / `shape_offset_pre` /
`shape_offset_spread` columns (auto-migrated: `_SCHEMA` + the `_migrate`
ADD COLUMN loop), computed at all three ingest paths and by
`backfill_event_shape.py`, exposed via `/db/events`.
Requires the shape/offset backfill on the prod store to populate existing events:
`python scripts/backfill_event_shape.py --db-path … --store-root …`.
## [Unreleased]
---
@@ -256,17 +39,9 @@ carried, and it found one real codec bug (below).
walk double-counts every binary — 127,035 paths are 63,535 distinct files. The
ASCII exports are *not* mirrored, so the 14,340 pair count is already distinct.)
**No prod backfill is required for this.** Verified after the fact: all four
recovered files are archive-only — none exists in the production store or the
events DB — and re-running stride detection over the production store's
**10,215** histogram binaries shows **0 files whose decode changes**. The fix
matters for future ingests of sub-minute histograms with a partial final block,
not for anything already stored.
(`TOOL_VERSION` moves with the release, so whenever a backfill *is* next run for
some other reason it will regenerate the whole store rather than skipping. That
is harmless — the output is byte-identical for every currently-stored file — but
it means the run takes its full ~2 hours on the NAS.)
⚠ Prod stores hold `.h5` files generated before this fix. Those 4 events stay
empty until `backfill_sidecars.py` is re-run — not worth a two-hour prod backfill
on its own; fold it into the next one.
- **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)
+26 -83
View File
@@ -2,11 +2,7 @@
Ground-up Python replacement for **Blastware**, Instantel's Windows-only software for
managing MiniMate Plus seismographs. Connects over direct RS-232 or cellular modem
(Sierra Wireless RV50 / RV55). Current version: **v0.30.0**.
Stack-level context — which repo owns what, and how the three project versions
pair — lives in `../terra-view/docs/tmi-stack.md`, which is also loaded as
`~/CLAUDE.md`.
(Sierra Wireless RV50 / RV55). Current version: **v0.27.0**.
---
@@ -24,43 +20,9 @@ Read this first when picking the project back up.
Independent corroboration of the 32000-count scale: 19,244 healthy
channel-events sit at a pre-trigger floor of exactly 0.000 (62.7%), 94.5%
within ±1 quantisation unit, median +0.0000 — no zero-point bias.
- **Series-4 (Thor / Micromate) is now verified per-sample (2026-09-10).**
**1,057,536 / 1,057,536** geo samples across all 153 genuine Thor waveform
files reproduce Thor's own CSV export exactly; IDFH peaks are within 2% on
858/858 (median -0.004%). The ground truth was in the corpus all along —
Thor writes `CSV/<name>.IDFW.csv` beside each binary with a **per-sample**
four-column block. Harness: `scratch/verify_thor_against_csv.py`.
Four bugs, all fixed: geo LSB was `0.0003` (display rounding of the real
`0.000310308`, so every sample read **3.3% low**); the IDFH segment
validator required a zero counter high byte, **capping every histogram at
250 intervals**; record mode `00 00` (raw int16) was unhandled, silently
dropping each channel's first 512 samples; and the body-offset search
matched `00 02 00` *inside* record headers, decoding a rotation-shifted
body. IDFW is no longer pinned to `decode_waveform_legacy`.
Series-3 re-verified unchanged at 14,338/14,338 after the shared-codec
change.
- **Mic-disabled (3-channel) units are a distinct shape (2026-09-10).**
Verified on a second corpus (`~/thor-csv-req`, UM11402/UM12947/UM20147):
**139/139** waveforms per-sample exact, **877/877** histograms within 2%.
Two structural differences: the shorter header puts the waveform record
chain head at `0x0dba` (below the old `_BODY_SCAN_FLOOR` of `0x0E00`, so it
was invisible and Vert came up exactly 512 short), and the histogram
interval record is **56 bytes, not 72** — `16 × n_channels + 8`, derived per
segment from the cumulative interval counter, never assumed.
- **`40 NN` blocks are not capped at NN=8 (2026-09-11).** `data_block_len()`
rejected `NN > 0x08`, a guard with no evidence behind it — the corpora
available when it was written only used NN ∈ {1,2,3,4,8}. Loud UM12947
events use NN up to 196, and since the walker stops at the first
unrecognised tag rather than raising, this surfaced as silently short
channels. Verified on 167 UM12947 waveforms: length mismatches 22 → 0,
1,476,242/1,476,249 samples exact.
- **Production IDFW is now 575/575** — zero truncations, zero decode
failures, median PPV error −0.0007% across 8 units (was 41 truncated + 1
failing, −3.3%). Across all three ground-truth corpora: **459 files,
3,807,158/3,807,165 samples exact**; the 7 stragglers differ by one
4th-decimal tick and are Thor's own rounding — no single linear LSB can
reproduce every printed value (the constraints are infeasible by 7e-5
relative), so do NOT retune `_GEO_LSB_IPS`.
- **Series-4 (Thor / Micromate) is NOT verified.** UM-series sits at ~48%
against device peaks with a ~1.7% systematic bias and a near-zero tail.
Thor IDFW is pinned to `decode_waveform_legacy` deliberately.
- **Open, not blocking:** 14 sensitive-range files show an exact 8x
(= 10.0/1.25) units discrepancy; `scripts/backfill_sidecars.py --force` also
inserts DB rows for store files that have none (one-time per store) and the
@@ -71,9 +33,8 @@ Read this first when picking the project back up.
(it gates regeneration). ⚠ On the office NAS this takes **~2 hours**
(~1.5 files/sec vs 85/sec on the dev box — gzip-4 in `sfm/event_hdf5.py`
against a Synology CPU). Budget it up front.
**v0.27.0 does NOT owe prod a backfill** — verified: the partial-final-block
fix changes 0 of the 10,215 histograms in the prod store (the 4 recovered
files are archive-only and were never ingested).
**v0.27.0 owes prod a backfill:** the partial-final-block fix recovers 4
histograms that are still empty in the store.
- **The "offset" hardware fault has its own journal** --
`docs/offset_investigation.md`. **5 of 45 units (11%)**, and the fault is
**persistent** — it stays until the geophone is serviced. Detect it with
@@ -154,34 +115,20 @@ should not import from `sfm/`, must not touch a DB, and have no I/O
beyond reading files passed as arguments. Keep them pure — both
tiers can then depend on them without circularity.
#### Thor IDF binary codec (updated 2026-09-10)
#### Thor IDF binary codec (2026-05-28)
`micromate/idf_file.read_idf_file()` decodes both Thor IDFW
(waveform) and IDFH (histogram) binaries. **Verified per-sample
against Thor's own CSV exports** — see
`scratch/verify_thor_against_csv.py`.
(waveform) and IDFH (histogram) binaries.
- **IDFW** uses the series-3 record-chain `decode_waveform_v2()`. The
body offset is **not** fixed: it is `<chain-head record> + 7`, found
by `_find_waveform_body_offset()` anchoring on record headers. All
**153/153** genuine Thor waveform files decode per-sample exact
(1,057,536/1,057,536 samples).
- **IDFH** segment header is `[len_be][0a 00 00 00][counter_be][05 3f]`,
where `counter` is a **uint16 cumulative interval index** — it must
not be constrained to a zero high byte (that capped histograms at 250
intervals). Intervals whose `min > max` on all channels are unwritten
slots carrying a ±full-scale seed and are skipped. 858/858 files land
within 2% of Thor's PPV (median -0.004%).
- **Geo LSB is `0.000310308` in/s per count** (full scale 10.0 in/s =
32226.05 counts). Series-3's 32000-count scale does NOT apply.
- **Record modes** are `02 00` deltas (14 B header), `01 00` absolute,
`00 03` raw 12-bit, and `00 00` **raw int16** (all 10 B headers).
`01 00` and `00 00` are also valid as the implicit segment-0 preamble.
⚠ **Thor's histogram PPV has a 0.0050 in/s display floor.** 41.4% of
prod IDFH sidecars report a component PPV exceeding their own vector
sum — impossible. On quiet files our decode is *more* accurate than
the reference; do not "fix" the decoder to match it.
- **IDFW** reuses `decode_waveform_v2()` on the body at fixed file
offset `0x0f1f`. Sample fidelity is 87–99% byte-exact on quiet
events; loud events hit the BW codec's known walker-stops-early
limitation.
- **IDFH** has its own segment-based decoder: `[len_be][0a 00 00 00]
[00 NN][05 3f]` + N × 72-byte interval records (4 × 16-byte
per-channel min/max/halfp). All 859 Thor IDFH corpus files
decode (181,071 intervals); peak matches sidecar within ~1.8%
(ADC quantization).
The two outlier `BE9439_*` files in the Thor example corpus are
actually Series III Blastware binaries that share the `.IDFW`/`.IDFH`
@@ -447,19 +394,15 @@ with zero mismatches. Before: 1 of 1196.
`BE12599/N599LPWJ.980W` @849, `BE9558/K558LOF2.820W` @1485.
(The series-3 histogram codec was fixed 2026-08-25 — see below.)
- ~~**Micromate (UM-series) IDF decode is ~1000× low**~~ — FIXED 2026-09-10.
`UM11402_20260406130113.IDFW` now decodes Tran 1.1168 / Vert 4.3220 /
Long 0.9135, matching the device report exactly. Root cause was the
body-offset search landing inside a record header plus the unhandled
`00 00` record mode, not anything UM-specific.
- ~~**Thor IDF per-count LSB**~~ — RESOLVED 2026-09-10. The 0.983 ratio was
exactly `0.0003 / 0.000310308`. Thor's geo LSB is **0.000310308 in/s per
count** (full scale 10.0 in/s = 32226.05 counts), pinned to ±6e-11 by
intersecting 991,415 rounding constraints from Thor's own exports and
corroborated by the ±full-scale seed (`±32226`) left in unwritten IDFH
interval slots. Series-3's 32000-count scale does **not** carry over.
Note `10.0/32226` is very slightly wrong — see
`docs/idf_protocol_reference.md`.
- **Micromate (UM-series) IDF decode is ~1000× low** — e.g.
`UM11402_20260406130113.IDFW` gives a Tran peak of 0.0009 in/s against
a device-reported 1.1168. The Thor IDF path decodes sanely, so this
is UM-specific.
- **Thor IDF per-count LSB** — after the 32000 geo full-scale
correction, series-4 Thor peaks sit at a median 0.983 of the
device-reported peak (was 0.960 under 32768). Closer but not exact;
Thor likely uses its own per-count LSB rather than the BW
16-count/0.005 in/s convention.
### Decoded sample counts (across the fixture bundle)
+1 -1
View File
@@ -1,4 +1,4 @@
# seismo-relay `v0.30.0`
# seismo-relay `v0.27.0`
A ground-up replacement for **Blastware** — Instantel's aging Windows-only
software for managing seismographs. Supports both the **MiniMate Plus
+1 -223
View File
@@ -6,15 +6,7 @@ Series IV event-file format. Sibling to
Series III "Rosetta Stone") — this doc holds what we know so far and
the open questions still to crack.
> ⚠ **The "Status (2026-05-28)" block below is SUPERSEDED.** Its geo LSB
> (0.0003), its IDFH scale (`/32768 × 10`), its fixed body offset (`0x0f1f`)
> and its "87–99% byte-exact / loud events truncate" caveat were all wrong or
> incomplete. See **[Verified against Thor's own exports
> (2026-09-10)](#verified-against-thors-own-exports-2026-09-10)** — the
> decoder is now per-sample exact on 1,057,536/1,057,536 samples. The block
> is kept only for the reverse-engineering trail.
**Status (2026-05-28, SUPERSEDED):** ASCII text sidecar fully decoded (1,014
**Status (2026-05-28):** ASCII text sidecar fully decoded (1,014
sample files round-trip). **Thor IDFW** binary now decodes via
`micromate.idf_file.read_idf_file()` — reuses the BW segment-rotated
block codec verbatim at fixed body offset `0x0f1f`; metadata (serial,
@@ -52,220 +44,6 @@ signature and raises `NotImplementedError` pointing callers at
time-of-peak); the two uint16 fields (probably PVS contributions);
8-byte interval tail (PVS data); mic dB(L) exact conversion constant.
## Verified against Thor's own exports (2026-09-10)
**The series-4 decoder is now per-sample exact.** 1,057,536 / 1,057,536
geophone samples across all 153 genuine Thor waveform files reproduce Thor's
own CSV export exactly; histogram peaks land within 2% on 858/858 files
(median error −0.004%).
### Ground truth — it was there all along
Thor writes `TXT/`, `CSV/`, `XML/` and `PDF/` exports beside every binary:
```
<serial dir>/UM13981_20220207084555.IDFW
<serial dir>/CSV/UM13981_20220207084555.IDFW.csv
```
The **CSV carries a per-sample block** — four columns (Tran, Vert, Long, Mic)
in in/s and psi, after the 2-column report header. That is the series-4
equivalent of Blastware's `_ASCII.TXT` exports, and it gives 1,012 paired
files (152 IDFW + 860 IDFH). Earlier notes in this file and in
`micromate/idf_file.py` asserted "Thor has no ASCII ground truth in the
corpus"; that was wrong, and it is why the decoder sat pinned to a
superseded walker with a scaling constant nobody could check.
Harness: `scratch/verify_thor_against_csv.py`.
### Geo LSB = 0.000310308 in/s per count (NOT 0.0003)
The old 0.0003 was read off the smallest non-zero sample in the exports —
but that is Thor's **4-decimal display rounding of the LSB, not the LSB**.
It read every series-4 geophone sample **3.3% low**. The quantisation
ladder gives it away: counts 1..6 export as 0.0003, 0.0006, 0.0009, 0.0012,
0.0016, 0.0019 — an LSB of exactly 0.0003 would end 0.0015, 0.0018.
Each exported sample constrains the LSB to the window that rounds to its
printed value. Intersecting 991,415 such constraints gives
```
LSB ∈ [0.000310307933, 0.000310308057] width 1.2e-10
```
so `_GEO_LSB_IPS = 0.000310308`, i.e. full scale 10.0 in/s = **32226.05
counts**. Corroboration: an IDFH interval that never recorded keeps its
min/max accumulator at its ±full-scale seed, and that seed is
`(min=+32226, max=-32226)`. ⚠ The tempting closed form `10.0/32226` is
very slightly wrong — it lands 4.5e-10 above the feasible window and loses
78 boundary samples while never winning one. **Series III uses 32000 counts
for the same 10.0 in/s, so the two generations do not share a scale.**
Independently confirmed on 8 production units (UM6047, UM11402, UM11719,
UM12947, UM13981, UM14133, UM20146, UM20147): every unit's median PPV error
against its device-reported peak moved from −3.3% to within ±0.03%. It is a
global constant, not a per-unit calibration.
### IDFH segment header: the counter is a uint16, and it is cumulative
```
[length_be 2B][0a 00 00 00][counter_be 2B][05 3f]
```
`counter` is the **0-based cumulative index of the last interval in the
segment** — 9, 19, 29, ... for the usual 10-intervals-per-segment layout
(`length` = 730).
The validator used to require `counter`'s high byte to be `0x00`. That
silently **capped every histogram at 250 intervals**: once the cumulative
counter passed 255 the high byte went non-zero and every later segment was
rejected. Any run longer than ~4 hours lost its tail — frequently the part
holding the event peak, so the file's PPV read low. **540 of 858 corpus
files were affected**; fixing it moved histogram peaks from 48.3% to 93.8%
within 0.5% of Thor's reported PPV.
### Unwritten interval slots carry a ±full-scale seed
An interval the device reserved but never wrote keeps `min = +32226`,
`max = -32226` on all four channels — `min > max`, impossible for real data.
Decoded naively it yields a 10.0 in/s peak on every channel and, being a
max-over-intervals, poisons the whole file's PPV. Rare but real: exactly 1
of 497,611 corpus intervals, and it inflated that file's Long PPV from
0.0081 to 10.0 in/s. The inversion is all-or-nothing across channels (0
partial cases), so requiring every channel to be inverted is a safe test.
### Record mode `00 00` — raw int16 absolute (MODE_RAW16)
The record chain's mode field at `off+8` takes a fourth value:
| mode | meaning | header |
|---|---|---|
| `02 00` | deltas + two int16 anchors | 14 B |
| `01 00` | absolute, tagged blocks | 10 B |
| `00 03` | raw 12-bit absolute, untagged | 10 B |
| **`00 00`** | **raw int16 BE absolute, untagged** | **10 B** |
A `MODE_RAW16` record with `length = 1032` carries exactly
`(1032 - 8) / 2 = 512` samples and reproduced Thor's export **512/512
exactly** on first test. Thor uses it for segment 0 (the pre-trigger
window) on some events. Before this mode existed the record fell through
the dispatch unhandled, so the channel silently lost its first 512 samples —
which is what produced the "loud events truncate" symptom.
`MODE_ABSOLUTE` is also valid as a **preamble** (the implicit segment-0 Tran
record); its tagged blocks start at `body[3]`, not `body[7]`, because its
header is 10 bytes rather than 14.
### Body offset is not fixed at 0x0f1f — and 0x0f1f is really a record + 7
A "body offset" is `<record start> + 7`, so that `body[0]` is the segment
index and `body[1:3]` is the mode. The canonical `0x0f1f` is simply the
record at `0x0f18`.
Searching for the literal preamble `00 02 00` finds only MODE_DELTA bodies,
and worse, it **matches the `[seg][mode]` bytes inside any record header**,
so the scan could pick a candidate part-way down the chain. That decodes a
plausible-looking but rotation-shifted body which drops each channel's
segment 0 — the real cause of the remaining truncations.
`_find_waveform_body_offset()` now anchors on record headers (the
`<channel_id> 00 00` signature at `+4`, validated with `is_record()`),
takes the **chain head** — a record no other record's length field points at
— and trial-decodes `head + 7`, preferring the candidate where all four
channels come out the same length.
⚠ Do **not** scan for candidate preambles instead: `MODE_RAW16` is
`00 00`, so every run of three zero bytes looks like a body start and each
costs a full trial decode (~0.5 s/file measured, vs 6 ms/file now).
### `40 NN` is not capped at NN=8 (2026-09-11)
`data_block_len()` rejected any `40 NN` int16 block with `NN > 0x08`. The cap
had no evidence behind it — every corpus available when it was written used
only NN ∈ {1, 2, 3, 4, 8}, so it was never exercised. Loud events use much
wider blocks:
| corpus | `40 NN` values | walker stops |
|---|---|---|
| first + 3-channel corpora | 1, 2, 3, 4, 8 | none |
| UM12947 2025-07..09 | 2, 4, 8, **12, 16, 20 … 196** | every value > 8 |
Because `walk_body`/`run` stop at the first unrecognised tag rather than
raising, this surfaced as **silently short channels** — e.g. Tran 1812 /
Vert 2132 / Long 2324 on a file whose export has 2324 for all three. The
real bound is the buffer (and the caller's record end), not a magic constant.
Verified against Thor's exports for UM12947 (2025-07-14 … 2025-09-25, 167
waveforms): length mismatches **22 → 0**, and **1,476,242 / 1,476,249**
samples exact.
⚠ These events are **not** truncated recordings, which was the competing
hypothesis — the exports carry the full sample count.
**The 7 residual samples are Thor's rounding, not ours.** Each differs by
exactly one 4th-decimal tick (e.g. decoded 3.3551 vs export 3.3550).
Intersecting the per-sample rounding constraints over this corpus is
**infeasible** — the binding pair (count 2013 → 0.6247, count 4351 → 1.3501)
contradict by 2.3e-11, i.e. 7e-5 relative. No single linear LSB can
reproduce every printed value, so Thor is not doing plain round-half-up on
`count × LSB`. Do not retune `_GEO_LSB_IPS` to chase these; it is already
pinned to ~1e-11.
### Mic-disabled units are a distinct shape (2026-09-10, second corpus)
Some units run with the microphone disabled — **3 channels, not 4** — and that
changes two structural things. Confirmed on the `9-10-26-csv-req` corpus
(UM11402, UM12947, UM20147): 139/139 waveforms and 877/877 histograms.
**Waveform: the body starts earlier.** A 3-channel unit has a shorter fixed
header and puts its record chain head at **`0x0dba`**, below the old
`_BODY_SCAN_FLOOR` of `0x0E00`. The head was therefore invisible to the scan,
which fell through to the *Vert* segment-0 record and decoded a body shifted
one position around the channel rotation. The signature is unmistakable:
```
Tran 3072 / Vert 2560 / Long 3072 / MicL 0 <- Vert exactly 512 short
```
46 of 139 files in that corpus were affected; all 46 became per-sample exact
once the floor dropped to `0x0C00`. Note the body-offset scoring also had to
stop requiring four channels — `len(lengths) >= 3`, not `== 4`, or `equal` is
permanently False for these events and the pick falls back to raw sample count.
**Histogram: the interval record is 56 bytes, not 72.**
```
interval_size = 16 × n_channels + 8 (72 for 4 channels, 56 for 3)
```
It is **not a constant**, and it cannot be inferred from `length` alone.
Derive the interval count from the segment counter — it is cumulative, so
`n = counter - previous_counter` — and then `stride = (length - 10) / n`.
`n_channels` follows from `(stride - 8) / 16`.
Assuming 72 read 7 intervals out of each 10-interval segment and then walked
off alignment into garbage that decoded as ~10 in/s peaks — inflating those
files' PPV by up to 191,000%. Fixing it moved the second corpus from 56.6% to
**100.0%** of histograms within 2% of Thor's reported PPV, and recovered 4
files that previously decoded no intervals at all.
### What is still open
- ~~23 of 575 production IDFW files~~ — **RESOLVED 2026-09-11.** Production
IDFW is now **575/575** with zero truncations and zero decode failures
(median PPV error −0.0007%). See "`40 NN` is not capped at NN=8" above.
- Mic → psi scale is still the rough `2.14e-6` regression, not derived.
- Per-channel `int16 field4` in the IDFH interval record (possibly
time-of-peak) and the 8-byte tail (PVS data) remain undecoded.
⚠ **Thor's histogram PPV has a display floor of 0.0050 in/s.** In the
production store 6,080 sidecar PPV values are exactly 0.0050 (next most
common value: 275 occurrences), and **41.4% of IDFH sidecars report a
component PPV larger than their own vector sum** — geometrically impossible.
On those quiet files the decoder's ~0.0025 in/s is *more* accurate than the
reference; do not "fix" the decoder to match it.
### Codec breakthroughs (2026-05-28)
- **Body offset is a fixed `0x0f1f`** across 151/154 corpus IDFW
+3 -476
View File
@@ -58,14 +58,6 @@ 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.
@@ -154,8 +146,8 @@ is the discriminator — and it is what the field experience predicts.
| runs of >=3 consecutive | — | 29 |
| runs of 1-2 events (noise) | — | 69 |
Units with a sustained pedestal: **BE9558, BA10895, BE11007, BE11529, BE12599,
BE13117, BE18003, BE18438**. BA10895 and BE18003 were invisible to v1.
Units with a sustained pedestal: **BE9558, BE10895, BE11007, BE11529, BE12599,
BE13117, BE18003, BE18438**. BE10895 and BE18003 were invisible to v1.
**The affected channel is most often Vert**, which v1 got wrong — it named
whichever axis had the largest peak. BE13117 and BE18438 are both Vert faults.
@@ -226,7 +218,7 @@ signal from a tuned one:
**BE9558, BE11529, BE12599, BE13117, BE18438.**
Unchanged across a 2x threshold range. BE11007 and BA10895 drop out — the
Unchanged across a 2x threshold range. BE11007 and BE10895 drop out — the
spread test identifies them as transients, not pedestals.
The 11% headline happens to match v1's, but the reasoning and the unit list
@@ -503,459 +495,6 @@ 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), BA10895|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.
### BA10895 — 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, BA15957, BE16754, BE16758, BE8081, BE8626, BA9229, 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.
---
## 8c. Mechanism — five hypotheses tested, all dead (2026-09-06)
**The mechanism is still unknown.** Five campaigns, ~105 effectively independent
tests, seven nominally significant results against **5.2 expected by chance**
under a global null. Every one died to its own confound analysis. What the
campaign bought is a set of *shape constraints* and a long list of dead ends.
### ⚠ Two things retracted from this journal
**1. "Polarity is perfectly consistent — 11 of 11, zero mixed cases."** That is
a **tautology of the spread gate**, not a property of the fault. `spread <= 0.02`
requires pre/mid/end to agree, which forces one sign. Amplitude-only at the same
0.025 threshold: **12 of 53 unit-channels are mixed**, including BE18438|Vert
(88+/1−) and BE9558|Vert (1+/35−). Withdrawn.
**2. "5 of 45 units, unchanged across a 2x threshold range."** The
threshold-insensitivity is also a property of the gate. Amplitude-only gives
**9 units at 0.020, 8 at 0.025** (adding BA10895, BE12844, BE18003), 5 at 0.040.
The fleet is **8–9 units, not 5**.
**3. "Persistent — it stays until the geophone is serviced."** Weakened, not
withdrawn. There are **23 recoveries after runs of >=3 flagged events, median
gap 6.03 days**, three inside ten minutes. BE18438|Vert reads `pre=mid=end=
+0.0000` on 2026-02-10, +0.185→+0.370 across 02-25/26, and `+0.0000` again on
2026-03-22 — identical Project, Seis Loc, calibration date, geo range and
trigger throughout. The one thing that cannot be excluded is a **field
autozero**: it is a button sequence at the unit and writes nothing into the
event header. So "persistent" may be "persistent unless somebody pressed the
buttons," and the archive cannot tell those apart.
### The one positive finding: onset is a RAMP, minutes to hours
Both onsets resolvable at minute cadence are ramps. **BE18438|Vert,
2026-02-20** — the histogram corpus collapses a 14 d 21 h waveform bracket to
**one minute**:
```
~14,200 consecutive quiet minutes at 0.000–0.005 (ten full daily files)
09:32 +0.005 09:39 +0.045 10:20 +0.125 16:00 +0.165
09:33 +0.010 09:42 +0.070 13:13 +0.150 20:17 +0.185 plateau
```
**50% of the excursion in 7 minutes**, the rest asymptotic over ~10 h, **>=25
distinct one-minute intermediates**. Validated **75/75** against Blastware's own
ASCII export. Its 2025-11-15 onset is the same shape over 2.7 h. BE13117 stage B
is a 91-minute monotone rise, +0.035 → +1.745 in/s over ~40 samples.
**This kills both poles of the original dichotomy** (journal Q1): not an
instantaneous latched step (a bad autozero, a stuck trim-DAC), and not slow
component degradation over days or weeks.
⚠ It rests on **2 of 45 instruments**. Clopper-Pearson on 4/4 resolved onsets
gives 95% CI [0.40, 1.00] — a mixed population with up to 60% true steps is not
excluded. BE13117 has zero paired ASCII, so its ramp rests on our decoder alone.
### The methodological corollary — more important than the finding
**A waveform-only bracket manufactures the appearance of a step, and the spread
gate is blind to onsets by construction.**
The offset is what fires the trigger, so no waveform event can exist until the
ramp has nearly reached the trigger level. BE18438's first event of each episode
sits at 0.280 against a 0.300 trigger, and 0.185 against 0.200. At daily cadence
against a 3 h ramp, P(catching an intermediate) = **0.125**.
And `spread <= 0.02` rejects any record in which the floor is *moving* — which
is exactly what an onset is. **The gate rejected the very BE18438 record where
the ramp is visible.** If the operational goal is catching a fault early, before
the unit floods the store with junk events, the current detector is the wrong
shape for the job.
### The surviving shape
An **electrical, reversible, two-time-constant settling process** (~10 min and
~hours), saturating at a ceiling, with occasional sub-3-minute discrete jumps
superposed (BE18438 2026-02-26: 13:24 pre +0.180 / mid +0.240 / end +0.255 →
13:27 +0.325, identical metadata). That is the signature of a **bias or leakage
path charging a high-impedance node** — the class of fault Instantel's autozero
recovers ~10% of the time, and what the X1/X8 gains measure.
**It is a shape constraint, not a mechanism. Do not write it up as one.**
### Dead — with the evidence, so none of this is re-derived
| Killed | Evidence |
|---|---|
| **Latched step at onset** | >=25 one-minute intermediates over ~10 h, ASCII-validated. Direct observation, not a test. |
| **Slow degradation over days/weeks** | Same observation — bulk of the excursion in 7 min to 2.7 h. |
| **Thermal driving of pedestal magnitude** | BE13117, 365-count pedestal, n=128: full-day modulation **−0.42% ± 0.42%**, 95% CI [−1.25%, +0.40%]. Healthy-fleet seasonal zero drift totals **~0.3 A/D counts** — 15x to 1200x too small. Best-powered result in the campaign. |
| **Ground-motion shock** | 30-day window-max percentile ranks 0.03/0.98/0.15/0.01/0.68/0.24, median **0.194** against a null of 0.5. **0 of 7 events >=9 in/s** was followed by an onset within 30 d. BE12599 hit 10.220 in/s (2023-11) and 10.005 (2025-04) and did not onset until 2026-08-14. |
| **Handling / redeployment** | **0 of 9** onsets had a Project/Client/Seis Loc change. Widened to 30 d: 2 observed vs 4.90 expected, P(X>=2)=0.995 — *depleted*, the wrong direction. The apparent gap effect (p=0.035) died on histogram coverage: BE18438's "59.7-day gap" contains 122 histogram files; true silence 0.52 d. |
| **Mechanical resonance / damping change** | BE18438|Vert at a 64-count pedestal (3x outside Instantel's ±21): ΔTest-Freq **CI [−0.090, +0.021]** against 0.127 Hz for a real calibration. Block permutation p=0.658. |
| **Accumulated-duty threshold** | ~4 clean units logged more monitoring than the largest positive onset dose; BE18193 logged **13.45M intervals, 6.2x**. A counterexample — no power argument weakens it. |
| **Firmware** | **14,338 of 14,340** exports read `V 10.72-8.17`. A constant cannot explain a variable. |
| **Unit age** | Serial rank-sum 118.0 vs null 115.0, p=0.549; unchanged on the 8-unit re-cut (p=0.586). Serial is a poor age proxy anyway (Spearman +0.113 against archive entry). |
| **Strong seasonal clustering** | 25 onsets, exposure-weighted permutation **p=0.59**. Excludes >=75%-in-one-season only; a 2x seasonal hazard is *not* excluded. |
Also retire two overstated bounds. H6's dose-response exclusion "|r| > 0.03" is
a **10x overstatement** once clustering is corrected — the honest bound is
|r| > 0.1–0.3, so a real r=0.2 is not excluded. And **any statistic quoted
per-event**: 512 flagged channel-events collapse to **4.9 effective independent
observations** (unequal-cluster design effect 104.6 at ICC=1), and **55% of the
flagged corpus is one instrument on two calendar days** (BE13117, 2023-05-03/04).
### Power — read every negative in this section as bounded
Fisher exact, 5 positives of 45, one-sided α=0.05, exposure a third of the fleet:
| relative risk | power |
|---|---|
| 1.5 | 0.059 |
| 2 | 0.112 |
| 3 | 0.231 |
| 6 | 0.497 |
| 15 | 0.753 |
80% power needs **RR ≈ 13–20**. Even a *perfect* split reaches p<0.05 only if
the exposed group is <=25 of 45 units. **This archive can detect only
near-deterministic unit-level causes.** Every negative above excludes a strong
effect, not a real one.
### What this archive can NEVER answer
- **The A/D zero and the X1/X8 gains.** The 2027–2069 numbers appear in no file,
header or decoded record. They exist only on a live device behind `SUB 0x0E`.
Q1 is structurally unanswerable from data.
- **Unit-level vs component-level cause.** **Zero of 14,340** exports carry a
geophone or sensor serial. Q4 is dead — there is no way to know whether the
same physical geophone came back after service.
- **Service history.** The only service-adjacent field is `Calibration: <date>`
— 30 distinct dates fleet-wide, none before 2023, ASCII corpus entirely
2025–26. BE9558's 2020 and BE13117's 2023 episodes have no calibration record.
- **Temperature.** Zero exports carry it. Battery Level is a verified coarse
thermometer (+0.204 V winter over summer, 20/20 unit-years, p=9.5e−7, matching
lead-acid tempco) but quantised at 0.1 V ≈ 10 °C — useless within a day. The
archive can *bound* thermal; it can never *test* it.
- **BE13117 specifically** — 55% of the flagged corpus, the largest pedestal at
1.92 in/s, **zero** ASCII exports, histogram record ending eight months before
its episode. The most informative case in the archive is permanently outside
every metadata test.
- **The mild-offset rate**, and therefore the base rate's denominator. Event
files only see offsets large enough to dominate the trace.
### The experiment to run — `SUB 0x0E`, one afternoon
Point Blastware at `bridges/ach_mitm.py` and run **Unit Channel Test** against
(1) a faulting unit, (2) a known-good control, (3) the same unit before and
after an autozero. BW's sequence is `0x0E x8 → 0x98 x2 → 0x0E x8`, the second
pass carrying live ADC. Eight 10-byte payloads with expected values near 2048 is
a very constrained puzzle.
- **Proves:** whether the X1/X8 gains are readable over the wire, and whether
the fault sits at or upstream of the ADC zero reference. Gains walk out of
2027–2069 with the pedestal → the fault *is* the zero reference, Q1 answered.
Gains hold while the trace moves → the fault is downstream, look at the front
end.
- **§8c hands it a falsifiable time course:** poll at ~1-minute cadence and the
numbers should **ramp over minutes-to-hours, not step**. If they step while
the trace ramps, the two are decoupled.
- **Payoff:** converts the 10%/90% ship-it-or-not gamble into a decision made
before packing a box, remotely, for the whole fleet.
- ⚠ In the MITM topology filenames are reversed — `raw_s3_*.bin` holds
Blastware's bytes.
**Second: swap the geophone** between a faulted base and a healthy one. Fault
follows the sensor → element or cable. Fault stays with the base → front-end
board. One afternoon, zero code, and it settles the one question the archive is
permanently blind to.
**Third: log a faulting unit for 72 h untouched.** Every recovery we have is
confounded by a possible field autozero. A shelf and a logger settles whether
the fault genuinely self-reverses.
**Fourth, free: re-cut the fleet label** — drop the spread gate, re-score
amplitude-only, screen the 11 unscreened thin-coverage units. Might reach 9–10
positives. Be honest about the gain: power against "older half carries 3x the
hazard" rises only 0.23 → 0.30.
**Highest-value item overall, and not an experiment: the RMA/repair records.**
Which unit went back, when, what was done (autozero vs geophone replaced vs
board), and the geophone serial fitted. "Same channel after a documented
geophone *replacement*" is component-level-negative in one observation.
---
### 8d. The non-motion test — Brian's "it doesn't cross zero" (2026-09-07)
Looking at BE12599's 2026-08-09 event, Brian noted it reports no ZC frequency
**because the trace never crosses zero**. That observation is the best detector
in this investigation, and it comes from physics rather than a threshold.
A geophone is a velocity sensor with no DC response, so its output over a record
must integrate to ~zero — the ground does not relocate. Real motion therefore
sits roughly half below zero. Anything electrical is one-sided.
mp = |mean| / peak ~0 for motion, ~1 for a fault
frac_neg = share of samples < 0
`scratch/nonmotion_scan.py`, all 6,577 waveforms, 19,731 channel-rows.
Restricted to peak >= 0.05 in/s (n = 12,068), the distribution is **bimodal
with an empty middle**:
| mp band | channel-events |
|---|---|
| 0.0–0.1 | 11,384 |
| 0.1–0.2 | 293 |
| **0.15–0.85 (dead zone)** | **131 = 1.09%** |
| 0.9–1.0 | 278 |
At `mp >= 0.8` with >=3 events it returns **exactly the five confirmed units** —
BE9558, BE11529, BE12599, BE13117, BE18438 — stable from 0.5 to 0.9. Two
detectors on entirely different principles agreeing on the unit list is the
strongest corroboration that list has.
**BE11007 is settled: NOT an offset.** It reaches mp 0.75–0.89, but with
`frac_neg = 0.99` at peaks of **7.4–9.4 in/s** — parked *negative* during a
near-full-scale blast. §4's guess was right. `mp` alone cannot separate a
pedestal from a large one-sided blast; pair it with a peak ceiling or with
sign-consistency across events.
⚠ **Not a rediscovery of the retracted v1 detector.** v1 scored only the
largest-peak axis and used the mean as a *baseline estimator* where the median
was required. Here the mean is the signal itself, per channel — that is what the
physics licenses.
**Correction to §8c.** That section says the spread gate is "blind to onsets by
construction." Too strong: of 87 BE18438|Vert events at mp >= 0.5 the gate
rejected **one** — the transitional record. It does not lose onsets
systematically; it loses the transition specifically.
### 8e. BE12599 — a connector, not a geophone (2026-09-07)
Waveform shapes across its August episode, measured rather than eyeballed:
| date | channel | shape |
|---|---|---|
| Aug 09 05:29 | Long | **unipolar +**, 0/2304 samples below zero, decay tau **26 ms** |
| Aug 09 05:35 | Long | unipolar +, 3 spikes at irregular gaps (744, 1032 ms), tau **38 ms** |
| Aug 14 05:00 | Long | single lobe, bipolar, tau **118 ms** |
| Aug 17–23 | Tran | **flat DC pedestal**, sd/level 0.015–0.020, 0 zero crossings |
**Unipolar impulses with an RC tail are not mechanical.** Fast rise, exponential
decay, one polarity, irregular timing — that is charge dumped into a
capacitively-coupled input and draining through the input resistance. The
progression 26 ms -> 118 ms -> never recovers, over 14 days, is a leakage path
worsening.
**And the fault moved channels** — Long on Aug 9/14, Tran on Aug 17–23, Long
again on Aug 21 (1.065 in/s) while Tran held its pedestal. Vert stayed clean
throughout. **A failing geophone element cannot hop channels. A connector can.**
That single fact explains what had been puzzling:
- **The sensor self-check keeps passing** (7.4/7.5/7.6 Hz, ratios 3.6–4.2, all
four channels Passed, on the very events where Long throws 0.5 in/s spikes).
The swing test drives the element; the element is fine. The fault is in the
wiring to it.
- **Why Instantel's autozero fixes only ~10%** — it cannot fix a connector.
- **Why onset "ramps" over minutes to hours** — contact resistance drifting.
All seven Aug 17–23 events are stamped **05:00:14**, the same second, and their
filename extensions run `8E → WE → KE → 8E → WE → KE → 8E` — the documented
3-day cycle for a fixed daily time. Clock-scheduled, not physically triggered:
the modem powers up, draws a surge, and a marginal connection responds.
**Field action: inspect and photograph the geophone connector BEFORE reseating
anything** — an intermittent contact clears the moment it is disturbed.
⚠ Scoped to BE12599. BE18438's onset was a smooth 7-minute ramp with no spikes,
which looks like a different failure mode wearing the same signature.
---
### ⚠ Serial prefixes — four of these units are BlastMates, not MiniMates
Corrected 2026-09-06, after Brian queried "BA10895?" against a report that
said BE10895. He was right. The BW filename encodes the serial **number
only** — `L895` -> 10895 — and every offset scanner synthesised the family
prefix as `"BE"`. Four of the 43 archive units are **BA** (BlastMate, the
MiniMate Plus's bigger sibling; same Series III, byte-identical data):
**BA9229, BA10060, BA10895, BA15957.**
Read off the file bodies, which carry the serial verbatim. No analysis
changed — grouping was always on the numeric part, and no unit number maps
to two serials — but every earlier reference to "BE10895" and the other
three is a label error and has been corrected throughout this document.
The same assumption was live in two production sites and is fixed
(`sfm/waveform_store.py`, `minimateplus/client.py`): the store would have
filed a BlastMate under a unit that does not exist, and the monitor-log
decoder lost the geo threshold along with the serial. See commit `9ceff65`.
---
## 9. Chronology
| date | event |
@@ -973,15 +512,3 @@ decoder lost the geo threshold along with the serial. See commit `9ceff65`.
| 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 | **BA10895 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. |
| 2026-09-06 | **Four units relabelled BA, not BE** — BA9229, BA10060, BA10895, BA15957 are BlastMates. The BW filename carries only the serial number; the family prefix must be read from the file body. Fixed in the scanners and in two production sites. |
| 2026-09-06 | **Mechanism campaign — five hypotheses, all dead.** Thermal, ground-motion shock, handling/redeployment, accumulated duty, unit age, firmware and a mechanical element fault are each refuted or bounded. 7 nominally significant results against 5.2 expected by chance. |
| 2026-09-06 | **Onset is a RAMP of minutes-to-hours, not a step** — BE18438 Vert resolved to one-minute cadence, 50% of the excursion in 7 min, >=25 intermediates, ASCII-validated 75/75. Kills both a latched digital step AND slow component degradation. Surviving shape: a reversible two-time-constant settling process — a bias/leakage path charging a high-impedance node. |
| 2026-09-06 | **Polarity consistency RETRACTED** (a tautology of the spread gate; amplitude-only gives 12 of 53 unit-channels mixed) and the fleet **re-cut to 8–9 units, not 5**. "Persistent until serviced" weakened: 23 recoveries, median gap 6 days — though a field autozero cannot be excluded. |
| 2026-09-06 | The spread gate is **blind to onsets by construction** — it rejects a moving floor, which is what an onset is. It rejected the very record in which the ramp is visible. |
| 2026-09-07 | **The non-motion test** (Brian: "it doesn't cross zero"). `\|mean\|/peak` is bimodal with a 1.09% dead zone and returns exactly the 5 confirmed units from physics, not a threshold. Independent corroboration of the unit list. **BE11007 settled as NOT an offset** — a one-sided 9 in/s blast. |
| 2026-09-07 | **BE12599 is a connector fault, not a geophone fault.** Unipolar spikes with a 26→118 ms RC tail progressing to a flat pedestal, and the fault MOVES between Long and Tran while the sensor self-check passes on every event. An element cannot hop channels; a connector can. Inspect before reseating. |
+40 -221
View File
@@ -47,24 +47,19 @@ from dataclasses import dataclass
from pathlib import Path
from typing import Optional, Union
# Thor IDFW bodies use the series-3 record-chain decoder.
# Thor IDFW bodies are pinned to the SUPERSEDED tag-dispatch decoder.
#
# This was previously pinned to the SUPERSEDED tag-dispatch walker
# (`decode_waveform_legacy`) on the stated grounds that "Thor has no ASCII
# ground truth in the corpus and its geo scaling is separately suspect".
# Both premises were false: Thor writes a per-sample CSV export next to every
# binary (see scratch/verify_thor_against_csv.py), and the scaling is now
# resolved (see _GEO_LSB_IPS). Measured against that ground truth on
# 2026-09-10, the record chain beats the legacy walker outright:
#
# channel truncation 55/153 files -> 3/153
# files exact 98/153 -> 150/153
# per-sample exact 99.781% -> 99.854%
#
# The legacy walker stops at the first unrecognised tag and returns whatever
# channels it had, so its failure mode is silent short channels rather than an
# error. Do not re-pin it.
from minimateplus.waveform_codec import _MODES, decode_waveform_v2, is_record
# _find_waveform_body_offset() trial-decodes every candidate offset and keeps
# whichever yields the most samples. The series-3 record-chain decoder
# correctly returns None where the legacy walker returned garbage, which
# changes that heuristic's winner on 33 of 577 files. The net effect measured
# 2026-08-25 was positive (all-channels-equal 8/577 -> 506/577, mean abs PPV
# error 0.228 -> 0.173 in/s) but Thor has no ASCII ground truth in the corpus
# and its geo scaling is separately suspect, so the switch is deferred until
# the body-offset search is reworked to use the record chain directly.
from minimateplus.waveform_codec import (
decode_waveform_legacy as decode_waveform_v2,
)
from .models import IdfEvent, IdfPeaks, IdfReport
@@ -94,70 +89,23 @@ _BODY_MAGIC = b"\x00\x02\x00"
# fixed-header region where the same magic legitimately appears inside
# channel-test records and the compliance block (offsets 0x015d, 0x091c,
# 0x0ae2, 0x0d30 in observed events).
# Lowered from 0x0E00 to 0x0C00 (2026-09-10). Three-channel events -- mic
# disabled -- have a shorter fixed header and put their record chain head at
# 0x0dba, below the old floor. The head was therefore invisible to the scan,
# which fell through to the *Vert* segment-0 record and decoded a body shifted
# one position around the channel rotation. 46 of 139 files in the
# 9-10-26-csv-req corpus were affected; all 46 became per-sample exact once
# the head was reachable. The floor still skips the fixed-header region,
# where `is_record()` can match channel-test records (0x015d, 0x091c, 0x0ae2).
_BODY_SCAN_FLOOR = 0x0C00
_BODY_SCAN_FLOOR = 0x0E00
# Cap on trial decodes per file. Chain-head detection normally yields one
# or two candidates; the cap only bounds the worst case on a corrupt file.
_MAX_BODY_CANDIDATES = 16
# Geophone count → in/s.
#
# The old value 0.0003 was read off the smallest non-zero sample in the
# sidecar corpus, but that sample is Thor's *4-decimal display rounding* of
# the true LSB, not the LSB itself. It read every series-4 geophone sample
# 3.3% low. The quantisation ladder gives it away: counts 1..6 export as
# 0.0003, 0.0006, 0.0009, 0.0012, 0.0016, 0.0019 — an LSB of exactly 0.0003
# would end 0.0015, 0.0018.
#
# The value below maximises exact 4-dp agreement over 1,046,016 paired
# samples (454 channel-events, 2 units) at 99.854%, versus 50.7% for 0.0003.
# It is a global constant, not a per-unit calibration: all 8 UM units in the
# production store independently agree to within ±0.07% on their
# device-reported PPV. 1/LSB = 3222.6 counts per in/s.
#
# The value is pinned, not guessed. Each exported sample constrains the LSB
# to the window that rounds to the printed 4-dp figure; intersecting 991,415
# such constraints (clean channel-events only) gives
#
# LSB in [0.000310307933, 0.000310308057] width 1.2e-10
#
# 0.000310308 sits at the centre of that window. Equivalent full scale is
# 10.0 in/s / 0.000310308 = 32226.05 counts.
#
# Corroboration from the device: an IDFH interval that never recorded keeps
# its min/max accumulator at its ±full-scale seed, and that seed is
# (min=+32226, max=-32226) — the same magnitude, independently. Note the
# tempting closed form 10.0/32226 is very slightly WRONG: it lands 4.5e-10
# above the feasible window and loses 78 boundary samples to the literal
# value while never winning one. Series-3 uses 32000 counts for the same
# 10.0 in/s, so the two generations do NOT share a scale.
#
# Ground truth + harness: scratch/verify_thor_against_csv.py
_GEO_LSB_IPS = 0.000310308
# Geophone count → in/s, derived from sidecar ground truth: the smallest
# non-zero sample in 1,014-file corpus is 0.0003 in/s.
_GEO_LSB_IPS = 0.0003
# Microphone count → psi, derived from sidecar regression on 50 sample
# pairs from UM11719_20231219162723.IDFW (mic-heavy event).
_MIC_LSB_PSI = 2.14e-6
# IDFH histogram constants.
# Bytes per interval record = 16 per channel + an 8-byte tail, so a
# 4-channel unit uses 72 and a mic-disabled 3-channel unit uses 56. It is
# NOT a constant: derive it per segment from the interval counter (see
# decode_idfh_body). This value survives only as the 4-channel default.
_IDFH_INTERVAL_SIZE = 72 # bytes per per-interval record (4 channels)
_IDFH_CHANNEL_BLOCK = 16 # bytes per channel inside an interval record
_IDFH_INTERVAL_TAIL = 8 # bytes after the per-channel blocks
_IDFH_INTERVAL_SIZE = 72 # bytes per per-interval record
_IDFH_SEGMENT_HEADER = 10 # bytes: [len_be 2B][0a 00 00 00 4B][00 NN 2B][05 3f 2B]
_IDFH_SEGMENT_TAIL = 2 # bytes after the interval data block, before next marker
_IDFH_HALFP_FREQ_NUM = 512.0 # freq_hz = NUM / halfp; halfp ≤ 5 means ">100 Hz" sentinel
_IDFH_GEO_FULL_SCALE = 10.0 # in/s — Normal range
_IDFH_INT16_FS = 32768.0
_IDFH_CHANNELS = ("Tran", "Vert", "Long", "MicL")
@@ -275,67 +223,26 @@ def _find_waveform_body_offset(buf: bytes) -> Optional[int]:
"""
if len(buf) < _BODY_SCAN_FLOOR + 8:
return None
# 1. Locate every plausible per-channel record header. A header carries
# [len 2B][channel_id][00][00] at +2..+6, so anchor the search on the
# three-byte ``<cid> 00 00`` signature and validate with is_record().
# Scanning candidate *preambles* instead is not viable: MODE_RAW16 is
# ``00 00``, so every run of three zero bytes would look like a body
# start and each would cost a full trial decode (~0.5 s/file measured).
floor = max(0, _BODY_SCAN_FLOOR - 7)
starts: list = []
for cid in (0x46, 0x47, 0x48, 0x49):
sig = bytes((cid, 0x00, 0x00))
i = floor
best: Optional[tuple[int, int]] = None # (total_samples, offset)
i = _BODY_SCAN_FLOOR
while True:
j = buf.find(sig, i)
j = buf.find(_BODY_MAGIC, i)
if j < 0:
break
i = j + 1
q = j - 4
if q >= floor and is_record(buf, q):
starts.append(q)
if not starts:
return None
starts.sort()
# 2. A body begins at the head of a record chain -- a record that no other
# record's length field points at. The head's own payload is the
# implicit segment-0 Tran record, and the body offset is head + 7 (past
# [len 2B][cid][00][00][seg]) so that body[1:3] lands on the mode.
ends = {q + 2 + int.from_bytes(buf[q + 2 : q + 4], "big") for q in starts}
heads = [q for q in starts if q not in ends] or starts[:1]
# 3. Trial-decode each head and keep the best. Prefer a candidate where
# all four channels come out the same length: scoring on raw sample
# count alone picks false positives sitting *inside* a record header,
# which decode a plausible-looking but rotation-shifted body that
# silently drops each channel's segment 0.
best = None
best_off = None
for head in heads[:_MAX_BODY_CANDIDATES]:
j = head + 7
if j + 3 > len(buf) or (buf[j + 1], buf[j + 2]) not in _MODES:
continue
try:
decoded = decode_waveform_v2(buf[j:])
except Exception:
continue
if not decoded:
continue
lengths = [len(v) for v in decoded.values() if v]
total = sum(len(v) for v in decoded.values())
# A "real" body has more than just the 2-sample preamble.
if total <= 2:
continue
# >= 3 rather than == 4: a mic-disabled event has only the three geo
# channels, and demanding four made `equal` permanently False for
# them, leaving the pick to raw sample count alone.
equal = len(lengths) >= 3 and len(set(lengths)) == 1
score = (equal, total)
if best is None or score > best:
best, best_off = score, j
return best_off
if best is None or total > best[0]:
best = (total, j)
return best[1] if best else None
def _decode_waveform_samples(buf: bytes) -> Optional[dict]:
@@ -392,12 +299,6 @@ class IdfhInterval:
micl_min: int
micl_max: int
micl_halfp: int
# 4 on a normal unit; 3 when the microphone is disabled, in which case the
# micl_* fields are absent from the record and read as zero.
n_channels: int = 4
def has_channel(self, channel: str) -> bool:
return channel != "MicL" or self.n_channels >= 4
def peak_count(self, channel: str) -> int:
mn = getattr(self, f"{channel.lower()}_min")
@@ -406,11 +307,7 @@ class IdfhInterval:
def peak_ips(self, channel: str) -> float:
"""Convert peak count to in/s (geo channels only)."""
# Same geo LSB as the waveform path — verified independently against
# the IDFH exports: as peak magnitude rises (and 4-dp quantisation
# noise falls) the implied LSB converges on 0.0003103, matching
# _GEO_LSB_IPS. The old 10.0/32768 read histogram peaks 1.7% low.
return self.peak_count(channel) * _GEO_LSB_IPS
return self.peak_count(channel) / _IDFH_INT16_FS * _IDFH_GEO_FULL_SCALE
def freq_hz(self, channel: str) -> Optional[float]:
halfp = getattr(self, f"{channel.lower()}_halfp")
@@ -419,46 +316,11 @@ class IdfhInterval:
return _IDFH_HALFP_FREQ_NUM / halfp
def _is_unwritten_interval(interval: "IdfhInterval") -> bool:
"""True for an interval slot the device reserved but never wrote.
Thor seeds each interval's per-channel accumulators at ``min = +full
scale`` and ``max = -full scale`` and then narrows them as samples
arrive. A slot that never recorded keeps that seed, so ``min > max`` —
impossible for real data. Such a record decodes to a full-scale
10.0 in/s peak on every channel and, being a max-over-intervals, poisons
the whole file's PPV.
Rare but real: exactly 1 of 497,611 corpus intervals, and it inflated
that file's Long PPV from 0.0081 to 10.0 in/s. The inversion is always
all-or-nothing across channels (0 partial cases in the corpus), so
requiring every channel to be inverted keeps this from ever firing on
genuine data.
"""
pairs = [
(interval.tran_min, interval.tran_max),
(interval.vert_min, interval.vert_max),
(interval.long_min, interval.long_max),
]
if interval.has_channel("MicL"):
pairs.append((interval.micl_min, interval.micl_max))
return all(mn > mx for mn, mx in pairs)
def _decode_idfh_interval(buf72: bytes, offset: int,
n_channels: int = 4) -> IdfhInterval:
"""Decode one interval record into per-channel min/max/halfp.
The record is ``n_channels`` × 16-byte blocks plus an 8-byte tail, so it
is 72 bytes on a normal unit and 56 when the microphone is disabled.
Missing channels read as zero.
"""
def _decode_idfh_interval(buf72: bytes, offset: int) -> IdfhInterval:
"""Decode one 72-byte interval record into per-channel min/max/halfp."""
import struct
fields = []
for i in range(4):
if i >= n_channels:
fields.extend([0, 0, 0])
continue
block = buf72[i * 16 : (i + 1) * 16]
mn = struct.unpack_from(">h", block, 0)[0]
mx = struct.unpack_from(">h", block, 2)[0]
@@ -474,7 +336,6 @@ def _decode_idfh_interval(buf72: bytes, offset: int,
vert_min=fields[3], vert_max=fields[4], vert_halfp=fields[5],
long_min=fields[6], long_max=fields[7], long_halfp=fields[8],
micl_min=fields[9], micl_max=fields[10], micl_halfp=fields[11],
n_channels=n_channels,
)
@@ -482,73 +343,36 @@ def decode_idfh_body(buf: bytes) -> list:
"""Walk an IDFH file and decode every interval record.
The body has one or more segments; each segment header is 12 bytes:
``[length_be 2B][0a 00 00 00][counter_be 2B][05 3f]`` where ``length``
``[length_be 2B][0a 00 00 00][00 NN_counter][05 3f]`` where ``length``
is bytes from the magic through the end of the interval block
(= 10 + 72 × n_intervals). Segments are separated by a 2-byte tail
+ next-segment 2-byte prefix (the bytes before the next length field).
``counter`` is a **uint16 BE cumulative interval index** — the 0-based
index of the LAST interval in this segment. Segments carry 10
intervals each, so it runs 9, 19, 29, ... across the file.
⚠ This validator used to require ``buf[j + 4] == 0x00``, i.e. that the
counter's high byte was zero. That silently capped every histogram at
**250 intervals**: the moment the cumulative counter passed 255 the high
byte went non-zero and every later segment was rejected, so any
monitoring run longer than ~4 hours lost its tail — frequently the part
holding the event peak, which is why those files' PPV read low. 540 of
858 corpus files were affected. Do not reinstate that check.
Confirmed against the 859-file corpus (181,071 intervals decoded; 1
failure is the sig-B BE9439 file).
"""
intervals: list = []
i = 0
prev_counter = -1 # so the first segment's n = counter + 1
while True:
j = buf.find(b"\x0a\x00\x00\x00", i)
if j < 0 or j < 2:
break
# Validate: [length_be][0a 00 00 00][counter_be][05 3f]. The counter
# is deliberately NOT constrained — see the note above.
if buf[j + 6 : j + 8] != b"\x05\x3f":
# Validate: [length_be][0a 00 00 00][00 NN][05 3f]
if buf[j + 4] != 0x00 or buf[j + 6 : j + 8] != b"\x05\x3f":
i = j + 1
continue
length = int.from_bytes(buf[j - 2 : j], "big")
counter = int.from_bytes(buf[j + 4 : j + 6], "big")
header_start = j - 2
if length < _IDFH_SEGMENT_HEADER or header_start + length > len(buf):
# Truncated / bogus length — not a real segment header.
i = j + 1
continue
# The counter is the cumulative index of this segment's LAST interval,
# so the interval count is its delta from the previous segment. That
# gives the record stride, which is NOT fixed: 16 bytes per channel
# plus an 8-byte tail, so 72 for a 4-channel unit and 56 for a
# mic-disabled 3-channel one. Assuming 72 unconditionally made every
# 3-channel histogram read 7 intervals per 10-interval segment,
# walking off alignment into garbage that decoded as ~10 in/s peaks.
n = counter - prev_counter
n = (length - _IDFH_SEGMENT_HEADER) // _IDFH_INTERVAL_SIZE
if n <= 0:
i = j + 1
continue
stride = (length - _IDFH_SEGMENT_HEADER) // n
n_channels, remainder = divmod(stride - _IDFH_INTERVAL_TAIL,
_IDFH_CHANNEL_BLOCK)
if remainder or not (1 <= n_channels <= 4):
i = j + 1
continue
header_start = j - 2
interval_start = header_start + _IDFH_SEGMENT_HEADER
for k in range(n):
off = interval_start + k * stride
if off + stride > len(buf):
off = interval_start + k * _IDFH_INTERVAL_SIZE
if off + _IDFH_INTERVAL_SIZE > len(buf):
break
chunk = buf[off : off + stride]
interval = _decode_idfh_interval(chunk, off, n_channels)
if _is_unwritten_interval(interval):
# Reserved-but-never-recorded slot: the min/max accumulators
# still hold their ±full-scale seed. Counting it would
# fabricate a 10.0 in/s peak on every channel.
continue
intervals.append(interval)
prev_counter = counter
chunk = buf[off : off + _IDFH_INTERVAL_SIZE]
intervals.append(_decode_idfh_interval(chunk, off))
# Advance past this segment + the 2-byte tail.
i = header_start + length + _IDFH_SEGMENT_TAIL
return intervals
@@ -628,12 +452,7 @@ def read_idf_file(
peak_long = max((iv.peak_ips("Long") for iv in intervals), default=0.0)
# Mic peak in psi — Thor stores per-interval mic ADC counts in the
# binary; convert the max count to psi via the per-count factor.
# Skip on a mic-disabled (3-channel) unit: those records carry no mic
# block at all, so peak_count("MicL") would report a synthetic zero.
mic_peak_count = max(
(iv.peak_count("MicL") for iv in intervals if iv.has_channel("MicL")),
default=0,
)
mic_peak_count = max((iv.peak_count("MicL") for iv in intervals), default=0)
mic_peak_psi = mic_count_to_psi(mic_peak_count) if mic_peak_count else None
rep = IdfReport(
serial_number=md.serial,
+1 -9
View File
@@ -30,7 +30,6 @@ from __future__ import annotations
import datetime
import logging
import re
import struct
from typing import Optional
@@ -2533,17 +2532,10 @@ def _decode_0a_partial_header(raw_data: bytes, index: int, key4: bytes) -> Optio
ts2 = try_ts(raw_data[ts1_end + 1:ts1_end + 1 + ts_size])
# Extract serial and geo threshold from "BE11529\0" and "Geo: X.XXX in/s\0".
#
# Match any two-letter family prefix, not a literal "BE" — a BlastMate
# reports "BA10895", and the old `find(b"BE")` returned -1 on one. That
# skipped this whole block, so the geo threshold went missing along with
# the serial. Requiring the NUL terminator in the pattern also makes the
# match stricter than the bare two-byte search it replaces.
serial: Optional[str] = None
geo_ips: Optional[float] = None
serial_match = re.search(rb"[A-Z]{2}\d{3,6}(?=\x00)", raw_data)
serial_pos = serial_match.start() if serial_match else -1
serial_pos = raw_data.find(b"BE")
if serial_pos >= 0:
# Read null-terminated serial starting at serial_pos.
null_pos = raw_data.find(b"\x00", serial_pos)
+1 -1
View File
@@ -50,7 +50,7 @@ SIDECAR_KIND = "sfm.event"
# bumped without a `pip install` re-run — leading to confusing stale
# version stamps in sidecars. Bump this constant and CHANGELOG.md
# together at release time.
TOOL_VERSION = "0.30.0"
TOOL_VERSION = "0.27.0"
try:
# Best-effort: prefer the installed metadata when it's NEWER than the
+7 -44
View File
@@ -722,18 +722,7 @@ STREAM_END_ID = 0x06
MODE_DELTA = (0x02, 0x00)
MODE_ABSOLUTE = (0x01, 0x00)
MODE_RAW12 = (0x00, 0x03)
# Raw int16 BE absolute samples, 10-byte header, no tags — the same shape as
# MODE_RAW12 but two bytes per sample instead of 1.5. Found on Thor/Micromate
# segment-0 records (2026-09-10): a `len=1032` record carries exactly
# (1032 - 8) / 2 = 512 samples and reproduces Thor's own export 512/512
# exactly. Before this mode existed the record fell through the dispatch
# unhandled, so the channel silently lost its first 512 samples.
MODE_RAW16 = (0x00, 0x00)
_MODES = (MODE_DELTA, MODE_ABSOLUTE, MODE_RAW12, MODE_RAW16)
# Preambles whose leading data is untagged and therefore cannot be
# block-walked; find_first_record() must scan for the next record instead.
_UNTAGGED_MODES = (MODE_RAW12, MODE_RAW16)
_MODES = (MODE_DELTA, MODE_ABSOLUTE, MODE_RAW12)
def _u16(b: bytes, p: int) -> int:
@@ -758,18 +747,7 @@ def data_block_len(body: bytes, p: int) -> Tuple[Optional[int], Optional[int]]:
hi = t0 & 0xF0
nn = ((t0 & 0x0F) << 8) | t1
if hi == 0x40: # int16 BE data block
# NN was capped at 0x08 until 2026-09-11. That cap had no basis: the
# two corpora available at the time only ever used NN in {1,2,3,4,8},
# so it was never exercised. Loud UM12947 events use NN of 12, 16,
# 20 ... up to 196, and every value above 8 halted the walk, which
# surfaced as silently short channels (walk_body/run stop at the first
# unrecognised tag rather than raising). Verified against Thor's own
# exports: 22 length-mismatched files -> 0, and the affected corpus
# went to 1,476,242/1,476,249 samples exact. The real bound is the
# buffer; the caller additionally clamps to the record end.
if nn == 0 or p + 2 * nn + 2 > len(body):
return None, None
return 2 * nn + 2, nn
return (None, None) if (nn == 0 or nn > 0x08) else (2 * nn + 2, nn)
if nn == 0 or nn % 4:
return None, None
if hi == 0x00:
@@ -783,11 +761,6 @@ def data_block_len(body: bytes, p: int) -> Tuple[Optional[int], Optional[int]]:
return None, None
def unpack16(data: bytes) -> List[int]:
"""Raw int16 BE absolute samples (MODE_RAW16)."""
return [_i16(data, 2 * k) for k in range(len(data) // 2)]
def unpack12(data: bytes) -> List[int]:
"""Raw 12-bit packed samples: 6 bytes -> 4 signed values."""
out: List[int] = []
@@ -812,17 +785,13 @@ def find_first_record(body: bytes) -> Optional[int]:
"""Offset of the first record, or None.
Under the normal ``00 02 00`` preamble the leading bytes are segment-0's
Tran blocks, so walk them. Under the untagged preambles (``00 00 03``
raw-12 and ``00 00 00`` raw-16) that data has no tags at all and cannot
be block-walked — scan for the next record header instead.
Tran blocks, so walk them. Under the ``00 00 03`` preamble that data is
raw 12-bit with no tags at all and cannot be block-walked — scan instead.
"""
if len(body) >= 3 and (body[1], body[2]) in _UNTAGGED_MODES:
if len(body) >= 3 and (body[1], body[2]) == MODE_RAW12:
scan_from = 3
else:
# Tagged preamble. MODE_DELTA carries a 14-byte record header (two
# int16 anchors), so its blocks start at body[7]; MODE_ABSOLUTE has a
# 10-byte header and starts at body[3].
i = 3 if (len(body) >= 3 and (body[1], body[2]) == MODE_ABSOLUTE) else 7
i = 7
while i < len(body):
if is_record(body, i):
nxt = i + 2 + _u16(body, i + 2)
@@ -881,7 +850,7 @@ def decode_waveform_v2(body: bytes) -> Optional[dict]:
if len(body) < 8 or body[0] != 0x00:
return None
preamble = (body[1], body[2])
if preamble not in (MODE_DELTA, MODE_ABSOLUTE, MODE_RAW12, MODE_RAW16):
if preamble not in (MODE_DELTA, MODE_RAW12):
return None
first = find_first_record(body)
if first is None:
@@ -926,10 +895,6 @@ def decode_waveform_v2(body: bytes) -> Optional[dict]:
if preamble == MODE_DELTA:
out["Tran"].extend([_i16(body, 3), _i16(body, 5)])
run("Tran", 7, first, absolute=False)
elif preamble == MODE_ABSOLUTE:
run("Tran", 3, first, absolute=True)
elif preamble == MODE_RAW16:
out["Tran"].extend(unpack16(body[3:first]))
else:
out["Tran"].extend(unpack12(body[3:first]))
@@ -943,6 +908,4 @@ def decode_waveform_v2(body: bytes) -> Optional[dict]:
run(ch, off + 10, end, absolute=True)
elif mode == MODE_RAW12:
out[ch].extend(unpack12(body[off + 10:end]))
elif mode == MODE_RAW16:
out[ch].extend(unpack16(body[off + 10:end]))
return out
+1 -1
View File
@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
[project]
name = "seismo-relay"
version = "0.30.0"
version = "0.27.0"
description = "Python client and REST server for MiniMate Plus seismographs"
requires-python = ">=3.10"
dependencies = [
-91
View File
@@ -1,91 +0,0 @@
#!/usr/bin/env python3
"""Detect NON-MOTION on a geophone channel: |mean| / peak.
A geophone is a velocity sensor with no DC response, so its output over a
record must integrate to ~zero — the ground does not relocate. Real motion
therefore sits roughly half above and half below zero. Anything electrical —
a charge-injection spike, a step, a parked pedestal — is one-sided.
mp = |mean| / peak ~0 for motion, ~1 for a pedestal
frac_neg = share of samples < 0 ~0.3-0.5 for motion, ~0 for a fault
Why this beats the pre-trigger floor (`offset_scan3.py`): that detector's
`spread <= 0.02` gate rejects any record whose floor is MOVING, which is
exactly what an onset is — it discarded the one BE18438 record in which the
ramp was visible. This test is indifferent to whether the fault is a spike,
a ramp or a flat pedestal; none of them cross zero.
⚠ Not a rediscovery of the retracted v1 detector. v1 scored only the
largest-peak axis and used the mean as a BASELINE estimator, where the median
was required. Here the mean is the signal itself, per channel, and that is
what the physics licenses.
"""
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
_WAVE=re.compile(r"\.[A-Za-z0-9]{2}0[Ww]$"); _STEM=re.compile(r"^([B-Z])(\d{3})")
_SER=re.compile(rb"[A-Z]{2}\d{3,6}")
def serial_of(name, path=None):
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 _SER.findall(Path(path).read_bytes()):
s=s.decode()
if s[2:].lstrip("0")==str(num): return s
except Exception: pass
return f"BE{num}"
def scan(ps):
import logging; logging.disable(logging.WARNING)
p=Path(ps)
try: ev=read_blastware_file(p)
except Exception: return None
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 ""
ser=serial_of(p.name,p); out=[]
for ch in GEO:
a=[x*K for x in s[ch]]
pk=max(abs(x) for x in a)
if pk<=0: continue
out.append({"serial":ser,"timestamp":stamp,"filename":p.name,"channel":ch,
"peak":round(pk,4),
"mean":round(statistics.fmean(a),4),
"mp":round(abs(statistics.fmean(a))/pk,4),
"frac_neg":round(sum(1 for x in a if x<0)/len(a),4),
"n":len(a)})
return out
COLS=["serial","timestamp","filename","channel","peak","mean","mp","frac_neg","n"]
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)
a=ap.parse_args()
seen=set(); files=[]
for q in sorted(Path(a.dir).rglob("*")):
if q.is_file() and _WAVE.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 i,f in enumerate(as_completed([ex.submit(scan,p) for p in files]),1):
r=f.result()
if r: rows.extend(r)
if i%1000==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)
print(f"\nwrote {a.out} ({len(rows)} channel-rows)")
if __name__=="__main__": main()
-234
View File
@@ -1,234 +0,0 @@
#!/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()
+4 -26
View File
@@ -28,31 +28,9 @@ from minimateplus.event_file_io import read_blastware_file
GEO=("Tran","Vert","Long"); K=10.0/32000.0
_WAVE=re.compile(r"\.[A-Za-z0-9]{2}0[Ww]$"); _STEM=re.compile(r"^([B-Z])(\d{3})")
_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 serial_of(n):
m=_STEM.match(n)
return f"BE{(ord(m.group(1))-ord('B'))*1000+int(m.group(2))}" if m else "?"
def scan(ps):
p=Path(ps)
@@ -70,7 +48,7 @@ def scan(ps):
mid, end = a[t:2*t], a[2*t:]
if not pre or not mid or not end: continue
v=[statistics.median(x)*K for x in (pre,mid,end)]
out.append({"serial":serial_of(p.name, p),"timestamp":stamp,
out.append({"serial":serial_of(p.name),"timestamp":stamp,
"filename":p.name,"channel":ch,
"pretrig_n": pre_n or 0,
"pre":round(v[0],4),"mid":round(v[1],4),"end":round(v[2],4),
-228
View File
@@ -1,228 +0,0 @@
#!/usr/bin/env python3
"""Verify the Thor / Micromate (series-4) IDF decoder against Thor's own exports.
Sister harness to ``scratch/verify_against_ascii.py`` (series-3 / Blastware).
Ground truth is the ``.IDFW.csv`` / ``.IDFH.csv`` file Thor writes next to each
binary, under a sibling ``CSV/`` directory:
<dir>/UM13981_20220207084555.IDFW
<dir>/CSV/UM13981_20220207084555.IDFW.csv
For waveforms the CSV carries a per-sample block of four columns
(Tran, Vert, Long, Mic) in in/s and psi -- i.e. true per-sample ground truth,
exactly what the BW ASCII exports give us for series-3. The leading 2-column
rows are the report header (PPV, sample rate, geo range, ...).
Usage:
python scratch/verify_thor_against_csv.py [--root DIR] [--lsb FLOAT]
[--limit N] [--kind idfw|idfh|both]
"""
from __future__ import annotations
import argparse
import csv
import os
import statistics
import sys
from collections import Counter, defaultdict
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
from micromate import idf_file as M
DEFAULT_ROOT = "/home/serversdown/thor-watcher/example-data"
GEO = ("Tran", "Vert", "Long")
def parse_export(path):
"""Return (header_dict, sample_rows) from a Thor CSV export."""
hdr, rows = {}, []
with open(path, newline="", encoding="utf-8", errors="replace") as fh:
for rec in csv.reader(fh):
if len(rec) == 2:
hdr[rec[0].strip()] = rec[1].strip()
elif len(rec) >= 3:
try:
rows.append([float(x) for x in rec])
except ValueError:
pass
return hdr, rows
def index_corpus(root):
"""Map BASENAME.IDFW -> (binary_path, csv_path) for every paired file."""
exports, binaries = {}, {}
for dirpath, _dirs, files in os.walk(root):
for name in files:
up = name.upper()
full = os.path.join(dirpath, name)
if up.endswith(".IDFW.CSV") or up.endswith(".IDFH.CSV"):
exports.setdefault(name[:-4].upper(), full)
elif up.endswith(".IDFW") or up.endswith(".IDFH"):
binaries.setdefault(up, full)
return {k: (binaries[k], exports[k]) for k in binaries.keys() & exports.keys()}
def hdr_float(hdr, key):
raw = hdr.get(key)
if not raw:
return None
try:
return float(raw.split()[0])
except (ValueError, IndexError):
return None
def verify_waveform(binpath, csvpath, lsb):
"""Compare one IDFW against its export. Returns a result dict."""
out = {"file": os.path.basename(binpath), "status": "ok"}
try:
res = M.read_idf_file(binpath)
except NotImplementedError:
out["status"] = "not-thor"
return out
except Exception as exc: # noqa: BLE001 - harness reports, never raises
out["status"] = "decode-error"
out["detail"] = f"{type(exc).__name__}: {exc}"
return out
hdr, rows = parse_export(csvpath)
if not rows:
out["status"] = "no-gt-samples"
return out
gt = {ch: [r[i] for r in rows] for i, ch in enumerate(GEO)}
out["gt_len"] = len(rows)
out["geo_range"] = hdr.get("GeoRange")
exact = total = 0
lens, chan_status = {}, {}
ppv_err = {}
for ch in GEO:
arr = res.samples.get(ch, [])
ref = gt[ch]
lens[ch] = len(arr)
if len(arr) != len(ref):
chan_status[ch] = "length"
continue
if not arr:
chan_status[ch] = "empty"
continue
hits = sum(1 for c, v in zip(arr, ref) if abs(c * lsb - v) < 5e-5)
exact += hits
total += len(arr)
chan_status[ch] = "exact" if hits == len(arr) else "value"
gp = hdr_float(hdr, f"{ch}PPV")
if gp:
ppv_err[ch] = (max(abs(c) for c in arr) * lsb - gp) / gp
out["lens"] = lens
out["chan_status"] = chan_status
out["exact"] = exact
out["total"] = total
out["ppv_err"] = ppv_err
if all(v == "exact" for v in chan_status.values()):
out["status"] = "exact"
elif any(v == "length" for v in chan_status.values()):
out["status"] = "length-mismatch"
else:
out["status"] = "value-mismatch"
return out
def verify_histogram(binpath, csvpath, lsb):
out = {"file": os.path.basename(binpath), "status": "ok"}
try:
res = M.read_idf_file(binpath)
except NotImplementedError:
out["status"] = "not-thor"
return out
except Exception as exc: # noqa: BLE001
out["status"] = "decode-error"
out["detail"] = f"{type(exc).__name__}: {exc}"
return out
hdr, _rows = parse_export(csvpath)
out["n_intervals"] = len(res.intervals or [])
errs = {}
for ch, attr in (("Tran", "transverse_ips"), ("Vert", "vertical_ips"),
("Long", "longitudinal_ips")):
gp = hdr_float(hdr, f"{ch}PPV")
dv = getattr(res.event.peaks, attr, None)
if gp and dv:
errs[ch] = (dv - gp) / gp
out["ppv_err"] = errs
out["status"] = "peaks" if errs else "no-gt-peaks"
return out
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--root", default=DEFAULT_ROOT)
ap.add_argument("--lsb", type=float, default=M._GEO_LSB_IPS)
ap.add_argument("--limit", type=int, default=0)
ap.add_argument("--kind", choices=("idfw", "idfh", "both"), default="both")
ap.add_argument("--show", type=int, default=15, help="worst-N detail rows")
args = ap.parse_args()
pairs = index_corpus(args.root)
keys = sorted(pairs)
if args.kind != "both":
keys = [k for k in keys if k.endswith(args.kind.upper())]
if args.limit:
keys = keys[: args.limit]
print(f"root: {args.root}")
print(f"geo LSB under test: {args.lsb!r} in/s per count")
print(f"paired files: {len(keys)}\n")
wf, hg = [], []
for k in keys:
binpath, csvpath = pairs[k]
if k.endswith(".IDFW"):
wf.append(verify_waveform(binpath, csvpath, args.lsb))
else:
hg.append(verify_histogram(binpath, csvpath, args.lsb))
if wf:
st = Counter(r["status"] for r in wf)
ex = sum(r.get("exact", 0) for r in wf)
tot = sum(r.get("total", 0) for r in wf)
print("=" * 68)
print(f"WAVEFORM (IDFW): {len(wf)} files")
for s, n in st.most_common():
print(f" {s:16} {n:5d} ({100*n/len(wf):5.1f}%)")
if tot:
print(f" per-sample exact: {ex}/{tot} = {100*ex/tot:.3f}%")
errs = [e for r in wf for e in r.get("ppv_err", {}).values()]
if errs:
print(f" PPV rel-error: median {statistics.median(errs):+.4%} "
f"mean {statistics.mean(errs):+.4%} "
f"max|.| {max(abs(e) for e in errs):.4%}")
bad = [r for r in wf if r["status"] not in ("exact",)]
if bad:
print(f"\n worst {min(args.show, len(bad))} of {len(bad)} non-exact:")
for r in bad[: args.show]:
print(f" {r['file']:42} {r['status']:16} "
f"lens={r.get('lens')} gt={r.get('gt_len')} "
f"{r.get('detail','')}")
if hg:
st = Counter(r["status"] for r in hg)
print("=" * 68)
print(f"HISTOGRAM (IDFH): {len(hg)} files")
for s, n in st.most_common():
print(f" {s:16} {n:5d} ({100*n/len(hg):5.1f}%)")
errs = [e for r in hg for e in r.get("ppv_err", {}).values()]
if errs:
print(f" PPV rel-error: median {statistics.median(errs):+.4%} "
f"mean {statistics.mean(errs):+.4%} "
f"max|.| {max(abs(e) for e in errs):.4%}")
within = lambda t: 100*sum(1 for e in errs if abs(e) <= t)/len(errs)
print(f" within 0.5%: {within(0.005):.1f}% "
f"within 2%: {within(0.02):.1f}% within 5%: {within(0.05):.1f}%")
return 0
if __name__ == "__main__":
raise SystemExit(main())
+6 -17
View File
@@ -1,12 +1,12 @@
#!/usr/bin/env python3
"""Backfill events.shape_* and shape_offset_* from each event's .h5 samples. Idempotent."""
"""Backfill events.shape_* from each event's .h5 waveform samples. Idempotent."""
from __future__ import annotations
import argparse, logging, sys
from pathlib import Path
sys.path.insert(0, str(Path(__file__).resolve().parent.parent))
from sfm.database import SeismoDb
from sfm.waveform_store import WaveformStore
from sfm.shape_metrics import shape_from_h5, offset_from_h5
from sfm.shape_metrics import shape_from_h5
log = logging.getLogger("backfill_event_shape")
@@ -21,7 +21,6 @@ def backfill_shape(db: SeismoDb, store: WaveformStore, *, dry_run: bool = False)
if not h5_path.exists():
counts["skipped_no_h5"] += 1; continue
shape = shape_from_h5(h5_path)
offset = offset_from_h5(h5_path)
if shape is None:
# The .h5 can no longer yield a shape (fewer than 2 samples, or a
# flat trace). Clear any previously stored value rather than
@@ -29,32 +28,22 @@ def backfill_shape(db: SeismoDb, store: WaveformStore, *, dry_run: bool = False)
# from and silently feeds the false-trigger detector. Seen after
# a decoder fix shrinks an event: 493 rows in the prod snapshot
# were carrying metrics from a superseded decode (2026-08-25).
if (row.get("shape_crest_factor") is not None
or row.get("shape_offset") is not None):
if row.get("shape_crest_factor") is not None:
if not dry_run:
with db._connect() as conn:
conn.execute(
"UPDATE events SET shape_crest_factor=NULL, "
"shape_near_peak_count=NULL, shape_sample_count=NULL, "
"shape_axis=NULL, shape_offset=NULL, shape_offset_axis=NULL, "
"shape_offset_pre=NULL, shape_offset_spread=NULL WHERE id=?",
(row["id"],))
"shape_axis=NULL WHERE id=?", (row["id"],))
counts["cleared_stale"] += 1
counts["skipped_no_samples"] += 1; continue
if not dry_run:
with db._connect() as conn:
conn.execute(
"UPDATE events SET shape_crest_factor=?, shape_near_peak_count=?, "
"shape_sample_count=?, shape_axis=?, shape_offset=?, "
"shape_offset_axis=?, shape_offset_pre=?, shape_offset_spread=? "
"WHERE id=?",
"shape_sample_count=?, shape_axis=? WHERE id=?",
(shape["crest_factor"], shape["near_peak_count"],
shape["sample_count"], shape["axis"],
(1 if offset["offset"] else 0) if offset else None,
offset["axis"] if offset else None,
offset["pre"] if offset else None,
offset["spread"] if offset else None,
row["id"]))
shape["sample_count"], shape["axis"], row["id"]))
counts["updated"] += 1
log.info("backfill_shape: %s", counts)
return counts
+10 -47
View File
@@ -82,7 +82,6 @@ CREATE TABLE IF NOT EXISTS events (
record_type TEXT, -- "single_shot" | "continuous"
false_trigger INTEGER NOT NULL DEFAULT 0, -- 0=no, 1=yes (manual flag)
reviewed_real INTEGER NOT NULL DEFAULT 0, -- 0=no, 1=operator-confirmed real (mutually exclusive with false_trigger)
false_trigger_reason TEXT, -- optional FT cause ("offset", ...); NULL = none. Only meaningful when false_trigger=1.
blastware_filename TEXT, -- event file within waveform store; extension is per-event (AB0T encodes timestamp)
blastware_filesize INTEGER, -- bytes; NULL if no event file saved
a5_pickle_filename TEXT, -- "<filename>.a5.pkl" sidecar
@@ -100,10 +99,6 @@ CREATE TABLE IF NOT EXISTS events (
shape_near_peak_count INTEGER, -- samples >= 0.5 * peak (FT: few; real: many)
shape_sample_count INTEGER, -- total samples (to normalize near_peak_count)
shape_axis TEXT, -- geophone channel measured ("Tran"/"Vert"/"Long")
shape_offset INTEGER, -- 1 = DC-offset false trigger (pre-trigger baseline off zero + flat). Meaningful for waveforms only.
shape_offset_axis TEXT, -- geo channel the offset was measured on
shape_offset_pre REAL, -- pre-trigger baseline median (in/s)
shape_offset_spread REAL, -- max(pre,mid,end) - min(...) in in/s; small = constant/DC
created_at TEXT NOT NULL DEFAULT (strftime('%Y-%m-%dT%H:%M:%SZ', 'now')),
UNIQUE(serial, timestamp)
);
@@ -230,12 +225,7 @@ class SeismoDb:
("shape_near_peak_count", "INTEGER"),
("shape_sample_count", "INTEGER"),
("shape_axis", "TEXT"),
("shape_offset", "INTEGER"),
("shape_offset_axis", "TEXT"),
("shape_offset_pre", "REAL"),
("shape_offset_spread", "REAL"),
("reviewed_real", "INTEGER NOT NULL DEFAULT 0"),
("false_trigger_reason", "TEXT"),
):
if col not in existing_cols:
log.info("_migrate: events ADD COLUMN %s %s", col, ddl)
@@ -440,11 +430,9 @@ class SeismoDb:
tran_zc_above_range, vert_zc_above_range,
long_zc_above_range, mic_zc_above_range,
shape_crest_factor, shape_near_peak_count,
shape_sample_count, shape_axis,
shape_offset, shape_offset_axis,
shape_offset_pre, shape_offset_spread)
shape_sample_count, shape_axis)
VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?,
?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?)
?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?)
""",
(
self._new_id(), serial, key, session_id, ts,
@@ -476,10 +464,6 @@ class SeismoDb:
rec.get("shape_near_peak_count"),
rec.get("shape_sample_count"),
rec.get("shape_axis"),
rec.get("shape_offset"),
rec.get("shape_offset_axis"),
rec.get("shape_offset_pre"),
rec.get("shape_offset_spread"),
),
)
inserted += 1
@@ -533,11 +517,7 @@ class SeismoDb:
shape_crest_factor = COALESCE(?, shape_crest_factor),
shape_near_peak_count = COALESCE(?, shape_near_peak_count),
shape_sample_count = COALESCE(?, shape_sample_count),
shape_axis = COALESCE(?, shape_axis),
shape_offset = COALESCE(?, shape_offset),
shape_offset_axis = COALESCE(?, shape_offset_axis),
shape_offset_pre = COALESCE(?, shape_offset_pre),
shape_offset_spread = COALESCE(?, shape_offset_spread)
shape_axis = COALESCE(?, shape_axis)
WHERE serial = ? AND timestamp = ?
""",
(
@@ -569,10 +549,6 @@ class SeismoDb:
rec.get("shape_near_peak_count") if rec else None,
rec.get("shape_sample_count") if rec else None,
rec.get("shape_axis") if rec else None,
rec.get("shape_offset") if rec else None,
rec.get("shape_offset_axis") if rec else None,
rec.get("shape_offset_pre") if rec else None,
rec.get("shape_offset_spread") if rec else None,
serial,
ts,
),
@@ -715,9 +691,9 @@ class SeismoDb:
def propagate_review_to_twins(self, event_id: str, *, window_seconds: int | None = None) -> list[str]:
"""
Copy this event's `false_trigger`/`reviewed_real`/`false_trigger_reason`
columns onto each of its histogram/waveform twins (see `find_twins`), so
flagging one twin flags both. Returns the list of twin ids updated.
Copy this event's `false_trigger`/`reviewed_real` columns onto each
of its histogram/waveform twins (see `find_twins`), so flagging one
twin flags both. Returns the list of twin ids updated.
``window_seconds`` is accepted for backward compatibility but ignored;
twin matching is now interval-based (see `find_twins`).
@@ -727,16 +703,12 @@ class SeismoDb:
return []
ft = 1 if row.get("false_trigger") else 0
real = 1 if row.get("reviewed_real") else 0
# The reason is a subtype of the FT flag — carry it only when the source
# is actually a false trigger, so a confirmed-real twin never keeps one.
reason = row.get("false_trigger_reason") if ft else None
twins = self.find_twins(event_id)
moved = []
with self._connect() as conn:
for tw in twins:
conn.execute(
"UPDATE events SET false_trigger=?, reviewed_real=?, false_trigger_reason=? WHERE id=?",
(ft, real, reason, tw["id"]))
conn.execute("UPDATE events SET false_trigger=?, reviewed_real=? WHERE id=?",
(ft, real, tw["id"]))
moved.append(tw["id"])
return moved
@@ -757,7 +729,7 @@ class SeismoDb:
)
else:
cur = conn.execute(
"UPDATE events SET false_trigger=0, false_trigger_reason=NULL WHERE id=?",
"UPDATE events SET false_trigger=0 WHERE id=?",
(event_id,),
)
return cur.rowcount > 0
@@ -851,8 +823,7 @@ class SeismoDb:
return False
has_ft = "false_trigger" in review
has_real = "reviewed_real" in review
has_reason = "false_trigger_reason" in review
if not has_ft and not has_real and not has_reason:
if not has_ft and not has_real:
# Nothing derived to update; just confirm the row exists.
with self._connect() as conn:
row = conn.execute(
@@ -865,19 +836,11 @@ class SeismoDb:
sets["false_trigger"] = 1 if review.get("false_trigger") else 0
if has_real:
sets["reviewed_real"] = 1 if review.get("reviewed_real") else 0
if has_reason:
reason = review.get("false_trigger_reason") or None
sets["false_trigger_reason"] = reason
if reason: # a reason is a subtype of FT → implies FT
sets["false_trigger"] = 1
# mutual exclusivity: a true in one forces the other column to 0
if sets.get("false_trigger") == 1:
sets["reviewed_real"] = 0
if sets.get("reviewed_real") == 1:
sets["false_trigger"] = 0
# the reason is only meaningful while flagged FT — clear it if FT ends up 0
if sets.get("false_trigger") == 0:
sets["false_trigger_reason"] = None
assign = ", ".join(f"{k}=?" for k in sets)
params = list(sets.values()) + [event_id]
with self._connect() as conn:
+10 -19
View File
@@ -777,19 +777,6 @@ def _draw_waveform_subplot(fig, gridspec_cell, rd: ReportData) -> None:
dt_s = (rd.dt_ms or (1000.0 / sr)) / 1000.0
t0_s = (rd.t0_ms if rd.t0_ms is not None else 0.0) / 1000.0
# Shared geo scale across Long/Vert/Tran (matches the event modal + BW's
# single amp/div): all three geo lanes use ONE Y scale = the max |sample|
# across them (padded, floored), so relative amplitudes stay honest instead
# of each lane auto-zooming to its own peak. Mic keeps its own (psi) scale.
GEO_FLOOR_INS = 0.05
_geo_amax = 0.0
for _gch in ("Long", "Vert", "Tran"):
for _x in (rd.channels.get(_gch) or []):
_a = abs(_x)
if _a > _geo_amax:
_geo_amax = _a
geo_shared = max(_geo_amax * 1.10, GEO_FLOOR_INS)
last_idx = len(order) - 1
for i, ch in enumerate(order):
ax = fig.add_subplot(inner[i])
@@ -799,10 +786,10 @@ def _draw_waveform_subplot(fig, gridspec_cell, rd: ReportData) -> None:
if values:
color = _channel_axis_color(ch)
ax.plot(times, values, color=color, linewidth=0.5)
# Geo: one shared symmetric scale (honest relative amplitudes).
# Mic: symmetric on its own psi scale (different unit).
# Symmetric y-axis for geo; zero-anchored for mic.
if ch != "MicL":
ax.set_ylim(-geo_shared, geo_shared)
amax = max((abs(v) for v in values), default=0.001)
ax.set_ylim(-amax * 1.10, amax * 1.10)
else:
amax = max((abs(v) for v in values), default=0.001)
ax.set_ylim(-amax * 1.10, amax * 1.10)
@@ -837,9 +824,13 @@ def _draw_waveform_subplot(fig, gridspec_cell, rd: ReportData) -> None:
# and find peak geo amplitude for the geo amp/div setting.
total_s = times[-1] - times[0] if values else 0
div_s = total_s / 10 if total_s > 0 else 0
# Footer div value reflects the SHARED geo scale (so it's correct for all
# three lanes, not just whichever one happened to be checked first).
geo_amp_div = f"{(geo_shared * 2) / 10:.3f}" if _geo_amax > 0 else "—"
geo_amp_div = "—"
for ch in ("Tran", "Vert", "Long"):
v = rd.channels.get(ch) or []
if v:
amax = max(abs(x) for x in v)
geo_amp_div = f"{(amax * 1.1 * 2) / 10:.3f}"
break
fig.text(
0.11, 0.030,
f"Time(Seconds) {div_s:.2f} sec/div Amplitude Geo: {geo_amp_div} in/s/div Mic: 0.001 psi(L)/div",
-69
View File
@@ -47,60 +47,6 @@ def shape_from_samples(chans: dict) -> dict | None:
return s
# ── Offset (DC-baseline) detection ────────────────────────────────────────────
# A DC offset is a false trigger where the geophone baseline sits at a constant
# non-zero floor (sensor bumped / settled / drifted) instead of oscillating
# around zero. Brian's method (validated in scratch/offset_scan3.py): the
# pre-trigger window is definitionally quiet, so a true offset shows |pre| off
# zero AND stays flat across the record (pre ≈ mid ≈ end). A transient moves one
# third relative to the others and is rejected by the spread test.
# Thresholds are in in/s (the .h5 samples are already range-scaled); validated at
# Normal range (10 in/s) — the only range in the fleet.
OFFSET_FLOOR = 0.025 # |pre| at/above this reads as an off-zero baseline (5 A/D counts)
OFFSET_MAX_SPREAD = 0.02 # max(pre,mid,end) - min(...) at/below this reads as flat/constant
def _channel_offset(x, pretrig_n):
"""Return (pre, spread, is_offset) for one channel, or None if unusable."""
x = np.asarray(x, dtype=float)
n = x.size
if n < 3:
return None
t = n // 3
pre = x[:pretrig_n] if (pretrig_n and 0 < pretrig_n < n) else x[:t]
mid, end = x[t:2 * t], x[2 * t:]
if pre.size == 0 or mid.size == 0 or end.size == 0:
return None
vals = [float(np.median(seg)) for seg in (pre, mid, end)]
spread = max(vals) - min(vals)
is_offset = abs(vals[0]) >= OFFSET_FLOOR and spread <= OFFSET_MAX_SPREAD
return vals[0], spread, is_offset
def offset_from_samples(chans: dict, pretrig_n) -> dict | None:
"""Detect a DC-offset false trigger across the geophone channels.
An event is offset if ANY geo channel's pre-trigger baseline is off zero and
flat across the record. Reports the tripping axis (or, if none trips, the
most-offset-like axis) with its ``pre``/``spread`` for transparency + tuning.
Returns None when no geo channel is usable.
"""
results = []
for ax in _GEO_CHANNELS:
x = chans.get(ax)
if x is None:
continue
r = _channel_offset(x, pretrig_n)
if r is not None:
results.append((ax, r[0], r[1], r[2]))
if not results:
return None
offenders = [r for r in results if r[3]]
ax, pre, spread, _ = max(offenders or results, key=lambda r: abs(r[1]))
return {"offset": bool(offenders), "axis": ax,
"pre": round(pre, 6), "spread": round(spread, 6)}
def shape_from_h5(path) -> dict | None:
import h5py
try:
@@ -110,18 +56,3 @@ def shape_from_h5(path) -> dict | None:
except Exception:
return None
return shape_from_samples(chans)
def offset_from_h5(path) -> dict | None:
"""offset_from_samples fed from an event's .h5 (float32 in/s geo samples +
the pretrig_samples attribute)."""
import h5py
try:
with h5py.File(path, "r") as f:
chans = {ax: f[f"samples/{ax}"][:] for ax in _GEO_CHANNELS
if f"samples/{ax}" in f}
pretrig_n = f.attrs.get("pretrig_samples")
except Exception:
return None
pretrig_n = int(pretrig_n) if pretrig_n is not None else 0
return offset_from_samples(chans, pretrig_n)
+13 -99
View File
@@ -32,7 +32,6 @@ from __future__ import annotations
import datetime
import logging
import pickle
import re
import shutil
from pathlib import Path
from typing import Optional, Union
@@ -42,7 +41,7 @@ from minimateplus.blastware_file import blastware_filename, write_blastware_file
from minimateplus.framing import S3Frame
from minimateplus.models import Event
from sfm import event_hdf5
from sfm.shape_metrics import shape_from_h5, offset_from_h5
from sfm.shape_metrics import shape_from_h5
log = logging.getLogger("sfm.waveform_store")
@@ -271,13 +270,6 @@ class WaveformStore:
"shape_sample_count": _shape["sample_count"],
"shape_axis": _shape["axis"],
} if _shape else {}
_offset = offset_from_h5(hdf5_path) if hdf5_filename else None
_offset_rec = {
"shape_offset": 1 if _offset["offset"] else 0,
"shape_offset_axis": _offset["axis"],
"shape_offset_pre": _offset["pre"],
"shape_offset_spread": _offset["spread"],
} if _offset else {}
return {
"filename": filename,
"filesize": filesize,
@@ -286,7 +278,6 @@ class WaveformStore:
"hdf5_filename": hdf5_filename,
"sidecar_filename": sidecar_path.name,
**_shape_rec,
**_offset_rec,
}
def save_imported_bw(
@@ -380,16 +371,8 @@ class WaveformStore:
# Resolve serial. blastware_filename derives a 4-char prefix from
# the numeric serial (e.g. BE11529 → M529); we go the other way
# if a hint wasn't given. The filename carries only the NUMBER,
# so read the family prefix out of the body first — a BlastMate
# ("BA") filed as "BE" is a unit that does not exist. The
# filename-only decoder stays as the last resort.
serial = (
serial_hint
or _serial_from_bw_bytes(bw_bytes, source_path.name)
or _serial_from_bw_filename(source_path.name)
or "UNKNOWN"
)
# via the source filename if a hint wasn't given.
serial = serial_hint or _serial_from_bw_filename(source_path.name) or "UNKNOWN"
# Use the source filename verbatim — it already encodes timestamp
# + record type per BW's AB0T scheme, and we want to preserve it
@@ -478,13 +461,6 @@ class WaveformStore:
"shape_sample_count": _shape["sample_count"],
"shape_axis": _shape["axis"],
} if _shape else {}
_offset = offset_from_h5(hdf5_path) if hdf5_filename else None
_offset_rec = {
"shape_offset": 1 if _offset["offset"] else 0,
"shape_offset_axis": _offset["axis"],
"shape_offset_pre": _offset["pre"],
"shape_offset_spread": _offset["spread"],
} if _offset else {}
return ev, {
"filename": filename,
"filesize": filesize,
@@ -494,7 +470,6 @@ class WaveformStore:
"sidecar_filename": sidecar_path.name,
"serial": serial,
**_shape_rec,
**_offset_rec,
}
def save_imported_idf(
@@ -595,19 +570,8 @@ class WaveformStore:
)
# Binary-derived peaks fill in when the .txt didn't supply them.
#
# The old justification for this precedence -- "binary peaks are ~3%
# low vs the .txt" -- was a decoder bug (geo LSB 0.0003 instead of
# 0.000310308) and was fixed 2026-09-10; the binary now agrees with
# Thor's own export per-sample. The .txt still wins when present
# because it is what the operator sees in Thor's report.
#
# ⚠ One case where the .txt is the *less* accurate of the two:
# Thor floors displayed histogram PPV at 0.0050 in/s, so on quiet
# IDFH events the .txt reports 0.0050 while the binary decodes the
# true ~0.0025. 41.4% of prod IDFH sidecars carry a component PPV
# larger than their own vector sum because of it. Left as-is
# deliberately, so stored peaks keep matching Thor's report.
# They're ~3% low vs the device-authoritative .txt values (residual
# codec drift), so .txt always wins when present.
if binary_peaks is not None:
if binary_peaks.transverse_ips and not report_dict.get("tran_ppv"):
report_dict["tran_ppv"] = binary_peaks.transverse_ips
@@ -787,13 +751,6 @@ class WaveformStore:
"shape_sample_count": _shape["sample_count"],
"shape_axis": _shape["axis"],
} if _shape else {}
_offset = offset_from_h5(hdf5_path) if hdf5_filename else None
_offset_rec = {
"shape_offset": 1 if _offset["offset"] else 0,
"shape_offset_axis": _offset["axis"],
"shape_offset_pre": _offset["pre"],
"shape_offset_spread": _offset["spread"],
} if _offset else {}
return ev, {
"filename": filename,
"filesize": filesize,
@@ -803,7 +760,6 @@ class WaveformStore:
"sidecar_filename": sidecar_path.name,
"serial": serial,
**_shape_rec,
**_offset_rec,
}
def load_a5(self, serial: str, filename: str) -> Optional[list[S3Frame]]:
@@ -860,24 +816,20 @@ class WaveformStore:
# ── helpers ─────────────────────────────────────────────────────────────────────
def _serial_number_from_bw_filename(name: str) -> Optional[int]:
def _serial_from_bw_filename(name: str) -> Optional[str]:
"""
Reverse of `blastware_filename`'s serial-prefix encoding — the NUMBER only.
Reverse of `blastware_filename`'s serial-prefix encoding.
BW filename format (V10.72): `<P><serial3><stem4>.<ext>`
where P = chr(ord('B') + floor(serial // 1000))
and serial3 = f"{serial % 1000:03d}".
Examples (from CLAUDE.md verification archive):
P036... → 14036 H907... → 6907
M529... → 11529 T003... → 18003
L895... → 10895
P036... → BE14036 H907... → BE6907
M529... → BE11529 T003... → BE18003
⚠ The filename encodes **only the number**. The two-letter family
prefix is NOT in it — "BE" is a MiniMate Plus, "BA" a BlastMate — so
the prefix has to come from the file body (`_serial_from_bw_bytes`)
or from an explicit hint. Returns None when the filename doesn't
match the expected pattern.
Returns the inferred BE-prefix serial (e.g. "BE11529") or None when
the filename doesn't match the expected pattern.
"""
if not name:
return None
@@ -890,43 +842,5 @@ def _serial_number_from_bw_filename(name: str) -> Optional[int]:
if prefix_letter < "B":
return None
thousands = ord(prefix_letter) - ord("B")
return thousands * 1000 + int(base[1:4])
_BW_SERIAL_RE = re.compile(rb"[A-Z]{2}\d{3,6}")
def _serial_from_bw_bytes(data: bytes, name: str) -> Optional[str]:
"""
Read the real serial — prefix included — out of a BW file body.
The body carries the serial as a plain ASCII string ("BE9558",
"BA10895"). We accept a candidate only when its numeric part matches
the number the filename encodes, which keeps a stray byte sequence in
the sample stream from being mistaken for a serial.
Returns None when the filename number can't be derived or no
candidate in the body agrees with it — the caller then falls back.
"""
num = _serial_number_from_bw_filename(name)
if num is None or not data:
return None
for match in _BW_SERIAL_RE.findall(data):
candidate = match.decode("ascii", errors="replace")
if candidate[2:].lstrip("0") == str(num):
return candidate
return None
def _serial_from_bw_filename(name: str) -> Optional[str]:
"""
Best-effort serial from the filename alone.
⚠ The family prefix is a **guess** — the filename does not carry it.
"BE" is right for every MiniMate Plus but wrong for a BlastMate, whose
serials start "BA". Prefer `_serial_from_bw_bytes` whenever the file
body is at hand; this exists for callers that only have a name
(log lines, dry-run output).
"""
num = _serial_number_from_bw_filename(name)
return None if num is None else f"BE{num}"
serial_num = thousands * 1000 + int(base[1:4])
return f"BE{serial_num}"
-103
View File
@@ -1,103 +0,0 @@
import sqlite3
from sfm.database import SeismoDb
from minimateplus.models import Event, Timestamp
def _ev(db, key="0111aaaa", serial="BE1"):
ev = Event(index=0)
ev._waveform_key = bytes.fromhex(key)
ev.timestamp = Timestamp(raw=b"", flag=0x10, year=2026, unknown_byte=0,
month=6, day=25, hour=8, minute=0, second=0)
ev.record_type = "Waveform"
db.insert_events([ev], serial=serial)
return [r for r in db.query_events(serial=serial) if r["waveform_key"] == key][0]["id"]
def test_flag_offset_reason_implies_ft(tmp_path):
db = SeismoDb(tmp_path / "s.db")
eid = _ev(db)
db.update_event_review(eid, {"false_trigger_reason": "offset"})
row = db.get_event(eid)
assert row["false_trigger"] == 1 # a reason is a subtype of FT
assert row["false_trigger_reason"] == "offset"
assert row["reviewed_real"] == 0
def test_plain_ft_leaves_reason_null(tmp_path):
# Reason is OPTIONAL — flagging FT without one records no reason.
db = SeismoDb(tmp_path / "s.db")
eid = _ev(db)
db.update_event_review(eid, {"false_trigger": True})
row = db.get_event(eid)
assert row["false_trigger"] == 1
assert row["false_trigger_reason"] is None
def test_confirm_real_clears_reason(tmp_path):
db = SeismoDb(tmp_path / "s.db")
eid = _ev(db)
db.update_event_review(eid, {"false_trigger_reason": "offset"})
db.update_event_review(eid, {"reviewed_real": True})
row = db.get_event(eid)
assert row["reviewed_real"] == 1
assert row["false_trigger"] == 0
assert row["false_trigger_reason"] is None
def test_clear_ft_clears_reason(tmp_path):
db = SeismoDb(tmp_path / "s.db")
eid = _ev(db)
db.update_event_review(eid, {"false_trigger_reason": "offset"})
db.update_event_review(eid, {"false_trigger": False})
row = db.get_event(eid)
assert row["false_trigger"] == 0
assert row["false_trigger_reason"] is None
def test_set_false_trigger_false_clears_reason(tmp_path):
db = SeismoDb(tmp_path / "s.db")
eid = _ev(db)
db.update_event_review(eid, {"false_trigger_reason": "offset"})
assert db.set_false_trigger(eid, False) is True
row = db.get_event(eid)
assert row["false_trigger"] == 0
assert row["false_trigger_reason"] is None
def test_reason_can_be_cleared_without_clearing_ft(tmp_path):
# Setting reason to None removes the reason but leaves the FT flag intact.
db = SeismoDb(tmp_path / "s.db")
eid = _ev(db)
db.update_event_review(eid, {"false_trigger_reason": "offset"})
db.update_event_review(eid, {"false_trigger_reason": None})
row = db.get_event(eid)
assert row["false_trigger"] == 1
assert row["false_trigger_reason"] is None
def _ts(h, m, d=25):
return Timestamp(raw=b"", flag=0x10, year=2026, unknown_byte=0,
month=2, day=d, hour=h, minute=m, second=0)
def test_offset_reason_propagates_to_twin(tmp_path):
# Flag a waveform as offset → its histogram twin also becomes FT with reason=offset.
db = SeismoDb(tmp_path / "s.db")
def ins(key, ts, rt):
ev = Event(index=0); ev._waveform_key = bytes.fromhex(key); ev.timestamp = ts
db.insert_events([ev], serial="BE1")
rid = [r for r in db.query_events(serial="BE1") if r["waveform_key"] == key][0]["id"]
with sqlite3.connect(db.db_path) as c:
c.execute("UPDATE events SET peak_vector_sum=0.4763, record_type=? WHERE id=?", (rt, rid))
return rid
hist = ins("01110001", _ts(19, 31), "Histogram") # interval start
wave = ins("01110002", _ts(20, 46), "Waveform") # trigger inside the interval
db.update_event_review(wave, {"false_trigger_reason": "offset"})
db.propagate_review_to_twins(wave)
row = db.get_event(hist)
assert row["false_trigger"] == 1
assert row["false_trigger_reason"] == "offset"
-322
View File
@@ -1,322 +0,0 @@
"""Per-sample verification of the Thor / Micromate (series-4) IDF binary codec.
Ground truth is Thor's own CSV export, written next to each binary by the
Thor desktop application. For waveforms the export carries a per-sample
block of four columns (Tran, Vert, Long, Mic) in in/s and psi -- the
series-4 equivalent of Blastware's ``_ASCII.TXT`` exports.
The full-corpus harness is ``scratch/verify_thor_against_csv.py``; these
tests pin the two constants that harness established so they cannot
regress silently.
"""
from __future__ import annotations
import csv
import os
import sys
from pathlib import Path
import pytest
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
from micromate.idf_file import (
_GEO_LSB_IPS,
geo_count_to_ips,
read_idf_file,
)
FIXTURES = Path(__file__).parent / "fixtures" / "thor-idf"
IDFW = FIXTURES / "UM11719_20231219162723.IDFW"
IDFH = FIXTURES / "UM11719_20231219162648.IDFH"
GEO_CHANNELS = ("Tran", "Vert", "Long")
# tests/fixtures/ is gitignored, so a fresh checkout has no sample data.
# Skip rather than fail, matching test_idf_ascii_report.py. To populate:
#
# B="<thor-watcher>/example-data/THORDATA_example/THORDATA_example/UPMC Presby"
# mkdir -p tests/fixtures/thor-idf
# for f in UM11719/UM11719_20231219162723.IDFW \
# UM11719/UM11719_20231219162648.IDFH \
# UM13981/UM13981_20220207084555.IDFW \
# UM13981/UM13981_20220207183102.IDFH \
# UM13981/UM13981_20221202063059.IDFH; do
# cp "$B/$f" tests/fixtures/thor-idf/
# cp "$B/$(dirname $f)/CSV/$(basename $f).csv" tests/fixtures/thor-idf/
# done
pytestmark = pytest.mark.skipif(
not FIXTURES.is_dir() or not any(FIXTURES.glob("*.IDFW")),
reason=f"Thor IDF fixtures not present under {FIXTURES}",
)
def _parse_export(path: Path):
"""Split a Thor CSV export into (header dict, per-sample rows)."""
header, rows = {}, []
with path.open(newline="", encoding="utf-8", errors="replace") as fh:
for rec in csv.reader(fh):
if len(rec) == 2:
header[rec[0].strip()] = rec[1].strip()
elif len(rec) >= 3:
try:
rows.append([float(x) for x in rec])
except ValueError:
pass
return header, rows
def _header_float(header, key):
return float(header[key].split()[0])
@pytest.fixture(scope="module")
def idfw_export():
return _parse_export(IDFW.with_suffix(".IDFW.csv"))
# ─── The geo scale constant ────────────────────────────────────────────────
def test_geo_lsb_matches_thor_quantisation():
"""Thor's own export quantises geo samples to this LSB.
Derived by maximising exact-match count over 1,046,016 paired samples
(454 channel-events, 2 units); independently corroborated on 8
production units via their device-reported PPV. The historical value
0.0003 read every series-4 geophone sample 3.3% low.
"""
assert _GEO_LSB_IPS == pytest.approx(0.000310308, rel=1e-6)
def test_geo_lsb_is_not_the_legacy_value():
# Guards against a revert to the truncated 0.0003 constant.
assert abs(_GEO_LSB_IPS - 0.0003) > 1e-6
# ─── Per-sample fidelity ───────────────────────────────────────────────────
def test_waveform_channel_lengths_match_export(idfw_export):
_header, rows = idfw_export
result = read_idf_file(IDFW)
for channel in GEO_CHANNELS:
assert len(result.samples[channel]) == len(rows), (
f"{channel} truncated: decoded {len(result.samples[channel])} "
f"samples, export has {len(rows)}"
)
def test_waveform_samples_match_export_exactly(idfw_export):
"""Every geo sample must reproduce Thor's exported value to 4 dp."""
_header, rows = idfw_export
result = read_idf_file(IDFW)
for index, channel in enumerate(GEO_CHANNELS):
decoded = result.samples[channel]
expected = [row[index] for row in rows]
mismatches = [
(i, geo_count_to_ips(c), v)
for i, (c, v) in enumerate(zip(decoded, expected))
if abs(geo_count_to_ips(c) - v) >= 5e-5
]
assert not mismatches, (
f"{channel}: {len(mismatches)} of {len(expected)} samples differ; "
f"first three {mismatches[:3]}"
)
def test_waveform_ppv_matches_export(idfw_export):
header, _rows = idfw_export
result = read_idf_file(IDFW)
for channel, attr in (
("Tran", "transverse_ips"),
("Vert", "vertical_ips"),
("Long", "longitudinal_ips"),
):
decoded = getattr(result.event.peaks, attr)
assert decoded == pytest.approx(
_header_float(header, f"{channel}PPV"), abs=5e-5
), f"{channel} PPV disagrees with Thor's export"
# ─── Histogram path shares the same scale ──────────────────────────────────
def test_histogram_peaks_match_export():
header, _rows = _parse_export(IDFH.with_suffix(".IDFH.csv"))
result = read_idf_file(IDFH)
assert result.intervals, "IDFH decoded no intervals"
for channel, attr in (
("Tran", "transverse_ips"),
("Vert", "vertical_ips"),
("Long", "longitudinal_ips"),
):
decoded = getattr(result.event.peaks, attr)
expected = _header_float(header, f"{channel}PPV")
# Histogram peaks are stored per-interval, so the export's PPV is
# reproduced within one quantisation step rather than exactly.
assert decoded == pytest.approx(expected, abs=2 * _GEO_LSB_IPS), (
f"{channel} histogram peak {decoded} vs export {expected}"
)
# ─── Regressions found 2026-09-10 ──────────────────────────────────────────
IDFH_LONG = FIXTURES / "UM13981_20220207183102.IDFH" # 719 intervals
IDFH_SENTINEL = FIXTURES / "UM13981_20221202063059.IDFH" # holds an unwritten slot
IDFW_RAW16 = FIXTURES / "UM13981_20220207084555.IDFW" # segment 0 is MODE_RAW16
def test_histogram_decodes_past_250_intervals():
"""The segment validator must not require a zero counter high byte.
The interval counter is a uint16 cumulative index. Requiring its high
byte to be zero rejected every segment past interval 255, capping each
histogram at 250 intervals and truncating any run longer than ~4 hours —
frequently discarding the part that held the peak.
"""
result = read_idf_file(IDFH_LONG)
header, _rows = _parse_export(IDFH_LONG.with_suffix(".IDFH.csv"))
expected = float(header["NumberOfIntervals"])
assert len(result.intervals) == 719
assert len(result.intervals) == pytest.approx(expected, abs=1.0)
def test_histogram_ignores_unwritten_interval_slot():
"""A never-written interval keeps its ±full-scale seed and must be dropped.
Counting it fabricates a 10.0 in/s peak on every channel, which then wins
the max-over-intervals and poisons the whole file's PPV.
"""
header, _rows = _parse_export(IDFH_SENTINEL.with_suffix(".IDFH.csv"))
result = read_idf_file(IDFH_SENTINEL)
for channel, attr in (
("Tran", "transverse_ips"),
("Vert", "vertical_ips"),
("Long", "longitudinal_ips"),
):
decoded = getattr(result.event.peaks, attr)
assert decoded < 1.0, f"{channel} peak {decoded} looks like the ±FS seed"
assert decoded == pytest.approx(
_header_float(header, f"{channel}PPV"), abs=2 * _GEO_LSB_IPS
)
def test_waveform_raw16_segment_zero_is_decoded():
"""Segment-0 records can be raw int16 (MODE_RAW16, 10-byte header).
That mode was absent from the dispatch, so the record fell through
unhandled and the channel silently lost its first 512 samples.
"""
rows = _parse_export(IDFW_RAW16.with_suffix(".IDFW.csv"))[1]
result = read_idf_file(IDFW_RAW16)
for index, channel in enumerate(GEO_CHANNELS):
decoded = result.samples[channel]
assert len(decoded) == len(rows), f"{channel} lost segment 0"
expected = [row[index] for row in rows]
bad = sum(
1 for c, v in zip(decoded, expected)
if abs(geo_count_to_ips(c) - v) >= 5e-5
)
assert bad == 0, f"{channel}: {bad} samples differ from Thor's export"
def test_body_offset_search_is_not_quadratic():
"""The body scan must stay cheap enough for bulk ingest.
MODE_RAW16 is (0x00, 0x00), so scanning for candidate *preambles* treats
every run of three zero bytes as a body start and trial-decodes each one
(~0.5 s/file measured). The search anchors on record headers instead.
"""
import time
start = time.perf_counter()
for _ in range(3):
read_idf_file(IDFW_RAW16)
elapsed = (time.perf_counter() - start) / 3
assert elapsed < 0.15, f"body-offset search took {elapsed*1000:.0f} ms/file"
# ─── Mic-disabled (3-channel) units, found 2026-09-10 ──────────────────────
IDFW_3CH = FIXTURES / "UM20147_20250531135901.IDFW" # body head below old floor
IDFH_3CH = FIXTURES / "UM20147_20250330070110.IDFH" # 56-byte interval records
def test_three_channel_waveform_decodes_all_geo_channels():
"""A mic-disabled unit's shorter header moves the record chain head.
Its head sits at 0x0dba, below the old ``_BODY_SCAN_FLOOR`` of 0x0E00, so
the scan could not see it and fell through to the *Vert* segment-0 record
— decoding a body shifted one position around the channel rotation, which
surfaced as Vert being exactly 512 samples short.
"""
rows = _parse_export(IDFW_3CH.with_suffix(".IDFW.csv"))[1]
result = read_idf_file(IDFW_3CH)
for index, channel in enumerate(GEO_CHANNELS):
decoded = result.samples[channel]
assert len(decoded) == len(rows), (
f"{channel}: {len(decoded)} samples, export has {len(rows)}"
)
expected = [row[index] for row in rows]
bad = sum(
1 for c, v in zip(decoded, expected)
if abs(geo_count_to_ips(c) - v) >= 5e-5
)
assert bad == 0, f"{channel}: {bad} samples differ from Thor's export"
# Mic is genuinely absent on these units, not merely undecoded.
assert not result.samples.get("MicL")
def test_three_channel_histogram_uses_56_byte_intervals():
"""Interval stride is 16 bytes per channel + an 8-byte tail, not a constant.
A mic-disabled unit packs 56-byte records, so assuming 72 read 7 intervals
out of every 10-interval segment and then walked off alignment into
garbage, which decoded as ~10 in/s peaks. The true count comes from the
segment's cumulative interval counter.
"""
header, _rows = _parse_export(IDFH_3CH.with_suffix(".IDFH.csv"))
result = read_idf_file(IDFH_3CH)
expected_intervals = float(header["NumberOfIntervals"])
assert len(result.intervals) == pytest.approx(expected_intervals, abs=1.0)
assert {iv.n_channels for iv in result.intervals} == {3}
for channel, attr in (
("Tran", "transverse_ips"),
("Vert", "vertical_ips"),
("Long", "longitudinal_ips"),
):
decoded = getattr(result.event.peaks, attr)
assert decoded < 1.0, f"{channel} peak {decoded} looks like walked-off garbage"
assert decoded == pytest.approx(
_header_float(header, f"{channel}PPV"), rel=0.02
)
# ─── `40 NN` blocks with NN > 8, verified 2026-09-11 ───────────────────────
IDFW_WIDE40 = FIXTURES / "UM12947_20250806134504.IDFW"
def test_wide_forty_nn_block_does_not_truncate_channels():
"""Loud events use `40 NN` blocks with NN well above the old cap of 8.
``data_block_len()`` rejected NN > 0x08, which halted the block walk
part-way through a record. The walker stops at the first unrecognised
tag instead of raising, so this surfaced as silently short channels —
here Tran 1812 / Vert 2132 / Long 2324 where the export has 2324 for all
three. The affected files use NN of 12, 16, 20 ... up to 196.
"""
rows = _parse_export(IDFW_WIDE40.with_suffix(".IDFW.csv"))[1]
result = read_idf_file(IDFW_WIDE40)
for index, channel in enumerate(GEO_CHANNELS):
decoded = result.samples[channel]
assert len(decoded) == len(rows), (
f"{channel}: {len(decoded)} samples, export has {len(rows)}"
)
expected = [row[index] for row in rows]
bad = sum(
1 for c, v in zip(decoded, expected)
if abs(geo_count_to_ips(c) - v) >= 5e-5
)
assert bad == 0, f"{channel}: {bad} samples differ from Thor's export"
-94
View File
@@ -1,94 +0,0 @@
import numpy as np
import h5py
from sfm.shape_metrics import offset_from_samples, offset_from_h5
def test_flags_constant_dc_floor():
# A geophone channel sitting at a constant +0.05 in/s across the whole record
# is a DC offset: baseline off zero AND flat across pre/mid/end thirds.
n = 300
chans = {"Tran": np.full(n, 0.05), "Vert": np.zeros(n), "Long": np.zeros(n)}
r = offset_from_samples(chans, pretrig_n=50)
assert r["offset"] is True
assert r["axis"] == "Tran"
assert abs(r["pre"] - 0.05) < 1e-6
assert r["spread"] < 0.02
def test_transient_rejected_by_spread():
# Off-zero pre-trigger but the baseline SETTLES back over the record — a
# transient, not a constant offset. The spread test must reject it.
x = np.concatenate([np.full(100, 0.05), np.full(100, 0.025), np.zeros(100)])
chans = {"Tran": x, "Vert": np.zeros(300), "Long": np.zeros(300)}
r = offset_from_samples(chans, pretrig_n=100)
assert r["offset"] is False
def test_clean_oscillation_not_offset():
t = np.arange(300)
x = 0.4 * np.sin(2 * np.pi * t / 20) # oscillates around zero — baseline IS zero
chans = {"Tran": x, "Vert": np.zeros(300), "Long": np.zeros(300)}
r = offset_from_samples(chans, pretrig_n=50)
assert r["offset"] is False
def test_below_floor_not_offset_but_reports_pre():
# A flat baseline below the floor is not an offset; still report the axis/pre
# for tuning transparency.
n = 300
chans = {"Tran": np.full(n, 0.01), "Vert": np.zeros(n), "Long": np.zeros(n)}
r = offset_from_samples(chans, pretrig_n=50)
assert r["offset"] is False
assert r["axis"] == "Tran"
assert abs(r["pre"] - 0.01) < 1e-6
def test_none_when_no_geo_channels():
assert offset_from_samples({"MicL": np.full(300, 0.05)}, pretrig_n=50) is None
def test_pretrig_fallback_when_invalid():
# pretrig_n of 0 (missing/unusable) falls back to the first third.
n = 300
chans = {"Tran": np.full(n, 0.05), "Vert": np.zeros(n), "Long": np.zeros(n)}
r = offset_from_samples(chans, pretrig_n=0)
assert r["offset"] is True
def test_flags_offset_on_any_axis():
# Offset on Vert alone still flags the event, and Vert is reported.
n = 300
chans = {"Tran": np.zeros(n), "Vert": np.full(n, -0.06), "Long": np.zeros(n)}
r = offset_from_samples(chans, pretrig_n=50)
assert r["offset"] is True
assert r["axis"] == "Vert"
def _write_h5(path, chans, pretrig_n):
with h5py.File(path, "w") as f:
g = f.create_group("samples")
for k, v in chans.items():
g.create_dataset(k, data=np.asarray(v, dtype="float32"))
if pretrig_n is not None:
f.attrs["pretrig_samples"] = pretrig_n
def test_offset_from_h5_reads_pretrig_attr(tmp_path):
p = tmp_path / "ev.h5"
n = 300
_write_h5(p, {"Tran": np.full(n, 0.05), "Vert": np.zeros(n), "Long": np.zeros(n)},
pretrig_n=50)
r = offset_from_h5(str(p))
assert r["offset"] is True and r["axis"] == "Tran"
def test_offset_from_h5_missing_pretrig_attr_falls_back(tmp_path):
p = tmp_path / "noattr.h5"
n = 300
_write_h5(p, {"Tran": np.full(n, 0.05), "Vert": np.zeros(n), "Long": np.zeros(n)},
pretrig_n=None)
assert offset_from_h5(str(p))["offset"] is True # falls back to first-third
def test_offset_from_h5_missing_file_is_none(tmp_path):
assert offset_from_h5(str(tmp_path / "nope.h5")) is None
-64
View File
@@ -1,64 +0,0 @@
from __future__ import annotations
from pathlib import Path
import numpy as np, h5py
from sfm.database import SeismoDb
from sfm.waveform_store import WaveformStore
from scripts.backfill_event_shape import backfill_shape
from minimateplus.models import Event, Timestamp, PeakValues
_FIX = Path(__file__).parent / "fixtures/histogram-extension-re/events-5-21-26/K558LL8B.7I0W"
def _event(waveform_key="0111abcd"):
ev = Event(index=0)
ev._waveform_key = bytes.fromhex(waveform_key)
ev.timestamp = Timestamp(raw=b"", flag=0x10, year=2026, unknown_byte=0,
month=6, day=25, hour=8, minute=50, second=0)
ev.record_type = "Waveform"
ev.peak_values = PeakValues(tran=0.075, vert=0.220, long=0.045,
peak_vector_sum=0.231, micl=0.01)
return ev
def test_insert_stores_offset_from_record(tmp_path: Path):
db = SeismoDb(tmp_path / "s.db")
ev = _event()
rec = {ev._waveform_key.hex(): {
"filename": "F.CE0W", "filesize": 10,
"shape_offset": 1, "shape_offset_axis": "Tran",
"shape_offset_pre": 0.05, "shape_offset_spread": 0.001}}
db.insert_events([ev], serial="BE1", waveform_records=rec)
row = db.query_events(serial="BE1")[0]
assert row["shape_offset"] == 1
assert row["shape_offset_axis"] == "Tran"
assert abs(row["shape_offset_pre"] - 0.05) < 1e-6
assert abs(row["shape_offset_spread"] - 0.001) < 1e-6
def test_save_imported_bw_attaches_offset(tmp_path: Path):
store = WaveformStore(tmp_path / "waveforms")
ev, rec = store.save_imported_bw(_FIX.read_bytes(), source_path=_FIX, serial_hint="BE9558")
assert rec["shape_offset"] in (0, 1)
assert rec["shape_offset_axis"] in ("Tran", "Vert", "Long")
assert "shape_offset_pre" in rec and "shape_offset_spread" in rec
def test_backfill_updates_offset(tmp_path: Path):
db = SeismoDb(tmp_path / "s.db")
store = WaveformStore(tmp_path / "waveforms")
ev = Event(index=0); ev._waveform_key = bytes.fromhex("0111abcd")
db.insert_events([ev], serial="BE1",
waveform_records={ev._waveform_key.hex(): {"filename": "F.CE0W", "filesize": 10}})
p = store.hdf5_path_for("BE1", "F.CE0W")
with h5py.File(p, "w") as f:
g = f.create_group("samples")
g.create_dataset("Tran", data=np.full(300, 0.05, "float32"))
g.create_dataset("Vert", data=np.zeros(300, "float32"))
g.create_dataset("Long", data=np.zeros(300, "float32"))
f.attrs["pretrig_samples"] = 50
backfill_shape(db, store)
row = db.query_events(serial="BE1")[0]
assert row["shape_offset"] == 1
assert row["shape_offset_axis"] == "Tran"
-61
View File
@@ -1,61 +0,0 @@
"""The event-report PDF must draw the three geo channels on ONE shared Y scale
(max |sample| across Long/Vert/Tran, floored), not each trace auto-zoomed to its
own peak — so relative amplitudes are honest and a small channel doesn't fill its
lane looking as big as a large one. Mirrors the event-modal waveform behaviour.
"""
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import pytest
from sfm.report_pdf import ReportData, _draw_waveform_subplot
def _draw(channels):
rd = ReportData(
channels=channels,
sample_rate_sps=1024,
dt_ms=1000.0 / 1024,
t0_ms=0.0,
)
fig = plt.figure()
cell = fig.add_gridspec(1, 1)[0, 0]
_draw_waveform_subplot(fig, cell, rd)
by_label = {ax.get_ylabel(): ax for ax in fig.axes}
try:
yield_ = {k: by_label[k].get_ylim() for k in ("Long", "Vert", "Tran", "MicL")}
finally:
plt.close(fig)
return yield_
def test_geo_traces_share_one_y_scale():
# Tran is the biggest geo channel (0.35); Long 0.10, Vert 0.02.
ylims = _draw({
"Long": [0.10, -0.10, 0.0],
"Vert": [0.02, -0.02, 0.0],
"Tran": [0.35, -0.35, 0.0],
"MicL": [0.0005, -0.0005, 0.0],
})
# Shared scale = max(0.35 * 1.10, floor 0.05) = 0.385, symmetric.
expected = pytest.approx(0.385, rel=1e-6)
for ch in ("Long", "Vert", "Tran"):
lo, hi = ylims[ch]
assert hi == expected, f"{ch} top ylim {hi} != shared 0.385"
assert lo == pytest.approx(-0.385, rel=1e-6), f"{ch} bottom ylim {lo}"
# All three geo lanes identical.
assert ylims["Long"] == ylims["Vert"] == ylims["Tran"]
# Mic keeps its own (much smaller) scale — not lumped into the geo max.
assert ylims["MicL"][1] < 0.01
def test_geo_shared_scale_has_floor():
# A tiny event (all geo well under the floor) clamps to the 0.05 floor.
ylims = _draw({
"Long": [0.008, -0.008, 0.0],
"Vert": [0.006, -0.006, 0.0],
"Tran": [0.010, -0.010, 0.0],
"MicL": [0.0001, -0.0001, 0.0],
})
for ch in ("Long", "Vert", "Tran"):
assert ylims[ch][1] == pytest.approx(0.05, rel=1e-6), f"{ch} not floored"
-101
View File
@@ -1,101 +0,0 @@
"""The BW filename encodes the serial NUMBER, never the family prefix.
"BE" is a MiniMate Plus; "BA" is a BlastMate. Both are Series III and their
files are byte-compatible — the whole archive's 1,493 BlastMate binaries
decode through the same codec at 100% — so the only thing that distinguishes
them downstream is the serial string, and that lives in the file body.
Synthesising the prefix as "BE" files a BlastMate under a unit that does not
exist. Four units in the DL2 archive are affected: BA9229, BA10060, BA10895
and BA15957.
"""
from __future__ import annotations
import pytest
from minimateplus.client import _decode_0a_partial_header
from sfm.waveform_store import (
_serial_from_bw_bytes,
_serial_from_bw_filename,
_serial_number_from_bw_filename,
)
# ── the filename gives a number, and only a number ──────────────────────────
@pytest.mark.parametrize("name,num", [
("P036L318.C80H", 14036), # BE14036
("H907KWRK.WB0H", 6907), # BE6907
("M529LKIQ.G10", 11529), # BE11529
("T003LQ9K.OE0H", 18003), # BE18003
("L895K63F.GE0W", 10895), # BA10895 — a BlastMate
("K229HGQI.XO0W", 9229), # BA9229 — a BlastMate
])
def test_number_from_filename(name, num):
assert _serial_number_from_bw_filename(name) == num
@pytest.mark.parametrize("name", ["", "not_a_bw_file.bin", "AB12", "1234ABCD.XX0W"])
def test_number_from_filename_rejects_junk(name):
assert _serial_number_from_bw_filename(name) is None
def test_filename_only_decoder_is_a_guess():
"""It still answers "BE" — that is why it must not be the first choice."""
assert _serial_from_bw_filename("L895K63F.GE0W") == "BE10895"
assert _serial_from_bw_filename("M529LKIQ.G10") == "BE11529"
assert _serial_from_bw_filename("nonsense") is None
# ── the body carries the truth ──────────────────────────────────────────────
def _body(serial: bytes) -> bytes:
return b"\x00" * 32 + b"STRT" + b"\xff\xfe" + serial + b"\x00Geo: 0.254 in/s\x00"
def test_body_wins_for_a_blastmate():
assert _serial_from_bw_bytes(_body(b"BA10895"), "L895K63F.GE0W") == "BA10895"
def test_body_wins_for_a_minimate():
assert _serial_from_bw_bytes(_body(b"BE11529"), "M529LKIQ.G10") == "BE11529"
def test_body_candidate_must_match_the_filename_number():
"""A serial-shaped byte run that disagrees with the filename is ignored."""
assert _serial_from_bw_bytes(_body(b"XX99999"), "L895K63F.GE0W") is None
def test_body_tolerates_a_leading_zero():
assert _serial_from_bw_bytes(_body(b"BA09229"), "K229HGQI.XO0W") == "BA09229"
@pytest.mark.parametrize("data,name", [
(b"", "L895K63F.GE0W"), # no bytes
(_body(b"BA10895"), "junk.bin"), # no derivable number
])
def test_body_returns_none_when_it_cannot_decide(data, name):
assert _serial_from_bw_bytes(data, name) is None
# ── the live monitor-log path ───────────────────────────────────────────────
def _partial_record(serial: bytes) -> bytes:
"""0x2C partial record: type, prefix, two 9-byte timestamps, then ASCII."""
ts = bytes([11, 0x10, 4, 0x07, 0xE9, 0, 16, 2, 0]) # 2025-04-11 16:02:00
return (bytes([0x2C]) + b"\x00" * 10 + ts + ts
+ b"\x00\x00\x00\x00" + serial + b"\x00Geo: 0.254 in/s\x00")
@pytest.mark.parametrize("serial", [b"BE11529", b"BA10895", b"UM11719"])
def test_monitor_log_reads_any_family_prefix(serial):
entry = _decode_0a_partial_header(_partial_record(serial), 0, b"\x01\x11\x00\x00")
assert entry is not None
assert entry.serial == serial.decode()
def test_monitor_log_geo_threshold_survives_a_blastmate():
"""The old find(b"BE") skipped the whole block, losing geo too."""
entry = _decode_0a_partial_header(_partial_record(b"BA10895"), 0, b"\x01\x11\x00\x00")
assert entry is not None
assert entry.geo_threshold_ips == pytest.approx(0.254)
+2 -21
View File
@@ -712,27 +712,8 @@ def test_forty_nn_is_a_data_block_not_a_segment_header():
"""
assert data_block_len(b"\x40\x02\x00\x01\x00\x02", 0) == (6, 2)
assert data_block_len(b"\x40\x08" + bytes(16), 0) == (18, 8)
def test_forty_nn_is_not_capped_at_eight():
"""NN > 8 is a perfectly ordinary `40 NN` block.
This test previously asserted the opposite (`40 0c` -> (None, None)),
codifying a guard that had no evidence behind it: the only corpora
available then used NN in {1,2,3,4,8}, so the cap was never exercised.
Loud UM12947 events use NN of 12, 16, 20 ... up to 196, and rejecting
them halted the block walk mid-record — surfacing as silently short
channels, since the walker stops at the first unrecognised tag rather
than raising. Lifting the cap took that corpus from 22 length-mismatched
files to 0, and 1,476,242 of 1,476,249 samples now reproduce Thor's own
CSV export exactly (the 7 stragglers differ by one 4th-decimal tick).
Verified 2026-09-11; see docs/idf_protocol_reference.md.
"""
assert data_block_len(b"\x40\x0c" + bytes(24), 0) == (26, 12)
assert data_block_len(b"\x40\xc4" + bytes(392), 0) == (394, 196)
# The real bound is the buffer: a block that cannot fit is not a block.
assert data_block_len(b"\x40\xc4" + bytes(8), 0) == (None, None)
assert data_block_len(b"\x40\x00" + bytes(8), 0) == (None, None)
# NN > 8 is not a data block
assert data_block_len(b"\x40\x0c" + bytes(24), 0) == (None, None)
def test_record_chain_is_followed_by_length_not_by_tag_sniffing():