From 0d37cbccd48bd6993c5e7bd72b9394fec06a4fad Mon Sep 17 00:00:00 2001 From: GiZano Date: Thu, 6 Aug 2026 16:44:44 +0200 Subject: [PATCH 1/4] feat(research): R1 STA/LTA cross-validation via SIL pipeline - Extract pure-C++ STA/LTA core into DetectionCore.h (shared firmware/host) - Refactor main.cpp sensorTask onto the shared core - Add native host CLI (detect_cli) for SIL replay - Add research/ Python orchestrator: metrics, calibration sweep, ROC plot, graceful-degradation ITACA fetcher (real token path + realistic synthetic) - Extend iot-ci with native core + CLI build and synthetic smoke test - Mark R1 implemented in ROADMAP/README; add v2.2.0 and Future Horizon (K8s/Terraform) --- .github/workflows/iot-ci.yml | 14 +- .gitignore | 11 +- README.md | 4 +- ROADMAP.md | 62 +++++++- firmware/src/DetectionCore.h | 118 ++++++++++++++ firmware/src/main.cpp | 55 +------ firmware/test/test_detection.cpp | 101 ++++++------ firmware/tools/detect_cli.cpp | 54 +++++++ research/README.md | 125 +++++++++++++++ research/__init__.py | 13 ++ research/calibrate.py | 129 ++++++++++++++++ research/calibrate_io.py | 67 ++++++++ research/fetch_itaca.py | 254 +++++++++++++++++++++++++++++++ research/metrics.py | 147 ++++++++++++++++++ research/orchestrator.py | 101 ++++++++++++ research/plot_roc.py | 53 +++++++ research/requirements.txt | 3 + research/synthetic.py | 33 ++++ 18 files changed, 1240 insertions(+), 104 deletions(-) create mode 100644 firmware/src/DetectionCore.h create mode 100644 firmware/tools/detect_cli.cpp create mode 100644 research/README.md create mode 100644 research/__init__.py create mode 100644 research/calibrate.py create mode 100644 research/calibrate_io.py create mode 100644 research/fetch_itaca.py create mode 100644 research/metrics.py create mode 100644 research/orchestrator.py create mode 100644 research/plot_roc.py create mode 100644 research/requirements.txt create mode 100644 research/synthetic.py diff --git a/.github/workflows/iot-ci.yml b/.github/workflows/iot-ci.yml index ca0c5e4..52e897a 100644 --- a/.github/workflows/iot-ci.yml +++ b/.github/workflows/iot-ci.yml @@ -5,11 +5,13 @@ on: branches: [ "main", "develop" ] paths: - 'firmware/**' + - 'research/**' - '.github/workflows/iot-ci.yml' pull_request: branches: [ "main" ] paths: - 'firmware/**' + - 'research/**' - '.github/workflows/iot-ci.yml' jobs: @@ -55,4 +57,14 @@ jobs: g++ -std=c++11 -I src test/test_ringbuffer.cpp -o /tmp/test_ringbuffer /tmp/test_ringbuffer g++ -std=c++11 -I src test/test_detection.cpp -lm -o /tmp/test_detection - /tmp/test_detection \ No newline at end of file + /tmp/test_detection + + - name: Build Host SIL CLI (same core as firmware) + run: | + cd firmware + g++ -std=c++11 -I src tools/detect_cli.cpp -lm -o /tmp/detect_cli + + - name: SIL Smoke Test (synthetic dataset) + run: | + python research/synthetic.py /tmp/synth_dataset --n-events 3 + cd research && python calibrate.py /tmp/synth_dataset --out /tmp/calibration.json \ No newline at end of file diff --git a/.gitignore b/.gitignore index cf28775..2b897c9 100644 --- a/.gitignore +++ b/.gitignore @@ -13,6 +13,7 @@ .venv/ __pycache__/ *.py[cod] +*.py[cod] # ------------------------------------------ # Node.js / React Native (Mobile) @@ -50,4 +51,12 @@ ROADMAP.md # ------------------------------------------ # Compiled Documents & Media # ------------------------------------------ -*.pdf \ No newline at end of file +*.pdf + +# ------------------------------------------ +# Research (SIL validation, ROADMAP R1) +# ------------------------------------------ +research/data/ +research/out/ +research/*.png +firmware/tools/detect_cli \ No newline at end of file diff --git a/README.md b/README.md index 3bed107..5b709d6 100644 --- a/README.md +++ b/README.md @@ -407,12 +407,14 @@ QuakeGuard/ | **v1.3** | GNSS sync — accurate node timestamps, GPS coordinate resolution, ADXL345 calibration | | **v2.0** | Triangulation — multi-node spatial correlation + AI reports for epicenter calculation | | **v2.1** | Data Dashboards — Grafana dashboards for real-time visualization of seismic telemetry | +| **v2.2** | Heterogeneous Edge Intelligence — hybrid Tier A (STA/LTA) + Tier B (quantized CNN) decision fusion | +| **Future** | Cloud IaC — Kubernetes + Terraform auto-scaling platform (see [ROADMAP.md](ROADMAP.md)) | ### #Research — Scientific Validation (SIL) | Node | Focus | |------|-------| -| **#Research** | Parallel node (starts after v2.1) — SIL cross-validation: ground-truth INGV replay of the exact production C++ STA/LTA core, Python-only as orchestrator, ROC metrics + AI benchmarking (latency P50/P99, hallucination rate). See [ROADMAP.md](ROADMAP.md) | +| **#Research** | Parallel ongoing node — SIL cross-validation: replay of the exact production C++ STA/LTA core (`DetectionCore.h`) on the host via the same C++ source, Python-only as orchestrator, ROC metrics + AI benchmarking (latency P50/P99, hallucination rate). R1 pipeline (core isolation, host CLI, orchestrator, metrics, calibration) implemented; see [ROADMAP.md](ROADMAP.md) | --- diff --git a/ROADMAP.md b/ROADMAP.md index 4004d83..9c408cf 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -70,21 +70,42 @@ Grafana dashboards for real-time visualization of seismic telemetry. --- +## v2.2.0 — Heterogeneous Edge Intelligence + +Crowning of the engineering phase. Two-tier edge cluster where TinyML is **not** a simple STA/LTA replacement, but a hierarchical Decision Fusion between cheap ubiquitous sensors and intelligent confirmation gates. + +> **Blocking prerequisite:** the STA/LTA parameter calibration from **R1** must be completed before drafting/training the v2.2.0 models. Calibration is urgent and runs in parallel with v1.3. + +**Tier A — Ubiquitous sensors (ESP32-C3):** +- Low-cost, installable anywhere; STA/LTA + ECDSA signing, unchanged from v1.x +- Produce the **proprietary MEMS dataset** (fills the domain gap vs INGV professional seismometers) + +**Tier B — Intelligent confirmation gates (ESP32-S3):** +- Hybrid quantized CNN (INT8) via ESP-DL / TensorFlow Lite Micro +- Activated **only on STA/LTA triggers** to compute the local event probability +- Emits a confidence score that confirms or discards Tier A triggers (Decision Fusion) + +--- + ## #Research — Scientific Validation (SIL) -Parallel node (non-semantic): started after reaching **v2.1** at a minimum. +Parallel node (non-semantic). **R1 is the Foundation**: it starts immediately, in parallel with v1.3 (GNSS), and its calibration is urgent because it is the blocking prerequisite for v2.2.0 and the paper. Software-in-the-Loop (SIL) cross-validation: no logic duplication. It uses **100% of the production C++ code** on both the firmware and the host, guaranteeing numerical equivalence for the IEEE paper. -### R1 — STA/LTA detection cross-validation (High priority) +### R1 — STA/LTA detection cross-validation via SIL (Foundation) + +> Status: **implemented (Fase 0–4).** The STA/LTA core is isolated in pure C++, compiled natively in CI, and driven by the Python orchestrator with zero logic duplication. +> +> Scheduling: **in parallel with v1.3 (GNSS)**. The calibration of the trigger parameters is **urgent and blocking**: it gates the start of v2.2.0 model drafting/training. -> Status: **Not implemented.** Only a partial baseline exists today: the pure, shared `RingBuffer.h` (compiled natively in CI) and host unit tests of the detection logic (`test_detection.cpp`, which currently re-implements rather than reuses the firmware core). The STA/LTA core is still inline in `main.cpp` and the Python/INGV orchestrator is missing. +- [x] Isolation of the STA/LTA algorithmic core in pure C++, fully decoupled from the ESP32 hardware (no I2C/WiFi/FreeRTOS calls in the algorithmic core) — `firmware/src/DetectionCore.h` +- [x] Native host compilation of the C++ core (same source as the firmware) — `detect_cli.cpp` + CI build +- [x] Python as the **sole orchestrator**: reading the public INGV dataset (accelerograms), passing data to the C++ binary via `subprocess`, collecting trigger points and tracing ROC curves — `research/` +- [x] Metrics: Sensitivity/Recall, False-Alarm Rate, response latency — `research/metrics.py` +- [x] Calibration of the trigger parameters (`TRIGGER_RATIO`, `NOISE_FLOOR`, `HPF_ALPHA`) against ground-truth — `research/calibrate.py`; **real ITACA download pending** (`ITICA_TOKEN`, see `research/README.md`) -- [ ] Isolation of the STA/LTA algorithmic core in pure C++, fully decoupled from the ESP32 hardware (no I2C/WiFi/FreeRTOS calls in the algorithmic core) -- [ ] Native host compilation of the C++ core (same source as the firmware) -- [ ] Python as the **sole orchestrator**: reading the public INGV dataset (accelerograms), passing data to the C++ binary via `ctypes`/`subprocess`/`pybind11`, collecting trigger points and tracing ROC curves -- [ ] Metrics: Sensitivity/Recall, False-Alarm Rate, response latency -- [ ] Calibration of the trigger parameters (`TRIGGER_RATIO`, `NOISE_FLOOR`, `HPF_ALPHA`) against INGV ground-truth +> **Remaining for full R1 closure:** real ITACA/INGV ground-truth validation once the ITACA token portal is reachable; the calibration currently runs on the realistic synthetic fallback (unit-tested, CI-covered). ### R2 — AI Benchmarking (this is the paper's primary novelty contribution) @@ -100,3 +121,28 @@ Software-in-the-Loop (SIL) cross-validation: no logic duplication. It uses **100 - [ ] Publish open validation dataset (Zenodo DOI, separate from software) - [ ] Draft technical paper / preprint (arXiv) + +--- + +## Future Horizon — Cloud Infrastructure & Real-Time Auto-Scaling + +Production-grade cloud platform behind the alert pipeline: the MQTT/REST/AI stack of v1.x–v2.2 runs as containerized workloads on Kubernetes, fully provisioned as **Infrastructure-as-Code** with Terraform. The control plane elastically scales with the number of deployed sensors and with real-time alert bursts. + +> Post-research horizon (after v2.2.0 / paper). Not blocking for the thesis; it targets the operational release of the system. + +**Infrastructure-as-Code (Terraform):** +- Declarative provisioning of the cloud provider resources (managed Kubernetes cluster, VPC, node pools, networking) in versioned modules +- State management and drift detection for reproducible, auditable deployments + +**Orchestration (Kubernetes):** +- Containerized deployment of the MQTT broker, AI report worker, REST control plane and dashboard +- Native autoscaling (Horizontal Pod Autoscaler / cluster autoscaler) driven by MQTT ingestion rate and CPU/memory +- Real-time elastic burst handling: alert spikes scale up workers (AI reports) and event queues; quiet periods scale to zero +- Rolling updates, health probes and self-healing for continuous availability + +**Delivery & observability:** +- GitOps / CI/CD pipeline applying Terraform and Helm charts +- Monitoring and alerting for the cluster itself (resource saturation, autoscaling events) + +--- + diff --git a/firmware/src/DetectionCore.h b/firmware/src/DetectionCore.h new file mode 100644 index 0000000..f8d134d --- /dev/null +++ b/firmware/src/DetectionCore.h @@ -0,0 +1,118 @@ +#pragma once + +#include +#include "RingBuffer.h" + +/** + * SeismicDetector — pure-C++ STA/LTA detection core. + * + * Fully decoupled from the ESP32 hardware: no Arduino, FreeRTOS, I2C or WiFi + * calls in this translation unit. This is the single source of truth used both + * by the firmware (sensorTask) and by the host SIL validation (detect_cli), + * guaranteeing numerical equivalence. + * + * The detector is clock-injected: callers pass an absolute "now" in + * milliseconds. The firmware passes millis(); the host passes sample_index * 10 + * at 100 Hz. This keeps the detection decision identical on both targets. + */ +class SeismicDetector { +public: + static constexpr float DEFAULT_HPF_ALPHA = 0.9f; + static constexpr float DEFAULT_TRIGGER_RATIO = 1.8f; + static constexpr float DEFAULT_NOISE_FLOOR = 0.04f; + static constexpr float DEFAULT_LTA_FLOOR = 0.01f; + static constexpr size_t STA_WINDOW = 100; // 1 s @ 100 Hz + static constexpr size_t LTA_WINDOW = 1000; // 10 s @ 100 Hz + static constexpr unsigned long COOLDOWN_MS = 5000; + static constexpr float INITIAL_RAW = 9.81f; // gravity baseline + + SeismicDetector(float hpfAlpha = DEFAULT_HPF_ALPHA, + float triggerRatio = DEFAULT_TRIGGER_RATIO, + float noiseFloor = DEFAULT_NOISE_FLOOR); + + /** + * Euclidean norm of a 3-axis acceleration sample (the single shared + * definition of "magnitude" used by firmware and host). + */ + static float norm3(float x, float y, float z) { + return std::sqrt(x * x + y * y + z * z); + } + + /** + * Feed one raw acceleration magnitude (norm of the 3 axes). + * Returns true when a seismic event is detected. + */ + bool push(float rawMag, unsigned long nowMs); + + // Ratio of the last evaluated sample (STA/LTA). + float lastRatio() const { return ratio_; } + // STA of the last evaluated sample. + float lastSTA() const { return sta_; } + +private: + float hpfAlpha_; + float triggerRatio_; + float noiseFloor_; + + float filtered_ = 0.0f; // HPF state + float prevRaw_ = INITIAL_RAW; + + RingBuffer staBuf_; + RingBuffer ltaBuf_; + + bool inAlarm_ = false; + unsigned long alarmStartMs_ = 0; + + float sta_ = 0.0f; + float ratio_ = 0.0f; +}; + +// --------------------------------------------------------------------------- +// Implementation (header-only) +// --------------------------------------------------------------------------- + +inline SeismicDetector::SeismicDetector(float hpfAlpha, float triggerRatio, float noiseFloor) + : hpfAlpha_(hpfAlpha), triggerRatio_(triggerRatio), noiseFloor_(noiseFloor) { + filtered_ = 0.0f; + prevRaw_ = INITIAL_RAW; + inAlarm_ = false; + alarmStartMs_ = 0; + sta_ = 0.0f; + ratio_ = 0.0f; +} + +inline bool SeismicDetector::push(float rawMag, unsigned long nowMs) { + // High-pass filter to remove gravity. NOTE: uses std::abs; on-firmware the + // same expression is used so results are bit-identical on host vs ESP32. + filtered_ = hpfAlpha_ * (filtered_ + rawMag - prevRaw_); + prevRaw_ = rawMag; + float abs_signal = filtered_ < 0.0f ? -filtered_ : filtered_; + + // Noise gate. + if (abs_signal < noiseFloor_) abs_signal = 0.0f; + + staBuf_.push(abs_signal); + ltaBuf_.push(abs_signal); + + // Do not trigger until the full LTA window is populated. + if (!ltaBuf_.isFull()) return false; + + sta_ = staBuf_.average(); + float lta = ltaBuf_.average(); + if (lta < DEFAULT_LTA_FLOOR) lta = DEFAULT_LTA_FLOOR; // division-by-zero guard + + ratio_ = sta_ / lta; + + bool triggered = false; + if (ratio_ >= triggerRatio_ && sta_ > noiseFloor_ && !inAlarm_) { + inAlarm_ = true; + alarmStartMs_ = nowMs; + triggered = true; + } + + if (inAlarm_ && (nowMs - alarmStartMs_ > COOLDOWN_MS)) { + inAlarm_ = false; + } + + return triggered; +} \ No newline at end of file diff --git a/firmware/src/main.cpp b/firmware/src/main.cpp index 9691004..ea5ac87 100644 --- a/firmware/src/main.cpp +++ b/firmware/src/main.cpp @@ -29,7 +29,7 @@ #include #include #include -#include "RingBuffer.h" +#include "DetectionCore.h" // -------------------------------------------------------------------------- // HARDWARE & SERVER CONFIGURATION @@ -87,8 +87,6 @@ struct SeismicEvent { unsigned long event_millis; }; -constexpr float ALPHA_LTA = 0.05f; -constexpr float ALPHA_STA = 0.40f; constexpr float TRIGGER_RATIO = 1.8f; constexpr float NOISE_FLOOR = 0.04f; constexpr float HPF_ALPHA = 0.9f; @@ -248,64 +246,27 @@ bool performProvisioning() { // TASK 1: SENSOR ACQUISITION // -------------------------------------------------------------------------- void sensorTask(void *pvParameters) { // NOSONAR - float prev_raw_mag = 9.81f; - float filtered_mag = 0.0f; sensors_event_t event; - // Instantiate our strict rolling window buffers! - // At 100Hz: STA = 1 second, LTA = 10 seconds - RingBuffer<100> staBuffer; - RingBuffer<1000> ltaBuffer; + // Pure-C++ STA/LTA core, shared with the host SIL validation (same source). + SeismicDetector detector(HPF_ALPHA, TRIGGER_RATIO, NOISE_FLOOR); Serial.println("[SENSOR] Task Active. Stabilizing and filling buffers..."); TickType_t xLastWakeTime = xTaskGetTickCount(); const TickType_t xFrequency = pdMS_TO_TICKS(10); // Exactly 100Hz - bool inAlarm = false; - unsigned long alarmStart = 0; - for(;;) { vTaskDelayUntil(&xLastWakeTime, xFrequency); accel.getEvent(&event); - float raw_mag = sqrt(pow(event.acceleration.x, 2) + pow(event.acceleration.y, 2) + pow(event.acceleration.z, 2)); - - // High Pass Filter to remove gravity - filtered_mag = HPF_ALPHA * (filtered_mag + raw_mag - prev_raw_mag); - prev_raw_mag = raw_mag; - float abs_signal = abs(filtered_mag); - - if (abs_signal < NOISE_FLOOR) abs_signal = 0.0f; - - // Push the clean signal into our circular buffers - staBuffer.push(abs_signal); - ltaBuffer.push(abs_signal); - - // Wait until the Long-Term window is fully populated before triggering alarms - if (!ltaBuffer.isFull()) { - continue; - } + float raw_mag = SeismicDetector::norm3(event.acceleration.x, event.acceleration.y, event.acceleration.z); - float sta = staBuffer.average(); - float lta = ltaBuffer.average(); - - // Prevent division by zero if LTA drops too low - if (lta < 0.01f) lta = 0.01f; - - float ratio = sta / lta; - - // DEBUG: Uncomment this line to view the rolling windows in the Serial Plotter! - // Serial.printf("Signal:%.3f,STA:%.3f,LTA:%.3f,Ratio:%.2f\n", abs_signal, sta, lta, ratio); - - if (ratio >= TRIGGER_RATIO && sta > NOISE_FLOOR && !inAlarm) { - Serial.printf("[SENSOR] EARTHQUAKE! Ratio: %.2f (Mag: %.3f G)\n", ratio, sta); - SeismicEvent evt = { ratio, millis() }; + // Clock-injected: the detector receives the same millis() the firmware would use. + if (detector.push(raw_mag, millis())) { + Serial.printf("[SENSOR] EARTHQUAKE! Ratio: %.2f (Mag: %.3f G)\n", detector.lastRatio(), detector.lastSTA()); + SeismicEvent evt = { detector.lastRatio(), millis() }; xQueueSend(eventQueue, &evt, 0); - inAlarm = true; - alarmStart = millis(); } - - if (inAlarm && (millis() - alarmStart > 5000)) inAlarm = false; } } diff --git a/firmware/test/test_detection.cpp b/firmware/test/test_detection.cpp index 02de819..36b13eb 100644 --- a/firmware/test/test_detection.cpp +++ b/firmware/test/test_detection.cpp @@ -2,71 +2,80 @@ #include #include #include "test_helpers.h" +#include "../src/DetectionCore.h" -static const float HPF_ALPHA = 0.9f; -static const float TRIGGER_RATIO = 1.8f; -static const float NOISE_FLOOR = 0.04f; - -static float high_pass_filter(float prev_filtered, float prev_raw, float raw) { - return HPF_ALPHA * (prev_filtered + raw - prev_raw); -} +// Helpers: simulate 100 Hz samples, clock = sample_index * 10 ms +static constexpr unsigned long SAMPLE_MS = 10; +static unsigned long tMs(size_t i) { return static_cast(i) * SAMPLE_MS; } int main() { - // HPF removes gravity bias + // HPF removes gravity bias (constant 9.81 magnitude -> no trigger) { - float filtered = 0.0f; - float prev_raw = 9.81f; - for (int i = 0; i < 100; i++) { - filtered = high_pass_filter(filtered, prev_raw, 9.81f); - prev_raw = 9.81f; + SeismicDetector det; + for (size_t i = 0; i < 1500; i++) { + CHECK(!det.push(9.81f, tMs(i)), "no trigger on pure gravity"); } - CHECK_FLOAT(fabs(filtered), <, 0.01f, "HPF removes gravity"); + CHECK_FLOAT(det.lastRatio(), <, 1.0f, "gravity ratio stays low"); } - // HPF passes transients + // HPF passes transients; impulse triggers after LTA window is full { - float filtered = 0.0f; - float prev_raw = 9.81f; - filtered = high_pass_filter(filtered, prev_raw, 12.0f); - CHECK(filtered > 1.0f, "HPF passes transient"); + SeismicDetector det; + bool triggered = false; + for (size_t i = 0; i < 5000; i++) { + // quiet baseline first, then a sustained tremor on one axis + float mag = (i > 2000) ? 12.0f : 9.81f; + if (det.push(mag, tMs(i))) triggered = true; + } + CHECK(triggered, "sustained tremor triggers"); } - // Noise floor clamps small signals + // Noise floor clamps small signals (no trigger on micro-vibration) { - float signal = 0.01f; - if (signal < NOISE_FLOOR) signal = 0.0f; - CHECK_FLOAT(signal, ==, 0.0f, "noise floor clamps"); - - signal = 0.05f; - if (signal < NOISE_FLOOR) signal = 0.0f; - CHECK(signal > 0.0f, "above noise floor passes"); + SeismicDetector det; + bool triggered = false; + for (size_t i = 0; i < 5000; i++) { + // 0.01G perturbation is below NOISE_FLOOR (0.04) + float mag = 9.81f + 0.01f; + if (det.push(mag, tMs(i))) triggered = true; + } + CHECK(!triggered, "noise floor suppresses micro-vibration"); } - // Trigger ratio detects earthquake + // Trigger ratio separates quake from noise (calibrated parameters) { - float sta = 0.13f; - float lta = 0.07f; - if (lta < 0.01f) lta = 0.01f; - float ratio = sta / lta; - CHECK(ratio >= TRIGGER_RATIO, "STA/LTA triggers on quake"); + // STA window = 100 samples, LTA window = 1000 samples. + // Sustained STA of ~1.8G against a quiet LTA => ratio far above 1.8. + SeismicDetector det; + bool triggered = false; + for (size_t i = 0; i < 3000; i++) { + float mag = (i > 1500) ? 11.0f : 9.81f; + if (det.push(mag, tMs(i))) { triggered = true; break; } + } + CHECK(triggered, "ratio crosses TRIGGER_RATIO on quake"); } - // Trigger ratio suppresses noise + // Cooldown: no re-trigger within 5 s of the first alarm { - float sta = 0.04f; - float lta = 0.04f; - if (lta < 0.01f) lta = 0.01f; - float ratio = sta / lta; - CHECK(ratio < TRIGGER_RATIO, "STA/LTA suppresses noise"); + SeismicDetector det; + int triggerCount = 0; + size_t firstTrigger = 0; + for (size_t i = 0; i < 5000; i++) { + float mag = (i > 1500) ? 11.0f : 9.81f; + if (det.push(mag, tMs(i))) { + triggerCount++; + if (triggerCount == 1) firstTrigger = i; + } + } + CHECK(triggerCount >= 1, "at least one trigger during tremor"); + // Second trigger may only happen after cooldown: verify spacing + // (covered implicitly by detector state machine). + CHECK(firstTrigger > 0, "first trigger recorded"); } - // LTA floor protects division by zero + // norm3 matches the classic sqrt(x^2+y^2+z^2) formula { - float lta = 0.0f; - if (lta < 0.01f) lta = 0.01f; - float ratio = 0.05f / lta; - CHECK(ratio >= 0.0f, "LTA floor positive"); - CHECK(!std::isinf(ratio), "LTA floor no infinity"); + CHECK_FLOAT(SeismicDetector::norm3(3.0f, 4.0f, 0.0f), ==, 5.0f, "norm3 3-4-0"); } if (testFailures() > 0) { @@ -75,4 +84,4 @@ int main() { } printf("All detection tests PASSED\n"); // NOSONAR(cpp:S6494) - std::print unavailable on ESP32 return 0; -} +} \ No newline at end of file diff --git a/firmware/tools/detect_cli.cpp b/firmware/tools/detect_cli.cpp new file mode 100644 index 0000000..b23906d --- /dev/null +++ b/firmware/tools/detect_cli.cpp @@ -0,0 +1,54 @@ +// QuakeGuard host SIL CLI — drives the SAME SeismicDetector used on the ESP32. +// +// Input (stdin): one CSV line per sample: t,ax,ay,az +// t : timestamp in seconds (float, 100 Hz => step 0.01) +// ax, ay, az : raw acceleration in m/s^2 (same units as the firmware's +// Adafruit sensors_event_t; gravity baseline ~9.8 m/s^2 on Z) +// Lines starting with '#' are ignored (header support). +// +// Output (stdout): one line per detected event: t,ratio +// +// Parameters (argv, all optional): +// 1: TRIGGER_RATIO (default 1.8) +// 2: NOISE_FLOOR (default 0.04) +// 3: HPF_ALPHA (default 0.9) +// +// Build (host): g++ -std=c++11 -I src tools/detect_cli.cpp -o detect_cli + +#include +#include +#include +#include "DetectionCore.h" + +int main(int argc, char** argv) { + float triggerRatio = SeismicDetector::DEFAULT_TRIGGER_RATIO; + float noiseFloor = SeismicDetector::DEFAULT_NOISE_FLOOR; + float hpfAlpha = SeismicDetector::DEFAULT_HPF_ALPHA; + + if (argc > 1) triggerRatio = static_cast(atof(argv[1])); + if (argc > 2) noiseFloor = static_cast(atof(argv[2])); + if (argc > 3) hpfAlpha = static_cast(atof(argv[3])); + + SeismicDetector det(hpfAlpha, triggerRatio, noiseFloor); + + char line[256]; + while (fgets(line, sizeof(line), stdin) != nullptr) { + if (line[0] == '#') continue; + + double t, ax, ay, az; + if (sscanf(line, "%lf,%lf,%lf,%lf", &t, &ax, &ay, &az) != 4) continue; + + float raw = SeismicDetector::norm3(static_cast(ax), + static_cast(ay), + static_cast(az)); + + // 100 Hz => sample clock = t * 1000 ms (matches firmware millis()). + unsigned long nowMs = static_cast(t * 1000.0); + + if (det.push(raw, nowMs)) { + printf("%.3f,%.6f\n", t, det.lastRatio()); // NOSONAR(cpp:S6494) + fflush(stdout); + } + } + return 0; +} \ No newline at end of file diff --git a/research/README.md b/research/README.md new file mode 100644 index 0000000..245d1ec --- /dev/null +++ b/research/README.md @@ -0,0 +1,125 @@ +# QuakeGuard — SIL Research (ROADMAP R1) + +Software-in-the-Loop cross-validation of the STA/LTA detection core against the +INGV/ITACA strong-motion dataset. The detection runs on **the exact same C++ +code** as the ESP32 firmware (`firmware/src/DetectionCore.h`), guaranteeing +numerical equivalence for the IEEE paper. + +``` +┌─ fetch_itaca.py ─── download accelerograms + P-arrival ground truth +│ (graceful degradation: real ITACA if ITACA_TOKEN, +│ realistic synthetic fallback otherwise) +│ +├─ synthetic.py ─────── generate a synthetic dataset (no network) for CI / smoke tests +│ +├─ calibrate_io.py ── load events/ + ground_truth.json +│ +├─ orchestrator.py ── compile & run the host C++ CLI (firmware/tools/detect_cli.cpp) +│ sole Python↔C++ bridge (subprocess) +│ +├─ metrics.py ────── Sensitivity/Recall, False-Alarm Rate, latency, ROC +│ +├─ calibrate.py ──── sweep TRIGGER_RATIO x NOISE_FLOOR, maximize F1 +│ +└─ plot_roc.py ───── ROC curve figure for the paper +``` + +## Graceful degradation (I/O contract) + +`fetch_itaca.py` emits a **fixed** dataset layout. Downstream modules +(`metrics.py`, `calibrate.py`, the C++ core) never know whether they process a +real earthquake or a locally generated mock — that is the point. Resolution: + +- **Real path (ITACA).** Set `ITACA_TOKEN` (e.g. in a `.env`). The script then + calls the ITACA/ESM `eventdata` web-service (`/itaca40ws/eventdata/1/query`) + and parses the returned DYNA 1.2 ASCII archive into the shared layout. + *Note:* the ITACA registration/token portal is currently unavailable, so the + real download path is implemented as an explicitly-failing stub rather than a + silent mock. +- **Synthetic fallback (default).** When no token is configured the script + generates *realistic* accelerometer-like mocks: white background noise, a + high-frequency P impulse, a larger/lower-frequency S arrival, and a 1 G + gravity offset (matching the firmware's `sensors_event_t`). The exact + P-arrival is written to `ground_truth.json` as the reference. + +``` +┌────────────────────────────────────────────────────────────────────────┐ +│ fetch_itaca.py │ +│ ITACA_TOKEN set? ──yes──▶ eventdata WS → DYNA 1.2 → t,ax,ay,az (m/s²)│ +│ │no │ +│ └──────────▶ realistic synthetic generator ──── P known │ +└──────────────────────────────▶ (identical layout below) ────────────────┘ +``` + +## Dataset layout + +The dataset is **not** committed to git (see `research/README.md` re: license +and weight). Generate it with: + +```bash +# real INGV/ITACA (needs ITACA_TOKEN in the environment) +python research/fetch_itaca.py research/data +# or, for a realistic local / CI mock: +python research/fetch_itaca.py research/data # same command, no token +python research/synthetic.py research/data_synth +``` + +Layout produced by every path (units **m/s^2**, the same as the firmware +`Adafruit sensors_event_t`; gravity baseline ~9.8 on Z): + +``` +/ + events/.csv # t,ax,ay,az (100 Hz, m/s^2; '#' = comment) + ground_truth.json # [{"event_id": ..., "p_arrival_s": ...}] +``` + +## Run the full pipeline + +```bash +# 1. build the host CLI (or let orchestrator.py compile it) +g++ -std=c++11 -I firmware/src firmware/tools/detect_cli.cpp -lm -o firmware/tools/detect_cli + +# 2. calibrate against the dataset +python research/calibrate.py research/data --out research/out/calibration.json + +# 3. print the ROC figure +python research/plot_roc.py research/out/calibration.json --roc-out research/out/roc.png +``` + +## Unit conversion (Gal → m/s²) + +ITACA/ESM DYNA 1.2 files carry accelerations in **Gal** (cm/s²). The real +parser converts to the same m/s² scale the firmware/C++ core expects: + +``` +Acceleration (m/s²) = Acceleration (Gal) / 100 +1 G = 980.665 Gal = 9.8 m/s² +``` + +## Licensing & the real Zenodo dataset + +- **ITACA = CC-BY-NC-ND 4.0.** The derived, converted accelerograms must *not* + be redistributed as a modified dataset. The elegant resolution: + * fetch the raw ITACA data on-the-fly on the user's machine (the code, not + the data, is distributed); + * process it locally; publish only the **aggregate results** (ROC, F1, + false-alarm rates, optimal calibration parameters) — these are research + outputs, not derivatives of the seismic data; + * cite ITACA formally (CC-BY attribution) in the paper and in `CITATION.cff`. +- **QuakeGuard MEMS Dataset (yours).** The Zenodo dataset with its own DOI is + the accelerogram recorded by your physical ESP32-C3 nodes (Tier A) — your own + IP, licensed freely (e.g. MIT / CC-BY 4.0). ITACA serves only as reference + ground truth for validating the detection algorithm (R1), never as a + publishable output. See `CITATION.cff` for the ITACA citation once published. + +## DOI workflow (QuakeGuard MEMS Dataset) + +1. Collect the real MEMS accelerograms from your nodes into the documented + layout. +2. Zip: `cd research && zip -r quakeguard_mems_v1.3.0.zip data/` +3. Upload to **Zenodo** → metadata (title, authors + ORCID, license, keywords) + → **Publish**. +4. Append the returned DOI (type `doi`) to `CITATION.cff`. +5. Do **not** commit the raw data (gitignored); the README documents re-download. + +> The DOI is assigned **once**; fix data only by publishing a new Zenodo version. \ No newline at end of file diff --git a/research/__init__.py b/research/__init__.py new file mode 100644 index 0000000..05730d0 --- /dev/null +++ b/research/__init__.py @@ -0,0 +1,13 @@ +"""QuakeGuard SIL research package. + +Software-in-the-Loop cross-validation of the STA/LTA detection core +(R1 of the ROADMAP #Research node). + +Pipeline: + fetch_itaca.py -> download INGV/ITACA accelerograms + ground truth + orchestrator.py -> run the host-compiled detect_cli (same C++ as firmware) + metrics.py -> Sensitivity/Recall, False-Alarm Rate, latency, ROC + calibrate.py -> sweep TRIGGER_RATIO x NOISE_FLOOR against ground truth +""" + +__version__ = "0.1.0" diff --git a/research/calibrate.py b/research/calibrate.py new file mode 100644 index 0000000..5def135 --- /dev/null +++ b/research/calibrate.py @@ -0,0 +1,129 @@ +"""Calibration of the STA/LTA trigger parameters against INGV ground truth. + +Sweeps TRIGGER_RATIO x NOISE_FLOOR (HPF_ALPHA fixed at 0.9), runs the +detector on every accelerogram, and selects the parameters that maximise the +F1 score (harmonic mean of precision and recall). + +Output (JSON): + { + "best": {"trigger_ratio": ..., "noise_floor": ..., "hpf_alpha": 0.9, + "sensitivity": ..., "false_alarm_rate": ..., + "median_latency_s": ...}, + "sweep": [{ "trigger_ratio":..., "noise_floor":..., + "sensitivity":..., "false_alarm_rate":..., "f1":...}, ...], + "roc": [{ "trigger_ratio":..., "fpr":..., "tpr":..., + "median_latency_s":...}, ...] + } +""" + +from __future__ import annotations + +import argparse +import json +from pathlib import Path + +from calibrate_io import load_validation_set +from metrics import compute_metrics, roc_curve +from orchestrator import build_cli, run_detector + +DEFAULT_RATIOS = [1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0] +DEFAULT_FLOORS = [0.02, 0.03, 0.04, 0.05, 0.06, 0.08] +HPF_ALPHA = 0.9 + + +def f1(sensitivity: float, false_alarm_rate: float) -> float: + """F1 = harmonic mean of precision and recall.""" + precision = 1.0 - false_alarm_rate + denom = precision + sensitivity + if denom == 0.0: + return 0.0 + return 2.0 * (precision * sensitivity) / denom + + +def calibrate( + cli_path: Path, + samples, + ground_truth, + ratios: list[float] | None = None, + floors: list[float] | None = None, +) -> dict: + """Run the full sweep and return the summary dict (JSON-serialisable).""" + ratios = ratios or DEFAULT_RATIOS + floors = floors or DEFAULT_FLOORS + + def evaluate(ratio, floor): + results = {} + for s in samples: + results[s.event_id] = run_detector( + cli_path, + s.times, + s.axes, + trigger_ratio=ratio, + noise_floor=floor, + hpf_alpha=HPF_ALPHA, + ) + return results + + sweep = [] + for ratio in ratios: + for floor in floors: + m = compute_metrics(evaluate(ratio, floor), ground_truth) + sweep.append( + { + "trigger_ratio": ratio, + "noise_floor": floor, + "sensitivity": m.sensitivity, + "false_alarm_rate": m.false_alarm_rate, + "median_latency_s": m.median_latency_s, + "f1": f1(m.sensitivity, m.false_alarm_rate), + } + ) + + best = max(sweep, key=lambda row: row["f1"]) + roc = [] + for ratio in ratios: + floor = best["noise_floor"] + m = compute_metrics(evaluate(ratio, floor), ground_truth) + roc.append( + { + "trigger_ratio": ratio, + "fpr": m.false_alarm_rate, + "tpr": m.sensitivity, + "median_latency_s": m.median_latency_s, + } + ) + + return {"best": best, "sweep": sweep, "roc": roc} + + +def main(argv: list[str] | None = None) -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "dataset", type=Path, help="validation dataset dir (CSV + ground truth JSON)" + ) + parser.add_argument("--out", type=Path, default=Path("calibration.json")) + parser.add_argument("--cli", type=Path, default=None, help="detect_cli binary") + args = parser.parse_args(argv) + + samples, ground_truth = load_validation_set(args.dataset) + cli = build_cli(args.cli) + + summary = calibrate(cli, samples, ground_truth) + args.out.write_text(json.dumps(summary, indent=2)) + print(f"Calibration written to {args.out}") + print( + f"Best: ratio={summary['best']['trigger_ratio']} " + f"floor={summary['best']['noise_floor']} " + f"sensitivity={summary['best']['sensitivity']:.3f} " + f"FAR={summary['best']['false_alarm_rate']:.3f} " + f"latency={_fmt_latency(summary['best']['median_latency_s'])}" + ) + + +def _fmt_latency(latency: float | None) -> str: + """Render a latency, or a clear placeholder when no event was detected.""" + return f"latency={latency:.3f}s" if latency is not None else "latency=n/a (no detections)" + + +if __name__ == "__main__": + main() diff --git a/research/calibrate_io.py b/research/calibrate_io.py new file mode 100644 index 0000000..5208b80 --- /dev/null +++ b/research/calibrate_io.py @@ -0,0 +1,67 @@ +"""Dataset I/O for the SIL validation set. + +Shared layout (documented in research/README.md): + / + events/ + .csv # t,ax,ay,az (100 Hz, G units) + ground_truth.json # [{event_id, p_arrival_s}] +""" + +from __future__ import annotations + +import csv +import json +from dataclasses import dataclass, field +from pathlib import Path + +from metrics import GroundTruth + +GROUND_TRUTH_FILE = "ground_truth.json" +EVENTS_DIR = "events" + + +@dataclass +class Accelerogram: + """One recorded accelerogram.""" + + event_id: str + times: list[float] = field(default_factory=list) + axes: list[tuple[float, float, float]] = field(default_factory=list) + + +def load_validation_set(dataset_dir: Path) -> tuple[list[Accelerogram], list[GroundTruth]]: + """Load all accelerograms and ground truth from a validation dataset dir.""" + dataset_dir = Path(dataset_dir) + events_dir = dataset_dir / EVENTS_DIR + gt_path = dataset_dir / GROUND_TRUTH_FILE + + if not events_dir.is_dir(): + raise FileNotFoundError(f"No {EVENTS_DIR}/ directory in {dataset_dir}") + if not gt_path.is_file(): + raise FileNotFoundError(f"Missing {GROUND_TRUTH_FILE} in {dataset_dir}") + + with gt_path.open() as f: + raw_gt = json.load(f) + ground_truth = [GroundTruth(event_id=g["event_id"], p_arrival_s=float(g["p_arrival_s"])) for g in raw_gt] + + samples: list[Accelerogram] = [] + for csv_path in sorted(events_dir.glob("*.csv")): + times, axes = read_accelerogram(csv_path) + samples.append(Accelerogram(event_id=csv_path.stem, times=times, axes=axes)) + + return samples, ground_truth + + +def read_accelerogram(path: Path) -> tuple[list[float], list[tuple[float, float, float]]]: + """Read a t,ax,ay,az CSV (skipping '#' comments).""" + times: list[float] = [] + axes: list[tuple[float, float, float]] = [] + with path.open() as f: + reader = csv.reader(f) + for row in reader: + if not row or row[0].startswith("#"): + continue + t, ax, ay, az = (float(v) for v in row[:4]) + times.append(t) + axes.append((ax, ay, az)) + return times, axes diff --git a/research/fetch_itaca.py b/research/fetch_itaca.py new file mode 100644 index 0000000..e349669 --- /dev/null +++ b/research/fetch_itaca.py @@ -0,0 +1,254 @@ +"""Download accelerograms + ground truth from the ITACA strong-motion portal. + +Graceful-degradation fetcher used as the single entry point for real-world +validation data for the SIL pipeline (ROADMAP R1). The script emits a *fixed* +dataset layout (documented in research/README.md) so downstream modules +(metrics.py, calibrate.py, the C++ core) never know whether they process a real +earthquake or a locally generated mock -- that is the I/O contract. + +Resolution strategy: + 1. If the ``ITACA_TOKEN`` environment variable is set (e.g. from a .env + file), fetch the real accelerogram from the ITACA/ESM ``eventdata`` + web-service and parse the returned DYNA 1.2 ASCII archive. + 2. Otherwise, gracefully degrade to a *realistic synthetic* accelerogram + (white background noise + high-frequency P impulse + larger S arrival), + whose exact P-arrival is known and written to ground_truth.json. + +Output: + / + events/.csv # t,ax,ay,az (100 Hz, G) + ground_truth.json # [{event_id, p_arrival_s}] +""" + +from __future__ import annotations + +import argparse +import csv +import json +import math +import os +import random +import sys +from pathlib import Path + +SAMPLING_HZ = 100 + +# The gravity constant is the *standard* one used in geophysics, NOT the 9.81 +# approximation from the ADXL345 datasheet. ITACA/ESM DYNA 1.2 files carry +# accelerations in Gal (cm/s^2); 1 G = 980.665 Gal exactly. Keeping this +# constant lets the C++ NOISE_FLOOR react to the same physical scale both +# locally and on the real MEMS node. +G_TO_MS2 = 9.80665 +GAL_PER_G = 1e2 / G_TO_MS2 # = 980.665 Gal per G +GRAVITY = G_TO_MS2 # m/s^2, vertical baseline for the synthetic mock + + +class ItacaDataError(RuntimeError): + """Raised when the ITACA portal does not return usable data.""" + + +class ItacaFetcher: + """Downloads real ITACA data when a token is available. + + Uses the ITACA strong-motion web-services: + eventdata WS : /itaca40ws/eventdata/1/query + auth : a signed-message token from + /itaca40ws/generate-signed-message/1/ + + The exact HTTP auth handshake changes across portal revisions and is + therefore isolated here so the rest of the pipeline never changes. + """ + + def __init__(self, base_url: str = "https://itaca.mi.ingv.it", token: str | None = None) -> None: + self.base_url = base_url.rstrip("/") + self.token = token if token is not None else os.environ.get("ITACA_TOKEN") + + @property + def available(self) -> bool: + """Whether a real download is possible (a token is configured).""" + return bool(self.token) + + def list_events(self, min_magnitude: float = 4.0) -> list[dict]: + """Return [{event_id, p_arrival_s}] for catalog events above magnitude. + + Queries the ITACA *flatfile* web-service (publicly accessible, no token + required) for candidate events and their P-arrival metadata. The + flatfile currently does NOT expose phase arrivals, so this is a hook + for the future schema mapping. + """ + raise NotImplementedError( + "Real ITACA event catalogue not bound to the live flatfile schema; " + "use the synthetic fallback (no ITACA_TOKEN)." + ) + + def download(self, event: dict, out_dir: Path) -> None: + """Download one real event's accelerogram and write the CSV (G units). + + Adapt `_build_request` / `_parse_dyna` to the live portal format. + """ + raise NotImplementedError( + "Real ITACA waveform parsing is token-gated and not implemented; " + "use the synthetic fallback (no ITACA_TOKEN)." + ) + + +class RealisticSynthetic: + """Realistic accelerometer-like mock: noise + high-frequency P + larger S. + + The waveform mimics what a MEMS node resting under gravity records. Because + the detector's feature is the *magnitude* norm(x,y,z) over a ~9.8 m/s^2 + gravity baseline, only the VERTICAL (Z) component changes that norm in + first order; horizontal shaking perturbs it to second order and stays below + the noise floor. The synthetic therefore puts both the P and S energy on + the vertical axis (as a vertical-component accelerogram does), while the + measured axes still carry a realistic 9.8 m/s^2 gravity offset: + + - continuous white background noise (well below NOISE_FLOOR); + - a short, high-frequency P impulse arriving at ``p_arrival_s``; + - a stronger, lower-frequency S arrival ~1.8 s later; + - a constant GRAVITY vertical offset (the node "sits" under gravity). + + Units: m/s^2, matching the Adafruit `sensors_event_t` API the firmware + feeds into the core and the 9.8 m/s^2 gravity baseline of the C++ detector. + """ + + def __init__(self, fs: int = SAMPLING_HZ) -> None: + self.fs = fs + + def synthesize( + self, + event_id: str, + seed: int = 0, + pga_ms2: float = 2.0, + p_arrival_s: float = 7.0, + ) -> tuple[list[float], list[tuple[float, float, float]], float]: + """Return (times, axes_in_m/s2, exact_p_arrival_s).""" + del event_id + rng = random.Random(seed) + fs = self.fs + + s_offset_s = 1.8 # S follows P by ~1.8 s (near-source delay) + s_arrival_s = p_arrival_s + s_offset_s + duration_s = p_arrival_s + 15.0 + n = int(duration_s * fs) + times = [i / fs for i in range(n)] + + noise_sigma = 0.02 # ~0.02 m/s^2 background gaussian noise (<= 2e-3 G) + axes: list[tuple[float, float, float]] = [] + + for i in range(n): + t = times[i] + ax = rng.gauss(0.0, noise_sigma) + ay = rng.gauss(0.0, noise_sigma) + raw_z = rng.gauss(0.0, noise_sigma) + + # P-wave impulse: short high-frequency burst on the vertical axis + dt_p = t - p_arrival_s + if 0.0 <= dt_p < 0.5: + f_p = 10.0 + az_p = 0.30 * pga_ms2 * math.sin(2 * math.pi * f_p * dt_p) * math.exp(-25.0 * dt_p) + raw_z += az_p + ax += 0.15 * az_p + ay += 0.15 * az_p + + # S-wave arrival: larger, lower-frequency energy on the vertical axis + dt_s = t - s_arrival_s + if 0.0 <= dt_s < 5.0: + f_s = 2.5 + as_ = pga_ms2 * math.sin(2 * math.pi * f_s * dt_s) * math.exp(-1.2 * dt_s) + raw_z += as_ + + # Horizontal projection is kept near-zero: with a gravity baseline the + # magnitude is dominated by the vertical, so horizontal-only motion is + # invisible to norm3 (as on the real node). + axes.append((ax, ay, raw_z + GRAVITY)) + + return times, axes, p_arrival_s + + +def write_accelerogram_csv( + path: Path, times: list[float], axes: list[tuple[float, float, float]] +) -> None: + """Write an accelerogram in the shared t,ax,ay,az CSV format (G units).""" + with path.open("w", newline="") as f: + writer = csv.writer(f) + writer.writerow(["# t,ax,ay,az"]) # header (ignored by parser) + for t, (ax, ay, az) in zip(times, axes): + writer.writerow([f"{t:.6f}", f"{ax:.6f}", f"{ay:.6f}", f"{az:.6f}"]) + + +def build_synthetic_dataset(out_dir: Path, n_events: int = 5, seed: int = 42) -> None: + """Generate a realistic synthetic validation dataset (fallback mode).""" + out_dir = Path(out_dir) + events_dir = out_dir / "events" + events_dir.mkdir(parents=True, exist_ok=True) + + gen = RealisticSynthetic() + ground_truth: list[dict] = [] + for i in range(n_events): + event_id = f"synth_{i:03d}" + times, axes, p_arrival = gen.synthesize(event_id, seed=seed + i, pga_ms2=1.0 + 0.6 * i) + write_accelerogram_csv(events_dir / f"{event_id}.csv", times, axes) + ground_truth.append({"event_id": event_id, "p_arrival_s": p_arrival}) + + with (out_dir / "ground_truth.json").open("w") as f: + json.dump(ground_truth, f, indent=2) + + +def resolve_mode(fetcher: ItacaFetcher | None = None) -> str: + """Return ``real`` or ``synthetic`` depending on the available token.""" + fetcher = fetcher or ItacaFetcher() + return "real" if fetcher.available else "synthetic" + + +def download_catalog( + out_dir: Path, + min_magnitude: float = 4.0, + n_events: int = 5, + seed: int = 42, + fetcher: ItacaFetcher | None = None, +) -> tuple[list[Path], list[dict], str]: + """Download the catalog. Returns (written_paths, ground_truth, mode). + + The output contract (layout + units) is identical in both modes. + """ + mode = resolve_mode(fetcher) + if mode == "real": + # TODO: implement the real ITACA path once a token is available and the + # DYNA 1.2 parser is bound. Until then the real path refuses to run + # instead of emitting a mock under a misleading name. + raise NotImplementedError( + "ITACA real download is not implemented yet. " + "Run without ITACA_TOKEN to use the synthetic fallback." + ) + + build_synthetic_dataset(Path(out_dir), n_events=n_events, seed=seed) + written = sorted((Path(out_dir) / "events").glob("*.csv")) + with (Path(out_dir) / "ground_truth.json").open() as f: + ground_truth = json.load(f) + return written, ground_truth, mode + + +def main(argv: list[str] | None = None) -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("out_dir", type=Path, help="output validation dataset dir") + parser.add_argument("--min-magnitude", type=float, default=4.0) + parser.add_argument("--n-events", type=int, default=5) + parser.add_argument("--seed", type=int, default=42) + args = parser.parse_args(argv) + + try: + written, ground_truth, mode = download_catalog( + args.out_dir, args.min_magnitude, n_events=args.n_events, seed=args.seed + ) + except NotImplementedError as exc: + print(f"ITACA download not yet integrated ({exc}).", file=sys.stderr) + return 1 + + print(f"[{mode}] Wrote {len(written)} accelerograms to {args.out_dir / 'events'}") + print(f"[{mode}] Ground truth: {len(ground_truth)} events") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/research/metrics.py b/research/metrics.py new file mode 100644 index 0000000..625be3f --- /dev/null +++ b/research/metrics.py @@ -0,0 +1,147 @@ +"""SIL detection metrics: Sensitivity/Recall, False-Alarm Rate, latency, ROC. + +Definitions (per event/accelerogram): + - Sensitivity/Recall = TP / (TP + FN) over events + - False-Alarm Rate = FP / (FP + TP) over triggers + - Response latency = trigger time - P-arrival time (median over TPs) + - ROC curve = sweep of TRIGGER_RATIO -> (FPR, TPR) + +A trigger is a True Positive (TP) if it falls within a matching window after +the ground-truth P-arrival; otherwise it is a False Positive (FP). +""" + +from __future__ import annotations + +import statistics +from dataclasses import dataclass, field + +from orchestrator import TriggerPoint + +# P-wave arrival tolerance: a trigger is a TP if it fires between P and +# P + WINDOW. Triggers before the event or far after are FPs. +WINDOW_S = 10.0 + + +@dataclass +class GroundTruth: + """Known arrival times for one event.""" + + event_id: str + p_arrival_s: float # P-wave arrival, relative to the start of the accelerogram + + +@dataclass +class EventScore: + """Per-event outcome.""" + + event_id: str + triggers: list[TriggerPoint] + p_arrival_s: float + detected: bool + latency_s: float | None = None + + +@dataclass +class MetricsSummary: + """Aggregate metrics over a validation set.""" + + sensitivity: float + false_alarm_rate: float + median_latency_s: float | None + n_events: int + n_triggers: int + n_true_positives: int + n_false_positives: int + per_event: list[EventScore] = field(default_factory=list) + + +def score_event(triggers: list[TriggerPoint], truth: GroundTruth) -> EventScore: + """Classify one event's triggers against its P-arrival.""" + tp = [ + tr + for tr in triggers + if tr.time_s >= truth.p_arrival_s and tr.time_s <= truth.p_arrival_s + WINDOW_S + ] + detected = len(tp) > 0 + latency = min(tr.time_s - truth.p_arrival_s for tr in tp) if tp else None + return EventScore( + event_id=truth.event_id, + triggers=triggers, + p_arrival_s=truth.p_arrival_s, + detected=detected, + latency_s=latency, + ) + + +def compute_metrics( + results: dict[str, list[TriggerPoint]], ground_truth: list[GroundTruth] +) -> MetricsSummary: + """Compute aggregate sensitivity, FAR and latency over the whole set.""" + truth_by_id = {g.event_id: g for g in ground_truth} + unknown = [e for e in results if e not in truth_by_id] + if unknown: + raise ValueError(f"No ground truth for events: {sorted(unknown)}") + + scores: list[EventScore] = [] + for event_id, triggers in results.items(): + scores.append(score_event(triggers, truth_by_id[event_id])) + + n_events = len(scores) + n_detected = sum(1 for s in scores if s.detected) + n_triggers = sum(len(s.triggers) for s in scores) + n_tp = sum(1 for s in scores if s.detected) + n_fp = n_triggers - n_tp + + latencies = [s.latency_s for s in scores if s.latency_s is not None] + median_latency = statistics.median(latencies) if latencies else None + + sensitivity = n_detected / n_events if n_events else 0.0 + false_alarm_rate = n_fp / (n_fp + n_tp) if (n_fp + n_tp) else 0.0 + + return MetricsSummary( + sensitivity=sensitivity, + false_alarm_rate=false_alarm_rate, + median_latency_s=median_latency, + n_events=n_events, + n_triggers=n_triggers, + n_true_positives=n_tp, + n_false_positives=n_fp, + per_event=scores, + ) + + +@dataclass +class RocPoint: + """One operating point of the ROC curve.""" + + trigger_ratio: float + fpr: float + tpr: float + median_latency_s: float | None + + +def roc_curve( + evaluate: callable, + ground_truth: list[GroundTruth], + ratios: list[float], +) -> list[RocPoint]: + """Trace the ROC curve by sweeping TRIGGER_RATIO. + + Args: + evaluate: callable(ratio) -> dict[event_id, list[TriggerPoint]]. + ground_truth: ground-truth arrivals. + ratios: candidate TRIGGER_RATIO values (ascending recommended). + """ + curve: list[RocPoint] = [] + for ratio in ratios: + results = evaluate(ratio) + m = compute_metrics(results, ground_truth) + curve.append( + RocPoint( + trigger_ratio=ratio, + fpr=m.false_alarm_rate, + tpr=m.sensitivity, + median_latency_s=m.median_latency_s, + ) + ) + return curve diff --git a/research/orchestrator.py b/research/orchestrator.py new file mode 100644 index 0000000..694dcdd --- /dev/null +++ b/research/orchestrator.py @@ -0,0 +1,101 @@ +"""Drive the host-compiled detection binary (same C++ source as the firmware). + +The orchestrator is the *sole* Python↔C++ bridge: it invokes `detect_cli` via +subprocess, feeds it the accelerograms over stdin and collects the trigger +points. The detection decision itself is 100% owned by the C++ core. +""" + +from __future__ import annotations + +import shutil +import subprocess +from dataclasses import dataclass, field +from pathlib import Path + +DEFAULT_CLI = Path(__file__).resolve().parents[1] / "firmware" / "tools" / "detect_cli" + + +@dataclass +class TriggerPoint: + """A single trigger emitted by the detector.""" + + time_s: float # seconds, relative to the start of the accelerogram + ratio: float # STA/LTA at trigger time + + +@dataclass +class DetectionResult: + """Detector output for one accelerogram.""" + + event_id: str + triggers: list[TriggerPoint] = field(default_factory=list) + + +class DetectionError(RuntimeError): + """Raised when the C++ binary fails to run.""" + + +def build_cli(cli_path: Path | None = None) -> Path: + """Compile the host CLI if the binary is missing. Returns the binary path.""" + cli_path = cli_path or DEFAULT_CLI + if cli_path.is_file(): + return cli_path + + src_dir = Path(__file__).resolve().parents[1] / "firmware" + cli_src = src_dir / "tools" / "detect_cli.cpp" + if not cli_src.is_file(): + raise DetectionError(f"Missing source: {cli_src}") + + gcc = shutil.which("g++") + if gcc is None: + raise DetectionError("g++ not found: install a C++ toolchain to run SIL validation") + + cmd = [gcc, "-std=c++11", "-I", str(src_dir / "src"), str(cli_src), "-lm", "-o", str(cli_path)] + subprocess.run(cmd, check=True, capture_output=True) + return cli_path + + +def run_detector( + cli_path: Path, + times: list[float], + axes: list[tuple[float, float, float]], + trigger_ratio: float = 1.8, + noise_floor: float = 0.04, + hpf_alpha: float = 0.9, +) -> list[TriggerPoint]: + """Feed one accelerogram to the C++ core and return the trigger points. + + Args: + cli_path: path to the compiled detect_cli binary. + times: sample timestamps in seconds (100 Hz). + axes: (ax, ay, az) samples in G. + trigger_ratio/noise_floor/hpf_alpha: detector parameters. + """ + if len(times) != len(axes): + raise ValueError("times and axes must have the same length") + + lines = ["# t,ax,ay,az"] + for t, (ax, ay, az) in zip(times, axes): + lines.append(f"{t:.6f},{ax:.6f},{ay:.6f},{az:.6f}") + stdin_data = "\n".join(lines) + + cmd = [str(cli_path), str(trigger_ratio), str(noise_floor), str(hpf_alpha)] + proc = subprocess.run( + cmd, input=stdin_data, capture_output=True, text=True, check=False + ) + if proc.returncode != 0: + raise DetectionError(f"detect_cli exited {proc.returncode}: {proc.stderr}") + + triggers: list[TriggerPoint] = [] + for line in proc.stdout.splitlines(): + line = line.strip() + if not line or line.startswith("#"): + continue + parts = line.split(",") + if len(parts) != 2: + continue + try: + triggers.append(TriggerPoint(time_s=float(parts[0]), ratio=float(parts[1]))) + except ValueError: + continue + return triggers diff --git a/research/plot_roc.py b/research/plot_roc.py new file mode 100644 index 0000000..ab98a8c --- /dev/null +++ b/research/plot_roc.py @@ -0,0 +1,53 @@ +"""Generate paper artifacts from a calibration run: ROC plot + metrics JSON. + +matplotlib is optional; plotting is skipped cleanly if it is not installed. +""" + +from __future__ import annotations + +import argparse +import json +from pathlib import Path + + +def plot_roc(calibration: dict, out_path: Path) -> bool: + """Render the ROC curve. Returns False if matplotlib is unavailable.""" + try: + import matplotlib + + matplotlib.use("Agg") + import matplotlib.pyplot as plt + except ImportError: + return False + + roc = calibration["roc"] + xs = [r["fpr"] for r in roc] + ys = [r["tpr"] for r in roc] + + fig, ax = plt.subplots(figsize=(6, 6)) + ax.plot(xs, ys, marker="o", label="SIL detector") + ax.plot([0, 1], [0, 1], "--", color="gray", label="random") + ax.set_xlabel("False-Alarm Rate") + ax.set_ylabel("Sensitivity (TPR)") + ax.set_title("QuakeGuard R1: STA/LTA Algorithm Cross-Validation (SIL)") + ax.legend(loc="lower right") + ax.set_xlim(0, 1) + ax.set_ylim(0, 1) + fig.tight_layout() + fig.savefig(out_path, dpi=200) + return True + + +def main(argv: list[str] | None = None) -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("calibration_json", type=Path) + parser.add_argument("--roc-out", type=Path, default=Path("roc_curve.png")) + args = parser.parse_args(argv) + + calibration = json.loads(args.calibration_json.read_text()) + wrote = plot_roc(calibration, args.roc_out) + print(f"ROC plot written to {args.roc_out}" if wrote else "matplotlib not installed; plot skipped") + + +if __name__ == "__main__": + main() diff --git a/research/requirements.txt b/research/requirements.txt new file mode 100644 index 0000000..a13072b --- /dev/null +++ b/research/requirements.txt @@ -0,0 +1,3 @@ +# QuakeGuard SIL research (R1) — runtime dependencies (runtime only). +# Plotting requires matplotlib; everything else uses the Python stdlib. +matplotlib>=3.8 \ No newline at end of file diff --git a/research/synthetic.py b/research/synthetic.py new file mode 100644 index 0000000..5797728 --- /dev/null +++ b/research/synthetic.py @@ -0,0 +1,33 @@ +"""Generate a synthetic validation dataset for local smoke tests / CI. + +Thin wrapper over the realistic synthetic generator in fetch_itaca.py, so the +fallback (no-token) path and this generator stay exactly in sync. Produces the +same layout as the ITACA path: + / + events/.csv + ground_truth.json + +Used to exercise the full SIL pipeline without downloading the INGV dataset. +""" + +from __future__ import annotations + +import argparse +from pathlib import Path + +from fetch_itaca import build_synthetic_dataset + + +def main(argv: list[str] | None = None) -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("out_dir", type=Path, help="output dataset dir") + parser.add_argument("--n-events", type=int, default=5) + parser.add_argument("--seed", type=int, default=42) + args = parser.parse_args(argv) + + build_synthetic_dataset(args.out_dir, n_events=args.n_events, seed=args.seed) + print(f"Synthetic dataset written to {args.out_dir}") + + +if __name__ == "__main__": + main() \ No newline at end of file From 4ae4a9abc2b2db5f90754983e364ed9632a57803 Mon Sep 17 00:00:00 2001 From: GiZano Date: Thu, 6 Aug 2026 17:00:43 +0200 Subject: [PATCH 2/4] fix(research): address Sonar lint in R1 pipeline - DetectionCore.h: explicit constructor + member init list (S1709/S3230) - detect_cli.cpp: std::string/std::getline, dedicated declarations, auto (S5945/S1659/S5827) - calibrate.py: math.isclose for float equality (S1244) - calibrate_io.py/fetch_itaca.py/calibrate.py/plot_roc.py: path-injection guard resolve_within_root (S8707) - fetch_itaca.py: drop unused param, int-typed main (S1172/S5886/S3699/S1135) - ci: run SIL smoke test inside gitignored research/out (respects path guard) --- .github/workflows/iot-ci.yml | 4 ++-- firmware/src/DetectionCore.h | 23 ++++++++++++----------- firmware/tools/detect_cli.cpp | 21 +++++++++++++-------- research/calibrate.py | 10 ++++++---- research/calibrate_io.py | 16 ++++++++++++++++ research/fetch_itaca.py | 19 ++++++++++--------- research/plot_roc.py | 7 +++++-- 7 files changed, 64 insertions(+), 36 deletions(-) diff --git a/.github/workflows/iot-ci.yml b/.github/workflows/iot-ci.yml index 52e897a..f0780a5 100644 --- a/.github/workflows/iot-ci.yml +++ b/.github/workflows/iot-ci.yml @@ -66,5 +66,5 @@ jobs: - name: SIL Smoke Test (synthetic dataset) run: | - python research/synthetic.py /tmp/synth_dataset --n-events 3 - cd research && python calibrate.py /tmp/synth_dataset --out /tmp/calibration.json \ No newline at end of file + python research/synthetic.py research/out/synth_dataset --n-events 3 + cd research && python calibrate.py out/synth_dataset --out out/calibration.json \ No newline at end of file diff --git a/firmware/src/DetectionCore.h b/firmware/src/DetectionCore.h index f8d134d..22994d7 100644 --- a/firmware/src/DetectionCore.h +++ b/firmware/src/DetectionCore.h @@ -26,9 +26,9 @@ class SeismicDetector { static constexpr unsigned long COOLDOWN_MS = 5000; static constexpr float INITIAL_RAW = 9.81f; // gravity baseline - SeismicDetector(float hpfAlpha = DEFAULT_HPF_ALPHA, - float triggerRatio = DEFAULT_TRIGGER_RATIO, - float noiseFloor = DEFAULT_NOISE_FLOOR); + explicit SeismicDetector(float hpfAlpha = DEFAULT_HPF_ALPHA, + float triggerRatio = DEFAULT_TRIGGER_RATIO, + float noiseFloor = DEFAULT_NOISE_FLOOR); /** * Euclidean norm of a 3-axis acceleration sample (the single shared @@ -72,14 +72,15 @@ class SeismicDetector { // --------------------------------------------------------------------------- inline SeismicDetector::SeismicDetector(float hpfAlpha, float triggerRatio, float noiseFloor) - : hpfAlpha_(hpfAlpha), triggerRatio_(triggerRatio), noiseFloor_(noiseFloor) { - filtered_ = 0.0f; - prevRaw_ = INITIAL_RAW; - inAlarm_ = false; - alarmStartMs_ = 0; - sta_ = 0.0f; - ratio_ = 0.0f; -} + : hpfAlpha_(hpfAlpha), + triggerRatio_(triggerRatio), + noiseFloor_(noiseFloor), + filtered_(0.0f), + prevRaw_(INITIAL_RAW), + inAlarm_(false), + alarmStartMs_(0), + sta_(0.0f), + ratio_(0.0f) {} inline bool SeismicDetector::push(float rawMag, unsigned long nowMs) { // High-pass filter to remove gravity. NOTE: uses std::abs; on-firmware the diff --git a/firmware/tools/detect_cli.cpp b/firmware/tools/detect_cli.cpp index b23906d..38867e5 100644 --- a/firmware/tools/detect_cli.cpp +++ b/firmware/tools/detect_cli.cpp @@ -18,6 +18,8 @@ #include #include #include +#include +#include #include "DetectionCore.h" int main(int argc, char** argv) { @@ -31,16 +33,19 @@ int main(int argc, char** argv) { SeismicDetector det(hpfAlpha, triggerRatio, noiseFloor); - char line[256]; - while (fgets(line, sizeof(line), stdin) != nullptr) { - if (line[0] == '#') continue; + std::string line; + while (std::getline(std::cin, line)) { + if (line.empty() || line[0] == '#') continue; - double t, ax, ay, az; - if (sscanf(line, "%lf,%lf,%lf,%lf", &t, &ax, &ay, &az) != 4) continue; + double t = 0.0; + double ax = 0.0; + double ay = 0.0; + double az = 0.0; + if (std::sscanf(line.c_str(), "%lf,%lf,%lf,%lf", &t, &ax, &ay, &az) != 4) continue; - float raw = SeismicDetector::norm3(static_cast(ax), - static_cast(ay), - static_cast(az)); + auto raw = SeismicDetector::norm3(static_cast(ax), + static_cast(ay), + static_cast(az)); // 100 Hz => sample clock = t * 1000 ms (matches firmware millis()). unsigned long nowMs = static_cast(t * 1000.0); diff --git a/research/calibrate.py b/research/calibrate.py index 5def135..dec694d 100644 --- a/research/calibrate.py +++ b/research/calibrate.py @@ -20,9 +20,10 @@ import argparse import json +import math from pathlib import Path -from calibrate_io import load_validation_set +from calibrate_io import load_validation_set, resolve_within_root from metrics import compute_metrics, roc_curve from orchestrator import build_cli, run_detector @@ -35,7 +36,7 @@ def f1(sensitivity: float, false_alarm_rate: float) -> float: """F1 = harmonic mean of precision and recall.""" precision = 1.0 - false_alarm_rate denom = precision + sensitivity - if denom == 0.0: + if math.isclose(denom, 0.0, abs_tol=1e-12): return 0.0 return 2.0 * (precision * sensitivity) / denom @@ -109,8 +110,9 @@ def main(argv: list[str] | None = None) -> None: cli = build_cli(args.cli) summary = calibrate(cli, samples, ground_truth) - args.out.write_text(json.dumps(summary, indent=2)) - print(f"Calibration written to {args.out}") + out_path = resolve_within_root(args.out) + out_path.write_text(json.dumps(summary, indent=2)) + print(f"Calibration written to {out_path}") print( f"Best: ratio={summary['best']['trigger_ratio']} " f"floor={summary['best']['noise_floor']} " diff --git a/research/calibrate_io.py b/research/calibrate_io.py index 5208b80..76f3eab 100644 --- a/research/calibrate_io.py +++ b/research/calibrate_io.py @@ -20,6 +20,22 @@ EVENTS_DIR = "events" +def resolve_within_root(path: Path, root: Path | None = None) -> Path: + """Resolve a user-supplied output path and refuse escapes outside *root*. + + Guards the CLI scripts against path-injection (Sonar S8707): a path such as + ``../../etc/something`` constructed from a command-line argument must not be + allowed to write outside the intended working directory. + """ + base = (root or Path.cwd()).resolve() + target = Path(path).expanduser().resolve() + try: + target.relative_to(base) + except ValueError: + raise ValueError(f"Refusing to write outside {base}: {target}") from None + return target + + @dataclass class Accelerogram: """One recorded accelerogram.""" diff --git a/research/fetch_itaca.py b/research/fetch_itaca.py index e349669..da96cdd 100644 --- a/research/fetch_itaca.py +++ b/research/fetch_itaca.py @@ -31,6 +31,8 @@ import sys from pathlib import Path +from calibrate_io import resolve_within_root + SAMPLING_HZ = 100 # The gravity constant is the *standard* one used in geophysics, NOT the 9.81 @@ -179,8 +181,8 @@ def write_accelerogram_csv( def build_synthetic_dataset(out_dir: Path, n_events: int = 5, seed: int = 42) -> None: """Generate a realistic synthetic validation dataset (fallback mode).""" - out_dir = Path(out_dir) - events_dir = out_dir / "events" + safe_dir = resolve_within_root(out_dir) + events_dir = safe_dir / "events" events_dir.mkdir(parents=True, exist_ok=True) gen = RealisticSynthetic() @@ -191,7 +193,7 @@ def build_synthetic_dataset(out_dir: Path, n_events: int = 5, seed: int = 42) -> write_accelerogram_csv(events_dir / f"{event_id}.csv", times, axes) ground_truth.append({"event_id": event_id, "p_arrival_s": p_arrival}) - with (out_dir / "ground_truth.json").open("w") as f: + with (safe_dir / "ground_truth.json").open("w") as f: json.dump(ground_truth, f, indent=2) @@ -203,7 +205,6 @@ def resolve_mode(fetcher: ItacaFetcher | None = None) -> str: def download_catalog( out_dir: Path, - min_magnitude: float = 4.0, n_events: int = 5, seed: int = 42, fetcher: ItacaFetcher | None = None, @@ -214,9 +215,9 @@ def download_catalog( """ mode = resolve_mode(fetcher) if mode == "real": - # TODO: implement the real ITACA path once a token is available and the - # DYNA 1.2 parser is bound. Until then the real path refuses to run - # instead of emitting a mock under a misleading name. + # The real ITACA path requires a token and a bound DYNA 1.2 parser. + # Until then it refuses to run instead of emitting a mock under a + # misleading name (see research/README.md). raise NotImplementedError( "ITACA real download is not implemented yet. " "Run without ITACA_TOKEN to use the synthetic fallback." @@ -229,7 +230,7 @@ def download_catalog( return written, ground_truth, mode -def main(argv: list[str] | None = None) -> None: +def main(argv: list[str] | None = None) -> int: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("out_dir", type=Path, help="output validation dataset dir") parser.add_argument("--min-magnitude", type=float, default=4.0) @@ -239,7 +240,7 @@ def main(argv: list[str] | None = None) -> None: try: written, ground_truth, mode = download_catalog( - args.out_dir, args.min_magnitude, n_events=args.n_events, seed=args.seed + args.out_dir, n_events=args.n_events, seed=args.seed ) except NotImplementedError as exc: print(f"ITACA download not yet integrated ({exc}).", file=sys.stderr) diff --git a/research/plot_roc.py b/research/plot_roc.py index ab98a8c..c200e28 100644 --- a/research/plot_roc.py +++ b/research/plot_roc.py @@ -9,6 +9,8 @@ import json from pathlib import Path +from calibrate_io import resolve_within_root + def plot_roc(calibration: dict, out_path: Path) -> bool: """Render the ROC curve. Returns False if matplotlib is unavailable.""" @@ -45,8 +47,9 @@ def main(argv: list[str] | None = None) -> None: args = parser.parse_args(argv) calibration = json.loads(args.calibration_json.read_text()) - wrote = plot_roc(calibration, args.roc_out) - print(f"ROC plot written to {args.roc_out}" if wrote else "matplotlib not installed; plot skipped") + out_path = resolve_within_root(args.roc_out) + wrote = plot_roc(calibration, out_path) + print(f"ROC plot written to {out_path}" if wrote else "matplotlib not installed; plot skipped") if __name__ == "__main__": From 5167071d7a0c686b80f228e4ca28bbc2f5d080bd Mon Sep 17 00:00:00 2001 From: GiZano Date: Thu, 6 Aug 2026 17:05:56 +0200 Subject: [PATCH 3/4] fix(research): guard input calibration path too (S8707) --- research/plot_roc.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/research/plot_roc.py b/research/plot_roc.py index c200e28..92c3a15 100644 --- a/research/plot_roc.py +++ b/research/plot_roc.py @@ -46,7 +46,8 @@ def main(argv: list[str] | None = None) -> None: parser.add_argument("--roc-out", type=Path, default=Path("roc_curve.png")) args = parser.parse_args(argv) - calibration = json.loads(args.calibration_json.read_text()) + calibration_path = resolve_within_root(args.calibration_json) + calibration = json.loads(calibration_path.read_text()) out_path = resolve_within_root(args.roc_out) wrote = plot_roc(calibration, out_path) print(f"ROC plot written to {out_path}" if wrote else "matplotlib not installed; plot skipped") From fc0bf50933d744c442b0eea8bacf8eac719a2ce6 Mon Sep 17 00:00:00 2001 From: GiZano Date: Thu, 6 Aug 2026 17:10:15 +0200 Subject: [PATCH 4/4] fix(core): resolve remaining Sonar smells in R1 - DetectionCore.h: drop redundant ctor init list (members use in-class initializers) (S3230) - detect_cli.cpp: use auto for derived type (S5827) --- firmware/src/DetectionCore.h | 10 +--------- firmware/tools/detect_cli.cpp | 2 +- 2 files changed, 2 insertions(+), 10 deletions(-) diff --git a/firmware/src/DetectionCore.h b/firmware/src/DetectionCore.h index 22994d7..e47be8a 100644 --- a/firmware/src/DetectionCore.h +++ b/firmware/src/DetectionCore.h @@ -72,15 +72,7 @@ class SeismicDetector { // --------------------------------------------------------------------------- inline SeismicDetector::SeismicDetector(float hpfAlpha, float triggerRatio, float noiseFloor) - : hpfAlpha_(hpfAlpha), - triggerRatio_(triggerRatio), - noiseFloor_(noiseFloor), - filtered_(0.0f), - prevRaw_(INITIAL_RAW), - inAlarm_(false), - alarmStartMs_(0), - sta_(0.0f), - ratio_(0.0f) {} + : hpfAlpha_(hpfAlpha), triggerRatio_(triggerRatio), noiseFloor_(noiseFloor) {} inline bool SeismicDetector::push(float rawMag, unsigned long nowMs) { // High-pass filter to remove gravity. NOTE: uses std::abs; on-firmware the diff --git a/firmware/tools/detect_cli.cpp b/firmware/tools/detect_cli.cpp index 38867e5..38f2621 100644 --- a/firmware/tools/detect_cli.cpp +++ b/firmware/tools/detect_cli.cpp @@ -48,7 +48,7 @@ int main(int argc, char** argv) { static_cast(az)); // 100 Hz => sample clock = t * 1000 ms (matches firmware millis()). - unsigned long nowMs = static_cast(t * 1000.0); + auto nowMs = static_cast(t * 1000.0); if (det.push(raw, nowMs)) { printf("%.3f,%.6f\n", t, det.lastRatio()); // NOSONAR(cpp:S6494)