diff --git a/CHANGELOG.md b/CHANGELOG.md index d0ec631..e3d1232 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,275 @@ All notable changes to seismo-relay are documented here. --- +## v0.31.0 — 2026-09-18 + +**Report parity, and a second way to rescue a runaway unit.** Two threads. + +The first closes out Blastware Event/FFT-Report parity: the FFT, the USBM +RI8507 compliance chart and the sensor self-check now render on the event +report, reverse-engineered against BE12844 (MiniMate Plus) and UM (Thor) +events. The sensor check is decoded for **both** series and standardized into +the `.h5` (schema **v2**, a new `/sensor_check` group), so SFM serves it +device-agnostically rather than decoding at report time. The Inspector — an +annotated hex reader for series-3 binaries — is what made the trailing-block +structure findable, and it earned its keep by *ruling out* a stored FFT block +and proving Blastware computes it from the samples. + +The second came out of a field emergency. BE12599's connector fault drove its +Tran channel to its trigger level, so the unit recorded back-to-back and dialed +the office ACH server every ~75 s, unreachable the whole time. +`bridges/ach_server.py` gained `--stop-monitoring` / `--disable-ach` / +`--rescue`, which **invert** the recovery: instead of racing a Stop into the +gaps between dial-outs, point the modem's Destination at our own ACH server and +answer the call. Proven in production the same night — the stop landed on the +first call-in and held. See `docs/runbooks/wedged_unit_recovery.md`. + +⚠ **This release owes prod a backfill** — see Migration below. + +### Added +- **Rescue-on-connect for `bridges/ach_server.py`** — `--stop-monitoring` + (SUB 0x97), `--disable-ach` (SUB 0x2C read → 0x7E write → 0x7F confirm) and + `--rescue` (both). They fire immediately after the startup handshake and + **before** the event walk, so a unit that is recording back-to-back on a + stuck-triggered geophone is quieted as early in the session as possible. + Each action is independently guarded — a failure does not abort the download + — and the outcome is written to `rescue.json` in the session directory. + + This inverts the `docs/runbooks/wedged_unit_recovery.md` approach. That + runbook reaches the unit *inbound* and clears the modem's Destination Address + to stop it dialing. When the device is instead wedged mid-modem-init — ALEOS + logs `tcpmode trying to send to invalid socket` and re-runs `Initialize Auto + answer` every ~75 s, orphaning any held inbound session — inbound cannot win. + Pointing the modem's Destination at an `ach_server` and letting the unit call + *us* gives a device-initiated session the modem bridges properly. + + ⚠ Prefer `--stop-monitoring` alone on first contact. `--disable-ach` stops + the unit calling, which is the only channel to a unit in this state; stopping + the recording ends the call-home loop on its own when ACH is + "after event recorded". + +- **Blastware-compatible channel FFT (`waveform_fft`).** Reproduces Blastware's + FFT Report: DC-removed, no window, zero-padded to 4096 (0.25 Hz bins at + 1024 sps), single-sided `2/N` amplitude. Matches Blastware's dominant + frequency to the exact bin and the amplitude to report precision across all + 28 channels of the 7-event BE12844 oracle set. `channel_spectrum()` / + `dominant_frequency()`; tests in `tests/test_waveform_fft.py`. + +- **USBM RI8507 / OSMRE compliance chart on the event-report PDF + (`sfm/compliance.py`).** The velocity-vs-frequency blasting-compliance + scatter Blastware draws in the upper-right of its Event Report: each channel's + significant cycles as `(frequency, peak velocity)` points (zero-crossing + method, so each channel's cloud tops out at its PPV) plotted against the + RI8507 Drywall (0.75 in/s) and plaster (0.50 in/s) limit curves, drawn + continuous (constant-displacement bounds meeting the plateaus — no vertical + steps). Sized and positioned to match a Blastware report, measured off the + reference PDF. A technical breakdown of the curve is in + `docs/ri8507_compliance_curve.md`. + +- **Sensor self-check waveforms decoded and drawn — both series.** The little + "Sensor Check" traces (geophone ring-downs — the transducer's damped impulse + response — plus a MicL pulse train, the mic's known-signal gain check) are the + unit's proof its sensors were healthy when it recorded the event. + - **Series-3** (`minimateplus.sensor_check`): four records (`0x3c`–`0x3f`) in + the binary's trailing block, same delta-block codec as the main waveform. + Verified against all 7 BE12844 reports (mic zero-crossing = 20.1 Hz exact; + geophone ring-downs ~7.5 Hz, overswing ~3.5). + - **Series-4** (`micromate.sensor_check`): the same self-test in the Thor IDFW + fixed header — four `01 0e 3c/3d/3e/3f` records (same channel ids) storing + raw int16 traces; three-channel (mic-disabled) units carry only the three + geophones. Validated by shape + cross-event consistency. + - **Standardized into the `.h5`** (`/sensor_check`, schema v2): each series' + decoder attaches the traces to the event at decode, the writer persists + them, and `gather_report_data` reads them back — so SFM renders the strip + (flush against the waveform panel) plus the **Sensor Check → Frequency / + Overswing Ratio** sub-rows without knowing the source instrument. + - Tests: `tests/test_sensor_check.py`, `tests/test_sensor_check_idf.py`, + `tests/test_event_hdf5_sensor_check.py`. + +- **Inspector tab in `seismo_lab.py` — annotated hex reader for series-3 + binaries (`minimateplus/binary_annotate.py`).** Tiles a raw Blastware file + into labeled spans (header / STRT / body record-chain / trailing metadata + + calibration + sensor-check records / footer) so a binary can be combed by eye. + +### Fixed +- **Event-report waveform panel — stacked-lane y-tick collision.** The lanes + touch, so each lane's bottom `-1.0` overprinted the next lane's top `1.0` at + the shared boundary. Prune the extreme ticks so each lane shows clean interior + ticks only. +- **Event-report header — serial+firmware line ran off the page.** The long + `BE##### V ##.##-#.## MiniMate Plus` string overflowed the right margin; + tighter right-column indent + BW's slightly smaller header size so it fits. + +--- + +### Migration + +⚠ **The sensor-check needs a backfill.** Existing `.h5` files are schema v1 +and carry no `/sensor_check` group, so their reports show no sensor-check strip +until regenerated. `TOOL_VERSION` is bumped to **0.31.0**, so the standard +backfill regenerates every event and picks up the traces with **no `--force`**: +`scripts/backfill_thor_events.py` for series-4 (it already owed a v0.30.0 Thor +backfill — this rides along) and the series-3 sidecar/shape backfill for +MiniMate events. Purely additive — no decoded value changes, and v1 `.h5` +files read fine until then (empty strip). DB backup first, as always. + +⚠ Budget **~2 h on the NAS** — ~1.5 files/sec there versus ~85/sec on the dev +box (gzip-4 in `sfm/event_hdf5.py` against a Synology CPU). + +Everything else in this release owes nothing: the FFT, the USBM compliance +chart and the `ach_server` rescue flags are additive and read data already on +disk — no schema change, no DB migration. + +--- + +## 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/.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 diff --git a/CLAUDE.md b/CLAUDE.md index 6d4c6ea..67dee54 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -2,7 +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.29.0**. +(Sierra Wireless RV50 / RV55). Current version: **v0.31.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 @@ -24,9 +24,43 @@ 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 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. +- **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/.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`. - **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 @@ -55,6 +89,47 @@ When new information about the protocol is discovered, please update the instant --- +## Changelog & release convention + +**Feature branches do NOT touch `CHANGELOG.md`. Write the entry on `dev`, as +part of finishing the merge, under `## Unreleased`. Cut the version on `dev` in a +dedicated release commit when you are ready to ship to `main`.** + +- **The changelog is written on `dev`, never on a feature branch.** With + several branches in flight they all edit the same few lines at the top of + the file and conflict every time. Writing it once, after the merge, also + lets it describe what actually *landed* — including anything that changed + during conflict resolution. +- ⚠ **The merge is not finished until `## Unreleased` is updated.** Same sitting, + not "later" — that is the one failure mode of writing it after the fact. + Reconstruct from the branch's own commit messages: + `git log --oneline dev..` before you merge, or + `git log --oneline ..` after. +- **No preamble under `## Unreleased`** — just the `### Added` / `### Changed` / + `### Fixed` lists. The themed opening paragraph gets written at release + time, when the whole release is visible and can be named honestly. A theme + written when the first item landed is stale by the third. +- ⚠ **State the operational consequence** on any entry touching the codec, the + waveform store, or the DB — **including when it is "none."** "requires + `backfill_sidecars.py` + `backfill_event_shape.py`, ~2 h on the NAS", + "`TOOL_VERSION` bumped", "no schema change, no migration". Silence is + ambiguous; "none" is information. This repo's changelog is how future-you + learns whether a deploy costs two hours. +- **Releases are cut on judgement, not on a schedule or a merge.** `Unreleased` + is the staging area for whatever is going into the next release; when enough + has accumulated to be worth shipping, it gets a number and a date. Nothing + about a merge to `dev` triggers a release. +- **Cutting a release** is its own `chore(release): vX.Y.Z — ` commit on + `dev`, renaming `## Unreleased` → `## vX.Y.Z — YYYY-MM-DD` and touching: + `CHANGELOG.md`, `pyproject.toml`, the version line in `CLAUDE.md` and + `README.md`, and `minimateplus/event_file_io.py` (`TOOL_VERSION`) **when the + codec changed** — that constant gates `.h5` regeneration. +- **`main` carries only released versions.** No `## Unreleased` section there; + it lands via the `dev` → `main` PR. `main` lagging `dev` by a version is + normal. + +--- + ## Architecture: three-tier conceptual model seismo-relay is a **suite of cooperating components**, not a single app. @@ -120,20 +195,34 @@ 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 (2026-05-28) +#### Thor IDF binary codec (updated 2026-09-10) `micromate/idf_file.read_idf_file()` decodes both Thor IDFW -(waveform) and IDFH (histogram) binaries. +(waveform) and IDFH (histogram) binaries. **Verified per-sample +against Thor's own CSV exports** — see +`scratch/verify_thor_against_csv.py`. -- **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). +- **IDFW** uses the series-3 record-chain `decode_waveform_v2()`. The + body offset is **not** fixed: it is ` + 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. The two outlier `BE9439_*` files in the Thor example corpus are actually Series III Blastware binaries that share the `.IDFW`/`.IDFH` @@ -399,15 +488,19 @@ 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** — 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. +- ~~**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`. ### Decoded sample counts (across the fixture bundle) diff --git a/README.md b/README.md index dfe95e0..7f384f2 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -# seismo-relay `v0.29.0` +# seismo-relay `v0.31.0` A ground-up replacement for **Blastware** — Instantel's aging Windows-only software for managing seismographs. Supports both the **MiniMate Plus diff --git a/bridges/ach_server.py b/bridges/ach_server.py index 0bfd8af..bbc72b8 100644 --- a/bridges/ach_server.py +++ b/bridges/ach_server.py @@ -177,6 +177,8 @@ class AchSession: store: "WaveformStore", clear_after_download: bool = False, restart_monitoring: bool = False, + rescue_stop_monitoring: bool = False, + rescue_disable_ach: bool = False, force_redownload: bool = False, ) -> None: self.sock = sock @@ -190,6 +192,9 @@ class AchSession: self.store = store self.clear_after_download = clear_after_download self.restart_monitoring = restart_monitoring + # Rescue actions for a runaway unit — fired before the event walk. + self.rescue_stop_monitoring = rescue_stop_monitoring + self.rescue_disable_ach = rescue_disable_ach # `force_redownload` tells this session to ignore ach_state and # re-download every event currently on the device, regardless of any # (key, timestamp) match. Useful as a manual override when state has @@ -290,6 +295,41 @@ class AchSession: root_logger.addHandler(fh) try: + # ── Step 1.5: rescue actions ────────────────────────────────────── + # Fired BEFORE the event walk so a runaway unit is quieted as early + # in the session as possible. A unit whose geophone sits above the + # trigger threshold records back-to-back and, with ACH set to "after + # event recorded", re-dials every time — saturating its own firmware + # so it never services inbound requests. See + # docs/runbooks/wedged_unit_recovery.md. + # + # Each action is independently guarded: a failure here must not + # abort the download that follows. + if self.rescue_stop_monitoring or self.rescue_disable_ach: + rescue: dict = {"peer": self.peer, "ts": ts} + + if self.rescue_stop_monitoring: + log.info("Step 1.5: RESCUE — stop monitoring (SUB 0x97)") + try: + client.stop_monitoring() + rescue["stop_monitoring"] = "ok" + log.info(" stop monitoring OK — device should stop recording") + except Exception as exc: + rescue["stop_monitoring"] = f"failed: {exc}" + log.error(" stop monitoring FAILED: %s", exc) + + if self.rescue_disable_ach: + log.info("Step 1.5: RESCUE — disable auto call home (SUB 0x2C/0x7E/0x7F)") + try: + client.set_call_home_config(auto_call_home_enabled=False) + rescue["disable_ach"] = "ok" + log.info(" disable ACH OK — unit should stop calling home") + except Exception as exc: + rescue["disable_ach"] = f"failed: {exc}" + log.error(" disable ACH FAILED: %s", exc) + + _save_json(session_dir / "rescue.json", rescue) + # ── Step 2: device info ─────────────────────────────────────────── device_info = None if not self.events_only: @@ -747,6 +787,13 @@ def serve(args: argparse.Namespace) -> None: print(f" Max events per session: {max_ev if max_ev else 'unlimited'}") print(f" Clear device after download: {'YES' if args.clear_after_download else 'no'}") print(f" Restart monitoring after download: {'YES' if args.restart_monitoring else 'no'}") + _stop_mon = args.stop_monitoring or args.rescue + _dis_ach = args.disable_ach or args.rescue + print(f" RESCUE stop monitoring on connect: {'YES' if _stop_mon else 'no'}") + print(f" RESCUE disable auto call home: {'YES' if _dis_ach else 'no'}") + if _stop_mon and args.restart_monitoring: + print(" !! --restart-monitoring will re-start the unit after download,") + print(" undoing --stop-monitoring. Drop one of them.") print(f" Force re-download all (ignore state): {'YES' if args.force_redownload_all else 'no'}") print(f"{'='*60}") print(f"\n Point your test unit's ACEmanager call-home settings to:") @@ -788,6 +835,8 @@ def serve(args: argparse.Namespace) -> None: store=store, clear_after_download=args.clear_after_download, restart_monitoring=args.restart_monitoring, + rescue_stop_monitoring=args.stop_monitoring or args.rescue, + rescue_disable_ach=args.disable_ach or args.rescue, force_redownload=args.force_redownload_all, ) t = threading.Thread(target=session.run, daemon=True, name=f"ach-{peer}") @@ -862,6 +911,32 @@ def parse_args() -> argparse.Namespace: "DCD on disconnect — without this the unit stays idle after a call-home." ), ) + p.add_argument( + "--stop-monitoring", + action="store_true", + default=False, + help=( + "RESCUE: send SUB 0x97 (stop monitoring) immediately after the " + "handshake, before any event download. Use on a unit that is " + "recording back-to-back because of a stuck-triggered geophone." + ), + ) + p.add_argument( + "--disable-ach", + action="store_true", + default=False, + help=( + "RESCUE: disable Auto Call Home on the device (SUB 0x2C read → " + "0x7E write → 0x7F confirm) immediately after the handshake. The " + "unit stops dialing out until ACH is explicitly re-enabled." + ), + ) + p.add_argument( + "--rescue", + action="store_true", + default=False, + help="Shorthand for --stop-monitoring --disable-ach.", + ) p.add_argument( "--clear-after-download", action="store_true", diff --git a/docs/idf_protocol_reference.md b/docs/idf_protocol_reference.md index aef3c69..2fdc035 100644 --- a/docs/idf_protocol_reference.md +++ b/docs/idf_protocol_reference.md @@ -6,7 +6,15 @@ 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. -**Status (2026-05-28):** ASCII text sidecar fully decoded (1,014 +> ⚠ **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 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, @@ -44,6 +52,220 @@ 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: + +``` +/UM13981_20220207084555.IDFW +/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 ` + 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 +` 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 diff --git a/docs/ri8507_compliance_curve.md b/docs/ri8507_compliance_curve.md new file mode 100644 index 0000000..befb1a3 --- /dev/null +++ b/docs/ri8507_compliance_curve.md @@ -0,0 +1,135 @@ +# USBM RI8507 / OSMRE Blasting Compliance Curve — Reference + +Reference for the **velocity-vs-frequency blasting compliance chart** Blastware +draws on its Event Report ("USBM RI8507 And OSMRE"), and how seismo-relay +reproduces it. Implemented in [`sfm/compliance.py`](../sfm/compliance.py); the +spectral (FFT) side lives in [`waveform_fft.py`](../waveform_fft.py). + +Reverse-engineered 2026-09-14 against 7 BE12844 (MiniMate Plus) events, each +with a Blastware Event Report + FFT Report as ground truth. Curve values from +USBM RI8507 Appendix B and 30 CFR 816.67. + +--- + +## What it is + +Two closely-related sources for the same limit curve: + +- **USBM RI8507** — Bureau of Mines *Report of Investigations 8507* (Siskind + et al., 1980), *"Structure Response and Damage Produced by Ground Vibration + From Surface Mine Blasting."* The curve is **Figure B-1**, Appendix B + ("Alternative Blasting Level Criteria"), p.73–74. +- **OSMRE / OSM** — the Office of Surface Mining Reclamation and Enforcement + codified it as **30 CFR 816.67, Figure 1**. "CFR" = the U.S. Code of Federal + Regulations. Same curve, regulatory force. + +The chart plots each geophone channel's significant vibration cycles as +`(frequency, peak velocity)` points against this limit. A point **below** the +line passes; **above** fails. + +--- + +## The limit curve + +A structure has a resonance band (~4–12 Hz for whole structures) where it is +most vulnerable, so the safe velocity is **lower** at those frequencies and +**higher** away from them. The curve captures this by alternating two kinds of +bound: + +- **Constant-velocity** segments — a flat horizontal line at a fixed PPV. +- **Constant-displacement** segments — a fixed peak *displacement* `d`. For + simple harmonic motion, peak velocity `v = 2πf·d`, so on a velocity-vs- + frequency **log-log** plot this is a straight line of slope +1 (velocity rises + with frequency). This is why the low- and high-frequency bounds are sloped. + +### Two lines — structure type + +RI8507 gives two lines for two interior-wall constructions (Table 13, p.67): + +| line | construction | plateau PPV | +|---|---|---| +| **Drywall** (solid) | modern gypsum wallboard | **0.75 in/s** | +| **Plaster** (dashed) | older plaster on wood lath | **0.50 in/s** | + +Plaster-on-lath is more damage-prone, hence the lower limit. You apply **one** +line depending on the monitored structure. + +### The four segments (Figure B-1, p.74) + +Going low → high frequency, each line is: + +1. **Ultimate low-frequency bound** — constant displacement **0.030 in** + (`v = 2πf·0.030`). Only relevant below ~4 Hz. +2. **Plateau** — constant velocity **0.75** (Drywall) / **0.50** (plaster) in/s. +3. **Rising diagonal** — constant displacement **0.008 in** (`v = 2πf·0.008`), + climbing from the plateau up to the high-frequency cap. +4. **High-frequency cap** — constant velocity **2.0 in/s** above ~40 Hz. + +The segments are drawn **continuous**: each bound is used over the frequency +range where it is the binding (lowest) limit, and consecutive bounds meet where +they are equal — so there are no vertical steps. Transition frequencies come +straight from the values (`f = V / (2π·d)`): + +| transition | formula | Drywall | Plaster | +|---|---|---|---| +| 0.030 in → plateau | `V_mid / (2π·0.030)` | 3.98 Hz | 2.65 Hz | +| plateau → 0.008 in | `V_mid / (2π·0.008)` | 14.92 Hz | 9.95 Hz | +| 0.008 in → 2.0 in/s | `2.0 / (2π·0.008)` | 39.79 Hz | 39.79 Hz | + +Because both lines share the same **0.008 in** rising diagonal, above ~15 Hz +they lie on the *same* line (both reach 2.0 in/s at ~40 Hz) — RI8507's literal +construction merges them there. Blastware renders the dashed line as a separate +parallel diagonal, but that is cosmetic: above ~15 Hz both structure types carry +the identical limit, so compliance is unaffected. + +> ⚠ RI8507's *Table 13* is a simpler two-range criterion with a **sharp +> discontinuity at 40 Hz** (flat plateau, then a jump to 2.0). Figure B-1 is the +> **smoothed** version that adds the 0.008 in transition — that is the one drawn +> on reports and implemented here. + +--- + +## The compliance scatter (the points) + +The cloud is **not** the FFT spectrum. It is a per-cycle, time-domain measure by +the **zero-crossing method** (`channel_compliance_points`): + +- Split the channel's waveform at its zero crossings. +- Each half-cycle contributes one point: **frequency** `= 1 / (2 · half-period)` + (from the samples between the two crossings), **velocity** `= peak |amplitude|` + in that half-cycle. + +This yields ~90–110 points per channel, and — by construction — each channel's +**highest** point equals that channel's PPV. Verified against Blastware: the +cloud shape, density, and ceiling all match. + +### Why not the FFT? + +A broadband blast spreads its energy across many FFT bins, so no single bin +reaches the time-domain peak — the FFT amplitudes come out ~10× below the +compliance-chart velocities. The compliance chart is a *per-cycle peak* view; +the **FFT** is a separate analysis (Blastware's *FFT Report*), reproduced by +[`waveform_fft.py`](../waveform_fft.py) and used for the dominant-frequency +readout and the #10 FFT view — not for this scatter. + +--- + +## Implementation + +- `sfm/compliance.py` + - `limit_at(freq, curve)` — the limit PPV at a frequency (`curve` = `"Drywall"` + or `"Plaster"`); curves are data in `_CURVES`, so more standards can be added. + - `channel_compliance_points(samples, sps)` — the zero-crossing scatter. + - `draw_compliance_chart(ax, channels, sps)` — matplotlib rendering (both + limit lines + per-channel scatter, Blastware's tick scales and channel + markers: Tran `+` red, Vert `×` green, Long `o` blue). +- Tests: `tests/test_compliance.py`. + +--- + +## Sources + +- USBM **RI8507** (Siskind, Stagg, Kopp, Dowding, 1980), Appendix B / Figure B-1, + p.73–74; Table 13, p.67. (`ref-stuff/usbm-ri8507-ground_vibration.pdf`.) +- **30 CFR 816.67**, "Use of explosives: Control of adverse effects," Figure 1 — + diff --git a/docs/runbooks/wedged_unit_recovery.md b/docs/runbooks/wedged_unit_recovery.md index 8d27dd0..bf427a0 100644 --- a/docs/runbooks/wedged_unit_recovery.md +++ b/docs/runbooks/wedged_unit_recovery.md @@ -1,6 +1,7 @@ # Runbook — Recovering a wedged unit stuck in a call-home loop -**Original incident:** BE9558H at `166.246.130.1:9034`, recovered 2026-05-17. +**Incidents:** BE9558H at `166.246.130.1:9034`, 2026-05-17 (Method B) · +BE12599 at `166.246.64.226:9034`, 2026-09-16 (Method A). A field unit with a stuck-triggered geophone (or any hardware fault causing constant event triggering) will record events back-to-back, and if Auto Call @@ -14,6 +15,33 @@ This runbook describes how to break the loop and recover control. --- +## ⚠ Two cures for one disease — intercept first + +Both incidents below are the **same failure**: a geophone offset crosses the +trigger level, the unit records back-to-back, ACH set to "after event +recorded" dials continuously, and the unit becomes unreachable because its +modem is in client mode almost all of the time. + +There are two ways to get a Stop Monitoring command into it. + +| | **A — intercept the call** (preferred) | **B — catch it between calls** (original) | +|---|---|---| +| Idea | Be the server it dials. Point the modem's Destination at our own ACH server and answer it. | Clear the Destination so it stops dialing, then race a Stop into the gap. | +| Needs inbound? | **No — the unit calls us** | Yes: working inbound TCP to the modem | +| Determinism | Deterministic — it dials every ~75 s, we only have to be listening | A race. BE9558H took ~7 h of attempts before one landed. | +| Tool | `bridges/ach_server.py --stop-monitoring` | `scripts/slow_drip.sh` | +| Proven on | BE12599, 2026-09-16 | BE9558H, 2026-05-17 | + +**Method A is the standard procedure now.** The unit won't answer us because +it is on the phone — so stop dialing it and be the one it calls. It rings, +we pick up, take its data, and tell it to stop calling here. + +Method B is kept because it is proven, and because A needs a listener the +modem can actually reach (public IP + forwarded port). When you have that, +don't race it — intercept it. + +--- + ## Symptoms - Terra-View / SFM `/device/info` either hangs or fails on `count_events()`. @@ -31,9 +59,85 @@ If you see *all* of these, the unit is in this exact failure mode. --- -## Quick reference — how to recover +## Method A (preferred) — intercept the call -You need **ACEmanager access** to the unit's modem. +You need **ACEmanager access** and a host the modem can dial: public IP with +the listener's port forwarded to it. + +### A1 — start the listener BEFORE touching the modem + +```bash +cd /home/serversdown/seismo-relay +tmux new -s rescue +.venv/bin/python -u bridges/ach_server.py --port 12345 \ + -o bridges/captures/-diag --stop-monitoring -v +``` + +⚠ **Listener first, always.** A Destination pointed at a dead port is the +worst state available — the device still dials, the modem still flips to +client mode, inbound stays blocked, and nothing gets delivered. + +Do **not** add `--events-only` (it silently breaks dedup — see gotchas), and +do **not** add `--disable-ach` yet (see A4). + +### A2 — point the modem at it + +ACEmanager → **Serial → Port Configuration**: + +| Field | Set to | +|---|---| +| **Destination Address** | the listener's public IP | +| **Destination Port** | the listener's port (e.g. `12345`) | + +Apply. The modem auto-dials its Destination whenever serial data arrives +while the serial port is closed — so the unit's own retry cycle now lands on +you instead of nowhere. + +### A3 — answer, and stop the bleeding + +Within ~75 s you should see a call-in. `--stop-monitoring` fires SUB 0x97 at +step 1.5 — after the handshake, **before** the event walk — so the recording +halts at the earliest possible moment in the session. Confirm via +`rescue.json` in the session directory: + +```json +{"peer": "166.246.64.226:60921", "stop_monitoring": "ok"} +``` + +That is the bleeding stopped. Everything after this is cleanup. + +### A4 — drain the backlog, THEN disable ACH + +⚠ **Order matters, and it is counter-intuitive.** Stopping monitoring also +removes your call-in trigger: ACH fires on "after event recorded", so with +recording stopped the unit has no reason to dial again. The backlog sitting +in its memory does **not** re-arm it. + +So if the stored events are worth keeping — and on a fault unit they usually +are, they're the evidence — drain them across however many call-ins it takes +*before* you silence it. Only then add `--disable-ach` (or use +`scripts/rescue_device.sh --no-erase`). + +If the unit has gone quiet and you still need it, cycling the modem produces +a call-in, and a unit with a scheduled daily call will dial at its configured +time regardless. + +### A5 — restore the Destination, and confirm you did + +Put `Destination Address` back to `0.0.0.0` (or the office Instantel ACH +server) once you are finished, and only stop the listener after that is done. + +### A6 — do NOT re-enable ACH until the hardware fault is repaired + +Otherwise the loop restarts the moment monitoring resumes and you run this +runbook again. + +--- + +## Method B (fallback) — catch it between calls + +The original 2026-05 procedure. Use when you cannot stand up a listener the +modem can reach. You need **ACEmanager access** to the unit's modem. ### Step 1: stop the modem's mode-flipping @@ -253,3 +357,223 @@ service). Total time from "i was wondering if its possible to" first attempt to recovery: ~7 hours of intermittent debugging across one evening. + +--- + +# Second incident — BE12599, 2026-09-16/17 + +**Unit:** BE12599 at `166.246.64.226:9034`, RV50, job *I-80 North Fork Bridge +— Abut 1 West* (Fay Company). Same job as BE9558H, which is a coincidence. + +**Fault:** the connector fault documented in `docs/offset_investigation.md` +§8e progressed until the Tran pedestal reached **0.400 in/s** — its trigger +level. Constant triggering → constant recording → ACH "after event recorded" +→ continuous dialing. Same disease as BE9558H. + +**Same disease, inverted cure.** Method B's Step 1 *did* work — clearing the +Destination stopped the dial-outs, confirmed in the ALEOS log. It was Step 2 +that didn't land, and rather than keep racing we turned the rescue around: +gave the unit a different server to call, and answered it. + +Total time ≈ 5 h, of which ~90 min went to two red herrings documented below. +Much of the rest was rediscovering the May procedure, which is why the +"two cures" table now sits at the top of this file. + +--- + +## Turn on ALEOS_SERIAL debug FIRST + +This is the single highest-value diagnostic and it should be step zero on any +future incident. ACEmanager → **Admin → Log → ALEOS_SERIAL log level → +DEBUG**, then view the serial log. + +It is the only thing that tells you what the *device* is actually saying. +Everything before we did this was guesswork. + +## What the log showed — the unit is on the phone + +Every ~75 seconds, verbatim: + +``` +ALEOS_SERIAL_HIF: 29 byte(s) in buffer: 'ATQ1^MATE0^MATS0=2^M^MRADIO RING^M' +ALEOS_SERIAL_HMC: TCP recvhost fd 65535 len 29 state TCPMode::kClosed +ALEOS_SERIAL_HMC: tcpmode trying to send to invalid socket +ALEOS_SERIAL_HMC: Connect to IP: 0.0.0.0 Port 0 +ALEOS_SERIAL_HMC: Initialize Auto answer on port 9034 +ALEOS_SERIAL_HMC: Cannot connect to 0.0.0.0 +``` + +Read that carefully: + +- `ATQ1` (quiet) / `ATE0` (echo off) / `ATS0=2` (auto-answer after 2 rings). + **There is no `ATD`.** The device is not dialing — it is trying to + *configure* its modem. +- The modem's serial port is in TCP data mode, so it never interprets these + as AT commands. It treats them as payload and tries to ship them to a TCP + socket that does not exist. +- The device therefore never receives `OK`, never progresses, and **retries + the identical 29 bytes forever**. + +**While it is in this state it is busy placing a call, not listening for +us.** This is almost certainly what BE9558H was doing too — we simply never +turned on ALEOS_SERIAL debug in May to look. It is not a different disease; +it is the same one, seen properly for the first time. + +It is also the argument for Method A in one picture: the unit is mid-dial +every ~75 s, and our inbound Stop has to thread the gaps between those +attempts. Give it somewhere to dial and the problem inverts into a +deterministic one. + +### Why `slow_drip` lied + +`slow_drip` returned the *success* signature except for the one field that +mattered: + +```json +{"duration_s":120.0,"drips_sent":38,"bytes_sent":920, + "bytes_received":0,"send_error":null} +``` + +Full duration, no broken pipe — but zero bytes back. Cause is in the log +above: each 75 s cycle re-runs `Initialize Auto answer on port 9034`, which +orphans the held session (`data in for unknown reason 3 removing from +select`, `OnMsg recv error: 107 - Transport endpoint is not connected`). Our +local TCP stayed open so `sendall` never raised — but the modem stopped +bridging after the first re-init, so every drip after that went into a socket +nobody was reading. + +⚠ **`send_error: null` + full duration is NOT success. Only +`bytes_received > 0` is success.** + +⚠ **In fairness to slow_drip: it got exactly one attempt here**, run ~90 s +after a modem reboot, with a dead session visible in the log at 20:19:17 in +that same window. BE9558H took hours of attempts before one landed. Method B +was not ruled out on BE12599 so much as abandoned in favour of something that +doesn't need luck. + +--- + +## ⚠ Two red herrings that cost ~90 minutes + +### 1. The trusted-IP whitelist (this was the real reason inbound never worked) + +The RV50s run with **Security → Trusted IPs (Friends List) enabled**. A +source IP that is not on the list is dropped **silently** — inbound presents +as `Connection error: timed out`, never a refusal. + +Brian's dev-box public IP is **dynamic** and had changed, so `tmi-dev` was no +longer whitelisted. Every inbound attempt failed identically across four +different modem and device states, which looked exactly like the BE9558H +mode-flipping symptom and sent us chasing modem configuration for over an +hour. + +**Check this before diagnosing anything else.** Note that SFM in Docker +egresses via the *host's public IP*, not its LAN IP. + +### 2. A 502 from SFM does not mean TCP connected + +`sfm/server.py` raises **502 for both** failure classes: + +```python +raise HTTPException(status_code=502, detail=f"Protocol error: {exc}") +raise HTTPException(status_code=502, detail=f"Connection error: {exc}") +``` + +We read an early 502 as "TCP connected, modem bridged, device mute" and built +a whole theory on it. It was almost certainly a connect timeout. +**Always read the `detail` string** — "connect failed" and "device didn't +answer" are completely different problems and the status code will not +separate them. + +--- + +## What actually worked — invert the direction + +The key observation is in the log above: + +> `TCP recvhost ... state TCPMode::kClosed` → `Connect to IP: 0.0.0.0 Port 0` + +**The modem auto-dials its Destination whenever serial data arrives while +closed.** So instead of fighting for inbound, give it somewhere to dial: +point `Destination Address` at our own `ach_server` and the device's own +75-second attempts become **device-initiated sessions the modem bridges +correctly**. No race, no contention, worst case a 75-second wait. + +### Procedure + +1. **Run the rescue server** on a host the modem can reach (public IP + + forwarded port): + + ```bash + cd /home/serversdown/seismo-relay + .venv/bin/python -u bridges/ach_server.py --port 12345 \ + -o bridges/captures/-diag --stop-monitoring -v + ``` + +2. **Point the modem at it** — ACEmanager → Serial → Port Configuration → + `Destination Address` = your public IP, `Destination Port` = 12345. + +3. **Wait for the call-in.** `--stop-monitoring` fires SUB 0x97 at step 1.5, + after the handshake and *before* the event walk. Confirm via + `rescue.json` in the session directory: + + ```json + {"peer": "166.246.64.226:60921", "stop_monitoring": "ok"} + ``` + +4. **Restore the modem's Destination** once you are done, then finish the + device side (disable ACH, erase) through whichever channel works. + +On BE12599 the first call-in landed at 20:58:11 and reported +`stop_monitoring: ok`; a second at 20:58:20 confirmed it. `is_monitoring: +false` was still true **6½ hours later** — the fix is durable. + +--- + +## Hard-won gotchas (do not re-derive) + +- **Never leave the Destination pointed at a host with nothing listening.** + That is the worst state available: the device still dials, the modem still + flips, inbound stays blocked, and nothing is delivered. An 8-minute gap + with the listener down produced a spurious inbound timeout that cost + another round of misdiagnosis. + +- **Stopping monitoring removes your call-in channel.** ACH is "after event + recorded"; no new events means no new dials. The backlog sitting in memory + does *not* re-arm it. After a successful stop the unit goes quiet and you + need the modem cycled (works — produced a call-in), the scheduled daily call + (BE12599 calls at **05:00:14 device-local**, per §8e), or working inbound. + **Plan the order before you fire the stop.** + +- **`--events-only` silently breaks dedup.** It skips the device-info step, + so the serial is never read; `ach_state.json` then keys on + `peer:ephemeral_port`, which is unique per connection. Every session looks + like a new unit, starts from key 0, and re-downloads the same event. Four + sessions on BE12599 downloaded the identical event four times and made zero + progress on the backlog. Events also file as `serial=UNKNOWN` with a + `M000…` BW filename (serial_numeric 0) instead of `N599…`. + **Do not use `--events-only` when you intend to download anything.** + +- **`/device/events/index` reported `lifetime_count: 0`** on a unit with years + of history. Suspected decode bug in the SUB 0x08 field offset — do not + trust that number. The 88-byte payload is preserved in the `raw_hex` field + if someone wants to chase it. + +- **Memory used cross-checks the event keys exactly:** + `last_key − buffer_start = memory_total − memory_free`. On BE12599: + `0x011230ec − 0x01110000 = 78,060` and `983,028 − 904,968 = 78,060`. + Useful sanity check that you are reading the keys right. + +--- + +## Final state (2026-09-17 ~01:30 local) + +- `is_monitoring: false`, held 6½ hours +- Battery 6.76 V +- Memory 78,060 / 983,028 bytes used (8%) +- `first_key 01121728`, `last_key 011230ec` — ~6.6 KB of addressable event + chain, roughly 3 events +- ACH still **enabled** — to be disabled after the backlog is preserved +- Modem Destination still pointed at tmi-dev — to be restored +- ⚠ **Do not re-enable ACH until the connector is serviced.** Tran is still + sitting at 0.400 and the loop restarts the moment monitoring resumes. diff --git a/docs/superpowers/plans/2026-09-17-rescue-listener.md b/docs/superpowers/plans/2026-09-17-rescue-listener.md new file mode 100644 index 0000000..58a5de4 --- /dev/null +++ b/docs/superpowers/plans/2026-09-17-rescue-listener.md @@ -0,0 +1,134 @@ +# Plan — "Rescue Listener": a first-class tool for the inverted rescue + +**Status:** proposal, not started. Written 2026-09-17 ~01:40 local, straight +off the BE12599 incident. Open questions at the bottom need Brian's answer +before anything is built. + +**Background:** `docs/runbooks/wedged_unit_recovery.md`, "Second incident — +BE12599". The manual version of this worked; this plan is about making it a +tool instead of a sequence of remembered steps at 1 AM. + +--- + +## The problem, stated plainly + +When a unit is wedged in the BE12599 mode — geophone offset above trigger, +recording back-to-back, ACH dialing constantly, device stuck repeating an AT +modem-init string and therefore **deaf to S3 over inbound** — the only channel +that works is the one the *device* opens. + +Recovering it currently means: + +1. Remember that `bridges/ach_server.py` exists and takes the right flags +2. Start it by hand on a box the modem can reach, with a public port forwarded +3. Go into ACEmanager and repoint the modem's Destination +4. Watch a terminal for a call-in +5. Read `rescue.json` to find out whether it worked +6. Go back into ACEmanager and repoint the modem to where it belongs +7. **Not forget step 6**, because leaving the Destination pointed at a dead + listener is worse than never having started + +That is six manual steps and one landmine, executed under pressure while a +unit floods the office server. + +## What the tool should be + +**A "rescue listener" an operator can start for one unit, which handles +whatever that unit says when it calls in, and refuses to go away until the +operator confirms the modem has been pointed back.** + +Lifecycle: + +1. **Start** — operator names the target unit and starts a rescue listener. + The tool reports the exact address/port to enter in ACEmanager, plus the + actions it will take. +2. **Operator repoints the modem** to that address. +3. **Wait** — listener sits there. Live status: "waiting for call-in", + elapsed, last-seen. +4. **Act** — on call-in, run the configured rescue actions automatically, + in a safe order, each independently guarded. Report per-action outcome. +5. **Hold** — the listener **stays up** and keeps reporting, because the + modem is still pointed at it. +6. **Confirm & stop** — the operator explicitly confirms the Destination has + been restored (to `0.0.0.0`, or to the office Instantel ACH server). + Only then does the listener shut down. + +Step 6 is the whole point of making this a tool. It is the step that is +easiest to skip and most expensive to skip. + +## Default action set + +Ordered deliberately — see "order matters" below. + +| # | Action | Default | Why | +|---|---|---|---| +| 1 | **Stop monitoring** (SUB 0x97) | ✅ on | Halts recording; ends the trigger→record→dial loop at its source. Already implemented as `--stop-monitoring`. | +| 2 | **Drain events** to a diagnostics store | ⚙ configurable | The backlog is usually evidence, not garbage — see the BE12599 offset investigation. Must NOT land in the prod SFM DB. | +| 3 | **Disable ACH** (SUB 0x2C/0x7E/0x7F) | ❌ off by default | Stops the dialing — **and stops your only channel**. Opt-in, and ideally gated on step 1 having succeeded. | +| 4 | **Erase events** | ❌ off by default | Destructive. Only after a verified drain. | + +### Order matters — the lesson from BE12599 + +Stopping monitoring *removes the call-in trigger*. ACH fires on "after event +recorded"; with recording stopped, the unit has no reason to dial again, even +though the backlog is still sitting in its memory. So a naive +"stop + disable + erase, all at once" rescue can silence the unit before +you've collected anything, leaving you with no channel and a device full of +evidence. + +The tool should either sequence around this or warn loudly about it. My +instinct is: **stop monitoring immediately** (it's the bleeding), then drain +across however many call-ins it takes, and treat disable-ACH/erase as a +separate, explicit "finish" action once the operator is satisfied. + +## Where it should live — open question, with a proposal + +The natural tier is **SFM** (device-side, per the three-tier model in +CLAUDE.md). But the rescue listener must be reachable *from the cellular +network*, which is a deployment constraint SFM's usual profile doesn't have. + +**Proposal worth considering:** run it at the office, beside the real Instantel +ACH server, on a **different port** (e.g. 12346 while Instantel holds 12345). +Then the ACEmanager change is a **port change, not an IP change** — smaller, +faster, less to get wrong, and trivially reversible. It also means the office +public IP (already stable and known) is the destination, rather than whatever +Brian's dynamic home IP happens to be that week. + +The tmi-dev approach used on BE12599 worked, but required a router forward and +ran into the dynamic-IP problem in the same session. + +## Open questions + +1. **Where does it run?** Office beside Instantel ACH (port swap), SFM on the + NAS, or ad-hoc on tmi-dev? Affects everything else. +2. **What drives it?** Terra-View admin page (fits "operator UI"), an SFM + endpoint pair (`POST /device/rescue_listener/start` + `/stop` + `/status`), + or a CLI wrapper? A long-lived listener doesn't fit the request/response + endpoint shape well — probably needs a background task with a status poll. +3. **How does it identify the unit?** It can't know the serial until the + device calls in and the handshake reads it. Allowlist by modem IP? Accept + anything and report what showed up? +4. **Where do drained events go?** A per-incident diagnostics store + (`bridges/captures/-diag`) seems right — explicitly *not* the prod + SFM DB. Does that store need to be a first-class thing with its own + retention, or is a directory fine? +5. **How is "confirm the modem is repointed" verified?** Operator attestation + (a button), or can we actually probe it? If the listener stops seeing + call-ins that's weak evidence; if inbound to the unit starts working that's + stronger. +6. **Multi-unit?** One listener per incident, or one listener that handles any + unit that dials in? Probably the former for safety. +7. **Timeout / abandonment policy.** If nobody ever confirms, does it run + forever? Alert after N hours? + +## What already exists + +- `bridges/ach_server.py` — the listener itself, with `--stop-monitoring`, + `--disable-ach`, `--rescue` (added on `feat/ach-rescue-on-connect`, commit + `9f1050b`), `--clear-after-download`, `--max-events`, `--allow-ip`. +- Per-session `rescue.json` recording per-action outcomes. +- Isolated per-output-dir SQLite + waveform store, so a diagnostics capture is + already separate from prod by construction. + +So the gap is not protocol work — it's lifecycle, operator surface, and the +confirmation gate. Most of the risk is in questions 1 and 2. diff --git a/micromate/idf_file.py b/micromate/idf_file.py index 60937b8..b862e03 100644 --- a/micromate/idf_file.py +++ b/micromate/idf_file.py @@ -47,19 +47,24 @@ from dataclasses import dataclass from pathlib import Path from typing import Optional, Union -# Thor IDFW bodies are pinned to the SUPERSEDED tag-dispatch decoder. +# Thor IDFW bodies use the series-3 record-chain decoder. # -# _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, -) +# 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 from .models import IdfEvent, IdfPeaks, IdfReport @@ -89,23 +94,70 @@ _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). -_BODY_SCAN_FLOOR = 0x0E00 +# 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 -# 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 +# 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 # 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. -_IDFH_INTERVAL_SIZE = 72 # bytes per per-interval record +# 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_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") @@ -223,26 +275,67 @@ def _find_waveform_body_offset(buf: bytes) -> Optional[int]: """ if len(buf) < _BODY_SCAN_FLOOR + 8: return None - best: Optional[tuple[int, int]] = None # (total_samples, offset) - i = _BODY_SCAN_FLOOR - while True: - j = buf.find(_BODY_MAGIC, i) - if j < 0: - break - i = j + 1 + + # 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 `` 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 + while True: + j = buf.find(sig, 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 - if best is None or total > best[0]: - best = (total, j) - return best[1] if best else None + # >= 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 def _decode_waveform_samples(buf: bytes) -> Optional[dict]: @@ -299,6 +392,12 @@ 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") @@ -307,7 +406,11 @@ class IdfhInterval: def peak_ips(self, channel: str) -> float: """Convert peak count to in/s (geo channels only).""" - return self.peak_count(channel) / _IDFH_INT16_FS * _IDFH_GEO_FULL_SCALE + # 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 def freq_hz(self, channel: str) -> Optional[float]: halfp = getattr(self, f"{channel.lower()}_halfp") @@ -316,11 +419,46 @@ class IdfhInterval: return _IDFH_HALFP_FREQ_NUM / halfp -def _decode_idfh_interval(buf72: bytes, offset: int) -> IdfhInterval: - """Decode one 72-byte interval record into per-channel min/max/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. + """ 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] @@ -336,6 +474,7 @@ def _decode_idfh_interval(buf72: bytes, offset: int) -> IdfhInterval: 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, ) @@ -343,36 +482,73 @@ 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][00 NN_counter][05 3f]`` where ``length`` + ``[length_be 2B][0a 00 00 00][counter_be 2B][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). - Confirmed against the 859-file corpus (181,071 intervals decoded; 1 - failure is the sig-B BE9439 file). + + ``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. """ 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][00 NN][05 3f] - if buf[j + 4] != 0x00 or buf[j + 6 : j + 8] != b"\x05\x3f": + # 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": i = j + 1 continue length = int.from_bytes(buf[j - 2 : j], "big") - n = (length - _IDFH_SEGMENT_HEADER) // _IDFH_INTERVAL_SIZE + 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 if n <= 0: i = j + 1 continue - header_start = j - 2 + 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 interval_start = header_start + _IDFH_SEGMENT_HEADER for k in range(n): - off = interval_start + k * _IDFH_INTERVAL_SIZE - if off + _IDFH_INTERVAL_SIZE > len(buf): + off = interval_start + k * stride + if off + stride > len(buf): break - chunk = buf[off : off + _IDFH_INTERVAL_SIZE] - intervals.append(_decode_idfh_interval(chunk, off)) + 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 # Advance past this segment + the 2-byte tail. i = header_start + length + _IDFH_SEGMENT_TAIL return intervals @@ -452,7 +628,12 @@ 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. - mic_peak_count = max((iv.peak_count("MicL") for iv in intervals), default=0) + # 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_psi = mic_count_to_psi(mic_peak_count) if mic_peak_count else None rep = IdfReport( serial_number=md.serial, diff --git a/micromate/sensor_check.py b/micromate/sensor_check.py new file mode 100644 index 0000000..dea67e6 --- /dev/null +++ b/micromate/sensor_check.py @@ -0,0 +1,89 @@ +r"""Decode the Thor / Micromate (series-4) sensor self-check waveforms from an +IDFW event binary. + +Reverse-engineered 2026-09-15 against 4 UM (Thor) oracle events. The IDFW +binary carries the sensor self-check in its fixed-header region (before the +waveform body), as up to four records tagged ``01 0e 3c/3d/3e/3f`` — the SAME +channel ids as the series-3 MiniMate Plus (Tran / Vert / Long / MicL), which is +the physical self-test: + + * 3c / 3d / 3e = Tran / Vert / Long geophone ring-downs (a damped impulse + response — resonant frequency + damping). + * 3f = MicL pulse train (the mic's known-signal gain check). Absent + on three-channel (mic-disabled) units. + +Record framing (per record):: + + 01 0e [id:1] [flags:3] [count:2 BE] [pad:10] [int16-BE samples × count] + \___ 18-byte header ___/ + +Unlike series-3's delta-coded trailing block, series-4 stores each trace as a +raw int16 big-endian array. ``count`` (the 2-byte field at header offset +8) +is the sample count; the record is padded to a fixed stride after that. +""" +from __future__ import annotations + +import struct +from typing import Dict, List + +# Record id → channel. Same ids/order as series-3 (minimateplus.sensor_check). +_ID_TO_CHANNEL = {0x3C: "Tran", 0x3D: "Vert", 0x3E: "Long", 0x3F: "MicL"} +_CHAIN_IDS = (0x3C, 0x3D, 0x3E, 0x3F) + +_MARKER = b"\x01\x0e" # precedes the 1-byte channel id +_HEADER_LEN = 18 # bytes from the marker start to the first sample +_COUNT_OFF = 8 # 2-byte BE sample count, from the marker start +_MAX_COUNT = 4000 # sanity cap (traces are ~70-200 samples) + + +def _find_chain(raw: bytes): + """Locate the sensor-check record chain. Returns a list of + ``(offset, id, count)`` for the first run of markers whose ids run + 3c, 3d, 3e[, 3f] in order, or ``[]``. + + Records are padded to a fixed stride, so the next marker is not at + ``header + count*2``; instead collect every ``01 0e [id]`` marker with a + sane count and take the first id-ordered run. Validating the id sequence + (not a lone ``01 0e 3c``) keeps a stray marker in the waveform body from + matching — the real chain sits in the fixed header, ahead of the body. + """ + n = len(raw) + markers = [] + for p in range(n - _HEADER_LEN): + if raw[p:p + 2] == _MARKER and raw[p + 2] in _ID_TO_CHANNEL: + count = int.from_bytes(raw[p + _COUNT_OFF:p + _COUNT_OFF + 2], "big") + if 0 < count <= _MAX_COUNT: + markers.append((p, raw[p + 2], count)) + + for i, (off, rid, _c) in enumerate(markers): + if rid != 0x3C: + continue + run = [markers[i]] + for m in markers[i + 1:]: + if len(run) < len(_CHAIN_IDS) and m[1] == _CHAIN_IDS[len(run)]: + run.append(m) + else: + break + if len(run) >= 3: # 3-channel (mic-disabled) units are valid + return run + return [] + + +def decode_idf_sensor_check(raw: bytes) -> Dict[str, List[int]]: + """Decode the sensor self-check traces from a Thor/Micromate IDFW binary. + + Returns ``{"Tran": [...], "Vert": [...], "Long": [...], "MicL": [...]}`` in + raw int16 ADC counts (MicL omitted on 3-channel units), or ``{}`` if the + binary carries no sensor-check chain (a non-IDF file, or an IDFH histogram). + """ + chain = _find_chain(raw) + if not chain: + return {} + out: Dict[str, List[int]] = {} + for off, rid, count in chain: + start = off + _HEADER_LEN + blob = raw[start:start + count * 2] + if len(blob) < count * 2: + continue + out[_ID_TO_CHANNEL[rid]] = list(struct.unpack(">%dh" % count, blob)) + return out diff --git a/minimateplus/binary_annotate.py b/minimateplus/binary_annotate.py new file mode 100644 index 0000000..25e164d --- /dev/null +++ b/minimateplus/binary_annotate.py @@ -0,0 +1,75 @@ +"""Structural annotation of a Series-3 Blastware waveform binary. + +Pure, no I/O: takes the raw file bytes and returns a flat, gap-free tiling of +labelled :class:`Span` regions for a hex viewer to paint. Every byte is +covered — anything the decoder can't account for becomes an ``unknown`` span, +so undecoded regions (e.g. a stored spectral/FFT block, if one exists) stand +out instead of hiding. + +File layout (see ``blastware_file.py``): ``[header][21B STRT][body][26B footer]``. +The body is the record chain walked by :func:`waveform_codec.walk_records`. +""" +from __future__ import annotations + +from dataclasses import dataclass +from typing import List + +from .waveform_codec import walk_records + +_STRT_LEN = 21 +_FOOTER_LEN = 26 + + +@dataclass +class Span: + start: int # inclusive byte offset + end: int # exclusive byte offset + label: str # human-readable description + kind: str # 'header' | 'strt' | 'sample' | 'footer' | 'unknown' + + +def _tile(known: List[Span], total: int) -> List[Span]: + """Sort *known* spans and fill every gap with an ``unknown`` span, so the + result is a contiguous, non-overlapping tiling of ``[0, total)``. Overlaps + are resolved by clamping to the running position (first writer wins).""" + out: List[Span] = [] + pos = 0 + for s in sorted(known, key=lambda x: (x.start, x.end)): + if s.end <= pos: + continue # fully behind — dropped overlap + start = max(s.start, pos) + if start > pos: + out.append(Span(pos, start, "unknown", "unknown")) + out.append(s if start == s.start else Span(start, s.end, s.label, s.kind)) + pos = s.end + if pos < total: + out.append(Span(pos, total, "unknown", "unknown")) + return out + + +def annotate_blastware_binary(raw: bytes) -> List[Span]: + """Annotate a Series-3 waveform binary into a gap-free list of spans.""" + total = len(raw) + strt_pos = raw.find(b"STRT") + if strt_pos < 0: + return [Span(0, total, "unrecognized — no STRT record", "unknown")] + + known: List[Span] = [] + if strt_pos > 0: + known.append(Span(0, strt_pos, "File header", "header")) + known.append(Span(strt_pos, strt_pos + _STRT_LEN, "STRT record", "strt")) + + body_start = strt_pos + _STRT_LEN + footer_start = total - _FOOTER_LEN + if footer_start >= body_start: + known.append(Span(footer_start, total, "File footer", "footer")) + else: + footer_start = total # file too short for a footer + + body = raw[body_start:footer_start] + for rec in walk_records(body): + hi, lo = rec["mode"] + label = f"{rec['channel']} record (seg {rec['segment_index']}, mode {hi:02x} {lo:02x})" + known.append(Span(body_start + rec["offset"], body_start + rec["end"], label, "sample")) + + return _tile(known, total) diff --git a/minimateplus/event_file_io.py b/minimateplus/event_file_io.py index d23fc6d..188df85 100644 --- a/minimateplus/event_file_io.py +++ b/minimateplus/event_file_io.py @@ -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.29.0" +TOOL_VERSION = "0.31.0" # +/sensor_check group (schema v2); gates the backfill regen try: # Best-effort: prefer the installed metadata when it's NEWER than the @@ -960,6 +960,11 @@ def read_blastware_file(path: Union[str, Path]) -> Event: project=project, client=client, operator=user, sensor_location=seisloc, ) ev.raw_samples = samples + # Sensor self-check traces from the binary's trailing block (waveform + # events only; returns {} for histograms / when absent). Carried on the + # Event so the .h5 writer persists them device-agnostically. + from minimateplus.sensor_check import decode_sensor_check + ev.sensor_check = decode_sensor_check(raw) or None # Only compute peaks from samples when we actually have samples. # For events the codec couldn't decode (histogram-mode bodies, until # the §7.6.2 histogram codec is wired in), samples is an empty dict diff --git a/minimateplus/models.py b/minimateplus/models.py index 48fd326..9e7865d 100644 --- a/minimateplus/models.py +++ b/minimateplus/models.py @@ -544,6 +544,15 @@ class Event: pretrig_samples: Optional[int] = None # from STRT record: pre-trigger sample count rectime_seconds: Optional[int] = None # from STRT record: record duration (seconds) + # Sensor self-check traces keyed by channel label — the short diagnostic + # waveforms the unit records when it pulses each sensor before monitoring + # (geophone ring-downs + a mic pulse train). Decoded from the binary by + # the per-series decoder (minimateplus.sensor_check / micromate.sensor_check) + # and carried here so the .h5 writer can persist them device-agnostically. + # Raw ADC counts; the source series' scale differs but the trace is a + # shape diagnostic (rendered fit-to-box). None when absent. + sensor_check: Optional[dict] = None # {"Tran": [...], ..., "MicL": [...]} + # ── Debug / introspection ───────────────────────────────────────────────── # Raw 210-byte waveform record bytes, set when debug mode is active. # Exposed by the SFM server via ?debug=true so field layouts can be verified. diff --git a/minimateplus/sensor_check.py b/minimateplus/sensor_check.py new file mode 100644 index 0000000..7432ea8 --- /dev/null +++ b/minimateplus/sensor_check.py @@ -0,0 +1,146 @@ +r"""Decode the Blastware sensor self-check waveforms from a series-3 event binary. + +Reverse-engineered 2026-09-15 against 7 BE12844 (MiniMate Plus) oracle events. +After the main waveform record-chain and the trailing metadata / per-channel +calibration records, the binary carries four length-prefixed records tagged +0x3c-0x3f: the sensor self-check traces the unit records when it pulses each +sensor before monitoring. Blastware draws these as the little waveforms in the +"Sensor Check" strip on the right of the Event Report. + + * 0x3c / 0x3d / 0x3e = Tran / Vert / Long geophone ring-downs (a damped + oscillation at the geophone's resonance, ~7-8 Hz at 1024 sps). + * 0x3f = MicL, a pulse train at the mic self-test frequency + (~20 Hz), whose zero-crossing frequency is BW's mic "Channel Test" freq. + +Record framing (per record, all four chained by their length prefix):: + + [len:2 BE][id:1][00 00][Nchan:1][12-byte header][delta stream][40 02][6B] + \_________________ payload (len bytes) _______________________________/ + +The delta stream is ``payload[20 : len-8]`` (the ``40 02`` terminator sits at +``len-8``, followed by 6 trailing bytes). It uses the exact same 10/20/30/00 +delta-block tags as the main waveform codec +(:mod:`minimateplus.waveform_codec`), decoded here from an implicit anchor of 0 +— so the traces come out in the same 16-count raw units as the main waveform +(LSB = 0.005 in/s at Normal range for the geophones). +""" +from __future__ import annotations + +from typing import Dict, List + +from minimateplus.waveform_codec import walk_body + +# Record id → channel. Order mirrors the trailing per-channel calibration +# records (Tran / Vert / Long / MicL), confirmed against BW's sensor-check +# frequencies on all 7 oracle events. +_ID_TO_CHANNEL = {0x3C: "Tran", 0x3D: "Vert", 0x3E: "Long", 0x3F: "MicL"} +_CHAIN_IDS = (0x3C, 0x3D, 0x3E, 0x3F) + +_HEADER_LEN = 20 # payload bytes before the delta stream +_TRAILER_LEN = 8 # 40 02 terminator + 6 trailing bytes after the stream + + +def _s4(nib: int) -> int: + """Sign-extend a 4-bit nibble delta.""" + return nib - 16 if nib >= 8 else nib + + +def _i8(byte: int) -> int: + """Sign-extend an 8-bit int delta.""" + return byte - 256 if byte >= 128 else byte + + +def _decode_delta_stream(buf: bytes) -> List[int]: + """Accumulate a 10/20/30/00 delta-block stream from an anchor of 0, + stopping at the 0x40 terminator. + + Mirrors the block semantics in + :func:`minimateplus.waveform_codec.decode_waveform_v2` (fully decoded & + byte-exact as of 2026-05-11); see that module for the format details. + """ + out: List[int] = [] + cur = 0 + for blk in walk_body(buf, 0): + fam = blk.tag_hi & 0xF0 + if fam == 0x10: + # nibble deltas, high nibble first + for byte in blk.data: + for nib in ((byte >> 4) & 0xF, byte & 0xF): + cur += _s4(nib) + out.append(cur) + elif fam == 0x20: + # int8 deltas + for byte in blk.data: + cur += _i8(byte) + out.append(cur) + elif fam == 0x30: + # 12-bit signed deltas, packed as tag_lo/4 groups of 6 bytes + for g in range(blk.tag_lo // 4): + grp = blk.data[g * 6:(g + 1) * 6] + if len(grp) < 6: + break + high_word = (grp[0] << 8) | grp[1] + for k in range(4): + nib = (high_word >> (12 - 4 * k)) & 0xF + v = (nib << 8) | grp[2 + k] + if v >= 0x800: + v -= 0x1000 + cur += v + out.append(cur) + elif fam == 0x00: + # RLE zero-delta run (wide form carries the high nibble in the tag) + run = ((blk.tag_hi & 0x0F) << 8) | blk.tag_lo + out.extend([cur] * run) + elif fam == 0x40: + # segment / record terminator + break + return out + + +def _find_chain(body: bytes): + """Locate the four length-prefixed sensor-check records. + + Returns a list of ``(offset, id, length)`` or ``None``. The chain is + validated by walking the ids 0x3c → 0x3d → 0x3e → 0x3f via their own length + prefixes, so a stray 0x3c byte in the waveform data cannot match. + """ + for p in range(len(body) - 6): + if body[p + 2] == 0x3C and body[p + 3] == 0 and body[p + 4] == 0: + q = p + recs = [] + ok = True + for expect in _CHAIN_IDS: + if q + 3 > len(body) or body[q + 2] != expect: + ok = False + break + length = int.from_bytes(body[q:q + 2], "big") + recs.append((q, expect, length)) + q = q + 2 + length + if ok and len(recs) == 4: + return recs + return None + + +def decode_sensor_check(raw: bytes) -> Dict[str, List[int]]: + """Decode the four sensor self-check traces from a series-3 event binary. + + Returns ``{"Tran": [...], "Vert": [...], "Long": [...], "MicL": [...]}`` in + raw decode units (same 16-count LSB as the main waveform), or ``{}`` if the + binary carries no sensor-check block (a histogram event, a non-series-3 + file, or a unit/firmware that doesn't store it). + """ + strt = raw.find(b"STRT") + if strt < 0 or len(raw) < strt + 21 + 26: + return {} + body = raw[strt + 21: len(raw) - 26] + chain = _find_chain(body) + if not chain: + return {} + out: Dict[str, List[int]] = {} + for off, rid, length in chain: + payload = body[off + 2: off + 2 + length] + if len(payload) < _HEADER_LEN + _TRAILER_LEN: + continue + stream = payload[_HEADER_LEN: length - _TRAILER_LEN] + out[_ID_TO_CHANNEL[rid]] = _decode_delta_stream(stream) + return out diff --git a/minimateplus/waveform_codec.py b/minimateplus/waveform_codec.py index 9a74780..e8d04dc 100644 --- a/minimateplus/waveform_codec.py +++ b/minimateplus/waveform_codec.py @@ -722,7 +722,18 @@ STREAM_END_ID = 0x06 MODE_DELTA = (0x02, 0x00) MODE_ABSOLUTE = (0x01, 0x00) MODE_RAW12 = (0x00, 0x03) -_MODES = (MODE_DELTA, MODE_ABSOLUTE, MODE_RAW12) +# 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) def _u16(b: bytes, p: int) -> int: @@ -747,7 +758,18 @@ 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 - return (None, None) if (nn == 0 or nn > 0x08) else (2 * nn + 2, nn) + # 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 if nn == 0 or nn % 4: return None, None if hi == 0x00: @@ -761,6 +783,11 @@ 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] = [] @@ -785,13 +812,17 @@ 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 ``00 00 03`` preamble that data is - raw 12-bit with no tags at all and cannot be block-walked — scan instead. + 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. """ - if len(body) >= 3 and (body[1], body[2]) == MODE_RAW12: + if len(body) >= 3 and (body[1], body[2]) in _UNTAGGED_MODES: scan_from = 3 else: - i = 7 + # 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 while i < len(body): if is_record(body, i): nxt = i + 2 + _u16(body, i + 2) @@ -850,7 +881,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_RAW12): + if preamble not in (MODE_DELTA, MODE_ABSOLUTE, MODE_RAW12, MODE_RAW16): return None first = find_first_record(body) if first is None: @@ -895,6 +926,10 @@ 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])) @@ -908,4 +943,6 @@ 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 diff --git a/pyproject.toml b/pyproject.toml index cd9281b..10ea7e7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "seismo-relay" -version = "0.29.0" +version = "0.31.0" description = "Python client and REST server for MiniMate Plus seismographs" requires-python = ">=3.10" dependencies = [ diff --git a/scratch/verify_thor_against_csv.py b/scratch/verify_thor_against_csv.py new file mode 100644 index 0000000..4b8c40c --- /dev/null +++ b/scratch/verify_thor_against_csv.py @@ -0,0 +1,228 @@ +#!/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: + + /UM13981_20220207084555.IDFW + /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()) diff --git a/scripts/backfill_thor_events.py b/scripts/backfill_thor_events.py index 41e7935..581b729 100644 --- a/scripts/backfill_thor_events.py +++ b/scripts/backfill_thor_events.py @@ -305,6 +305,11 @@ def main(argv=None) -> int: default=0, ) ev.total_samples = ev.total_samples or n_samp + # Sensor self-check traces from the IDFW fixed + # header, so regenerated .h5 files gain the v2 + # /sensor_check group (mirrors save_imported_idf). + from micromate.sensor_check import decode_idf_sensor_check + ev.sensor_check = decode_idf_sensor_check(binary_bytes) or None event_hdf5.write_event_hdf5( hdf5_path, ev, diff --git a/seismo_lab.py b/seismo_lab.py index 1986127..5cad795 100644 --- a/seismo_lab.py +++ b/seismo_lab.py @@ -54,6 +54,7 @@ from s3_analyzer import ( # noqa: E402 write_claude_export, ) from frame_db import FrameDB # noqa: E402 +from minimateplus.binary_annotate import annotate_blastware_binary # noqa: E402 # ── colour palette ──────────────────────────────────────────────────────────── BG = "#1e1e1e" @@ -2675,6 +2676,95 @@ class DownloadPanel(tk.Frame): self._on_capture_ready(bw_path, s3_path, label) +# ───────────────────────────────────────────────────────────────────────────── +# Inspector panel — annotated hex view of a Series-3 binary +# ───────────────────────────────────────────────────────────────────────────── + +class InspectorPanel(tk.Frame): + """Load any Series-3 waveform binary and read it as an annotated hex dump. + + Regions the decoder understands (header, STRT, per-channel sample records, + footer) are labelled and colour-coded; everything the decoder cannot account + for is flagged UNKNOWN, so undecoded bytes stand out for hand-inspection. + """ + + _KIND_COLOR = { + "header": ACCENT, + "strt": YELLOW, + "sample": COL_S3, + "footer": FG_DIM, + "unknown": RED, + } + + def __init__(self, parent: tk.Widget, initialdir=None, **kw) -> None: + super().__init__(parent, bg=BG, **kw) + self._path = None + self._initialdir = initialdir + self._build() + + def _build(self) -> None: + bar = tk.Frame(self, bg=BG2) + bar.pack(side=tk.TOP, fill=tk.X) + tk.Button(bar, text="Open binary…", command=self._open, bg=BG3, fg=FG, + relief=tk.FLAT, font=MONO, activebackground=ACCENT).pack(side=tk.LEFT, padx=6, pady=6) + self._path_var = tk.StringVar(value="(no file loaded)") + tk.Label(bar, textvariable=self._path_var, bg=BG2, fg=FG_DIM, font=MONO).pack(side=tk.LEFT, padx=6) + self._summary_var = tk.StringVar(value="") + tk.Label(bar, textvariable=self._summary_var, bg=BG2, fg=FG, font=MONO).pack(side=tk.RIGHT, padx=10) + + legend = tk.Frame(self, bg=BG2) + legend.pack(side=tk.TOP, fill=tk.X) + tk.Label(legend, text="legend:", bg=BG2, fg=FG_DIM, font=MONO).pack(side=tk.LEFT, padx=(8, 2)) + for kind, color in self._KIND_COLOR.items(): + tk.Label(legend, text=f"■ {kind}", bg=BG2, fg=color, font=MONO).pack(side=tk.LEFT, padx=5, pady=2) + + self._text = scrolledtext.ScrolledText( + self, bg=BG, fg=FG, insertbackground=FG, font=MONO, wrap=tk.NONE, borderwidth=0) + self._text.pack(side=tk.TOP, fill=tk.BOTH, expand=True) + for kind, color in self._KIND_COLOR.items(): + self._text.tag_configure(kind, foreground=color) + self._text.tag_configure("label", foreground="#ffffff", font=("Consolas", 9, "bold")) + self._text.tag_configure("dim", foreground=FG_DIM) + self._text.configure(state=tk.DISABLED) + + def _open(self) -> None: + p = filedialog.askopenfilename(title="Open a Series-3 binary", initialdir=self._initialdir) + if p: + self.load(Path(p)) + + def load(self, path: Path) -> None: + try: + raw = path.read_bytes() + spans = annotate_blastware_binary(raw) + except Exception as e: # noqa: BLE001 — surface any read/annotate failure to the user + messagebox.showerror("Inspector", f"Failed to read/annotate:\n{path}\n\n{e}") + return + self._path = path + self._path_var.set(str(path)) + self._render(raw, spans) + + def _render(self, raw: bytes, spans) -> None: + t = self._text + t.configure(state=tk.NORMAL) + t.delete("1.0", tk.END) + unknown = sum(s.end - s.start for s in spans if s.kind == "unknown") + pct = 100 * unknown / max(1, len(raw)) + self._summary_var.set(f"{len(raw)} B · {len(spans)} regions · {pct:.1f}% unknown") + for s in spans: + t.insert(tk.END, f"\n── {s.label} [0x{s.start:04x}:0x{s.end:04x}] {s.end - s.start} B ──\n", ("label",)) + self._insert_hex(t, raw, s.start, s.end, s.kind) + t.configure(state=tk.DISABLED) + + def _insert_hex(self, t: tk.Text, raw: bytes, start: int, end: int, kind: str) -> None: + for off in range(start, end, 16): + row = raw[off:min(off + 16, end)] + hx = " ".join(f"{b:02x}" for b in row).ljust(16 * 3 - 1) + txt = "".join(chr(b) if 32 <= b < 127 else "." for b in row) + t.insert(tk.END, f" 0x{off:04x} ", ("dim",)) + t.insert(tk.END, hx, (kind,)) + t.insert(tk.END, f" {txt}\n", ("dim",)) + + # ───────────────────────────────────────────────────────────────────────────── # Main application window # ───────────────────────────────────────────────────────────────────────────── @@ -2730,6 +2820,9 @@ class SeismoLab(tk.Tk): ) nb.add(self._download_panel, text=" Download ") + self._inspector_panel = InspectorPanel(nb) + nb.add(self._inspector_panel, text=" Inspector ") + self._nb = nb self.protocol("WM_DELETE_WINDOW", self._on_close) diff --git a/sfm/compliance.py b/sfm/compliance.py new file mode 100644 index 0000000..bdc4692 --- /dev/null +++ b/sfm/compliance.py @@ -0,0 +1,133 @@ +"""USBM RI8507 / OSMRE blasting compliance chart. + +Renders the velocity-vs-frequency compliance scatter Blastware draws on its Event +Report: each channel's significant waveform cycles as ``(frequency, peak +velocity)`` points on log-log axes against the regulatory limit curve(s). A point +below the curve passes; above fails. + +Two pieces, kept separate so both can be reused/extended: + * ``limit_at`` / ``limit_curve`` — the regulatory limit curve(s), as data. + * ``channel_compliance_points`` — the per-cycle (freq, velocity) scatter, by + the zero-crossing method (matches Blastware: each channel's cloud tops out + at that channel's PPV). + +Limit curves (USBM RI8507 Figure B-1 / OSM 30 CFR 816.67), drawn CONTINUOUS — a +constant-displacement bound (sloped, ``v = 2πf·d``) meets a constant-velocity +plateau at the frequency where they're equal, so there are no vertical steps +(matching how Blastware draws it). Two lines: + * **Drywall** (modern gypsum board) — 0.75 in/s plateau (solid). + * **Plaster** on wood lath (older homes) — 0.50 in/s plateau (dashed). +Both use a 0.030 in low-frequency displacement bound and rise through a 0.010 in +displacement bound to a 2.0 in/s high-frequency plateau. Values from USBM RI8507 +(Appendix B) / 30 CFR 816.67; ⚠ confirm the exact shape against a Blastware +report before trusting for compliance. +""" +from __future__ import annotations + +import math +from typing import Dict, Sequence, Tuple + +import numpy as np +from matplotlib.ticker import FixedLocator, NullLocator + +# curve name → (low-freq "ultimate" displacement in, mid velocity plateau in/s, +# high-freq displacement in, high-freq velocity plateau in/s). +# RI8507 Fig B-1 (p.74): ultimate max displacement 0.030 in (< ~4 Hz), plateau +# 0.75 (Drywall) / 0.50 (plaster), rising diagonal at 0.008 in displacement up to +# a 2.0 in/s plateau reached at ~40 Hz. +_CURVES: Dict[str, Tuple[float, float, float, float]] = { + "Drywall": (0.030, 0.75, 0.008, 2.00), + "Plaster": (0.030, 0.50, 0.008, 2.00), +} +# how each curve is stroked on the chart +_CURVE_STYLE = {"Drywall": {"ls": "-", "lw": 1.0}, "Plaster": {"ls": "--", "lw": 0.9}} + +STANDARDS = tuple(_CURVES) + +# Blastware's channel markers/colours on the compliance chart. +_CHANNEL_STYLE = { + "Tran": ("+", "#d62728"), # red + + "Vert": ("x", "#2ca02c"), # green x + "Long": ("o", "#1f77b4"), # blue o +} + + +def limit_at(freq_hz: float, curve: str = "Drywall") -> float: + """Max allowed PPV (in/s) at ``freq_hz`` for ``curve`` (continuous).""" + d_low, v_mid, d_high, v_high = _CURVES[curve] + f = max(freq_hz, 1.0) + f_a = v_mid / (2.0 * math.pi * d_low) # disp_low → vel_mid + f_b = v_mid / (2.0 * math.pi * d_high) # vel_mid → disp_high + f_c = v_high / (2.0 * math.pi * d_high) # disp_high → vel_high + if f <= f_a: + return 2.0 * math.pi * f * d_low + if f <= f_b: + return v_mid + if f <= f_c: + return 2.0 * math.pi * f * d_high + return v_high + + +def limit_curve(curve: str = "Drywall", fmin: float = 1.0, fmax: float = 100.0, n: int = 400): + """(freqs, limits) sampled across the band for plotting one curve.""" + freqs = np.logspace(np.log10(fmin), np.log10(fmax), n) + return freqs, np.array([limit_at(f, curve) for f in freqs]) + + +def channel_compliance_points( + samples: Sequence[float], sps: float, fmin: float = 1.0, fmax: float = 100.0, + vmin: float = 0.0, +) -> Tuple[np.ndarray, np.ndarray]: + """Per-cycle (frequency, peak velocity) scatter for one channel. + + Zero-crossing method: split the trace at sign changes; each half-cycle + contributes one point at ``(1/(2·half_period), max|amplitude|)``. Matches + Blastware — the cloud's ceiling is the channel PPV. ``samples`` must be in the + velocity unit you want plotted (in/s). Points outside ``[fmin, fmax]`` or at + or below ``vmin`` are dropped. + """ + x = np.asarray(samples, dtype=float) + if x.size < 3: + return np.empty(0), np.empty(0) + zc = np.where(np.diff(np.signbit(x)))[0] + freqs, vels = [], [] + for a, b in zip(zc[:-1], zc[1:]): + half_period = (b - a) / sps + if half_period <= 0: + continue + freqs.append(1.0 / (2.0 * half_period)) + vels.append(float(np.abs(x[a:b + 1]).max())) + f = np.array(freqs) + v = np.array(vels) + keep = (f >= fmin) & (f <= fmax) & (v > vmin) + return f[keep], v[keep] + + +def draw_compliance_chart(ax, channels: Dict[str, Sequence[float]], sps: float) -> None: + """Draw the compliance chart (both limit curves + per-channel scatter).""" + for name, style in _CURVE_STYLE.items(): + cf, cv = limit_curve(name) + ax.plot(cf, cv, color="#333", zorder=3, **style) + + for ch, (marker, color) in _CHANNEL_STYLE.items(): + samples = channels.get(ch) + if samples is None or len(samples) == 0: + continue + f, v = channel_compliance_points(samples, sps) + ax.scatter(f, v, marker=marker, s=12, c=color, linewidths=0.7, zorder=4, label=ch) + + ax.set_xscale("log") + ax.set_yscale("log") + ax.set_xlim(1, 100) + ax.set_ylim(0.0394, 10) + ax.set_box_aspect(1) # square plot box (log-log compliance charts are square) + xt = [1, 2, 5, 10, 20, 50, 100] + yt = [0.0394, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10] + ax.xaxis.set_major_locator(FixedLocator(xt)); ax.xaxis.set_minor_locator(NullLocator()) + ax.yaxis.set_major_locator(FixedLocator(yt)); ax.yaxis.set_minor_locator(NullLocator()) + ax.set_xticklabels([str(v) for v in xt]) + ax.set_yticklabels([("%g" % v) for v in yt]) + ax.set_xlabel("Frequency (Hz)", fontsize=7) + ax.set_ylabel("Velocity (in/s)", fontsize=7) + ax.tick_params(labelsize=6) + ax.grid(True, which="both", ls=":", lw=0.4, color="#ccc") diff --git a/sfm/event_hdf5.py b/sfm/event_hdf5.py index a25b34d..c63d993 100644 --- a/sfm/event_hdf5.py +++ b/sfm/event_hdf5.py @@ -12,8 +12,11 @@ Layout written to `.h5`: ├─ samples_int16/ (optional) │ ├─ Tran (int16, raw ADC counts) shape: (N,) │ └─ ... per channel (only when present in the source) + ├─ sensor_check/ (optional, schema v2+) + │ ├─ Tran (int32, raw counts) shape: (M,) M ≪ N + │ └─ ... per channel present in the source (MicL absent on 3-channel units) └─ root attrs (event metadata): - schema_version int = 1 + schema_version int = 2 kind str = "sfm.event.hdf5" serial str waveform_key str (8-hex) @@ -64,7 +67,7 @@ from minimateplus.models import Event log = logging.getLogger(__name__) -SCHEMA_VERSION = 1 +SCHEMA_VERSION = 2 # v2 adds the optional /sensor_check group HDF5_KIND = "sfm.event.hdf5" # Geophone full-scale velocity per range (in/s). Confirmed in CLAUDE.md @@ -270,6 +273,22 @@ def write_event_hdf5( ) igrp.attrs["mic_psi_per_count"] = float(mic_factor) + # /sensor_check — optional short diagnostic self-check traces (schema + # v2+). Raw ADC counts (a shape diagnostic; the per-series count scale + # differs, and the renderer fits each trace to its box). Only channels + # the decoder found are written — 3-channel units carry no MicL. + sc = event.sensor_check or {} + if sc: + scgrp = f.create_group("sensor_check") + for ch in ("Tran", "Vert", "Long", "MicL"): + vals = sc.get(ch) + if vals: + scgrp.create_dataset( + ch, data=np.asarray(vals, dtype=np.int32), + compression="gzip", compression_opts=4, shuffle=True, + ) + scgrp.attrs["units"] = "raw_counts" + import os os.replace(tmp, path) @@ -334,6 +353,16 @@ def read_event_hdf5(path: Union[str, Path]) -> dict: if mic_attr is not None: mic_psi = float(mic_attr) + # /sensor_check — optional (schema v2+); absent on older files. + sensor_check = None + scgrp = f.get("sensor_check") + if scgrp is not None: + sensor_check = {} + for ch in ("Tran", "Vert", "Long", "MicL"): + ds = scgrp.get(ch) + if ds is not None: + sensor_check[ch] = np.asarray(ds[()]) + return { "schema_version": sv, "kind": attrs.get("kind"), @@ -341,6 +370,7 @@ def read_event_hdf5(path: Union[str, Path]) -> dict: "samples": samples, "samples_int16": samples_int16, "mic_psi_per_count": mic_psi, + "sensor_check": sensor_check, } @@ -431,11 +461,16 @@ def plot_json_from_hdf5( event_id: Optional[str] = None, index: Optional[int] = None, ) -> dict: - """Build a `sfm.plot.v1` JSON dict from a stored .h5 file.""" + """Build a `sfm.plot.v1` JSON dict from a stored .h5 file. + + The dict also carries a top-level ``sensor_check`` key (the raw self-check + traces as ``{ch: [int]}``, or None) beyond the plot schema, so report + generation can read the traces from the same single .h5 load. + """ data = read_event_hdf5(path) a = data["attrs"] s = data["samples"] - return _build_plot_dict( + out = _build_plot_dict( n_samples=len(s["Tran"]) if "Tran" in s else 0, sample_rate=int(a.get("sample_rate", 1024) or 1024), pretrig_samples=int(a.get("pretrig_samples", 0) or 0), @@ -463,6 +498,11 @@ def plot_json_from_hdf5( event_id=event_id, index=index, ) + scd = data.get("sensor_check") + out["sensor_check"] = ( + {ch: v.tolist() for ch, v in scd.items()} if scd else None + ) + return out def _build_plot_dict( diff --git a/sfm/report_pdf.py b/sfm/report_pdf.py index 25859d1..d4d60e7 100644 --- a/sfm/report_pdf.py +++ b/sfm/report_pdf.py @@ -121,6 +121,13 @@ class ReportData: t0_ms: Optional[float] = None dt_ms: Optional[float] = None + # Sensor self-check traces — {ch: [samples]} in raw counts, read from the + # standardized .h5 (/sensor_check group, schema v2+) where the per-series + # decoder stored them at ingest. The little diagnostic waveforms BW draws + # in its "Sensor Check" strip. Empty when absent (pre-v2 .h5, histogram, + # or 3-channel unit's MicL). + sensor_check_waveforms: dict = field(default_factory=dict) + # Record-type discriminator record_type: Optional[str] = None is_histogram: bool = False @@ -246,6 +253,8 @@ def gather_report_data( "peak_accel_g": ch.get("peak_accel_g"), "peak_disp_in": ch.get("peak_disp_in"), "sensor_check": sc_ch.get("result"), + "sc_freq_hz": sc_ch.get("freq_hz"), + "sc_ratio": sc_ch.get("ratio"), "peak_date": peak_date, "peak_time": peak_time, }) @@ -287,6 +296,12 @@ def gather_report_data( rd.pretrig_samples = ta.get("pretrig_samples") rd.t0_ms = ta.get("t0_ms") rd.dt_ms = ta.get("dt_ms") + # Sensor self-check traces — read from the standardized .h5 (schema + # v2+). Device-agnostic: whichever decoder produced the event + # stored them at ingest, so SFM reads them here without knowing or + # caring about the source instrument series. Empty on pre-v2 files + # (until backfilled) and on 3-channel / histogram events. + rd.sensor_check_waveforms = wf.get("sensor_check") or {} except Exception as exc: log.warning("gather_report_data: hdf5 read failed: %s", exc) @@ -396,9 +411,34 @@ def _render_waveform_layout(fig, rd: ReportData) -> None: ax_stats = fig.add_subplot(gs[2]); ax_stats.axis("off") _draw_channel_stats_waveform(ax_stats, rd) + _draw_compliance_panel(fig, rd) _draw_waveform_subplot(fig, gs[3], rd) +# Compliance-chart placement, in figure fractions. Measured directly off a +# Blastware Event Report PDF (ref-stuff/n844lqhbzt0w_bw_pdf.pdf) so the chart +# matches BW's size and position: it spans from just under the header down +# through the stats band, hard against the right page margin. The left edge +# leaves room for the y-axis tick labels + "Velocity (in/s)" title, which the +# compacted stats table (see _draw_channel_stats_waveform) is sized to clear. +_COMPLIANCE_BOX = (0.489, 0.502, 0.951, 0.867) # x0, y0, x1, y1 + + +def _draw_compliance_panel(fig, rd: ReportData) -> None: + """Large USBM RI8507 compliance chart in the upper-right, sized and + positioned to match Blastware's Event Report (see _COMPLIANCE_BOX).""" + x0, y0, x1, y1 = _COMPLIANCE_BOX + fig.text((x0 + x1) / 2, y1 + 0.006, "USBM RI8507 And OSMRE", fontsize=9, + weight="bold", color="#333", ha="center", va="bottom") + if rd.channels and rd.sample_rate_sps: + from sfm.compliance import draw_compliance_chart + ax = fig.add_axes([x0, y0, x1 - x0, y1 - y0]) + draw_compliance_chart(ax, rd.channels, rd.sample_rate_sps) + else: + fig.text((x0 + x1) / 2, (y0 + y1) / 2, "(no waveform data)", fontsize=8, + color="#bbb", ha="center", va="center", style="italic") + + def _render_histogram_layout(fig, rd: ReportData) -> None: """Histogram layout: header / mic-only / per-channel stats / bar plot. @@ -477,11 +517,11 @@ def _split_iso_to_date_time(iso: Optional[str]) -> tuple[Optional[str], Optional return (None, None) -def _kv(ax, x, y, label, value, *, label_w=0.18): +def _kv(ax, x, y, label, value, *, label_w=0.18, fontsize=8): """Render a 'Label Value' row at axes-coordinates (x, y).""" - ax.text(x, y, label, fontsize=8, color="#555", ha="left", va="top", + ax.text(x, y, label, fontsize=fontsize, color="#555", ha="left", va="top", transform=ax.transAxes) - ax.text(x + label_w, y, _fmt(value), fontsize=8, ha="left", va="top", + ax.text(x + label_w, y, _fmt(value), fontsize=fontsize, ha="left", va="top", transform=ax.transAxes, family="monospace") @@ -544,14 +584,17 @@ def _draw_header_columns(ax, rows_left, rd: ReportData) -> None: ("File Name", rd.file_name), ("Post Event Notes", rd.post_event_notes), ] + # fontsize 7.5 (BW's header is a touch smaller than our body text) + a + # tighter right-column value indent so the long serial+firmware line + # ("BE##### V ##.##-#.## MiniMate Plus") fits without running off the page. y = 0.95 dy = 0.095 for label, value in rows_left: - _kv(ax, 0.0, y, label, value, label_w=0.18) + _kv(ax, 0.0, y, label, value, label_w=0.18, fontsize=7.5) y -= dy y = 0.95 for label, value in rows_right: - _kv(ax, 0.55, y, label, value, label_w=0.20) + _kv(ax, 0.55, y, label, value, label_w=0.14, fontsize=7.5) y -= dy @@ -574,19 +617,14 @@ def _draw_mic_and_usbm(ax, rd: ReportData) -> None: transform=ax.transAxes, va="top") rows = _mic_rows(rd) y = 0.80 + # Tighter label indent + slightly smaller font so the long "Channel Test + # Passed (Freq = … Amp = … mv)" line clears the enlarged compliance chart's + # left edge (_COMPLIANCE_BOX) instead of running behind it. for label, value in rows: - _kv(ax, 0.0, y, label, value, label_w=0.18) + _kv(ax, 0.0, y, label, value, label_w=0.13, fontsize=7) y -= 0.15 - - # USBM chart placeholder — upper-right. Real piecewise compliance - # curves are a separate work item; for now this just shows the title - # + a "see report" message so the layout is correct. - ax.text(0.72, 0.97, "USBM RI8507 And OSMRE", - fontsize=9, weight="bold", color="#333", ha="center", va="top", - transform=ax.transAxes) - ax.text(0.72, 0.50, "[compliance chart\ncoming soon]", - fontsize=8, color="#bbb", ha="center", va="center", - transform=ax.transAxes, style="italic") + # The USBM compliance chart is drawn as its own large square panel spanning + # the mic + stats rows on the right — see _draw_compliance_panel(). def _mic_rows(rd: ReportData) -> list[tuple[str, Optional[str]]]: @@ -636,8 +674,18 @@ def _draw_channel_stats_waveform(ax, rd: ReportData) -> None: ("Peak Acceleration", "peak_accel_g", "g"), ("Peak Displacement", "peak_disp_in", "in"), ("Sensor Check", "sensor_check", ""), + # Sensor-check sub-rows (indented under "Sensor Check", like BW): the + # geophone ring-down frequency + overswing ratio from the self-check. + (" Frequency", "sc_freq_hz", "Hz"), + (" Overswing Ratio", "sc_ratio", ""), ] - _draw_stats_table(ax, rd, rows_spec) + # Compacted to the left half so the enlarged compliance chart (BW-sized, + # right against the page margin) has room — see _COMPLIANCE_BOX. + _draw_stats_table( + ax, rd, rows_spec, + bbox_width=0.42, fontsize=7.5, + col_widths=[0.185, 0.065, 0.065, 0.065, 0.040], + ) _draw_pvs_summary(ax, rd, n_data_rows=len(rows_spec)) @@ -698,19 +746,39 @@ def _draw_pvs_summary( table_bottom_y = getattr(ax, "_stats_table_bottom", -0.10) pvs_y = table_bottom_y - 0.04 # small gap below the table border - # Centered for visual balance — looks intentional rather than offset. - # The original BW-replica had a "NA: Not Applicable" caption below - # this line; dropped because we use "—" for missing values and the - # legend was always squished against the PVS line. - ax.text(0.5, pvs_y, line, fontsize=9, weight="bold", - ha="center", va="top", transform=ax.transAxes) + # Centered under the stats table for visual balance — looks intentional + # rather than offset. When the table is compacted (waveform layout), it + # occupies only the left portion of the axes, so center on the table's + # width rather than the full axes (which would push the line under the + # compliance chart). The original BW-replica had a "NA: Not Applicable" + # caption below this line; dropped because we use "—" for missing values. + table_w = getattr(ax, "_stats_table_width", 0.80) + if table_w < 0.79: + # Compacted (waveform) layout: left-align under the table, one point + # smaller, so the line clears the enlarged compliance chart's + # bottom-left tick labels on the right. + ax.text(0.0, pvs_y, line, fontsize=8, weight="bold", + ha="left", va="top", transform=ax.transAxes) + else: + ax.text(0.5, pvs_y, line, fontsize=9, weight="bold", + ha="center", va="top", transform=ax.transAxes) -def _draw_stats_table(ax, rd: ReportData, rows_spec: list[tuple[str, str, str]]) -> None: +def _draw_stats_table( + ax, rd: ReportData, rows_spec: list[tuple[str, str, str]], + *, bbox_width: float = 0.80, fontsize: float = 8, + col_widths: Optional[list[float]] = None, +) -> None: """Render a per-channel stats table (Tran/Vert/Long). rows_spec: list of (label, field_name_in_channel_stats, unit_string) + + ``bbox_width`` / ``col_widths`` / ``fontsize`` let a caller compact the + table (the waveform layout packs it into the left half to clear the + compliance chart; the histogram layout keeps the wider defaults). """ + if col_widths is None: + col_widths = [0.28, 0.14, 0.14, 0.14, 0.10] headers = ["", "Tran", "Vert", "Long", ""] ch_lookup = {c["name"]: c for c in rd.channel_stats} @@ -726,6 +794,8 @@ def _draw_stats_table(ax, rd: ReportData, rows_spec: list[tuple[str, str, str]]) if field == "zc_freq_hz": prefix = ">" if ch_rec.get("zc_freq_above_range") else "" return f"{prefix}{val:.0f}" + if field in ("sc_freq_hz", "sc_ratio"): + return f"{val:.1f}" # BW shows 1 decimal (7.5 Hz, 3.6) return f"{val:.3f}" return str(val) @@ -750,16 +820,17 @@ def _draw_stats_table(ax, rd: ReportData, rows_spec: list[tuple[str, str, str]]) table_bottom = 1.0 - table_height tbl = ax.table( cellText=table_data, - colWidths=[0.28, 0.14, 0.14, 0.14, 0.10], + colWidths=col_widths, cellLoc="left", edges="open", - bbox=[0.0, table_bottom, 0.80, table_height], + bbox=[0.0, table_bottom, bbox_width, table_height], ) tbl.auto_set_font_size(False) - tbl.set_fontsize(8) + tbl.set_fontsize(fontsize) for j in range(5): tbl[(0, j)].set_text_props(weight="bold", color="#555") - # Stash the bottom Y so _draw_pvs_summary can position itself below. + # Stash the bottom Y + width so _draw_pvs_summary can position itself. ax._stats_table_bottom = table_bottom + ax._stats_table_width = bbox_width def _channel_axis_color(ch: str) -> str: @@ -769,27 +840,59 @@ def _channel_axis_color(ch: str) -> str: def _draw_waveform_subplot(fig, gridspec_cell, rd: ReportData) -> None: """4-channel stacked waveform plot — Instantel printout order (MicL on top, Tran on bottom), shared x-axis in SECONDS, trigger - triangle markers at t=0, '0.0' baseline label on right of each.""" - inner = gridspec_cell.subgridspec(4, 1, hspace=0.0) + triangle markers at t=0, '0.0' baseline label on right of each. + + When sensor self-check traces are present (rd.sensor_check_waveforms), a + narrow "Sensor Check" strip of per-channel mini-plots is drawn to the right, + aligned to the lanes — matching Blastware's Event Report. + """ + from matplotlib.ticker import MaxNLocator + order = ["MicL", "Long", "Vert", "Tran"] + has_sc = bool(rd.sensor_check_waveforms) + if has_sc: + # main lanes + a narrow sensor-check strip column, flush against the + # main panel (BW shares the border — no gap), with the "0.0" baseline + # labels moved to the right of the strip. Proportions match BW's + # Event Report (main ~0.75 / strip ~0.10 of the panel width). + inner = gridspec_cell.subgridspec(4, 2, width_ratios=[1.0, 0.13], + wspace=0.0, hspace=0.0) + else: + inner = gridspec_cell.subgridspec(4, 1, hspace=0.0) sr = rd.sample_rate_sps or 1024 # Convert ms-based time axis to seconds for the x-axis 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) + + main_axes = [] + sc_axes = [] last_idx = len(order) - 1 for i, ch in enumerate(order): - ax = fig.add_subplot(inner[i]) + ax = fig.add_subplot(inner[i, 0] if has_sc else inner[i]) + main_axes.append(ax) values = rd.channels.get(ch) or [] times = [t0_s + j * dt_s for j in range(len(values))] if values: color = _channel_axis_color(ch) ax.plot(times, values, color=color, linewidth=0.5) - # Symmetric y-axis for geo; zero-anchored for mic. + # Geo: one shared symmetric scale (honest relative amplitudes). + # Mic: symmetric on its own psi scale (different unit). if ch != "MicL": - amax = max((abs(v) for v in values), default=0.001) - ax.set_ylim(-amax * 1.10, amax * 1.10) + ax.set_ylim(-geo_shared, geo_shared) else: amax = max((abs(v) for v in values), default=0.001) ax.set_ylim(-amax * 1.10, amax * 1.10) @@ -797,9 +900,12 @@ def _draw_waveform_subplot(fig, gridspec_cell, rd: ReportData) -> None: # Channel label on the LEFT (matches BW) ax.set_ylabel(ch, fontsize=8, rotation=0, ha="right", va="center", color=_channel_axis_color(ch), weight="bold", labelpad=14) - # "0.0" on the RIGHT (BW convention) - ax.text(1.005, 0.5, "0.0", transform=ax.transAxes, - fontsize=7, color="#555", va="center", ha="left") + # "0.0" baseline label on the RIGHT (BW convention). With the sensor- + # check strip attached, it goes to the right of the STRIP (drawn below); + # otherwise just outside the main lane. + if not has_sc: + ax.text(1.005, 0.5, "0.0", transform=ax.transAxes, + fontsize=7, color="#555", va="center", ha="left") ax.grid(True, linestyle="--", linewidth=0.3, color="#bbb", alpha=0.6) # Vertical dashed trigger line at t=0 @@ -814,23 +920,53 @@ def _draw_waveform_subplot(fig, gridspec_cell, rd: ReportData) -> None: else: ax.tick_params(axis="x", labelsize=7) ax.tick_params(axis="y", labelsize=6) + # Stacked lanes touch, so the top/bottom y-tick labels of adjacent lanes + # would overprint at the shared boundary. Prune the extreme ticks so + # each boundary shows clean interior ticks (0.5 / 0.0 / -0.5) only. + ax.yaxis.set_major_locator(MaxNLocator(nbins=4, prune="both")) + + # Sensor self-check mini-plot in the right strip (aligned to this lane). + if has_sc: + scx = fig.add_subplot(inner[i, 1]) + sc_axes.append(scx) + sc_vals = rd.sensor_check_waveforms.get(ch) or [] + if sc_vals: + _col = _channel_axis_color(ch) + # Faint zero baseline (BW draws the channel baseline through the + # strip) — reference for the one-sided geophone ring-downs. + scx.axhline(0.0, color=_col, linewidth=0.3, alpha=0.4) + scx.plot(range(len(sc_vals)), sc_vals, color=_col, linewidth=0.5) + # Fit the trace to the box (BW-style) rather than a symmetric + # scale: the geo self-checks are one-sided dips, so a symmetric + # scale would strand them in the bottom half with an empty top. + _lo, _hi = min(sc_vals), max(sc_vals) + _pad = 0.10 * ((_hi - _lo) or 1.0) + scx.set_ylim(_lo - _pad, _hi + _pad) + scx.set_xticks([]); scx.set_yticks([]) + for _s in scx.spines.values(): + _s.set_linewidth(0.4); _s.set_color("#999") + # "0.0" baseline label to the RIGHT of the strip (BW convention) + scx.text(1.10, 0.5, "0.0", transform=scx.transAxes, + fontsize=7, color="#555", va="center", ha="left") # Trigger triangle marker ▼ above the top channel at t=0 - top_ax = fig.axes[-4] # MicL is the first added in this gridspec + top_ax = main_axes[0] # MicL top_ax.plot([0], [top_ax.get_ylim()[1]], marker="v", color="black", markersize=8, clip_on=False, zorder=10) + # "Sensor Check" caption under the strip (BW convention) + if has_sc and sc_axes: + pos = sc_axes[-1].get_position() + fig.text((pos.x0 + pos.x1) / 2, pos.y0 - 0.012, "Sensor Check", + fontsize=7, color="#555", ha="center", va="top") + # Compute scale-per-division for the footer (10 divs across the chart) # 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 - 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 + # 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 "—" 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", diff --git a/sfm/waveform_store.py b/sfm/waveform_store.py index 8373e8d..8c2e26c 100644 --- a/sfm/waveform_store.py +++ b/sfm/waveform_store.py @@ -595,8 +595,19 @@ class WaveformStore: ) # Binary-derived peaks fill in when the .txt didn't supply them. - # They're ~3% low vs the device-authoritative .txt values (residual - # codec drift), so .txt always wins when present. + # + # 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. 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 @@ -651,6 +662,11 @@ class WaveformStore: ev.raw_samples = idf_samples n_samples = max((len(idf_samples.get(ch, [])) for ch in ("Tran", "Vert", "Long", "MicL")), default=0) ev.total_samples = ev.total_samples or n_samples + # Sensor self-check traces from the IDFW fixed header (waveform + # events only; {} on histograms / when absent). Carried on the + # bridged Event so the .h5 writer persists them like series-3. + from micromate.sensor_check import decode_idf_sensor_check + ev.sensor_check = decode_idf_sensor_check(idf_bytes) or None # For IDFH histograms there are no per-sample waveform arrays — the # device stores one peak ADC count per interval per channel. Synthesise diff --git a/tests/fixtures/fft-oracle-2026-09-14/N844LPGH.VV0W b/tests/fixtures/fft-oracle-2026-09-14/N844LPGH.VV0W new file mode 100644 index 0000000..785721f Binary files /dev/null and b/tests/fixtures/fft-oracle-2026-09-14/N844LPGH.VV0W differ diff --git a/tests/fixtures/fft-oracle-2026-09-14/N844LPPR.3S0W b/tests/fixtures/fft-oracle-2026-09-14/N844LPPR.3S0W new file mode 100644 index 0000000..b4395a1 Binary files /dev/null and b/tests/fixtures/fft-oracle-2026-09-14/N844LPPR.3S0W differ diff --git a/tests/fixtures/fft-oracle-2026-09-14/N844LQHB.ZT0W b/tests/fixtures/fft-oracle-2026-09-14/N844LQHB.ZT0W new file mode 100644 index 0000000..6c7afb6 Binary files /dev/null and b/tests/fixtures/fft-oracle-2026-09-14/N844LQHB.ZT0W differ diff --git a/tests/fixtures/fft-oracle-2026-09-14/N844LQUE.T50W b/tests/fixtures/fft-oracle-2026-09-14/N844LQUE.T50W new file mode 100644 index 0000000..fd9a7cd Binary files /dev/null and b/tests/fixtures/fft-oracle-2026-09-14/N844LQUE.T50W differ diff --git a/tests/fixtures/fft-oracle-2026-09-14/N844LR8W.790W b/tests/fixtures/fft-oracle-2026-09-14/N844LR8W.790W new file mode 100644 index 0000000..8726fdb Binary files /dev/null and b/tests/fixtures/fft-oracle-2026-09-14/N844LR8W.790W differ diff --git a/tests/fixtures/fft-oracle-2026-09-14/N844LRCO.G60W b/tests/fixtures/fft-oracle-2026-09-14/N844LRCO.G60W new file mode 100644 index 0000000..9e4fa8f Binary files /dev/null and b/tests/fixtures/fft-oracle-2026-09-14/N844LRCO.G60W differ diff --git a/tests/fixtures/fft-oracle-2026-09-14/N844LRCW.F30W b/tests/fixtures/fft-oracle-2026-09-14/N844LRCW.F30W new file mode 100644 index 0000000..6e9ef75 Binary files /dev/null and b/tests/fixtures/fft-oracle-2026-09-14/N844LRCW.F30W differ diff --git a/tests/fixtures/thor-idf-sc/UM11719_20231219162723.IDFW b/tests/fixtures/thor-idf-sc/UM11719_20231219162723.IDFW new file mode 100644 index 0000000..14620e8 Binary files /dev/null and b/tests/fixtures/thor-idf-sc/UM11719_20231219162723.IDFW differ diff --git a/tests/fixtures/thor-idf-sc/UM12947_20250806134504.IDFW b/tests/fixtures/thor-idf-sc/UM12947_20250806134504.IDFW new file mode 100644 index 0000000..a16b414 Binary files /dev/null and b/tests/fixtures/thor-idf-sc/UM12947_20250806134504.IDFW differ diff --git a/tests/fixtures/thor-idf-sc/UM13981_20220207084555.IDFW b/tests/fixtures/thor-idf-sc/UM13981_20220207084555.IDFW new file mode 100644 index 0000000..4a14035 Binary files /dev/null and b/tests/fixtures/thor-idf-sc/UM13981_20220207084555.IDFW differ diff --git a/tests/fixtures/thor-idf-sc/UM20147_20250531135901.IDFW b/tests/fixtures/thor-idf-sc/UM20147_20250531135901.IDFW new file mode 100644 index 0000000..476ddc9 Binary files /dev/null and b/tests/fixtures/thor-idf-sc/UM20147_20250531135901.IDFW differ diff --git a/tests/test_binary_annotate.py b/tests/test_binary_annotate.py new file mode 100644 index 0000000..3748041 --- /dev/null +++ b/tests/test_binary_annotate.py @@ -0,0 +1,50 @@ +"""Structural annotation of a Series-3 Blastware binary (for the seismo_lab +Binary Inspector). The annotator maps byte ranges to labelled spans; anything +the decoder can't account for is a first-class ``unknown`` span, so the whole +file is tiled and the gaps (candidate FFT/spectral data) are visible. +""" +from pathlib import Path + +from minimateplus.binary_annotate import annotate_blastware_binary, Span + +# A known-good full-3-channel Series-3 waveform binary (the V70 cracking fixture). +FIXTURE = Path(__file__).parent / "fixtures" / "5-11-26" / "M529LL1L.V70" + + +def _raw() -> bytes: + return FIXTURE.read_bytes() + + +def test_spans_tile_the_whole_file(): + raw = _raw() + spans = annotate_blastware_binary(raw) + assert spans, "expected at least one span" + assert spans[0].start == 0 + assert spans[-1].end == len(raw) + for a, b in zip(spans, spans[1:]): + assert a.end == b.start, f"gap/overlap between {a!r} and {b!r}" + for s in spans: + assert s.start < s.end, f"empty/negative span {s!r}" + + +def test_strt_record_is_located(): + raw = _raw() + spans = annotate_blastware_binary(raw) + strt = [s for s in spans if s.kind == "strt"] + assert strt, "expected a STRT region" + assert raw[strt[0].start : strt[0].start + 4] == b"STRT" + + +def test_geo_sample_records_annotated(): + raw = _raw() + spans = annotate_blastware_binary(raw) + chans = {s.label.split()[0] for s in spans if s.kind == "sample"} + # V70 is a full three-geo-channel event. + assert {"Tran", "Vert", "Long"} <= chans, f"expected geo records, got {chans}" + + +def test_footer_is_last(): + raw = _raw() + spans = annotate_blastware_binary(raw) + assert spans[-1].kind == "footer" + assert spans[-1].end - spans[-1].start == 26 diff --git a/tests/test_compliance.py b/tests/test_compliance.py new file mode 100644 index 0000000..c102e9c --- /dev/null +++ b/tests/test_compliance.py @@ -0,0 +1,35 @@ +"""USBM/OSMRE compliance curve + scatter logic (sfm.compliance). +Rendering is verified visually against Blastware reports.""" +import math + +import numpy as np + +from sfm.compliance import limit_at, channel_compliance_points + + +def test_osmre_velocity_segments(): + assert abs(limit_at(6.0) - 0.75) < 1e-9 # 3.5–12 Hz flat + assert abs(limit_at(50.0) - 2.00) < 1e-9 # 30–100 Hz flat + + +def test_displacement_segments(): + assert abs(limit_at(2.0) - 2 * math.pi * 2.0 * 0.030) < 1e-9 # low-freq 0.030 in + assert abs(limit_at(20.0) - 2 * math.pi * 20.0 * 0.008) < 1e-9 # rising diagonal 0.008 in + + +def test_limit_clamps_below_1hz(): + assert limit_at(0.1) == limit_at(1.0) + + +def test_scatter_ceiling_is_ppv_at_dominant_freq(): + # ~27 Hz blast-like trace whose energy peaks mid-record (inside full cycles, + # as a real event does): the scatter cloud's ceiling is the trace PPV and the + # top point sits near the dominant frequency. + sps, n = 1024.0, 3328 + t = np.arange(n) / sps + env = np.exp(-((t - 1.5) ** 2) / (2 * 0.3 ** 2)) + x = 0.9 * env * np.sin(2 * np.pi * 27.0 * t) + f, v = channel_compliance_points(x, sps) + assert len(f) > 20 + assert v.max() >= 0.99 * np.abs(x).max() + assert 20.0 < f[int(np.argmax(v))] < 35.0 diff --git a/tests/test_event_hdf5_sensor_check.py b/tests/test_event_hdf5_sensor_check.py new file mode 100644 index 0000000..fda0dc8 --- /dev/null +++ b/tests/test_event_hdf5_sensor_check.py @@ -0,0 +1,71 @@ +"""The event .h5 carries the sensor self-check traces (schema v2). + +The sensor check is decoded by the per-series decoder and attached to the +standardized Event, so the .h5 writer persists it device-agnostically and SFM +reads it back without knowing which instrument produced it. Old v1 files (no +sensor_check group) must still read cleanly. +""" +import tempfile +from pathlib import Path + +import numpy as np + +from minimateplus.models import Event +from minimateplus.event_file_io import read_blastware_file +from sfm import event_hdf5 + +S3_FIX = Path(__file__).parent / "fixtures" / "fft-oracle-2026-09-14" / "N844LQHB.ZT0W" + + +def _write(ev, **kw): + d = Path(tempfile.mkdtemp()) + p = d / "e.h5" + event_hdf5.write_event_hdf5(p, ev, serial="BE12844", **kw) + return p + + +def test_sensor_check_roundtrips_through_hdf5(): + ev = Event(index=0) + ev.raw_samples = {"Tran": [1, 2, -3], "Vert": [0, 1], "Long": [2], "MicL": [5, -5]} + ev.sample_rate = 1024 + sc = {"Tran": [0, -990, -500, -100], "Vert": [0, -980, -480], + "Long": [0, -986, -470], "MicL": [0, -1800, 1800, -1800]} + ev.sensor_check = sc + + r = event_hdf5.read_event_hdf5(_write(ev)) + assert r["schema_version"] == 2 + assert set(r["sensor_check"]) == {"Tran", "Vert", "Long", "MicL"} + for ch, vals in sc.items(): + assert r["sensor_check"][ch].tolist() == vals + + +def test_plot_json_carries_sensor_check(): + ev = Event(index=0) + ev.raw_samples = {"Tran": [1, 2, 3]} + ev.sample_rate = 1024 + ev.sensor_check = {"Tran": [0, -990, -500], "Vert": [0, -980], + "Long": [0, -986]} # 3-channel: no MicL + pj = event_hdf5.plot_json_from_hdf5(_write(ev)) + assert pj["sensor_check"] is not None + assert "MicL" not in pj["sensor_check"] + assert pj["sensor_check"]["Tran"] == [0, -990, -500] + + +def test_event_without_sensor_check_still_reads_as_v2(): + ev = Event(index=0) + ev.raw_samples = {"Tran": [1, 2, 3]} + ev.sample_rate = 1024 + r = event_hdf5.read_event_hdf5(_write(ev)) + assert r["schema_version"] == 2 + assert r["sensor_check"] is None + assert event_hdf5.plot_json_from_hdf5(_write(ev))["sensor_check"] is None + + +def test_series3_decode_populates_event_sensor_check(): + # The real series-3 decoder attaches the traces to the Event, so the + # ingest/backfill .h5 write picks them up with no extra plumbing. + ev = read_blastware_file(S3_FIX) + assert ev.sensor_check is not None + assert set(ev.sensor_check) == {"Tran", "Vert", "Long", "MicL"} + tran = np.asarray(ev.sensor_check["Tran"], dtype=float) + assert tran.min() < -800 # the geophone ring-down deflection diff --git a/tests/test_idf_binary_codec.py b/tests/test_idf_binary_codec.py new file mode 100644 index 0000000..06e7c1f --- /dev/null +++ b/tests/test_idf_binary_codec.py @@ -0,0 +1,322 @@ +"""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="/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" diff --git a/tests/test_report_pdf_geo_scale.py b/tests/test_report_pdf_geo_scale.py new file mode 100644 index 0000000..6511f03 --- /dev/null +++ b/tests/test_report_pdf_geo_scale.py @@ -0,0 +1,61 @@ +"""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" diff --git a/tests/test_sensor_check.py b/tests/test_sensor_check.py new file mode 100644 index 0000000..b39924a --- /dev/null +++ b/tests/test_sensor_check.py @@ -0,0 +1,67 @@ +"""Blastware sensor self-check waveform decode (minimateplus.sensor_check). + +Reverse-engineered 2026-09-15 against 7 BE12844 (MiniMate Plus) oracle events. +After the main waveform record-chain and the trailing metadata / per-channel +calibration records, a series-3 binary carries four length-prefixed records +tagged 0x3c-0x3f: the sensor self-check traces the unit records when it pulses +each sensor before monitoring (Blastware draws these as the little waveforms in +the "Sensor Check" strip on the right of the Event Report). + + * 0x3c / 0x3d / 0x3e = Tran / Vert / Long geophone ring-downs. + * 0x3f = MicL, a pulse train at the mic self-test frequency. + +The self-check injects a fixed pulse, so the response is near-identical across +events — asserted here as an invariant shape (damped one-sided ring-down for +the geophones, a multi-pulse train for the mic). +""" +from pathlib import Path + +import numpy as np + +from minimateplus.sensor_check import decode_sensor_check + +FIXDIR = Path(__file__).parent / "fixtures" / "fft-oracle-2026-09-14" +EVENTS = sorted(p.name for p in FIXDIR.iterdir()) # 7 BE12844 event binaries + + +def _decode(name): + return decode_sensor_check((FIXDIR / name).read_bytes()) + + +def test_all_four_channels_present(): + for name in EVENTS: + sc = _decode(name) + assert set(sc) == {"Tran", "Vert", "Long", "MicL"}, name + + +def test_geo_channels_are_damped_ringdowns(): + # Each geophone self-check is a large one-sided deflection (~-990 raw) that + # rings back and damps toward a settled value well above the trough. + for name in EVENTS: + sc = _decode(name) + for ch in ("Tran", "Vert", "Long"): + tr = np.asarray(sc[ch], dtype=float) + assert 240 <= len(tr) <= 260, f"{name}:{ch} n={len(tr)}" + assert abs(tr[:3].mean()) < 50, f"{name}:{ch} starts off-baseline" + assert tr.min() < -800, f"{name}:{ch} min {tr.min()}" + assert tr.max() < 60, f"{name}:{ch} unexpected positive swing {tr.max()}" + # damped: settles between the trough and zero, well above the trough + assert tr.min() < tr[-1] < 0, f"{name}:{ch} end {tr[-1]} not between trough and 0" + assert abs(tr[-1]) < 0.6 * abs(tr.min()), f"{name}:{ch} not damped, end {tr[-1]}" + + +def test_mic_channel_is_a_pulse_train(): + for name in EVENTS: + tr = np.asarray(_decode(name)["MicL"], dtype=float) + assert 235 <= len(tr) <= 255, f"{name} mic n={len(tr)}" + # larger dynamic range than the geo ring-down, and swings both ways + assert tr.min() < -1500, f"{name} mic min {tr.min()}" + assert tr.max() > 100, f"{name} mic max {tr.max()}" + # multiple pulses: several deep local minima + deep = (tr[1:-1] < tr[:-2]) & (tr[1:-1] < tr[2:]) & (tr[1:-1] < -800) + assert int(deep.sum()) >= 4, f"{name} mic pulses {int(deep.sum())}" + + +def test_returns_empty_when_no_sensor_check_block(): + assert decode_sensor_check(b"not a blastware file") == {} + assert decode_sensor_check(b"") == {} diff --git a/tests/test_sensor_check_idf.py b/tests/test_sensor_check_idf.py new file mode 100644 index 0000000..00074b2 --- /dev/null +++ b/tests/test_sensor_check_idf.py @@ -0,0 +1,67 @@ +"""Series-4 (Thor / Micromate IDFW) sensor self-check waveform decode. + +Reverse-engineered 2026-09-15 against 4 UM (Thor) oracle events. The IDFW +binary carries the sensor self-check in its fixed-header region (before the +waveform body) as up to four records tagged ``01 0e 3c/3d/3e/3f`` — the SAME +channel ids as series-3 (Tran/Vert/Long/MicL). Unlike series-3's delta-coded +trailing block, series-4 stores each trace as a raw int16-BE array after an +18-byte record header whose sample count is a 2-byte field at offset +8. + +Three-channel (mic-disabled) Thor units carry only 3c/3d/3e — no MicL record. + +Validated by shape (geophone ring-down / mic pulse train) and cross-event +consistency, since there's no Thor Event-Report strip to exact-match against. +""" +from pathlib import Path + +import numpy as np + +from micromate.sensor_check import decode_idf_sensor_check + +FIXDIR = Path(__file__).parent / "fixtures" / "thor-idf-sc" +EVENTS = sorted(p.name for p in FIXDIR.glob("*.IDFW")) + + +def _decode(name): + return decode_idf_sensor_check((FIXDIR / name).read_bytes()) + + +def test_geo_channels_present_and_ringdown_shaped(): + # Every IDFW event has the three geophone self-checks; each is a large + # one-sided deflection (~15000 raw counts) that rings back — the geophone's + # damped impulse response. + for name in EVENTS: + sc = _decode(name) + for ch in ("Tran", "Vert", "Long"): + assert ch in sc, f"{name} missing {ch}" + tr = np.asarray(sc[ch], dtype=float) + tr = tr - tr[:4].mean() # reference to the pre-trigger baseline + assert 40 <= len(tr) <= 300, f"{name}:{ch} n={len(tr)}" + assert tr.min() < -8000, f"{name}:{ch} min {tr.min()}" + # deflects one way and rings back toward / past the baseline + assert tr.max() < abs(tr.min()), f"{name}:{ch} not one-sided" + + +def test_mic_present_only_on_four_channel_units(): + # UM11719 / UM12947 record a mic; UM13981 / UM20147 are 3-channel + # (mic-disabled) units and carry no MicL self-check. + got = {name: ("MicL" in _decode(name)) for name in EVENTS} + assert any(got.values()), "expected at least one 4-channel unit" + assert not all(got.values()), "expected at least one 3-channel unit" + for name, has_mic in got.items(): + if has_mic: + tr = np.asarray(_decode(name)["MicL"], dtype=float) + tr = tr - tr[:4].mean() + # mic self-check is a bipolar pulse train — swings both ways, wide range + assert tr.max() > 5000 and tr.min() < -5000, f"{name} mic not bipolar" + + +def test_channel_ids_and_order(): + # ids decode to the canonical channel names, geo always in Tran/Vert/Long order + sc = _decode(EVENTS[0]) + assert [c for c in ("Tran", "Vert", "Long") if c in sc] == ["Tran", "Vert", "Long"] + + +def test_returns_empty_on_non_idf_input(): + assert decode_idf_sensor_check(b"not an IDF file") == {} + assert decode_idf_sensor_check(b"") == {} diff --git a/tests/test_waveform_codec.py b/tests/test_waveform_codec.py index eebcf9d..d2c97b3 100644 --- a/tests/test_waveform_codec.py +++ b/tests/test_waveform_codec.py @@ -712,8 +712,27 @@ 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) - # NN > 8 is not a data block - assert data_block_len(b"\x40\x0c" + bytes(24), 0) == (None, None) + + +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) def test_record_chain_is_followed_by_length_not_by_tag_sniffing(): diff --git a/tests/test_waveform_fft.py b/tests/test_waveform_fft.py new file mode 100644 index 0000000..e1212fa --- /dev/null +++ b/tests/test_waveform_fft.py @@ -0,0 +1,85 @@ +"""Blastware-compatible channel FFT (waveform_fft). + +Reverse-engineered 2026-09-14 against 7 BE12844 (MiniMate Plus) events, each with +a Blastware FFT report as ground truth. The recipe (DC-remove, no window, +zero-pad to 4096 → 0.25 Hz bins, single-sided 2/N amplitude) reproduces +Blastware's dominant frequency to the exact bin on all 28 channels and the +amplitude to report precision. +""" +from pathlib import Path + +import numpy as np + +from waveform_fft import channel_spectrum, dominant_frequency +from minimateplus.waveform_codec import decode_waveform_v2 + +FIXDIR = Path(__file__).parent / "fixtures" / "fft-oracle-2026-09-14" +GEO_LSB = 0.005 # 1 decode unit = 16 ADC counts = 0.005 in/s (series-3 Normal range) + +# Blastware FFT-report ground truth: file → {channel: (dominant_hz, amplitude_ips)}. +# amplitude is None where the channel is at the noise floor (report amp 0.000/0.001) +# — the dominant frequency still matches exactly, but the amplitude isn't meaningful. +ORACLE = { + "N844LPGH.VV0W": {"Tran": (27.00, 0.018), "Vert": (26.75, 0.009), "Long": (26.50, 0.021), "MicL": (2.000, None)}, + "N844LPPR.3S0W": {"Tran": (30.75, None), "Vert": (46.75, None), "Long": (26.75, None), "MicL": (49.50, None)}, + "N844LQHB.ZT0W": {"Tran": (19.75, 0.040), "Vert": (26.50, 0.018), "Long": (26.50, 0.083), "MicL": (2.750, None)}, + "N844LQUE.T50W": {"Tran": (21.50, 0.080), "Vert": (14.25, 0.028), "Long": (28.50, 0.046), "MicL": (5.750, None)}, + "N844LR8W.790W": {"Tran": (31.00, None), "Vert": (31.00, None), "Long": (34.00, None), "MicL": (66.25, None)}, + "N844LRCO.G60W": {"Tran": (32.25, 0.009), "Vert": (32.00, 0.005), "Long": (32.00, 0.008), "MicL": (32.00, None)}, + "N844LRCW.F30W": {"Tran": (21.25, 0.010), "Vert": (42.25, 0.002), "Long": (21.25, 0.014), "MicL": (21.25, None)}, +} + + +def test_pure_sine_frequency_and_amplitude(): + # A pure sine at a bin-centre frequency (128 cycles over 4096 samples) has no + # leakage, so the single-sided 2/N normalisation returns the amplitude exactly. + sps, n, f0, amp = 1024.0, 4096, 32.0, 0.5 + x = amp * np.sin(2 * np.pi * f0 * np.arange(n) / sps) + freqs, amps = channel_spectrum(x, sps=sps, nfft=4096) + fpk, apk = dominant_frequency(freqs, amps) + assert fpk == 32.0 + assert abs(apk - amp) < 1e-3 + + +def test_bin_resolution_is_quarter_hz(): + freqs, _ = channel_spectrum(np.zeros(3328), sps=1024.0, nfft=4096) + assert abs((freqs[1] - freqs[0]) - 0.25) < 1e-9 + + +def test_empty_input(): + freqs, amps = channel_spectrum([]) + assert len(freqs) == 0 and len(amps) == 0 + + +def _spectra(fname): + raw = (FIXDIR / fname).read_bytes() + dec = decode_waveform_v2(raw[raw.find(b"STRT") + 21:]) + out = {} + for ch, samples in dec.items(): + ips = np.asarray(samples, float) * GEO_LSB + out[ch] = channel_spectrum(ips, sps=1024.0) + return out + + +def test_dominant_frequency_matches_blastware_exactly(): + misses = [] + for fname, chans in ORACLE.items(): + spectra = _spectra(fname) + for ch, (want_hz, _) in chans.items(): + got_hz, _ = dominant_frequency(*spectra[ch]) + if abs(got_hz - want_hz) > 0.25: + misses.append(f"{fname}:{ch} got {got_hz} want {want_hz}") + assert not misses, "dominant-frequency mismatches:\n" + "\n".join(misses) + + +def test_amplitude_matches_blastware(): + misses = [] + for fname, chans in ORACLE.items(): + spectra = _spectra(fname) + for ch, (_, want_amp) in chans.items(): + if want_amp is None: + continue + _, got_amp = dominant_frequency(*spectra[ch]) + if abs(got_amp - want_amp) > 0.0015: + misses.append(f"{fname}:{ch} got {got_amp:.4f} want {want_amp:.3f}") + assert not misses, "amplitude mismatches:\n" + "\n".join(misses) diff --git a/waveform_fft.py b/waveform_fft.py new file mode 100644 index 0000000..7e84565 --- /dev/null +++ b/waveform_fft.py @@ -0,0 +1,66 @@ +"""Blastware-compatible FFT of a decoded seismograph channel. + +Pure numpy; no I/O, no device or DB dependencies. Feed it a channel's decoded +samples **in the unit you want the amplitudes in** (e.g. in/s) and it returns the +single-sided amplitude spectrum that Blastware's *FFT Report* draws. + +Reverse-engineered 2026-09-14 against 7 BE12844 (MiniMate Plus) events with +Blastware FFT reports as ground truth. The recipe reproduces Blastware's +**dominant frequency to the exact 0.25 Hz bin on all 28 channels** and the +amplitude to report precision: + + 1. remove the DC component (subtract the mean); **no window** — a window + smears the peak and measurably worsens the match, + 2. zero-pad to ``nfft`` (4096 → 0.25 Hz bins at 1024 sps — Blastware's + resolution), + 3. single-sided amplitude ``A[k] = 2·|X[k]| / N`` where ``N`` is the real + sample count (not ``nfft``). + +The compliance chart (USBM RI8507 / OSMRE) is this spectrum's ``(freq, amp)`` +points plotted against the regulatory limit curve; the #10 FFT view is the +spectrum itself. +""" +from __future__ import annotations + +import numpy as np + +BW_NFFT = 4096 # 0.25 Hz bins at 1024 sps — Blastware's FFT resolution +BW_FMIN = 2.0 # dominant-frequency search floor (Hz) +BW_FMAX = 250.0 # dominant-frequency search ceiling (Hz) + + +def channel_spectrum(samples, sps: float = 1024.0, nfft: int = BW_NFFT): + """Single-sided amplitude spectrum of one channel, Blastware-compatible. + + ``samples`` is a 1-D sequence in the desired amplitude unit (in/s). Returns + ``(freqs, amps)`` numpy arrays covering ``0 .. sps/2`` in ``sps/nfft`` steps. + + Records longer than ``nfft`` are truncated by the transform — untested + against Blastware for that case (real MiniMate Plus records are ≤ ~3.3 s, + well under 4096 samples at 1024 sps). + """ + x = np.asarray(samples, dtype=float) + n = x.size + if n == 0: + return np.empty(0), np.empty(0) + x = x - x.mean() # DC removal, no window + mag = np.abs(np.fft.rfft(x, nfft)) + freqs = np.fft.rfftfreq(nfft, 1.0 / sps) + amps = (2.0 / n) * mag # single-sided amplitude + return freqs, amps + + +def dominant_frequency(freqs, amps, fmin: float = BW_FMIN, fmax: float = BW_FMAX): + """Peak ``(frequency_hz, amplitude)`` of a spectrum within ``[fmin, fmax)``. + + Matches Blastware's "Dominant Frequency" — the largest spectral bin in the + reportable band (below 2 Hz is baseline/DC drift, above 250 Hz is noise). + """ + freqs = np.asarray(freqs) + amps = np.asarray(amps) + lo = int(np.searchsorted(freqs, fmin)) + hi = int(np.searchsorted(freqs, fmax)) + if hi <= lo: + return 0.0, 0.0 + k = lo + int(np.argmax(amps[lo:hi])) + return float(freqs[k]), float(amps[k])