diff --git a/CHANGELOG.md b/CHANGELOG.md index 6b4fad0..d0ec631 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,7 +4,76 @@ All notable changes to seismo-relay are documented here. --- -## [Unreleased] +## v0.29.0 — 2026-09-04 + +First release to reach prod since **v0.27.0**, so it ships **both** the +`false_trigger_reason` column below *and* the v0.28.0 offset (DC-baseline) +detector: v0.28.0 was version-bumped in-tree (`TOOL_VERSION`, CHANGELOG) but +never tagged or deployed, so 0.29.0 is the first build to carry either to prod. +Pairs with Terra-View ≥ 0.24.0. The `false_trigger_reason` column auto-migrates +on startup; the offset detector still needs the shape backfill on the prod store +(`scripts/backfill_event_shape.py`) to populate `shape_offset*` on existing rows. + +### Added +- **`events.false_trigger_reason` — optional FT cause.** A nullable `TEXT` + column recording *why* an event is a false trigger (e.g. `"offset"`), as a + subtype of the FT flag: setting a reason via the sidecar review PATCH implies + `false_trigger=1`, and the reason is cleared whenever FT ends up 0 + (confirm-real, clear-FT, `set_false_trigger(false)`). `propagate_review_to_twins` + carries the reason to the histogram/waveform twin alongside the flag. + Auto-migrated (`_SCHEMA` + `_migrate` ADD COLUMN — not the Migration-1 + rebuild); exposed via `/db/events`. Terra-View surfaces it as a manual + "Flag as offset" action + an `FT · offset` badge. + +### Fixed +- **BlastMate serials — the family prefix is read from the file, not guessed.** + The Blastware filename encodes only the serial *number* (`L895…` → 10895); + the two-letter prefix is not in it. `waveform_store` synthesised `"BE"`, so + an imported **BlastMate** (serials `BA…`) was filed under a MiniMate Plus + serial that does not exist — silently, and Terra-View read it straight + through. `save_imported_bw` now resolves serial as hint → file body → + filename guess, via a new `_serial_from_bw_bytes` that accepts a candidate + only when its numeric part matches the filename. `client._decode_0a_partial_header` + likewise matched a literal `b"BE"` in monitor-log partial records; on a + BlastMate that returned −1 and skipped the whole block, losing the **geo + threshold** along with the serial. It now matches any two-letter prefix and + requires the NUL terminator — stricter than the search it replaces. + + BlastMate is the MiniMate Plus's larger Series III sibling and its files are + byte-compatible: all 1,493 in the DL2 archive decode through the existing + codec at 100%, same four channels. **The serial string was the only thing + blocking BlastMate support in SFM.** Four archive units were affected — + BA9229, BA10060, BA10895, BA15957. + + **No backfill and no `TOOL_VERSION` bump**: this changes which serial an + *import* is filed under, not any decoded value, so existing sidecars and + `.h5` files are untouched. **No migration either** — prod holds no BlastMate + events (the archive's BA units last recorded 2018-10 through 2023-11; the + prod backfill reaches back only to ~May 2025). + +--- + +## v0.28.0 — 2026-09-02 + +**Offset (DC-baseline) false-trigger detector.** Productionizes the validated +pre-trigger detector: a geophone event whose baseline sits off zero and stays +flat across the record (sensor bumped / settled / drifted) is now flagged and +surfaced in Terra-View as an `offset` false-trigger reason — catching offsets the +crest/near-peak spike rule misses (an offset is low-crest and flat). + +### Added +- `shape_metrics.offset_from_samples` / `offset_from_h5`: per geophone channel, + `|median(pre-trigger)| ≥ 0.025 in/s` AND `pre/mid/end spread ≤ 0.02` → offset; + the consistency test rejects transients (a real event moves one third). Reads + the `.h5` samples + the `pretrig_samples` attr, range-aware via the in/s float + samples. Constants `OFFSET_FLOOR` / `OFFSET_MAX_SPREAD` are tunable. +- `events.shape_offset` / `shape_offset_axis` / `shape_offset_pre` / + `shape_offset_spread` columns (auto-migrated: `_SCHEMA` + the `_migrate` + ADD COLUMN loop), computed at all three ingest paths and by + `backfill_event_shape.py`, exposed via `/db/events`. + +Requires the shape/offset backfill on the prod store to populate existing events: +`python scripts/backfill_event_shape.py --db-path … --store-root …`. --- diff --git a/CLAUDE.md b/CLAUDE.md index 71616d9..6d4c6ea 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -2,7 +2,7 @@ Ground-up Python replacement for **Blastware**, Instantel's Windows-only software for managing MiniMate Plus seismographs. Connects over direct RS-232 or cellular modem -(Sierra Wireless RV50 / RV55). Current version: **v0.27.0**. +(Sierra Wireless RV50 / RV55). Current version: **v0.29.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 diff --git a/README.md b/README.md index 4cb9676..dfe95e0 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -# seismo-relay `v0.27.0` +# seismo-relay `v0.29.0` A ground-up replacement for **Blastware** — Instantel's aging Windows-only software for managing seismographs. Supports both the **MiniMate Plus diff --git a/docs/offset_investigation.md b/docs/offset_investigation.md index 20da227..d69f0a1 100644 --- a/docs/offset_investigation.md +++ b/docs/offset_investigation.md @@ -58,6 +58,14 @@ Companion material: sensor-check failures. Do not try to use it as a screen. - **Cause is still unsettled.** Instantel's autozero fixes the minority of cases; the rest are hardware. We cannot yet tell which is which remotely. +- **The histogram corpus (63,535 files, 9.7x the waveforms) is now scanned too** — + see §8b. It independently confirms BE18438 and BE9558 with a clean 2.5x + separation, but detects only **2 of the 5** confirmed units, cannot attribute a + channel, and resolves time to ~a month. **A negative histogram result is not + evidence of health** — DC leakage into the interval peak varies 45x between units. +- **`offset_scan3.py` has a label defect** (§8b): its spread gate discards 18.8% of + high-|pre| rows onto units currently counted as clean. Re-cut before quoting any + precision number again. - **Best open lead:** `SUB 0x0E` (channel sensor data, 8 channels × 10 bytes, unimplemented) may carry the very numbers Instantel says to check against **2027–2069**. Untested. @@ -146,8 +154,8 @@ is the discriminator — and it is what the field experience predicts. | runs of >=3 consecutive | — | 29 | | runs of 1-2 events (noise) | — | 69 | -Units with a sustained pedestal: **BE9558, BE10895, BE11007, BE11529, BE12599, -BE13117, BE18003, BE18438**. BE10895 and BE18003 were invisible to v1. +Units with a sustained pedestal: **BE9558, BA10895, BE11007, BE11529, BE12599, +BE13117, BE18003, BE18438**. BA10895 and BE18003 were invisible to v1. **The affected channel is most often Vert**, which v1 got wrong — it named whichever axis had the largest peak. BE13117 and BE18438 are both Vert faults. @@ -218,7 +226,7 @@ signal from a tuned one: **BE9558, BE11529, BE12599, BE13117, BE18438.** -Unchanged across a 2x threshold range. BE11007 and BE10895 drop out — the +Unchanged across a 2x threshold range. BE11007 and BA10895 drop out — the spread test identifies them as transients, not pedestals. The 11% headline happens to match v1's, but the reasoning and the unit list @@ -495,6 +503,459 @@ doubled two reported figures before it was caught. --- +## 8b. The histogram corpus — the other 90% of the archive (2026-09-04) + +Every result above §8 comes from **waveform** files. `offset_scan3.py` filters on +`\.[A-Za-z0-9]{2}0[Ww]$`, so the corpus it scanned is 6,577 unique binaries. The +archive also holds **63,535 unique histograms** — 9.7x more files — which the +pre-trigger method cannot touch, because a histogram carries no samples: only a +per-interval, per-channel peak and half-period. + +`scratch/offset_hist_scan.py` scans them. **63,505 of 63,535 decoded (99.95%), +43 units, 77.9M intervals.** Two of the 45 units have no histograms at all. +Output: `/home/serversdown/dl2-archive/offset_hist.csv` (190,515 channel-rows). + +### The premise, and how far it actually holds + +A histogram file is hours of continuous monitoring, so most of its intervals are +definitionally quiet, and a channel parked off zero cannot report a peak below +its own displacement. The signal is real — two within-unit contrasts, siblings +unmoved in both: + +| unit | channel | in-episode floor | outside | waveform \|pre\| same window | +|---|---|---|---|---| +| BE18438 | Vert | 0.0350 | 0.0050 | +0.18 .. +0.37 | +| BE12599 | Tran | 0.0250 | 0.0050 | +0.03 .. +0.49 | + +But the **leakage from a waveform pedestal into the histogram floor is bimodal, +not merely partial**: measured ratio ~0.9 on BE18438 Vert, ~0.7 on BE9558, +**~0.02 on BE12599** — two orders of magnitude on one instrument. The device +evidently measures each interval peak against a running baseline, and how much +DC survives that varies per unit. **Consequence: a negative histogram result +carries almost no information.** Do not read "clean in the histograms" as clean. + +### The detector that survived + + dmin(file, ch) = min[ch] - min over the other two geo channels, SAME file + gates (both hard): n_intervals >= 60 AND mic_p5 <= 5 raw counts + day statistic: median of dmin over that day's qualifying files + flag day at dmin >= 0.020 in/s (4 A/D counts) + episode at >= 3 CONSECUTIVE observed days + +**Result: BE18438|Vert, BE9558|Tran, BE9558|Long.** Threshold-insensitive — +the journal's own test for a real signal against a tuned one — and this is the +first operating point in the investigation that passes it cleanly. The identical +answer holds across: statistic `min` or `p5`; length gate 10/30/60/120/300; mic +gate 3/5/8/10; threshold 0.015–0.035 (a 2.3x span); persistence K = 2,3,4,5,7. + +Separation, ranked by highest floor sustained over 3 consecutive gated days +across all 135 unit-channels: + +| unit-channel | best3 | +|---|---| +| BE18438 Vert | 0.1650 | +| BE9558 Long | 0.0350 | +| BE9558 Tran | 0.0250 | +| *(2.5x gap)* | | +| BE7145 Tran | 0.0100 | +| entire rest of fleet | <= 0.0050 (one quantisation count) | + +Day-level false alarm: **37 of 99,432 gated unit-channel-days = 0.037%.** + +### What it does NOT do — read this before trusting it + +- **It finds 2 of the 5 confirmed units, not 5.** The site-quiet gate is what + makes it work and it is also what costs BE11529 and BE12599. BE11529's + four-day single-axis ramp (Tran 0.025 -> 0.055, both siblings pinned at 0.005) + is the most offset-shaped thing in the corpus outside the two detections, and + the gate discards it. +- **The positive class is two units.** Every threshold here is fitted to + BE18438 and BE9558, which contribute 22 of the 37 flagged days in the entire + corpus. No cross-validation is possible at n=2. +- **Per-channel attribution is NOT established.** Rotating the three geo channel + labels within each file — preserving every value, file and day, destroying + only channel identity — reproduces the episode *count* with p = 0.769 and the + label agreement at p = 0.038–0.077. Report a **unit and a window**; do not + name a geophone axis on the strength of this detector alone. +- **Timing resolution is ~1 month, not ~1 day.** A 30-day label shift still + scores 2 of 9 episode hits; the signal dies only past ~60 days. The day-level + series look far crisper than they are. +- **Ground truth here is a sibling detector, not a service record.** Agreement + between the two corpora is corroboration of a shared method. Nothing in this + section has been checked against an actual repair, calibration or RMA. + +### Dead ends — keep these dead + +- **Absolute floor (min / p1 / p5 / p10 / p25, thresholded alone) — RETIRED.** + Not fleet-comparable and mostly not about the channel. Scoring each cell using + *only the other two channels* — a statistic containing zero information about + the suspect channel — reaches AUC 0.746 against the same labels, versus 0.872 + for the absolute floor itself. **66% of its apparent discrimination is "that + day was noisy at that site."** Interval size alone moves its p99 7x (0.0350 at + 1 min vs 0.0050 at 2 s). And of all files with any channel above 0.025, 56.5% + have **all three** channels above it — common-mode, i.e. the wrong physics. +- **Zero-fraction — STRUCTURALLY IMPOSSIBLE, not merely weak.** The device never + reports a zero histogram interval peak. The value is a max over hundreds of + samples of a channel that always carries at least 1 count of noise, so it is + clamped at 1 A/D count (0.005 in/s). There is no zero to count. +- **Interval size, sample rate, geo range, firmware — refuted as confounds for + the differential.** All four are *file-level scalars*: they move all three geo + channels together, so they cannot produce a single-channel lift and the + within-file differential is immune to them by construction. Geo range is + identical across the three geo channels in **63,535 of 63,535** binaries. + (Interval size remains fatal to the *absolute*-floor version, above.) + +### Two findings that are independent of the histogram detector + +**1. `offset_scan3.py`'s `spread <= 0.02` gate is discarding real signal.** +It rejects **113 of the 600 channel-rows with |pre| >= 0.025 (18.8%)**, and the +rejections are not random — 92 of them fall across 41 unit-channels currently +labelled NEGATIVE. Four would become sustained positives under an +amplitude-only >=3-consecutive rule: **BE12599|Long (run of 8), BE18003|Vert +(4), BA10895|Vert (3), BE12844|Tran (3).** Until this is re-cut, the fleet label +is **three-state — POSITIVE / NEGATIVE / SPREAD-REJECTED(unknown)** — and the +third state should be excluded from both TP and FP counts rather than silently +scored as healthy. Every precision figure computed against the two-state label, +in this section and in §3, is affected. + +**2. The waveform corpus sees ~7% of the days a unit was deployed.** 2,627 +(unit, day) observations against the histogram corpus's 35,105 — 13.4x — with a +per-unit median ratio of 0.070. BE12599, a confirmed unit, is waveform-observed +on 39 of its 1,666 histogram-observed days (**2.3%**). Any statement of the form +"the fault was absent before date X" that rests on waveform coverage alone is +much weaker than its event count suggests. + +### BA10895 — reclassified (see also §4) + +Previously dismissed as a transient. The histogram record shows its **Vert** +quiet-minute floor at 0.005 on 62/62 qualifying files from 2023-07-07, then +0.010–0.015 on 48/58 files from 2023-08-03 to 08-27, while Tran moves on 2/58 +and Long on 9/58 and the site mic floor never leaves 1–3 counts. Independently, +**42 of its 85 waveform events (49.4%) are single-axis-dominant** — one geo peak +>= 10x both siblings and >= 0.05 in/s — the **highest rate in the 45-unit +fleet** (BE13117 36.1%, BE18438 29.4%), and **100% of it on Vert**. Vert +excursions of 0.1–1.5 in/s with Tran/Long at 0.005–0.035 are not ground motion. + +This is a genuine Vert-channel hardware fault, but **not the classic pedestal** — +the differential is only one A/D count. Caveat: its entire histogram record is a +single 52-day deployment ending 2023-08-27, so nothing says whether it +persisted, was serviced, or resolved. + +The other six marginal units — BE11007, BE17354, BE18004, BE18104, BE9557, +BE18003 — are **clean**. All seven cap at +0.005 to +0.007 (one A/D count) +lifetime under the quiet-site gate, against +0.175 for BE18438 Vert and +0.062 +for BE9558 Long. Three individual waveform flags fall in windows with **zero** +histogram coverage and are NO-DATA, not clean: BE18004|Tran 2024-10-16, +BE9557|Tran 2021-06-28, BE9557|Vert 2025-06-12. + +### Still open in this section + +- **The 11 thin-coverage units were not screened** (BE10202, BE11462, BE13779, + BE15760, BA15957, BE16754, BE16758, BE8081, BE8626, BA9229, BE9887 — each + under 20 waveform events, several with hundreds of histograms). This is the + population most likely to hold a previously unknown offset, and it is the one + slice of the plan that did not run. BE11462 was incidentally scored clean by + the full-archive pass; BE10202 has no histogram files at all. +- **No completeness audit was run** over the above. +- Re-cutting the ground truth three-state (finding 1) and re-scoring everything + against it. + +--- + +## 8c. Mechanism — five hypotheses tested, all dead (2026-09-06) + +**The mechanism is still unknown.** Five campaigns, ~105 effectively independent +tests, seven nominally significant results against **5.2 expected by chance** +under a global null. Every one died to its own confound analysis. What the +campaign bought is a set of *shape constraints* and a long list of dead ends. + +### ⚠ Two things retracted from this journal + +**1. "Polarity is perfectly consistent — 11 of 11, zero mixed cases."** That is +a **tautology of the spread gate**, not a property of the fault. `spread <= 0.02` +requires pre/mid/end to agree, which forces one sign. Amplitude-only at the same +0.025 threshold: **12 of 53 unit-channels are mixed**, including BE18438|Vert +(88+/1−) and BE9558|Vert (1+/35−). Withdrawn. + +**2. "5 of 45 units, unchanged across a 2x threshold range."** The +threshold-insensitivity is also a property of the gate. Amplitude-only gives +**9 units at 0.020, 8 at 0.025** (adding BA10895, BE12844, BE18003), 5 at 0.040. +The fleet is **8–9 units, not 5**. + +**3. "Persistent — it stays until the geophone is serviced."** Weakened, not +withdrawn. There are **23 recoveries after runs of >=3 flagged events, median +gap 6.03 days**, three inside ten minutes. BE18438|Vert reads `pre=mid=end= ++0.0000` on 2026-02-10, +0.185→+0.370 across 02-25/26, and `+0.0000` again on +2026-03-22 — identical Project, Seis Loc, calibration date, geo range and +trigger throughout. The one thing that cannot be excluded is a **field +autozero**: it is a button sequence at the unit and writes nothing into the +event header. So "persistent" may be "persistent unless somebody pressed the +buttons," and the archive cannot tell those apart. + +### The one positive finding: onset is a RAMP, minutes to hours + +Both onsets resolvable at minute cadence are ramps. **BE18438|Vert, +2026-02-20** — the histogram corpus collapses a 14 d 21 h waveform bracket to +**one minute**: + +``` +~14,200 consecutive quiet minutes at 0.000–0.005 (ten full daily files) +09:32 +0.005 09:39 +0.045 10:20 +0.125 16:00 +0.165 +09:33 +0.010 09:42 +0.070 13:13 +0.150 20:17 +0.185 plateau +``` + +**50% of the excursion in 7 minutes**, the rest asymptotic over ~10 h, **>=25 +distinct one-minute intermediates**. Validated **75/75** against Blastware's own +ASCII export. Its 2025-11-15 onset is the same shape over 2.7 h. BE13117 stage B +is a 91-minute monotone rise, +0.035 → +1.745 in/s over ~40 samples. + +**This kills both poles of the original dichotomy** (journal Q1): not an +instantaneous latched step (a bad autozero, a stuck trim-DAC), and not slow +component degradation over days or weeks. + +⚠ It rests on **2 of 45 instruments**. Clopper-Pearson on 4/4 resolved onsets +gives 95% CI [0.40, 1.00] — a mixed population with up to 60% true steps is not +excluded. BE13117 has zero paired ASCII, so its ramp rests on our decoder alone. + +### The methodological corollary — more important than the finding + +**A waveform-only bracket manufactures the appearance of a step, and the spread +gate is blind to onsets by construction.** + +The offset is what fires the trigger, so no waveform event can exist until the +ramp has nearly reached the trigger level. BE18438's first event of each episode +sits at 0.280 against a 0.300 trigger, and 0.185 against 0.200. At daily cadence +against a 3 h ramp, P(catching an intermediate) = **0.125**. + +And `spread <= 0.02` rejects any record in which the floor is *moving* — which +is exactly what an onset is. **The gate rejected the very BE18438 record where +the ramp is visible.** If the operational goal is catching a fault early, before +the unit floods the store with junk events, the current detector is the wrong +shape for the job. + +### The surviving shape + +An **electrical, reversible, two-time-constant settling process** (~10 min and +~hours), saturating at a ceiling, with occasional sub-3-minute discrete jumps +superposed (BE18438 2026-02-26: 13:24 pre +0.180 / mid +0.240 / end +0.255 → +13:27 +0.325, identical metadata). That is the signature of a **bias or leakage +path charging a high-impedance node** — the class of fault Instantel's autozero +recovers ~10% of the time, and what the X1/X8 gains measure. + +**It is a shape constraint, not a mechanism. Do not write it up as one.** + +### Dead — with the evidence, so none of this is re-derived + +| Killed | Evidence | +|---|---| +| **Latched step at onset** | >=25 one-minute intermediates over ~10 h, ASCII-validated. Direct observation, not a test. | +| **Slow degradation over days/weeks** | Same observation — bulk of the excursion in 7 min to 2.7 h. | +| **Thermal driving of pedestal magnitude** | BE13117, 365-count pedestal, n=128: full-day modulation **−0.42% ± 0.42%**, 95% CI [−1.25%, +0.40%]. Healthy-fleet seasonal zero drift totals **~0.3 A/D counts** — 15x to 1200x too small. Best-powered result in the campaign. | +| **Ground-motion shock** | 30-day window-max percentile ranks 0.03/0.98/0.15/0.01/0.68/0.24, median **0.194** against a null of 0.5. **0 of 7 events >=9 in/s** was followed by an onset within 30 d. BE12599 hit 10.220 in/s (2023-11) and 10.005 (2025-04) and did not onset until 2026-08-14. | +| **Handling / redeployment** | **0 of 9** onsets had a Project/Client/Seis Loc change. Widened to 30 d: 2 observed vs 4.90 expected, P(X>=2)=0.995 — *depleted*, the wrong direction. The apparent gap effect (p=0.035) died on histogram coverage: BE18438's "59.7-day gap" contains 122 histogram files; true silence 0.52 d. | +| **Mechanical resonance / damping change** | BE18438|Vert at a 64-count pedestal (3x outside Instantel's ±21): ΔTest-Freq **CI [−0.090, +0.021]** against 0.127 Hz for a real calibration. Block permutation p=0.658. | +| **Accumulated-duty threshold** | ~4 clean units logged more monitoring than the largest positive onset dose; BE18193 logged **13.45M intervals, 6.2x**. A counterexample — no power argument weakens it. | +| **Firmware** | **14,338 of 14,340** exports read `V 10.72-8.17`. A constant cannot explain a variable. | +| **Unit age** | Serial rank-sum 118.0 vs null 115.0, p=0.549; unchanged on the 8-unit re-cut (p=0.586). Serial is a poor age proxy anyway (Spearman +0.113 against archive entry). | +| **Strong seasonal clustering** | 25 onsets, exposure-weighted permutation **p=0.59**. Excludes >=75%-in-one-season only; a 2x seasonal hazard is *not* excluded. | + +Also retire two overstated bounds. H6's dose-response exclusion "|r| > 0.03" is +a **10x overstatement** once clustering is corrected — the honest bound is +|r| > 0.1–0.3, so a real r=0.2 is not excluded. And **any statistic quoted +per-event**: 512 flagged channel-events collapse to **4.9 effective independent +observations** (unequal-cluster design effect 104.6 at ICC=1), and **55% of the +flagged corpus is one instrument on two calendar days** (BE13117, 2023-05-03/04). + +### Power — read every negative in this section as bounded + +Fisher exact, 5 positives of 45, one-sided α=0.05, exposure a third of the fleet: + +| relative risk | power | +|---|---| +| 1.5 | 0.059 | +| 2 | 0.112 | +| 3 | 0.231 | +| 6 | 0.497 | +| 15 | 0.753 | + +80% power needs **RR ≈ 13–20**. Even a *perfect* split reaches p<0.05 only if +the exposed group is <=25 of 45 units. **This archive can detect only +near-deterministic unit-level causes.** Every negative above excludes a strong +effect, not a real one. + +### What this archive can NEVER answer + +- **The A/D zero and the X1/X8 gains.** The 2027–2069 numbers appear in no file, + header or decoded record. They exist only on a live device behind `SUB 0x0E`. + Q1 is structurally unanswerable from data. +- **Unit-level vs component-level cause.** **Zero of 14,340** exports carry a + geophone or sensor serial. Q4 is dead — there is no way to know whether the + same physical geophone came back after service. +- **Service history.** The only service-adjacent field is `Calibration: ` + — 30 distinct dates fleet-wide, none before 2023, ASCII corpus entirely + 2025–26. BE9558's 2020 and BE13117's 2023 episodes have no calibration record. +- **Temperature.** Zero exports carry it. Battery Level is a verified coarse + thermometer (+0.204 V winter over summer, 20/20 unit-years, p=9.5e−7, matching + lead-acid tempco) but quantised at 0.1 V ≈ 10 °C — useless within a day. The + archive can *bound* thermal; it can never *test* it. +- **BE13117 specifically** — 55% of the flagged corpus, the largest pedestal at + 1.92 in/s, **zero** ASCII exports, histogram record ending eight months before + its episode. The most informative case in the archive is permanently outside + every metadata test. +- **The mild-offset rate**, and therefore the base rate's denominator. Event + files only see offsets large enough to dominate the trace. + +### The experiment to run — `SUB 0x0E`, one afternoon + +Point Blastware at `bridges/ach_mitm.py` and run **Unit Channel Test** against +(1) a faulting unit, (2) a known-good control, (3) the same unit before and +after an autozero. BW's sequence is `0x0E x8 → 0x98 x2 → 0x0E x8`, the second +pass carrying live ADC. Eight 10-byte payloads with expected values near 2048 is +a very constrained puzzle. + +- **Proves:** whether the X1/X8 gains are readable over the wire, and whether + the fault sits at or upstream of the ADC zero reference. Gains walk out of + 2027–2069 with the pedestal → the fault *is* the zero reference, Q1 answered. + Gains hold while the trace moves → the fault is downstream, look at the front + end. +- **§8c hands it a falsifiable time course:** poll at ~1-minute cadence and the + numbers should **ramp over minutes-to-hours, not step**. If they step while + the trace ramps, the two are decoupled. +- **Payoff:** converts the 10%/90% ship-it-or-not gamble into a decision made + before packing a box, remotely, for the whole fleet. +- ⚠ In the MITM topology filenames are reversed — `raw_s3_*.bin` holds + Blastware's bytes. + +**Second: swap the geophone** between a faulted base and a healthy one. Fault +follows the sensor → element or cable. Fault stays with the base → front-end +board. One afternoon, zero code, and it settles the one question the archive is +permanently blind to. + +**Third: log a faulting unit for 72 h untouched.** Every recovery we have is +confounded by a possible field autozero. A shelf and a logger settles whether +the fault genuinely self-reverses. + +**Fourth, free: re-cut the fleet label** — drop the spread gate, re-score +amplitude-only, screen the 11 unscreened thin-coverage units. Might reach 9–10 +positives. Be honest about the gain: power against "older half carries 3x the +hazard" rises only 0.23 → 0.30. + +**Highest-value item overall, and not an experiment: the RMA/repair records.** +Which unit went back, when, what was done (autozero vs geophone replaced vs +board), and the geophone serial fitted. "Same channel after a documented +geophone *replacement*" is component-level-negative in one observation. + +--- + +### 8d. The non-motion test — Brian's "it doesn't cross zero" (2026-09-07) + +Looking at BE12599's 2026-08-09 event, Brian noted it reports no ZC frequency +**because the trace never crosses zero**. That observation is the best detector +in this investigation, and it comes from physics rather than a threshold. + +A geophone is a velocity sensor with no DC response, so its output over a record +must integrate to ~zero — the ground does not relocate. Real motion therefore +sits roughly half below zero. Anything electrical is one-sided. + + mp = |mean| / peak ~0 for motion, ~1 for a fault + frac_neg = share of samples < 0 + +`scratch/nonmotion_scan.py`, all 6,577 waveforms, 19,731 channel-rows. +Restricted to peak >= 0.05 in/s (n = 12,068), the distribution is **bimodal +with an empty middle**: + +| mp band | channel-events | +|---|---| +| 0.0–0.1 | 11,384 | +| 0.1–0.2 | 293 | +| **0.15–0.85 (dead zone)** | **131 = 1.09%** | +| 0.9–1.0 | 278 | + +At `mp >= 0.8` with >=3 events it returns **exactly the five confirmed units** — +BE9558, BE11529, BE12599, BE13117, BE18438 — stable from 0.5 to 0.9. Two +detectors on entirely different principles agreeing on the unit list is the +strongest corroboration that list has. + +**BE11007 is settled: NOT an offset.** It reaches mp 0.75–0.89, but with +`frac_neg = 0.99` at peaks of **7.4–9.4 in/s** — parked *negative* during a +near-full-scale blast. §4's guess was right. `mp` alone cannot separate a +pedestal from a large one-sided blast; pair it with a peak ceiling or with +sign-consistency across events. + +⚠ **Not a rediscovery of the retracted v1 detector.** v1 scored only the +largest-peak axis and used the mean as a *baseline estimator* where the median +was required. Here the mean is the signal itself, per channel — that is what the +physics licenses. + +**Correction to §8c.** That section says the spread gate is "blind to onsets by +construction." Too strong: of 87 BE18438|Vert events at mp >= 0.5 the gate +rejected **one** — the transitional record. It does not lose onsets +systematically; it loses the transition specifically. + +### 8e. BE12599 — a connector, not a geophone (2026-09-07) + +Waveform shapes across its August episode, measured rather than eyeballed: + +| date | channel | shape | +|---|---|---| +| Aug 09 05:29 | Long | **unipolar +**, 0/2304 samples below zero, decay tau **26 ms** | +| Aug 09 05:35 | Long | unipolar +, 3 spikes at irregular gaps (744, 1032 ms), tau **38 ms** | +| Aug 14 05:00 | Long | single lobe, bipolar, tau **118 ms** | +| Aug 17–23 | Tran | **flat DC pedestal**, sd/level 0.015–0.020, 0 zero crossings | + +**Unipolar impulses with an RC tail are not mechanical.** Fast rise, exponential +decay, one polarity, irregular timing — that is charge dumped into a +capacitively-coupled input and draining through the input resistance. The +progression 26 ms -> 118 ms -> never recovers, over 14 days, is a leakage path +worsening. + +**And the fault moved channels** — Long on Aug 9/14, Tran on Aug 17–23, Long +again on Aug 21 (1.065 in/s) while Tran held its pedestal. Vert stayed clean +throughout. **A failing geophone element cannot hop channels. A connector can.** + +That single fact explains what had been puzzling: +- **The sensor self-check keeps passing** (7.4/7.5/7.6 Hz, ratios 3.6–4.2, all + four channels Passed, on the very events where Long throws 0.5 in/s spikes). + The swing test drives the element; the element is fine. The fault is in the + wiring to it. +- **Why Instantel's autozero fixes only ~10%** — it cannot fix a connector. +- **Why onset "ramps" over minutes to hours** — contact resistance drifting. + +All seven Aug 17–23 events are stamped **05:00:14**, the same second, and their +filename extensions run `8E → WE → KE → 8E → WE → KE → 8E` — the documented +3-day cycle for a fixed daily time. Clock-scheduled, not physically triggered: +the modem powers up, draws a surge, and a marginal connection responds. + +**Field action: inspect and photograph the geophone connector BEFORE reseating +anything** — an intermittent contact clears the moment it is disturbed. + +⚠ Scoped to BE12599. BE18438's onset was a smooth 7-minute ramp with no spikes, +which looks like a different failure mode wearing the same signature. + +--- + +### ⚠ Serial prefixes — four of these units are BlastMates, not MiniMates + +Corrected 2026-09-06, after Brian queried "BA10895?" against a report that +said BE10895. He was right. The BW filename encodes the serial **number +only** — `L895` -> 10895 — and every offset scanner synthesised the family +prefix as `"BE"`. Four of the 43 archive units are **BA** (BlastMate, the +MiniMate Plus's bigger sibling; same Series III, byte-identical data): + +**BA9229, BA10060, BA10895, BA15957.** + +Read off the file bodies, which carry the serial verbatim. No analysis +changed — grouping was always on the numeric part, and no unit number maps +to two serials — but every earlier reference to "BE10895" and the other +three is a label error and has been corrected throughout this document. + +The same assumption was live in two production sites and is fixed +(`sfm/waveform_store.py`, `minimateplus/client.py`): the store would have +filed a BlastMate under a unit that does not exist, and the monitor-log +decoder lost the geo threshold along with the serial. See commit `9ceff65`. + +--- + ## 9. Chronology | date | event | @@ -512,3 +973,15 @@ doubled two reported figures before it was caught. | 2026-08-28 | Bimodality established; sensor check proven **blind** to offsets; `SUB 0x0E` identified as the best open lead. | | 2026-08-28 | **v1 detector retracted.** Brian challenged the "come and go" finding against field experience. Two flaws found: dominant-axis-only scoring and mean-instead-of-median. Corrected detector shows persistent pedestals on **8 of 45 units**, and the gaps are service windows. | | 2026-08-28 | **Detector v3 (Brian's method):** pre-trigger floor + pre/mid/end consistency. Healthy channels proven to sit at 0.000 +/-1 unit (94.5%), confirming no decoder zero-point bias. Final: **5 of 45 units (11%)**, threshold-insensitive. | +| 2026-09-04 | **Histogram corpus scanned** — 63,505 of 63,535 files, 43 units, 77.9M intervals (9.7x the waveform corpus). `scratch/offset_hist_scan.py`. | +| 2026-09-04 | Absolute-floor statistic **retired**: 66% of its discrimination is a day/site confound (other-channels-only AUC 0.746 vs 0.872). Zero-fraction shown **structurally impossible** — the device clamps every interval peak at >= 1 count. | +| 2026-09-04 | Site-quiet-gated cross-channel differential established: **BE18438 Vert, BE9558 Tran+Long**, threshold-insensitive over a 2.3x span. Finds only **2 of the 5** confirmed units — leakage into the histogram floor is bimodal (0.9 to 0.02), so a negative result carries almost no information. Per-channel attribution **not** established (channel-scramble p = 0.769). | +| 2026-09-04 | **BA10895 reclassified** from transient to a genuine Vert fault of a different subtype — 49.4% single-axis-dominant events, the highest in the fleet, 100% on Vert. The other six marginal units are clean. | +| 2026-09-04 | **Defect found in `offset_scan3.py`**: its `spread <= 0.02` gate discards 18.8% of rows with \|pre\| >= 0.025, concentrated on 41 negative unit-channels; 4 would be sustained positives without it. The fleet label is three-state, not two. | +| 2026-09-06 | **Four units relabelled BA, not BE** — BA9229, BA10060, BA10895, BA15957 are BlastMates. The BW filename carries only the serial number; the family prefix must be read from the file body. Fixed in the scanners and in two production sites. | +| 2026-09-06 | **Mechanism campaign — five hypotheses, all dead.** Thermal, ground-motion shock, handling/redeployment, accumulated duty, unit age, firmware and a mechanical element fault are each refuted or bounded. 7 nominally significant results against 5.2 expected by chance. | +| 2026-09-06 | **Onset is a RAMP of minutes-to-hours, not a step** — BE18438 Vert resolved to one-minute cadence, 50% of the excursion in 7 min, >=25 intermediates, ASCII-validated 75/75. Kills both a latched digital step AND slow component degradation. Surviving shape: a reversible two-time-constant settling process — a bias/leakage path charging a high-impedance node. | +| 2026-09-06 | **Polarity consistency RETRACTED** (a tautology of the spread gate; amplitude-only gives 12 of 53 unit-channels mixed) and the fleet **re-cut to 8–9 units, not 5**. "Persistent until serviced" weakened: 23 recoveries, median gap 6 days — though a field autozero cannot be excluded. | +| 2026-09-06 | The spread gate is **blind to onsets by construction** — it rejects a moving floor, which is what an onset is. It rejected the very record in which the ramp is visible. | +| 2026-09-07 | **The non-motion test** (Brian: "it doesn't cross zero"). `\|mean\|/peak` is bimodal with a 1.09% dead zone and returns exactly the 5 confirmed units from physics, not a threshold. Independent corroboration of the unit list. **BE11007 settled as NOT an offset** — a one-sided 9 in/s blast. | +| 2026-09-07 | **BE12599 is a connector fault, not a geophone fault.** Unipolar spikes with a 26→118 ms RC tail progressing to a flat pedestal, and the fault MOVES between Long and Tran while the sensor self-check passes on every event. An element cannot hop channels; a connector can. Inspect before reseating. | diff --git a/minimateplus/client.py b/minimateplus/client.py index d82604d..49f3e83 100644 --- a/minimateplus/client.py +++ b/minimateplus/client.py @@ -30,6 +30,7 @@ from __future__ import annotations import datetime import logging +import re import struct from typing import Optional @@ -2532,10 +2533,17 @@ def _decode_0a_partial_header(raw_data: bytes, index: int, key4: bytes) -> Optio ts2 = try_ts(raw_data[ts1_end + 1:ts1_end + 1 + ts_size]) # Extract serial and geo threshold from "BE11529\0" and "Geo: X.XXX in/s\0". + # + # Match any two-letter family prefix, not a literal "BE" — a BlastMate + # reports "BA10895", and the old `find(b"BE")` returned -1 on one. That + # skipped this whole block, so the geo threshold went missing along with + # the serial. Requiring the NUL terminator in the pattern also makes the + # match stricter than the bare two-byte search it replaces. serial: Optional[str] = None geo_ips: Optional[float] = None - serial_pos = raw_data.find(b"BE") + serial_match = re.search(rb"[A-Z]{2}\d{3,6}(?=\x00)", raw_data) + serial_pos = serial_match.start() if serial_match else -1 if serial_pos >= 0: # Read null-terminated serial starting at serial_pos. null_pos = raw_data.find(b"\x00", serial_pos) diff --git a/minimateplus/event_file_io.py b/minimateplus/event_file_io.py index 09a5bf1..d23fc6d 100644 --- a/minimateplus/event_file_io.py +++ b/minimateplus/event_file_io.py @@ -50,7 +50,7 @@ SIDECAR_KIND = "sfm.event" # bumped without a `pip install` re-run — leading to confusing stale # version stamps in sidecars. Bump this constant and CHANGELOG.md # together at release time. -TOOL_VERSION = "0.27.0" +TOOL_VERSION = "0.29.0" try: # Best-effort: prefer the installed metadata when it's NEWER than the diff --git a/pyproject.toml b/pyproject.toml index 6153c0f..cd9281b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "seismo-relay" -version = "0.27.0" +version = "0.29.0" description = "Python client and REST server for MiniMate Plus seismographs" requires-python = ">=3.10" dependencies = [ diff --git a/scratch/nonmotion_scan.py b/scratch/nonmotion_scan.py new file mode 100644 index 0000000..f6621a0 --- /dev/null +++ b/scratch/nonmotion_scan.py @@ -0,0 +1,91 @@ +#!/usr/bin/env python3 +"""Detect NON-MOTION on a geophone channel: |mean| / peak. + +A geophone is a velocity sensor with no DC response, so its output over a +record must integrate to ~zero — the ground does not relocate. Real motion +therefore sits roughly half above and half below zero. Anything electrical — +a charge-injection spike, a step, a parked pedestal — is one-sided. + + mp = |mean| / peak ~0 for motion, ~1 for a pedestal + frac_neg = share of samples < 0 ~0.3-0.5 for motion, ~0 for a fault + +Why this beats the pre-trigger floor (`offset_scan3.py`): that detector's +`spread <= 0.02` gate rejects any record whose floor is MOVING, which is +exactly what an onset is — it discarded the one BE18438 record in which the +ramp was visible. This test is indifferent to whether the fault is a spike, +a ramp or a flat pedestal; none of them cross zero. + +⚠ Not a rediscovery of the retracted v1 detector. v1 scored only the +largest-peak axis and used the mean as a BASELINE estimator, where the median +was required. Here the mean is the signal itself, per channel, and that is +what the physics licenses. +""" +from __future__ import annotations +import argparse, csv, re, statistics, sys +from concurrent.futures import ProcessPoolExecutor, as_completed +from pathlib import Path +sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) +from minimateplus.event_file_io import read_blastware_file + +GEO=("Tran","Vert","Long"); K=10.0/32000.0 +_WAVE=re.compile(r"\.[A-Za-z0-9]{2}0[Ww]$"); _STEM=re.compile(r"^([B-Z])(\d{3})") +_SER=re.compile(rb"[A-Z]{2}\d{3,6}") + +def serial_of(name, path=None): + m=_STEM.match(name) + if not m: return "?" + num=(ord(m.group(1))-ord("B"))*1000+int(m.group(2)) + if path is not None: + try: + for s in _SER.findall(Path(path).read_bytes()): + s=s.decode() + if s[2:].lstrip("0")==str(num): return s + except Exception: pass + return f"BE{num}" + +def scan(ps): + import logging; logging.disable(logging.WARNING) + p=Path(ps) + try: ev=read_blastware_file(p) + except Exception: return None + s=ev.raw_samples or {} + if not all(s.get(c) for c in GEO): return None + ts=ev.timestamp + stamp=(f"{ts.year:04d}-{ts.month:02d}-{ts.day:02d}T" + f"{ts.hour:02d}:{ts.minute:02d}:{ts.second:02d}") if ts else "" + ser=serial_of(p.name,p); out=[] + for ch in GEO: + a=[x*K for x in s[ch]] + pk=max(abs(x) for x in a) + if pk<=0: continue + out.append({"serial":ser,"timestamp":stamp,"filename":p.name,"channel":ch, + "peak":round(pk,4), + "mean":round(statistics.fmean(a),4), + "mp":round(abs(statistics.fmean(a))/pk,4), + "frac_neg":round(sum(1 for x in a if x<0)/len(a),4), + "n":len(a)}) + return out + +COLS=["serial","timestamp","filename","channel","peak","mean","mp","frac_neg","n"] + +def main(): + ap=argparse.ArgumentParser() + ap.add_argument("--dir",required=True); ap.add_argument("--out",required=True) + ap.add_argument("--jobs",type=int,default=4) + a=ap.parse_args() + seen=set(); files=[] + for q in sorted(Path(a.dir).rglob("*")): + if q.is_file() and _WAVE.search(q.name) and q.name not in seen: + seen.add(q.name); files.append(str(q)) + print(f"unique waveform binaries: {len(files)}",flush=True) + rows=[] + with ProcessPoolExecutor(max_workers=a.jobs) as ex: + for i,f in enumerate(as_completed([ex.submit(scan,p) for p in files]),1): + r=f.result() + if r: rows.extend(r) + if i%1000==0: print(f" {i}/{len(files)}",flush=True) + with open(a.out,"w",newline="") as fh: + w=csv.DictWriter(fh,fieldnames=COLS); w.writeheader(); w.writerows(rows) + print(f"\nwrote {a.out} ({len(rows)} channel-rows)") + +if __name__=="__main__": main() diff --git a/scratch/offset_hist_scan.py b/scratch/offset_hist_scan.py new file mode 100644 index 0000000..a43a097 --- /dev/null +++ b/scratch/offset_hist_scan.py @@ -0,0 +1,234 @@ +#!/usr/bin/env python3 +"""Offset detector — HISTOGRAM corpus (the other 90% of the archive). + +`offset_scan3.py` measures the pre-trigger floor in *waveform* samples. That +covers 6,577 of the archive's 70,112 unique series-3 files; the remaining +63,535 are **histograms**, which carry no samples — only a per-interval, +per-channel peak + half-period. So the pre-trigger method cannot run on them. + +The histogram analogue of "the resting floor" is the **low percentile of the +per-interval peaks**. A histogram file is typically hours of continuous +monitoring, so the great majority of its intervals are definitionally quiet; +the bottom of that distribution is what the channel reads when nothing is +happening. A healthy channel bottoms out at 0.000-0.005 in/s. A channel +parked off zero cannot report a peak below its own displacement, so its floor +is pinned up. + +⚠ The DC leakage into the histogram peak is PARTIAL. Measured within-unit +against episodes already established from the waveform scan: + + BE18438 Vert in-episode 0.0350 vs 0.0050 outside (waveform pre = +0.18..+0.37) + BE12599 Tran in-episode 0.0250 vs 0.0050 outside (waveform pre = +0.03..+0.49) + +so the device's per-interval peak is evidently measured against a running / +AC-coupled baseline that removes most, but not all, of the DC. The residual +is real and channel-specific, but the margin is ~5 quantisation counts rather +than the ~70 the waveform detector enjoys. Do not carry the waveform +detector's 0.025 in/s floor across unexamined — calibrate on the CSV. + +Because the absolute floor also moves with site noise (traffic, wind, a +generator), the statistic that matters most is the **cross-channel +differential**: a channel's floor minus the quietest of the other two geo +channels in the same file. Site noise lifts all three together and cancels; +a DC offset lifts one. + +This script does not decide anything. It emits every candidate statistic per +(file, channel) so thresholds can be calibrated against the waveform-derived +ground truth in `offset_v3.csv` rather than guessed. + +Usage: + python scratch/offset_hist_scan.py --dir /home/serversdown/dl2-archive/files \ + --out /home/serversdown/dl2-archive/offset_hist.csv --jobs 4 +""" +from __future__ import annotations + +import argparse +import csv +import datetime +import logging +import re +import statistics +import sys +from concurrent.futures import ProcessPoolExecutor, as_completed +from pathlib import Path + +sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) +from minimateplus.event_file_io import read_blastware_file # noqa: E402 + +GEO = ("Tran", "Vert", "Long") +K = 10.0 / 32000.0 # ADC count -> in/s (see CLAUDE.md: full scale 32000) +_HIST = re.compile(r"\.[A-Za-z0-9]{2}0[Hh]$") +_STEM = re.compile(r"^([B-Z])(\d{3})") +_B36 = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ" + + +_SERIAL_RE = re.compile(rb"\b([A-Z]{2}\d{3,6})\b") + + +def serial_of(name: str, path=None) -> str: + """Real serial for a BW file. + + The filename encodes only the NUMBER: `<3 digits>` where + letter = chr(ord('B') + serial // 1000). The two-letter family prefix + ("BE", "BA", ...) is **not** in the filename, so it must be read out of + the file body. Four units in the DL2 archive are BA, not BE — assuming + "BE" mislabels BA9229, BA10060, BA10895 and BA15957. + """ + m = _STEM.match(name) + if not m: + return "?" + num = (ord(m.group(1)) - ord("B")) * 1000 + int(m.group(2)) + if path is not None: + try: + for s in _SERIAL_RE.findall(Path(path).read_bytes()): + s = s.decode() + if s[2:].lstrip("0") == str(num): + return s + except Exception: + pass + return f"BE{num}" # last-resort fallback; prefix unverified + + +def stem_time(name: str): + """Decode the filename's base-36 timestamp. Epoch 1985-01-01, 1296 s/tick. + + Preferred over the file's own footer timestamp only because it costs + nothing; the caller falls back to the decoded event when this fails. + """ + try: + base, ext = name.rsplit(".", 1) + n = 0 + for c in base[4:8].upper(): + n = n * 36 + _B36.index(c) + ab = _B36.index(ext[0].upper()) * 36 + _B36.index(ext[1].upper()) + return datetime.datetime(1985, 1, 1) + datetime.timedelta(seconds=n * 1296 + ab) + except Exception: + return None + + +def _pct(sorted_vals, q): + """Nearest-rank percentile on an already-sorted list.""" + if not sorted_vals: + return None + i = min(len(sorted_vals) - 1, max(0, int(len(sorted_vals) * q / 100.0))) + return sorted_vals[i] + + +def scan(path_str: str): + logging.disable(logging.WARNING) # per-worker: the codec warns on undecodables + p = Path(path_str) + try: + ev = read_blastware_file(p) + except Exception: + return None + s = ev.raw_samples or {} + if not any(s.get(c) for c in GEO): + return None + + ts = stem_time(p.name) or ev.timestamp + stamp = "" + if ts is not None: + stamp = (f"{ts.year:04d}-{ts.month:02d}-{ts.day:02d}T" + f"{ts.hour:02d}:{ts.minute:02d}:{ts.second:02d}") + + # Per-channel floor candidates, in in/s. + stats = {} + for ch in GEO: + v = sorted(s.get(ch) or []) + if not v: + continue + stats[ch] = { + "n": len(v), + "min": v[0] * K, + "p1": _pct(v, 1) * K, + "p5": _pct(v, 5) * K, + "p10": _pct(v, 10) * K, + "p25": _pct(v, 25) * K, + "med": statistics.median(v) * K, + "peak": v[-1] * K, + "zeros": sum(1 for x in v if x == 0) / len(v), + } + if len(stats) < 2: # need at least one sibling channel for the differential + return None + + # Mic floor as a site-noise proxy (raw counts; the dB conversion is not + # needed — only its relative movement matters here). + mic = sorted(s.get("MicL") or []) + mic_p5 = _pct(mic, 5) if mic else "" + + rows = [] + for ch, st in stats.items(): + others = [stats[o]["p5"] for o in stats if o != ch] + rows.append({ + "serial": serial_of(p.name, p), + "timestamp": stamp, + "filename": p.name, + "channel": ch, + "n_intervals": st["n"], + "min": round(st["min"], 4), + "p1": round(st["p1"], 4), + "p5": round(st["p5"], 4), + "p10": round(st["p10"], 4), + "p25": round(st["p25"], 4), + "median": round(st["med"], 4), + "peak": round(st["peak"], 4), + "frac_zero": round(st["zeros"], 4), + # the site-noise-cancelling statistic: this channel's floor above + # the quietest sibling geo channel in the same file + "diff_p5": round(st["p5"] - min(others), 4), + "mic_p5": mic_p5, + }) + return rows + + +COLS = ["serial", "timestamp", "filename", "channel", "n_intervals", + "min", "p1", "p5", "p10", "p25", "median", "peak", "frac_zero", + "diff_p5", "mic_p5"] + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("--dir", required=True) + ap.add_argument("--out", required=True) + ap.add_argument("--jobs", type=int, default=4) + ap.add_argument("--limit", type=int, default=0, help="stop after N files (smoke test)") + a = ap.parse_args() + + # Dedupe by basename — the DL2 export keeps a byte-identical `Sent/` + # mirror of its root, which doubled two figures before it was caught. + seen, files = set(), [] + for q in sorted(Path(a.dir).rglob("*")): + if q.is_file() and _HIST.search(q.name) and q.name not in seen: + seen.add(q.name) + files.append(str(q)) + if a.limit: + files = files[:a.limit] + print(f"unique histogram binaries: {len(files)}", flush=True) + + rows, undecodable = [], 0 + with ProcessPoolExecutor(max_workers=a.jobs) as ex: + futs = [ex.submit(scan, f) for f in files] + for i, fut in enumerate(as_completed(futs), 1): + r = fut.result() + if r: + rows.extend(r) + else: + undecodable += 1 + if i % 5000 == 0: + print(f" {i}/{len(files)}", flush=True) + + with open(a.out, "w", newline="") as fh: + w = csv.DictWriter(fh, fieldnames=COLS) + w.writeheader() + w.writerows(rows) + + files_ok = len({r["filename"] for r in rows}) + units = len({r["serial"] for r in rows}) + ivals = sum(r["n_intervals"] for r in rows) // 3 + print(f"\ndecoded {files_ok}/{len(files)} files " + f"({undecodable} undecodable), {units} units, ~{ivals/1e6:.1f}M intervals") + print(f"wrote {a.out} ({len(rows)} channel-rows)") + + +if __name__ == "__main__": + main() diff --git a/scratch/offset_scan3.py b/scratch/offset_scan3.py index 2c2b953..3a3cde6 100644 --- a/scratch/offset_scan3.py +++ b/scratch/offset_scan3.py @@ -28,9 +28,31 @@ from minimateplus.event_file_io import read_blastware_file GEO=("Tran","Vert","Long"); K=10.0/32000.0 _WAVE=re.compile(r"\.[A-Za-z0-9]{2}0[Ww]$"); _STEM=re.compile(r"^([B-Z])(\d{3})") -def serial_of(n): - m=_STEM.match(n) - return f"BE{(ord(m.group(1))-ord('B'))*1000+int(m.group(2))}" if m else "?" +_SERIAL_RE = re.compile(rb"\b([A-Z]{2}\d{3,6})\b") + + +def serial_of(name: str, path=None) -> str: + """Real serial for a BW file. + + The filename encodes only the NUMBER: `<3 digits>` where + letter = chr(ord('B') + serial // 1000). The two-letter family prefix + ("BE", "BA", ...) is **not** in the filename, so it must be read out of + the file body. Four units in the DL2 archive are BA, not BE — assuming + "BE" mislabels BA9229, BA10060, BA10895 and BA15957. + """ + m = _STEM.match(name) + if not m: + return "?" + num = (ord(m.group(1)) - ord("B")) * 1000 + int(m.group(2)) + if path is not None: + try: + for s in _SERIAL_RE.findall(Path(path).read_bytes()): + s = s.decode() + if s[2:].lstrip("0") == str(num): + return s + except Exception: + pass + return f"BE{num}" # last-resort fallback; prefix unverified def scan(ps): p=Path(ps) @@ -48,7 +70,7 @@ def scan(ps): mid, end = a[t:2*t], a[2*t:] if not pre or not mid or not end: continue v=[statistics.median(x)*K for x in (pre,mid,end)] - out.append({"serial":serial_of(p.name),"timestamp":stamp, + out.append({"serial":serial_of(p.name, p),"timestamp":stamp, "filename":p.name,"channel":ch, "pretrig_n": pre_n or 0, "pre":round(v[0],4),"mid":round(v[1],4),"end":round(v[2],4), diff --git a/scripts/backfill_event_shape.py b/scripts/backfill_event_shape.py index 975c716..4b7c308 100644 --- a/scripts/backfill_event_shape.py +++ b/scripts/backfill_event_shape.py @@ -1,12 +1,12 @@ #!/usr/bin/env python3 -"""Backfill events.shape_* from each event's .h5 waveform samples. Idempotent.""" +"""Backfill events.shape_* and shape_offset_* from each event's .h5 samples. Idempotent.""" from __future__ import annotations import argparse, logging, sys from pathlib import Path sys.path.insert(0, str(Path(__file__).resolve().parent.parent)) from sfm.database import SeismoDb from sfm.waveform_store import WaveformStore -from sfm.shape_metrics import shape_from_h5 +from sfm.shape_metrics import shape_from_h5, offset_from_h5 log = logging.getLogger("backfill_event_shape") @@ -21,6 +21,7 @@ def backfill_shape(db: SeismoDb, store: WaveformStore, *, dry_run: bool = False) if not h5_path.exists(): counts["skipped_no_h5"] += 1; continue shape = shape_from_h5(h5_path) + offset = offset_from_h5(h5_path) if shape is None: # The .h5 can no longer yield a shape (fewer than 2 samples, or a # flat trace). Clear any previously stored value rather than @@ -28,22 +29,32 @@ def backfill_shape(db: SeismoDb, store: WaveformStore, *, dry_run: bool = False) # from and silently feeds the false-trigger detector. Seen after # a decoder fix shrinks an event: 493 rows in the prod snapshot # were carrying metrics from a superseded decode (2026-08-25). - if row.get("shape_crest_factor") is not None: + if (row.get("shape_crest_factor") is not None + or row.get("shape_offset") is not None): if not dry_run: with db._connect() as conn: conn.execute( "UPDATE events SET shape_crest_factor=NULL, " "shape_near_peak_count=NULL, shape_sample_count=NULL, " - "shape_axis=NULL WHERE id=?", (row["id"],)) + "shape_axis=NULL, shape_offset=NULL, shape_offset_axis=NULL, " + "shape_offset_pre=NULL, shape_offset_spread=NULL WHERE id=?", + (row["id"],)) counts["cleared_stale"] += 1 counts["skipped_no_samples"] += 1; continue if not dry_run: with db._connect() as conn: conn.execute( "UPDATE events SET shape_crest_factor=?, shape_near_peak_count=?, " - "shape_sample_count=?, shape_axis=? WHERE id=?", + "shape_sample_count=?, shape_axis=?, shape_offset=?, " + "shape_offset_axis=?, shape_offset_pre=?, shape_offset_spread=? " + "WHERE id=?", (shape["crest_factor"], shape["near_peak_count"], - shape["sample_count"], shape["axis"], row["id"])) + shape["sample_count"], shape["axis"], + (1 if offset["offset"] else 0) if offset else None, + offset["axis"] if offset else None, + offset["pre"] if offset else None, + offset["spread"] if offset else None, + row["id"])) counts["updated"] += 1 log.info("backfill_shape: %s", counts) return counts diff --git a/sfm/database.py b/sfm/database.py index 04c5d24..b2ebf08 100644 --- a/sfm/database.py +++ b/sfm/database.py @@ -82,6 +82,7 @@ CREATE TABLE IF NOT EXISTS events ( record_type TEXT, -- "single_shot" | "continuous" false_trigger INTEGER NOT NULL DEFAULT 0, -- 0=no, 1=yes (manual flag) reviewed_real INTEGER NOT NULL DEFAULT 0, -- 0=no, 1=operator-confirmed real (mutually exclusive with false_trigger) + false_trigger_reason TEXT, -- optional FT cause ("offset", ...); NULL = none. Only meaningful when false_trigger=1. blastware_filename TEXT, -- event file within waveform store; extension is per-event (AB0T encodes timestamp) blastware_filesize INTEGER, -- bytes; NULL if no event file saved a5_pickle_filename TEXT, -- ".a5.pkl" sidecar @@ -99,6 +100,10 @@ CREATE TABLE IF NOT EXISTS events ( shape_near_peak_count INTEGER, -- samples >= 0.5 * peak (FT: few; real: many) shape_sample_count INTEGER, -- total samples (to normalize near_peak_count) shape_axis TEXT, -- geophone channel measured ("Tran"/"Vert"/"Long") + shape_offset INTEGER, -- 1 = DC-offset false trigger (pre-trigger baseline off zero + flat). Meaningful for waveforms only. + shape_offset_axis TEXT, -- geo channel the offset was measured on + shape_offset_pre REAL, -- pre-trigger baseline median (in/s) + shape_offset_spread REAL, -- max(pre,mid,end) - min(...) in in/s; small = constant/DC created_at TEXT NOT NULL DEFAULT (strftime('%Y-%m-%dT%H:%M:%SZ', 'now')), UNIQUE(serial, timestamp) ); @@ -225,7 +230,12 @@ class SeismoDb: ("shape_near_peak_count", "INTEGER"), ("shape_sample_count", "INTEGER"), ("shape_axis", "TEXT"), + ("shape_offset", "INTEGER"), + ("shape_offset_axis", "TEXT"), + ("shape_offset_pre", "REAL"), + ("shape_offset_spread", "REAL"), ("reviewed_real", "INTEGER NOT NULL DEFAULT 0"), + ("false_trigger_reason", "TEXT"), ): if col not in existing_cols: log.info("_migrate: events ADD COLUMN %s %s", col, ddl) @@ -430,9 +440,11 @@ class SeismoDb: tran_zc_above_range, vert_zc_above_range, long_zc_above_range, mic_zc_above_range, shape_crest_factor, shape_near_peak_count, - shape_sample_count, shape_axis) + shape_sample_count, shape_axis, + shape_offset, shape_offset_axis, + shape_offset_pre, shape_offset_spread) VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, - ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?) + ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?) """, ( self._new_id(), serial, key, session_id, ts, @@ -464,6 +476,10 @@ class SeismoDb: rec.get("shape_near_peak_count"), rec.get("shape_sample_count"), rec.get("shape_axis"), + rec.get("shape_offset"), + rec.get("shape_offset_axis"), + rec.get("shape_offset_pre"), + rec.get("shape_offset_spread"), ), ) inserted += 1 @@ -517,7 +533,11 @@ class SeismoDb: shape_crest_factor = COALESCE(?, shape_crest_factor), shape_near_peak_count = COALESCE(?, shape_near_peak_count), shape_sample_count = COALESCE(?, shape_sample_count), - shape_axis = COALESCE(?, shape_axis) + shape_axis = COALESCE(?, shape_axis), + shape_offset = COALESCE(?, shape_offset), + shape_offset_axis = COALESCE(?, shape_offset_axis), + shape_offset_pre = COALESCE(?, shape_offset_pre), + shape_offset_spread = COALESCE(?, shape_offset_spread) WHERE serial = ? AND timestamp = ? """, ( @@ -549,6 +569,10 @@ class SeismoDb: rec.get("shape_near_peak_count") if rec else None, rec.get("shape_sample_count") if rec else None, rec.get("shape_axis") if rec else None, + rec.get("shape_offset") if rec else None, + rec.get("shape_offset_axis") if rec else None, + rec.get("shape_offset_pre") if rec else None, + rec.get("shape_offset_spread") if rec else None, serial, ts, ), @@ -691,9 +715,9 @@ class SeismoDb: def propagate_review_to_twins(self, event_id: str, *, window_seconds: int | None = None) -> list[str]: """ - Copy this event's `false_trigger`/`reviewed_real` columns onto each - of its histogram/waveform twins (see `find_twins`), so flagging one - twin flags both. Returns the list of twin ids updated. + Copy this event's `false_trigger`/`reviewed_real`/`false_trigger_reason` + columns onto each of its histogram/waveform twins (see `find_twins`), so + flagging one twin flags both. Returns the list of twin ids updated. ``window_seconds`` is accepted for backward compatibility but ignored; twin matching is now interval-based (see `find_twins`). @@ -703,12 +727,16 @@ class SeismoDb: return [] ft = 1 if row.get("false_trigger") else 0 real = 1 if row.get("reviewed_real") else 0 + # The reason is a subtype of the FT flag — carry it only when the source + # is actually a false trigger, so a confirmed-real twin never keeps one. + reason = row.get("false_trigger_reason") if ft else None twins = self.find_twins(event_id) moved = [] with self._connect() as conn: for tw in twins: - conn.execute("UPDATE events SET false_trigger=?, reviewed_real=? WHERE id=?", - (ft, real, tw["id"])) + conn.execute( + "UPDATE events SET false_trigger=?, reviewed_real=?, false_trigger_reason=? WHERE id=?", + (ft, real, reason, tw["id"])) moved.append(tw["id"]) return moved @@ -729,7 +757,7 @@ class SeismoDb: ) else: cur = conn.execute( - "UPDATE events SET false_trigger=0 WHERE id=?", + "UPDATE events SET false_trigger=0, false_trigger_reason=NULL WHERE id=?", (event_id,), ) return cur.rowcount > 0 @@ -823,7 +851,8 @@ class SeismoDb: return False has_ft = "false_trigger" in review has_real = "reviewed_real" in review - if not has_ft and not has_real: + has_reason = "false_trigger_reason" in review + if not has_ft and not has_real and not has_reason: # Nothing derived to update; just confirm the row exists. with self._connect() as conn: row = conn.execute( @@ -836,11 +865,19 @@ class SeismoDb: sets["false_trigger"] = 1 if review.get("false_trigger") else 0 if has_real: sets["reviewed_real"] = 1 if review.get("reviewed_real") else 0 + if has_reason: + reason = review.get("false_trigger_reason") or None + sets["false_trigger_reason"] = reason + if reason: # a reason is a subtype of FT → implies FT + sets["false_trigger"] = 1 # mutual exclusivity: a true in one forces the other column to 0 if sets.get("false_trigger") == 1: sets["reviewed_real"] = 0 if sets.get("reviewed_real") == 1: sets["false_trigger"] = 0 + # the reason is only meaningful while flagged FT — clear it if FT ends up 0 + if sets.get("false_trigger") == 0: + sets["false_trigger_reason"] = None assign = ", ".join(f"{k}=?" for k in sets) params = list(sets.values()) + [event_id] with self._connect() as conn: diff --git a/sfm/shape_metrics.py b/sfm/shape_metrics.py index 05c9dd5..cd02b80 100644 --- a/sfm/shape_metrics.py +++ b/sfm/shape_metrics.py @@ -47,6 +47,60 @@ def shape_from_samples(chans: dict) -> dict | None: return s +# ── Offset (DC-baseline) detection ──────────────────────────────────────────── +# A DC offset is a false trigger where the geophone baseline sits at a constant +# non-zero floor (sensor bumped / settled / drifted) instead of oscillating +# around zero. Brian's method (validated in scratch/offset_scan3.py): the +# pre-trigger window is definitionally quiet, so a true offset shows |pre| off +# zero AND stays flat across the record (pre ≈ mid ≈ end). A transient moves one +# third relative to the others and is rejected by the spread test. +# Thresholds are in in/s (the .h5 samples are already range-scaled); validated at +# Normal range (10 in/s) — the only range in the fleet. +OFFSET_FLOOR = 0.025 # |pre| at/above this reads as an off-zero baseline (5 A/D counts) +OFFSET_MAX_SPREAD = 0.02 # max(pre,mid,end) - min(...) at/below this reads as flat/constant + + +def _channel_offset(x, pretrig_n): + """Return (pre, spread, is_offset) for one channel, or None if unusable.""" + x = np.asarray(x, dtype=float) + n = x.size + if n < 3: + return None + t = n // 3 + pre = x[:pretrig_n] if (pretrig_n and 0 < pretrig_n < n) else x[:t] + mid, end = x[t:2 * t], x[2 * t:] + if pre.size == 0 or mid.size == 0 or end.size == 0: + return None + vals = [float(np.median(seg)) for seg in (pre, mid, end)] + spread = max(vals) - min(vals) + is_offset = abs(vals[0]) >= OFFSET_FLOOR and spread <= OFFSET_MAX_SPREAD + return vals[0], spread, is_offset + + +def offset_from_samples(chans: dict, pretrig_n) -> dict | None: + """Detect a DC-offset false trigger across the geophone channels. + + An event is offset if ANY geo channel's pre-trigger baseline is off zero and + flat across the record. Reports the tripping axis (or, if none trips, the + most-offset-like axis) with its ``pre``/``spread`` for transparency + tuning. + Returns None when no geo channel is usable. + """ + results = [] + for ax in _GEO_CHANNELS: + x = chans.get(ax) + if x is None: + continue + r = _channel_offset(x, pretrig_n) + if r is not None: + results.append((ax, r[0], r[1], r[2])) + if not results: + return None + offenders = [r for r in results if r[3]] + ax, pre, spread, _ = max(offenders or results, key=lambda r: abs(r[1])) + return {"offset": bool(offenders), "axis": ax, + "pre": round(pre, 6), "spread": round(spread, 6)} + + def shape_from_h5(path) -> dict | None: import h5py try: @@ -56,3 +110,18 @@ def shape_from_h5(path) -> dict | None: except Exception: return None return shape_from_samples(chans) + + +def offset_from_h5(path) -> dict | None: + """offset_from_samples fed from an event's .h5 (float32 in/s geo samples + + the pretrig_samples attribute).""" + import h5py + try: + with h5py.File(path, "r") as f: + chans = {ax: f[f"samples/{ax}"][:] for ax in _GEO_CHANNELS + if f"samples/{ax}" in f} + pretrig_n = f.attrs.get("pretrig_samples") + except Exception: + return None + pretrig_n = int(pretrig_n) if pretrig_n is not None else 0 + return offset_from_samples(chans, pretrig_n) diff --git a/sfm/waveform_store.py b/sfm/waveform_store.py index 6b25e14..8373e8d 100644 --- a/sfm/waveform_store.py +++ b/sfm/waveform_store.py @@ -32,6 +32,7 @@ from __future__ import annotations import datetime import logging import pickle +import re import shutil from pathlib import Path from typing import Optional, Union @@ -41,7 +42,7 @@ from minimateplus.blastware_file import blastware_filename, write_blastware_file from minimateplus.framing import S3Frame from minimateplus.models import Event from sfm import event_hdf5 -from sfm.shape_metrics import shape_from_h5 +from sfm.shape_metrics import shape_from_h5, offset_from_h5 log = logging.getLogger("sfm.waveform_store") @@ -270,6 +271,13 @@ class WaveformStore: "shape_sample_count": _shape["sample_count"], "shape_axis": _shape["axis"], } if _shape else {} + _offset = offset_from_h5(hdf5_path) if hdf5_filename else None + _offset_rec = { + "shape_offset": 1 if _offset["offset"] else 0, + "shape_offset_axis": _offset["axis"], + "shape_offset_pre": _offset["pre"], + "shape_offset_spread": _offset["spread"], + } if _offset else {} return { "filename": filename, "filesize": filesize, @@ -278,6 +286,7 @@ class WaveformStore: "hdf5_filename": hdf5_filename, "sidecar_filename": sidecar_path.name, **_shape_rec, + **_offset_rec, } def save_imported_bw( @@ -371,8 +380,16 @@ class WaveformStore: # Resolve serial. blastware_filename derives a 4-char prefix from # the numeric serial (e.g. BE11529 → M529); we go the other way - # via the source filename if a hint wasn't given. - serial = serial_hint or _serial_from_bw_filename(source_path.name) or "UNKNOWN" + # if a hint wasn't given. The filename carries only the NUMBER, + # so read the family prefix out of the body first — a BlastMate + # ("BA") filed as "BE" is a unit that does not exist. The + # filename-only decoder stays as the last resort. + serial = ( + serial_hint + or _serial_from_bw_bytes(bw_bytes, source_path.name) + or _serial_from_bw_filename(source_path.name) + or "UNKNOWN" + ) # Use the source filename verbatim — it already encodes timestamp # + record type per BW's AB0T scheme, and we want to preserve it @@ -461,6 +478,13 @@ class WaveformStore: "shape_sample_count": _shape["sample_count"], "shape_axis": _shape["axis"], } if _shape else {} + _offset = offset_from_h5(hdf5_path) if hdf5_filename else None + _offset_rec = { + "shape_offset": 1 if _offset["offset"] else 0, + "shape_offset_axis": _offset["axis"], + "shape_offset_pre": _offset["pre"], + "shape_offset_spread": _offset["spread"], + } if _offset else {} return ev, { "filename": filename, "filesize": filesize, @@ -470,6 +494,7 @@ class WaveformStore: "sidecar_filename": sidecar_path.name, "serial": serial, **_shape_rec, + **_offset_rec, } def save_imported_idf( @@ -751,6 +776,13 @@ class WaveformStore: "shape_sample_count": _shape["sample_count"], "shape_axis": _shape["axis"], } if _shape else {} + _offset = offset_from_h5(hdf5_path) if hdf5_filename else None + _offset_rec = { + "shape_offset": 1 if _offset["offset"] else 0, + "shape_offset_axis": _offset["axis"], + "shape_offset_pre": _offset["pre"], + "shape_offset_spread": _offset["spread"], + } if _offset else {} return ev, { "filename": filename, "filesize": filesize, @@ -760,6 +792,7 @@ class WaveformStore: "sidecar_filename": sidecar_path.name, "serial": serial, **_shape_rec, + **_offset_rec, } def load_a5(self, serial: str, filename: str) -> Optional[list[S3Frame]]: @@ -816,20 +849,24 @@ class WaveformStore: # ── helpers ───────────────────────────────────────────────────────────────────── -def _serial_from_bw_filename(name: str) -> Optional[str]: +def _serial_number_from_bw_filename(name: str) -> Optional[int]: """ - Reverse of `blastware_filename`'s serial-prefix encoding. + Reverse of `blastware_filename`'s serial-prefix encoding — the NUMBER only. BW filename format (V10.72): `

