Compare commits
24
Commits
v0.29.0
...
0408c37866
+229
@@ -4,6 +4,235 @@ All notable changes to seismo-relay are documented here.
|
||||
|
||||
---
|
||||
|
||||
## Unreleased
|
||||
|
||||
### 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 (`minimateplus.sensor_check`).**
|
||||
The "Sensor Check" traces Blastware shows to the right of the waveform panel
|
||||
live in the series-3 binary's trailing block as four length-prefixed records
|
||||
(`0x3c`–`0x3f`) using the same delta-block codec as the main waveform:
|
||||
Tran/Vert/Long geophone ring-downs (the transducer's damped impulse response —
|
||||
resonant frequency + overswing/damping) and a MicL pulse train (the mic's
|
||||
known-signal gain check). `gather_report_data` decodes them from the retained
|
||||
BW binary at report time; the report renders them as a strip flush against the
|
||||
waveform panel plus the **Sensor Check → Frequency / Overswing Ratio** sub-rows
|
||||
in the stats table. Verified against the reports on all 7 oracle events (mic
|
||||
zero-crossing frequency = 20.1 Hz exact; geophone ring-downs consistent
|
||||
~7.5 Hz with overswing ~3.5). Tests in `tests/test_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
|
||||
|
||||
**None.** Every change here is additive and reads from data already on disk —
|
||||
the `.h5` samples and the retained raw BW binary. No `.h5`/DB change, no
|
||||
schema change, no migration, no backfill, and **no `TOOL_VERSION` bump**: a
|
||||
report regenerated for an existing event simply gains the new panels, and the
|
||||
`ach_server` rescue flags don't touch the codec.
|
||||
|
||||
---
|
||||
|
||||
## v0.30.0 — 2026-09-12
|
||||
|
||||
**The series-4 correctness release** — the Thor / Micromate counterpart to
|
||||
v0.26.0's series-3 work. The decoder is now verified per-sample against
|
||||
Thor's own CSV exports: **459 waveform files, 3,807,158 / 3,807,165 samples
|
||||
exact** across three independent ground-truth corpora, and production IDFW is
|
||||
**575/575** with zero truncations and zero decode failures. Series-3
|
||||
re-verified **unchanged at 14,338/14,338** after every shared-codec change.
|
||||
|
||||
⚠ **This release owes the prod store a Thor backfill.** Every stored
|
||||
series-4 geophone value is **3.3% low**, and histogram peaks from monitoring
|
||||
runs longer than ~4 hours can be far worse (the interval cap discarded the
|
||||
tail, frequently the part holding the peak). Run
|
||||
`scripts/backfill_thor_events.py` — `TOOL_VERSION` is bumped to `0.30.0`, so
|
||||
regeneration is gated correctly and **no `--force` is needed**. DB backup
|
||||
first. Series-3 events are untouched by this release and do not need
|
||||
re-running.
|
||||
|
||||
⚠ **Terra-View displays these values.** Series-4 geophone readings will rise
|
||||
~3.3% after the backfill, and some histogram PPVs will rise a great deal more.
|
||||
That is a correction, not a regression.
|
||||
|
||||
|
||||
### Fixed — event-report PDF used a per-trace geo Y scale
|
||||
|
||||
The waveform plot scaled each geo lane to its own peak, so a small channel
|
||||
filled its lane and looked as large as a big one, and the `Geo: X in/s/div`
|
||||
footer reflected only whichever channel was measured first — wrong for the
|
||||
other two. All three geo lanes now share one symmetric scale (max |sample|
|
||||
across them, padded, 0.05 in/s floor), matching the event modal and BW's
|
||||
single amp/div; the footer reflects that shared scale. Mic keeps its own psi
|
||||
scale. Large events are unchanged.
|
||||
|
||||
### Fixed — series-4 (Thor / Micromate) decoder is now per-sample exact
|
||||
|
||||
Verified against **Thor's own CSV exports**, which carry a per-sample
|
||||
four-column block beside every binary (`CSV/<name>.IDFW.csv`) — 1,012 paired
|
||||
files that had been sitting in the corpus unused. Previous notes asserted
|
||||
"Thor has no ASCII ground truth", which is why the decoder stayed pinned to a
|
||||
superseded walker with an unverifiable scale factor.
|
||||
|
||||
| metric | before | after |
|
||||
|---|---|---|
|
||||
| IDFW per-sample exact | 39.1% | **100.000%** (1,057,536/1,057,536) |
|
||||
| IDFW files fully exact | 0/153 | **153/153** |
|
||||
| IDFW PPV median error | −3.32% | **−0.002%** |
|
||||
| IDFH within 2% of Thor PPV | 51.1% | **100.0%** (858/858) |
|
||||
| prod IDFW PPV median error (8 units) | −3.3% | **−0.001%** |
|
||||
| decode cost | — | 6 ms/file |
|
||||
|
||||
Four independent root causes:
|
||||
|
||||
- **Geo LSB was `0.0003`, should be `0.000310308`** — the old value was Thor's
|
||||
4-decimal *display rounding* of the LSB mistaken for the LSB, so every
|
||||
series-4 geophone sample read **3.3% low**. Pinned to ±6e-11 by
|
||||
intersecting 991,415 rounding constraints; corroborated by the ±full-scale
|
||||
seed (`±32226`) in unwritten IDFH slots. Applies to IDFH too, which had a
|
||||
separate (also wrong) `10.0/32768`.
|
||||
- **IDFH histograms were capped at 250 intervals** — the segment validator
|
||||
required the interval counter's high byte to be zero, but the counter is a
|
||||
uint16 cumulative index, so every segment past interval 255 was rejected.
|
||||
Any run over ~4 hours lost its tail, often the part holding the peak.
|
||||
540/858 corpus files affected.
|
||||
- **Record mode `00 00` (raw int16, 10-byte header) was unhandled** — the
|
||||
record fell through the dispatch, silently dropping each channel's first
|
||||
512 samples. This produced the long-standing "loud events truncate"
|
||||
symptom. `MODE_ABSOLUTE` is now also accepted as a segment-0 preamble.
|
||||
- **Body-offset search matched `00 02 00` inside record headers** — picking a
|
||||
candidate part-way down the chain, which decodes a rotation-shifted body
|
||||
that drops each channel's segment 0. The search now anchors on record
|
||||
headers and takes the chain head.
|
||||
|
||||
Also fixes the separately-tracked "UM-series decodes ~1000× low" bug
|
||||
(`UM11402_20260406130113.IDFW` now matches its device report exactly).
|
||||
|
||||
Series-3 re-verified **unchanged at 14,338/14,338 exact** after the shared
|
||||
`waveform_codec` change.
|
||||
|
||||
⚠ **This is a codec change: the Thor store owes a regeneration.** Run
|
||||
`scripts/backfill_thor_events.py` (bump `TOOL_VERSION` first, or pass
|
||||
`--force`), DB backup first. All stored series-4 `.h5`/sidecar peaks are
|
||||
currently ~3.3% low, and histogram peaks for runs over ~4 hours may be
|
||||
badly low.
|
||||
|
||||
⚠ **Thor's histogram PPV has a 0.0050 in/s display floor** — 41.4% of prod
|
||||
IDFH sidecars report a component PPV larger than their own vector sum. On
|
||||
quiet files the decoder is now *more* accurate than that reference.
|
||||
|
||||
New: `scratch/verify_thor_against_csv.py`, `tests/test_idf_binary_codec.py`
|
||||
(10 tests, fixtures under `tests/fixtures/thor-idf/`).
|
||||
|
||||
### Fixed — mic-disabled (3-channel) units
|
||||
|
||||
Verified on a second corpus (`9-10-26-csv-req`: UM11402, UM12947, UM20147) —
|
||||
**139/139 waveforms per-sample exact (1,273,380 samples), 877/877 histograms
|
||||
within 2%** (was 66.9% and 56.6%).
|
||||
|
||||
- **Waveform body head sat below the scan floor.** A 3-channel unit's shorter
|
||||
header puts the record chain head at `0x0dba`, under the old
|
||||
`_BODY_SCAN_FLOOR` of `0x0E00`. The scan couldn't see it and fell through
|
||||
to the Vert segment-0 record, decoding a body shifted one position around
|
||||
the channel rotation — Vert came up exactly 512 samples short. Floor
|
||||
lowered to `0x0C00`; body-offset scoring now accepts 3 channels as "equal"
|
||||
instead of demanding 4.
|
||||
- **Histogram interval record is 56 bytes, not 72.** It is
|
||||
`16 × n_channels + 8`, so mic-disabled units pack 56. Assuming 72 read 7
|
||||
intervals out of every 10-interval segment then walked off alignment into
|
||||
garbage decoding as ~10 in/s peaks (errors up to +191,000%). The interval
|
||||
count now comes from the segment's cumulative counter and the stride is
|
||||
derived from it; also recovers 4 files that decoded no intervals at all.
|
||||
|
||||
Combined across both corpora: **292/292 waveform files, 2,330,916/2,330,916
|
||||
samples exact.** Production IDFW truncations 41 → 22.
|
||||
|
||||
### Fixed — `40 NN` int16 blocks with NN > 8
|
||||
|
||||
`data_block_len()` rejected any `40 NN` block with `NN > 0x08`. The cap had
|
||||
no evidence behind it: every corpus available when it was written used only
|
||||
NN ∈ {1,2,3,4,8}, so it was never exercised. Loud UM12947 events use NN of
|
||||
12, 16, 20 … up to 196, and because the block walker stops at the first
|
||||
unrecognised tag rather than raising, rejecting them surfaced as **silently
|
||||
short channels** (e.g. Tran 1812 / Vert 2132 / Long 2324 on a file whose
|
||||
export has 2324 for all three). The bound is the buffer, not a constant.
|
||||
|
||||
Verified against Thor exports for UM12947 (2025-07-14 … 09-25, 167
|
||||
waveforms): length mismatches **22 → 0**, **1,476,242/1,476,249** samples
|
||||
exact. These are not truncated recordings — the exports carry full sample
|
||||
counts.
|
||||
|
||||
`tests/test_waveform_codec.py` asserted the cap as intended behaviour; that
|
||||
assertion was wrong and has been replaced with one pinning the opposite,
|
||||
carrying the evidence.
|
||||
|
||||
### Result across all three ground-truth corpora
|
||||
|
||||
**459 waveform files, 3,807,158 / 3,807,165 samples exact.** Production
|
||||
IDFW: **575/575**, zero truncations, zero decode failures, median PPV error
|
||||
−0.0007% across 8 units. Series-3 re-verified **unchanged at 14,338/14,338**
|
||||
after every shared-codec change.
|
||||
|
||||
The 7 residual samples each differ by one 4th-decimal tick and are **Thor's
|
||||
own rounding**: intersecting the per-sample rounding constraints over that
|
||||
corpus is infeasible (the binding pair contradict by 2.3e-11, 7e-5 relative),
|
||||
so no single linear LSB reproduces every printed value. `_GEO_LSB_IPS` is
|
||||
already pinned to ~1e-11 — do not retune it to chase these.
|
||||
|
||||
---
|
||||
|
||||
## v0.29.0 — 2026-09-04
|
||||
|
||||
First release to reach prod since **v0.27.0**, so it ships **both** the
|
||||
|
||||
@@ -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.30.0**.
|
||||
|
||||
Stack-level context — which repo owns what, and how the three project versions
|
||||
pair — lives in `../terra-view/docs/tmi-stack.md`, which is also loaded as
|
||||
@@ -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/<name>.IDFW.csv` beside each binary with a **per-sample**
|
||||
four-column block. Harness: `scratch/verify_thor_against_csv.py`.
|
||||
Four bugs, all fixed: geo LSB was `0.0003` (display rounding of the real
|
||||
`0.000310308`, so every sample read **3.3% low**); the IDFH segment
|
||||
validator required a zero counter high byte, **capping every histogram at
|
||||
250 intervals**; record mode `00 00` (raw int16) was unhandled, silently
|
||||
dropping each channel's first 512 samples; and the body-offset search
|
||||
matched `00 02 00` *inside* record headers, decoding a rotation-shifted
|
||||
body. IDFW is no longer pinned to `decode_waveform_legacy`.
|
||||
Series-3 re-verified unchanged at 14,338/14,338 after the shared-codec
|
||||
change.
|
||||
- **Mic-disabled (3-channel) units are a distinct shape (2026-09-10).**
|
||||
Verified on a second corpus (`~/thor-csv-req`, UM11402/UM12947/UM20147):
|
||||
**139/139** waveforms per-sample exact, **877/877** histograms within 2%.
|
||||
Two structural differences: the shorter header puts the waveform record
|
||||
chain head at `0x0dba` (below the old `_BODY_SCAN_FLOOR` of `0x0E00`, so it
|
||||
was invisible and Vert came up exactly 512 short), and the histogram
|
||||
interval record is **56 bytes, not 72** — `16 × n_channels + 8`, derived per
|
||||
segment from the cumulative interval counter, never assumed.
|
||||
- **`40 NN` blocks are not capped at NN=8 (2026-09-11).** `data_block_len()`
|
||||
rejected `NN > 0x08`, a guard with no evidence behind it — the corpora
|
||||
available when it was written only used NN ∈ {1,2,3,4,8}. Loud UM12947
|
||||
events use NN up to 196, and since the walker stops at the first
|
||||
unrecognised tag rather than raising, this surfaced as silently short
|
||||
channels. Verified on 167 UM12947 waveforms: length mismatches 22 → 0,
|
||||
1,476,242/1,476,249 samples exact.
|
||||
- **Production IDFW is now 575/575** — zero truncations, zero decode
|
||||
failures, median PPV error −0.0007% across 8 units (was 41 truncated + 1
|
||||
failing, −3.3%). Across all three ground-truth corpora: **459 files,
|
||||
3,807,158/3,807,165 samples exact**; the 7 stragglers differ by one
|
||||
4th-decimal tick and are Thor's own rounding — no single linear LSB can
|
||||
reproduce every printed value (the constraints are infeasible by 7e-5
|
||||
relative), so do NOT retune `_GEO_LSB_IPS`.
|
||||
- **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..<branch>` before you merge, or
|
||||
`git log --oneline <merge-base>..<branch>` 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 — <theme>` 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 `<chain-head record> + 7`, found
|
||||
by `_find_waveform_body_offset()` anchoring on record headers. All
|
||||
**153/153** genuine Thor waveform files decode per-sample exact
|
||||
(1,057,536/1,057,536 samples).
|
||||
- **IDFH** segment header is `[len_be][0a 00 00 00][counter_be][05 3f]`,
|
||||
where `counter` is a **uint16 cumulative interval index** — it must
|
||||
not be constrained to a zero high byte (that capped histograms at 250
|
||||
intervals). Intervals whose `min > max` on all channels are unwritten
|
||||
slots carrying a ±full-scale seed and are skipped. 858/858 files land
|
||||
within 2% of Thor's PPV (median -0.004%).
|
||||
- **Geo LSB is `0.000310308` in/s per count** (full scale 10.0 in/s =
|
||||
32226.05 counts). Series-3's 32000-count scale does NOT apply.
|
||||
- **Record modes** are `02 00` deltas (14 B header), `01 00` absolute,
|
||||
`00 03` raw 12-bit, and `00 00` **raw int16** (all 10 B headers).
|
||||
`01 00` and `00 00` are also valid as the implicit segment-0 preamble.
|
||||
|
||||
⚠ **Thor's histogram PPV has a 0.0050 in/s display floor.** 41.4% of
|
||||
prod IDFH sidecars report a component PPV exceeding their own vector
|
||||
sum — impossible. On quiet files our decode is *more* accurate than
|
||||
the reference; do not "fix" the decoder to match it.
|
||||
|
||||
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)
|
||||
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# seismo-relay `v0.29.0`
|
||||
# seismo-relay `v0.30.0`
|
||||
|
||||
A ground-up replacement for **Blastware** — Instantel's aging Windows-only
|
||||
software for managing seismographs. Supports both the **MiniMate Plus
|
||||
|
||||
@@ -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",
|
||||
|
||||
@@ -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:
|
||||
|
||||
```
|
||||
<serial dir>/UM13981_20220207084555.IDFW
|
||||
<serial dir>/CSV/UM13981_20220207084555.IDFW.csv
|
||||
```
|
||||
|
||||
The **CSV carries a per-sample block** — four columns (Tran, Vert, Long, Mic)
|
||||
in in/s and psi, after the 2-column report header. That is the series-4
|
||||
equivalent of Blastware's `_ASCII.TXT` exports, and it gives 1,012 paired
|
||||
files (152 IDFW + 860 IDFH). Earlier notes in this file and in
|
||||
`micromate/idf_file.py` asserted "Thor has no ASCII ground truth in the
|
||||
corpus"; that was wrong, and it is why the decoder sat pinned to a
|
||||
superseded walker with a scaling constant nobody could check.
|
||||
|
||||
Harness: `scratch/verify_thor_against_csv.py`.
|
||||
|
||||
### Geo LSB = 0.000310308 in/s per count (NOT 0.0003)
|
||||
|
||||
The old 0.0003 was read off the smallest non-zero sample in the exports —
|
||||
but that is Thor's **4-decimal display rounding of the LSB, not the LSB**.
|
||||
It read every series-4 geophone sample **3.3% low**. The quantisation
|
||||
ladder gives it away: counts 1..6 export as 0.0003, 0.0006, 0.0009, 0.0012,
|
||||
0.0016, 0.0019 — an LSB of exactly 0.0003 would end 0.0015, 0.0018.
|
||||
|
||||
Each exported sample constrains the LSB to the window that rounds to its
|
||||
printed value. Intersecting 991,415 such constraints gives
|
||||
|
||||
```
|
||||
LSB ∈ [0.000310307933, 0.000310308057] width 1.2e-10
|
||||
```
|
||||
|
||||
so `_GEO_LSB_IPS = 0.000310308`, i.e. full scale 10.0 in/s = **32226.05
|
||||
counts**. Corroboration: an IDFH interval that never recorded keeps its
|
||||
min/max accumulator at its ±full-scale seed, and that seed is
|
||||
`(min=+32226, max=-32226)`. ⚠ The tempting closed form `10.0/32226` is
|
||||
very slightly wrong — it lands 4.5e-10 above the feasible window and loses
|
||||
78 boundary samples while never winning one. **Series III uses 32000 counts
|
||||
for the same 10.0 in/s, so the two generations do not share a scale.**
|
||||
|
||||
Independently confirmed on 8 production units (UM6047, UM11402, UM11719,
|
||||
UM12947, UM13981, UM14133, UM20146, UM20147): every unit's median PPV error
|
||||
against its device-reported peak moved from −3.3% to within ±0.03%. It is a
|
||||
global constant, not a per-unit calibration.
|
||||
|
||||
### IDFH segment header: the counter is a uint16, and it is cumulative
|
||||
|
||||
```
|
||||
[length_be 2B][0a 00 00 00][counter_be 2B][05 3f]
|
||||
```
|
||||
|
||||
`counter` is the **0-based cumulative index of the last interval in the
|
||||
segment** — 9, 19, 29, ... for the usual 10-intervals-per-segment layout
|
||||
(`length` = 730).
|
||||
|
||||
The validator used to require `counter`'s high byte to be `0x00`. That
|
||||
silently **capped every histogram at 250 intervals**: once the cumulative
|
||||
counter passed 255 the high byte went non-zero and every later segment was
|
||||
rejected. Any run longer than ~4 hours lost its tail — frequently the part
|
||||
holding the event peak, so the file's PPV read low. **540 of 858 corpus
|
||||
files were affected**; fixing it moved histogram peaks from 48.3% to 93.8%
|
||||
within 0.5% of Thor's reported PPV.
|
||||
|
||||
### Unwritten interval slots carry a ±full-scale seed
|
||||
|
||||
An interval the device reserved but never wrote keeps `min = +32226`,
|
||||
`max = -32226` on all four channels — `min > max`, impossible for real data.
|
||||
Decoded naively it yields a 10.0 in/s peak on every channel and, being a
|
||||
max-over-intervals, poisons the whole file's PPV. Rare but real: exactly 1
|
||||
of 497,611 corpus intervals, and it inflated that file's Long PPV from
|
||||
0.0081 to 10.0 in/s. The inversion is all-or-nothing across channels (0
|
||||
partial cases), so requiring every channel to be inverted is a safe test.
|
||||
|
||||
### Record mode `00 00` — raw int16 absolute (MODE_RAW16)
|
||||
|
||||
The record chain's mode field at `off+8` takes a fourth value:
|
||||
|
||||
| mode | meaning | header |
|
||||
|---|---|---|
|
||||
| `02 00` | deltas + two int16 anchors | 14 B |
|
||||
| `01 00` | absolute, tagged blocks | 10 B |
|
||||
| `00 03` | raw 12-bit absolute, untagged | 10 B |
|
||||
| **`00 00`** | **raw int16 BE absolute, untagged** | **10 B** |
|
||||
|
||||
A `MODE_RAW16` record with `length = 1032` carries exactly
|
||||
`(1032 - 8) / 2 = 512` samples and reproduced Thor's export **512/512
|
||||
exactly** on first test. Thor uses it for segment 0 (the pre-trigger
|
||||
window) on some events. Before this mode existed the record fell through
|
||||
the dispatch unhandled, so the channel silently lost its first 512 samples —
|
||||
which is what produced the "loud events truncate" symptom.
|
||||
|
||||
`MODE_ABSOLUTE` is also valid as a **preamble** (the implicit segment-0 Tran
|
||||
record); its tagged blocks start at `body[3]`, not `body[7]`, because its
|
||||
header is 10 bytes rather than 14.
|
||||
|
||||
### Body offset is not fixed at 0x0f1f — and 0x0f1f is really a record + 7
|
||||
|
||||
A "body offset" is `<record start> + 7`, so that `body[0]` is the segment
|
||||
index and `body[1:3]` is the mode. The canonical `0x0f1f` is simply the
|
||||
record at `0x0f18`.
|
||||
|
||||
Searching for the literal preamble `00 02 00` finds only MODE_DELTA bodies,
|
||||
and worse, it **matches the `[seg][mode]` bytes inside any record header**,
|
||||
so the scan could pick a candidate part-way down the chain. That decodes a
|
||||
plausible-looking but rotation-shifted body which drops each channel's
|
||||
segment 0 — the real cause of the remaining truncations.
|
||||
|
||||
`_find_waveform_body_offset()` now anchors on record headers (the
|
||||
`<channel_id> 00 00` signature at `+4`, validated with `is_record()`),
|
||||
takes the **chain head** — a record no other record's length field points at
|
||||
— and trial-decodes `head + 7`, preferring the candidate where all four
|
||||
channels come out the same length.
|
||||
|
||||
⚠ Do **not** scan for candidate preambles instead: `MODE_RAW16` is
|
||||
`00 00`, so every run of three zero bytes looks like a body start and each
|
||||
costs a full trial decode (~0.5 s/file measured, vs 6 ms/file now).
|
||||
|
||||
### `40 NN` is not capped at NN=8 (2026-09-11)
|
||||
|
||||
`data_block_len()` rejected any `40 NN` int16 block with `NN > 0x08`. The cap
|
||||
had no evidence behind it — every corpus available when it was written used
|
||||
only NN ∈ {1, 2, 3, 4, 8}, so it was never exercised. Loud events use much
|
||||
wider blocks:
|
||||
|
||||
| corpus | `40 NN` values | walker stops |
|
||||
|---|---|---|
|
||||
| first + 3-channel corpora | 1, 2, 3, 4, 8 | none |
|
||||
| UM12947 2025-07..09 | 2, 4, 8, **12, 16, 20 … 196** | every value > 8 |
|
||||
|
||||
Because `walk_body`/`run` stop at the first unrecognised tag rather than
|
||||
raising, this surfaced as **silently short channels** — e.g. Tran 1812 /
|
||||
Vert 2132 / Long 2324 on a file whose export has 2324 for all three. The
|
||||
real bound is the buffer (and the caller's record end), not a magic constant.
|
||||
|
||||
Verified against Thor's exports for UM12947 (2025-07-14 … 2025-09-25, 167
|
||||
waveforms): length mismatches **22 → 0**, and **1,476,242 / 1,476,249**
|
||||
samples exact.
|
||||
|
||||
⚠ These events are **not** truncated recordings, which was the competing
|
||||
hypothesis — the exports carry the full sample count.
|
||||
|
||||
**The 7 residual samples are Thor's rounding, not ours.** Each differs by
|
||||
exactly one 4th-decimal tick (e.g. decoded 3.3551 vs export 3.3550).
|
||||
Intersecting the per-sample rounding constraints over this corpus is
|
||||
**infeasible** — the binding pair (count 2013 → 0.6247, count 4351 → 1.3501)
|
||||
contradict by 2.3e-11, i.e. 7e-5 relative. No single linear LSB can
|
||||
reproduce every printed value, so Thor is not doing plain round-half-up on
|
||||
`count × LSB`. Do not retune `_GEO_LSB_IPS` to chase these; it is already
|
||||
pinned to ~1e-11.
|
||||
|
||||
### Mic-disabled units are a distinct shape (2026-09-10, second corpus)
|
||||
|
||||
Some units run with the microphone disabled — **3 channels, not 4** — and that
|
||||
changes two structural things. Confirmed on the `9-10-26-csv-req` corpus
|
||||
(UM11402, UM12947, UM20147): 139/139 waveforms and 877/877 histograms.
|
||||
|
||||
**Waveform: the body starts earlier.** A 3-channel unit has a shorter fixed
|
||||
header and puts its record chain head at **`0x0dba`**, below the old
|
||||
`_BODY_SCAN_FLOOR` of `0x0E00`. The head was therefore invisible to the scan,
|
||||
which fell through to the *Vert* segment-0 record and decoded a body shifted
|
||||
one position around the channel rotation. The signature is unmistakable:
|
||||
|
||||
```
|
||||
Tran 3072 / Vert 2560 / Long 3072 / MicL 0 <- Vert exactly 512 short
|
||||
```
|
||||
|
||||
46 of 139 files in that corpus were affected; all 46 became per-sample exact
|
||||
once the floor dropped to `0x0C00`. Note the body-offset scoring also had to
|
||||
stop requiring four channels — `len(lengths) >= 3`, not `== 4`, or `equal` is
|
||||
permanently False for these events and the pick falls back to raw sample count.
|
||||
|
||||
**Histogram: the interval record is 56 bytes, not 72.**
|
||||
|
||||
```
|
||||
interval_size = 16 × n_channels + 8 (72 for 4 channels, 56 for 3)
|
||||
```
|
||||
|
||||
It is **not a constant**, and it cannot be inferred from `length` alone.
|
||||
Derive the interval count from the segment counter — it is cumulative, so
|
||||
`n = counter - previous_counter` — and then `stride = (length - 10) / n`.
|
||||
`n_channels` follows from `(stride - 8) / 16`.
|
||||
|
||||
Assuming 72 read 7 intervals out of each 10-interval segment and then walked
|
||||
off alignment into garbage that decoded as ~10 in/s peaks — inflating those
|
||||
files' PPV by up to 191,000%. Fixing it moved the second corpus from 56.6% to
|
||||
**100.0%** of histograms within 2% of Thor's reported PPV, and recovered 4
|
||||
files that previously decoded no intervals at all.
|
||||
|
||||
### What is still open
|
||||
|
||||
- ~~23 of 575 production IDFW files~~ — **RESOLVED 2026-09-11.** Production
|
||||
IDFW is now **575/575** with zero truncations and zero decode failures
|
||||
(median PPV error −0.0007%). See "`40 NN` is not capped at NN=8" above.
|
||||
|
||||
- Mic → psi scale is still the rough `2.14e-6` regression, not derived.
|
||||
- Per-channel `int16 field4` in the IDFH interval record (possibly
|
||||
time-of-peak) and the 8-byte tail (PVS data) remain undecoded.
|
||||
|
||||
⚠ **Thor's histogram PPV has a display floor of 0.0050 in/s.** In the
|
||||
production store 6,080 sidecar PPV values are exactly 0.0050 (next most
|
||||
common value: 275 occurrences), and **41.4% of IDFH sidecars report a
|
||||
component PPV larger than their own vector sum** — geometrically impossible.
|
||||
On those quiet files the decoder's ~0.0025 in/s is *more* accurate than the
|
||||
reference; do not "fix" the decoder to match it.
|
||||
|
||||
### Codec breakthroughs (2026-05-28)
|
||||
|
||||
- **Body offset is a fixed `0x0f1f`** across 151/154 corpus IDFW
|
||||
|
||||
@@ -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 —
|
||||
<https://www.ecfr.gov/current/title-30/chapter-VII/subchapter-K/part-816/section-816.67>
|
||||
@@ -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/<unit>-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 <host> <port> --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/<unit>-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.
|
||||
|
||||
@@ -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/<unit>-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.
|
||||
+221
-40
@@ -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
|
||||
|
||||
# 1. Locate every plausible per-channel record header. A header carries
|
||||
# [len 2B][channel_id][00][00] at +2..+6, so anchor the search on the
|
||||
# three-byte ``<cid> 00 00`` signature and validate with is_record().
|
||||
# Scanning candidate *preambles* instead is not viable: MODE_RAW16 is
|
||||
# ``00 00``, so every run of three zero bytes would look like a body
|
||||
# start and each would cost a full trial decode (~0.5 s/file measured).
|
||||
floor = max(0, _BODY_SCAN_FLOOR - 7)
|
||||
starts: list = []
|
||||
for cid in (0x46, 0x47, 0x48, 0x49):
|
||||
sig = bytes((cid, 0x00, 0x00))
|
||||
i = floor
|
||||
while True:
|
||||
j = buf.find(_BODY_MAGIC, i)
|
||||
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,
|
||||
|
||||
@@ -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)
|
||||
@@ -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.30.0"
|
||||
|
||||
try:
|
||||
# Best-effort: prefer the installed metadata when it's NEWER than the
|
||||
|
||||
@@ -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
|
||||
@@ -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
|
||||
|
||||
+1
-1
@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
|
||||
|
||||
[project]
|
||||
name = "seismo-relay"
|
||||
version = "0.29.0"
|
||||
version = "0.30.0"
|
||||
description = "Python client and REST server for MiniMate Plus seismographs"
|
||||
requires-python = ">=3.10"
|
||||
dependencies = [
|
||||
|
||||
@@ -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:
|
||||
|
||||
<dir>/UM13981_20220207084555.IDFW
|
||||
<dir>/CSV/UM13981_20220207084555.IDFW.csv
|
||||
|
||||
For waveforms the CSV carries a per-sample block of four columns
|
||||
(Tran, Vert, Long, Mic) in in/s and psi -- i.e. true per-sample ground truth,
|
||||
exactly what the BW ASCII exports give us for series-3. The leading 2-column
|
||||
rows are the report header (PPV, sample rate, geo range, ...).
|
||||
|
||||
Usage:
|
||||
python scratch/verify_thor_against_csv.py [--root DIR] [--lsb FLOAT]
|
||||
[--limit N] [--kind idfw|idfh|both]
|
||||
"""
|
||||
from __future__ import annotations
|
||||
|
||||
import argparse
|
||||
import csv
|
||||
import os
|
||||
import statistics
|
||||
import sys
|
||||
from collections import Counter, defaultdict
|
||||
|
||||
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
|
||||
|
||||
from micromate import idf_file as M
|
||||
|
||||
DEFAULT_ROOT = "/home/serversdown/thor-watcher/example-data"
|
||||
GEO = ("Tran", "Vert", "Long")
|
||||
|
||||
|
||||
def parse_export(path):
|
||||
"""Return (header_dict, sample_rows) from a Thor CSV export."""
|
||||
hdr, rows = {}, []
|
||||
with open(path, newline="", encoding="utf-8", errors="replace") as fh:
|
||||
for rec in csv.reader(fh):
|
||||
if len(rec) == 2:
|
||||
hdr[rec[0].strip()] = rec[1].strip()
|
||||
elif len(rec) >= 3:
|
||||
try:
|
||||
rows.append([float(x) for x in rec])
|
||||
except ValueError:
|
||||
pass
|
||||
return hdr, rows
|
||||
|
||||
|
||||
def index_corpus(root):
|
||||
"""Map BASENAME.IDFW -> (binary_path, csv_path) for every paired file."""
|
||||
exports, binaries = {}, {}
|
||||
for dirpath, _dirs, files in os.walk(root):
|
||||
for name in files:
|
||||
up = name.upper()
|
||||
full = os.path.join(dirpath, name)
|
||||
if up.endswith(".IDFW.CSV") or up.endswith(".IDFH.CSV"):
|
||||
exports.setdefault(name[:-4].upper(), full)
|
||||
elif up.endswith(".IDFW") or up.endswith(".IDFH"):
|
||||
binaries.setdefault(up, full)
|
||||
return {k: (binaries[k], exports[k]) for k in binaries.keys() & exports.keys()}
|
||||
|
||||
|
||||
def hdr_float(hdr, key):
|
||||
raw = hdr.get(key)
|
||||
if not raw:
|
||||
return None
|
||||
try:
|
||||
return float(raw.split()[0])
|
||||
except (ValueError, IndexError):
|
||||
return None
|
||||
|
||||
|
||||
def verify_waveform(binpath, csvpath, lsb):
|
||||
"""Compare one IDFW against its export. Returns a result dict."""
|
||||
out = {"file": os.path.basename(binpath), "status": "ok"}
|
||||
try:
|
||||
res = M.read_idf_file(binpath)
|
||||
except NotImplementedError:
|
||||
out["status"] = "not-thor"
|
||||
return out
|
||||
except Exception as exc: # noqa: BLE001 - harness reports, never raises
|
||||
out["status"] = "decode-error"
|
||||
out["detail"] = f"{type(exc).__name__}: {exc}"
|
||||
return out
|
||||
|
||||
hdr, rows = parse_export(csvpath)
|
||||
if not rows:
|
||||
out["status"] = "no-gt-samples"
|
||||
return out
|
||||
|
||||
gt = {ch: [r[i] for r in rows] for i, ch in enumerate(GEO)}
|
||||
out["gt_len"] = len(rows)
|
||||
out["geo_range"] = hdr.get("GeoRange")
|
||||
|
||||
exact = total = 0
|
||||
lens, chan_status = {}, {}
|
||||
ppv_err = {}
|
||||
for ch in GEO:
|
||||
arr = res.samples.get(ch, [])
|
||||
ref = gt[ch]
|
||||
lens[ch] = len(arr)
|
||||
if len(arr) != len(ref):
|
||||
chan_status[ch] = "length"
|
||||
continue
|
||||
if not arr:
|
||||
chan_status[ch] = "empty"
|
||||
continue
|
||||
hits = sum(1 for c, v in zip(arr, ref) if abs(c * lsb - v) < 5e-5)
|
||||
exact += hits
|
||||
total += len(arr)
|
||||
chan_status[ch] = "exact" if hits == len(arr) else "value"
|
||||
gp = hdr_float(hdr, f"{ch}PPV")
|
||||
if gp:
|
||||
ppv_err[ch] = (max(abs(c) for c in arr) * lsb - gp) / gp
|
||||
|
||||
out["lens"] = lens
|
||||
out["chan_status"] = chan_status
|
||||
out["exact"] = exact
|
||||
out["total"] = total
|
||||
out["ppv_err"] = ppv_err
|
||||
if all(v == "exact" for v in chan_status.values()):
|
||||
out["status"] = "exact"
|
||||
elif any(v == "length" for v in chan_status.values()):
|
||||
out["status"] = "length-mismatch"
|
||||
else:
|
||||
out["status"] = "value-mismatch"
|
||||
return out
|
||||
|
||||
|
||||
def verify_histogram(binpath, csvpath, lsb):
|
||||
out = {"file": os.path.basename(binpath), "status": "ok"}
|
||||
try:
|
||||
res = M.read_idf_file(binpath)
|
||||
except NotImplementedError:
|
||||
out["status"] = "not-thor"
|
||||
return out
|
||||
except Exception as exc: # noqa: BLE001
|
||||
out["status"] = "decode-error"
|
||||
out["detail"] = f"{type(exc).__name__}: {exc}"
|
||||
return out
|
||||
hdr, _rows = parse_export(csvpath)
|
||||
out["n_intervals"] = len(res.intervals or [])
|
||||
errs = {}
|
||||
for ch, attr in (("Tran", "transverse_ips"), ("Vert", "vertical_ips"),
|
||||
("Long", "longitudinal_ips")):
|
||||
gp = hdr_float(hdr, f"{ch}PPV")
|
||||
dv = getattr(res.event.peaks, attr, None)
|
||||
if gp and dv:
|
||||
errs[ch] = (dv - gp) / gp
|
||||
out["ppv_err"] = errs
|
||||
out["status"] = "peaks" if errs else "no-gt-peaks"
|
||||
return out
|
||||
|
||||
|
||||
def main():
|
||||
ap = argparse.ArgumentParser()
|
||||
ap.add_argument("--root", default=DEFAULT_ROOT)
|
||||
ap.add_argument("--lsb", type=float, default=M._GEO_LSB_IPS)
|
||||
ap.add_argument("--limit", type=int, default=0)
|
||||
ap.add_argument("--kind", choices=("idfw", "idfh", "both"), default="both")
|
||||
ap.add_argument("--show", type=int, default=15, help="worst-N detail rows")
|
||||
args = ap.parse_args()
|
||||
|
||||
pairs = index_corpus(args.root)
|
||||
keys = sorted(pairs)
|
||||
if args.kind != "both":
|
||||
keys = [k for k in keys if k.endswith(args.kind.upper())]
|
||||
if args.limit:
|
||||
keys = keys[: args.limit]
|
||||
|
||||
print(f"root: {args.root}")
|
||||
print(f"geo LSB under test: {args.lsb!r} in/s per count")
|
||||
print(f"paired files: {len(keys)}\n")
|
||||
|
||||
wf, hg = [], []
|
||||
for k in keys:
|
||||
binpath, csvpath = pairs[k]
|
||||
if k.endswith(".IDFW"):
|
||||
wf.append(verify_waveform(binpath, csvpath, args.lsb))
|
||||
else:
|
||||
hg.append(verify_histogram(binpath, csvpath, args.lsb))
|
||||
|
||||
if wf:
|
||||
st = Counter(r["status"] for r in wf)
|
||||
ex = sum(r.get("exact", 0) for r in wf)
|
||||
tot = sum(r.get("total", 0) for r in wf)
|
||||
print("=" * 68)
|
||||
print(f"WAVEFORM (IDFW): {len(wf)} files")
|
||||
for s, n in st.most_common():
|
||||
print(f" {s:16} {n:5d} ({100*n/len(wf):5.1f}%)")
|
||||
if tot:
|
||||
print(f" per-sample exact: {ex}/{tot} = {100*ex/tot:.3f}%")
|
||||
errs = [e for r in wf for e in r.get("ppv_err", {}).values()]
|
||||
if errs:
|
||||
print(f" PPV rel-error: median {statistics.median(errs):+.4%} "
|
||||
f"mean {statistics.mean(errs):+.4%} "
|
||||
f"max|.| {max(abs(e) for e in errs):.4%}")
|
||||
bad = [r for r in wf if r["status"] not in ("exact",)]
|
||||
if bad:
|
||||
print(f"\n worst {min(args.show, len(bad))} of {len(bad)} non-exact:")
|
||||
for r in bad[: args.show]:
|
||||
print(f" {r['file']:42} {r['status']:16} "
|
||||
f"lens={r.get('lens')} gt={r.get('gt_len')} "
|
||||
f"{r.get('detail','')}")
|
||||
|
||||
if hg:
|
||||
st = Counter(r["status"] for r in hg)
|
||||
print("=" * 68)
|
||||
print(f"HISTOGRAM (IDFH): {len(hg)} files")
|
||||
for s, n in st.most_common():
|
||||
print(f" {s:16} {n:5d} ({100*n/len(hg):5.1f}%)")
|
||||
errs = [e for r in hg for e in r.get("ppv_err", {}).values()]
|
||||
if errs:
|
||||
print(f" PPV rel-error: median {statistics.median(errs):+.4%} "
|
||||
f"mean {statistics.mean(errs):+.4%} "
|
||||
f"max|.| {max(abs(e) for e in errs):.4%}")
|
||||
within = lambda t: 100*sum(1 for e in errs if abs(e) <= t)/len(errs)
|
||||
print(f" within 0.5%: {within(0.005):.1f}% "
|
||||
f"within 2%: {within(0.02):.1f}% within 5%: {within(0.05):.1f}%")
|
||||
return 0
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
raise SystemExit(main())
|
||||
@@ -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)
|
||||
|
||||
|
||||
@@ -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")
|
||||
+182
-41
@@ -121,6 +121,11 @@ class ReportData:
|
||||
t0_ms: Optional[float] = None
|
||||
dt_ms: Optional[float] = None
|
||||
|
||||
# Sensor self-check traces — {ch: [samples]} in raw decode units, decoded
|
||||
# from the binary's trailing block (see minimateplus.sensor_check). The
|
||||
# little waveforms BW draws in its "Sensor Check" strip. Empty when absent.
|
||||
sensor_check_waveforms: dict = field(default_factory=dict)
|
||||
|
||||
# Record-type discriminator
|
||||
record_type: Optional[str] = None
|
||||
is_histogram: bool = False
|
||||
@@ -246,6 +251,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,
|
||||
})
|
||||
@@ -290,6 +297,19 @@ def gather_report_data(
|
||||
except Exception as exc:
|
||||
log.warning("gather_report_data: hdf5 read failed: %s", exc)
|
||||
|
||||
# ── Sensor self-check traces — decoded from the retained raw binary ──
|
||||
# The .h5 holds only the main waveform; the sensor-check traces live in the
|
||||
# binary's trailing block, so decode them straight from the kept BW file.
|
||||
# Waveform events only (histograms have no sensor-check strip).
|
||||
if not rd.is_histogram:
|
||||
try:
|
||||
from minimateplus.sensor_check import decode_sensor_check
|
||||
bw_path, _a5 = store.paths_for(serial, filename)
|
||||
if bw_path.exists():
|
||||
rd.sensor_check_waveforms = decode_sensor_check(bw_path.read_bytes())
|
||||
except Exception as exc:
|
||||
log.warning("gather_report_data: sensor-check decode failed: %s", exc)
|
||||
|
||||
# ── Histogram aggregation ──
|
||||
# Codec emits ~N per-block samples (typically 1/sec); BW reports
|
||||
# one bar per configured interval (1 min / 5 min / etc.). When
|
||||
@@ -396,9 +416,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 +522,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 +589,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 +622,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 +679,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 +751,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.
|
||||
# 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 +799,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 +825,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 +845,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,7 +905,10 @@ 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)
|
||||
# "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")
|
||||
|
||||
@@ -814,23 +925,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",
|
||||
|
||||
+13
-2
@@ -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
|
||||
|
||||
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
@@ -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
|
||||
@@ -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
|
||||
@@ -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="<thor-watcher>/example-data/THORDATA_example/THORDATA_example/UPMC Presby"
|
||||
# mkdir -p tests/fixtures/thor-idf
|
||||
# for f in UM11719/UM11719_20231219162723.IDFW \
|
||||
# UM11719/UM11719_20231219162648.IDFH \
|
||||
# UM13981/UM13981_20220207084555.IDFW \
|
||||
# UM13981/UM13981_20220207183102.IDFH \
|
||||
# UM13981/UM13981_20221202063059.IDFH; do
|
||||
# cp "$B/$f" tests/fixtures/thor-idf/
|
||||
# cp "$B/$(dirname $f)/CSV/$(basename $f).csv" tests/fixtures/thor-idf/
|
||||
# done
|
||||
pytestmark = pytest.mark.skipif(
|
||||
not FIXTURES.is_dir() or not any(FIXTURES.glob("*.IDFW")),
|
||||
reason=f"Thor IDF fixtures not present under {FIXTURES}",
|
||||
)
|
||||
|
||||
|
||||
def _parse_export(path: Path):
|
||||
"""Split a Thor CSV export into (header dict, per-sample rows)."""
|
||||
header, rows = {}, []
|
||||
with path.open(newline="", encoding="utf-8", errors="replace") as fh:
|
||||
for rec in csv.reader(fh):
|
||||
if len(rec) == 2:
|
||||
header[rec[0].strip()] = rec[1].strip()
|
||||
elif len(rec) >= 3:
|
||||
try:
|
||||
rows.append([float(x) for x in rec])
|
||||
except ValueError:
|
||||
pass
|
||||
return header, rows
|
||||
|
||||
|
||||
def _header_float(header, key):
|
||||
return float(header[key].split()[0])
|
||||
|
||||
|
||||
@pytest.fixture(scope="module")
|
||||
def idfw_export():
|
||||
return _parse_export(IDFW.with_suffix(".IDFW.csv"))
|
||||
|
||||
|
||||
# ─── The geo scale constant ────────────────────────────────────────────────
|
||||
|
||||
|
||||
def test_geo_lsb_matches_thor_quantisation():
|
||||
"""Thor's own export quantises geo samples to this LSB.
|
||||
|
||||
Derived by maximising exact-match count over 1,046,016 paired samples
|
||||
(454 channel-events, 2 units); independently corroborated on 8
|
||||
production units via their device-reported PPV. The historical value
|
||||
0.0003 read every series-4 geophone sample 3.3% low.
|
||||
"""
|
||||
assert _GEO_LSB_IPS == pytest.approx(0.000310308, rel=1e-6)
|
||||
|
||||
|
||||
def test_geo_lsb_is_not_the_legacy_value():
|
||||
# Guards against a revert to the truncated 0.0003 constant.
|
||||
assert abs(_GEO_LSB_IPS - 0.0003) > 1e-6
|
||||
|
||||
|
||||
# ─── Per-sample fidelity ───────────────────────────────────────────────────
|
||||
|
||||
|
||||
def test_waveform_channel_lengths_match_export(idfw_export):
|
||||
_header, rows = idfw_export
|
||||
result = read_idf_file(IDFW)
|
||||
for channel in GEO_CHANNELS:
|
||||
assert len(result.samples[channel]) == len(rows), (
|
||||
f"{channel} truncated: decoded {len(result.samples[channel])} "
|
||||
f"samples, export has {len(rows)}"
|
||||
)
|
||||
|
||||
|
||||
def test_waveform_samples_match_export_exactly(idfw_export):
|
||||
"""Every geo sample must reproduce Thor's exported value to 4 dp."""
|
||||
_header, rows = idfw_export
|
||||
result = read_idf_file(IDFW)
|
||||
for index, channel in enumerate(GEO_CHANNELS):
|
||||
decoded = result.samples[channel]
|
||||
expected = [row[index] for row in rows]
|
||||
mismatches = [
|
||||
(i, geo_count_to_ips(c), v)
|
||||
for i, (c, v) in enumerate(zip(decoded, expected))
|
||||
if abs(geo_count_to_ips(c) - v) >= 5e-5
|
||||
]
|
||||
assert not mismatches, (
|
||||
f"{channel}: {len(mismatches)} of {len(expected)} samples differ; "
|
||||
f"first three {mismatches[:3]}"
|
||||
)
|
||||
|
||||
|
||||
def test_waveform_ppv_matches_export(idfw_export):
|
||||
header, _rows = idfw_export
|
||||
result = read_idf_file(IDFW)
|
||||
for channel, attr in (
|
||||
("Tran", "transverse_ips"),
|
||||
("Vert", "vertical_ips"),
|
||||
("Long", "longitudinal_ips"),
|
||||
):
|
||||
decoded = getattr(result.event.peaks, attr)
|
||||
assert decoded == pytest.approx(
|
||||
_header_float(header, f"{channel}PPV"), abs=5e-5
|
||||
), f"{channel} PPV disagrees with Thor's export"
|
||||
|
||||
|
||||
# ─── Histogram path shares the same scale ──────────────────────────────────
|
||||
|
||||
|
||||
def test_histogram_peaks_match_export():
|
||||
header, _rows = _parse_export(IDFH.with_suffix(".IDFH.csv"))
|
||||
result = read_idf_file(IDFH)
|
||||
assert result.intervals, "IDFH decoded no intervals"
|
||||
for channel, attr in (
|
||||
("Tran", "transverse_ips"),
|
||||
("Vert", "vertical_ips"),
|
||||
("Long", "longitudinal_ips"),
|
||||
):
|
||||
decoded = getattr(result.event.peaks, attr)
|
||||
expected = _header_float(header, f"{channel}PPV")
|
||||
# Histogram peaks are stored per-interval, so the export's PPV is
|
||||
# reproduced within one quantisation step rather than exactly.
|
||||
assert decoded == pytest.approx(expected, abs=2 * _GEO_LSB_IPS), (
|
||||
f"{channel} histogram peak {decoded} vs export {expected}"
|
||||
)
|
||||
|
||||
|
||||
# ─── Regressions found 2026-09-10 ──────────────────────────────────────────
|
||||
|
||||
IDFH_LONG = FIXTURES / "UM13981_20220207183102.IDFH" # 719 intervals
|
||||
IDFH_SENTINEL = FIXTURES / "UM13981_20221202063059.IDFH" # holds an unwritten slot
|
||||
IDFW_RAW16 = FIXTURES / "UM13981_20220207084555.IDFW" # segment 0 is MODE_RAW16
|
||||
|
||||
|
||||
def test_histogram_decodes_past_250_intervals():
|
||||
"""The segment validator must not require a zero counter high byte.
|
||||
|
||||
The interval counter is a uint16 cumulative index. Requiring its high
|
||||
byte to be zero rejected every segment past interval 255, capping each
|
||||
histogram at 250 intervals and truncating any run longer than ~4 hours —
|
||||
frequently discarding the part that held the peak.
|
||||
"""
|
||||
result = read_idf_file(IDFH_LONG)
|
||||
header, _rows = _parse_export(IDFH_LONG.with_suffix(".IDFH.csv"))
|
||||
expected = float(header["NumberOfIntervals"])
|
||||
assert len(result.intervals) == 719
|
||||
assert len(result.intervals) == pytest.approx(expected, abs=1.0)
|
||||
|
||||
|
||||
def test_histogram_ignores_unwritten_interval_slot():
|
||||
"""A never-written interval keeps its ±full-scale seed and must be dropped.
|
||||
|
||||
Counting it fabricates a 10.0 in/s peak on every channel, which then wins
|
||||
the max-over-intervals and poisons the whole file's PPV.
|
||||
"""
|
||||
header, _rows = _parse_export(IDFH_SENTINEL.with_suffix(".IDFH.csv"))
|
||||
result = read_idf_file(IDFH_SENTINEL)
|
||||
for channel, attr in (
|
||||
("Tran", "transverse_ips"),
|
||||
("Vert", "vertical_ips"),
|
||||
("Long", "longitudinal_ips"),
|
||||
):
|
||||
decoded = getattr(result.event.peaks, attr)
|
||||
assert decoded < 1.0, f"{channel} peak {decoded} looks like the ±FS seed"
|
||||
assert decoded == pytest.approx(
|
||||
_header_float(header, f"{channel}PPV"), abs=2 * _GEO_LSB_IPS
|
||||
)
|
||||
|
||||
|
||||
def test_waveform_raw16_segment_zero_is_decoded():
|
||||
"""Segment-0 records can be raw int16 (MODE_RAW16, 10-byte header).
|
||||
|
||||
That mode was absent from the dispatch, so the record fell through
|
||||
unhandled and the channel silently lost its first 512 samples.
|
||||
"""
|
||||
rows = _parse_export(IDFW_RAW16.with_suffix(".IDFW.csv"))[1]
|
||||
result = read_idf_file(IDFW_RAW16)
|
||||
for index, channel in enumerate(GEO_CHANNELS):
|
||||
decoded = result.samples[channel]
|
||||
assert len(decoded) == len(rows), f"{channel} lost segment 0"
|
||||
expected = [row[index] for row in rows]
|
||||
bad = sum(
|
||||
1 for c, v in zip(decoded, expected)
|
||||
if abs(geo_count_to_ips(c) - v) >= 5e-5
|
||||
)
|
||||
assert bad == 0, f"{channel}: {bad} samples differ from Thor's export"
|
||||
|
||||
|
||||
def test_body_offset_search_is_not_quadratic():
|
||||
"""The body scan must stay cheap enough for bulk ingest.
|
||||
|
||||
MODE_RAW16 is (0x00, 0x00), so scanning for candidate *preambles* treats
|
||||
every run of three zero bytes as a body start and trial-decodes each one
|
||||
(~0.5 s/file measured). The search anchors on record headers instead.
|
||||
"""
|
||||
import time
|
||||
|
||||
start = time.perf_counter()
|
||||
for _ in range(3):
|
||||
read_idf_file(IDFW_RAW16)
|
||||
elapsed = (time.perf_counter() - start) / 3
|
||||
assert elapsed < 0.15, f"body-offset search took {elapsed*1000:.0f} ms/file"
|
||||
|
||||
|
||||
# ─── Mic-disabled (3-channel) units, found 2026-09-10 ──────────────────────
|
||||
|
||||
IDFW_3CH = FIXTURES / "UM20147_20250531135901.IDFW" # body head below old floor
|
||||
IDFH_3CH = FIXTURES / "UM20147_20250330070110.IDFH" # 56-byte interval records
|
||||
|
||||
|
||||
def test_three_channel_waveform_decodes_all_geo_channels():
|
||||
"""A mic-disabled unit's shorter header moves the record chain head.
|
||||
|
||||
Its head sits at 0x0dba, below the old ``_BODY_SCAN_FLOOR`` of 0x0E00, so
|
||||
the scan could not see it and fell through to the *Vert* segment-0 record
|
||||
— decoding a body shifted one position around the channel rotation, which
|
||||
surfaced as Vert being exactly 512 samples short.
|
||||
"""
|
||||
rows = _parse_export(IDFW_3CH.with_suffix(".IDFW.csv"))[1]
|
||||
result = read_idf_file(IDFW_3CH)
|
||||
for index, channel in enumerate(GEO_CHANNELS):
|
||||
decoded = result.samples[channel]
|
||||
assert len(decoded) == len(rows), (
|
||||
f"{channel}: {len(decoded)} samples, export has {len(rows)}"
|
||||
)
|
||||
expected = [row[index] for row in rows]
|
||||
bad = sum(
|
||||
1 for c, v in zip(decoded, expected)
|
||||
if abs(geo_count_to_ips(c) - v) >= 5e-5
|
||||
)
|
||||
assert bad == 0, f"{channel}: {bad} samples differ from Thor's export"
|
||||
# Mic is genuinely absent on these units, not merely undecoded.
|
||||
assert not result.samples.get("MicL")
|
||||
|
||||
|
||||
def test_three_channel_histogram_uses_56_byte_intervals():
|
||||
"""Interval stride is 16 bytes per channel + an 8-byte tail, not a constant.
|
||||
|
||||
A mic-disabled unit packs 56-byte records, so assuming 72 read 7 intervals
|
||||
out of every 10-interval segment and then walked off alignment into
|
||||
garbage, which decoded as ~10 in/s peaks. The true count comes from the
|
||||
segment's cumulative interval counter.
|
||||
"""
|
||||
header, _rows = _parse_export(IDFH_3CH.with_suffix(".IDFH.csv"))
|
||||
result = read_idf_file(IDFH_3CH)
|
||||
expected_intervals = float(header["NumberOfIntervals"])
|
||||
assert len(result.intervals) == pytest.approx(expected_intervals, abs=1.0)
|
||||
assert {iv.n_channels for iv in result.intervals} == {3}
|
||||
for channel, attr in (
|
||||
("Tran", "transverse_ips"),
|
||||
("Vert", "vertical_ips"),
|
||||
("Long", "longitudinal_ips"),
|
||||
):
|
||||
decoded = getattr(result.event.peaks, attr)
|
||||
assert decoded < 1.0, f"{channel} peak {decoded} looks like walked-off garbage"
|
||||
assert decoded == pytest.approx(
|
||||
_header_float(header, f"{channel}PPV"), rel=0.02
|
||||
)
|
||||
|
||||
|
||||
# ─── `40 NN` blocks with NN > 8, verified 2026-09-11 ───────────────────────
|
||||
|
||||
IDFW_WIDE40 = FIXTURES / "UM12947_20250806134504.IDFW"
|
||||
|
||||
|
||||
def test_wide_forty_nn_block_does_not_truncate_channels():
|
||||
"""Loud events use `40 NN` blocks with NN well above the old cap of 8.
|
||||
|
||||
``data_block_len()`` rejected NN > 0x08, which halted the block walk
|
||||
part-way through a record. The walker stops at the first unrecognised
|
||||
tag instead of raising, so this surfaced as silently short channels —
|
||||
here Tran 1812 / Vert 2132 / Long 2324 where the export has 2324 for all
|
||||
three. The affected files use NN of 12, 16, 20 ... up to 196.
|
||||
"""
|
||||
rows = _parse_export(IDFW_WIDE40.with_suffix(".IDFW.csv"))[1]
|
||||
result = read_idf_file(IDFW_WIDE40)
|
||||
for index, channel in enumerate(GEO_CHANNELS):
|
||||
decoded = result.samples[channel]
|
||||
assert len(decoded) == len(rows), (
|
||||
f"{channel}: {len(decoded)} samples, export has {len(rows)}"
|
||||
)
|
||||
expected = [row[index] for row in rows]
|
||||
bad = sum(
|
||||
1 for c, v in zip(decoded, expected)
|
||||
if abs(geo_count_to_ips(c) - v) >= 5e-5
|
||||
)
|
||||
assert bad == 0, f"{channel}: {bad} samples differ from Thor's export"
|
||||
@@ -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"
|
||||
@@ -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"") == {}
|
||||
@@ -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():
|
||||
|
||||
@@ -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)
|
||||
@@ -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])
|
||||
Reference in New Issue
Block a user