diff --git a/.agent/work-plans/issue-92/progress.md b/.agent/work-plans/issue-92/progress.md new file mode 100644 index 0000000..63b2232 --- /dev/null +++ b/.agent/work-plans/issue-92/progress.md @@ -0,0 +1,71 @@ +--- +issue: 92 +--- + +# Issue #92 — Bound resident-tile RAM in the offline importer (persist-then-drop eviction) + +## Local Review (Pre-Push) +**Status**: complete +**When**: 2026-06-29 08:53 +0000 +**By**: Claude Code Agent (Claude Opus) +**Verdict**: changes-requested + +**Branch**: feature/issue-92 at `65dfd9f` +**Mode**: pre-push +**Depth**: Deep (reason: 1502/-319 over 17 files; cross-cutting resource-management + numerical + binary-serialization) +**Must-fix**: 3 | **Suggestions**: 2 +**Round**: 1 | **Ship**: continue — three must-fixes, all mechanical doc/comment/ADR reconciliation (no algorithmic defect); fast convergence expected after addressing. + +### Findings +- [x] (must-fix) `extractNodeRecord` docstring still says intensity is "emitted UNCORRECTED" and "per-beam raw set is retained so it is re-derivable when #15" — both false after the record-time correction + Welford change — `cube_bathymetry/include/cube_bathymetry/node.h:195-204` +- [x] (must-fix) ADR-0007 (and ADR-0001's re-derivability notes) document raw `{raw_intensity, beam_angle}` sufficient-statistic retention + correct-at-extract for cube#15 re-derivability; this PR reverses that decision (record-time correction, O(1) Welford, raw samples discarded) without updating/superseding the ADR — `cube_bathymetry/docs/decisions/0007-mbes-backscatter-store-addendum-phase-b-transition.md` +- [x] (must-fix) Reload-failure path silently drops this batch's soundings on the failed tile and the comment claims it "retries on the next revisit" (only *future* batches retry; these beams are lost permanently offline) — contradicts the headline "lossless" guarantee; fix the comment and log the dropped count (buffer+re-add optional) — `cube_bathymetry/src/store_import.cpp:1503-1535` (cross-confirmed by both adversarial passes) +- [x] (suggestion) `Node::setSettledIntensityWelford` writes to `chooseHypothesis()` without the `number_of_samples > 0` gate its siblings (`chosenIntensityWelford`, `extractNodeRecord`) use — harmless today, add for symmetry/robustness — `cube_bathymetry/src/node.cpp:354-362` +- [x] (suggestion) Scratch-dir name `cube_import_spill__` is not collision-safe across concurrent `import_bag` processes (counter resets per process; ns is the only cross-process distinguisher) → potential silent spill cross-corruption; add `getpid()`/`mkdtemp` as the test helper already does — `cube_bathymetry/src/store_import.cpp:1291-1300` + +## Implementation +**Status**: complete +**When**: 2026-06-29 09:15 +0000 +**By**: Claude Code Agent (Claude Opus) + +**Branch**: feature/issue-92 at `6b9a7d7` +**Addressed**: Local Review (Pre-Push) — When 2026-06-29 08:53 +0000, at `65dfd9f` (Verdict: changes-requested, 3 must-fix + 2 suggestions) +**Commits**: `31ced91`, `4b4caff`, `50c8b56`, `5334319`, `6b9a7d7` + +### Actions +- [x] (must-fix) `extractNodeRecord` docstring no longer claims "emitted UNCORRECTED" / raw-set retained — rewritten for record-time correction + streaming Welford; the matching `intensity`-field docstring (node.h:56-59) corrected too — `cube_bathymetry/include/cube_bathymetry/node.h` (`31ced91`) +- [x] (must-fix) ADR reconciled — added a "Phase B.2 addendum (cube#92/#93)" to ADR-0007 documenting the record-time-correction + Welford reversal and the lost in-place re-derivability; marked the now-false Context/Consequences/Tier-2 claims superseded inline, and updated ADR-0001's backscatter re-derivability trade-off note — `cube_bathymetry/docs/decisions/0007-…md`, `0001-…md` (`4b4caff`) +- [x] (must-fix) Reload-failure path no longer mislabels the loss — `addBatch` now counts soundings centred in un-reloadable tiles and emits a `WARNING` before dropping; the "retries on the next revisit" comments/cerr reworded to state this batch's soundings are not re-processed offline; added `GeoMapSheet::gridIndexForSounding` to support the count — `cube_bathymetry/src/store_import.cpp`, `geo_map_sheet.{h,cpp}` (`50c8b56`) +- [x] (suggestion) `setSettledIntensityWelford` now applies the `number_of_samples > 0` gate its siblings use — `cube_bathymetry/src/node.cpp:354-362` (`5334319`) +- [x] (suggestion) Spill scratch dir now created with `mkdtemp` (atomic, retries on collision) — collision-safe across concurrent `import_bag` processes; dropped the steady-clock-ns + per-process counter naming and the now-unused `` include — `cube_bathymetry/src/store_import.cpp:1291-1300` (`6b9a7d7`) + +**Quick checks**: `ament_cpplint` + `ament_uncrustify` clean on all five touched C++/header files. A full `colcon build` was **not** possible in this worktree — the lower layers (`underlay_ws`/`core_ws`/… providing `marine_autonomy`/`gggs.h`) have no install trees here, so the package cannot be configured. Changes are mechanical (doc/comment + a gate + an additive helper + a `mkdtemp` swap) and were reviewed by hand against current source. + +### Next step +Lifecycle: **Implementation** → **review-code** (re-review the fixes). Hand off to a fresh-context sub-agent: + + .agent/scripts/dispatch_subagent.sh --mode in-process --issue 92 --skill review-code + +## Local Review (Pre-Push) +**Status**: complete +**When**: 2026-06-29 09:30 +0000 +**By**: Claude Code Agent (Claude Opus) +**Verdict**: approved + +**Branch**: feature/issue-92 at `dbe7b21` +**Mode**: pre-push +**Depth**: Deep (reason: ~2000 lines / 20 files; cross-cutting resource-management + numerical estimator + binary serialization) +**Must-fix**: 0 | **Suggestions**: 5 +**Round**: 2 | **Ship**: recommended — re-review of round-1 fixes; all 3 must-fix + 2 suggestions confirmed addressed, no new must-fix, remaining items are doc-accuracy / observability robustness. + +### Findings +- [ ] (suggestion) `restoreSpilledSamples` is silent on a missing/corrupt spill while `persistBackscatterTile` plain-overwrites the tile → out-of-band spill loss (e.g. /tmp reaper on a multi-day import) silently drops pre-eviction backscatter; add a WARN mirroring the bathy reload-failure WARN — `cube_bathymetry/src/store_import.cpp` (`restoreSpilledSamples`/`persistBackscatterTile`) +- [ ] (suggestion) Comment claims "Atomic temp-then-rename via tile_io::saveTile"; verified false — `marine_tiled_raster_store/src/tile_io.cpp:132` writes directly to the final path via GDAL `Create` (no temp/rename). Runtime is still safe (throw keeps tile resident) but the crash-safety reasoning is wrong — fix the comment — `cube_bathymetry/src/store_import.cpp` (`persistBathyTile`) +- [ ] (suggestion) Bolded "Backscatter is lossless under eviction." reads unconditionally; the real guarantee is scoped to a consistent re-survey (a depth-disambiguation divergence between visits could divert backscatter) — tighten the claim — `cube_bathymetry/include/cube_bathymetry/store_import.h` +- [ ] (suggestion) Registries written only at finalize; a mid-pass crash leaves evicted tiles with a source_index but no registry.json — likely acceptable (re-runnable single pass) but undocumented at the eviction site; add a one-line note — `cube_bathymetry/src/store_import.cpp` (`persistBathyTile`/`finalize`) +- [ ] (suggestion) `reloadEvictedTile`'s `loadWindow(sw,ne)` over the tile's own corners relies on loadWindow inclusivity to reload that exact tile (correct per source); add a targeted single-tile reload round-trip test to lock it in — `cube_bathymetry/src/store_import.cpp` (`reloadEvictedTile`) + +**Static analysis**: `ament_cpplint` + `ament_uncrustify` clean on all changed C++/headers/tests. A full `colcon build` was not possible here (lower layers provide `marine_autonomy`/`gggs.h` only via install trees absent in this worktree); the two highest-risk correctness paths were verified by hand against source: (1) the reload window equals the add window exactly (shared `boundsForSoundings` + identical `GridAreaIterator`), so no touched evicted tile is missed; (2) a reload-seeded hypothesis has `number_of_samples == 1`, so the Welford-restore gate never silently skips the restore. + +### Next step +Lifecycle: **Local Review** (approved) → push / open PR → **triage-reviews**. The 5 suggestions are advisory (doc/observability); the operator/host decides whether to apply them before or after push. diff --git a/cube_bathymetry/CMakeLists.txt b/cube_bathymetry/CMakeLists.txt index 1916393..9c932af 100644 --- a/cube_bathymetry/CMakeLists.txt +++ b/cube_bathymetry/CMakeLists.txt @@ -292,6 +292,17 @@ if(BUILD_TESTING) cube_bathymetry_store_import marine_bathymetry_store::marine_bathymetry_store ) + + # Bounded == unbounded equivalence for the offline import accumulator (#92): + # an eviction-forcing budget reproduces the unbounded bathy + backscatter store + # exactly (lossless persist-then-drop + reload-on-revisit), and a persist + # failure never drops a tile. + ament_add_gtest(test_import_eviction test/test_import_eviction.cpp) + target_link_libraries(test_import_eviction + cube_bathymetry_store_import + marine_bathymetry_store::marine_bathymetry_store + marine_mbes_backscatter_store::marine_mbes_backscatter_store + ) endif() # The installed public header store_import.h (#include'd by downstream consumers diff --git a/cube_bathymetry/docs/decisions/0001-tile-eviction-and-incremental-publish.md b/cube_bathymetry/docs/decisions/0001-tile-eviction-and-incremental-publish.md index bd5f00e..0b09ed5 100644 --- a/cube_bathymetry/docs/decisions/0001-tile-eviction-and-incremental-publish.md +++ b/cube_bathymetry/docs/decisions/0001-tile-eviction-and-incremental-publish.md @@ -190,6 +190,10 @@ dependency) rather than retained. the full multi-hypothesis / pre-filter / backscatter-sample state (out of scope; backscatter is re-derivable per cube#15). A revisited cell continues from the persisted best estimate, which is the intended behavior. + **(Updated by cube#92/#93 — see ADR-0007 § Phase B.2 addendum: the corrected + backscatter Welford `(n, mean, M2)` IS now spilled and restored losslessly on + eviction rather than dropped/re-derived; raw per-beam samples are no longer + kept, so a richer correction requires re-importing the bag.)** - A misconfiguration where `max_resident_tiles` < CA-window tile span degrades the live view (over-lethal holes). Mitigated by the on_configure WARN and a default pair chosen to satisfy the constraint. diff --git a/cube_bathymetry/docs/decisions/0007-mbes-backscatter-store-addendum-phase-b-transition.md b/cube_bathymetry/docs/decisions/0007-mbes-backscatter-store-addendum-phase-b-transition.md index 73e7e97..6e6bfe6 100644 --- a/cube_bathymetry/docs/decisions/0007-mbes-backscatter-store-addendum-phase-b-transition.md +++ b/cube_bathymetry/docs/decisions/0007-mbes-backscatter-store-addendum-phase-b-transition.md @@ -23,6 +23,9 @@ now the node-output correction was a **Phase B no-op**: `corrected == raw` for every beam, with the per-beam pairs kept so the surfaced backscatter stays re-derivable once a real correction lands. +> **Superseded (cube#92/#93):** the correction was later moved to record time and +> the per-beam raw set is no longer retained — see the *Phase B.2 addendum* below. + An empirical study of the `m3_dryrun` Portsmouth Harbor survey (`/detections`, 2025 pings, ~445k beams) showed the earlier candidate model — a cos²θ Lambert correction — is wrong for this data: it explains only ~¼ of the observed @@ -83,6 +86,9 @@ sonar. The `rx_angles` sign/zero convention was verified against reflect "uncorrected by default". - The per-beam `{raw_intensity, beam_angle}` retention is unchanged, so the node value stays fully re-derivable when a richer correction lands. + **(Superseded by the Phase B.2 addendum, cube#92/#93: raw samples are no longer + retained — only an O(1) corrected Welford — so a richer correction now requires + re-importing the bag rather than in-place re-derivation.)** - Tests: 6 ARA cases at `extractNodeRecord()` (interpolation, nadir identity, beyond-max identity, NaN-angle identity, mode-None identity, port/starboard symmetry) plus curve-loader/mode-parse tests; the 3 pre-existing default-path @@ -120,7 +126,10 @@ Design decisions: verbatim and never recomputes α, so Python and C++ apply an **identical** TL by construction (a divergence would silently corrupt the correction). At 500 kHz / 24 °C / fresh water α ≈ 0.049 dB/m. -3. **Per-beam slant range R is a new sufficient statistic**, threaded exactly like +3. **Per-beam slant range R is a new sufficient statistic** (Superseded by the + Phase B.2 addendum, cube#92/#93: R is now consumed by `correctBeamIntensity` at + record time and folded into the Welford, not retained per beam), threaded + exactly like `beam_angle`: `Sounding::slant_range` → `DepthAndUncertainty::range` (the pack(1) raster struct grows to 20 bytes; layout-safe because every raster read strides by `sizeof(DepthAndUncertainty)`) → `BeamIntensitySample::range`, @@ -139,6 +148,46 @@ beam-pattern), which is the deferred tier-3. The TVG/absorption/frequency state ultimately belongs in `SonarInfo` (unh_marine_autonomy#240); α is computed locally meanwhile. +## Phase B.2 addendum — record-time correction + streaming Welford (cube_bathymetry#92/#93) + +cube_bathymetry#93 moves the per-beam correction from **extract** to **record** +time and replaces the per-beam `{raw_intensity, beam_angle, range}` retention with +a streaming **Welford** `(n, mean, M2)` of the *corrected* intensity, kept on the +winning depth hypothesis. cube_bathymetry#92 then makes tile eviction **lossless** +by spilling and restoring that triplet. + +This **reverses** two earlier decisions recorded above: + +1. **Correction timing.** Decision 1 located the empirical-ARA (and the tier-2 TL) + correction at `Node::extractNodeRecord()`. It now runs at record + (`Hypothesis::recordBeam` → `correctBeamIntensity`, the per-beam math moved + verbatim to `angular_response_curve.cpp`); extract only reads the triplet. The + math is unchanged, so the node-output mean + estimate variance are **identical** + to the old correct-at-extract path — only the timing moved. + +2. **Raw-sample retention / re-derivability.** The Context and Consequences above + state the per-beam `{raw_intensity, beam_angle}` set is retained so the surfaced + backscatter "stays fully re-derivable when a richer correction lands." That is + **no longer true**: only the O(1) corrected Welford is kept (the cube#93 + Massabesic OOM fix — per-cell intensity memory used to grow with total beam + count). A future richer correction (the deferred tier-3 GeoCoder / + cube_bathymetry#15) therefore requires **re-importing the source bag**, not + in-place re-derivation from stored raw samples. Per-beam slant range `R` + (Tier-2) is likewise folded into the correction at record and not retained. + +Why this is acceptable: the empirical ARA + tier-2 TL correction the store needs +*today* is applied before the Welford folds each beam, so the persisted value is +the final corrected backscatter; applying a *different* correction was always going +to be a re-processing operation, and bounding live/offline RAM (the cube#92/#93 +motivation) is the harder constraint. The Welford triplet is a perfect sufficient +statistic for the mean + estimate variance, so eviction spill/reload is +bit-identical to never-evicting. + +This addendum **supersedes** the "retention is unchanged / re-derivable" statements +in the Context, Consequences, and Tier-2 sections above, and the backscatter +"out of scope; re-derivable per cube#15" trade-off note in +`0001-tile-eviction-and-incremental-publish.md` (§ Negative / trade-offs). + ## Deferred (explicit follow-ups) - **Full radiometric GeoCoder (tier-3)** — insonified-area, beam-pattern, TVG diff --git a/cube_bathymetry/include/cube_bathymetry/angular_response_curve.h b/cube_bathymetry/include/cube_bathymetry/angular_response_curve.h index c8e2ba6..9f8e642 100644 --- a/cube_bathymetry/include/cube_bathymetry/angular_response_curve.h +++ b/cube_bathymetry/include/cube_bathymetry/angular_response_curve.h @@ -74,6 +74,27 @@ namespace cube /// provenance). Equivalent to `loadAngularResponseCurveWithHeader(path).points`. std::vector < std::pair < float, float >> loadAngularResponseCurve(const std::string & path); +/// @brief Apply the node-output backscatter correction to one raw beam +/// (ADR-0007 D3, cube_bathymetry#81/#87). +/// +/// Returns the corrected backscatter (dB) for a single beam, computed at RECORD +/// time (cube_bathymetry#93) so the CUBE hypothesis can stream a Welford of the +/// CORRECTED intensity instead of retaining every raw beam. This is the exact +/// per-beam math that used to run inside `Node::extractNodeRecord`, factored out +/// verbatim so correct-at-record reproduces correct-at-extract bit-for-bit: +/// - tier-2 TL add-back (cube#87): when @c backscatter_tl_removed and the range +/// is finite and > 0, add `40*log10(R) + 2*alpha*R` (alpha == +/// @c backscatter_absorption_db_per_m). A NaN / non-positive range skips it. +/// - angular-response residual (cube#81): when the mode is Empirical with a +/// non-empty curve and the beam angle is not NaN, subtract +/// `curve(|beam_angle| in degrees)` (linearly interpolated, edge-clamped). +/// Identity (returns @p raw_intensity) when the mode is None or the curve is +/// empty, and for the angle term when @p beam_angle is NaN. The caller gates a +/// NaN @p raw_intensity (a NaN beam is never a sample). Computed and returned in +/// `double` to match the estimator's accumulation precision. + double correctBeamIntensity( + float raw_intensity, float beam_angle, float range, const Parameters & parameters); + } // namespace cube #endif // CUBE_BATHYMETRY__ANGULAR_RESPONSE_CURVE_H_ diff --git a/cube_bathymetry/include/cube_bathymetry/geo_grid.h b/cube_bathymetry/include/cube_bathymetry/geo_grid.h index 411011a..a45fe0c 100644 --- a/cube_bathymetry/include/cube_bathymetry/geo_grid.h +++ b/cube_bathymetry/include/cube_bathymetry/geo_grid.h @@ -127,6 +127,29 @@ namespace cube /// beam-angle/GeoCoder correction is deferred to cube#81). std::vector < NodeRecord > nodeRecords() const; + /// @brief Corrected-intensity Welford of every cell's winning hypothesis, keyed + /// by `gggs::CellIndex` (cube_bathymetry#92/#93 lossless eviction spill). + /// + /// Iterates the SPARSE node map (not all 960x960 cells), flushes each node's + /// median pre-filter (`queueFlush`, like @ref nodeRecords) and emits the + /// @ref Node::chosenIntensityWelford for cells with intensity-bearing beams + /// (n > 0). The result is what must be persisted to a scratch spill so a tile + /// evicted from RAM can be reloaded with its backscatter accumulator intact + /// (the bathy tile stores only the depth summary; the intensity Welford lives + /// only here). Three numbers per cell -- O(1), the cube#93 memory fix. + std::map < gggs::CellIndex, IntensityWelford > nodeIntensityWelford() const; + + /// @brief Restore a corrected-intensity Welford onto @p cell's winning + /// hypothesis (cube_bathymetry#92/#93 lossless eviction reload). + /// + /// Forwards to @ref Node::setSettledIntensityWelford on the node at @p cell. A + /// no-op if no node exists there (the cell was not reseeded by the depth + /// reload). Must run AFTER @ref setSettledDepthAt (which lazy-creates the node) + /// and BEFORE the revisit's soundings are added, so those beams then continue + /// the Welford on the same reloaded hypothesis (bit-identical to never-evicting). + void setSettledIntensityWelfordAt( + const gggs::CellIndex & cell, const IntensityWelford & intensity); + private: gggs::GridIndex index_; diff --git a/cube_bathymetry/include/cube_bathymetry/geo_map_sheet.h b/cube_bathymetry/include/cube_bathymetry/geo_map_sheet.h index 4976b6e..669f8b9 100644 --- a/cube_bathymetry/include/cube_bathymetry/geo_map_sheet.h +++ b/cube_bathymetry/include/cube_bathymetry/geo_map_sheet.h @@ -50,6 +50,27 @@ namespace cube const std::vector < GeoSounding > &soundings, std::chrono::steady_clock::time_point time = std::chrono::steady_clock::now()); + /// @brief The grid indices @ref addSoundings would touch for @p soundings, + /// WITHOUT creating them (cube_bathymetry#92 reload-before-add). + /// + /// Computes the same expanded bounds @ref addSoundings uses (the sounding extent + /// grown by one cell so a near-seam sounding reaches its neighbour tile) and + /// returns every `gggs::GridIndex` in that window. The bounded-RAM importer calls + /// this BEFORE addSoundings to reload any evicted tile the batch is about to + /// touch, so the new soundings accrete onto the reloaded hypotheses (lossless + /// blend) instead of forming a fresh, partial tile. Empty for empty input. + std::vector < gggs::GridIndex > gridIndicesForSoundings( + const std::vector < GeoSounding > &soundings) const; + + /// @brief The grid index a single sounding's centre falls in (no creation). + /// + /// Maps the sounding position to its home tile via the sheet's grid level. The + /// bounded-RAM importer uses this to count the soundings lost when a tile that + /// failed to reload is dropped. A sounding's influence radius can spread it into + /// neighbour tiles too (see GeoGrid::insert), so this home-tile count is a + /// conservative floor on the affected soundings, not an exact cell tally. + gggs::GridIndex gridIndexForSounding(const GeoSounding & sounding) const; + /// @brief Configure the per-beam backscatter angular-response correction /// applied at node-output (ADR-0007 D3, cube_bathymetry#81). /// @@ -105,6 +126,18 @@ namespace cube /// startup prime. Does NOT mark the grid dirty (reproduces persisted data). void setSettledDepthAt(const gggs::CellIndex & cell, float depth, float uncertainty); + /// @brief Restore a corrected-intensity Welford onto @p cell's winning + /// hypothesis (cube_bathymetry#92/#93 lossless eviction reload). + /// + /// Forwards to `GeoGrid::setSettledIntensityWelfordAt` on the grid owning + /// @p cell (lazy-creates the grid; a no-op on the node if @ref setSettledDepthAt + /// did not seed it). Must run AFTER the depth reload and BEFORE the revisit's + /// soundings are added, so those beams continue the Welford on the same reloaded + /// hypothesis (bit-identical to never-evicting). Does NOT mark the grid dirty + /// (it reproduces already-persisted data). + void setSettledIntensityWelfordAt( + const gggs::CellIndex & cell, const IntensityWelford & intensity); + /// @brief Grid indices touched (returning true from insert) since the last /// clearDirtyGrids(). Returned by value -- safe to iterate while saving. std::set < gggs::GridIndex > dirtyGrids() const; diff --git a/cube_bathymetry/include/cube_bathymetry/hypothesis.h b/cube_bathymetry/include/cube_bathymetry/hypothesis.h index 4468aa3..dbf233b 100644 --- a/cube_bathymetry/include/cube_bathymetry/hypothesis.h +++ b/cube_bathymetry/include/cube_bathymetry/hypothesis.h @@ -32,32 +32,29 @@ namespace cube { -/// Per-beam backscatter sufficient statistics retained on a hypothesis -/// (ADR-0007 D3). Each contributing beam contributes its RAW (uncorrected) -/// intensity and the per-beam receive/steering (beam/incidence) angle, so the -/// radiometric (empirical-ARA / GeoCoder) correction can be applied later -- at -/// node-output, once depth and local slope have settled -- rather than baking -/// an angle-corrupted, non-re-correctable average into the hypothesis. - struct BeamIntensitySample +/// @brief Streaming Welford accumulator of the CORRECTED per-beam backscatter on +/// a hypothesis (ADR-0007 D3/D4, cube_bathymetry#93). +/// +/// Replaces the unbounded `std::vector` raw-sample retention: +/// the angular-response + tier-2 TL correction is now applied at RECORD time +/// (`correctBeamIntensity`), and only the running `(n, mean, M2)` of the corrected +/// dB value is kept -- O(1) per cell regardless of beam count, the fix for the +/// Massabesic OOM (per-cell intensity memory used to grow with total beams). +/// +/// `mean` and `m2` (sum of squared deviations from the mean) are `double` for +/// numerical stability over the ~10^9 beams a multi-day survey accumulates. +/// Node-output reads `intensity = mean` and the estimate variance +/// `intensity_var = (m2/(n-1)) / n` (variance of the mean, ADR-0007 D4) -- the +/// SAME values the retain-and-correct-at-extract path produced. + struct IntensityWelford { - /// Raw, uncorrected per-beam intensity (e.g. M3 reflectivity in dB). - float raw_intensity; - - /// Per-beam receive/steering (beam/incidence) angle from nadir, in radians - /// (NOT a true grazing angle of 90 deg - theta; nadir = 0). May be NaN when - /// the source did not report an angle for the beam; the output correction - /// treats a NaN angle as "no angle correction available" and emits that beam - /// uncorrected. - float beam_angle; - - /// Per-beam slant range R from the sonar head to the touchdown, in meters, - /// for the tier-2 backscatter 2-way transmission-loss correction - /// (cube_bathymetry#87): `corrected = raw + (40*log10(R) + 2*alpha*R) - - /// residualCurve(|angle|)` (the TL is ADDED BACK to compensate the loss). NaN - /// (or non-positive) when not reported; the TL - /// term is then skipped (identity) so the beam is corrected by the residual - /// angular-response curve alone (tier-1 behavior). - float range = std::numeric_limits < float > ::quiet_NaN(); + /// Number of corrected intensity samples folded in (intensity-bearing beams). + uint32_t n = 0; + /// Running mean of the corrected intensity (dB). + double mean = 0.0; + /// Running sum of squared deviations from the mean (Welford M2); m2/(n-1) is + /// the unbiased sample variance. + double m2 = 0.0; }; /// Depth hypothesis structure used to maintain a current track on the depth @@ -117,20 +114,23 @@ namespace cube /// the hypothesis represents (i.e., an intervention is required). bool update(float depth, float variance, const Parameters & parameters); - /// Record one beam's backscatter sufficient statistics on this hypothesis - /// (ADR-0007 D3). Called by Node::update() only for a beam whose depth was - /// accepted by (or used to seed) THIS hypothesis, so the backscatter - /// association mirrors the depth association exactly -- depth-geometry - /// outliers the W&H monitor rejects never enter this set. + /// Record one beam's backscatter on this hypothesis (ADR-0007 D3/D4, + /// cube_bathymetry#93). Called by Node::update() only for a beam whose depth was + /// accepted by (or used to seed) THIS hypothesis, so the backscatter association + /// mirrors the depth association exactly -- depth-geometry outliers the W&H + /// monitor rejects never enter the accumulator. /// - /// A NaN raw_intensity is skipped (a source that omits intensities must never - /// inject a phantom sample); a NaN beam_angle is retained (the beam is - /// still a valid intensity sample, merely uncorrectable for angle). A NaN / - /// non-positive range is also retained (the beam is a valid intensity sample, - /// merely uncorrectable for the tier-2 TL term, cube_bathymetry#87). + /// Applies the angular-response + tier-2 TL correction (`correctBeamIntensity`, + /// using @p parameters) at RECORD time, then folds the corrected value into the + /// streaming Welford `intensity` (no raw-sample retention -- O(1) per cell). + /// + /// A NaN @p raw_intensity is skipped (a source that omits intensities must never + /// inject a phantom sample); a NaN @p beam_angle is still counted (a valid + /// intensity sample, merely uncorrected for angle); a NaN / non-positive + /// @p range still counted (the tier-2 TL term is skipped for it, cube#87). void recordBeam( - float raw_intensity, float beam_angle, - float range = std::numeric_limits < float > ::quiet_NaN()); + float raw_intensity, float beam_angle, float range, + const Parameters & parameters); /// Current depth mean estimate double current_estimate; @@ -171,20 +171,14 @@ namespace cube /// This tracks the maximum of the two estimates. float maximum_of_input_and_predicted_variance = 0.0; - /// Per-beam backscatter sufficient statistics for the beams associated with - /// this hypothesis (ADR-0007 D3). Bounded by hypothesis membership -- one - /// entry per contributing beam with a non-NaN intensity, a few bytes each. - /// The radiometric correction (D3 GeoCoder) is applied per element at - /// node-output (Node::extractNodeRecord); these raw pairs are retained so the - /// node value stays re-derivable when slope (cube_bathymetry#15) lands. - /// - /// Memory budget: sizeof(BeamIntensitySample) == 8 bytes (two floats). Growth - /// class is identical to number_of_samples (one entry per accepted beam, for - /// the survey lifetime of the hypothesis). Worst-case is ~8 bytes/beam/node. - /// Once cube_bathymetry#15's correction model settles this can be reduced to - /// pure sufficient statistics (mean, M2, count) if the per-beam retention is - /// no longer needed for re-derivation. - std::vector < BeamIntensitySample > intensity_samples; + /// Streaming Welford of the CORRECTED backscatter for the beams associated with + /// this hypothesis (ADR-0007 D3/D4, cube_bathymetry#93). O(1) -- 16 bytes, + /// independent of beam count (the cube#93 OOM fix). The correction is applied at + /// record (`recordBeam` -> `correctBeamIntensity`), so this is already the + /// node-output value; `Node::extractNodeRecord` reads it without re-correcting. + /// The tile-eviction spill/reload persists this triplet (a perfect sufficient + /// statistic, so reload + continue is bit-identical to never-evicting, #92). + IntensityWelford intensity; }; } // namespace cube diff --git a/cube_bathymetry/include/cube_bathymetry/node.h b/cube_bathymetry/include/cube_bathymetry/node.h index deff492..4ca8692 100644 --- a/cube_bathymetry/include/cube_bathymetry/node.h +++ b/cube_bathymetry/include/cube_bathymetry/node.h @@ -54,9 +54,10 @@ namespace cube float depth_var = std::nan(""); /// Co-estimated backscatter intensity: the mean of the per-beam corrected - /// intensities on the winning hypothesis. NaN when no intensity-bearing beams. - /// (Phase B: correction is currently the identity -- emitted UNCORRECTED - /// pending cube_bathymetry#15; see extractNodeRecord().) + /// intensities (streaming Welford) on the winning hypothesis. NaN when no + /// intensity-bearing beams. The correction (empirical ARA + tier-2 TL) is + /// applied at record time and is the identity only when no curve is configured; + /// see extractNodeRecord(). float intensity = std::nan(""); /// Intensity ESTIMATE variance (variance of the mean, shrinks with n_samples; @@ -196,14 +197,39 @@ namespace cube /// co-estimated backscatter intensity/uncertainty/sample-count of the winning /// hypothesis. The depth fields match extractDepthAndUncertainty() exactly, /// including the nominated_hypothesis_ priority path, so the two outputs never - /// disagree for the same node state. Intensity is computed from the winning - /// hypothesis's per-beam {raw, angle} set: each beam is angle-corrected (D3), - /// then the corrected values are combined into a mean + ESTIMATE variance - /// (D4). The angle correction is currently the identity (no-op) pending - /// cube_bathymetry#15 -- intensity is emitted UNCORRECTED, but the per-beam - /// raw set is retained so it is re-derivable when #15 provides slope. + /// disagree for the same node state. Intensity is the winning hypothesis's + /// streaming Welford of the CORRECTED per-beam backscatter: each beam is + /// angle-/TL-corrected at RECORD time (recordBeam -> correctBeamIntensity, + /// D3/cube_bathymetry#93), so extract just reads the triplet -- mean as the + /// intensity and M2/(n-1)/n as the ESTIMATE variance (D4), with no re-correction. + /// This O(1) (n, mean, M2) replaces the old per-beam {raw, angle} retention (the + /// cube_bathymetry#93 OOM fix): raw samples are no longer kept, so a future + /// richer correction (cube_bathymetry#15) requires re-importing the bag rather + /// than re-deriving in place. NodeRecord extractNodeRecord(const Parameters & parameters); + /// @brief The corrected-intensity Welford of the WINNING hypothesis + /// (cube_bathymetry#92/#93 lossless eviction spill). + /// + /// Returns the `(n, mean, M2)` of exactly the hypothesis @ref extractNodeRecord + /// would pick (the `nominated_hypothesis_` priority path, else + /// @ref chooseHypothesis), so spilling and restoring the triplet reproduces the + /// same node-output intensity. A default (n == 0) triplet when there is no valid + /// hypothesis or it carries no intensity-bearing beams. Does NOT flush the queue + /// -- the caller flushes first (GeoGrid does, via queueFlush). + IntensityWelford chosenIntensityWelford(); + + /// @brief Restore a corrected-intensity Welford onto the WINNING hypothesis + /// (cube_bathymetry#92/#93 lossless eviction reload). + /// + /// Sets @p intensity as the Welford of the chosen hypothesis (same selection as + /// @ref chosenIntensityWelford / @ref extractNodeRecord). Used by the tile- + /// eviction reload, which runs BEFORE the revisit's soundings are added so those + /// beams then CONTINUE the Welford on the same hypothesis. Because `(n, mean, + /// M2)` is a perfect sufficient statistic, restore-then-continue is bit-identical + /// to never-evicting. A no-op when the node has no hypothesis. + void setSettledIntensityWelford(const IntensityWelford & intensity); + /* Routine: cube_node_choose_hypothesis * Purpose: Choose the current best hypothesis for the node in question * Inputs: *list Pointer to the list of hypotheses diff --git a/cube_bathymetry/include/cube_bathymetry/store_import.h b/cube_bathymetry/include/cube_bathymetry/store_import.h index 5070e3c..f093e89 100644 --- a/cube_bathymetry/include/cube_bathymetry/store_import.h +++ b/cube_bathymetry/include/cube_bathymetry/store_import.h @@ -23,16 +23,23 @@ #ifndef CUBE_BATHYMETRY__STORE_IMPORT_H_ #define CUBE_BATHYMETRY__STORE_IMPORT_H_ +#include #include #include +#include +#include +#include #include "cube_bathymetry/geo_grid.h" #include "cube_bathymetry/geo_map_sheet.h" +#include "cube_bathymetry/geo_sounding.h" #include "marine_autonomy/gggs.h" #include "marine_bathymetry_store/bathymetry_store.hpp" #include "marine_bathymetry_store/bathymetry_tile.hpp" #include "marine_bathymetry_store/bathy_cell.hpp" +#include "marine_bathymetry_store/registry.hpp" #include "marine_mbes_backscatter_store/mbes_cell.hpp" +#include "marine_mbes_backscatter_store/registry.hpp" namespace cube { @@ -182,6 +189,133 @@ namespace cube GeoMapSheet & map_sheet, bool seed_settled = true); +/// @brief Configuration for @ref ImportAccumulator (cube_bathymetry#92). + struct ImportAccumulatorConfig + { + /// Output bathymetry-store directory (the importer's `-o`). Evicted and + /// final-resident bathy tiles are written here under the @ref bathy_layer + /// subdirectory; a revisited tile is reloaded from here. Empty = no persistence + /// (eviction is then disabled to stay lossless — there is nowhere to drop to). + std::string store_dir; + /// Bathy store layer the import writes (`Processed` default, `Draft` opt-in, #85). + marine_bathymetry_store::SourceLayer bathy_layer = + marine_bathymetry_store::SourceLayer::Processed; + /// Registry source index stamped into every persisted bathy cell. + uint16_t source_index = 0; + /// Deterministic per-cell acquisition timestamp (ns since the Unix epoch). + int64_t timestamp_ns = 0; + /// Survey nominal cell size (m); fixes the GGGS level of the scratch stores used + /// to reload an evicted tile and to merge backscatter. Must equal the + /// GeoMapSheet's `nominalCellSizeMeters()` so the levels line up. + float cell_size_m = 1.0f; + /// Optional MBES backscatter store directory (the importer's `--bs-store`). + /// Empty disables backscatter co-persistence. + std::string bs_store_dir; + /// Registry source index stamped into every persisted backscatter cell. + uint16_t bs_source_index = 0; + /// Maximum resident GeoGrid tiles before persist-then-drop eviction runs. + /// 0 = unbounded (never evict — the pre-#92 whole-survey-in-RAM behavior). + std::size_t max_resident_tiles = 0; + }; + +/// @brief Bounded-RAM offline import accumulator (cube_bathymetry#92). +/// +/// Ports the live node's cube#70 persist-then-drop eviction + lossless +/// reload-on-revisit (ADR-0001) to the offline `import_bag` path, so resident RAM +/// is bounded by tile COUNT rather than surveyed AREA (a multi-day survey used to +/// OOM the importer because `grids_` grew with coverage). The host feeds one +/// batch of soundings per call (one ping in `import_bag`, a synthetic region in +/// tests); the accumulator drives the underlying @ref GeoMapSheet and persists to +/// the configured stores. +/// +/// **Persist-then-drop (lossless):** when the resident tile count exceeds the +/// budget, each cold tile is written to disk — bathy via @ref geoGridToTile + +/// `saveTile`, backscatter SUMMARY via `saveTile`, AND its per-cell intensity +/// Welford `(n, mean, M2)` spilled to a temporary scratch — BEFORE it is dropped +/// from RAM. A tile whose persist throws is left resident (never dropped), so a +/// disk failure costs transient RAM, never data. +/// +/// **Reload-before-add (lossless blend):** the bathy tile stores only the depth +/// SUMMARY and the backscatter SUMMARY is the corrected mean — neither retains the +/// running Welford CUBE needs to keep blending. So eviction also spills each cell's +/// intensity Welford to a scratch file (cube#93: a 3-number sufficient statistic of +/// the CORRECTED intensity, since correction now runs at record), and the reload +/// runs BEFORE a batch's soundings are added: it computes the grids the batch will +/// touch (@ref GeoMapSheet::gridIndicesForSoundings), `loadWindow` + +/// @ref primeFromTile restores each cell's settled depth as one CUBE hypothesis, +/// then the spilled Welford is restored onto that hypothesis. The batch's new beams +/// then continue the Welford on the SAME reloaded hypothesis, so the node-output +/// intensity is the FULL pre+post-eviction blend — bit-for-bit equal to a +/// never-evicted build for a consistent re-survey (a Welford triplet is a perfect +/// sufficient statistic). **Backscatter is lossless under eviction.** +/// +/// **Remaining approximation (bathy uncertainty only):** the reload reconstructs a +/// SINGLE depth hypothesis (a Bayesian prior from the stored depth + variance, +/// ADR-0001), seeded with one sample rather than the original count. The depth +/// VALUE is faithful (the prior is refined by the revisit), but the re-derived +/// depth UNCERTAINTY drifts slightly from the unbounded build for a tile evicted +/// mid-disambiguation — a hypothesis-state collapse, not data loss. The backscatter +/// estimate does not depend on the depth sample count, so it is unaffected. +/// +/// **Transient disk cost:** the spill holds a fixed 24-byte Welford record per +/// surveyed cell for every currently-evicted tile (NOT per beam — cube#93 made +/// per-cell intensity state O(1)), in a scratch dir deleted in @ref finalize (and +/// by the destructor on an exception). Scratch-only and local/fast. + class ImportAccumulator + { +public: + /// @param sheet The map sheet to accumulate into (lifetime must outlast this). + /// @param config Persistence + budget configuration. + ImportAccumulator(GeoMapSheet & sheet, ImportAccumulatorConfig config); + + /// @brief Clean up the scratch spill directory (RAII safety net for finalize). + ~ImportAccumulator(); + + /// @brief Reload-before-add any evicted tile this batch touches, accumulate the + /// batch, then evict cold tiles back to budget. + void addBatch( + const std::vector < GeoSounding > &soundings, + std::chrono::steady_clock::time_point time = + std::chrono::steady_clock::now()); + + /// @brief Persist every still-resident tile (bathy + backscatter), write both + /// registries, and delete the scratch spill. Evicted tiles are already + /// durable. Call once at end-of-stream. + void finalize( + const marine_bathymetry_store::SourceRegistry & bathy_registry, + const marine_mbes_backscatter_store::SourceRegistry * bs_registry = nullptr); + + /// @brief Tiles currently resident in RAM. + std::size_t residentTileCount() const; + /// @brief Indices evicted to disk and not yet reloaded. + const std::set < gggs::GridIndex > & evictedIndices() const { + return evicted_; + } + /// @brief Cumulative bathy tile writes (eviction + finalize). + std::size_t bathyTilesPersisted() const {return bathy_persisted_;} + /// @brief Cumulative backscatter tile writes (eviction + finalize). + std::size_t backscatterTilesPersisted() const {return bs_persisted_;} + /// @brief The scratch spill directory (empty until the first eviction). Exposed + /// for tests that assert it is cleaned up after @ref finalize. + const std::string & scratchDir() const {return scratch_dir_;} + +private: + void evictColdTiles(); + void persistBathyTile(const gggs::GridIndex & index); + void persistBackscatterTile(const gggs::GridIndex & index); + void spillIntensitySamples(const gggs::GridIndex & index); + void restoreSpilledSamples(const gggs::GridIndex & index); + bool reloadEvictedTile(const gggs::GridIndex & index); + void cleanupScratch(); + + GeoMapSheet & sheet_; + ImportAccumulatorConfig cfg_; + std::set < gggs::GridIndex > evicted_; + std::string scratch_dir_; // lazily created on first eviction; "" = none + std::size_t bathy_persisted_ = 0; + std::size_t bs_persisted_ = 0; + }; + } // namespace cube #endif // CUBE_BATHYMETRY__STORE_IMPORT_H_ diff --git a/cube_bathymetry/src/angular_response_curve.cpp b/cube_bathymetry/src/angular_response_curve.cpp index 346fb79..202e8fa 100644 --- a/cube_bathymetry/src/angular_response_curve.cpp +++ b/cube_bathymetry/src/angular_response_curve.cpp @@ -24,6 +24,8 @@ #include #include +#include +#include #include #include #include @@ -33,6 +35,84 @@ namespace cube { +namespace +{ +/// Empirical angular-response lookup: the curve's db_relative_to_nadir at +/// @p abs_angle_deg, linearly interpolated between adjacent bin centres. The +/// curve is ascending {abs_angle_deg, db_relative_to_nadir} pairs (loaded by +/// loadAngularResponseCurve). Returns 0 (identity) when the curve is empty or +/// @p abs_angle_deg lies beyond the curve's max angle; clamps to the first bin +/// below the curve's min angle (nadir bin ~0 -> ~identity near nadir). +/// +/// Moved verbatim from node.cpp (cube_bathymetry#93): the correction now runs at +/// record (recordBeam) instead of extract, so this helper lives next to +/// correctBeamIntensity. Unchanged logic -> identical corrected values. +double curveRelativeDb( + const std::vector> & curve, double abs_angle_deg) +{ + if (curve.empty()) { + return 0.0; + } + // Below the first bin centre: clamp to it (the nadir bin's value, ~0). + if (abs_angle_deg <= curve.front().first) { + return curve.front().second; + } + // Beyond the last bin centre: clamp to the outermost bin's correction rather + // than jump to identity. The empirical curve is bounded (unlike a cos/log model + // that diverges near grazing), so continuing the edge value keeps the correction + // continuous and avoids a swath-edge discontinuity / bright ring (#81 review). + if (abs_angle_deg > curve.back().first) { + return curve.back().second; + } + // Find the bracketing pair [lo, hi] and linearly interpolate. + for (std::size_t i = 1; i < curve.size(); ++i) { + if (abs_angle_deg <= curve[i].first) { + const double a0 = curve[i - 1].first; + const double d0 = curve[i - 1].second; + const double a1 = curve[i].first; + const double d1 = curve[i].second; + const double span = a1 - a0; + if (span <= 0.0) { + return d1; // duplicate angle: take the upper bin's value + } + const double t = (abs_angle_deg - a0) / span; + return d0 + t * (d1 - d0); + } + } + return 0.0; // unreachable (guarded by the back() check above) +} +} // namespace + +double correctBeamIntensity( + float raw_intensity, float beam_angle, float range, const Parameters & parameters) +{ + // Exact per-beam math formerly in Node::extractNodeRecord (cube#80/#81/#87), + // moved to record time (cube#93). Apply order is TL add-back THEN angular + // residual, accumulated in double -- identical to the old extract. + const bool apply_ara = + parameters.backscatter_angle_correction == + BackscatterAngleCorrection::Empirical && + !parameters.angular_response_curve.empty(); + const bool apply_tl = apply_ara && parameters.backscatter_tl_removed; + const double alpha = parameters.backscatter_absorption_db_per_m; + + double corrected = raw_intensity; + if (apply_tl) { + const double range_d = static_cast(range); + // Skip the TL term for a missing / non-positive range (log10 undefined); the + // beam is still corrected by the residual angular-response curve below. + if (std::isfinite(range_d) && range_d > 0.0) { + corrected += 40.0 * std::log10(range_d) + 2.0 * alpha * range_d; + } + } + if (apply_ara && !std::isnan(beam_angle)) { + const double abs_angle_deg = + std::abs(static_cast(beam_angle)) * 180.0 / M_PI; + corrected -= curveRelativeDb(parameters.angular_response_curve, abs_angle_deg); + } + return corrected; +} + bool parseBackscatterAngleCorrection( const std::string & text, BackscatterAngleCorrection & out) { diff --git a/cube_bathymetry/src/geo_grid.cpp b/cube_bathymetry/src/geo_grid.cpp index 35b3c68..254c2a3 100644 --- a/cube_bathymetry/src/geo_grid.cpp +++ b/cube_bathymetry/src/geo_grid.cpp @@ -185,4 +185,34 @@ std::vector GeoGrid::nodeRecords() const return ret; } +std::map +GeoGrid::nodeIntensityWelford() const +{ + std::map ret; + // Iterate the SPARSE node map (only created cells), not a CellAreaIterator over + // all 960x960 cells: a tile has far fewer surveyed nodes than cells, so this is + // both faster and avoids visiting ~900k empty cells. + for (const auto & entry : nodes_) { + if(!entry.second) { + continue; + } + entry.second->queueFlush(parameters_); + const IntensityWelford w = entry.second->chosenIntensityWelford(); + if(w.n > 0) { + ret.emplace(entry.first, w); + } + } + return ret; +} + +void GeoGrid::setSettledIntensityWelfordAt( + const gggs::CellIndex & cell, const IntensityWelford & intensity) +{ + auto it = nodes_.find(cell); + if(it == nodes_.end() || !it->second) { + return; // no node here (depth reload did not seed this cell) -- nothing to do + } + it->second->setSettledIntensityWelford(intensity); +} + } // namespace cube diff --git a/cube_bathymetry/src/geo_map_sheet.cpp b/cube_bathymetry/src/geo_map_sheet.cpp index a6554ef..f1dd664 100644 --- a/cube_bathymetry/src/geo_map_sheet.cpp +++ b/cube_bathymetry/src/geo_map_sheet.cpp @@ -30,6 +30,32 @@ namespace cube { +namespace +{ +// Expanded geographic bounds covering a sounding batch: the sounding extent grown +// by one cell on every side so a sounding near a tile seam still reaches its +// neighbour tile. Shared by addSoundings (which then creates the grids) and +// gridIndicesForSoundings (which only enumerates them) so the two never drift. +gz4d::BoundsDegrees boundsForSoundings( + const std::vector & soundings, const gggs::Level & grid_level) +{ + gz4d::BoundsDegrees bounds; + for (const auto & s : soundings) { + bounds.expand(s); + } + // quick hack to make sure to go a bit beyond the outer soundings + const auto angular_span = grid_level.cellAngularSpan(); + auto min = bounds.minimum(); + min = gz4d::PositionDegrees(min.latitude - angular_span, min.longitude - angular_span); + bounds.expand(min); + + auto max = bounds.maximum(); + max = gz4d::PositionDegrees(max.latitude + angular_span, max.longitude + angular_span); + bounds.expand(max); + return bounds; +} +} // namespace + GeoMapSheet::GeoMapSheet(float cell_size, std::string iho_order) :parameters_(CellSizes(cell_size), iho_order), grid_level_(gggs::Level::fromCellSize(cell_size)) { @@ -56,22 +82,7 @@ void GeoMapSheet::addSoundings( return; } - gz4d::BoundsDegrees bounds; - for (const auto & s : soundings) { - bounds.expand(s); - } - - // quick hack to make sure to go a bit beyond the outer soundings - auto angular_span = grid_level_.cellAngularSpan(); - auto min = bounds.minimum(); - min = gz4d::PositionDegrees(min.latitude - angular_span, min.longitude - angular_span); - bounds.expand(min); - - auto max = bounds.maximum(); - max = gz4d::PositionDegrees(max.latitude + angular_span, max.longitude + angular_span); - bounds.expand(max); - - auto grids = getOrCreateGridsIn(bounds); + auto grids = getOrCreateGridsIn(boundsForSoundings(soundings, grid_level_)); for (auto g : grids) { if(g->insert(soundings)) { last_update_time_ = time; @@ -84,6 +95,38 @@ void GeoMapSheet::addSoundings( } } +std::vector GeoMapSheet::gridIndicesForSoundings( + const std::vector & soundings) const +{ + std::vector ret; + if(soundings.empty()) { + return ret; + } + const gz4d::BoundsDegrees bounds = boundsForSoundings(soundings, grid_level_); + gggs::GridAreaIterator i( + grid_level_.gridIndex(bounds.minimum().latitude, bounds.minimum().longitude), + grid_level_.gridIndex(bounds.maximum().latitude, bounds.maximum().longitude)); + while(i.valid()) { + ret.push_back(*i); + i.next(); + } + return ret; +} + +gggs::GridIndex GeoMapSheet::gridIndexForSounding(const GeoSounding & sounding) const +{ + return grid_level_.gridIndex(sounding.latitude, sounding.longitude); +} + +void GeoMapSheet::setSettledIntensityWelfordAt( + const gggs::CellIndex & cell, const IntensityWelford & intensity) +{ + // Lazy-create the grid (mirrors setSettledDepthAt) so the call is safe even if + // the grid is absent; GeoGrid::setSettledIntensityWelfordAt is a no-op when the + // node was not seeded. Does NOT mark dirty (reproduces persisted data). + getOrCreateGrid(cell.grid())->setSettledIntensityWelfordAt(cell, intensity); +} + std::vector> GeoMapSheet::getOrCreateGridsIn( const gz4d::BoundsDegrees & bounds) { diff --git a/cube_bathymetry/src/hypothesis.cpp b/cube_bathymetry/src/hypothesis.cpp index 722935a..7ace651 100644 --- a/cube_bathymetry/src/hypothesis.cpp +++ b/cube_bathymetry/src/hypothesis.cpp @@ -24,6 +24,8 @@ #include +#include "cube_bathymetry/angular_response_curve.h" + namespace cube { @@ -115,18 +117,31 @@ bool Hypothesis::update(float depth, float variance, const Parameters & paramete return true; } -void Hypothesis::recordBeam(float raw_intensity, float beam_angle, float range) +void Hypothesis::recordBeam( + float raw_intensity, float beam_angle, float range, const Parameters & parameters) { // Skip beams with no reported intensity -- a NaN must never be mistaken for a // real backscatter sample (ADR-0007: missing intensity != zero backscatter). // A NaN beam_angle is allowed through: the beam is still a valid intensity - // measurement, it merely cannot be angle-corrected at output. A NaN / non- - // positive range is likewise allowed through (the tier-2 TL term is skipped - // for that beam, cube_bathymetry#87). + // measurement, it merely cannot be angle-corrected. A NaN / non-positive range + // is likewise allowed through (the tier-2 TL term is skipped, cube#87). if (std::isnan(raw_intensity)) { return; } - intensity_samples.push_back(BeamIntensitySample{raw_intensity, beam_angle, range}); + + // Correct at RECORD (cube#93): apply the angular-response + tier-2 TL correction + // now, then stream the corrected value into the Welford -- O(1) per cell instead + // of retaining every raw beam. correctBeamIntensity is the exact per-beam math + // that used to run at extract, so the node-output mean/variance are unchanged. + const double corrected = + correctBeamIntensity(raw_intensity, beam_angle, range, parameters); + + // Welford online update (n, mean, M2) in double precision. + ++intensity.n; + const double delta = corrected - intensity.mean; + intensity.mean += delta / intensity.n; + const double delta2 = corrected - intensity.mean; + intensity.m2 += delta * delta2; } diff --git a/cube_bathymetry/src/import_bag_main.cpp b/cube_bathymetry/src/import_bag_main.cpp index aaf6c40..74f082e 100644 --- a/cube_bathymetry/src/import_bag_main.cpp +++ b/cube_bathymetry/src/import_bag_main.cpp @@ -98,6 +98,16 @@ "--backscatter-correction empirical; empty -> correction is a no-op. A tier-2 " "curve (header '# tl_removed: true' + '# absorption_db_per_m: ') also makes " "the estimator remove per-beam 2-way TL 40*log10(R)+2*alpha*R (cube#87)\n"; + std::cout << " --max-resident-tiles : bound resident-tile RAM (cube#92). When " + "the in-memory GGGS tile count exceeds N, the coldest tiles are persisted to " + "the -o store (and their backscatter to --bs-store), their per-cell intensity " + "Welford spilled to a temp scratch, and dropped from RAM. A revisit reloads the " + "tile (settled depth + restored Welford) BEFORE the new soundings, so " + "backscatter blends LOSSLESSLY with the pre-eviction beams; only the bathy " + "depth UNCERTAINTY is re-derived (the depth value stays faithful). Default 256 " + "(generous; disk is local/fast offline). 0 = unbounded (whole survey in RAM). " + "NOTE: per-cell intensity memory is O(1) regardless of beam count (cube#93), so " + "a small heavily-oversampled survey no longer OOMs even without eviction\n"; std::cout << " -l : Stop after this many pings (debugging)\n"; std::cout << " --source-id : Registry source id recorded for every cell " "(default cube-replay)\n"; @@ -291,6 +301,11 @@ int main(int argc, char * argv[]) double resolution = 1.0; std::string iho_order = "order1a"; int ping_count_limit = 0; + // Bounded-RAM eviction budget (cube#92). Default 256 resident tiles: generous + // for offline (each ~960x960-cell GeoGrid + CUBE state is the heavy object), so + // small/medium surveys never evict (identical output to the old path) while a + // large multi-day survey stays bounded. 0 = unbounded (old behavior). + std::size_t max_resident_tiles = 256; // Backscatter angular-response correction (cube#81). Default none = identity. std::string backscatter_correction_str = "none"; std::string backscatter_curve_file; @@ -342,6 +357,13 @@ int main(int argc, char * argv[]) backscatter_correction_str = next_value("--backscatter-correction"); } else if (*arg == "--backscatter-curve") { backscatter_curve_file = next_value("--backscatter-curve"); + } else if (*arg == "--max-resident-tiles") { + const int64_t v = std::stoll(next_value("--max-resident-tiles")); + if (v < 0) { + std::cerr << "error: --max-resident-tiles must be >= 0 (0 = unbounded)\n"; + usage(); + } + max_resident_tiles = static_cast(v); } else if (*arg == "-l") { ping_count_limit = std::stoi(next_value("-l")); } else if (*arg == "--source-id") { @@ -572,6 +594,51 @@ int main(int argc, char * argv[]) << " (blunder rejection active, #89)." << std::endl; } + // Bathy store provenance + bounded-RAM accumulator (cube#92). The registry/ + // source index are created up front (not at the end) so evicted tiles persisted + // mid-pass already carry the right source index; registry.json is written once + // at finalize(). The accumulator owns the persist-then-drop eviction + lossless + // reload-on-revisit that bounds resident RAM by tile COUNT, not surveyed AREA. + marine_bathymetry_store::SourceRegistry registry; + const uint16_t source_index = registry.registerSource(source_record); + + // Backscatter store provenance (cube#80), registered up front for the same + // reason. Only used when --bs-store is set. + marine_mbes_backscatter_store::SourceRegistry bs_registry; + uint16_t bs_source_index = 0; + if (!bs_store_dir.empty()) { + marine_mbes_backscatter_store::SourceRecord bs_source_record; + bs_source_record.source_id = source_record.source_id; + bs_source_record.platform = source_record.platform; + bs_source_record.sensor = source_record.sensor; + bs_source_record.sensor_class = "mbes-backscatter"; + bs_source_record.campaign = source_record.campaign; + bs_source_index = bs_registry.registerSource(bs_source_record); + } + + cube::ImportAccumulatorConfig accumulator_config; + accumulator_config.store_dir = store_dir; + accumulator_config.bathy_layer = bathy_layer; + accumulator_config.source_index = source_index; + accumulator_config.timestamp_ns = cell_timestamp_ns; + // Match the store level to the sheet's actual (GGGS-snapped) cell size so the + // reload/merge scratch stores tile identically. + accumulator_config.cell_size_m = + static_cast(geo_map_sheet.nominalCellSizeMeters()); + accumulator_config.bs_store_dir = bs_store_dir; + accumulator_config.bs_source_index = bs_source_index; + accumulator_config.max_resident_tiles = max_resident_tiles; + cube::ImportAccumulator accumulator(geo_map_sheet, accumulator_config); + + if (max_resident_tiles > 0) { + std::cout << "Bounded resident tiles: " << max_resident_tiles + << " (cold tiles persist to the -o store and drop from RAM; " + "revisits reload losslessly, cube#92)." << std::endl; + } else { + std::cout << "Unbounded resident tiles (--max-resident-tiles 0): the whole " + "survey stays in RAM." << std::endl; + } + std::cout << "reading messages..." << std::endl; int ping_count = 0; @@ -687,7 +754,11 @@ int main(int argc, char * argv[]) gs.sounding.slant_range = s.slant_range; soundings.push_back(gs); } - geo_map_sheet.addSoundings(soundings); + // Accumulate through the bounded-RAM accumulator (cube#92): adds the + // batch, reloads any evicted tile this ping revisits, then evicts cold + // tiles back to the budget. With --max-resident-tiles 0 this is a plain + // addSoundings (no eviction). + accumulator.addBatch(soundings); ping_count++; } catch (const tf2::TransformException & e) { // A ping with no earth transform in the (bounded) buffer at its stamp -- @@ -818,85 +889,41 @@ int main(int argc, char * argv[]) std::cout << "Building store tiles..." << std::endl; - // The GeoMapSheet picks a GGGS level from the requested cell size; build the - // store at the matching default level. importTiles is multi-level, so the - // GridIndex on each tile carries the authoritative level regardless. - marine_bathymetry_store::BathymetryStore store = - marine_bathymetry_store::BathymetryStore::fromCellSize( - static_cast(geo_map_sheet.nominalCellSizeMeters())); - - marine_bathymetry_store::SourceRegistry registry; - const uint16_t source_index = registry.registerSource(source_record); - - auto tiles = cube::mapSheetToTiles(geo_map_sheet, cell_timestamp_ns, source_index); - std::cout << "Tiles with data: " << tiles.size() << " (build: " << phase_secs() << "s)" - << std::endl; - - if (tiles.empty()) { + // Persist the still-resident tiles and write the registries. Tiles evicted + // during the pass were already written to disk (bathy in the -o store, their + // backscatter merged into --bs-store newest-finite-wins); finalize() writes + // whatever is still in RAM, so the on-disk store is the union of evicted + + // resident — identical to an unbounded build (cube#92). The off-boat full-bag + // CUBE replay is the authoritative product, so bathy defaults to the `Processed` + // layer (#85; --bathy-layer overrides); the live node writes `Draft`. Single + // fused grid per layer (unh_marine_autonomy#221 — newest value wins per cell). + const std::size_t resident_before_final = accumulator.residentTileCount(); + const std::size_t evicted_count = accumulator.evictedIndices().size(); + accumulator.finalize( + registry, bs_store_dir.empty() ? nullptr : &bs_registry); + std::cout << "Persisted " << accumulator.bathyTilesPersisted() + << " bathy tile(s) to " << store_dir << " (" << bathy_layer_str + << " layer; " << evicted_count << " evicted mid-pass, " + << resident_before_final << " resident at end; build: " + << phase_secs() << "s)." << std::endl; + + if (accumulator.bathyTilesPersisted() == 0) { std::cerr << "WARNING: no tiles had finite data -- nothing imported. Check " "the projector frame overrides and the detections topic." << std::endl; } - // The off-boat full-bag CUBE replay is the authoritative product, so by default it - // lands in the `Processed` layer (#85; --bathy-layer overrides). The live node - // writes the `Draft` layer instead. Single fused grid per layer - // (unh_marine_autonomy#221 — no per-day epochs, newest value wins per cell). - store.importTiles(bathy_layer, std::move(tiles)); - - std::size_t written = marine_bathymetry_store::save(store, store_dir, ®istry); - std::cout << "Saved " << written << " tiles to " << store_dir << " (" - << bathy_layer_str << " layer)." << std::endl; - - // Optional: surface the co-estimated backscatter into an MBES backscatter store - // layer from the SAME CUBE pass (#80), written to the Processed layer (the - // off-boat CUBE re-run is the authoritative product). By default the value is - // UNCORRECTED; --backscatter-correction empirical (with a --backscatter-curve) - // applies the per-beam angular-response correction at node-output (cube#81), - // which corrects both this offline layer and the live tile (#78). Bathy (above) - // defaults to the same Processed layer (#85). + // The co-estimated backscatter was surfaced into the --bs-store layer (#80) from + // the SAME CUBE pass, incrementally under eviction (newest-finite-wins merge, + // cube#92). By default UNCORRECTED; --backscatter-correction empirical applies + // the per-beam angular-response correction at node-output (cube#81). if (!bs_store_dir.empty()) { - std::cout << "Building backscatter store tiles..." << std::endl; - - // Match the bathy store's GGGS level so both products tile identically. - marine_mbes_backscatter_store::MbesBackscatterStore bs_store = - marine_mbes_backscatter_store::MbesBackscatterStore::fromCellSize( - static_cast(geo_map_sheet.nominalCellSizeMeters())); - - // Provenance: register the same physical source in the backscatter registry, - // tagged with the backscatter sensor class. Every cell carries this index + - // the import timestamp, so the Processed product is not source/time-blank. - marine_mbes_backscatter_store::SourceRegistry bs_registry; - marine_mbes_backscatter_store::SourceRecord bs_source_record; - bs_source_record.source_id = source_record.source_id; - bs_source_record.platform = source_record.platform; - bs_source_record.sensor = source_record.sensor; - bs_source_record.sensor_class = "mbes-backscatter"; - bs_source_record.campaign = source_record.campaign; - const uint16_t bs_source_index = bs_registry.registerSource(bs_source_record); - - // Surface the co-estimated intensity cell-by-cell (no bulk-import API is added - // to the separate marine_mbes_backscatter_store package; #80 stays in-repo). - const std::map - bs_cells = cube::mapSheetToBackscatterCells( - geo_map_sheet, cell_timestamp_ns, bs_source_index); - std::cout << "Backscatter cells with data: " << bs_cells.size() - << " (build: " << phase_secs() << "s)" << std::endl; - - if (bs_cells.empty()) { + std::cout << "Persisted " << accumulator.backscatterTilesPersisted() + << " backscatter tile(s) to " << bs_store_dir << "." << std::endl; + if (accumulator.backscatterTilesPersisted() == 0) { std::cerr << "WARNING: no cells had finite backscatter -- nothing written to " "the backscatter store. Check that the detections carry intensities." << std::endl; } - - for (const auto & cell : bs_cells) { - bs_store.set( - marine_mbes_backscatter_store::SourceLayer::Processed, cell.first, cell.second); - } - - std::size_t bs_written = - marine_mbes_backscatter_store::save(bs_store, bs_store_dir, &bs_registry); - std::cout << "Saved " << bs_written << " backscatter tiles to " << bs_store_dir - << "." << std::endl; } std::cout << "done!" << std::endl; diff --git a/cube_bathymetry/src/node.cpp b/cube_bathymetry/src/node.cpp index 00296cf..8468422 100644 --- a/cube_bathymetry/src/node.cpp +++ b/cube_bathymetry/src/node.cpp @@ -21,6 +21,7 @@ #include "cube_bathymetry/node.h" +#include #include #include #include @@ -29,52 +30,6 @@ namespace cube { -namespace -{ - -/// Empirical angular-response lookup: the curve's db_relative_to_nadir at -/// @p abs_angle_deg, linearly interpolated between adjacent bin centres. The -/// curve is ascending {abs_angle_deg, db_relative_to_nadir} pairs (loaded by -/// loadAngularResponseCurve). Returns 0 (identity) when the curve is empty or -/// @p abs_angle_deg lies beyond the curve's max angle; clamps to the first bin -/// below the curve's min angle (nadir bin ~0 -> ~identity near nadir). -double curveRelativeDb( - const std::vector> & curve, double abs_angle_deg) -{ - if (curve.empty()) { - return 0.0; - } - // Below the first bin centre: clamp to it (the nadir bin's value, ~0). - if (abs_angle_deg <= curve.front().first) { - return curve.front().second; - } - // Beyond the last bin centre: clamp to the outermost bin's correction rather - // than jump to identity. The empirical curve is bounded (unlike a cos/log model - // that diverges near grazing), so continuing the edge value keeps the correction - // continuous and avoids a swath-edge discontinuity / bright ring (#81 review). - if (abs_angle_deg > curve.back().first) { - return curve.back().second; - } - // Find the bracketing pair [lo, hi] and linearly interpolate. - for (std::size_t i = 1; i < curve.size(); ++i) { - if (abs_angle_deg <= curve[i].first) { - const double a0 = curve[i - 1].first; - const double d0 = curve[i - 1].second; - const double a1 = curve[i].first; - const double d1 = curve[i].second; - const double span = a1 - a0; - if (span <= 0.0) { - return d1; // duplicate angle: take the upper bin's value - } - const double t = (abs_angle_deg - a0) / span; - return d0 + t * (d1 - d0); - } - } - return 0.0; // unreachable (guarded by the back() check above) -} - -} // namespace - bool Node::addHypothesis(float depth, float variance) { auto new_hypothesis = std::make_shared(depth, variance); @@ -131,14 +86,14 @@ bool Node::update( addHypothesis(depth, variance); // This beam seeded (and so is a member of) the new hypothesis: record its // backscatter on it. Without this the first beam at every node is dropped. - depth_hypotheses_.back()->recordBeam(intensity, beam_angle, range); + depth_hypotheses_.back()->recordBeam(intensity, beam_angle, range, parameters); return true; } else { /* Update the best hypothesis with the current data */ if(best->update(depth, variance, parameters)) { // Accepted into the winning hypothesis -- accumulate its backscatter on // the SAME hypothesis (ADR-0007 D2: intensity rides the winning depth). - best->recordBeam(intensity, beam_angle, range); + best->recordBeam(intensity, beam_angle, range, parameters); } else { /* Failed update --- indicates an intervention, so that we need to * start a new hypothesis to capture the outlier/datum shift. @@ -148,7 +103,7 @@ bool Node::update( // The beam was REJECTED from `best` (depth-geometry outlier) and used to // seed the new hypothesis, so its backscatter belongs to the new // hypothesis, NOT to `best` (exclusion-on-intervention invariant). - depth_hypotheses_.back()->recordBeam(intensity, beam_angle, range); + depth_hypotheses_.back()->recordBeam(intensity, beam_angle, range, parameters); } } return true; @@ -344,91 +299,72 @@ NodeRecord Node::extractNodeRecord(const Parameters & parameters) return record; } - // Backscatter half (ADR-0007 D2/D3/D4). Apply the per-beam angular-response - // correction to each retained {raw_intensity, beam_angle} sample, THEN combine - // the corrected values into a mean and an ESTIMATE variance. - // - // Empirical ARA (cube_bathymetry#81): corrected_dB = raw_dB - curveRel(|beam_angle|), - // where curveRel is the per-sonar angular-response curve's db_relative_to_nadir - // column, linearly interpolated between bin centres by |beam_angle| in DEGREES. - // Nadir -> curveRel ~0 -> identity. Beyond the curve's max angle, a NaN - // beam_angle, mode None, or an empty curve all fall back to identity (raw). - // - // rx_angles sign/zero gate (verified against kongsberg_em_bridge): the producer - // negates the Kongsberg pointing angle (+PORT -> +STARBOARD; nadir = 0; radians; - // intensity in dB). The curve is keyed on |beam_angle|, so port/starboard beams - // of equal magnitude get the same correction. - // - // Tier-2 (cube_bathymetry#87): when the loaded curve is a TL-REMOVED residual - // (backscatter_tl_removed), compensate each beam's 2-way transmission loss by - // ADDING IT BACK -- a distant return lost more energy, so it is boosted to - // recover range-independent backscatter (TVG-style): - // TL(R) = 40*log10(R) + 2*alpha*R (R = per-beam slant range, m) - // so the residual curve is depth/range transferable: - // corrected = raw + TL(R) - residualCurve(|beam_angle|). - // alpha is read verbatim from the curve header (backscatter_absorption_db_per_m); - // the estimator NEVER recomputes the (Francois-Garrison) absorption -- the - // Python derive tool is the single source of truth, guaranteeing consistency. - // A NaN / non-positive R skips the TL term (no log of a non-positive range); - // tier-1 (backscatter_tl_removed == false) skips it entirely (unchanged). - // - // Deferred (NOT done here): the full radiometric GeoCoder chain -- insonified - // area, beam-pattern, TVG residual, and the depth/slope incidence term - // (ADR-0007 D3, cube_bathymetry#15/#59). The empirical curve absorbs the - // aggregate angular falloff for a flat bottom; per-beam {raw_intensity, - // beam_angle, range} stay retained so the node value is re-derivable when #15 - // lands. - const bool apply_ara = - parameters.backscatter_angle_correction == - BackscatterAngleCorrection::Empirical && - !parameters.angular_response_curve.empty(); - const bool apply_tl = apply_ara && parameters.backscatter_tl_removed; - const double alpha = parameters.backscatter_absorption_db_per_m; - double sum = 0.0; - double sum_sq = 0.0; - uint32_t n = 0; - for(const auto & sample : chosen->intensity_samples) { - // recordBeam() already excluded NaN-intensity beams, so raw_intensity is real. - double corrected = sample.raw_intensity; - if(apply_tl) { - const double range = static_cast(sample.range); - // Skip the TL term for a missing / non-positive range (log10 undefined); - // the beam is still corrected by the residual angular-response curve below. - if(std::isfinite(range) && range > 0.0) { - // Compensate (ADD BACK) the 2-way transmission loss: a distant return - // lost more energy, so boost it to recover range-independent backscatter - // (TVG-style). #87 sign fix. - corrected += 40.0 * std::log10(range) + 2.0 * alpha * range; - } - } - if(apply_ara && !std::isnan(sample.beam_angle)) { - const double abs_angle_deg = - std::abs(static_cast(sample.beam_angle)) * 180.0 / M_PI; - corrected -= curveRelativeDb(parameters.angular_response_curve, abs_angle_deg); + // Backscatter half (ADR-0007 D2/D3/D4). The angular-response + tier-2 TL + // correction is now applied at RECORD (recordBeam -> correctBeamIntensity, + // cube#93), and the winning hypothesis streams a Welford of the CORRECTED + // intensity. So extract just reads (n, mean, M2) -- NO re-correction (the curve + // params are not consulted here anymore). This produces the SAME node-output + // mean + estimate variance as the old retain-and-correct-at-extract path; the + // correction math itself is unchanged, only its timing moved (record <- extract). + const IntensityWelford & bs = chosen->intensity; + record.n_samples = bs.n; + if(bs.n > 0) { + record.intensity = static_cast(bs.mean); + if(bs.n >= 2) { + // Unbiased sample variance M2/(n-1), then divided by n to get the variance + // of the MEAN -- the estimate variance that shrinks with n (ADR-0007 D4), + // matching the bathy store's depth uncertainty. Clamp to 0 defensively (the + // Welford M2 is >= 0 by construction, so this is a no-op in practice). + const double sample_variance = std::max(0.0, bs.m2 / (bs.n - 1)); + record.intensity_var = static_cast(sample_variance / bs.n); } - sum += corrected; - sum_sq += corrected * corrected; - ++n; + // n == 1: intensity set, intensity_var stays NaN (no spread from one beam). } - record.n_samples = n; - if(n > 0) { - const double mean = sum / n; - record.intensity = static_cast(mean); - if(n >= 2) { - // Sample variance (unbiased), then divide by n to get the variance of the - // MEAN -- the estimate variance that shrinks with n (ADR-0007 D4), NOT the - // raw sample variance. Mirrors the bathy store's depth uncertainty. - // Clamp to 0: float rounding on the sum-of-squares form can produce a - // tiny negative sample_variance at low-dB means over many beams. - const double sample_variance = std::max( - 0.0, (sum_sq - sum * sum / n) / (n - 1)); - record.intensity_var = static_cast(sample_variance / n); + return record; +} + +IntensityWelford Node::chosenIntensityWelford() +{ + // Same hypothesis selection as extractNodeRecord(), so the spilled Welford and + // the persisted node-output intensity come from the SAME hypothesis (#92/#93). + std::shared_ptr chosen; + if(nominated_hypothesis_) { + chosen = nominated_hypothesis_; + } else { + auto h = chooseHypothesis(); + if(h && h->number_of_samples > 0) { + chosen = h; } - // n == 1: intensity set, intensity_var stays NaN (no spread from one beam). } + if(!chosen) { + return {}; + } + return chosen->intensity; +} - return record; +void Node::setSettledIntensityWelford(const IntensityWelford & intensity) +{ + // Target the same hypothesis chosenIntensityWelford()/extractNodeRecord() read. + // On the eviction-reload path the node is freshly reseeded (seedSettledDepth + // pushed exactly one hypothesis), so this is that hypothesis; the revisit's + // beams then continue the Welford on it via recordBeam(). Because (n, mean, M2) + // is a perfect sufficient statistic, restore-then-continue is bit-identical to + // never-evicting (#93). + std::shared_ptr chosen; + if(nominated_hypothesis_) { + chosen = nominated_hypothesis_; + } else { + // Same number_of_samples > 0 gate chosenIntensityWelford()/extractNodeRecord() + // apply, so the restore targets exactly the hypothesis the spill read from. + auto h = chooseHypothesis(); + if(h && h->number_of_samples > 0) { + chosen = h; + } + } + if(chosen) { + chosen->intensity = intensity; + } } std::shared_ptr Node::chooseHypothesis() diff --git a/cube_bathymetry/src/store_import.cpp b/cube_bathymetry/src/store_import.cpp index 5f9e392..7fab67a 100644 --- a/cube_bathymetry/src/store_import.cpp +++ b/cube_bathymetry/src/store_import.cpp @@ -23,12 +23,30 @@ #include "cube_bathymetry/store_import.h" #include +#include +#include #include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include #include #include "cube_bathymetry/common.h" +#include "cube_bathymetry/hypothesis.h" #include "marine_autonomy/gggs.h" +#include "marine_bathymetry_store/tile_io.hpp" +#include "marine_mbes_backscatter_store/mbes_store.hpp" +#include "marine_mbes_backscatter_store/mbes_tile.hpp" +#include "marine_mbes_backscatter_store/tile_io.hpp" namespace cube { @@ -212,4 +230,378 @@ void loadIntoSheet( } } +// =========================================================================== +// ImportAccumulator (cube_bathymetry#92): bounded-RAM offline import. +// =========================================================================== + +namespace +{ +// Per-tile intensity-Welford spill (cube#92/#93). A tiny binary file, scratch-only +// and single-machine, so native layout/endianness is fine. Format: +// uint32 magic | uint32 ncells | +// ncells * { uint16 row, uint16 col, uint32 n, double mean, double m2 } +// Each cell's record is a fixed 24 bytes -- O(1)/cell (cube#93 replaced the +// variable-length raw-sample spill with the 3-number corrected-intensity Welford). +constexpr uint32_t kSpillMagic = 0x43425332u; // "CBS2" (cube#93 triplet format) + +std::string spillFileName(const gggs::GridIndex & index) +{ + // Reuse the bathy tile naming (level_row_col.tif) for a unique per-tile stem. + return marine_bathymetry_store::tileFilename(index) + ".spill"; +} +} // namespace + +ImportAccumulator::ImportAccumulator(GeoMapSheet & sheet, ImportAccumulatorConfig config) +: sheet_(sheet), cfg_(std::move(config)) +{ +} + +ImportAccumulator::~ImportAccumulator() +{ + // RAII safety net: finalize() normally cleans up, but an exception on the import + // path must not leak the scratch spill. + cleanupScratch(); +} + +void ImportAccumulator::cleanupScratch() +{ + if (!scratch_dir_.empty()) { + std::error_code ec; + std::filesystem::remove_all(scratch_dir_, ec); // best-effort; never throws here + scratch_dir_.clear(); + } +} + +std::size_t ImportAccumulator::residentTileCount() const +{ + return sheet_.residentTileCount(); +} + +void ImportAccumulator::spillIntensitySamples(const gggs::GridIndex & index) +{ + if (cfg_.bs_store_dir.empty()) { + return; // no backscatter product -> no need to retain the intensity Welford + } + auto grid = sheet_.gridAt(index); + if (!grid) { + return; + } + const std::map welford = + grid->nodeIntensityWelford(); + + // Lazily create a unique scratch dir on first spill. mkdtemp atomically creates + // a directory with a name guaranteed unique against any other process (it retries + // internally on collision), so two concurrent import_bag processes can never + // share a spill dir and cross-corrupt each other's tiles. A steady_clock-ns + + // per-process atomic counter name could collide cross-process: the counter resets + // per process, leaving the ns the only distinguisher, and two processes can read + // the same coarse ns tick. + if (scratch_dir_.empty()) { + std::string tmpl = + (std::filesystem::temp_directory_path() / "cube_import_spill_XXXXXX").string(); + std::vector buf(tmpl.begin(), tmpl.end()); + buf.push_back('\0'); + if (::mkdtemp(buf.data()) == nullptr) { + throw std::runtime_error( + std::string("import_bag: cannot create spill scratch dir from template ") + + tmpl + ": " + std::strerror(errno)); + } + scratch_dir_ = std::string(buf.data()); + } + + const std::string path = scratch_dir_ + "/" + spillFileName(index); + if (welford.empty()) { + // Nothing to retain; drop any stale spill so a later reload can't restore + // outdated state for this tile. + std::error_code ec; + std::filesystem::remove(path, ec); + return; + } + + std::ofstream out(path, std::ios::binary | std::ios::trunc); + if (!out) { + throw std::runtime_error("import_bag: cannot open spill file " + path); + } + const uint32_t magic = kSpillMagic; + const uint32_t ncells = static_cast(welford.size()); + out.write(reinterpret_cast(&magic), sizeof(magic)); + out.write(reinterpret_cast(&ncells), sizeof(ncells)); + for (const auto & entry : welford) { + const uint16_t row = entry.first.row(); + const uint16_t col = entry.first.column(); + const uint32_t n = entry.second.n; + const double mean = entry.second.mean; + const double m2 = entry.second.m2; + out.write(reinterpret_cast(&row), sizeof(row)); + out.write(reinterpret_cast(&col), sizeof(col)); + out.write(reinterpret_cast(&n), sizeof(n)); + out.write(reinterpret_cast(&mean), sizeof(mean)); + out.write(reinterpret_cast(&m2), sizeof(m2)); + } + if (!out) { + throw std::runtime_error("import_bag: failed writing spill file " + path); + } +} + +void ImportAccumulator::persistBathyTile(const gggs::GridIndex & index) +{ + if (cfg_.store_dir.empty()) { + return; + } + auto grid = sheet_.gridAt(index); + if (!grid) { + return; + } + // geoGridToTile flushes the median pre-filter (values()) and writes only the + // finite cells; an all-no-data tile is not persisted (matches mapSheetToTiles). + marine_bathymetry_store::BathymetryTile tile = + geoGridToTile(*grid, cfg_.timestamp_ns, cfg_.source_index); + if (!tile.dirty()) { + return; + } + const std::string layer_dir = cfg_.store_dir + "/" + + marine_bathymetry_store::layerDirName(cfg_.bathy_layer); + std::filesystem::create_directories(layer_dir); + // Atomic temp-then-rename via tile_io::saveTile -- a partial write never + // corrupts the on-disk surface (same primitive the live node's eviction uses). + marine_bathymetry_store::saveTile( + tile, layer_dir + "/" + marine_bathymetry_store::tileFilename(index)); + ++bathy_persisted_; +} + +void ImportAccumulator::persistBackscatterTile(const gggs::GridIndex & index) +{ + if (cfg_.bs_store_dir.empty()) { + return; + } + auto grid = sheet_.gridAt(index); + if (!grid) { + return; + } + namespace mbs = marine_mbes_backscatter_store; + // Only the finite co-estimated cells (NaN-intensity cells are skipped upstream). + const std::map cells = + geoGridToBackscatterCells(*grid, cfg_.timestamp_ns, cfg_.bs_source_index); + if (cells.empty()) { + // No finite intensity this pass -> nothing to write. Do NOT write an empty + // tile over a populated one. With the sample spill/restore this is only hit + // for a genuinely intensity-less tile (it never loses an earlier write). + return; + } + const std::string layer_dir = cfg_.bs_store_dir + "/" + + mbs::layerDirName(mbs::SourceLayer::Processed); + const std::string path = layer_dir + "/" + mbs::tileFilename(index); + // Plain overwrite (no on-disk merge): the reload-before-add path restores the + // tile's raw intensity samples before a revisit accretes onto them, so at every + // eviction the in-RAM tile already holds the COMPLETE sample population for each + // cell -- its corrected summary is the complete value, and overwriting the + // on-disk tile reproduces it losslessly. (The old newest-finite-wins disk merge + // was only needed when reload could not restore intensity; #92 removes that gap.) + mbs::MbesTile tile(index); + for (const auto & entry : cells) { + tile.set(entry.first.row(), entry.first.column(), entry.second); + } + std::filesystem::create_directories(layer_dir); + mbs::saveTile(tile, path); + ++bs_persisted_; +} + +void ImportAccumulator::restoreSpilledSamples(const gggs::GridIndex & index) +{ + // Best-effort: a missing/corrupt/truncated spill only degrades backscatter for + // this tile to the post-reload samples; the depth reload already succeeded. + if (scratch_dir_.empty()) { + return; + } + const std::string path = scratch_dir_ + "/" + spillFileName(index); + if (!std::filesystem::is_regular_file(path)) { + return; + } + std::ifstream in(path, std::ios::binary); + if (!in) { + return; + } + uint32_t magic = 0; + uint32_t ncells = 0; + in.read(reinterpret_cast(&magic), sizeof(magic)); + in.read(reinterpret_cast(&ncells), sizeof(ncells)); + if (!in || magic != kSpillMagic) { + return; + } + for (uint32_t c = 0; c < ncells; ++c) { + uint16_t row = 0; + uint16_t col = 0; + IntensityWelford w; + in.read(reinterpret_cast(&row), sizeof(row)); + in.read(reinterpret_cast(&col), sizeof(col)); + in.read(reinterpret_cast(&w.n), sizeof(w.n)); + in.read(reinterpret_cast(&w.mean), sizeof(w.mean)); + in.read(reinterpret_cast(&w.m2), sizeof(w.m2)); + if (!in) { + return; // truncated record -> stop (what was restored already stands) + } + sheet_.setSettledIntensityWelfordAt(gggs::CellIndex(index, row, col), w); + } +} + +bool ImportAccumulator::reloadEvictedTile(const gggs::GridIndex & index) +{ + if (cfg_.store_dir.empty()) { + return true; + } + try { + marine_bathymetry_store::BathymetryStore scratch = + marine_bathymetry_store::BathymetryStore::fromCellSize(cfg_.cell_size_m); + const auto sw = index.southWestPosition(); + const auto ne = index.northEastPosition(); + marine_bathymetry_store::loadWindow(scratch, cfg_.store_dir, sw, ne, nullptr); + const auto & tiles = scratch.tiles(cfg_.bathy_layer); + auto it = tiles.find(index); + if (it != tiles.end()) { + // seed_settled=true: restore each cell's settled depth/uncertainty as a CUBE + // hypothesis so accumulation continues from the saved state (ADR-0001). + primeFromTile(it->second, sheet_); + } + } catch (const std::exception & e) { + // The on-disk surface is the real data: on a load error the caller drops the + // partial re-created grid rather than let a later save clobber the intact file. + // The tile stays evicted, so a later batch that revisits it retries the reload; + // but this batch's soundings on the tile are dropped now and not re-processed + // (the caller WARNs with the count). + std::cerr << "import_bag: could not reload evicted tile on revisit: " + << e.what() << " (dropping this batch's soundings on the tile to " + "protect the on-disk surface; a later batch revisiting the tile retries the " + "reload, but these dropped soundings are not re-processed)" << std::endl; + return false; + } + // Restore the raw intensity samples onto the just-reseeded hypotheses so the + // revisit's beams blend with the pre-eviction population (lossless backscatter). + // Runs AFTER the depth reseed and BEFORE the batch's soundings are added. + restoreSpilledSamples(index); + return true; +} + +void ImportAccumulator::evictColdTiles() +{ + if (cfg_.max_resident_tiles == 0) { + return; // unbounded: never evict (preserves the pre-#92 behavior) + } + if (cfg_.store_dir.empty()) { + return; // nowhere to persist -> never evict (dropping would lose data) + } + if (sheet_.residentTileCount() <= cfg_.max_resident_tiles) { + return; + } + // Persist-then-drop the coldest tiles beyond budget. Each is written to disk + // (bathy + backscatter summary) and its corrected-intensity Welford spilled to + // scratch BEFORE it is dropped, so eviction is lossless and the tile reloads + // (with its intensity sufficient statistic) on revisit. A tile whose persist/spill + // THROWS is left + // resident (never dropped): RAM stays transiently over budget rather than losing + // unsaved data; the next batch retries. + for (const auto & index : sheet_.coldTiles(cfg_.max_resident_tiles)) { + try { + persistBathyTile(index); + persistBackscatterTile(index); + spillIntensitySamples(index); + } catch (const std::exception & e) { + std::cerr << "import_bag: failed to persist cold tile for eviction: " + << e.what() << " (keeping it resident to avoid data loss)" + << std::endl; + continue; // do NOT drop -- lossless guarantee + } + sheet_.dropTile(index); + evicted_.insert(index); + } +} + +void ImportAccumulator::addBatch( + const std::vector & soundings, + std::chrono::steady_clock::time_point time) +{ + // Reload-BEFORE-add (cube#92 lossless blend): reload any evicted tile this batch + // is about to touch FIRST, so the batch's beams accrete onto the reloaded + // hypotheses (settled depth + restored corrected-intensity Welford) instead of + // forming a fresh, partial tile. gridIndicesForSoundings computes the SAME + // expanded window addSoundings will touch (including a near-seam neighbour tile), + // so no evicted tile is missed. A tile whose reload FAILS is recorded and, AFTER + // the add, has this batch's soundings on it dropped to protect the intact on-disk + // surface. The tile stays evicted, so a LATER batch that revisits it retries the + // reload -- but the soundings dropped here are NOT replayed (offline is a single + // pass over the bag), so this is a real, bounded loss, not the "lossless" + // eviction blend; it is counted and WARNed below rather than hidden. + std::vector reload_failed; + if (!cfg_.store_dir.empty() && !evicted_.empty()) { + for (const auto & idx : sheet_.gridIndicesForSoundings(soundings)) { + if (evicted_.count(idx)) { + if (reloadEvictedTile(idx)) { + evicted_.erase(idx); + } else { + reload_failed.push_back(idx); + } + } + } + } + + sheet_.addSoundings(soundings, time); + + if (!reload_failed.empty()) { + // Count and report the loss before dropping: the tiles that failed to reload + // are about to be discarded together with the soundings this batch put on them, + // and offline those soundings are never replayed. Count soundings whose centre + // falls in a failed tile (a sounding can also spread into neighbour tiles, so + // this is a conservative floor, not an exact cell tally -- see + // GeoMapSheet::gridIndexForSounding). + const std::set failed_set( + reload_failed.begin(), reload_failed.end()); + std::size_t dropped_soundings = 0; + for (const auto & s : soundings) { + if (failed_set.count(sheet_.gridIndexForSounding(s))) { + ++dropped_soundings; + } + } + std::cerr << "import_bag: WARNING permanently dropping ~" << dropped_soundings + << " sounding(s) centred in " << reload_failed.size() + << " un-reloadable tile(s) this batch (reload failed; the on-disk " + "surface is preserved, but these soundings are NOT re-processed " + "offline -- not the lossless eviction blend)" << std::endl; + for (const auto & idx : reload_failed) { + sheet_.dropTile(idx); // protect the intact on-disk surface + } + } + + evictColdTiles(); +} + +void ImportAccumulator::finalize( + const marine_bathymetry_store::SourceRegistry & bathy_registry, + const marine_mbes_backscatter_store::SourceRegistry * bs_registry) +{ + // Persist every still-resident tile (evicted tiles are already durable). Snapshot + // the indices first so the persist loop iterates a stable, deterministic order. + std::vector resident; + for (const auto & grid : sheet_.grids()) { + if (grid) { + resident.push_back(grid->index()); + } + } + for (const auto & index : resident) { + persistBathyTile(index); + persistBackscatterTile(index); + } + // Registries are store-wide sidecars: written once at the end so every cell's + // source_index resolves (the per-tile eviction writes carry no registry). + if (!cfg_.store_dir.empty()) { + std::filesystem::create_directories(cfg_.store_dir); + bathy_registry.saveRegistry(cfg_.store_dir); + } + if (!cfg_.bs_store_dir.empty() && bs_registry != nullptr) { + std::filesystem::create_directories(cfg_.bs_store_dir); + bs_registry->saveRegistry(cfg_.bs_store_dir); + } + // The spilled raw samples were only needed to reload an evicted tile mid-run; + // the import is complete, so delete the scratch dir. + cleanupScratch(); +} + } // namespace cube diff --git a/cube_bathymetry/test/test_hypothesis.cpp b/cube_bathymetry/test/test_hypothesis.cpp index 5603bad..2106510 100644 --- a/cube_bathymetry/test/test_hypothesis.cpp +++ b/cube_bathymetry/test/test_hypothesis.cpp @@ -201,46 +201,72 @@ TEST_F(HypothesisTest, MonitorIsSymmetricInErrorSign) EXPECT_EQ(h_plus.sequence_length, h_minus.sequence_length); } -// ADR-0007 D3: recordBeam appends the per-beam {raw intensity, grazing angle} -// sufficient-statistics pair, keeping the raw value re-correctable at output. -TEST_F(HypothesisTest, RecordBeamAppendsRawAndAngle) +// ADR-0007 D3/D4, cube#93: recordBeam folds each beam's CORRECTED intensity into +// the streaming Welford (n, mean, M2). With the default None correction the +// corrected value equals the raw, so two beams give the running mean of the raw +// dB. (White-box storage changed from a raw-sample vector to the O(1) Welford; +// the node-output mean/variance these feed are unchanged -- see test_node.) +TEST_F(HypothesisTest, RecordBeamFoldsCorrectedIntensityIntoWelford) { Hypothesis h(10.0f, 1.0f); - EXPECT_TRUE(h.intensity_samples.empty()); + EXPECT_EQ(h.intensity.n, 0u); - h.recordBeam(-30.0f, 0.1f); - h.recordBeam(-28.0f, 0.2f); + h.recordBeam(-30.0f, 0.1f, std::nan(""), params); // None mode -> corrected = raw + h.recordBeam(-28.0f, 0.2f, std::nan(""), params); - ASSERT_EQ(h.intensity_samples.size(), 2u); - EXPECT_FLOAT_EQ(h.intensity_samples[0].raw_intensity, -30.0f); - EXPECT_FLOAT_EQ(h.intensity_samples[0].beam_angle, 0.1f); - EXPECT_FLOAT_EQ(h.intensity_samples[1].raw_intensity, -28.0f); - EXPECT_FLOAT_EQ(h.intensity_samples[1].beam_angle, 0.2f); + EXPECT_EQ(h.intensity.n, 2u); + EXPECT_DOUBLE_EQ(h.intensity.mean, -29.0); // (-30 + -28) / 2 + // M2 = sum (x - mean)^2 = 1 + 1 = 2 -> sample variance M2/(n-1) = 2. + EXPECT_DOUBLE_EQ(h.intensity.m2, 2.0); } // A NaN raw intensity (source omitted intensities) must never become a phantom -// backscatter sample. +// backscatter sample -- it leaves the Welford count unchanged. TEST_F(HypothesisTest, RecordBeamSkipsNanIntensity) { Hypothesis h(10.0f, 1.0f); - h.recordBeam(std::nan(""), 0.1f); - EXPECT_TRUE(h.intensity_samples.empty()); + h.recordBeam(std::nan(""), 0.1f, std::nan(""), params); + EXPECT_EQ(h.intensity.n, 0u); - h.recordBeam(-25.0f, 0.3f); - EXPECT_EQ(h.intensity_samples.size(), 1u); + h.recordBeam(-25.0f, 0.3f, std::nan(""), params); + EXPECT_EQ(h.intensity.n, 1u); + EXPECT_DOUBLE_EQ(h.intensity.mean, -25.0); } -// A NaN beam angle is retained: the beam is still a valid intensity sample, -// it simply cannot be angle-corrected (emitted uncorrected at output). -TEST_F(HypothesisTest, RecordBeamRetainsNanAngle) +// A NaN beam angle is still a valid intensity sample: it is counted, just not +// angle-corrected (with None mode there is no angle term anyway -> corrected = raw). +TEST_F(HypothesisTest, RecordBeamCountsNanAngle) { Hypothesis h(10.0f, 1.0f); - h.recordBeam(-20.0f, std::nan("")); - ASSERT_EQ(h.intensity_samples.size(), 1u); - EXPECT_FLOAT_EQ(h.intensity_samples[0].raw_intensity, -20.0f); - EXPECT_TRUE(std::isnan(h.intensity_samples[0].beam_angle)); + h.recordBeam(-20.0f, std::nan(""), std::nan(""), params); + EXPECT_EQ(h.intensity.n, 1u); + EXPECT_DOUBLE_EQ(h.intensity.mean, -20.0); +} + +// cube#93 O(1) memory: per-cell intensity state is a fixed-size Welford triplet, +// NOT a per-beam vector, so a heavily-oversampled cell (the Massabesic OOM cause: +// ~10^9 beams over 47 tiles) does not grow per-cell RAM. The portable proxy: fold +// a large number of beams and assert the per-hypothesis footprint is unchanged +// (the old raw-sample vector would have grown ~12 bytes/beam without bound). +TEST_F(HypothesisTest, IntensityStateIsConstantSizeRegardlessOfBeamCount) +{ + const std::size_t hypothesis_size = sizeof(Hypothesis); + + Hypothesis h(10.0f, 1.0f); + const int many = 500000; + for (int i = 0; i < many; ++i) { + // Vary the value so it is a real running estimate, not a constant fold. + h.recordBeam(-30.0f - static_cast(i % 11), 0.0f, std::nan(""), params); + } + + // All beams folded into the O(1) accumulator... + EXPECT_EQ(h.intensity.n, static_cast(many)); + // ...with NO growth in the per-hypothesis footprint (the accumulator is a value + // member of fixed size; this is the cube#93 invariant the old vector violated). + EXPECT_EQ(sizeof(Hypothesis), hypothesis_size); + EXPECT_LE(sizeof(IntensityWelford), static_cast(32)); } } // namespace cube diff --git a/cube_bathymetry/test/test_import_eviction.cpp b/cube_bathymetry/test/test_import_eviction.cpp new file mode 100644 index 0000000..293794b --- /dev/null +++ b/cube_bathymetry/test/test_import_eviction.cpp @@ -0,0 +1,403 @@ +// Copyright 2026 Center for Coastal and Ocean Mapping & NOAA-UNH Joint +// Hydrographic Center, University of New Hampshire +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +// THE SOFTWARE. + +// Bounded == unbounded equivalence for the offline import accumulator (#92). +// +// The proof that the cube#70 persist-then-drop + reload-on-revisit scheme, ported +// from the live node to import_bag, is LOSSLESS: a survey run with a tiny resident +// budget (forcing eviction + reload) must produce the SAME on-disk bathy store +// (depth + uncertainty) AND backscatter store (intensity + variance) as the same +// survey run unbounded (everything in RAM, persisted only at the end). A +// data-losing eviction would drop or NaN-out cells and the comparison would fail. + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "cube_bathymetry/geo_map_sheet.h" +#include "cube_bathymetry/geo_sounding.h" +#include "cube_bathymetry/store_import.h" +#include "marine_autonomy/gggs.h" +#include "marine_bathymetry_store/bathymetry_store.hpp" +#include "marine_bathymetry_store/bathymetry_tile.hpp" +#include "marine_bathymetry_store/registry.hpp" +#include "marine_bathymetry_store/tile_io.hpp" +#include "marine_mbes_backscatter_store/mbes_store.hpp" +#include "marine_mbes_backscatter_store/registry.hpp" +#include "marine_mbes_backscatter_store/tile_io.hpp" + +namespace cube +{ +namespace +{ +constexpr int64_t kStamp = 1234567890123456789LL; +constexpr float kCellSize = 1.0f; + +std::string makeTempDir(const std::string & tag) +{ + const auto base = std::filesystem::temp_directory_path() / + ("cube_import_evict_" + tag + "_" + std::to_string(::getpid())); + std::filesystem::remove_all(base); + std::filesystem::create_directories(base); + return base.string(); +} + +// One cell's worth of survey: several soundings at one position so CUBE settles a +// finite depth AND co-estimates a finite intensity (a single sample may not +// resolve a variance). A constant intensity keeps the Welford mean exact. +std::vector surveyCell( + double lat, double lon, float depth, float intensity, int reps = 8) +{ + std::vector soundings; + for (int rep = 0; rep < reps; ++rep) { + gz4d::GeoPointLatLongDegrees point(lat, lon, -depth - rep * 0.001); + GeoSounding s(point); + s.sounding.vertical_error = 0.5f; + s.sounding.horizontal_error = 0.1f; + s.sounding.intensity = intensity; + s.sounding.beam_angle = 0.0f; + s.sounding.slant_range = depth; + soundings.push_back(s); + } + return soundings; +} + +ImportAccumulatorConfig makeConfig( + const std::string & store_dir, const std::string & bs_dir, std::size_t budget) +{ + ImportAccumulatorConfig cfg; + cfg.store_dir = store_dir; + cfg.bathy_layer = marine_bathymetry_store::SourceLayer::Processed; + cfg.source_index = 1; + cfg.timestamp_ns = kStamp; + cfg.cell_size_m = kCellSize; + cfg.bs_store_dir = bs_dir; + cfg.bs_source_index = 1; + cfg.max_resident_tiles = budget; + return cfg; +} + +// Drive one full import: feed every batch through the accumulator, then finalize. +void runImport( + const std::vector> & batches, + const std::string & store_dir, const std::string & bs_dir, std::size_t budget) +{ + GeoMapSheet sheet(kCellSize); + ImportAccumulator accumulator(sheet, makeConfig(store_dir, bs_dir, budget)); + for (const auto & batch : batches) { + accumulator.addBatch(batch); + } + marine_bathymetry_store::SourceRegistry bathy_reg; + marine_bathymetry_store::SourceRecord bathy_rec; + bathy_rec.source_id = "test"; + bathy_reg.registerSource(bathy_rec); + marine_mbes_backscatter_store::SourceRegistry bs_reg; + marine_mbes_backscatter_store::SourceRecord bs_rec; + bs_rec.source_id = "test"; + bs_reg.registerSource(bs_rec); + accumulator.finalize(bathy_reg, bs_dir.empty() ? nullptr : &bs_reg); +} + +// Load every finite bathy cell of a store layer into a comparable map. +std::map> loadBathyCells( + const std::string & store_dir) +{ + marine_bathymetry_store::BathymetryStore store = + marine_bathymetry_store::BathymetryStore::fromCellSize(kCellSize); + marine_bathymetry_store::SourceRegistry reg; + marine_bathymetry_store::load(store, store_dir, ®); + std::map> out; + for (const auto & grid_tile : store.tiles( + marine_bathymetry_store::SourceLayer::Processed)) + { + const auto & depth = grid_tile.second.depthBand(); + const auto & unc = grid_tile.second.uncertaintyBand(); + gggs::CellAreaIterator it(grid_tile.second.index()); + std::size_t k = 0; + for (; it.valid() && k < depth.size(); it.next(), ++k) { + if (std::isfinite(depth[k])) { + out.emplace(*it, std::make_pair(depth[k], unc[k])); + } + } + } + return out; +} + +// Load every finite backscatter cell of a store into a comparable map. +std::map> loadBackscatterCells( + const std::string & bs_dir) +{ + marine_mbes_backscatter_store::MbesBackscatterStore store = + marine_mbes_backscatter_store::MbesBackscatterStore::fromCellSize(kCellSize); + marine_mbes_backscatter_store::SourceRegistry reg; + marine_mbes_backscatter_store::load(store, bs_dir, ®); + std::map> out; + for (const auto & grid_tile : store.tiles( + marine_mbes_backscatter_store::SourceLayer::Processed)) + { + gggs::CellAreaIterator it(grid_tile.second.index()); + for (; it.valid(); it.next()) { + const marine_mbes_backscatter_store::MbesCell c = + grid_tile.second.get((*it).row(), (*it).column()); + if (c.hasData()) { + out.emplace(*it, std::make_pair(c.intensity, c.intensity_variance)); + } + } + } + return out; +} + +// Compare two store dirs cell-for-cell. Cell COUNT must match exactly (a dropped +// or NaN-ed cell is data loss). Depth and backscatter (intensity + variance) are +// compared tightly: eviction+reload+merge must reproduce them, not approximate. +// +// @param uncertainty_tol Bathy-uncertainty tolerance. For a cell that is evicted +// and never revisited (or never evicted), the uncertainty is bit-reproduced +// (pass 1e-3). For a cell reloaded on revisit, the warm-start restores the +// settled DEPTH losslessly but RE-DERIVES the uncertainty from the reseeded +// single hypothesis rather than the original sample history (ADR-0001 / #21: +// "does NOT reconstruct CUBE hypothesis, queue, or pre-filter state"), so its +// uncertainty drifts slightly — a documented, depth-preserving artifact, not +// data loss. Such tests pass a looser bound. +void expectStoresEqual( + const std::string & bounded_dir, const std::string & unbounded_dir, + const std::string & bounded_bs, const std::string & unbounded_bs, + double uncertainty_tol) +{ + const auto bathy_b = loadBathyCells(bounded_dir); + const auto bathy_u = loadBathyCells(unbounded_dir); + ASSERT_FALSE(bathy_u.empty()) << "unbounded build produced no bathy cells"; + ASSERT_EQ(bathy_b.size(), bathy_u.size()) + << "bounded build lost or gained bathy cells"; + for (const auto & [cell, du] : bathy_u) { + auto it = bathy_b.find(cell); + ASSERT_NE(it, bathy_b.end()) << "bounded build is missing a bathy cell"; + EXPECT_NEAR(it->second.first, du.first, 1e-3) << "depth differs"; + EXPECT_NEAR(it->second.second, du.second, uncertainty_tol) + << "uncertainty differs"; + } + + const auto bs_b = loadBackscatterCells(bounded_bs); + const auto bs_u = loadBackscatterCells(unbounded_bs); + ASSERT_FALSE(bs_u.empty()) << "unbounded build produced no backscatter cells"; + ASSERT_EQ(bs_b.size(), bs_u.size()) + << "bounded build lost or gained backscatter cells"; + for (const auto & [cell, iv] : bs_u) { + auto it = bs_b.find(cell); + ASSERT_NE(it, bs_b.end()) << "bounded build is missing a backscatter cell"; + EXPECT_NEAR(it->second.first, iv.first, 1e-3) << "intensity differs"; + EXPECT_NEAR(it->second.second, iv.second, 1e-3) << "intensity_variance differs"; + } +} +} // namespace + +// A disjoint multi-tile survey: a tiny budget forces every tile to be evicted +// once (and never revisited), so the bounded store is the union of the evicted +// tiles. It must equal the unbounded store exactly. +TEST(ImportEviction, BoundedEqualsUnboundedNoRevisit) +{ + std::vector> batches; + const int n_tiles = 12; + for (int i = 0; i < n_tiles; ++i) { + // 0.02 deg (~2.2 km) apart -> every batch lands on a distinct GGGS tile. + batches.push_back( + surveyCell(43.0 + 0.02 * i, -70.0, 10.0f + i, 30.0f + i)); + } + + const std::string root = makeTempDir("norevisit"); + const std::string b_dir = root + "/bounded"; + const std::string u_dir = root + "/unbounded"; + const std::string b_bs = root + "/bounded_bs"; + const std::string u_bs = root + "/unbounded_bs"; + + runImport(batches, b_dir, b_bs, /*budget=*/3); // forces eviction + runImport(batches, u_dir, u_bs, /*budget=*/0); // unbounded + + // No revisit: every cell is bit-reproduced, uncertainty included. + expectStoresEqual(b_dir, u_dir, b_bs, u_bs, /*uncertainty_tol=*/1e-3); + std::filesystem::remove_all(root); +} + +// A tile surveyed, evicted (by surveying many other tiles), then revisited at +// DIFFERENT cells within it. The reload-on-revisit must restore the first visit's +// settled cells (bathy) and the on-disk merge must preserve their backscatter, so +// the final store has BOTH visits' cells -- identical to the unbounded build. +TEST(ImportEviction, RevisitEqualsUnbounded) +{ + std::vector> batches; + // Visit 1: tile A (cell near its SW). + batches.push_back(surveyCell(43.00000, -70.00000, 12.0f, 40.0f)); + // Evict A by surveying 10 distinct far-away tiles. + for (int i = 1; i <= 10; ++i) { + batches.push_back(surveyCell(43.0 + 0.02 * i, -71.0, 15.0f + i, 25.0f + i)); + } + // Visit 2: revisit tile A at a DIFFERENT cell (~55 m N, same ~960 m tile). + batches.push_back(surveyCell(43.00050, -70.00000, 18.0f, 55.0f)); + + const std::string root = makeTempDir("revisit"); + const std::string b_dir = root + "/bounded"; + const std::string u_dir = root + "/unbounded"; + const std::string b_bs = root + "/bounded_bs"; + const std::string u_bs = root + "/unbounded_bs"; + + runImport(batches, b_dir, b_bs, /*budget=*/3); // A is evicted then reloaded + runImport(batches, u_dir, u_bs, /*budget=*/0); // A never leaves RAM + + // No cell is lost (count matches) and every depth + backscatter value is + // reproduced. The first visit's reloaded cells keep their depth losslessly but + // re-derive uncertainty from the warm-start reseed (ADR-0001), so uncertainty + // is compared with a looser bound. + expectStoresEqual(b_dir, u_dir, b_bs, u_bs, /*uncertainty_tol=*/0.05); + std::filesystem::remove_all(root); +} + +// Lossless backscatter blend (cube#92 Option A): a cell resurveyed AFTER its tile +// was evicted blends BOTH visits' backscatter, matching a never-evicted build. +// The raw per-beam intensity samples are spilled to scratch on eviction and +// restored onto the reloaded hypothesis BEFORE the revisit's beams accrete, so the +// node-output intensity is the full population, not just the newest visit. Both +// visits survey the SAME depth so they land on one hypothesis (a consistent +// re-survey -- the overlap case of crossing survey lines). +TEST(ImportEviction, ResurveyedBackscatterBlendsLosslessly) +{ + constexpr double kLat = 43.0; + constexpr double kLon = -70.0; + constexpr float kDepth = 12.0f; + constexpr float kI1 = 20.0f; // visit-1 intensity + constexpr float kI2 = 60.0f; // visit-2 intensity (same cell, same depth) + + std::vector> batches; + batches.push_back(surveyCell(kLat, kLon, kDepth, kI1)); // visit 1 + for (int i = 1; i <= 10; ++i) { // evict the tile + batches.push_back(surveyCell(43.0 + 0.02 * i, -71.0, 15.0f + i, 25.0f + i)); + } + batches.push_back(surveyCell(kLat, kLon, kDepth, kI2)); // visit 2, SAME cell + + const std::string root = makeTempDir("resurvey"); + const std::string b_dir = root + "/bounded"; + const std::string u_dir = root + "/unbounded"; + const std::string b_bs = root + "/bounded_bs"; + const std::string u_bs = root + "/unbounded_bs"; + runImport(batches, b_dir, b_bs, /*budget=*/3); // evict + reload-before-add + runImport(batches, u_dir, u_bs, /*budget=*/0); // never evicts + + // The resurveyed cell's CellIndex (deterministic for a fixed position). + marine_mbes_backscatter_store::MbesBackscatterStore probe = + marine_mbes_backscatter_store::MbesBackscatterStore::fromCellSize(kCellSize); + const gggs::CellIndex xcell = probe.cellIndex(kLat, kLon); + + const auto bs_b = loadBackscatterCells(b_bs); + const auto bs_u = loadBackscatterCells(u_bs); + ASSERT_TRUE(bs_b.count(xcell)) << "bounded build lost the resurveyed cell"; + ASSERT_TRUE(bs_u.count(xcell)) << "unbounded build lost the resurveyed cell"; + + const float bounded_i = bs_b.at(xcell).first; + const float unbounded_i = bs_u.at(xcell).first; + + // Bounded == unbounded == the blend of both equal-count visits ~ (I1 + I2) / 2. + EXPECT_NEAR(unbounded_i, 0.5f * (kI1 + kI2), 0.5f) << "unbounded should blend"; + EXPECT_NEAR(bounded_i, unbounded_i, 1e-3f) + << "bounded build must reproduce the unbounded blend (lossless backscatter)"; + // And the intensity VARIANCE (estimate variance, shrinks with sample count) + // matches too -- same sample population => same variance. + EXPECT_NEAR(bs_b.at(xcell).second, bs_u.at(xcell).second, 1e-3f) + << "bounded build must reproduce the unbounded intensity variance"; +} + +// The scratch spill directory is created during eviction and DELETED by finalize. +TEST(ImportEviction, ScratchDirCleanedUpAfterFinalize) +{ + const std::string root = makeTempDir("scratch"); + const std::string store_dir = root + "/store"; + const std::string bs_dir = root + "/bs"; + + GeoMapSheet sheet(kCellSize); + ImportAccumulator accumulator(sheet, makeConfig(store_dir, bs_dir, /*budget=*/3)); + for (int i = 0; i < 12; ++i) { // >> budget -> forces eviction (and a spill) + accumulator.addBatch(surveyCell(43.0 + 0.02 * i, -70.0, 10.0f + i, 30.0f + i)); + } + const std::string scratch = accumulator.scratchDir(); + ASSERT_FALSE(scratch.empty()) << "eviction should have created a scratch dir"; + EXPECT_TRUE(std::filesystem::exists(scratch)) << "scratch dir should exist mid-run"; + + marine_bathymetry_store::SourceRegistry bathy_reg; + marine_mbes_backscatter_store::SourceRegistry bs_reg; + accumulator.finalize(bathy_reg, &bs_reg); + + EXPECT_FALSE(std::filesystem::exists(scratch)) + << "finalize must delete the scratch spill dir"; + EXPECT_TRUE(accumulator.scratchDir().empty()) + << "finalize must clear the scratch path"; + + std::filesystem::remove_all(root); +} + +// Lossless guarantee: if a cold tile's persist FAILS, it must NOT be dropped +// (that would lose unsaved data). Force the failure by planting a regular file +// where the layer directory must be created, then survey past the budget. +TEST(ImportEviction, NeverDropsTileWhenPersistFails) +{ + const std::string root = makeTempDir("persistfail"); + const std::string store_dir = root + "/store"; + std::filesystem::create_directories(store_dir); + // layerDirName(Processed) == "processed": plant a FILE there so + // create_directories() throws and persistBathyTile cannot write. + { + std::ofstream blocker( + store_dir + "/" + + marine_bathymetry_store::layerDirName( + marine_bathymetry_store::SourceLayer::Processed)); + blocker << "not a directory"; + } + + GeoMapSheet sheet(kCellSize); + const std::size_t budget = 2; + ImportAccumulator accumulator(sheet, makeConfig(store_dir, "", budget)); + + const int n_tiles = 6; // >> budget + for (int i = 0; i < n_tiles; ++i) { + accumulator.addBatch(surveyCell(43.0 + 0.02 * i, -70.0, 10.0f + i, 30.0f)); + } + + // No bathy tile could be written (the layer dir can't be created), so no + // DATA-bearing tile was dropped: all surveyed (finite-data) tiles remain + // resident and RAM grew past the budget rather than losing unsaved data. + // (An empty tile a survey spilled into across a GGGS seam carries no data; the + // accumulator may drop such an empty tile harmlessly -- losing nothing -- so we + // assert on the data-bearing invariant, not an exact resident/evicted count.) + EXPECT_EQ(accumulator.bathyTilesPersisted(), static_cast(0)); + EXPECT_GE(accumulator.residentTileCount(), static_cast(n_tiles)); + EXPECT_GT(accumulator.residentTileCount(), budget); + + std::filesystem::remove_all(root); +} + +} // namespace cube diff --git a/cube_bathymetry/test/test_node.cpp b/cube_bathymetry/test/test_node.cpp index 6d61fe6..c5489e0 100644 --- a/cube_bathymetry/test/test_node.cpp +++ b/cube_bathymetry/test/test_node.cpp @@ -22,6 +22,8 @@ #include #include #include +#include +#include "cube_bathymetry/angular_response_curve.h" #include "cube_bathymetry/node.h" namespace cube @@ -439,13 +441,14 @@ TEST_F(NodeTest, ExclusionOnInterventionLeavesOriginalIntensityUnchanged) ASSERT_NE(deep, nullptr); ASSERT_NE(shallow, deep); - // Original (shallow) hypothesis keeps ONLY its own beam's intensity. - ASSERT_EQ(shallow->intensity_samples.size(), 1u); - EXPECT_FLOAT_EQ(shallow->intensity_samples[0].raw_intensity, -30.0f); + // Original (shallow) hypothesis keeps ONLY its own beam's intensity (Welford + // of one corrected sample == the raw, default None correction). + EXPECT_EQ(shallow->intensity.n, 1u); + EXPECT_DOUBLE_EQ(shallow->intensity.mean, -30.0); // The outlier's intensity is on the newly-seeded (deep) hypothesis. - ASSERT_EQ(deep->intensity_samples.size(), 1u); - EXPECT_FLOAT_EQ(deep->intensity_samples[0].raw_intensity, -10.0f); + EXPECT_EQ(deep->intensity.n, 1u); + EXPECT_DOUBLE_EQ(deep->intensity.mean, -10.0); } // extractNodeRecord() must honor the nominated_hypothesis_ priority path exactly @@ -462,8 +465,8 @@ TEST_F(NodeTest, NodeRecordHonorsNominatedHypothesis) // tell which path extractNodeRecord() took. auto nominated = std::make_shared(42.0f, 1.0f); nominated->input_sample_variance = 4.0f; - nominated->recordBeam(-12.0f, 0.0f); - nominated->recordBeam(-16.0f, 0.0f); + nominated->recordBeam(-12.0f, 0.0f, std::nan(""), params); // None -> corrected = raw + nominated->recordBeam(-16.0f, 0.0f, std::nan(""), params); NodeNominationTestAccess::nominate(n, nominated); auto record = n.extractNodeRecord(params); @@ -589,13 +592,14 @@ Parameters empiricalParams() } // Build a single-beam node via a nominated hypothesis so the surfaced intensity -// equals the corrected value of exactly that beam (no averaging). +// equals the corrected value of exactly that beam (no averaging). The correction +// is applied AT RECORD (cube#93), so the curve params are passed to recordBeam. std::shared_ptr singleBeamNominated( - Node & n, float raw_intensity, float beam_angle_rad) + Node & n, float raw_intensity, float beam_angle_rad, const Parameters & record_params) { auto h = std::make_shared(10.0f, 1.0f); h->input_sample_variance = 1.0f; - h->recordBeam(raw_intensity, beam_angle_rad); + h->recordBeam(raw_intensity, beam_angle_rad, std::nan(""), record_params); NodeNominationTestAccess::nominate(n, h); return h; } @@ -606,7 +610,7 @@ TEST_F(NodeTest, ARAInterpolation) { Parameters ara = empiricalParams(); Node n; - singleBeamNominated(n, -30.0f, deg2rad(30.0f)); + singleBeamNominated(n, -30.0f, deg2rad(30.0f), ara); auto record = n.extractNodeRecord(ara); ASSERT_EQ(record.n_samples, 1u); @@ -618,7 +622,7 @@ TEST_F(NodeTest, ARANadirIdentity) { Parameters ara = empiricalParams(); Node n; - singleBeamNominated(n, -30.0f, 0.0f); + singleBeamNominated(n, -30.0f, 0.0f, ara); auto record = n.extractNodeRecord(ara); ASSERT_EQ(record.n_samples, 1u); @@ -633,7 +637,7 @@ TEST_F(NodeTest, ARABeyondMaxAngleClampsToEdge) // corrected = raw - curveRel = -30 - (-12) = -18. Parameters ara = empiricalParams(); Node n; - singleBeamNominated(n, -30.0f, deg2rad(70.0f)); + singleBeamNominated(n, -30.0f, deg2rad(70.0f), ara); auto record = n.extractNodeRecord(ara); ASSERT_EQ(record.n_samples, 1u); @@ -645,7 +649,7 @@ TEST_F(NodeTest, ARANaNAngleIdentity) { Parameters ara = empiricalParams(); Node n; - singleBeamNominated(n, -30.0f, std::nan("")); + singleBeamNominated(n, -30.0f, std::nan(""), ara); auto record = n.extractNodeRecord(ara); ASSERT_EQ(record.n_samples, 1u); @@ -656,7 +660,7 @@ TEST_F(NodeTest, ARANaNAngleIdentity) TEST_F(NodeTest, ARAOffIsIdentity) { Node n; - singleBeamNominated(n, -30.0f, deg2rad(30.0f)); + singleBeamNominated(n, -30.0f, deg2rad(30.0f), params); // default None at record auto record = n.extractNodeRecord(params); // default None ASSERT_EQ(record.n_samples, 1u); @@ -672,8 +676,8 @@ TEST_F(NodeTest, ARAPortStarboardSymmetry) Node n; auto h = std::make_shared(10.0f, 1.0f); h->input_sample_variance = 1.0f; - h->recordBeam(-30.0f, deg2rad(30.0f)); // starboard - h->recordBeam(-30.0f, deg2rad(-30.0f)); // port + h->recordBeam(-30.0f, deg2rad(30.0f), std::nan(""), ara); // starboard + h->recordBeam(-30.0f, deg2rad(-30.0f), std::nan(""), ara); // port NodeNominationTestAccess::nominate(n, h); auto record = n.extractNodeRecord(ara); @@ -700,13 +704,15 @@ Parameters tier2Params(float alpha) } // Single-beam nominated hypothesis carrying a per-beam slant range, so the -// surfaced intensity equals the tier-2 corrected value of exactly that beam. +// surfaced intensity equals the tier-2 corrected value of exactly that beam. The +// correction is applied AT RECORD (cube#93), so the tier-2 params are passed in. std::shared_ptr singleBeamNominatedRange( - Node & n, float raw_intensity, float beam_angle_rad, float range) + Node & n, float raw_intensity, float beam_angle_rad, float range, + const Parameters & record_params) { auto h = std::make_shared(10.0f, 1.0f); h->input_sample_variance = 1.0f; - h->recordBeam(raw_intensity, beam_angle_rad, range); + h->recordBeam(raw_intensity, beam_angle_rad, range, record_params); NodeNominationTestAccess::nominate(n, h); return h; } @@ -729,7 +735,7 @@ TEST_F(NodeTest, Tier2RemovesTLAndResidual) const float range = 50.0f; Parameters t2 = tier2Params(alpha); Node n; - singleBeamNominatedRange(n, -30.0f, deg2rad(30.0f), range); + singleBeamNominatedRange(n, -30.0f, deg2rad(30.0f), range, t2); auto record = n.extractNodeRecord(t2); ASSERT_EQ(record.n_samples, 1u); @@ -745,7 +751,7 @@ TEST_F(NodeTest, Tier2NadirRemovesTLOnly) const float range = 35.0f; Parameters t2 = tier2Params(alpha); Node n; - singleBeamNominatedRange(n, -40.0f, 0.0f, range); + singleBeamNominatedRange(n, -40.0f, 0.0f, range, t2); auto record = n.extractNodeRecord(t2); ASSERT_EQ(record.n_samples, 1u); @@ -759,7 +765,7 @@ TEST_F(NodeTest, Tier2NaNRangeSkipsTL) { Parameters t2 = tier2Params(0.05f); Node n; - singleBeamNominatedRange(n, -30.0f, deg2rad(30.0f), std::nan("")); + singleBeamNominatedRange(n, -30.0f, deg2rad(30.0f), std::nan(""), t2); auto record = n.extractNodeRecord(t2); ASSERT_EQ(record.n_samples, 1u); @@ -772,7 +778,7 @@ TEST_F(NodeTest, Tier2NonPositiveRangeSkipsTL) { Parameters t2 = tier2Params(0.05f); Node n; - singleBeamNominatedRange(n, -30.0f, deg2rad(30.0f), -5.0f); + singleBeamNominatedRange(n, -30.0f, deg2rad(30.0f), -5.0f, t2); auto record = n.extractNodeRecord(t2); ASSERT_EQ(record.n_samples, 1u); @@ -785,7 +791,7 @@ TEST_F(NodeTest, Tier1IgnoresRange) { Parameters ara = empiricalParams(); // tl_removed defaults false Node n; - singleBeamNominatedRange(n, -30.0f, deg2rad(30.0f), 50.0f); + singleBeamNominatedRange(n, -30.0f, deg2rad(30.0f), 50.0f, ara); auto record = n.extractNodeRecord(ara); ASSERT_EQ(record.n_samples, 1u); @@ -810,4 +816,49 @@ TEST_F(NodeTest, Tier2RangeThreadsThroughUpdate) EXPECT_NEAR(record.intensity, static_cast(expected), 1e-2); } +// THE load-bearing cube#93 proof: correct-at-record (streaming Welford of the +// corrected intensity) reproduces correct-at-extract (retain raw + correct each at +// extract + naive mean/variance) for a fixed curve. Feed N varied beams through +// recordBeam (which corrects + folds), extract the node-output mean/variance, and +// compare to a reference computed the OLD way (correctBeamIntensity per beam, then +// naive sum / sum-of-squares -- the exact formula extractNodeRecord used before). +TEST_F(NodeTest, CorrectAtRecordMatchesCorrectAtExtract) +{ + Parameters ara = empiricalParams(); // Empirical, curve {{0,0},{60,-12}} + + // A spread of beams across the swath (mix of angles + dB), all on one nominated + // hypothesis so the record surfaces exactly this set. + const std::vector> beams = { // {raw_dB, angle_rad} + {-30.0f, deg2rad(0.0f)}, {-28.5f, deg2rad(12.0f)}, {-33.0f, deg2rad(25.0f)}, + {-26.0f, deg2rad(40.0f)}, {-31.5f, deg2rad(55.0f)}, {-29.0f, deg2rad(70.0f)}, + {-35.0f, deg2rad(33.0f)}, {-27.0f, deg2rad(8.0f)}}; + + Node n; + auto h = std::make_shared(10.0f, 1.0f); + h->input_sample_variance = 1.0f; + for (const auto & b : beams) { + h->recordBeam(b.first, b.second, std::nan(""), ara); // correct-at-record + } + NodeNominationTestAccess::nominate(n, h); + const auto record = n.extractNodeRecord(ara); + + // Reference: the OLD retain-and-correct-at-extract computation. + double sum = 0.0; + double sum_sq = 0.0; + const auto nN = static_cast(beams.size()); + for (const auto & b : beams) { + const double c = correctBeamIntensity(b.first, b.second, std::nan(""), ara); + sum += c; + sum_sq += c * c; + } + const double ref_mean = sum / nN; + const double ref_sample_var = (sum_sq - sum * sum / nN) / (nN - 1.0); + const double ref_var_of_mean = ref_sample_var / nN; + + ASSERT_EQ(record.n_samples, beams.size()); + // Welford vs naive agree to ~13+ digits; assert tight equivalence. + EXPECT_NEAR(record.intensity, static_cast(ref_mean), 1e-5); + EXPECT_NEAR(record.intensity_var, static_cast(ref_var_of_mean), 1e-5); +} + } // namespace cube