.` where P = chr(ord('B') + floor(serial // 1000)) and serial3 = f"{serial % 1000:03d}". Examples (from CLAUDE.md verification archive): - P036... → BE14036 H907... → BE6907 - M529... → BE11529 T003... → BE18003 + P036... → 14036 H907... → 6907 + M529... → 11529 T003... → 18003 + L895... → 10895 - Returns the inferred BE-prefix serial (e.g. "BE11529") or None when - the filename doesn't match the expected pattern. + ⚠ The filename encodes **only the number**. The two-letter family + prefix is NOT in it — "BE" is a MiniMate Plus, "BA" a BlastMate — so + the prefix has to come from the file body (`_serial_from_bw_bytes`) + or from an explicit hint. Returns None when the filename doesn't + match the expected pattern. """ if not name: return None @@ -842,5 +879,43 @@ def _serial_from_bw_filename(name: str) -> Optional[str]: if prefix_letter < "B": return None thousands = ord(prefix_letter) - ord("B") - serial_num = thousands * 1000 + int(base[1:4]) - return f"BE{serial_num}" + return thousands * 1000 + int(base[1:4]) + + +_BW_SERIAL_RE = re.compile(rb"[A-Z]{2}\d{3,6}") + + +def _serial_from_bw_bytes(data: bytes, name: str) -> Optional[str]: + """ + Read the real serial — prefix included — out of a BW file body. + + The body carries the serial as a plain ASCII string ("BE9558", + "BA10895"). We accept a candidate only when its numeric part matches + the number the filename encodes, which keeps a stray byte sequence in + the sample stream from being mistaken for a serial. + + Returns None when the filename number can't be derived or no + candidate in the body agrees with it — the caller then falls back. + """ + num = _serial_number_from_bw_filename(name) + if num is None or not data: + return None + for match in _BW_SERIAL_RE.findall(data): + candidate = match.decode("ascii", errors="replace") + if candidate[2:].lstrip("0") == str(num): + return candidate + return None + + +def _serial_from_bw_filename(name: str) -> Optional[str]: + """ + Best-effort serial from the filename alone. + + ⚠ The family prefix is a **guess** — the filename does not carry it. + "BE" is right for every MiniMate Plus but wrong for a BlastMate, whose + serials start "BA". Prefer `_serial_from_bw_bytes` whenever the file + body is at hand; this exists for callers that only have a name + (log lines, dry-run output). + """ + num = _serial_number_from_bw_filename(name) + return None if num is None else f"BE{num}" diff --git a/tests/test_ft_reason.py b/tests/test_ft_reason.py new file mode 100644 index 0000000..f0699d1 --- /dev/null +++ b/tests/test_ft_reason.py @@ -0,0 +1,103 @@ +import sqlite3 + +from sfm.database import SeismoDb +from minimateplus.models import Event, Timestamp + + +def _ev(db, key="0111aaaa", serial="BE1"): + ev = Event(index=0) + ev._waveform_key = bytes.fromhex(key) + ev.timestamp = Timestamp(raw=b"", flag=0x10, year=2026, unknown_byte=0, + month=6, day=25, hour=8, minute=0, second=0) + ev.record_type = "Waveform" + db.insert_events([ev], serial=serial) + return [r for r in db.query_events(serial=serial) if r["waveform_key"] == key][0]["id"] + + +def test_flag_offset_reason_implies_ft(tmp_path): + db = SeismoDb(tmp_path / "s.db") + eid = _ev(db) + db.update_event_review(eid, {"false_trigger_reason": "offset"}) + row = db.get_event(eid) + assert row["false_trigger"] == 1 # a reason is a subtype of FT + assert row["false_trigger_reason"] == "offset" + assert row["reviewed_real"] == 0 + + +def test_plain_ft_leaves_reason_null(tmp_path): + # Reason is OPTIONAL — flagging FT without one records no reason. + db = SeismoDb(tmp_path / "s.db") + eid = _ev(db) + db.update_event_review(eid, {"false_trigger": True}) + row = db.get_event(eid) + assert row["false_trigger"] == 1 + assert row["false_trigger_reason"] is None + + +def test_confirm_real_clears_reason(tmp_path): + db = SeismoDb(tmp_path / "s.db") + eid = _ev(db) + db.update_event_review(eid, {"false_trigger_reason": "offset"}) + db.update_event_review(eid, {"reviewed_real": True}) + row = db.get_event(eid) + assert row["reviewed_real"] == 1 + assert row["false_trigger"] == 0 + assert row["false_trigger_reason"] is None + + +def test_clear_ft_clears_reason(tmp_path): + db = SeismoDb(tmp_path / "s.db") + eid = _ev(db) + db.update_event_review(eid, {"false_trigger_reason": "offset"}) + db.update_event_review(eid, {"false_trigger": False}) + row = db.get_event(eid) + assert row["false_trigger"] == 0 + assert row["false_trigger_reason"] is None + + +def test_set_false_trigger_false_clears_reason(tmp_path): + db = SeismoDb(tmp_path / "s.db") + eid = _ev(db) + db.update_event_review(eid, {"false_trigger_reason": "offset"}) + assert db.set_false_trigger(eid, False) is True + row = db.get_event(eid) + assert row["false_trigger"] == 0 + assert row["false_trigger_reason"] is None + + +def test_reason_can_be_cleared_without_clearing_ft(tmp_path): + # Setting reason to None removes the reason but leaves the FT flag intact. + db = SeismoDb(tmp_path / "s.db") + eid = _ev(db) + db.update_event_review(eid, {"false_trigger_reason": "offset"}) + db.update_event_review(eid, {"false_trigger_reason": None}) + row = db.get_event(eid) + assert row["false_trigger"] == 1 + assert row["false_trigger_reason"] is None + + +def _ts(h, m, d=25): + return Timestamp(raw=b"", flag=0x10, year=2026, unknown_byte=0, + month=2, day=d, hour=h, minute=m, second=0) + + +def test_offset_reason_propagates_to_twin(tmp_path): + # Flag a waveform as offset → its histogram twin also becomes FT with reason=offset. + db = SeismoDb(tmp_path / "s.db") + + def ins(key, ts, rt): + ev = Event(index=0); ev._waveform_key = bytes.fromhex(key); ev.timestamp = ts + db.insert_events([ev], serial="BE1") + rid = [r for r in db.query_events(serial="BE1") if r["waveform_key"] == key][0]["id"] + with sqlite3.connect(db.db_path) as c: + c.execute("UPDATE events SET peak_vector_sum=0.4763, record_type=? WHERE id=?", (rt, rid)) + return rid + + hist = ins("01110001", _ts(19, 31), "Histogram") # interval start + wave = ins("01110002", _ts(20, 46), "Waveform") # trigger inside the interval + + db.update_event_review(wave, {"false_trigger_reason": "offset"}) + db.propagate_review_to_twins(wave) + row = db.get_event(hist) + assert row["false_trigger"] == 1 + assert row["false_trigger_reason"] == "offset" diff --git a/tests/test_offset_metrics.py b/tests/test_offset_metrics.py new file mode 100644 index 0000000..e349a82 --- /dev/null +++ b/tests/test_offset_metrics.py @@ -0,0 +1,94 @@ +import numpy as np +import h5py +from sfm.shape_metrics import offset_from_samples, offset_from_h5 + + +def test_flags_constant_dc_floor(): + # A geophone channel sitting at a constant +0.05 in/s across the whole record + # is a DC offset: baseline off zero AND flat across pre/mid/end thirds. + n = 300 + chans = {"Tran": np.full(n, 0.05), "Vert": np.zeros(n), "Long": np.zeros(n)} + r = offset_from_samples(chans, pretrig_n=50) + assert r["offset"] is True + assert r["axis"] == "Tran" + assert abs(r["pre"] - 0.05) < 1e-6 + assert r["spread"] < 0.02 + + +def test_transient_rejected_by_spread(): + # Off-zero pre-trigger but the baseline SETTLES back over the record — a + # transient, not a constant offset. The spread test must reject it. + x = np.concatenate([np.full(100, 0.05), np.full(100, 0.025), np.zeros(100)]) + chans = {"Tran": x, "Vert": np.zeros(300), "Long": np.zeros(300)} + r = offset_from_samples(chans, pretrig_n=100) + assert r["offset"] is False + + +def test_clean_oscillation_not_offset(): + t = np.arange(300) + x = 0.4 * np.sin(2 * np.pi * t / 20) # oscillates around zero — baseline IS zero + chans = {"Tran": x, "Vert": np.zeros(300), "Long": np.zeros(300)} + r = offset_from_samples(chans, pretrig_n=50) + assert r["offset"] is False + + +def test_below_floor_not_offset_but_reports_pre(): + # A flat baseline below the floor is not an offset; still report the axis/pre + # for tuning transparency. + n = 300 + chans = {"Tran": np.full(n, 0.01), "Vert": np.zeros(n), "Long": np.zeros(n)} + r = offset_from_samples(chans, pretrig_n=50) + assert r["offset"] is False + assert r["axis"] == "Tran" + assert abs(r["pre"] - 0.01) < 1e-6 + + +def test_none_when_no_geo_channels(): + assert offset_from_samples({"MicL": np.full(300, 0.05)}, pretrig_n=50) is None + + +def test_pretrig_fallback_when_invalid(): + # pretrig_n of 0 (missing/unusable) falls back to the first third. + n = 300 + chans = {"Tran": np.full(n, 0.05), "Vert": np.zeros(n), "Long": np.zeros(n)} + r = offset_from_samples(chans, pretrig_n=0) + assert r["offset"] is True + + +def test_flags_offset_on_any_axis(): + # Offset on Vert alone still flags the event, and Vert is reported. + n = 300 + chans = {"Tran": np.zeros(n), "Vert": np.full(n, -0.06), "Long": np.zeros(n)} + r = offset_from_samples(chans, pretrig_n=50) + assert r["offset"] is True + assert r["axis"] == "Vert" + + +def _write_h5(path, chans, pretrig_n): + with h5py.File(path, "w") as f: + g = f.create_group("samples") + for k, v in chans.items(): + g.create_dataset(k, data=np.asarray(v, dtype="float32")) + if pretrig_n is not None: + f.attrs["pretrig_samples"] = pretrig_n + + +def test_offset_from_h5_reads_pretrig_attr(tmp_path): + p = tmp_path / "ev.h5" + n = 300 + _write_h5(p, {"Tran": np.full(n, 0.05), "Vert": np.zeros(n), "Long": np.zeros(n)}, + pretrig_n=50) + r = offset_from_h5(str(p)) + assert r["offset"] is True and r["axis"] == "Tran" + + +def test_offset_from_h5_missing_pretrig_attr_falls_back(tmp_path): + p = tmp_path / "noattr.h5" + n = 300 + _write_h5(p, {"Tran": np.full(n, 0.05), "Vert": np.zeros(n), "Long": np.zeros(n)}, + pretrig_n=None) + assert offset_from_h5(str(p))["offset"] is True # falls back to first-third + + +def test_offset_from_h5_missing_file_is_none(tmp_path): + assert offset_from_h5(str(tmp_path / "nope.h5")) is None diff --git a/tests/test_offset_storage.py b/tests/test_offset_storage.py new file mode 100644 index 0000000..73cf069 --- /dev/null +++ b/tests/test_offset_storage.py @@ -0,0 +1,64 @@ +from __future__ import annotations +from pathlib import Path + +import numpy as np, h5py + +from sfm.database import SeismoDb +from sfm.waveform_store import WaveformStore +from scripts.backfill_event_shape import backfill_shape +from minimateplus.models import Event, Timestamp, PeakValues + +_FIX = Path(__file__).parent / "fixtures/histogram-extension-re/events-5-21-26/K558LL8B.7I0W" + + +def _event(waveform_key="0111abcd"): + ev = Event(index=0) + ev._waveform_key = bytes.fromhex(waveform_key) + ev.timestamp = Timestamp(raw=b"", flag=0x10, year=2026, unknown_byte=0, + month=6, day=25, hour=8, minute=50, second=0) + ev.record_type = "Waveform" + ev.peak_values = PeakValues(tran=0.075, vert=0.220, long=0.045, + peak_vector_sum=0.231, micl=0.01) + return ev + + +def test_insert_stores_offset_from_record(tmp_path: Path): + db = SeismoDb(tmp_path / "s.db") + ev = _event() + rec = {ev._waveform_key.hex(): { + "filename": "F.CE0W", "filesize": 10, + "shape_offset": 1, "shape_offset_axis": "Tran", + "shape_offset_pre": 0.05, "shape_offset_spread": 0.001}} + db.insert_events([ev], serial="BE1", waveform_records=rec) + row = db.query_events(serial="BE1")[0] + assert row["shape_offset"] == 1 + assert row["shape_offset_axis"] == "Tran" + assert abs(row["shape_offset_pre"] - 0.05) < 1e-6 + assert abs(row["shape_offset_spread"] - 0.001) < 1e-6 + + +def test_save_imported_bw_attaches_offset(tmp_path: Path): + store = WaveformStore(tmp_path / "waveforms") + ev, rec = store.save_imported_bw(_FIX.read_bytes(), source_path=_FIX, serial_hint="BE9558") + assert rec["shape_offset"] in (0, 1) + assert rec["shape_offset_axis"] in ("Tran", "Vert", "Long") + assert "shape_offset_pre" in rec and "shape_offset_spread" in rec + + +def test_backfill_updates_offset(tmp_path: Path): + db = SeismoDb(tmp_path / "s.db") + store = WaveformStore(tmp_path / "waveforms") + ev = Event(index=0); ev._waveform_key = bytes.fromhex("0111abcd") + db.insert_events([ev], serial="BE1", + waveform_records={ev._waveform_key.hex(): {"filename": "F.CE0W", "filesize": 10}}) + p = store.hdf5_path_for("BE1", "F.CE0W") + with h5py.File(p, "w") as f: + g = f.create_group("samples") + g.create_dataset("Tran", data=np.full(300, 0.05, "float32")) + g.create_dataset("Vert", data=np.zeros(300, "float32")) + g.create_dataset("Long", data=np.zeros(300, "float32")) + f.attrs["pretrig_samples"] = 50 + backfill_shape(db, store) + row = db.query_events(serial="BE1")[0] + assert row["shape_offset"] == 1 + assert row["shape_offset_axis"] == "Tran" diff --git a/tests/test_serial_prefix.py b/tests/test_serial_prefix.py new file mode 100644 index 0000000..9d2d111 --- /dev/null +++ b/tests/test_serial_prefix.py @@ -0,0 +1,101 @@ +"""The BW filename encodes the serial NUMBER, never the family prefix. + +"BE" is a MiniMate Plus; "BA" is a BlastMate. Both are Series III and their +files are byte-compatible — the whole archive's 1,493 BlastMate binaries +decode through the same codec at 100% — so the only thing that distinguishes +them downstream is the serial string, and that lives in the file body. + +Synthesising the prefix as "BE" files a BlastMate under a unit that does not +exist. Four units in the DL2 archive are affected: BA9229, BA10060, BA10895 +and BA15957. +""" +from __future__ import annotations + +import pytest + +from minimateplus.client import _decode_0a_partial_header +from sfm.waveform_store import ( + _serial_from_bw_bytes, + _serial_from_bw_filename, + _serial_number_from_bw_filename, +) + + +# ── the filename gives a number, and only a number ────────────────────────── + +@pytest.mark.parametrize("name,num", [ + ("P036L318.C80H", 14036), # BE14036 + ("H907KWRK.WB0H", 6907), # BE6907 + ("M529LKIQ.G10", 11529), # BE11529 + ("T003LQ9K.OE0H", 18003), # BE18003 + ("L895K63F.GE0W", 10895), # BA10895 — a BlastMate + ("K229HGQI.XO0W", 9229), # BA9229 — a BlastMate +]) +def test_number_from_filename(name, num): + assert _serial_number_from_bw_filename(name) == num + + +@pytest.mark.parametrize("name", ["", "not_a_bw_file.bin", "AB12", "1234ABCD.XX0W"]) +def test_number_from_filename_rejects_junk(name): + assert _serial_number_from_bw_filename(name) is None + + +def test_filename_only_decoder_is_a_guess(): + """It still answers "BE" — that is why it must not be the first choice.""" + assert _serial_from_bw_filename("L895K63F.GE0W") == "BE10895" + assert _serial_from_bw_filename("M529LKIQ.G10") == "BE11529" + assert _serial_from_bw_filename("nonsense") is None + + +# ── the body carries the truth ────────────────────────────────────────────── + +def _body(serial: bytes) -> bytes: + return b"\x00" * 32 + b"STRT" + b"\xff\xfe" + serial + b"\x00Geo: 0.254 in/s\x00" + + +def test_body_wins_for_a_blastmate(): + assert _serial_from_bw_bytes(_body(b"BA10895"), "L895K63F.GE0W") == "BA10895" + + +def test_body_wins_for_a_minimate(): + assert _serial_from_bw_bytes(_body(b"BE11529"), "M529LKIQ.G10") == "BE11529" + + +def test_body_candidate_must_match_the_filename_number(): + """A serial-shaped byte run that disagrees with the filename is ignored.""" + assert _serial_from_bw_bytes(_body(b"XX99999"), "L895K63F.GE0W") is None + + +def test_body_tolerates_a_leading_zero(): + assert _serial_from_bw_bytes(_body(b"BA09229"), "K229HGQI.XO0W") == "BA09229" + + +@pytest.mark.parametrize("data,name", [ + (b"", "L895K63F.GE0W"), # no bytes + (_body(b"BA10895"), "junk.bin"), # no derivable number +]) +def test_body_returns_none_when_it_cannot_decide(data, name): + assert _serial_from_bw_bytes(data, name) is None + + +# ── the live monitor-log path ─────────────────────────────────────────────── + +def _partial_record(serial: bytes) -> bytes: + """0x2C partial record: type, prefix, two 9-byte timestamps, then ASCII.""" + ts = bytes([11, 0x10, 4, 0x07, 0xE9, 0, 16, 2, 0]) # 2025-04-11 16:02:00 + return (bytes([0x2C]) + b"\x00" * 10 + ts + ts + + b"\x00\x00\x00\x00" + serial + b"\x00Geo: 0.254 in/s\x00") + + +@pytest.mark.parametrize("serial", [b"BE11529", b"BA10895", b"UM11719"]) +def test_monitor_log_reads_any_family_prefix(serial): + entry = _decode_0a_partial_header(_partial_record(serial), 0, b"\x01\x11\x00\x00") + assert entry is not None + assert entry.serial == serial.decode() + + +def test_monitor_log_geo_threshold_survives_a_blastmate(): + """The old find(b"BE") skipped the whole block, losing geo too.""" + entry = _decode_0a_partial_header(_partial_record(b"BA10895"), 0, b"\x01\x11\x00\x00") + assert entry is not None + assert entry.geo_threshold_ips == pytest.approx(0.254)