diff --git a/.agent/work-plans/issue-96/plan.md b/.agent/work-plans/issue-96/plan.md new file mode 100644 index 0000000..d138d0f --- /dev/null +++ b/.agent/work-plans/issue-96/plan.md @@ -0,0 +1,189 @@ +# Plan: Adopt simplified cube stores: seed precedence, batch-regen, sufficient-statistic backscatter + +## Issue + +https://github.com/rolker/cube_bathymetry/issues/96 + +## Context + +uma#248 has landed (jazzy), introducing: +- Two-layer store: `SourceLayer::Survey` + `SourceLayer::Reference` (replacing + `Draft`/`Processed`/`Chart`). +- 3-band `MbesCell {mean, standard_error, sample_sd}` sufficient statistics + (replacing `{intensity, intensity_variance, timestamp, source_index}`). +- `layerDirName()` now returns `"survey"` / `"reference"`. + +The current `cube_bathymetry` code does not compile against the new library. +This issue adds three features on top of the migration: + +1. **Seed precedence** — lazy per-tile seeding when a tile is first touched + (`survey` → settled, `reference` → predicted-only, else blank). +2. **Batch-regen** — scatter/gather path: no eviction → bit-exact output. +3. **Backscatter 3-band write/reload** — encode/reconstruct the Welford triplet + via the new `MbesCell` format. + +Greenfield: no backwards-compat required. + +## Approach + +### Step 1 — SourceLayer migration (compile fix for uma#248) + +Update every occurrence of the old enum values and MbesCell fields: + +| Old | New | +|-----|-----| +| `SourceLayer::Processed` | `SourceLayer::Survey` | +| `SourceLayer::Draft` | `SourceLayer::Survey` | +| `SourceLayer::Chart` | `SourceLayer::Reference` | +| `MbesCell::intensity` | `MbesCell::mean` | +| `MbesCell::intensity_variance` | computed from `standard_error`/`sample_sd` | + +Remove `--bathy-layer draft\|processed` CLI flag (both collapsed to `Survey`). +Remove `timestamp_ns`/`source_index` from all `MbesCell` construction sites. + +### Step 2 — Backscatter 3-band write (`geoGridToBackscatterCells`) + +Convert `NodeRecord {intensity, intensity_var, n_samples}` → `MbesCell`: + +- **n = 1 sentinel** (`intensity_var` is NaN): `{mean=intensity, standard_error=0, sample_sd=0}`. +- **n ≥ 2**: `sample_sd = sqrt(intensity_var * n_samples)`; + `standard_error = kScale * sqrt(intensity_var)` where `kScale = 1.96f`. + +### Step 3 — Backscatter Welford reconstruction (seed from store) + +New helper `welfordFromCell(MbesCell) -> IntensityWelford`: + +- `sample_sd == 0 && isfinite(mean)` → `{n=1, mean=mean, m2=0}` (n=1 sentinel). +- Else: `SE = standard_error / kScale`; `n = round((sample_sd / SE)^2)`; + `m2 = sample_sd^2 * (n - 1)`. + +**Edge case (review-plan suggestion): `sample_sd==0` for n≥2** — a cell whose +n≥2 samples are all *identical* has `M2=0` → `sample_sd=0`, indistinguishable +from the n=1 sentinel on reload (collapses to n=1). For continuous corrected-dB +data this is astronomically rare (exact float equality across ≥2 beams), and the +only consequence is a slightly under-counted `n` on a subsequent re-survey blend +of that one cell — the mean is still exact. **Documented as an accepted +limitation in the ADR-0007 addendum** (not guarded in code: a guard would need a +separate "n but zero-variance" encoding, not worth a band for a non-occurring +case). `kScale` is a single shared constant (matches the store's bathy +`stddev_to_confidence_interval_scale` convention) so write/read round-trip is +self-consistent. + +Used in Step 4 when reloading backscatter from the survey layer on first tile touch. + +### Step 4 — Seed precedence (`ImportAccumulator`) + +Add `std::set seeded_` to `ImportAccumulator`. Add private +`seedNewTile(index)`. + +`ImportAccumulatorConfig` gains `std::string reference_store_dir` (replaces +`--prior`; the survey seed always comes from `store_dir`). + +`addBatch` now: for each tile the batch touches — + +1. If in `evicted_`: run existing `reloadEvictedTile` + spill restore. +2. Else if NOT in `seeded_`: run `seedNewTile`: + a. Load bathy tile from `store_dir/survey/`; if found: `primeFromTile(..., seed_settled=true)`, + then reconstruct and restore each cell's `IntensityWelford` via `welfordFromCell`. + b. Else load bathy tile from `reference_store_dir/reference/`; if found: + `primeFromTile(..., seed_settled=false)` (no backscatter seed for reference). + c. Add index to `seeded_`. + +Remove the upfront whole-sheet `loadIntoSheet` call from `import_bag_main.cpp` +`--prior` handling. Replace `--prior` with `--reference-store ` that sets +`accumulator_config.reference_store_dir`. + +### Step 5 — Batch-regen (`batch_regen_bag` binary) + +New `src/batch_regen_main.cpp` and `CMakeLists.txt` entry. + +**Scatter phase**: same TF + detection-reading loop as `import_bag_main.cpp`. +For each projected ping, instead of CUBE accumulation: call +`sheet_.gridIndexForSounding(s)` to determine tile, then serialize the +`GeoSounding` to a per-tile binary bucket file (`/__.bin`). + +**Gather phase**: iterate tile bucket files. For each: +- Create a fresh `GeoMapSheet`; add all soundings from the bucket in one pass + (no eviction possible — single tile, bounded RAM). +- Seed from survey layer first if a tile already exists (same `seedNewTile` logic). +- Write the tile once via `ImportAccumulator::finalize` equivalent. + +**Test**: `test_batch_regen.cpp` — assert bit-exact depth/backscatter output vs +a single-pass unbounded `ImportAccumulator` run over the same soundings. + +## Files to Change + +| File | Change | +|------|--------| +| `include/cube_bathymetry/store_import.h` | `ImportAccumulatorConfig`: remove `bathy_layer`/`bs_source_index`, add `reference_store_dir`; update doc comments | +| `src/store_import.cpp` | Steps 2–4: `geoGridToBackscatterCells`, `persistBackscatterTile`, `reloadEvictedTile`, `addBatch`, new `seedNewTile`/`welfordFromCell` | +| `src/cube_bathymetry_node.cpp` | **Live-node SourceLayer migration (compile-breaker): `Draft`→`Survey`.** The live boat node now writes the `survey` layer directly (draft+processed collapsed; reference is read-only prior). Bounded/approximate boat-side product; batch-regen overwrites it as authoritative. Operator-confirmed behavior; document in the ADR-0001 addendum. | +| `test/test_persistence.cpp` | SourceLayer/MbesCell migration (compile-breaker: uses removed `Draft`/`Chart`) | +| `test/test_tile_eviction_rss.cpp` | SourceLayer/MbesCell migration (compile-breaker: uses removed `Draft`/`Chart`) | +| `src/import_bag_main.cpp` | Remove `--bathy-layer`, replace `--prior` with `--reference-store`, remove upfront loadIntoSheet block | +| `test/test_store_import.cpp` | SourceLayer + MbesCell field migration | +| `test/test_import_eviction.cpp` | SourceLayer + MbesCell migration; add reference-layer seed test (asserts no sample count increment) | +| `test/test_batch_regen.cpp` (new) | Bit-exact batch-regen equivalence test | +| `src/batch_regen_main.cpp` (new) | Batch-regen scatter/gather binary | +| `CMakeLists.txt` | Add `batch_regen_bag` executable target | +| `docs/decisions/0001-tile-eviction-and-incremental-publish.md` | Addendum: two-rung seed precedence semantics + `seed_settled` contract; batch-regen as the exact rebuild path | +| `docs/decisions/0007-mbes-backscatter-store-addendum-phase-b-transition.md` | Addendum: 3-band MbesCell write/reconstruct formulas, n=1 sentinel, Welford round-trip guarantee | +| `README.md` (or equivalent) | Document seed precedence, batch-regen mode, backscatter fidelity notes | + +## Principles Self-Check + +| Principle | Consideration | +|---|---| +| Human control and transparency | `--help` + README cover seed precedence rungs, batch-regen, backscatter notes | +| Enforcement over documentation | Reference-layer seed explicitly asserts `seed_settled=false`; test enforces no sample count increment | +| Capture decisions | ADR-0001 addendum documents two-rung precedence; ADR-0007 addendum documents 3-band formula | +| A change includes its consequences | `build_massabesic_store.sh` must be updated to remove `--append` grep guard; follow-up ops PR | +| Only what's needed | Batch-regen is a new binary, not an overload of the existing accumulator | + +## ADR Compliance + +The issue body's **"ADR-0002 (store)" citation is unh_marine_autonomy's ADR-0002** +(the bathy-store ADR) — already amended by #248 on the store side. This CUBE-side +work does **not** touch it; all decisions here are documented in cube_bathymetry's +own ADRs (ADR-0001 + ADR-0007 addenda) per the operator's Issue-Review resolution. + +| ADR | Triggered | How addressed | +|---|---|---| +| cube_bathymetry ADR-0001 | Yes | Addendum: two-rung seed precedence generalizes the reload/eviction prime paths; live-node now targets the `survey` layer | +| cube_bathymetry ADR-0007 addendum | Yes | Addendum: 3-band MbesCell sufficient statistics, reconstruction formula, n=1 sentinel + the n≥2 identical-sample limitation | + +## Consequences + +| If we change… | Also update… | Included? | +|---|---|---| +| Remove `--bathy-layer` | `build_massabesic_store.sh` (unh_echoboats_project11 #352) | No — follow-up ops PR (#353 update) | +| `MbesCell` 3-band format | Any consumer that reads `intensity`/`intensity_variance` | In scope: `test_import_eviction.cpp`; node.cpp already uses the new Welford triplet | +| Replace `--prior` with lazy seeding | Existing users of `--prior` flag | Document in `--help`; ops script update is follow-up | + +## Open Questions + +- None — uma#248 has landed; layer names (`survey`/`reference`) are finalized; all design decisions are specified in the issue header resolutions. + +## Estimated Scope + +Single PR (all three features share the `store_import.{h,cpp}` tile-seeding infrastructure and the uma#248 migration is a prerequisite for all three). + +## Deviations during implementation + +- **Registry migration (unplanned but required).** uma#248 also replaced the + per-cell `SourceRegistry`/`SourceRecord` interning table with a store-level + `StoreMetadata` sidecar. The plan's Files-to-Change was silent on this, but the + old types no longer compile. `ImportAccumulator::finalize` now takes + `const StoreMetadata *` pointers; `import_bag` builds `StoreMetadata` from + `--platform`/`--sensor`/`--campaign`(→survey). `--source-id`/`--sensor-class` + were dropped (greenfield; no per-cell source id exists anymore). +- **Dead conversion params removed.** Because `BathyCell` (2-band) and `MbesCell` + (3-band) no longer carry `timestamp`/`source_index`, those args were removed from + `geoGridToTile`/`mapSheetToTiles`/`geoGridToBackscatterCells`/ + `mapSheetToBackscatterCells`, and `source_index`/`timestamp_ns` were removed from + `ImportAccumulatorConfig` (plan only named `bathy_layer`/`bs_source_index`). +- **Batch-regen is a library class, not just a binary.** `BatchRegen` + (`batch_regen.{h,cpp}`) lives in the `cube_bathymetry_store_import` library so + `test_batch_regen` can drive it directly (a bag-only binary is not unit-testable). + Added `ImportAccumulator::persistResidentTile` for the single-tile gather write. + `batch_regen_main.cpp` + the `batch_regen_bag` target are as planned. diff --git a/.agent/work-plans/issue-96/progress.md b/.agent/work-plans/issue-96/progress.md new file mode 100644 index 0000000..b1094ce --- /dev/null +++ b/.agent/work-plans/issue-96/progress.md @@ -0,0 +1,289 @@ +--- +issue: 96 +--- + +# Issue #96 — Adopt simplified cube stores: seed precedence, batch-regen, sufficient-statistic backscatter + +## Issue Review +**Status**: complete +**When**: 2026-07-01 00:00 +00:00 +**By**: Claude Code Agent (Claude Sonnet) + +**Issue**: #96 +**Comment**: (best-effort post follows this entry; not recorded inline) +**Scope verdict**: well-scoped + +### Scope Assessment + +The three components (seed precedence, batch-regen, backscatter bands) are tightly +interdependent — all three consume the new uma#248 store format and share the +tile-reload path. Grouping them is appropriate. The batch-regen scatter/gather path +is the most novel and complex piece; it is explicitly scoped by an "bit-exact" +acceptance criterion that makes validation concrete. + +**Hard dependency: unh_marine_autonomy#248 must land first.** The two named store +layers (`survey` + `reference`, finalized names per the issue header note) and the +3-band backscatter sufficient-statistic format are defined there. Implementation +cannot begin until uma#248 is merged. + +**Right repo?** Yes — cube_bathymetry project repo. All three parts are CUBE-side +processing changes. + +### Principle Alignment + +| Principle | Status | Notes | +|---|---|---| +| Human control and transparency | OK | `--help` + README coverage is an acceptance criterion; seed precedence semantics clearly specified with per-rung behavior | +| Enforcement over documentation | Watch | Acceptance criteria imply testable behaviors (bit-exact regen, Welford round-trip, n=1 sentinel) but test coverage for seed-precedence rung semantics (reference layer must NOT be accumulated as measured data) is not explicitly stated as a test requirement | +| Capture decisions, not just implementations | Action needed | Issue references "ADR-0002 (store) — addendum expected" but no ADR-0002 exists in `docs/decisions/`; the existing ADR-0007 addendum file notes that the canonical `0007-mbes-backscatter-store.md` has not yet been authored. Implementation must resolve which ADR documents are authored/addended (see Actions) | +| A change includes its consequences | OK | Acceptance criteria cover ops script (`build_massabesic_store.sh`), help/README, backscatter round-trip test | +| Only what's needed | OK | Explicitly greenfield; scoped to what uma#248 enables; no backwards-compatibility required | +| Improve incrementally | Watch | Scope is substantial (seed precedence + batch-regen + backscatter), but all three parts share the same store-interface change and cannot be split without partial uma#248 adoption — grouping is justified | +| Test what breaks | Watch | "bit-exact vs never-evicted pass" and "Welford round-trip incl. n=1 sentinel" are good regression targets; explicitly add a test asserting reference-layer seed does not increment accumulated sample counts | +| Workspace vs. project separation | OK | All changes in cube_bathymetry project repo | + +### ADR Applicability + +| ADR | Triggered | Notes | +|---|---|---| +| Workspace ADR-0001 (adopt ADRs) | Yes | Design decisions being made; addendums expected as part of this work | +| Workspace ADR-0002 (worktree isolation) | Yes | Feature worktree already exists | +| cube_bathymetry ADR-0001 (lossless reload + eviction) | Yes | Seed precedence generalizes the startup prime and revisit-reload paths in ADR-0001 §1/§2; addendum required to document the two-rung precedence semantics and the `seed_settled` flag distinction | +| cube_bathymetry ADR-0007 addendum (MBES backscatter) | Yes | The sufficient-statistic backscatter bands (mean, confidence-scaled standard_error, sample_sd) extend Phase B.2 (Welford spill/restore); the Welford reconstruction formula (`n=(sample_sd/standard_error)^2`, M2, n=1 sentinel) is a new design decision that needs ADR coverage | + +**Gap**: "ADR-0002 (store)" cited in the issue is unresolvable as written — no `0002-*.md` exists in `docs/decisions/`. This is likely the not-yet-authored canonical `docs/decisions/0007-mbes-backscatter-store.md` mentioned in the existing ADR-0007 addendum file, or a new addendum to ADR-0001. The plan phase must decide and document which files are created/updated. + +### Consequences (from workspace map) + +- New store-layer names (`survey` + `reference`) replace prior placeholder names (`cube` + `pre-existing`) throughout — ensure all code paths, docs, and scripts use the finalized names. +- The ops script `build_massabesic_store.sh` (unh_echoboats_project11 #352/PR#353) is an explicit acceptance criterion dependency; changes to the importer CLI (`--seed` flags or equivalent) must keep that script working. +- Downstream: unh_marine_autonomy#247 (sidescan) is parked behind this; no action needed here, but plan/review should note the interface this issue establishes is a dependency for that work. + +### Actions +- [ ] Resolve the "ADR-0002 (store)" reference: decide whether to author `docs/decisions/0002-store-layer-design.md` (new) or add an addendum to ADR-0001 and the canonical ADR-0007, and state that decision in the plan. +- [ ] Add an explicit test for reference-layer seed precedence semantics: assert that a tile seeded from the `reference` layer does not contribute to accumulated sample counts (i.e., `seed_settled=false` path is not counted as measured data). +- [ ] Confirm that all identifiers, CLI flags, and store-path logic use the finalized layer names `survey` + `reference` (not the old `cube` + `pre-existing`). +- [ ] Verify uma#248 has landed before beginning implementation; gate the plan on this dependency. + +## Plan Authored +**Status**: complete +**When**: 2026-07-01 00:00 +00:00 +**By**: Claude Code Agent (Claude Sonnet) + +**Plan**: `.agent/work-plans/issue-96/plan.md` at `795b888` +**Branch**: feature/issue-96 at `795b888` +**Phases**: single + +### Open questions +- [ ] No open questions — plan is review-plan-ready. + +## Plan Review +**Status**: complete +**When**: 2026-07-01 11:54 +00:00 +**By**: Claude Code Agent (Claude Opus) + +**Plan**: `.agent/work-plans/issue-96/plan.md` at `795b888` +**PR**: PR-less (`--issue 96`) +**Verdict**: approve-with-suggestions + +### Findings +- [ ] (must-fix) File table omits three compile-breaking files that use removed `SourceLayer::Draft`/`Chart`: `src/cube_bathymetry_node.cpp`, `test/test_persistence.cpp`, `test/test_tile_eviction_rss.cpp` — `plan.md:104-116` +- [ ] (must-fix) Live-node `Draft`→`Survey` layer switch is an unstated behavioral consequence; confirm/document the live node's target layer — `plan.md:29-42` +- [ ] (suggestion) `welfordFromCell` collapses `n→1` when `sample_sd==0` for `n≥2` (identical samples); guard or document in ADR-0007 addendum — `plan.md:54-60` +- [ ] (suggestion) Explicitly resolve the issue's "ADR-0002 (store)" citation as addressed via ADR-0001/0007 addenda — `plan.md:129-134` + +## Implementation +**Status**: complete +**When**: 2026-07-01 00:00 +00:00 +**By**: Claude Opus + +**Branch**: feature/issue-96 (5 new commits: `ae5f68c`..`e022eb5`, on top of plan `6b92cb1`) + +### What was implemented (all plan steps) + +1. **uma#248 migration (compile fix).** `SourceLayer` Draft/Processed/Chart → + Survey/Reference; `BathyCell` 2-band; `MbesCell` 3-band; per-cell + `SourceRegistry` → store-level `StoreMetadata`. Touched `store_import.{h,cpp}`, + `cube_bathymetry_node.cpp` (live node now writes the **`survey`** layer, + operator-confirmed), `import_bag_main.cpp`, and all four affected tests. +2. **Backscatter 3-band write** (`geoGridToBackscatterCells`): n=1 sentinel + + n≥2 formulas, `kBackscatterConfidenceScale = 1.96f` shared constant. +3. **Backscatter Welford reconstruction** (`welfordFromCell`): exact inverse; + n≥2 identical-sample limitation documented in the ADR-0007 addendum (not + code-guarded, as planned). +4. **Seed precedence** (`ImportAccumulator::seedNewTile` + `seeded_` set): survey + warm-start (settled + backscatter Welford restore) > reference gate + (predicted-only, not measured, no backscatter) > blank. `import_bag` gains + `--reference-store` (replacing `--prior`); the upfront whole-sheet load is gone. +5. **Batch-regen**: `BatchRegen` (`batch_regen.{h,cpp}`, in the store_import + library) + `batch_regen_bag` binary (`batch_regen_main.cpp`) + CMake target. + Scatter to per-tile disk buckets → gather each tile in one unbounded pass → + write once via new `ImportAccumulator::persistResidentTile`. Scratch dir cleaned + up in `finalize`. +6. **Tests**: `test_batch_regen` (byte-exact vs single-pass unbounded, multi-tile + + revisit + scratch cleanup); `test_import_eviction.ReferenceSeedDoesNotAddMeasuredData`; + `test_store_import` Welford round-trips (`WelfordFromCellInvertsEncode`, + `BackscatterWelfordRoundTripMultiSample`, n=1 sentinel in `BackscatterCellsMatchGridRecords`). +7. **Docs**: ADR-0001 addendum (two-rung seed precedence + `seed_settled` contract + + batch-regen exact path + live-node→survey); ADR-0007 addendum (3-band + write/reconstruct formulas + n=1 sentinel + n≥2 identical-sample limitation); + README offline-import/rebuild + seed-precedence + backscatter-fidelity section. + +### Build / test status — GREEN (built in-container) + +The store libs (`marine_bathymetry_store`, `marine_mbes_backscatter_store` and deps) +were built from `layers/main/core_ws` source in-container, then `cube_bathymetry` +built and tested against them: + +- `cube_bathymetry`: **build GREEN**, **tests GREEN** — 441 tests, 0 failures, 58 + skipped (skips are unrelated pre-existing); lint (uncrustify + cpplint) clean on + all new/changed files. `batch_regen`'s 3 bit-exact-match tests pass. + +### Deviations (also noted in plan.md § Deviations) + +- Registry migration to `StoreMetadata` was required for compilation but unlisted + in the plan; `finalize()` signature changed to `StoreMetadata *` pointers. +- Removed dead `timestamp_ns`/`source_index` params from the conversion functions + and `ImportAccumulatorConfig` (the new cells don't carry them). +- Dropped `--source-id`/`--sensor-class` from `import_bag` (greenfield; no per-cell + source id). `BatchRegen` added as a library class (not just a binary) so it is + unit-testable. + +### Host must finish + +- **Out of scope (follow-up ops PR, different repo):** `build_massabesic_store.sh` + in unh_echoboats_project11 must drop `--bathy-layer`/`--prior`/`--source-id`/ + `--sensor-class` and switch to `--reference-store`; it will otherwise error on the + removed flags. Provenance now flows through `--platform`/`--sensor`/`--campaign`. +- Push + open PR (local-first; not pushed by this sub-agent). + +## Local Review (Pre-Push) +**Status**: complete +**When**: 2026-07-01 15:30 +00:00 +**By**: Claude Code Agent (Claude Opus) +**Verdict**: changes-requested + +**Branch**: feature/issue-96 at `0a83ad0` +**Mode**: pre-push +**Depth**: Deep (reason: +2688/-354 across 18 files, 2 new ADRs, store-format migration + new binary) +**Must-fix**: 2 | **Suggestions**: 6 +**Round**: 1 | **Ship**: continue — a real correctness gap in the headline "bit-exact" contract warrants a fix + re-review + +### Findings +- [x] (must-fix) Batch-regen not bit-exact near seams: scatter routes each sounding to its ±1-cell window, but `GeoGrid::insert` deposits up to `max_radius=CONF_99PC·√(horizontal_error)` (multiple cells for realistic TPU); single-pass gets far-radius cross-sounding deposits the gather omits. Tests use he=0.1 (sub-cell) so they mask it — `cube_bathymetry/src/batch_regen.cpp:162` +- [x] (must-fix) `welfordFromCell` UB: unguarded `se=standard_error/kScale` division when `standard_error==0` & `sample_sd!=0` → `lround(inf)` + out-of-range uint32 cast; add `se==0→n=1` guard (cross-confirmed both lenses) — `cube_bathymetry/src/store_import.cpp:199` +- [x] (suggestion) Re-running batch-regen into a non-empty `-o` store silently blends (seedNewTile warm-starts output-store survey tiles) instead of exact rebuild; warn/guard/`--append` — `cube_bathymetry/src/batch_regen.cpp` finalize +- [x] (suggestion) `closeAllStreams` doesn't check ofstream close-time flush; a disk-full final flush silently truncates a bucket → silently-wrong tile — `cube_bathymetry/src/batch_regen.cpp:181` +- [x] (suggestion) `seedNewTile` swallows a survey-seed load failure then overwrites the tile from scratch, discarding existing coverage on a transient read error (asymmetric with reload path) — `cube_bathymetry/src/store_import.cpp:584` +- [x] (suggestion) `batch_regen_main` `std::stod`/`std::stoi` unguarded → `std::terminate` on a bad flag value (mirrors import_bag) — `cube_bathymetry/src/batch_regen_main.cpp:343` +- [x] (suggestion) `-r` resolution not validated `> 0` before driving GGGS level — `cube_bathymetry/src/batch_regen_main.cpp:343` +- [x] (suggestion) `batch_regen_main.cpp` duplicates ~800 lines of `import_bag_main.cpp` (BagReaders/speedAt/projection pipeline); extract shared helper — `cube_bathymetry/src/batch_regen_main.cpp` (deferred: extracting BagReaders + speedAt + the projection pipeline into a shared helper is a substantial structural refactor with its own review/test risk, out of proportion to this fix pass; tracked as a follow-up) +- [x] (watch, pre-existing/out-of-scope) `marine_tiled_raster_store::saveTile` writes direct to the final path (no temp-then-rename); the "atomic … never corrupts the surface" comments (from #92) overstate the durability #96's gather/eviction leans on — follow-up to make the write atomic or fix the comments (partially addressed: corrected the overstated cube-side comments in `store_import.{cpp,h}` and `cube_bathymetry_node.cpp` to say saveTile is a direct, flush/close-checked, NOT crash-atomic write; deferred: making the write atomic is a `marine_tiled_raster_store` change in another repo — out of scope for this worktree, tracked as a follow-up) + +## Implementation +**Status**: complete +**When**: 2026-07-01 16:11 +00:00 +**By**: Claude Code Agent (Claude Opus) + +**Branch**: feature/issue-96 at `30136ff` +**Addressed**: Local Review (Pre-Push) — 2026-07-01 15:30 +00:00, branch feature/issue-96 at `0a83ad0` +**Commits**: `d68476f` (store_import), `83fc027` (batch_regen), `30136ff` (batch_regen_main) + +### Actions +- [x] (must-fix) Batch-regen bit-exact near seams — `addBatch` now routes the WHOLE batch to every tile in the batch's one-cell-expanded bounds, mirroring `GeoMapSheet::addSoundings` exactly, so a near-seam sounding's far-radius (`max_radius=CONF_99PC·√(horizontal_error)`) cross-sounding deposits are no longer dropped; the gather's `GeoGrid::insert` trims the superset bucket to the exact deposit set. Adds `SeamCrossingExactMatch` (multi-cell radius across a real GGGS seam) — `cube_bathymetry/src/batch_regen.cpp`, `test/test_batch_regen.cpp` +- [x] (must-fix) `welfordFromCell` `se<=0`/non-finite guard → `n=1` sentinel, avoiding `lround(inf/NaN)` UB + out-of-range uint32 cast — `cube_bathymetry/src/store_import.cpp` +- [x] (suggestion) Batch-regen exact rebuild: gather sets `skip_survey_seed` so it never blends onto tiles it rebuilds, and warns when the output survey layer is non-empty — `cube_bathymetry/src/batch_regen.cpp`, `include/cube_bathymetry/store_import.h`, `src/store_import.cpp` +- [x] (suggestion) Bucket close-time flush now checked — `flushOpenStreams` before gather + a flush/check on LRU eviction; a disk-full truncation is a hard error, not a silently-wrong tile — `cube_bathymetry/src/batch_regen.cpp` +- [x] (suggestion) `seedNewTile` protective on a survey-seed READ error — returns false so the caller drops the tile like a failed reload, protecting the intact-but-unreadable on-disk surface — `cube_bathymetry/src/store_import.cpp` +- [x] (suggestion) `batch_regen_main` numeric parsing guarded — `std::stod`/`std::stoi` via helpers that catch and route to `usage()` (and reject trailing junk) instead of `std::terminate` — `cube_bathymetry/src/batch_regen_main.cpp` +- [x] (suggestion) `batch_regen_main` rejects non-positive `-r` resolution up front — `cube_bathymetry/src/batch_regen_main.cpp` +- [x] (suggestion) `batch_regen_main`/`import_bag_main` ~800-line duplication — `cube_bathymetry/src/batch_regen_main.cpp` (deferred: substantial structural refactor with its own review/test risk, out of proportion to this fix pass; tracked follow-up) +- [x] (watch) `saveTile` "atomic … never corrupts" comments overstate durability — corrected the cube-side comments to say the write is direct + flush/close-checked but NOT crash-atomic — `src/store_import.cpp`, `include/cube_bathymetry/store_import.h`, `src/cube_bathymetry_node.cpp` (deferred: making the write crash-atomic is a `marine_tiled_raster_store` change in another repo, out of scope for this worktree) + +### Build / test — GREEN +Store libs (`marine_bathymetry_store`, `marine_mbes_backscatter_store`, +`marine_tiled_raster_store`) built from this worktree's `core_ws` source, then +`cube_bathymetry` built + tested against them: +- **build GREEN** (only pre-existing warnings: `tmpnam` in `test_angular_response_curve`, + unused-var in `test_tile_eviction_rss`). +- **tests GREEN** — full suite **442 tests, 0 failures, 58 skipped** (skips pre-existing). + `BatchRegen.SeamCrossingExactMatch` passes (asserts ≥2 tiles spanned and bit-exact + vs single-pass at a ~3-cell deposit radius); `test_store_import` (Welford round-trips) + and `test_import_eviction` green. +- **lint** — uncrustify + cpplint clean on all changed files. + +### Deferred (follow-ups, not done here) +- Extract a shared `BagReaders`/`speedAt`/projection helper from + `batch_regen_main.cpp` + `import_bag_main.cpp` (finding #8). +- Make `marine_tiled_raster_store::saveTile` crash-atomic (temp-then-rename) so the + "atomic" durability the gather/eviction leans on is real (finding #9, other repo). + +### 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 96 --skill review-code + +## Local Review (Pre-Push) +**Status**: complete +**When**: 2026-07-01 16:25 +00:00 +**By**: Claude Code Agent (Claude Opus) +**Verdict**: changes-requested + +**Branch**: feature/issue-96 at `16b0de1` +**Mode**: pre-push +**Depth**: Deep (reason: store-format migration + new binary + 2 ADR addenda, +3002/-354) +**Must-fix**: 1 | **Suggestions**: 3 +**Round**: 2 | **Ship**: recommended — one mechanical read-back guard; must-fix fell 2→1, not a design question + +### Findings +- [x] (must-fix) Gather read-back swallows I/O errors, breaking the bit-exact contract on the error path: read-open failure on a known-created bucket is `continue`d and `in.read` can't tell `badbit` from EOF → silently missing/truncated tile; guard symmetric with the hard-throwing write path — `cube_bathymetry/src/batch_regen.cpp:274-281` +- [x] (suggestion) `closeAllStreams` discards close-time errors, safe only because `flushOpenStreams` runs first; bind that ordering (assert/comment or check-and-throw) — `cube_bathymetry/src/batch_regen.cpp:224-228` +- [x] (suggestion) `import_bag_main` still uses raw `std::stod`/`std::stoi`/`std::stoll` → `std::terminate` on bad CLI value; mirror the guarded parse helpers added to the sibling `batch_regen_main` — `cube_bathymetry/src/import_bag_main.cpp:351-382` +- [x] (suggestion) `topic_info` loop var can be `const &` (cppcheck); folds into the tracked ~800-line dedup follow-up — `cube_bathymetry/src/batch_regen_main.cpp:183` + +### Notes +- Round 1's two must-fixes (near-seam bit-exactness, `welfordFromCell` division UB) verified correctly fixed under adversarial tracing (Lens A + Lens B) with new regression tests. +- Static analysis: cpplint + uncrustify clean; cppcheck 2 trivial style notes (finding above + a `useStlAlgorithm` nit dropped as below threshold). +- Governance Watch: `build_massabesic_store.sh` (unh_echoboats_project11) flag migration is a real breaking consequence, documented as an out-of-scope follow-up ops PR. +- Base `origin/jazzy` reviewed from local ref (remote fetch failed on host-key verification — possibly stale). + +### Next step +Lifecycle: **Local Review** → `address-findings` (verdict changes-requested) → re-run `review-code` → push / open PR → `triage-reviews`. +Ship advisory is **recommended** after the one mechanical must-fix lands; the host decides. + +## Implementation +**Status**: complete +**When**: 2026-07-01 16:44 +00:00 +**By**: Claude Code Agent (Claude Opus) + +**Branch**: feature/issue-96 at `b9d9a94` +**Addressed**: Local Review (Pre-Push) — 2026-07-01 16:25 +00:00, branch feature/issue-96 at `16b0de1` +**Commits**: `23dfe42` (batch_regen), `3ab4035` (import_bag_main), `b9d9a94` (batch_regen_main) + +### Actions +- [x] (must-fix) Gather read-back now symmetric with the hard-throwing scatter/flush path: a reopen failure on a bucket this run created **throws** (was `continue`d); after the read loop, `in.bad()` (a real I/O fault, not EOF) **throws**, and a non-zero `in.gcount()` at loop exit (a truncated partial trailing record) **throws** — no error path can silently drop or truncate a tile from the "bit-exact" rebuild — `cube_bathymetry/src/batch_regen.cpp` +- [x] (suggestion) `closeAllStreams` swallow-close-errors invariant bound with a comment: it is safe only because the success path (`finalize`) runs the checked `flushOpenStreams()` first; the only unflushed caller is `~BatchRegen` (exception mid-scatter), where throwing during unwinding is wrong and dropping soon-to-be-deleted scratch is correct — `cube_bathymetry/src/batch_regen.cpp` +- [x] (suggestion) `import_bag_main` numeric flag parsing guarded: `parse_double`/`parse_int`/`parse_long` helpers (mirroring the sibling `batch_regen_main`) catch and route bad `-r`/`--max-resident-tiles`/`-l`/`--minimum-range`/`--maximum-range` values to `usage()` (rejecting trailing junk) instead of `std::terminate` — `cube_bathymetry/src/import_bag_main.cpp` +- [x] (suggestion) `batch_regen_main` bag-topic lookup loop var taken by `const &` (cppcheck) — `cube_bathymetry/src/batch_regen_main.cpp` + +### Build / test — GREEN +Full layer stack rebuilt from source in-container (underlay_ws → core_ws store libs +`marine_autonomy`/`marine_bathymetry_store`/`marine_mbes_backscatter_store`/`marine_tiled_raster_store`), +then `cube_bathymetry` built + tested against them: +- **build GREEN** (only pre-existing warnings: `tmpnam` in `test_angular_response_curve`, + unused-var in `test_tile_eviction_rss`). +- **tests GREEN** — full suite **442 tests, 0 failures, 58 skipped** (skips pre-existing); + `BatchRegen` gather/seam tests, `test_store_import` Welford round-trips, `test_import_eviction` green. +- **lint** — uncrustify + cpplint clean on all changed files. + +### Deferred (unchanged from prior round — not re-opened here) +- Extract a shared `BagReaders`/`speedAt`/projection helper from + `batch_regen_main.cpp` + `import_bag_main.cpp` (the ~800-line dedup follow-up). +- Make `marine_tiled_raster_store::saveTile` crash-atomic (temp-then-rename), other repo. + +### 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 96 --skill review-code diff --git a/README.md b/README.md index 5f790ff..7c0b5b0 100644 --- a/README.md +++ b/README.md @@ -22,6 +22,49 @@ lookup `map_frame ← header.frame_id` at the ping stamp. Consumers read fields name** (`cube_bathymetry_node` and `bag_to_geotiff` both use named PointCloud2 iterators), so the field order is not load-bearing. +## Offline store import & rebuild + +Two offline tools replay a detections bag through CUBE and write a +`marine_bathymetry_store` (and optional `marine_mbes_backscatter_store`). Both +share `bag_to_geotiff`'s projection pipeline and write the **`survey`** layer +(unh_marine_autonomy#248 collapsed the old `draft`/`processed`/`chart` layers into +`survey` + `reference`; the off-boat re-run is authoritative). + +| Tool | RAM | Output | Use when | +|---|---|---|---| +| `import_bag` | bounded by `--max-resident-tiles` (persist-then-drop eviction, #92) | lossless, but a tile evicted mid-disambiguation has a slightly re-derived depth **uncertainty** (depth value faithful) | streaming / very large surveys where RAM is the constraint | +| `batch_regen_bag` | bounded by one tile's soundings | **bit-exact** vs a whole-survey-in-RAM build (depth, uncertainty, and backscatter) | the authoritative off-boat product | + +`batch_regen` scatters each projected sounding to a per-tile bucket on disk, then +gathers each tile in a single unbounded pass (no eviction) — so no tile is ever +evicted mid-disambiguation. Its scratch scatter directory is cleaned up at the end. + +### Seed precedence (`--reference-store`) + +On the first touch of each tile, both tools seed it with a two-rung precedence: + +1. **survey** — a `survey/` tile already in the output store (a pre-existing store, + or a tile written earlier this run) is warm-started as measured CUBE data + (settled depth **and** its backscatter Welford, restored losslessly). +2. **reference** — else, if `--reference-store ` is given, a `reference/` + (prior/contour) tile primes the **predicted surface only**: it turns CUBE's + blunder-rejection gate on so false-deep detections are dropped, but is **never + settled** as measured data and seeds **no** backscatter. Only tiles at the survey + GGGS level gate. (Replaces the pre-#96 `--prior` flag, which loaded the whole + prior into RAM up front and defeated eviction.) +3. else **blank**. + +### Backscatter fidelity + +The co-estimated per-cell backscatter is stored as a 3-band Welford sufficient +statistic `{mean, standard_error, sample_sd}` (uma#248), so the estimate +reconstructs losslessly on an off-boat re-run: a single-beam cell is the `n = 1` +sentinel (`sample_sd = 0`, finite mean), and a multi-beam cell round-trips +`n`/mean/variance exactly. Values are **uncorrected** by default; pass +`--backscatter-correction empirical --backscatter-curve ` for the per-beam +angular-response (and tier-2 TL) correction. See +`docs/decisions/0007-mbes-backscatter-store-addendum-phase-b-transition.md`. + ## Data flow & pose sourcing (design: #31) The package faithfully ports Calder's CUBE *algorithm* but originally diverged diff --git a/cube_bathymetry/CMakeLists.txt b/cube_bathymetry/CMakeLists.txt index 9c932af..633ccbc 100644 --- a/cube_bathymetry/CMakeLists.txt +++ b/cube_bathymetry/CMakeLists.txt @@ -91,7 +91,10 @@ install(DIRECTORY include/ DESTINATION include/${PROJECT_NAME}) # store). Linked by the offline import_bag tool and its test, AND -- since the # draft-tile persistence work (#21) -- by the runtime cube_bathymetry_node, which # uses geoGridToTile / loadIntoSheet to save and warm-start draft tiles. -add_library(cube_bathymetry_store_import src/store_import.cpp) +add_library(cube_bathymetry_store_import + src/store_import.cpp + src/batch_regen.cpp # exact scatter/gather store rebuild (#96) +) target_link_libraries(cube_bathymetry_store_import cube_bathymetry marine_bathymetry_store::marine_bathymetry_store @@ -229,6 +232,26 @@ install(TARGETS import_bag DESTINATION lib/${PROJECT_NAME} ) +# Exact, bounded-RAM store rebuild: scatter soundings to per-tile buckets, gather +# each tile in one unbounded pass (no eviction -> bit-exact output) (#96). +add_executable(batch_regen_bag src/batch_regen_main.cpp) +target_link_libraries(batch_regen_bag + cube_bathymetry_store_import # BatchRegen + cube_bathymetry + marine_bathymetry_store + marine_mbes_backscatter_store::marine_mbes_backscatter_store # --bs-store output (#80) + ${geometry_msgs_TARGETS} + ${marine_acoustic_msgs_TARGETS} + ${nav_msgs_TARGETS} + rclcpp::rclcpp + rosbag2_cpp::rosbag2_cpp + rosbag2_transport::rosbag2_transport + tf2_geometry_msgs::tf2_geometry_msgs + tf2_ros::tf2_ros +) + +install(TARGETS batch_regen_bag + DESTINATION lib/${PROJECT_NAME} +) + if(BUILD_TESTING) find_package(ament_lint_auto REQUIRED) ament_lint_auto_find_test_dependencies() @@ -303,6 +326,15 @@ if(BUILD_TESTING) marine_bathymetry_store::marine_bathymetry_store marine_mbes_backscatter_store::marine_mbes_backscatter_store ) + + # Batch-regen exactness (#96): the scatter/gather rebuild produces a byte-exact + # bathy + backscatter store vs a single-pass unbounded ImportAccumulator. + ament_add_gtest(test_batch_regen test/test_batch_regen.cpp) + target_link_libraries(test_batch_regen + 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 0b09ed5..d715dce 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 @@ -206,3 +206,59 @@ dependency) rather than retained. - Reload model aligns with cube#15 (cross-epoch self-improvement). - Consumer of the deferred `~/tiles` contract: rolker/unh_marine_autonomy#86, rolker/unh_marine_autonomy#250. + +## Addendum (cube_bathymetry#96) — two-rung seed precedence, `seed_settled`, batch-regen, live-node → survey + +unh_marine_autonomy#248 simplified the store taxonomy (`draft`/`processed`/`chart` +→ `survey`/`reference`) and the backscatter cell (3-band Welford sufficient +statistic). This addendum records the CUBE-side decisions #96 layered on top; the +backscatter-format decisions live in the ADR-0007 addendum. + +### Live node writes the `survey` layer + +The live `cube_bathymetry_node` now persists directly to the **`survey`** layer +(the old `draft` layer is gone; #248 collapsed draft+processed into one). The +boat-side product is bounded/approximate (the eviction-uncertainty artifact noted +in § Negative trade-offs); the off-boat **batch-regen** rebuild overwrites it as +the authoritative surface. The `draft_dir` ROS parameter name is retained (external +interface stability) but its tiles land under `/survey/`. +Operator-confirmed. + +### Two-rung seed precedence + the `seed_settled` contract + +The startup-prime (#21) and revisit-reload (#70) paths generalize into a single +**per-tile, first-touch seed precedence** in `ImportAccumulator::seedNewTile` +(reused by the batch-regen gather). When a tile is first touched: + +1. **survey** — a `survey/` bathy tile already on disk (a pre-existing store, or a + tile written earlier this run) is restored with `seed_settled=true`: settled + depth as a CUBE hypothesis (round-trips through `values()`, refines under new + soundings) **and** its per-cell backscatter Welford reconstructed from the + `survey/` backscatter tile (`welfordFromCell`, ADR-0007 addendum). Measured data. +2. **reference** — else a `reference/` prior tile is primed with + `seed_settled=false`: predicted-surface only (turns the blunder-rejection gate + on) but **never settled**, so it produces no `values()` output, seeds no + backscatter, and is **not counted as measured data**. Enforced by + `test_import_eviction.ReferenceSeedDoesNotAddMeasuredData`. +3. else **blank**. + +`seed_settled` is thus the contract boundary between *measured* (survey: settled + +backscatter) and *prior-only* (reference: gate-only). This replaces the pre-#96 +upfront whole-store `loadIntoSheet` of a `--prior` chart, which loaded the entire +prior into the sheet at once and defeated bounded-RAM eviction; the importer flag is +now `--reference-store` and priming is lazy, per tile. + +### Batch-regen — the exact rebuild path + +The eviction reload restores a *single reseeded hypothesis*, so a tile evicted +mid-disambiguation has a faithful depth **value** but a slightly re-derived depth +**uncertainty** (the § Negative trade-off above). **batch-regen** (`batch_regen_bag` +/ `BatchRegen`) removes even that artifact: it scatters every sounding to a per-tile +on-disk bucket (the sounding's one-cell-expanded window, matching `addSoundings`' +spread), then gathers each tile in a single **unbounded** pass over the complete set +of soundings that touch it — so no tile is ever evicted mid-disambiguation and the +output is **bit-exact** vs a whole-survey-in-RAM build (depth, uncertainty, and the +3-band backscatter). RAM is bounded by one tile's soundings, not surveyed area. +This is the authoritative off-boat product; the bounded-RAM `import_bag` remains the +live/streaming path. Enforced by `test_batch_regen` (byte-exact vs a single-pass +unbounded `ImportAccumulator`). 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 6e6bfe6..0738ed6 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 @@ -206,3 +206,60 @@ in the Context, Consequences, and Tier-2 sections above, and the backscatter correction assumes dB (a subtraction). Sonars reporting linear intensity would need a divide, or a domain conversion at ingest. **Follow-up to file:** model the intensity domain as a sonar property rather than assuming dB. + +## Addendum (cube_bathymetry#96) — 3-band `MbesCell` write/reconstruct + +unh_marine_autonomy#248 replaced the per-cell backscatter record +`{intensity, intensity_variance, timestamp, source_index}` with a **3-band +`Float32` Welford sufficient statistic** `{mean, standard_error, sample_sd}` (the +per-cell `timestamp`/`source_index` companions were dropped; coarse provenance is +now store-level `StoreMetadata`). The single layer is `survey` (draft/processed +collapsed). This addendum records how `cube_bathymetry` writes and reconstructs it; +the store round-trips the three floats bit-exactly and does no scaling itself. + +### Write (`geoGridToBackscatterCells`) + +From each node's `NodeRecord {intensity (running mean), intensity_var (estimate +variance = M2/(n−1)/n, NaN with <2 samples), n_samples}`: + +- **`mean = intensity`** always. +- **n = 1 sentinel** (`intensity_var` is NaN — a single beam, no dispersion): + `standard_error = 0`, `sample_sd = 0`. A finite `mean` with `sample_sd == 0` + reconstructs to `n = 1`; distinct from no-data (`mean = NaN`). +- **n ≥ 2**: `sample_sd = sqrt(intensity_var · n_samples)` (the beams' sample + standard deviation, `sqrt(M2/(n−1))`); `standard_error = kScale · + sqrt(intensity_var)` (the confidence-scaled standard error of the mean, true + `SE = sqrt(intensity_var) = sample_sd/sqrt(n)`). + +`kScale = kBackscatterConfidenceScale = 1.96f`, a single shared constant on both the +write and reconstruct sides (mirrors the bathy store's +`stddev_to_confidence_interval_scale` convention), so the round-trip is +self-consistent regardless of its numeric value. + +### Reconstruct (`welfordFromCell`) — the exact inverse + +Used to seed backscatter accumulation from a persisted `survey` tile on first tile +touch (seed precedence, ADR-0001 addendum) so an off-boat re-run blends with the +stored population rather than restarting it: + +- `!hasData()` (`mean` NaN) → empty Welford `{n = 0}`. +- `sample_sd == 0 && isfinite(mean)` → `{n = 1, mean, M2 = 0}` (n=1 sentinel). +- else `SE = standard_error / kScale`; `n = round((sample_sd / SE)²)`; + `M2 = sample_sd² · (n − 1)`. + +`round()` absorbs float round-trip error in `n`. Round-trip guarantee: for a node +with ≥2 samples of genuine dispersion the reconstructed `(n, mean, M2)` — hence the +estimate variance `(M2/(n−1))/n` — reproduces the original. Enforced by +`test_store_import.WelfordFromCellInvertsEncode` and +`BackscatterWelfordRoundTripMultiSample`. + +### Accepted limitation — n ≥ 2 identical samples + +A cell whose n ≥ 2 samples are all **exactly** identical has `M2 = 0` → +`sample_sd = 0`, indistinguishable on reload from the n=1 sentinel: it collapses to +`n = 1`. For continuous corrected-dB data this is astronomically rare (exact float +equality across ≥2 beams), and the only consequence is a slightly under-counted `n` +on a subsequent re-survey blend of that one cell — the **mean stays exact**. It is +**not code-guarded**: a guard would need a separate "n but zero-variance" encoding +(a fourth band, or a reserved sentinel), which is not worth it for a +non-occurring case. Documented here as an accepted limitation. diff --git a/cube_bathymetry/include/cube_bathymetry/batch_regen.h b/cube_bathymetry/include/cube_bathymetry/batch_regen.h new file mode 100644 index 0000000..1736a1d --- /dev/null +++ b/cube_bathymetry/include/cube_bathymetry/batch_regen.h @@ -0,0 +1,155 @@ +// 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. + + +#ifndef CUBE_BATHYMETRY__BATCH_REGEN_H_ +#define CUBE_BATHYMETRY__BATCH_REGEN_H_ + +#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/registry.hpp" +#include "marine_mbes_backscatter_store/registry.hpp" + +namespace cube +{ + +/// @brief Exact, bounded-RAM store rebuild via scatter/gather (cube_bathymetry#96). +/// +/// The authoritative off-boat rebuild path. Where @ref ImportAccumulator bounds RAM +/// by evicting cold tiles (which re-derives a revisited tile's depth UNCERTAINTY +/// from a single reseeded hypothesis — faithful depth, slightly drifted +/// uncertainty), batch-regen produces a **bit-exact** result: every tile is built +/// in a single pass over the COMPLETE set of soundings that touch it, so no tile is +/// ever evicted mid-disambiguation. +/// +/// **Scatter (streaming):** @ref addBatch appends each batch to the per-tile binary +/// bucket files on disk — one bucket per GGGS tile in the batch's one-cell-expanded +/// bounds, each receiving the WHOLE batch, exactly the tiles + soundings +/// @ref GeoMapSheet::addSoundings feeds each grid. Writing the whole batch (not just +/// the soundings whose own centre is near a tile) preserves the far-radius +/// cross-sounding deposits addSoundings makes near a seam (out to +/// max_radius = CONF_99PC·√(horizontal_error)), so the rebuild stays bit-exact for +/// any TPU. Nothing is accumulated in RAM during the scatter; resident memory is +/// bounded by the open bucket-stream cache, not by surveyed area. +/// +/// **Gather:** @ref finalize processes one bucket at a time. Each tile's soundings +/// are replayed, in global scatter order, into a FRESH @ref GeoMapSheet driven by a +/// single-tile @ref ImportAccumulator (unbounded, so it never evicts), which also +/// applies the two-rung seed precedence (survey warm-start / reference gate). Only +/// the target tile is written (@ref ImportAccumulator::persistResidentTile) — the +/// neighbour grids a seam sounding also created in the fresh sheet are each written +/// by their own bucket's gather. Because the tile's bucket holds exactly the +/// soundings a single unbounded pass would feed its grid, in the same order, the +/// output is bit-for-bit identical to that unbounded pass. +/// +/// The scatter scratch directory is deleted by @ref finalize (and by the destructor +/// on an exception). + class BatchRegen + { +public: + /// @brief Factory that builds a GeoMapSheet configured EXACTLY as the intended + /// single-pass run would (cell size, IHO order, backscatter correction). + /// + /// The gather builds one sheet per tile from this factory, and one long-lived + /// index sheet for tile routing, so every sheet shares identical CUBE parameters + /// — the precondition for bit-exact output. + using SheetFactory = std::function < std::unique_ptr < GeoMapSheet > () >; + + /// @param factory Builds a fully-configured, empty GeoMapSheet (see @ref SheetFactory). + /// @param config Store paths + seed precedence. `max_resident_tiles` is ignored + /// (the gather is single-tile and never evicts); the rest (store_dir, + /// reference_store_dir, bs_store_dir, cell_size_m) drive persistence and seeding. + BatchRegen(SheetFactory factory, ImportAccumulatorConfig config); + + /// @brief Delete the scratch scatter dir (RAII safety net for @ref finalize). + ~BatchRegen(); + + BatchRegen(const BatchRegen &) = delete; + BatchRegen & operator = (const BatchRegen &) = delete; + + /// @brief Scatter a batch of soundings to their per-tile bucket files. + void addBatch(const std::vector < GeoSounding > &soundings); + + /// @brief Gather every tile bucket, build + write each tile once, write the + /// store-level metadata, and delete the scratch scatter dir. + void finalize( + const marine_bathymetry_store::StoreMetadata * bathy_metadata = nullptr, + const marine_mbes_backscatter_store::StoreMetadata * bs_metadata = nullptr); + + /// @brief Cumulative bathy tile writes (valid after @ref finalize). + std::size_t bathyTilesPersisted() const {return bathy_persisted_;} + /// @brief Cumulative backscatter tile writes (valid after @ref finalize). + std::size_t backscatterTilesPersisted() const {return bs_persisted_;} + /// @brief Number of distinct tiles scattered to (bucket count). + std::size_t tileCount() const {return tiles_.size();} + /// @brief The scratch scatter directory (empty until the first scatter). Exposed + /// for tests that assert it is cleaned up after @ref finalize. + const std::string & scratchDir() const {return scratch_dir_;} + +private: + /// Return an append-open stream for @p index's bucket, lazily creating the + /// scratch dir and truncating the file the first time this tile is seen; keeps a + /// bounded LRU of open streams so the scatter never exhausts file descriptors. + std::ofstream & bucketStream(const gggs::GridIndex & index); + /// Flush every open bucket stream and throw if any flush failed (disk-full / I/O + /// error), so a truncated bucket is a hard error before the gather reads it rather + /// than a silently-wrong tile. Called by @ref finalize before @ref closeAllStreams. + void flushOpenStreams(); + void closeAllStreams(); + void cleanupScratch(); + std::string bucketPath(const gggs::GridIndex & index) const; + + SheetFactory factory_; + ImportAccumulatorConfig cfg_; + /// Tile routing only (no accumulation). + std::unique_ptr < GeoMapSheet > index_sheet_; + + /// Lazily created on first scatter; "" = none. + std::string scratch_dir_; + /// Every tile ever scattered to (its bucket file exists). + std::set < gggs::GridIndex > tiles_; + + // Bounded LRU of open append streams (front = most recently used). + std::list < gggs::GridIndex > lru_; + std::map < gggs::GridIndex, std::pair < std::ofstream, + std::list < gggs::GridIndex > ::iterator >> + open_streams_; + + std::size_t bathy_persisted_ = 0; + std::size_t bs_persisted_ = 0; + }; + +} // namespace cube + +#endif // CUBE_BATHYMETRY__BATCH_REGEN_H_ diff --git a/cube_bathymetry/include/cube_bathymetry/store_import.h b/cube_bathymetry/include/cube_bathymetry/store_import.h index f093e89..5f55f6e 100644 --- a/cube_bathymetry/include/cube_bathymetry/store_import.h +++ b/cube_bathymetry/include/cube_bathymetry/store_import.h @@ -33,6 +33,7 @@ #include "cube_bathymetry/geo_grid.h" #include "cube_bathymetry/geo_map_sheet.h" #include "cube_bathymetry/geo_sounding.h" +#include "cube_bathymetry/hypothesis.h" #include "marine_autonomy/gggs.h" #include "marine_bathymetry_store/bathymetry_store.hpp" #include "marine_bathymetry_store/bathymetry_tile.hpp" @@ -49,8 +50,9 @@ namespace cube /// The live `cube_bathymetry_node` accumulates into a geographic /// @ref GeoMapSheet (migrated from the Cartesian `MapSheet`, #21) so that /// `geoGridToTile` / `mapSheetToTiles` persist live data directly into the -/// `marine_bathymetry_store` `draft/` layer tiles via `tile_io` (atomic -/// temp-then-rename) with NO lossy Cartesian→geographic resample. The single +/// `marine_bathymetry_store` `draft/` layer tiles via `tile_io` (a direct, +/// flush/close-checked write — not crash-atomic) with NO lossy Cartesian→geographic +/// resample. The single /// fused `draft` grid (no per-day epochs, unh_marine_autonomy#221) accumulates /// newest-value-wins per cell. The costmap `bathymetry_layer` (#164) and the sim /// live loop (#77) then read exactly what CUBE writes. @@ -60,6 +62,16 @@ namespace cube /// queue, or pre-filter state (full Node deserialization is out of scope). CUBE /// continues accumulating from scratch on top of the seeded prediction surface. +/// @brief Confidence-interval scale applied to the backscatter standard error on +/// write and divided back out on reconstruction (uma#248 A.1). +/// +/// Mirrors the bathy store's `stddev_to_confidence_interval_scale` convention so +/// the stored `MbesCell::standard_error` is a confidence-scaled quantity, exactly +/// as the bathy uncertainty band is. A single shared constant on both the write +/// (@ref geoGridToBackscatterCells) and reconstruct (@ref welfordFromCell) sides +/// makes the Welford round-trip self-consistent regardless of its numeric value. + inline constexpr float kBackscatterConfidenceScale = 1.96f; + /// @brief Convert one CUBE @ref GeoGrid into a marine_bathymetry_store tile. /// /// `GeoGrid::values()` is a flat, positional `vector` in @@ -73,16 +85,14 @@ namespace cube /// /// `values()` returns `float`; the store cell is `double`, so each field is /// widened. `values()` mutates node state (it flushes the median pre-filter), so -/// it is called exactly **once** per grid here. +/// it is called exactly **once** per grid here. The pre-#248 per-cell +/// `timestamp`/`source_index` bands were dropped (`BathyCell` is now 2-band +/// `{depth, uncertainty}`); coarse provenance lives in the store-level +/// `StoreMetadata` sidecar instead. /// /// @param grid The CUBE grid to convert. -/// @param timestamp_ns Acquisition/import time written into every finite cell -/// (nanoseconds since the Unix epoch). A single value per -/// import keeps the result deterministic. -/// @param source_index Registry source index written into every finite cell. /// @return A `BathymetryTile` for `grid.index()` holding the finite cells. - marine_bathymetry_store::BathymetryTile geoGridToTile( - const GeoGrid & grid, int64_t timestamp_ns, uint16_t source_index); + marine_bathymetry_store::BathymetryTile geoGridToTile(const GeoGrid & grid); /// @brief Convert every grid of a @ref GeoMapSheet into a per-grid tile map. /// @@ -96,11 +106,8 @@ namespace cube /// `gggs::GridIndex` map order, and each grid converts deterministically. /// /// @param map_sheet The CUBE map sheet to convert. -/// @param timestamp_ns Acquisition/import time for every finite cell. -/// @param source_index Registry source index for every finite cell. std::map < gggs::GridIndex, marine_bathymetry_store::BathymetryTile > - mapSheetToTiles( - const GeoMapSheet & map_sheet, int64_t timestamp_ns, uint16_t source_index); + mapSheetToTiles(const GeoMapSheet & map_sheet); /// @brief Convert one CUBE @ref GeoGrid into its finite-backscatter cells (#80). /// @@ -115,35 +122,62 @@ namespace cube /// The returned map is keyed by `gggs::CellIndex` so the caller can write each /// cell with the store's public `set()` (no bulk import API is added to the /// separate `marine_mbes_backscatter_store` package; #80 stays within -/// `cube_bathymetry`). `intensity_var` (NaN with < 2 samples) maps to the cell's -/// `intensity_variance` quality band (ADR-0007 D6). `timestamp_ns`/`source_index` -/// are stamped into every emitted cell, exactly as the bathy path does — so the -/// Processed product carries provenance, not timestamp=0/source_index=0. +/// `cube_bathymetry`). /// -/// `nodeRecords()` flushes the median pre-filter, so it is called once per grid. +/// Each `NodeRecord {intensity (mean), intensity_var (variance of the mean), +/// n_samples}` is encoded into the 3-band `MbesCell {mean, standard_error, +/// sample_sd}` Welford sufficient statistic (uma#248 A.1): +/// - **n = 1 sentinel** (`intensity_var` is NaN, a single beam): `mean = +/// intensity`, `standard_error = 0`, `sample_sd = 0` (no dispersion; a finite +/// mean with `sample_sd == 0` reconstructs to `n = 1`). +/// - **n ≥ 2**: `sample_sd = sqrt(intensity_var * n_samples)` (the sample stddev +/// of the beams) and `standard_error = kBackscatterConfidenceScale * +/// sqrt(intensity_var)` (the confidence-scaled standard error of the mean). +/// +/// @ref welfordFromCell is the exact inverse, so the estimate round-trips +/// losslessly through the store on an off-boat re-run. `nodeRecords()` flushes the +/// median pre-filter, so it is called once per grid. /// /// @param grid The CUBE grid to convert. -/// @param timestamp_ns Acquisition/import time written into every emitted cell. -/// @param source_index Registry source index written into every emitted cell. /// @return A `gggs::CellIndex -> MbesCell` map of the finite-intensity cells. std::map < gggs::CellIndex, marine_mbes_backscatter_store::MbesCell > - geoGridToBackscatterCells( - const GeoGrid & grid, int64_t timestamp_ns, uint16_t source_index); + geoGridToBackscatterCells(const GeoGrid & grid); /// @brief Convert every grid of a @ref GeoMapSheet into one backscatter-cell map. /// /// Iterates `map_sheet.grids()` and merges each grid's /// @ref geoGridToBackscatterCells result. Grids cover disjoint GGGS cells, so the /// merge never collides. The caller writes the cells into a -/// `marine_mbes_backscatter_store::MbesBackscatterStore` via `set()` (Processed +/// `marine_mbes_backscatter_store::MbesBackscatterStore` via `set()` (Survey /// layer, #80). Deterministic for a fixed map sheet. /// /// @param map_sheet The CUBE map sheet to convert. -/// @param timestamp_ns Acquisition/import time for every emitted cell. -/// @param source_index Registry source index for every emitted cell. std::map < gggs::CellIndex, marine_mbes_backscatter_store::MbesCell > - mapSheetToBackscatterCells( - const GeoMapSheet & map_sheet, int64_t timestamp_ns, uint16_t source_index); + mapSheetToBackscatterCells(const GeoMapSheet & map_sheet); + +/// @brief Reconstruct a corrected-intensity @ref IntensityWelford from a stored +/// 3-band @ref MbesCell — the exact inverse of @ref geoGridToBackscatterCells. +/// +/// Used to seed backscatter accumulation from a persisted `survey` backscatter +/// tile on first tile touch (seed precedence, #96) so an off-boat re-run blends +/// with the stored population rather than restarting it. From the confidence- +/// scaled bands (dividing @ref kBackscatterConfidenceScale back out of +/// `standard_error` to recover the true `SE = sample_sd / sqrt(n)`): +/// - `sample_sd == 0 && isfinite(mean)` → `{n = 1, mean, m2 = 0}` (n=1 sentinel). +/// - else `SE = standard_error / kBackscatterConfidenceScale`; +/// `n = round((sample_sd / SE)^2)`; `m2 = sample_sd^2 * (n - 1)`. +/// +/// **Limitation (ADR-0007 addendum):** a cell whose n≥2 samples are all *exactly* +/// identical has `M2 = 0` → `sample_sd = 0`, indistinguishable from the n=1 +/// sentinel on reload (it collapses to `n = 1`). For continuous corrected-dB data +/// this is astronomically rare (exact float equality across ≥2 beams); the only +/// consequence is a slightly under-counted `n` on a later re-survey blend of that +/// one cell — the mean stays exact. Accepted, not code-guarded. +/// +/// @param cell A finite-mean (`hasData()`) backscatter cell. A no-data cell +/// (`mean` NaN) reconstructs to `{n = 0}` (empty Welford). + IntensityWelford welfordFromCell( + const marine_mbes_backscatter_store::MbesCell & cell); /// @brief Seed predicted depths in @p map_sheet from every finite cell of @p tile. /// @@ -163,10 +197,10 @@ namespace cube /// turns on `Node::insert`'s predicted-surface gate so a false-deep sounding /// below `target - blunder_scalar*sqrt(var)` is rejected) but does NOT /// fill/settle the cell -- no hypothesis, no `values()` output, the sheet stays -/// clean. This is the `Chart` (contour) prior path (cube#89): settling coarse -/// contour depths would contaminate the survey layer (and its co-estimated +/// clean. This is the `Reference` (contour/prior) path (cube#89, #96): settling +/// coarse prior depths would contaminate the survey layer (and its co-estimated /// backscatter) with non-measured fill, so the prior only gates, it does not -/// fill; survey-falls-through-to-chart gap-filling stays a query-time concern. +/// fill; survey-falls-through-to-reference gap-filling stays a query-time concern. void primeFromTile( const marine_bathymetry_store::BathymetryTile & tile, GeoMapSheet & map_sheet, bool seed_settled = true); @@ -182,40 +216,52 @@ namespace cube /// /// @param seed_settled Forwarded to @ref primeFromTile (default true = /// settled+predicted warm-start reload; false = predicted-only blunder-rejection -/// prior, e.g. priming the predicted surface from a `Chart` layer, cube#89). +/// prior, e.g. priming the predicted surface from a `Reference` layer, cube#89). void loadIntoSheet( const marine_bathymetry_store::BathymetryStore & store, marine_bathymetry_store::SourceLayer layer, GeoMapSheet & map_sheet, bool seed_settled = true); -/// @brief Configuration for @ref ImportAccumulator (cube_bathymetry#92). +/// @brief Configuration for @ref ImportAccumulator (cube_bathymetry#92, #96). 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). + /// final-resident bathy tiles are written here under the `survey/` layer + /// subdirectory; a revisited tile is reloaded from here, and a first-touched + /// tile is seeded from any pre-existing `survey/` tile (seed precedence #96). + /// Empty = no persistence (eviction is then disabled to stay lossless — there + /// is nowhere to drop to). The store always writes the `survey` layer (#248 + /// collapsed the draft/processed split; the off-boat CUBE re-run is authoritative). 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; + /// Optional reference-prior store directory (the importer's `--reference-store`, + /// replacing the pre-#96 `--prior`). When a tile is first touched and no `survey` + /// tile exists for it in @ref store_dir, a `reference/` tile here (if present) + /// seeds the CUBE predicted surface ONLY (blunder-rejection gate, seed_settled= + /// false): the coarse prior is never settled as measured data and never seeds + /// backscatter. Empty disables reference seeding. + std::string reference_store_dir; /// 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. + /// to reload an evicted tile, to seed a first-touched 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. + /// Empty disables backscatter co-persistence. Always writes the `survey` layer. 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; + /// When true, @ref ImportAccumulator::seedNewTile SKIPS the rung-1 survey + /// warm-start (it still applies the rung-2 reference gate). Set by @ref BatchRegen + /// for its gather: batch-regen replays a tile's COMPLETE sounding population in one + /// pass, so warm-starting that same tile from a pre-existing `survey/` tile in the + /// OUTPUT store would double-count it — the gather would blend onto data it is + /// about to fully reproduce, silently corrupting an exact rebuild. Forcing + /// from-scratch makes batch-regen a true rebuild regardless of the `-o` store's + /// prior contents. The live import path leaves this false (its warm-start is the + /// intended incremental-import behavior). + bool skip_survey_seed = false; }; /// @brief Bounded-RAM offline import accumulator (cube_bathymetry#92). @@ -279,11 +325,28 @@ namespace cube 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. + /// store-level metadata sidecars, and delete the scratch spill. Evicted + /// tiles are already durable. Call once at end-of-stream. + /// + /// @param bathy_metadata Optional store-level `registry.json` provenance for the + /// bathy store (uma#248 replaced the per-cell `SourceRegistry` interning table + /// with a single coarse `StoreMetadata` at the store root). Written only when + /// non-null and not `empty()`. + /// @param bs_metadata Optional store-level provenance for the backscatter store. void finalize( - const marine_bathymetry_store::SourceRegistry & bathy_registry, - const marine_mbes_backscatter_store::SourceRegistry * bs_registry = nullptr); + const marine_bathymetry_store::StoreMetadata * bathy_metadata = nullptr, + const marine_mbes_backscatter_store::StoreMetadata * bs_metadata = nullptr); + + /// @brief Persist ONE resident tile's bathy + backscatter, nothing else + /// (batch-regen gather, #96). + /// + /// Unlike @ref finalize this writes only @p index and no store-level metadata, + /// so the batch-regen gather can process one tile bucket in a fresh sheet and + /// write ONLY the target tile — never the neighbour grids a near-seam sounding + /// also created in that sheet. Runs the same @ref seedNewTile precedence via the + /// preceding @ref addBatch, so a survey/reference prior is honoured. Increments + /// the persisted-tile counters. + void persistResidentTile(const gggs::GridIndex & index); /// @brief Tiles currently resident in RAM. std::size_t residentTileCount() const; @@ -306,11 +369,33 @@ namespace cube void spillIntensitySamples(const gggs::GridIndex & index); void restoreSpilledSamples(const gggs::GridIndex & index); bool reloadEvictedTile(const gggs::GridIndex & index); + /// @brief Seed a tile the batch is touching for the FIRST time (seed precedence + /// #96), then mark it @ref seeded_. Two-rung precedence: + /// 1. survey — a `survey/` bathy tile in @ref store_dir (a pre-existing + /// store, or a resurvey of an already-written tile): prime settled + /// (`seed_settled=true`) AND restore each cell's backscatter Welford + /// via @ref welfordFromCell from the `survey/` backscatter tile. + /// 2. reference — else a `reference/` tile in @ref reference_store_dir: + /// prime predicted-only (`seed_settled=false`, blunder gate); NOT + /// counted as measured data, NO backscatter seed. + /// else blank (no prior). A no-op beyond marking @ref seeded_ when no + /// seed source is configured or found. + /// @return false if the rung-1 survey seed threw (the on-disk survey tile exists + /// but could not be read): the caller then drops this tile like a failed + /// reload, protecting the intact-but-unreadable surface from being + /// overwritten with from-scratch data, and leaves it UNseeded so a later + /// batch retries. true otherwise (seeded, or nothing to seed). + bool seedNewTile(const gggs::GridIndex & index); void cleanupScratch(); GeoMapSheet & sheet_; ImportAccumulatorConfig cfg_; std::set < gggs::GridIndex > evicted_; + /// Tiles already seeded (or confirmed blank) on first touch — seedNewTile runs + /// at most once per tile (seed precedence #96). Distinct from @ref evicted_: an + /// evicted tile was seeded, so it reloads (from its own written survey tile) + /// rather than re-seeding from the reference prior. + std::set < gggs::GridIndex > seeded_; std::string scratch_dir_; // lazily created on first eviction; "" = none std::size_t bathy_persisted_ = 0; std::size_t bs_persisted_ = 0; diff --git a/cube_bathymetry/src/batch_regen.cpp b/cube_bathymetry/src/batch_regen.cpp new file mode 100644 index 0000000..1de5060 --- /dev/null +++ b/cube_bathymetry/src/batch_regen.cpp @@ -0,0 +1,354 @@ +// 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. + + +#include "cube_bathymetry/batch_regen.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "marine_autonomy/gz4d_geo.h" +#include "marine_bathymetry_store/tile_io.hpp" + +namespace cube +{ + +namespace +{ +// Bound on simultaneously-open bucket streams so a wide survey never exhausts file +// descriptors. Scatter is spatially local (consecutive pings touch adjacent tiles), +// so a small LRU keeps nearly every write on an already-open stream. +constexpr std::size_t kMaxOpenStreams = 128; + +// One scattered sounding: the fields the CUBE insert path reads (geo_grid.cpp / +// node.cpp), so a reconstructed GeoSounding drives an identical estimate. Native +// layout — scratch-only, single-machine. `depth` is the position altitude (the +// GeoSounding ctor sets sounding.depth from point[2]). +#pragma pack(push, 1) +struct ScatterRecord +{ + double latitude; + double longitude; + double depth; + float vertical_error; + float horizontal_error; + float intensity; + float beam_angle; + float slant_range; + float predicted_depth_at_touchdown; +}; +#pragma pack(pop) + +ScatterRecord toRecord(const GeoSounding & s) +{ + ScatterRecord r; + r.latitude = s.latitude; + r.longitude = s.longitude; + r.depth = static_cast(s.sounding.depth); + r.vertical_error = s.sounding.vertical_error; + r.horizontal_error = s.sounding.horizontal_error; + r.intensity = s.sounding.intensity; + r.beam_angle = s.sounding.beam_angle; + r.slant_range = s.sounding.slant_range; + r.predicted_depth_at_touchdown = s.sounding.predicted_depth_at_touchdown; + return r; +} + +GeoSounding fromRecord(const ScatterRecord & r) +{ + GeoSounding s(gz4d::GeoPointLatLongDegrees(r.latitude, r.longitude, r.depth)); + s.sounding.vertical_error = r.vertical_error; + s.sounding.horizontal_error = r.horizontal_error; + s.sounding.intensity = r.intensity; + s.sounding.beam_angle = r.beam_angle; + s.sounding.slant_range = r.slant_range; + s.sounding.predicted_depth_at_touchdown = r.predicted_depth_at_touchdown; + return s; +} +} // namespace + +BatchRegen::BatchRegen(SheetFactory factory, ImportAccumulatorConfig config) +: factory_(std::move(factory)), cfg_(std::move(config)), index_sheet_(factory_()) +{ +} + +BatchRegen::~BatchRegen() +{ + // RAII safety net: finalize() normally cleans up, but an exception mid-scatter + // must not leak the scratch buckets. + closeAllStreams(); + cleanupScratch(); +} + +std::string BatchRegen::bucketPath(const gggs::GridIndex & index) const +{ + // Reuse the bathy tile naming for a unique per-tile stem. + return scratch_dir_ + "/" + marine_bathymetry_store::tileFilename(index) + ".bin"; +} + +std::ofstream & BatchRegen::bucketStream(const gggs::GridIndex & index) +{ + auto oit = open_streams_.find(index); + if (oit != open_streams_.end()) { + // Already open: move to the front of the LRU and reuse. + lru_.splice(lru_.begin(), lru_, oit->second.second); + return oit->second.first; + } + + // Lazily create the scratch dir on first scatter. mkdtemp atomically creates a + // uniquely-named directory (retries internally on collision), so two concurrent + // batch_regen processes can never share a scatter dir and cross-corrupt buckets. + if (scratch_dir_.empty()) { + std::string tmpl = + (std::filesystem::temp_directory_path() / "cube_regen_scatter_XXXXXX").string(); + std::vector buf(tmpl.begin(), tmpl.end()); + buf.push_back('\0'); + if (::mkdtemp(buf.data()) == nullptr) { + throw std::runtime_error( + std::string("batch_regen: cannot create scatter scratch dir from template ") + + tmpl + ": " + std::strerror(errno)); + } + scratch_dir_ = std::string(buf.data()); + } + + // First time we see this tile: truncate (fresh bucket). A later reopen after LRU + // eviction appends, so the bucket accumulates the full sounding sequence in order. + const bool first_time = tiles_.insert(index).second; + const std::ios::openmode mode = std::ios::binary | + (first_time ? std::ios::trunc : std::ios::app); + + // Evict the least-recently-used open stream(s) to stay under the FD budget. + while (open_streams_.size() >= kMaxOpenStreams && !lru_.empty()) { + const gggs::GridIndex victim = lru_.back(); + lru_.pop_back(); + auto vit = open_streams_.find(victim); + if (vit != open_streams_.end()) { + // Flush before the stream is closed+destroyed so a disk-full / I/O failure on + // the buffered data surfaces HERE. Letting the ofstream destructor swallow a + // close-time flush failure would truncate this bucket and yield a + // silently-wrong tile when the gather reads it back. + vit->second.first.flush(); + if (!vit->second.first) { + throw std::runtime_error( + "batch_regen: failed flushing scatter bucket " + bucketPath(victim) + + " on eviction (disk full?) -- aborting to avoid a silently-truncated tile"); + } + open_streams_.erase(vit); // closes the stream + } + } + + const std::string path = bucketPath(index); + std::ofstream out(path, mode); + if (!out) { + throw std::runtime_error("batch_regen: cannot open scatter bucket " + path); + } + lru_.push_front(index); + auto res = open_streams_.emplace( + index, std::make_pair(std::move(out), lru_.begin())); + return res.first->second.first; +} + +void BatchRegen::addBatch(const std::vector & soundings) +{ + if (soundings.empty()) { + return; + } + // Route the WHOLE batch into every tile addSoundings would create for it, so the + // scatter mirrors a single unbounded pass EXACTLY. GeoMapSheet::addSoundings takes + // gridIndicesForSoundings(batch) (the batch's one-cell-expanded bounds) and calls + // grid->insert(batch) on every grid in that set -- i.e. each such grid sees the + // ENTIRE batch, not just the soundings whose own centre is near it. A sounding + // near a tile seam deposits into a neighbour tile's cells out to + // max_radius = CONF_99PC*sqrt(horizontal_error), MULTIPLE cells for realistic TPU. + // Scattering each sounding only to its own one-cell window (the earlier approach) + // dropped those far-radius cross-sounding deposits, so the gather was not bit-exact + // near seams -- masked only because the tests use a sub-cell horizontal_error. + // Writing the whole batch to every grid in its expanded bounds reproduces + // addSoundings' touch set exactly; the gather's GeoGrid::insert re-applies the true + // per-cell radius test, so a bucketed sounding that does not actually reach the tile + // is harmlessly filtered (no false deposit) -- the bucket is a superset the radius + // test trims back to the exact single-pass deposit set. + for (const auto & idx : index_sheet_->gridIndicesForSoundings(soundings)) { + std::ofstream & out = bucketStream(idx); + for (const auto & s : soundings) { + const ScatterRecord rec = toRecord(s); + out.write(reinterpret_cast(&rec), sizeof(rec)); + if (!out) { + throw std::runtime_error("batch_regen: failed writing scatter bucket"); + } + } + } +} + +void BatchRegen::flushOpenStreams() +{ + // Force every still-open bucket's buffered data out and CHECK the result, so a + // disk-full / I/O error surfaces before the gather reads the buckets back rather + // than being swallowed by ofstream's destructor in closeAllStreams(). + for (auto & entry : open_streams_) { + std::ofstream & out = entry.second.first; + out.flush(); + if (!out) { + throw std::runtime_error( + "batch_regen: failed flushing scatter bucket " + bucketPath(entry.first) + + " (disk full?) -- aborting to avoid a silently-truncated tile"); + } + } +} + +void BatchRegen::closeAllStreams() +{ + // Deliberately swallows close-time flush errors (the ofstream destructor cannot + // report them). That is only safe because of a strict ordering invariant: on any + // SUCCESS path a checked flushOpenStreams() MUST run first (finalize() does), so every + // buffer is already empty and surfaced any disk-full/I/O fault. The only caller that + // reaches here with unflushed buffers is ~BatchRegen (an exception mid-scatter), where + // throwing would terminate during stack unwinding -- there, dropping the scratch data + // we were going to delete anyway is correct. Do NOT call this as the success-path close + // without a preceding flushOpenStreams(). + open_streams_.clear(); // destructor flushes + closes each ofstream + lru_.clear(); +} + +void BatchRegen::finalize( + const marine_bathymetry_store::StoreMetadata * bathy_metadata, + const marine_mbes_backscatter_store::StoreMetadata * bs_metadata) +{ + // Flush every bucket and CHECK the result before reading them back, so a disk-full + // truncation is a hard error here rather than a silently-wrong tile at gather; then + // close them all. + flushOpenStreams(); + closeAllStreams(); + + // A non-empty output survey layer means batch-regen is rebuilding over a populated + // store. The gather forces from-scratch (skip_survey_seed below) so it never blends + // onto the tiles it rebuilds, but tiles NOT touched by this run stay behind as + // stale survey data mixed with the fresh rebuild -- warn so the operator can point + // -o at an empty directory for a clean exact rebuild. + if (!cfg_.store_dir.empty()) { + const std::string survey_dir = cfg_.store_dir + "/" + + marine_bathymetry_store::layerDirName( + marine_bathymetry_store::SourceLayer::Survey); + std::error_code ec; + if (std::filesystem::is_directory(survey_dir, ec) && + !std::filesystem::is_empty(survey_dir, ec)) + { + std::cerr << "batch_regen: WARNING output survey layer '" << survey_dir + << "' is not empty; batch-regen rebuilds each touched tile from " + "scratch, but any pre-existing tile this run does NOT touch is left in " + "place (stale data mixed with the rebuild). Point -o at an empty directory " + "for a clean exact rebuild." << std::endl; + } + } + + // No eviction in the gather (each bucket is a single tile). Copy the config with + // the budget forced unbounded so the gather accumulator never drops a grid, and + // skip the rung-1 survey warm-start so the gather never double-counts a tile it is + // rebuilding from its complete bucket (exact rebuild, not a blend). + ImportAccumulatorConfig gather_cfg = cfg_; + gather_cfg.max_resident_tiles = 0; + gather_cfg.skip_survey_seed = true; + + // Gather one tile at a time, in deterministic (sorted GridIndex) order. + for (const auto & idx : tiles_) { + const std::string path = bucketPath(idx); + std::vector bucket; + { + // Every idx in tiles_ had its bucket created (and written) by bucketStream, so a + // reopen failure here is a real I/O fault, NOT a benign "missing bucket" -- treat + // it symmetrically with the hard-throwing scatter/flush path so it cannot silently + // drop a tile from the "bit-exact" rebuild. + std::ifstream in(path, std::ios::binary); + if (!in) { + throw std::runtime_error( + "batch_regen: cannot reopen scatter bucket " + path + + " for gather -- a bucket this run created is unreadable; aborting to avoid a " + "silently-missing tile"); + } + ScatterRecord rec; + while (in.read(reinterpret_cast(&rec), sizeof(rec))) { + bucket.push_back(fromRecord(rec)); + } + // read() sets eofbit+failbit at a clean end-of-file, but badbit on an actual I/O + // error -- distinguish them so a hardware read fault is not mistaken for EOF. A + // partial trailing record (gcount != 0 after the failed read) means the bucket was + // truncated to a non-record boundary; both are silently-wrong tiles, so throw. + if (in.bad()) { + throw std::runtime_error( + "batch_regen: I/O error reading scatter bucket " + path + + " during gather -- aborting to avoid a silently-truncated tile"); + } + if (in.gcount() != 0) { + throw std::runtime_error( + "batch_regen: scatter bucket " + path + + " ends with a partial record (truncated) -- aborting to avoid a " + "silently-wrong tile"); + } + } + if (bucket.empty()) { + continue; + } + // Fresh sheet + single-tile accumulator: replays the bucket in one pass (never + // evicts), applies the survey/reference seed precedence via addBatch, then + // writes ONLY the target tile (not the neighbour grids a seam sounding created). + std::unique_ptr sheet = factory_(); + ImportAccumulator acc(*sheet, gather_cfg); + acc.addBatch(bucket); + acc.persistResidentTile(idx); + bathy_persisted_ += acc.bathyTilesPersisted(); + bs_persisted_ += acc.backscatterTilesPersisted(); + } + + // Store-level provenance sidecars, written once (uma#248 StoreMetadata). + if (!cfg_.store_dir.empty() && bathy_metadata != nullptr && + !bathy_metadata->empty()) + { + std::filesystem::create_directories(cfg_.store_dir); + bathy_metadata->save(cfg_.store_dir); + } + if (!cfg_.bs_store_dir.empty() && bs_metadata != nullptr && + !bs_metadata->empty()) + { + std::filesystem::create_directories(cfg_.bs_store_dir); + bs_metadata->save(cfg_.bs_store_dir); + } + + // The scatter buckets were only needed to bucket-by-tile; the rebuild is done. + cleanupScratch(); +} + +void BatchRegen::cleanupScratch() +{ + if (!scratch_dir_.empty()) { + std::error_code ec; + std::filesystem::remove_all(scratch_dir_, ec); // best-effort; never throws here + scratch_dir_.clear(); + } +} + +} // namespace cube diff --git a/cube_bathymetry/src/batch_regen_main.cpp b/cube_bathymetry/src/batch_regen_main.cpp new file mode 100644 index 0000000..fb81889 --- /dev/null +++ b/cube_bathymetry/src/batch_regen_main.cpp @@ -0,0 +1,860 @@ +// 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. + +// batch_regen: offline detections-bag -> CUBE GeoMapSheet -> bathymetry-store +// EXACT rebuild (cube_bathymetry#96). Shares import_bag's projection pipeline but +// scatters each projected sounding to a per-tile bucket on disk, then gathers each +// tile in a single unbounded pass (no eviction) so the output is BIT-EXACT vs a +// whole-survey-in-RAM build -- the authoritative off-boat rebuild path. +// +// Mirrors the bag_to_geotiff `-d` offline-projection chain (rosbag2 +// SequentialReader -> tf2::BufferCore from /tf + /tf_static -> DetectionsProjector +// -> per-sounding lookupTransform("earth", frame_id, stamp) -> GeoSounding -> +// GeoMapSheet) but writes marine_bathymetry_store survey tiles instead of a GeoTIFF. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "cube_bathymetry/angular_response_curve.h" +#include "cube_bathymetry/batch_regen.h" +#include "cube_bathymetry/detections_projector.h" +#include "cube_bathymetry/geo_map_sheet.h" +#include "cube_bathymetry/geo_sounding.h" +#include "cube_bathymetry/store_import.h" +#include "geometry_msgs/msg/point_stamped.hpp" +#include "marine_acoustic_msgs/msg/sonar_detections.hpp" +#include "marine_autonomy/gz4d_geo.h" +#include "nav_msgs/msg/odometry.hpp" +#include "marine_bathymetry_store/bathymetry_store.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 "rclcpp/rclcpp.hpp" +#include "rosbag2_transport/reader_writer_factory.hpp" +#include "tf2_geometry_msgs/tf2_geometry_msgs.hpp" +#include "tf2_msgs/msg/tf_message.hpp" +#include "tf2/time.h" +#include "tf2_ros/buffer.h" + +[[noreturn]] void usage() +{ + std::cout << "usage: batch_regen [options] -o " + "-d [ ...]\n"; + std::cout << " Exact rebuild: scatters soundings to per-tile buckets on disk, " + "then gathers each tile in one unbounded pass (no eviction), so the output is " + "bit-exact vs a whole-survey-in-RAM build (cube#96). Use this for the " + "authoritative off-boat product; use import_bag for the bounded-RAM path.\n"; + std::cout << " -o : Output bathymetry-store directory (created if " + "needed). The off-boat CUBE re-run is authoritative, so it always writes the " + "`survey` layer (uma#248 collapsed the draft/processed split into one)\n"; + std::cout << " --reference-store : seed the CUBE predicted surface " + "lazily, per tile on first touch, from this store's `reference` (prior) layer " + "so blunder rejection drops false-deep detections (#89, #96). Only tiles at the " + "survey GGGS level gate. Predicted-only: the coarse prior is NEVER settled as " + "measured data and seeds no backscatter. NOTE: a coarse/shallow-biased prior " + "can also reject LEGITIMATE deeper-than-charted returns; the rejection margin " + "is tunable via the blunder_* params. (A `survey` tile already in -o takes " + "precedence and is warm-started as measured data instead.)\n"; + std::cout << " --bs-store : Also write an MBES backscatter store `survey` " + "layer from the same CUBE pass (optional; surfaces the co-estimated " + "intensity -- uncorrected by default, angle-corrected when " + "--backscatter-correction empirical is set, cube#81)\n"; + std::cout << " -d : marine_acoustic_msgs/SonarDetections " + "topic to replay through CUBE (required)\n"; + std::cout << " --odom-topic : nav_msgs/Odometry topic for per-ping " + "vessel speed-over-ground (optional; without it the speed-dependent " + "horizontal-TPU terms are floored to 0)\n"; + std::cout << " -r : Grid resolution (nominal; snapped to GGGS). " + "Default 1.0\n"; + std::cout << " --iho-order : CUBE IHO order (default order1a)\n"; + std::cout << " --backscatter-correction none|empirical: per-beam angular-response " + "correction at node-output (default none = identity). 'empirical' subtracts the " + "per-sonar curve from --backscatter-curve (cube#81)\n"; + std::cout << " --backscatter-curve : empirical angular-response curve CSV " + "(abs_angle_deg_center,mean_bs_db,n,db_relative_to_nadir). Required for " + "--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 << " -l : Stop after this many pings (debugging)\n"; + std::cout << " --platform / --sensor / --campaign : store-level provenance " + "written once to /registry.json (uma#248 StoreMetadata; --campaign maps " + "to the survey/campaign id). Per-cell source interning was retired for the " + "single-platform deployment\n"; + std::cout << " Frame/range overrides for the offline projector (must match " + "the bag's namespaced frames, or the grid comes out empty):\n"; + std::cout << " --base-link-frame, --level-frame, --tide-frame\n"; + std::cout << " --minimum-range, --maximum-range (meters)\n"; + exit(-1); +} + +// Speed-over-ground (m/s) nearest @p ns in the odometry timeline. Returns NaN +// when the timeline is empty OR the nearest sample is more than kSpeedMaxAgeNs +// away (an odom dropout) -- the error model then floors the speed-dependent +// horizontal-TPU terms to 0, the same as having no odometry at all, rather than +// applying an arbitrarily stale speed (cube_bathymetry#63 review). +float speedAt(const std::map & by_ns, int64_t ns) +{ + // Max staleness for a usable speed sample. Generous vs typical odom rates + // (tens of Hz) so normal jitter never trips it; trips only on a real dropout. + constexpr int64_t kSpeedMaxAgeNs = 5LL * 1000000000LL; + if (by_ns.empty()) { + return std::nanf(""); + } + auto it = by_ns.lower_bound(ns); + int64_t nearest_ns; + double nearest_speed; + if (it == by_ns.end()) { + nearest_ns = std::prev(it)->first; + nearest_speed = std::prev(it)->second; + } else if (it == by_ns.begin()) { + nearest_ns = it->first; + nearest_speed = it->second; + } else { + auto prev = std::prev(it); + const bool prev_closer = (ns - prev->first <= it->first - ns); + nearest_ns = prev_closer ? prev->first : it->first; + nearest_speed = prev_closer ? prev->second : it->second; + } + if (std::llabs(ns - nearest_ns) > kSpeedMaxAgeNs) { + return std::nanf(""); + } + return static_cast(nearest_speed); +} + +bool ends_with(const std::string & str, const std::string & suffix) +{ + if (str.length() >= suffix.length()) { + return 0 == str.compare(str.length() - suffix.length(), suffix.length(), suffix); + } + return false; +} + +/// Keep track of multiple bag files and retrieve messages in chronological order. +/// (Mirrors bag_to_geotiff's BagReaders.) +class BagReaders +{ +public: + /// A bag serialized message with the data type. + struct Message + { + std::string data_type; + rosbag2_storage::SerializedBagMessageSharedPtr message; + + using ConstPtr = std::shared_ptr; + + Message( + rosbag2_storage::SerializedBagMessageSharedPtr message, + const rosbag2_cpp::Reader & reader) + : message(message) + { + for (const auto & topic_info : reader.get_all_topics_and_types()) { + if (topic_info.name == message->topic_name) { + data_type = topic_info.type; + break; + } + } + } + }; + + explicit BagReaders(const std::vector & bagfile_names) + { + for (const auto & bagfile_name : bagfile_names) { + readers_[bagfile_name].open(bagfile_name); + } + } + + auto start_time() + { + auto start_time = readers_.begin()->second.reader->get_metadata().starting_time; + for (auto & reader : readers_) { + if (reader.second.reader->get_metadata().starting_time < start_time) { + start_time = reader.second.reader->get_metadata().starting_time; + } + } + return start_time; + } + + auto end_time() + { + auto end_time = readers_.begin()->second.reader->get_metadata().starting_time + + readers_.begin()->second.reader->get_metadata().duration; + for (auto & reader : readers_) { + if (reader.second.reader->get_metadata().starting_time + + reader.second.reader->get_metadata().duration > end_time) + { + end_time = reader.second.reader->get_metadata().starting_time + + reader.second.reader->get_metadata().duration; + } + } + return end_time; + } + + Message::ConstPtr next() + { + Bag * has_next = nullptr; + for (auto & reader : readers_) { + if (reader.second.has_next()) { + if (!has_next || + reader.second.peek_next()->message->send_timestamp < + has_next->peek_next()->message->send_timestamp) + { + has_next = &reader.second; + } + } + } + if (has_next) { + return has_next->pop_next(); + } + return {}; + } + +private: + struct Bag + { + std::unique_ptr reader; + Message::ConstPtr next_message; + + void open(const std::string & file_name) + { + rosbag2_storage::StorageOptions storage_options; + storage_options.uri = file_name; + reader = rosbag2_transport::ReaderWriterFactory::make_reader(storage_options); + reader->open(storage_options); + if (reader->has_next()) { + next_message = std::make_shared(reader->read_next(), *reader); + } + } + + bool has_next() + { + return static_cast(next_message); + } + + Message::ConstPtr peek_next() + { + return next_message; + } + + Message::ConstPtr pop_next() + { + auto return_value = next_message; + if (reader->has_next()) { + next_message = std::make_shared(reader->read_next(), *reader); + } else { + next_message.reset(); + } + return return_value; + } + }; + + std::map readers_; +}; + + +int main(int argc, char * argv[]) +{ + std::vector arguments(argv + 1, argv + argc); + if (arguments.empty()) { + usage(); + } + + std::vector bagfile_names; + std::string store_dir; + // optional: reference-prior store to seed the predicted surface lazily (#89, #96) + std::string reference_store_dir; + std::string bs_store_dir; // optional: MBES backscatter store output (#80) + std::string detections_topic; // required + std::string odom_topic; // optional: nav_msgs/Odometry for per-ping vessel speed + double resolution = 1.0; + std::string iho_order = "order1a"; + int ping_count_limit = 0; + // Backscatter angular-response correction (cube#81). Default none = identity. + std::string backscatter_correction_str = "none"; + std::string backscatter_curve_file; + + // Store-level provenance (uma#248 StoreMetadata, written once at finalize). + marine_bathymetry_store::StoreMetadata store_metadata; + + // DetectionsProjector configuration for the offline path. Defaults match + // detections_to_pointcloud's node defaults so a detections-only bag projects + // the same way the live node would. + cube::ProjectorParams projector_params; + + // Consume the value token following a value-taking flag, bounds-checked: a + // flag given as the last argument (no value after it) would otherwise advance + // `arg` past end() and dereference it (UB). usage() is [[noreturn]], so on the + // missing-value error this never falls through to the ++arg/deref. `arg` is + // declared before the lambda so the by-reference capture binds the loop cursor. + auto arg = arguments.begin(); + auto next_value = [&](const char * flag) -> const std::string & { + if (std::next(arg) == arguments.end()) { + std::cerr << "error: option '" << flag << "' requires a value\n"; + usage(); + } + ++arg; + return *arg; + }; + + // Guarded numeric parses: std::stod/std::stoi THROW on a non-numeric value and, + // uncaught in main, would std::terminate the process with an opaque message + // (mirrors the hardening import_bag still lacks). Catch and route to usage() with + // a clear diagnostic instead. usage() is [[noreturn]], so these never fall through. + auto parse_double = [&](const char * flag, const std::string & value) -> double { + try { + std::size_t consumed = 0; + const double parsed = std::stod(value, &consumed); + if (consumed != value.size()) { + throw std::invalid_argument("trailing characters"); + } + return parsed; + } catch (const std::exception &) { + std::cerr << "error: option '" << flag << "' expects a number, got '" + << value << "'\n"; + usage(); + } + }; + auto parse_int = [&](const char * flag, const std::string & value) -> int { + try { + std::size_t consumed = 0; + const int parsed = std::stoi(value, &consumed); + if (consumed != value.size()) { + throw std::invalid_argument("trailing characters"); + } + return parsed; + } catch (const std::exception &) { + std::cerr << "error: option '" << flag << "' expects an integer, got '" + << value << "'\n"; + usage(); + } + }; + + for (; arg != arguments.end(); arg++) { + if (*arg == "-h") { + usage(); + } else if (*arg == "-o") { + store_dir = next_value("-o"); + } else if (*arg == "--reference-store") { + reference_store_dir = next_value("--reference-store"); + } else if (*arg == "--bs-store") { + bs_store_dir = next_value("--bs-store"); + } else if (*arg == "-d") { + detections_topic = next_value("-d"); + } else if (*arg == "--odom-topic") { + odom_topic = next_value("--odom-topic"); + } else if (*arg == "-r") { + resolution = parse_double("-r", next_value("-r")); + } else if (*arg == "--iho-order") { + iho_order = next_value("--iho-order"); + } else if (*arg == "--backscatter-correction") { + backscatter_correction_str = next_value("--backscatter-correction"); + } else if (*arg == "--backscatter-curve") { + backscatter_curve_file = next_value("--backscatter-curve"); + } else if (*arg == "-l") { + ping_count_limit = parse_int("-l", next_value("-l")); + } else if (*arg == "--platform") { + store_metadata.platform = next_value("--platform"); + } else if (*arg == "--sensor") { + store_metadata.sensor = next_value("--sensor"); + } else if (*arg == "--campaign") { + store_metadata.survey = next_value("--campaign"); + } else if (*arg == "--base-link-frame") { + projector_params.base_link_frame = next_value("--base-link-frame"); + } else if (*arg == "--level-frame") { + projector_params.level_frame = next_value("--level-frame"); + } else if (*arg == "--tide-frame") { + projector_params.tide_frame = next_value("--tide-frame"); + } else if (*arg == "--minimum-range") { + projector_params.minimum_range = + parse_double("--minimum-range", next_value("--minimum-range")); + } else if (*arg == "--maximum-range") { + projector_params.maximum_range = + parse_double("--maximum-range", next_value("--maximum-range")); + } else if (!arg->empty() && (*arg)[0] == '-' && *arg != "-") { + // An unrecognized flag would otherwise be silently treated as a bag path + // and fail later with a confusing "cannot open bag". Reject it up front. + std::cerr << "error: unknown option '" << *arg << "'"; + if (*arg == "-e") { + std::cerr << " -- the -e argument was removed in " + "cube_bathymetry#69; the store no longer uses per-day epochs, so the " + "bag now imports into a single fused draft grid with no date label"; + } + std::cerr << "\n"; + usage(); + } else { + bagfile_names.push_back(*arg); + } + } + + if (store_dir.empty() || detections_topic.empty() || bagfile_names.empty()) { + std::cerr << "error: -o , -d , and " + "at least one bag are all required\n"; + usage(); + } + + // Resolution drives the GGGS level (Level::fromCellSize); a non-positive value is + // nonsensical and would produce a degenerate/garbage level rather than fail cleanly. + if (!(resolution > 0.0)) { + std::cerr << "error: -r resolution must be > 0 (got " << resolution << ")\n"; + usage(); + } + + std::cout << "Detections topic: " << detections_topic + << " (offline projection, vessel_speed = NaN)" << std::endl; + std::cout << "Store dir: " << store_dir << std::endl; + + cube::DetectionsProjector projector(projector_params); + + // Accumulated offline-projection diagnostics, surfaced at the end so a + // misconfigured-frames or over-tight-range run is diagnosable rather than a + // silently sparse/empty import (the failure mode #43 exists to kill). + size_t proj_pings = 0; + size_t proj_soundings = 0; + size_t proj_filtered_range = 0; + size_t proj_missing_attitude = 0; + size_t proj_missing_heave = 0; + size_t proj_dropped_georef = 0; // pings with no earth transform at their stamp + + BagReaders bag_readers(bagfile_names); + + std::cout << "calculating total time..." << std::endl; + auto begin_time = bag_readers.start_time(); + auto end_time = bag_readers.end_time(); + auto total_duration = end_time - begin_time; + + auto start_time_t = std::chrono::system_clock::to_time_t(begin_time); + std::cout << "start time: " + << std::put_time(std::gmtime(&start_time_t), "%Y-%m-%d %H:%M:%S") << std::endl; + auto end_time_t = std::chrono::system_clock::to_time_t(end_time); + std::cout << "end time: " + << std::put_time(std::gmtime(&end_time_t), "%Y-%m-%d %H:%M:%S") << std::endl; + std::cout << "total time: " + << std::chrono::duration_cast(total_duration).count() + << " seconds" << std::endl; + + // Single, deterministic per-cell timestamp: the bag's nominal start time in ns + // since the Unix epoch. A single value per import keeps the epoch deterministic + // (every replay of the same bag yields the same cell stamps). + const int64_t cell_timestamp_ns = + std::chrono::duration_cast( + begin_time.time_since_epoch()).count(); + + // Bounded TF cache: the offline projection only needs the transforms + // bracketing each ping's stamp, so a small rolling window keeps tf2's + // std::list-backed TimeCache short. Whole-bag buffering made every + // lookupTransform an O(n) linear list walk (~7 pings/s; cube_bathymetry#63). + // The guard is how far the TF frontier must advance past a ping before we + // project it, so both bracketing transforms are present (no frontier drops). + constexpr double kCacheWindowSec = 60.0; + constexpr double kGuardSec = 3.0; + const int64_t kGuardNs = static_cast(kGuardSec * 1e9); + + auto clock = std::make_shared(); + tf2_ros::Buffer tfBuffer(clock, tf2::durationFromSec(kCacheWindowSec)); + + cube::GeoMapSheet geo_map_sheet(resolution, iho_order); + std::cout << "requested resolution: " << resolution << " nominal used: " + << geo_map_sheet.nominalCellSizeMeters() << std::endl; + + // Backscatter angular-response correction (cube#81). The setter must run AFTER + // the sheet is constructed (its grids hold a const ref to the sheet Parameters). + cube::BackscatterAngleCorrection backscatter_mode = + cube::BackscatterAngleCorrection::None; + if (!cube::parseBackscatterAngleCorrection( + backscatter_correction_str, backscatter_mode)) + { + std::cerr << "error: --backscatter-correction must be 'none' or 'empirical' " + << "(got '" << backscatter_correction_str << "')\n"; + usage(); + } + cube::AngularResponseCurve backscatter_curve; + if (backscatter_mode == cube::BackscatterAngleCorrection::Empirical && + !backscatter_curve_file.empty()) + { + backscatter_curve = cube::loadAngularResponseCurveWithHeader(backscatter_curve_file); + } + if (backscatter_mode == cube::BackscatterAngleCorrection::Empirical && + backscatter_curve.points.empty()) + { + // Loud, not silent: enabled but no curve loaded -> correction is a no-op. + std::cerr << "warning: --backscatter-correction empirical but no curve was " + "loaded from --backscatter-curve '" << backscatter_curve_file + << "' -- the correction is ENABLED but a NO-OP (intensity emitted " + "uncorrected). Provide a valid curve CSV.\n"; + } else if (backscatter_mode == cube::BackscatterAngleCorrection::Empirical) { + std::cout << "Backscatter angular-response correction: empirical, " + << backscatter_curve.points.size() << "-point curve from " + << backscatter_curve_file; + if (backscatter_curve.tl_removed) { + // tier-2 (cube#87): the curve is a TL-removed residual; the estimator + // also removes 40*log10(R) + 2*alpha*R per beam. + std::cout << " [tier-2: TL-removed, alpha=" + << backscatter_curve.absorption_db_per_m << " dB/m]"; + } + std::cout << std::endl; + } + geo_map_sheet.setBackscatterCorrection( + backscatter_mode, backscatter_curve.points, + backscatter_curve.tl_removed, backscatter_curve.absorption_db_per_m); + + // Sheet factory (#96): batch-regen builds one routing sheet + one gather sheet + // per tile, all of which MUST be configured identically to this projection sheet + // (cell size, IHO order, backscatter correction) for the rebuild to be exact. + // Capture the correction settings by value so the factory can build many sheets. + cube::BatchRegen::SheetFactory make_sheet = + [resolution, iho_order, backscatter_mode, + curve_points = backscatter_curve.points, + tl_removed = backscatter_curve.tl_removed, + absorption = backscatter_curve.absorption_db_per_m]() { + auto sheet = std::make_unique(resolution, iho_order); + sheet->setBackscatterCorrection( + backscatter_mode, curve_points, tl_removed, absorption); + return sheet; + }; + + // Reference-prior seeding (#89, #96) is LAZY, per tile on first touch, driven by + // the gather accumulator's seedNewTile: a `reference/` tile primes the CUBE + // predicted surface only (seed_settled=false) so the blunder gate turns on WITHOUT + // settling coarse prior depths as measured data. Only tiles at the survey GGGS + // level coincide with a survey node and gate. + + // Store-level provenance (uma#248 StoreMetadata), written once at finalize. + // Backscatter provenance mirrors the bathy platform/sensor with an MBES-specific + // calibration ref (empty until a beam-pattern calibration exists). + marine_mbes_backscatter_store::StoreMetadata bs_metadata; + bs_metadata.platform = store_metadata.platform; + bs_metadata.sensor = store_metadata.sensor; + bs_metadata.survey = store_metadata.survey; + bs_metadata.date = store_metadata.date; + + cube::ImportAccumulatorConfig regen_config; + regen_config.store_dir = store_dir; + regen_config.reference_store_dir = reference_store_dir; + // Match the store level to the sheet's actual (GGGS-snapped) cell size so the + // scatter routing + gather scratch stores tile identically. + regen_config.cell_size_m = + static_cast(geo_map_sheet.nominalCellSizeMeters()); + regen_config.bs_store_dir = bs_store_dir; + cube::BatchRegen regen(make_sheet, regen_config); + + if (!reference_store_dir.empty()) { + std::cout << "Reference-prior seeding from " << reference_store_dir + << " (lazy per-tile; predicted-only blunder gate, #96)." << std::endl; + } + std::cout << "Exact rebuild: scatter to per-tile buckets, then gather each tile " + "in one unbounded pass (cube#96)." << std::endl; + + std::cout << "reading messages..." << std::endl; + + int ping_count = 0; + auto last_report_time = std::chrono::system_clock::now(); + const std::chrono::seconds report_interval(1); + + // SINGLE INTERLEAVED PASS over the chronological message stream. TF is fed + // into the bounded-window buffer as it arrives; each detection waits in a + // short FIFO until the TF frontier has advanced a guard interval past its + // stamp, so its bracketing transforms (both sides) are present before we + // project -- the correctness the old whole-bag 2-pass bought, but without the + // O(n) lookup cost of a 100k-entry-per-frame TimeCache (cube_bathymetry#63). + auto phase_tp = std::chrono::steady_clock::now(); + auto phase_secs = [&phase_tp]() { + auto now = std::chrono::steady_clock::now(); + double s = std::chrono::duration(now - phase_tp).count(); + phase_tp = now; + return s; + }; + + const int64_t begin_ns = cell_timestamp_ns; + const int64_t total_ns = + std::chrono::duration_cast(total_duration).count(); + + std::map speed_by_ns; + size_t tf_count = 0; + int64_t tf_frontier_ns = std::numeric_limits::min(); + + // Detections awaiting their bracketing TF (ping stamp ns -> message). + std::deque> pending; + + bool limit_reached = false; + + // Project + georeference one ping into the GeoMapSheet. + auto process_detection = + [&](const marine_acoustic_msgs::msg::SonarDetections & detections, int64_t ping_ns) { + // Real speed-over-ground from odom (Fix B) if available, else NaN -- the + // error model floors the speed-dependent horizontal-TPU terms to 0 (Fix A). + const float vessel_speed = speedAt(speed_by_ns, ping_ns); + + auto projection = projector.project(detections, tfBuffer, vessel_speed); + ++proj_pings; + proj_soundings += projection.soundings.size(); + proj_filtered_range += projection.diagnostics.filtered_range; + proj_missing_attitude += projection.diagnostics.missing_attitude; + proj_missing_heave += projection.diagnostics.missing_heave; + + try { + auto transform = tfBuffer.lookupTransform( + "earth", detections.header.frame_id, detections.header.stamp); + + // Per-ping geodetic reference: convert the sensor origin ECEF -> + // lat/long ONCE (the single iterative ECEF->geodetic solve per ping), + // then map each sounding's ECEF position through a local-ENU tangent + // plane (a cheap 4x4 matvec) + a first-order radii-of-curvature + // linearization. This collapses the per-sounding iterative solves into + // one per ping. At swath scale (tens of metres about the reference) the + // linearization error is well under a millimetre -- far below the GGGS + // cell size and the soundings' own TPU. + const auto & origin = transform.transform.translation; + gz4d::GeoPointLatLongDegrees ref_ll( + gz4d::GeoPointECEF(origin.x, origin.y, origin.z)); + gz4d::LocalENU enu(ref_ll); + + const double ref_lat_deg = ref_ll.latitude(); + const double ref_lon_deg = ref_ll.longitude(); + const double ref_height = ref_ll.altitude(); + constexpr double kDeg2Rad = M_PI / 180.0; + constexpr double kRad2Deg = 180.0 / M_PI; + constexpr double kA = 6378137.0; // WGS84 semi-major axis + constexpr double kF = 1.0 / 298.257223563; // WGS84 flattening + constexpr double kE2 = kF * (2.0 - kF); // first eccentricity^2 + const double lat0_rad = ref_lat_deg * kDeg2Rad; + const double sin_lat0 = std::sin(lat0_rad); + const double cos_lat0 = std::cos(lat0_rad); + const double w = std::sqrt(1.0 - kE2 * sin_lat0 * sin_lat0); + const double prime_vertical = kA / w; // N(lat0) + const double meridional = kA * (1.0 - kE2) / (w * w * w); // M(lat0) + + std::vector soundings; + soundings.reserve(projection.soundings.size()); + for (const auto & s : projection.soundings) { + geometry_msgs::msg::PointStamped sounding_re_sensor; + sounding_re_sensor.point.x = s.sonar_relative_position.x; + sounding_re_sensor.point.y = s.sonar_relative_position.y; + sounding_re_sensor.point.z = s.sonar_relative_position.z; + sounding_re_sensor.header = detections.header; + + geometry_msgs::msg::PointStamped sounding_ecef; + tf2::doTransform(sounding_re_sensor, sounding_ecef, transform); + + // ECEF -> local ENU (East, North, Up), then linearize ENU -> geodetic + // delta about the per-ping reference latitude. + const gz4d::Point local = enu.toLocal(gz4d::GeoPointECEF( + sounding_ecef.point.x, sounding_ecef.point.y, sounding_ecef.point.z)); + const double lat_deg = ref_lat_deg + (local[1] / meridional) * kRad2Deg; + const double lon_deg = + ref_lon_deg + (local[0] / (prime_vertical * cos_lat0)) * kRad2Deg; + const double height = ref_height + local[2]; + + cube::GeoSounding gs(gz4d::GeoPointLatLongDegrees(lat_deg, lon_deg, height)); + gs.sounding.vertical_error = s.vertical_error; + gs.sounding.horizontal_error = s.horizontal_error; + // Carry the {raw intensity, beam angle} sufficient-stats pair so the + // CUBE node co-estimates backscatter (#54) on the winning depth + // hypothesis -- without this every beam has NaN intensity and the + // backscatter store (--bs-store, #80) accumulates nothing. node.cpp + // emits the value UNCORRECTED; the angle correction is cube#81. + gs.sounding.intensity = s.intensity; + gs.sounding.beam_angle = s.beam_angle; + // Per-beam slant range R = twtt*c/2 (set in the Sounding detections + // ctor) for the tier-2 TL correction (cube#87). + gs.sounding.slant_range = s.slant_range; + soundings.push_back(gs); + } + // Scatter this ping's soundings to their per-tile buckets on disk (cube#96). + // Nothing accumulates in RAM here; the gather (finalize) builds each tile. + regen.addBatch(soundings); + ping_count++; + } catch (const tf2::TransformException & e) { + // A ping with no earth transform in the (bounded) buffer at its stamp -- + // e.g. before the first earth fix, or a TF gap wider than the cache + // window. Counted (not just logged) so an empty or sparse import is + // diagnosable rather than silently dropped. + ++proj_dropped_georef; + if (proj_dropped_georef <= 5) { + std::cerr << "Transform Exception: " << e.what() << std::endl; + } + } + }; + + // Drain detections whose bracketing TF is now present (frontier advanced a + // guard interval past the ping stamp). With flush=true, project whatever + // remains at end-of-stream against the available coverage. + auto drain_pending = [&](bool flush) { + while (!pending.empty() && !limit_reached) { + const int64_t ping_ns = pending.front().first; + if (!flush && (tf_frontier_ns == std::numeric_limits::min() || + ping_ns > tf_frontier_ns - kGuardNs)) + { + break; + } + process_detection(pending.front().second, ping_ns); + pending.pop_front(); + if (ping_count_limit > 0 && ping_count >= ping_count_limit) { + limit_reached = true; + } + } + }; + + std::cout << "projecting detections (single interleaved pass)..." << std::endl; + for (auto message = bag_readers.next(); message && !limit_reached; + message = bag_readers.next()) + { + const bool is_tf = message->data_type == "tf2_msgs/msg/TFMessage"; + const bool is_odom = !odom_topic.empty() && + message->data_type == "nav_msgs/msg/Odometry" && + message->message->topic_name == odom_topic; + const bool is_detection = message->data_type == "marine_acoustic_msgs/msg/SonarDetections" && + message->message->topic_name == detections_topic; + + if (is_tf) { + const bool is_static = ends_with(message->message->topic_name, "/tf_static"); + if (!is_static && !ends_with(message->message->topic_name, "/tf")) { + continue; + } + try { + rclcpp::SerializedMessage sm(*message->message->serialized_data); + tf2_msgs::msg::TFMessage tf_message; + rclcpp::Serialization().deserialize_message(&sm, &tf_message); + for (const auto & t : tf_message.transforms) { + tfBuffer.setTransform(t, "", is_static); + ++tf_count; + if (!is_static) { + const int64_t ns = static_cast(t.header.stamp.sec) * 1000000000LL + + t.header.stamp.nanosec; + tf_frontier_ns = std::max(tf_frontier_ns, ns); + } + } + } catch (const std::exception & e) { + std::cerr << e.what() << '\n'; + } + drain_pending(false); + } else if (is_odom) { + try { + rclcpp::SerializedMessage sm(*message->message->serialized_data); + nav_msgs::msg::Odometry odom; + rclcpp::Serialization().deserialize_message(&sm, &odom); + const int64_t ns = static_cast(odom.header.stamp.sec) * 1000000000LL + + odom.header.stamp.nanosec; + speed_by_ns[ns] = std::hypot(odom.twist.twist.linear.x, odom.twist.twist.linear.y); + } catch (const std::exception & e) { + std::cerr << e.what() << '\n'; + } + } else if (is_detection) { + try { + rclcpp::SerializedMessage sm(*message->message->serialized_data); + marine_acoustic_msgs::msg::SonarDetections detections; + rclcpp::Serialization().deserialize_message( + &sm, &detections); + const int64_t ping_ns = static_cast(detections.header.stamp.sec) * 1000000000LL + + detections.header.stamp.nanosec; + pending.emplace_back(ping_ns, std::move(detections)); + } catch (const std::exception & e) { + std::cerr << e.what() << '\n'; + } + } + + auto now = std::chrono::system_clock::now(); + if (now >= last_report_time + report_interval && total_ns > 0 && + tf_frontier_ns != std::numeric_limits::min()) + { + double progress = (tf_frontier_ns - begin_ns) / static_cast(total_ns); + std::cout << "\r " << static_cast(100 * progress) << "% " << ping_count + << " pings (" << proj_dropped_georef << " dropped, " << pending.size() + << " pending) " << std::flush; + last_report_time = now; + } + } + // Flush detections still pending at end-of-stream against available coverage. + drain_pending(true); + if (limit_reached) { + std::cout << "\nPing count limit of " << ping_count_limit << " reached" << std::endl; + } + std::cout << "\n buffered " << tf_count << " transforms"; + if (!odom_topic.empty()) {std::cout << ", " << speed_by_ns.size() << " odometry samples";} + std::cout << "; projected " << ping_count << " pings in " << phase_secs() << "s." << std::endl; + + std::cout << "\ndone." << std::endl; + std::cout << "Offline projection: " << proj_pings << " pings projected, " + << ping_count << " georeferenced into the grid, " << proj_dropped_georef + << " dropped (no earth TF); " << proj_soundings << " soundings (" + << proj_filtered_range << " range-filtered, " << proj_missing_attitude + << " missing attitude, " << proj_missing_heave << " missing heave)" << std::endl; + if (proj_pings > 0 && proj_soundings == 0) { + std::cerr << "WARNING: projected 0 soundings from " << proj_pings + << " pings -- check the --*-frame overrides match the bag's " + << "namespaced frames (see README 'Configuring frames per platform')." + << std::endl; + } + if (ping_count == 0 && proj_dropped_georef > 0) { + std::cerr << "WARNING: every ping was dropped for lack of an earth transform -- " + << "check that the bag has a localization chain to the 'earth' frame." + << std::endl; + } + + std::cout << "Gathering per-tile buckets (exact rebuild)..." << std::endl; + + // Gather every scattered tile bucket: each tile is rebuilt in a single unbounded + // pass over the complete set of soundings that touch it (no eviction), then + // written once to the `survey` layer (uma#248 collapsed draft/processed). Store- + // level metadata is written once here. The scatter scratch dir is cleaned up. + const std::size_t tile_buckets = regen.tileCount(); + regen.finalize( + store_metadata.empty() ? nullptr : &store_metadata, + (bs_store_dir.empty() || bs_metadata.empty()) ? nullptr : &bs_metadata); + std::cout << "Persisted " << regen.bathyTilesPersisted() + << " bathy tile(s) to " << store_dir << " (survey layer; " + << tile_buckets << " tile bucket(s) gathered; build: " + << phase_secs() << "s)." << std::endl; + + if (regen.bathyTilesPersisted() == 0) { + std::cerr << "WARNING: no tiles had finite data -- nothing imported. Check " + "the projector frame overrides and the detections topic." << std::endl; + } + + // The co-estimated backscatter was surfaced into the --bs-store `survey` layer + // (#80) from the SAME CUBE pass. By default UNCORRECTED; --backscatter-correction + // empirical applies the per-beam angular-response correction at node-output (#81). + if (!bs_store_dir.empty()) { + std::cout << "Persisted " << regen.backscatterTilesPersisted() + << " backscatter tile(s) to " << bs_store_dir << "." << std::endl; + if (regen.backscatterTilesPersisted() == 0) { + std::cerr << "WARNING: no cells had finite backscatter -- nothing written to " + "the backscatter store. Check that the detections carry intensities." + << std::endl; + } + } + + std::cout << "done!" << std::endl; + return 0; +} diff --git a/cube_bathymetry/src/cube_bathymetry_node.cpp b/cube_bathymetry/src/cube_bathymetry_node.cpp index fb0b431..3c99beb 100644 --- a/cube_bathymetry/src/cube_bathymetry_node.cpp +++ b/cube_bathymetry/src/cube_bathymetry_node.cpp @@ -229,13 +229,15 @@ class CubeBathymetry : public rclcpp_lifecycle::LifecycleNode marine_bathymetry_store::BathymetryStore::fromCellSize( static_cast(cell_size_)); marine_bathymetry_store::load(store, draft_dir_); + // The live node writes (and re-primes from) the `survey` layer since + // uma#248 collapsed the draft/processed split into one (ADR-0001 addendum). const auto & draft_tiles = - store.tiles(marine_bathymetry_store::SourceLayer::Draft); + store.tiles(marine_bathymetry_store::SourceLayer::Survey); if (!draft_tiles.empty()) { cube::loadIntoSheet( - store, marine_bathymetry_store::SourceLayer::Draft, *geo_map_sheet_); + store, marine_bathymetry_store::SourceLayer::Survey, *geo_map_sheet_); RCLCPP_INFO(get_logger(), - "Primed GeoMapSheet from %zu draft tiles under %s", + "Primed GeoMapSheet from %zu survey tiles under %s", draft_tiles.size(), draft_dir_.c_str()); // Bound the prime to the resident budget (must-fix): loadIntoSheet loads // the WHOLE store, so without this a restart mid-long-survey re-creates @@ -761,7 +763,7 @@ class CubeBathymetry : public rclcpp_lifecycle::LifecycleNode const auto ne = index.northEastPosition(); marine_bathymetry_store::loadWindow(scratch, draft_dir_, sw, ne, nullptr); const auto & tiles = - scratch.tiles(marine_bathymetry_store::SourceLayer::Draft); + scratch.tiles(marine_bathymetry_store::SourceLayer::Survey); auto it = tiles.find(index); if(it != tiles.end()) { cube::primeFromTile(it->second, *geo_map_sheet_); @@ -777,8 +779,9 @@ class CubeBathymetry : public rclcpp_lifecycle::LifecycleNode } // Persist every grid touched since the last save as a marine_bathymetry_store - // draft tile (atomic temp-then-rename via tile_io::saveTile), then clear the - // dirty set. A no-op when persistence is disabled or nothing changed. + // draft tile (via tile_io::saveTile -- a direct, flush/close-checked write, not + // crash-atomic), then clear the dirty set. A no-op when persistence is disabled or + // nothing changed. // // NOTE: geoGridToTile() calls GeoGrid::values(), which flushes the median // pre-filter -- the same flush the end-of-session export does, now happening @@ -795,14 +798,16 @@ class CubeBathymetry : public rclcpp_lifecycle::LifecycleNode return; } - const int64_t ts_ns = get_clock()->now().nanoseconds(); - // Single fused draft grid (unh_marine_autonomy#221): tiles go directly under - // /draft/ with no per-day epoch segment. Newest value wins per + // Single fused survey grid (unh_marine_autonomy#221, #248): tiles go directly + // under /survey/ with no per-day epoch segment. The live node now + // writes the `survey` layer directly (uma#248 collapsed the draft/processed + // split; ADR-0001 addendum): a bounded/approximate boat-side product that an + // off-boat batch-regen later overwrites as authoritative. Newest value wins per // cell, so successive saves and sessions accumulate into one grid. const std::string dir = draft_dir_ + "/" + marine_bathymetry_store::layerDirName( - marine_bathymetry_store::SourceLayer::Draft); + marine_bathymetry_store::SourceLayer::Survey); std::size_t written = 0; try { @@ -812,10 +817,11 @@ class CubeBathymetry : public rclcpp_lifecycle::LifecycleNode if (!grid_ptr) { continue; } - // source_index 0 = no registry wired yet (#21 scope); follow-on once - // marine_control device-control lands. + // BathyCell is 2-band since uma#248 (no per-cell timestamp/source_index); + // coarse provenance is store-level StoreMetadata (not wired in the live + // node yet, #21 scope). marine_bathymetry_store::BathymetryTile tile = - cube::geoGridToTile(*grid_ptr, ts_ns, /*source_index=*/0); + cube::geoGridToTile(*grid_ptr); if (!tile.dirty()) { // No finite cells (all queued/NaN) -- nothing to write. dirty() is the // value-raster flag geoGridToTile sets iff it wrote a finite cell diff --git a/cube_bathymetry/src/import_bag_main.cpp b/cube_bathymetry/src/import_bag_main.cpp index 74f082e..183a83c 100644 --- a/cube_bathymetry/src/import_bag_main.cpp +++ b/cube_bathymetry/src/import_bag_main.cpp @@ -31,6 +31,7 @@ #include #include #include +#include #include #include #include @@ -39,6 +40,7 @@ #include #include #include +#include #include #include #include @@ -68,18 +70,19 @@ { std::cout << "usage: import_bag [options] -o " "-d [ ...]\n"; - std::cout << " --bathy-layer draft|processed: bathymetry store layer for this " - "import (default processed: the off-boat CUBE re-run is authoritative; the live " - "node writes draft. #85)\n"; std::cout << " -o : Output bathymetry-store directory (created if " - "needed)\n"; - std::cout << " --prior : seed the CUBE predicted surface from this " - "store's Chart (contour) layer so blunder rejection drops false-deep " - "detections (#89). The chart must be at the survey GGGS level. NOTE: a " - "coarse/shallow-biased chart can also reject LEGITIMATE deeper-than-charted " - "returns; the rejection margin is tunable via the blunder_* params\n"; - std::cout << " --bs-store : Also write an MBES backscatter store layer " - "(Processed) from the same CUBE pass (optional; surfaces the co-estimated " + "needed). The off-boat CUBE re-run is authoritative, so it always writes the " + "`survey` layer (uma#248 collapsed the draft/processed split into one)\n"; + std::cout << " --reference-store : seed the CUBE predicted surface " + "lazily, per tile on first touch, from this store's `reference` (prior) layer " + "so blunder rejection drops false-deep detections (#89, #96). Only tiles at the " + "survey GGGS level gate. Predicted-only: the coarse prior is NEVER settled as " + "measured data and seeds no backscatter. NOTE: a coarse/shallow-biased prior " + "can also reject LEGITIMATE deeper-than-charted returns; the rejection margin " + "is tunable via the blunder_* params. (A `survey` tile already in -o takes " + "precedence and is warm-started as measured data instead.)\n"; + std::cout << " --bs-store : Also write an MBES backscatter store `survey` " + "layer from the same CUBE pass (optional; surfaces the co-estimated " "intensity -- uncorrected by default, angle-corrected when " "--backscatter-correction empirical is set, cube#81)\n"; std::cout << " -d : marine_acoustic_msgs/SonarDetections " @@ -109,10 +112,10 @@ "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"; - std::cout << " --platform / --sensor / --sensor-class / --campaign : " - "registry provenance fields\n"; + std::cout << " --platform / --sensor / --campaign : store-level provenance " + "written once to /registry.json (uma#248 StoreMetadata; --campaign maps " + "to the survey/campaign id). Per-cell source interning was retired for the " + "single-platform deployment\n"; std::cout << " Frame/range overrides for the offline projector (must match " "the bag's namespaced frames, or the grid comes out empty):\n"; std::cout << " --base-link-frame, --level-frame, --tide-frame\n"; @@ -293,9 +296,9 @@ int main(int argc, char * argv[]) std::vector bagfile_names; std::string store_dir; - std::string prior_dir; // optional: Chart-prior store to seed the predicted surface (#89) + // optional: reference-prior store to seed the predicted surface lazily (#89, #96) + std::string reference_store_dir; std::string bs_store_dir; // optional: MBES backscatter store output (#80) - std::string bathy_layer_str = "processed"; // bathy target layer (#85): draft|processed std::string detections_topic; // required std::string odom_topic; // optional: nav_msgs/Odometry for per-ping vessel speed double resolution = 1.0; @@ -310,9 +313,8 @@ int main(int argc, char * argv[]) std::string backscatter_correction_str = "none"; std::string backscatter_curve_file; - // Registry provenance fields for the imported cells. - marine_bathymetry_store::SourceRecord source_record; - source_record.source_id = "cube-replay"; + // Store-level provenance (uma#248 StoreMetadata, written once at finalize). + marine_bathymetry_store::StoreMetadata store_metadata; // DetectionsProjector configuration for the offline path. Defaults match // detections_to_pointcloud's node defaults so a detections-only bag projects @@ -334,15 +336,60 @@ int main(int argc, char * argv[]) return *arg; }; + // Guarded numeric parses: std::stod/std::stoi/std::stoll THROW on a non-numeric value + // and, uncaught in main, would std::terminate the process with an opaque message. + // Catch and route to usage() with a clear diagnostic instead (mirrors the sibling + // batch_regen_main). usage() is [[noreturn]], so these never fall through. + auto parse_double = [&](const char * flag, const std::string & value) -> double { + try { + std::size_t consumed = 0; + const double parsed = std::stod(value, &consumed); + if (consumed != value.size()) { + throw std::invalid_argument("trailing characters"); + } + return parsed; + } catch (const std::exception &) { + std::cerr << "error: option '" << flag << "' expects a number, got '" + << value << "'\n"; + usage(); + } + }; + auto parse_int = [&](const char * flag, const std::string & value) -> int { + try { + std::size_t consumed = 0; + const int parsed = std::stoi(value, &consumed); + if (consumed != value.size()) { + throw std::invalid_argument("trailing characters"); + } + return parsed; + } catch (const std::exception &) { + std::cerr << "error: option '" << flag << "' expects an integer, got '" + << value << "'\n"; + usage(); + } + }; + auto parse_long = [&](const char * flag, const std::string & value) -> int64_t { + try { + std::size_t consumed = 0; + const int64_t parsed = std::stoll(value, &consumed); + if (consumed != value.size()) { + throw std::invalid_argument("trailing characters"); + } + return parsed; + } catch (const std::exception &) { + std::cerr << "error: option '" << flag << "' expects an integer, got '" + << value << "'\n"; + usage(); + } + }; + for (; arg != arguments.end(); arg++) { if (*arg == "-h") { usage(); } else if (*arg == "-o") { store_dir = next_value("-o"); - } else if (*arg == "--prior") { - prior_dir = next_value("--prior"); - } else if (*arg == "--bathy-layer") { - bathy_layer_str = next_value("--bathy-layer"); + } else if (*arg == "--reference-store") { + reference_store_dir = next_value("--reference-store"); } else if (*arg == "--bs-store") { bs_store_dir = next_value("--bs-store"); } else if (*arg == "-d") { @@ -350,7 +397,7 @@ int main(int argc, char * argv[]) } else if (*arg == "--odom-topic") { odom_topic = next_value("--odom-topic"); } else if (*arg == "-r") { - resolution = std::stod(next_value("-r")); + resolution = parse_double("-r", next_value("-r")); } else if (*arg == "--iho-order") { iho_order = next_value("--iho-order"); } else if (*arg == "--backscatter-correction") { @@ -358,24 +405,20 @@ int main(int argc, char * argv[]) } 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")); + const int64_t v = parse_long("--max-resident-tiles", 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") { - source_record.source_id = next_value("--source-id"); + ping_count_limit = parse_int("-l", next_value("-l")); } else if (*arg == "--platform") { - source_record.platform = next_value("--platform"); + store_metadata.platform = next_value("--platform"); } else if (*arg == "--sensor") { - source_record.sensor = next_value("--sensor"); - } else if (*arg == "--sensor-class") { - source_record.sensor_class = next_value("--sensor-class"); + store_metadata.sensor = next_value("--sensor"); } else if (*arg == "--campaign") { - source_record.campaign = next_value("--campaign"); + store_metadata.survey = next_value("--campaign"); } else if (*arg == "--base-link-frame") { projector_params.base_link_frame = next_value("--base-link-frame"); } else if (*arg == "--level-frame") { @@ -383,9 +426,11 @@ int main(int argc, char * argv[]) } else if (*arg == "--tide-frame") { projector_params.tide_frame = next_value("--tide-frame"); } else if (*arg == "--minimum-range") { - projector_params.minimum_range = std::stod(next_value("--minimum-range")); + projector_params.minimum_range = parse_double("--minimum-range", + next_value("--minimum-range")); } else if (*arg == "--maximum-range") { - projector_params.maximum_range = std::stod(next_value("--maximum-range")); + projector_params.maximum_range = parse_double("--maximum-range", + next_value("--maximum-range")); } else if (!arg->empty() && (*arg)[0] == '-' && *arg != "-") { // An unrecognized flag would otherwise be silently treated as a bag path // and fail later with a confusing "cannot open bag". Reject it up front. @@ -408,19 +453,6 @@ int main(int argc, char * argv[]) usage(); } - // Bathy target layer (#85). The off-boat CUBE re-run is the authoritative product - // (the live node writes Draft), so this defaults to Processed; --bathy-layer draft - // overrides (e.g. to seed a draft from a bag). - marine_bathymetry_store::SourceLayer bathy_layer = - marine_bathymetry_store::SourceLayer::Processed; - if (bathy_layer_str == "draft") { - bathy_layer = marine_bathymetry_store::SourceLayer::Draft; - } else if (bathy_layer_str != "processed") { - std::cerr << "error: --bathy-layer must be 'draft' or 'processed' (got '" - << bathy_layer_str << "')\n"; - usage(); - } - std::cout << "Detections topic: " << detections_topic << " (offline projection, vessel_speed = NaN)" << std::endl; std::cout << "Store dir: " << store_dir << std::endl; @@ -519,117 +551,44 @@ int main(int argc, char * argv[]) backscatter_mode, std::move(backscatter_curve.points), backscatter_curve.tl_removed, backscatter_curve.absorption_db_per_m); - // Chart-prior prime (#89): seed the CUBE predicted surface from the prior store's - // Chart (contour) layer BEFORE any soundings are added. primeFromTile lazy-creates - // a node per chart cell carrying a predicted depth, which turns ON CUBE's - // predicted-surface blunder gate (Node::insert is a pass-through when - // predicted_depth_ is NaN); the survey soundings then hit nodes that already carry - // a predicted depth, so a false-deep detection below - // `target - blunder_scalar*sqrt(predicted_var)` is rejected. seed_settled=false: - // seed ONLY the predicted surface, never settle the coarse contour into the survey - // layer (that would contaminate both the bathy and co-estimated backscatter). - // - // ALIGNMENT: load() restores each Chart tile at the GGGS level encoded in its - // filename (the store is multi-level; the fromCellSize() arg below only sets the - // store's default cellIndex() level, it does NOT pin the level tiles load at). A - // chart cell only coincides with -- and thus gates -- a survey node when its level - // matches the survey level, i.e. the Chart layer was imported at the survey's - // resolution. A chart tile at a different level primes a node the survey soundings - // never land on, so it cannot gate. Rather than silently build an un-gated store - // (the operator passed --prior to reject outliers), we validate below and refuse - // when nothing would gate. A future refinement could resample the chart. - if (!prior_dir.empty()) { - marine_bathymetry_store::BathymetryStore prior = - marine_bathymetry_store::BathymetryStore::fromCellSize( - static_cast(geo_map_sheet.nominalCellSizeMeters())); - marine_bathymetry_store::SourceRegistry prior_reg; - try { - marine_bathymetry_store::load(prior, prior_dir, &prior_reg); - } catch (const std::exception & e) { - std::cerr << "ERROR: --prior store '" << prior_dir << "' failed to load: " - << e.what() << std::endl; - return 1; - } - const auto & chart_tiles = - prior.tiles(marine_bathymetry_store::SourceLayer::Chart); - const std::size_t prior_tiles = chart_tiles.size(); - if (prior_tiles == 0) { - std::cerr << "ERROR: --prior store '" << prior_dir << "' has no Chart-layer " - << "tiles; nothing would seed the predicted surface and blunder " - << "rejection would be INACTIVE. Check the path and that the store " - << "has a chart/ layer (refusing to build an un-gated store)." - << std::endl; - return 1; - } - // Only tiles at the survey level can coincide with a survey node and gate - // (see ALIGNMENT above). The prior store's default level is fromCellSize() of - // the survey cell size, so it equals the survey sheet's level. - const uint8_t survey_level = prior.level().level(); - std::size_t mismatched = 0; - for (const auto & grid_tile : chart_tiles) { - if (grid_tile.second.index().level() != survey_level) { - ++mismatched; - } - } - if (mismatched == prior_tiles) { - std::cerr << "ERROR: all " << prior_tiles << " Chart tile(s) in '" << prior_dir - << "' are at a GGGS level other than the survey level " - << static_cast(survey_level) << "; none would coincide with a " - << "survey node, so blunder rejection would be INACTIVE. Re-import " - << "the chart at the survey resolution (refusing to build an " - << "un-gated store)." << std::endl; - return 1; - } - if (mismatched > 0) { - std::cerr << "WARNING: " << mismatched << " of " << prior_tiles << " Chart " - << "tile(s) are not at the survey level " - << static_cast(survey_level) << " and will not gate; only the " - << "level-matched tiles seed the predicted surface." << std::endl; - } - cube::loadIntoSheet( - prior, marine_bathymetry_store::SourceLayer::Chart, geo_map_sheet, - /*seed_settled=*/false); - std::cout << "Primed CUBE predicted surface from " << (prior_tiles - mismatched) - << " of " << prior_tiles << " Chart tile(s) in " << prior_dir - << " (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); - } + // Reference-prior seeding (#89, #96) is now LAZY, per tile on first touch, driven + // by the accumulator's seedNewTile (see ImportAccumulatorConfig::reference_store_dir + // below): a `reference/` tile primes the CUBE predicted surface only (seed_settled= + // false) so the blunder gate turns on WITHOUT settling coarse prior depths as + // measured data. Only tiles at the survey GGGS level coincide with a survey node + // and gate; a reference tile at another level primes a node the soundings never + // land on (harmless no-op). The pre-#96 upfront whole-store loadIntoSheet was + // removed: it defeated the bounded-RAM eviction by loading the entire prior into + // the sheet at once. + + // Store-level provenance (uma#248 StoreMetadata) + bounded-RAM accumulator + // (cube#92). The accumulator owns the persist-then-drop eviction + lossless + // reload-on-revisit that bounds resident RAM by tile COUNT (not surveyed AREA) + // and the two-rung seed precedence (#96). registry.json is written once at + // finalize(). Backscatter provenance mirrors the bathy platform/sensor with an + // MBES-specific calibration ref (empty until a beam-pattern calibration exists). + marine_mbes_backscatter_store::StoreMetadata bs_metadata; + bs_metadata.platform = store_metadata.platform; + bs_metadata.sensor = store_metadata.sensor; + bs_metadata.survey = store_metadata.survey; + bs_metadata.date = store_metadata.date; 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; + accumulator_config.reference_store_dir = reference_store_dir; // Match the store level to the sheet's actual (GGGS-snapped) cell size so the - // reload/merge scratch stores tile identically. + // reload/seed/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 (!reference_store_dir.empty()) { + std::cout << "Reference-prior seeding from " << reference_store_dir + << " (lazy per-tile; predicted-only blunder gate, #96)." << std::endl; + } + if (max_resident_tiles > 0) { std::cout << "Bounded resident tiles: " << max_resident_tiles << " (cold tiles persist to the -o store and drop from RAM; " @@ -889,21 +848,21 @@ int main(int argc, char * argv[]) std::cout << "Building store tiles..." << std::endl; - // 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). + // Persist the still-resident tiles and write the store-level metadata. Tiles + // evicted during the pass were already written to disk (bathy in the -o store, + // their backscatter to --bs-store); 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 it always writes the `survey` layer (uma#248 + // collapsed draft/processed). Single fused grid per layer (uma#221). 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); + store_metadata.empty() ? nullptr : &store_metadata, + (bs_store_dir.empty() || bs_metadata.empty()) ? nullptr : &bs_metadata); std::cout << "Persisted " << accumulator.bathyTilesPersisted() - << " bathy tile(s) to " << store_dir << " (" << bathy_layer_str - << " layer; " << evicted_count << " evicted mid-pass, " + << " bathy tile(s) to " << store_dir << " (survey layer; " + << evicted_count << " evicted mid-pass, " << resident_before_final << " resident at end; build: " << phase_secs() << "s)." << std::endl; diff --git a/cube_bathymetry/src/store_import.cpp b/cube_bathymetry/src/store_import.cpp index 7fab67a..44908c4 100644 --- a/cube_bathymetry/src/store_import.cpp +++ b/cube_bathymetry/src/store_import.cpp @@ -61,8 +61,7 @@ namespace constexpr double kPrimeVarianceFloor = 1e-4; } // namespace -marine_bathymetry_store::BathymetryTile geoGridToTile( - const GeoGrid & grid, int64_t timestamp_ns, uint16_t source_index) +marine_bathymetry_store::BathymetryTile geoGridToTile(const GeoGrid & grid) { marine_bathymetry_store::BathymetryTile tile(grid.index()); @@ -81,21 +80,21 @@ marine_bathymetry_store::BathymetryTile geoGridToTile( if (std::isnan(v.depth)) { continue; // no estimate here; leave the tile's NaN no-data sentinel } + // BathyCell is 2-band {depth, uncertainty} since uma#248 (the per-cell + // timestamp/source bands were dropped; coarse provenance is store-level + // StoreMetadata). values() returns float; widen to the store's double. tile.set( (*it).row(), (*it).column(), marine_bathymetry_store::BathyCell{ static_cast(v.depth), - static_cast(v.uncertainty), - timestamp_ns, - source_index}); + static_cast(v.uncertainty)}); } return tile; } std::map -mapSheetToTiles( - const GeoMapSheet & map_sheet, int64_t timestamp_ns, uint16_t source_index) +mapSheetToTiles(const GeoMapSheet & map_sheet) { std::map tiles; @@ -103,8 +102,7 @@ mapSheetToTiles( if (!grid) { continue; } - marine_bathymetry_store::BathymetryTile tile = - geoGridToTile(*grid, timestamp_ns, source_index); + marine_bathymetry_store::BathymetryTile tile = geoGridToTile(*grid); // Drop a grid that yielded no finite cells -- an all-no-data tile would just // persist as an empty file. if (!tile.dirty()) { @@ -117,8 +115,7 @@ mapSheetToTiles( } std::map -geoGridToBackscatterCells( - const GeoGrid & grid, int64_t timestamp_ns, uint16_t source_index) +geoGridToBackscatterCells(const GeoGrid & grid) { std::map cells; @@ -136,20 +133,34 @@ geoGridToBackscatterCells( if (std::isnan(r.intensity)) { continue; // no co-estimated backscatter here -- skip (mirrors NaN-depth skip) } - // intensity_var is the estimate variance (NaN with < 2 samples); it rides - // into the quality band (ADR-0007 D6). timestamp/source stamp provenance. - cells.emplace( - *it, - marine_mbes_backscatter_store::MbesCell{ - r.intensity, r.intensity_var, timestamp_ns, source_index}); + // Encode the Welford sufficient statistic into the 3-band MbesCell (uma#248 + // A.1); welfordFromCell is the exact inverse. r.intensity is the running mean; + // r.intensity_var is the estimate variance (variance of the mean = M2/(n-1)/n), + // NaN with < 2 samples. + marine_mbes_backscatter_store::MbesCell cell; + cell.mean = r.intensity; + if (std::isnan(r.intensity_var)) { + // n = 1 sentinel: a single beam has no dispersion. sample_sd == 0 with a + // finite mean reconstructs to n = 1 (distinct from mean = NaN no-data). + cell.standard_error = 0.0f; + cell.sample_sd = 0.0f; + } else { + // n >= 2. sample_sd = sqrt(intensity_var * n) is the beams' sample stddev; + // standard_error = scale * sqrt(intensity_var) is the confidence-scaled SE + // of the mean (true SE = sqrt(intensity_var)). + const float n = static_cast(r.n_samples); + cell.sample_sd = std::sqrt(r.intensity_var * n); + cell.standard_error = + kBackscatterConfidenceScale * std::sqrt(r.intensity_var); + } + cells.emplace(*it, cell); } return cells; } std::map -mapSheetToBackscatterCells( - const GeoMapSheet & map_sheet, int64_t timestamp_ns, uint16_t source_index) +mapSheetToBackscatterCells(const GeoMapSheet & map_sheet) { std::map cells; @@ -158,7 +169,7 @@ mapSheetToBackscatterCells( continue; } std::map grid_cells = - geoGridToBackscatterCells(*grid, timestamp_ns, source_index); + geoGridToBackscatterCells(*grid); // Grids cover disjoint GGGS cells, so merge never collides. cells.merge(grid_cells); } @@ -166,6 +177,48 @@ mapSheetToBackscatterCells( return cells; } +IntensityWelford welfordFromCell( + const marine_mbes_backscatter_store::MbesCell & cell) +{ + IntensityWelford w; + if (!cell.hasData()) { + return w; // no-data cell (mean NaN) -> empty Welford {n=0} + } + w.mean = static_cast(cell.mean); + if (cell.sample_sd == 0.0f) { + // n = 1 sentinel (a single-beam node, or -- the accepted limitation -- an n>=2 + // cell whose samples were all exactly identical, M2 == 0). Reconstruct n = 1. + w.n = 1; + w.m2 = 0.0; + return w; + } + // Divide the confidence scale back out of standard_error to recover the true + // standard error of the mean, then invert n = (sample_sd / SE)^2 and + // M2 = sample_sd^2 * (n - 1). round() absorbs float round-trip error in n. + const double sample_sd = static_cast(cell.sample_sd); + const double se = + static_cast(cell.standard_error) / kBackscatterConfidenceScale; + if (!(se > 0.0)) { + // Defensive: a well-formed encode pairs a non-zero sample_sd with a non-zero + // standard_error (both are zero iff the variance is zero, handled above). A + // corrupt or hand-built tile with sample_sd != 0 but a zero / negative / + // non-finite standard_error would drive ratio = sample_sd / se to +/-inf or + // NaN, and std::lround(inf/NaN) is undefined behaviour with an out-of-range + // uint32_t cast. A zero SE carries no dispersion information, so fall back to + // the n = 1 sentinel rather than invoke UB. + w.n = 1; + w.m2 = 0.0; + return w; + } + const double ratio = sample_sd / se; + w.n = static_cast(std::lround(ratio * ratio)); + if (w.n < 2) { + w.n = 2; // sample_sd != 0 implies >= 2 samples; guard against round-to-1 + } + w.m2 = sample_sd * sample_sd * static_cast(w.n - 1); + return w; +} + void primeFromTile( const marine_bathymetry_store::BathymetryTile & tile, GeoMapSheet & map_sheet, bool seed_settled) @@ -277,6 +330,15 @@ std::size_t ImportAccumulator::residentTileCount() const return sheet_.residentTileCount(); } +void ImportAccumulator::persistResidentTile(const gggs::GridIndex & index) +{ + // Batch-regen gather (#96): write only this tile (bathy + backscatter). The + // gather sheet may also hold neighbour grids a near-seam sounding spilled into; + // those are each written by their OWN tile's gather, so we never persist them here. + persistBathyTile(index); + persistBackscatterTile(index); +} + void ImportAccumulator::spillIntensitySamples(const gggs::GridIndex & index) { if (cfg_.bs_store_dir.empty()) { @@ -354,16 +416,21 @@ void ImportAccumulator::persistBathyTile(const gggs::GridIndex & index) } // 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); + marine_bathymetry_store::BathymetryTile tile = geoGridToTile(*grid); if (!tile.dirty()) { return; } + // The off-boat CUBE re-run is the authoritative product: always the `survey` + // layer (uma#248 collapsed the draft/processed split into one). const std::string layer_dir = cfg_.store_dir + "/" + - marine_bathymetry_store::layerDirName(cfg_.bathy_layer); + marine_bathymetry_store::layerDirName( + marine_bathymetry_store::SourceLayer::Survey); 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). + // tile_io::saveTile writes the GTiff directly to the final path and checks the + // flush/close result (an I/O error or full disk throws), but it is NOT crash-atomic: + // it does not temp-then-rename, so a crash/kill mid-write can leave a partial tile + // at the final path. Making that write atomic is a tracked marine_tiled_raster_store + // follow-up (#96 review); a batch-regen re-run reproduces the tile exactly. marine_bathymetry_store::saveTile( tile, layer_dir + "/" + marine_bathymetry_store::tileFilename(index)); ++bathy_persisted_; @@ -381,7 +448,7 @@ void ImportAccumulator::persistBackscatterTile(const gggs::GridIndex & index) 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); + geoGridToBackscatterCells(*grid); 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 @@ -389,7 +456,7 @@ void ImportAccumulator::persistBackscatterTile(const gggs::GridIndex & index) return; } const std::string layer_dir = cfg_.bs_store_dir + "/" + - mbs::layerDirName(mbs::SourceLayer::Processed); + mbs::layerDirName(mbs::SourceLayer::Survey); 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 @@ -455,7 +522,8 @@ bool ImportAccumulator::reloadEvictedTile(const gggs::GridIndex & index) 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); + const auto & tiles = + scratch.tiles(marine_bathymetry_store::SourceLayer::Survey); auto it = tiles.find(index); if (it != tiles.end()) { // seed_settled=true: restore each cell's settled depth/uncertainty as a CUBE @@ -481,6 +549,100 @@ bool ImportAccumulator::reloadEvictedTile(const gggs::GridIndex & index) return true; } +bool ImportAccumulator::seedNewTile(const gggs::GridIndex & index) +{ + // Two-rung seed precedence (#96), run once per tile on first touch. survey wins + // over reference: a survey tile is measured CUBE data (settle it, seed its + // backscatter); a reference tile is a coarse read-only prior (gate only). + + // Rung 1 -- survey: a pre-existing survey bathy tile (an incremental import into + // an existing store, or -- via reloadEvictedTile -- an already-written tile). + // loadWindow silently returns 0 when store_dir has no survey/ layer yet (a fresh + // import), so this falls through to the reference rung with no error/warning. + // skip_survey_seed disables this rung for the batch-regen gather: replaying a + // tile's complete sounding population onto a warm-start from that same tile in the + // OUTPUT store would double-count it (a silent blend, not an exact rebuild). + if (!cfg_.store_dir.empty() && !cfg_.skip_survey_seed) { + 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(marine_bathymetry_store::SourceLayer::Survey); + auto it = tiles.find(index); + if (it != tiles.end()) { + // Settled warm-start: the survey layer round-trips as a CUBE hypothesis and + // refines under new soundings (ADR-0001). + primeFromTile(it->second, sheet_, /*seed_settled=*/true); + // Reconstruct each cell's corrected-intensity Welford from the survey + // backscatter tile so the re-run's beams blend with the stored population + // (lossless backscatter seed). welfordFromCell inverts the 3-band write. + if (!cfg_.bs_store_dir.empty()) { + namespace mbs = marine_mbes_backscatter_store; + const std::string bs_path = cfg_.bs_store_dir + "/" + + mbs::layerDirName(mbs::SourceLayer::Survey) + "/" + + mbs::tileFilename(index); + if (std::filesystem::is_regular_file(bs_path)) { + const gggs::Level level = gggs::Level::fromCellSize(cfg_.cell_size_m); + const mbs::MbesTile bs_tile = mbs::loadTile(bs_path, level); + gggs::CellAreaIterator cit(index); + for (; cit.valid(); cit.next()) { + const mbs::MbesCell c = bs_tile.get((*cit).row(), (*cit).column()); + if (c.hasData()) { + sheet_.setSettledIntensityWelfordAt(*cit, welfordFromCell(c)); + } + } + } + } + seeded_.insert(index); + return true; + } + } catch (const std::exception & e) { + // A survey-seed read error means the on-disk survey tile EXISTS but could not + // be loaded (a fresh import returns 0 tiles WITHOUT throwing, so it never + // reaches here). Accumulating from scratch and then persisting would overwrite + // that intact-but-unreadable tile with partial data -- the same data-loss the + // reload path guards against. Mirror reloadEvictedTile: signal failure so the + // caller drops this tile (and its soundings) to protect the on-disk surface, + // and leave it UNseeded so a later batch retries the seed. + std::cerr << "import_bag: could not survey-seed tile on first touch: " + << e.what() << " (dropping this batch's soundings on the tile to " + "protect the on-disk surface; a later batch retries the seed, but these " + "dropped soundings are not re-processed)" << std::endl; + return false; + } + } + + // Rung 2 -- reference: a coarse read-only prior. Predicted-only prime + // (seed_settled=false) turns the blunder gate on but does NOT settle the cell as + // measured data (no values() output, sheet stays clean) and seeds NO backscatter. + if (!cfg_.reference_store_dir.empty()) { + try { + marine_bathymetry_store::BathymetryStore ref = + marine_bathymetry_store::BathymetryStore::fromCellSize(cfg_.cell_size_m); + const auto sw = index.southWestPosition(); + const auto ne = index.northEastPosition(); + marine_bathymetry_store::loadWindow( + ref, cfg_.reference_store_dir, sw, ne, nullptr); + const auto & tiles = + ref.tiles(marine_bathymetry_store::SourceLayer::Reference); + auto it = tiles.find(index); + if (it != tiles.end()) { + primeFromTile(it->second, sheet_, /*seed_settled=*/false); + } + } catch (const std::exception & e) { + std::cerr << "import_bag: could not reference-seed tile on first touch: " + << e.what() << " (no prior gate for this tile)" << std::endl; + } + } + + // else blank -- nothing to seed; still mark it seeded so we do not retry. + seeded_.insert(index); + return true; +} + void ImportAccumulator::evictColdTiles() { if (cfg_.max_resident_tiles == 0) { @@ -530,15 +692,25 @@ void ImportAccumulator::addBatch( // 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. + // For each tile this batch is about to touch: reload it if it was evicted, else + // seed it if this is its first touch (seed precedence #96 -- survey warm-start, + // reference gate, or blank). Both run BEFORE the soundings are added so the new + // beams accrete onto the seeded/reloaded hypotheses (lossless blend). Iterating + // the window every batch is cheap; seedNewTile is a marked-once no-op thereafter. 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); - } + for (const auto & idx : sheet_.gridIndicesForSoundings(soundings)) { + if (evicted_.count(idx)) { + if (reloadEvictedTile(idx)) { + evicted_.erase(idx); + } else { + reload_failed.push_back(idx); + } + } else if (!seeded_.count(idx)) { + if (!seedNewTile(idx)) { + // Survey-seed read error on an existing tile: protect it exactly like a + // failed reload -- drop this batch's soundings on it after the add so the + // intact-but-unreadable on-disk surface is never overwritten from scratch. + reload_failed.push_back(idx); } } } @@ -562,9 +734,9 @@ void ImportAccumulator::addBatch( } 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; + << " un-reloadable/un-seedable tile(s) this batch (reload or " + "survey-seed 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 } @@ -574,8 +746,8 @@ void ImportAccumulator::addBatch( } void ImportAccumulator::finalize( - const marine_bathymetry_store::SourceRegistry & bathy_registry, - const marine_mbes_backscatter_store::SourceRegistry * bs_registry) + const marine_bathymetry_store::StoreMetadata * bathy_metadata, + const marine_mbes_backscatter_store::StoreMetadata * bs_metadata) { // Persist every still-resident tile (evicted tiles are already durable). Snapshot // the indices first so the persist loop iterates a stable, deterministic order. @@ -589,15 +761,20 @@ void ImportAccumulator::finalize( 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()) { + // Store-level provenance sidecars (uma#248 replaced the per-cell SourceRegistry + // interning table with one coarse StoreMetadata `registry.json` at the store + // root). Written once at the end; skipped when absent or empty. + if (!cfg_.store_dir.empty() && bathy_metadata != nullptr && + !bathy_metadata->empty()) + { std::filesystem::create_directories(cfg_.store_dir); - bathy_registry.saveRegistry(cfg_.store_dir); + bathy_metadata->save(cfg_.store_dir); } - if (!cfg_.bs_store_dir.empty() && bs_registry != nullptr) { + if (!cfg_.bs_store_dir.empty() && bs_metadata != nullptr && + !bs_metadata->empty()) + { std::filesystem::create_directories(cfg_.bs_store_dir); - bs_registry->saveRegistry(cfg_.bs_store_dir); + bs_metadata->save(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. diff --git a/cube_bathymetry/test/test_batch_regen.cpp b/cube_bathymetry/test/test_batch_regen.cpp new file mode 100644 index 0000000..0e07380 --- /dev/null +++ b/cube_bathymetry/test/test_batch_regen.cpp @@ -0,0 +1,360 @@ +// 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. + +// Batch-regen exactness (cube_bathymetry#96). The scatter/gather rebuild builds +// each tile in a single unbounded pass over the COMPLETE set of soundings that +// touch it, so its output is BIT-EXACT against a single-pass unbounded +// ImportAccumulator over the same soundings -- depth, uncertainty, and the 3-band +// backscatter all reproduced to the last bit (no eviction, no reload, no +// approximation). Also asserts the scatter scratch dir is cleaned up. + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "cube_bathymetry/batch_regen.h" +#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/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 +{ +namespace +{ +constexpr float kCellSize = 1.0f; + +std::string makeTempDir(const std::string & tag) +{ + const auto base = std::filesystem::temp_directory_path() / + ("cube_batch_regen_" + 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. Varying intensity gives real +// backscatter dispersion (so the 3-band statistic is non-trivial). +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 + rep * 0.7f; // real dispersion (n>=2) + s.sounding.beam_angle = 0.0f; + s.sounding.slant_range = depth; + soundings.push_back(s); + } + return soundings; +} + +BatchRegen::SheetFactory sheetFactory() +{ + return []() {return std::make_unique(kCellSize);}; +} + +ImportAccumulatorConfig makeConfig( + const std::string & store_dir, const std::string & bs_dir) +{ + ImportAccumulatorConfig cfg; + cfg.store_dir = store_dir; + cfg.cell_size_m = kCellSize; + cfg.bs_store_dir = bs_dir; + cfg.max_resident_tiles = 0; + return cfg; +} + +// Single-pass unbounded reference build. +void runSinglePass( + const std::vector> & batches, + const std::string & store_dir, const std::string & bs_dir) +{ + GeoMapSheet sheet(kCellSize); + ImportAccumulator acc(sheet, makeConfig(store_dir, bs_dir)); + for (const auto & batch : batches) { + acc.addBatch(batch); + } + acc.finalize(); +} + +// Batch-regen scatter/gather build. +void runBatchRegen( + const std::vector> & batches, + const std::string & store_dir, const std::string & bs_dir) +{ + BatchRegen regen(sheetFactory(), makeConfig(store_dir, bs_dir)); + for (const auto & batch : batches) { + regen.addBatch(batch); + } + regen.finalize(); +} + +std::map> loadBathyCells( + const std::string & store_dir) +{ + marine_bathymetry_store::BathymetryStore store = + marine_bathymetry_store::BathymetryStore::fromCellSize(kCellSize); + marine_bathymetry_store::load(store, store_dir); + std::map> out; + for (const auto & grid_tile : store.tiles( + marine_bathymetry_store::SourceLayer::Survey)) + { + 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; +} + +// The stored 3-band backscatter statistic per finite cell (mean, standard_error, +// sample_sd) -- compared byte-for-byte between the two builds. +struct BsBands +{ + float mean; + float standard_error; + float sample_sd; +}; + +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::load(store, bs_dir); + std::map out; + for (const auto & grid_tile : store.tiles( + marine_mbes_backscatter_store::SourceLayer::Survey)) + { + 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, BsBands{c.mean, c.standard_error, c.sample_sd}); + } + } + } + return out; +} + +// Assert batch-regen == single-pass, bit-exact, cell-for-cell. +void expectExactMatch( + const std::string & regen_dir, const std::string & single_dir, + const std::string & regen_bs, const std::string & single_bs) +{ + const auto bathy_r = loadBathyCells(regen_dir); + const auto bathy_s = loadBathyCells(single_dir); + ASSERT_FALSE(bathy_s.empty()) << "single-pass build produced no bathy cells"; + ASSERT_EQ(bathy_r.size(), bathy_s.size()) + << "batch-regen bathy cell count differs from single-pass"; + for (const auto & [cell, ds] : bathy_s) { + auto it = bathy_r.find(cell); + ASSERT_NE(it, bathy_r.end()) << "batch-regen is missing a bathy cell"; + EXPECT_DOUBLE_EQ(it->second.first, ds.first) << "depth differs (not bit-exact)"; + EXPECT_DOUBLE_EQ(it->second.second, ds.second) + << "uncertainty differs (not bit-exact)"; + } + + const auto bs_r = loadBackscatterCells(regen_bs); + const auto bs_s = loadBackscatterCells(single_bs); + ASSERT_FALSE(bs_s.empty()) << "single-pass build produced no backscatter cells"; + ASSERT_EQ(bs_r.size(), bs_s.size()) + << "batch-regen backscatter cell count differs from single-pass"; + for (const auto & [cell, bs] : bs_s) { + auto it = bs_r.find(cell); + ASSERT_NE(it, bs_r.end()) << "batch-regen is missing a backscatter cell"; + EXPECT_FLOAT_EQ(it->second.mean, bs.mean) << "backscatter mean differs"; + EXPECT_FLOAT_EQ(it->second.standard_error, bs.standard_error) + << "backscatter standard_error differs"; + EXPECT_FLOAT_EQ(it->second.sample_sd, bs.sample_sd) + << "backscatter sample_sd differs"; + } +} +} // namespace + +// A disjoint multi-tile survey: batch-regen must reproduce the single-pass store +// exactly. +TEST(BatchRegen, MultiTileExactMatch) +{ + std::vector> batches; + for (int i = 0; i < 12; ++i) { + // ~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("multitile"); + const std::string r_dir = root + "/regen"; + const std::string s_dir = root + "/single"; + const std::string r_bs = root + "/regen_bs"; + const std::string s_bs = root + "/single_bs"; + + runBatchRegen(batches, r_dir, r_bs); + runSinglePass(batches, s_dir, s_bs); + + expectExactMatch(r_dir, s_dir, r_bs, s_bs); + std::filesystem::remove_all(root); +} + +// A tile visited, then other tiles surveyed, then the first tile REVISITED at a +// different cell. In batch-regen both visits land in the same tile bucket and are +// gathered in one pass; the single-pass keeps the tile resident the whole time. +// Both must agree exactly. +TEST(BatchRegen, RevisitExactMatch) +{ + std::vector> batches; + batches.push_back(surveyCell(43.00000, -70.00000, 12.0f, 40.0f)); // visit 1, tile A + for (int i = 1; i <= 5; ++i) { // other tiles + batches.push_back(surveyCell(43.0 + 0.02 * i, -71.0, 15.0f + i, 25.0f + i)); + } + batches.push_back(surveyCell(43.00050, -70.00000, 18.0f, 55.0f)); // visit 2, tile A + + const std::string root = makeTempDir("revisit"); + const std::string r_dir = root + "/regen"; + const std::string s_dir = root + "/single"; + const std::string r_bs = root + "/regen_bs"; + const std::string s_bs = root + "/single_bs"; + + runBatchRegen(batches, r_dir, r_bs); + runSinglePass(batches, s_dir, s_bs); + + expectExactMatch(r_dir, s_dir, r_bs, s_bs); + std::filesystem::remove_all(root); +} + +// A survey straddling a GGGS tile seam with a MULTI-CELL deposit radius. This is +// the case the earlier per-sounding ±1-cell scatter got WRONG: GeoGrid::insert +// spreads a sounding out to radius ~= CONF_99PC*sqrt(horizontal_error) (several +// cells for realistic TPU), so a sounding two-plus cells inside one tile still +// deposits into the neighbour tile's cells across the seam. A single unbounded pass +// makes those cross-seam deposits (the neighbour grid receives the whole batch); +// per-sounding scatter dropped them because the sounding's own ±1-cell window never +// reached the neighbour tile. The other tests use a sub-cell horizontal_error (0.1), +// so they never cross a seam and cannot catch this -- this test pins the fix. +TEST(BatchRegen, SeamCrossingExactMatch) +{ + const gggs::Level level = gggs::Level::fromCellSize(kCellSize); + // East edge of the tile containing (43, -70) == west edge of its east neighbour. + const double seam_lon = level.gridIndex(43.0, -70.0).eastLongitude(); + const double lat = 43.0; + const double deg_per_m_lon = 1.0 / (111320.0 * std::cos(lat * M_PI / 180.0)); + + // horizontal_error = 4 m -> max_radius = CONF_99PC*sqrt(4) = 5.152 m. With a small + // vertical_error the effective radius (sqrt(ratio-1) - max_radius, clamped) settles + // near 3 m -- about three 1 m cells. vertical_error is chosen so ratio = + // maxVarianceAllowed(~12 m, order1a)/vertical_error ~= 67, i.e. sqrt(ratio-1) ~= + // 8.15 and radius ~= 3 m. Survey cells sit 1..4 m either side of the seam, so the + // ones >= 2 m from it deposit across the seam yet lie outside their own ±1-cell + // window -- exactly the deposits the old scatter dropped. + const float kHorizErr = 4.0f; + const float kVertErr = 0.00106f; + + std::vector batch; + for (int c = -4; c <= 4; ++c) { + if (c == 0) {continue;} + const double lon = seam_lon + c * 1.0 * deg_per_m_lon; + for (int rep = 0; rep < 8; ++rep) { + gz4d::GeoPointLatLongDegrees point(lat, lon, -12.0 - rep * 0.001); + GeoSounding s(point); + s.sounding.vertical_error = kVertErr; + s.sounding.horizontal_error = kHorizErr; + s.sounding.intensity = 30.0f + c + rep * 0.7f; // real dispersion (n>=2) + s.sounding.beam_angle = 0.0f; + s.sounding.slant_range = 12.0f; + batch.push_back(s); + } + } + const std::vector> batches{batch}; + + const std::string root = makeTempDir("seam"); + const std::string r_dir = root + "/regen"; + const std::string s_dir = root + "/single"; + const std::string r_bs = root + "/regen_bs"; + const std::string s_bs = root + "/single_bs"; + + runBatchRegen(batches, r_dir, r_bs); + runSinglePass(batches, s_dir, s_bs); + + // The survey must actually span both tiles, else the seam was never crossed and + // the test would be vacuous. + const auto cells = loadBathyCells(s_dir); + std::set grids; + for (const auto & kv : cells) { + grids.insert(kv.first.grid()); + } + ASSERT_GE(grids.size(), 2u) + << "seam survey must span >= 2 tiles to exercise cross-seam deposits"; + + expectExactMatch(r_dir, s_dir, r_bs, s_bs); + std::filesystem::remove_all(root); +} + +// The scatter scratch directory is created during scatter and DELETED by finalize. +TEST(BatchRegen, ScratchDirCleanedUpAfterFinalize) +{ + const std::string root = makeTempDir("scratch"); + const std::string store_dir = root + "/store"; + const std::string bs_dir = root + "/bs"; + + BatchRegen regen(sheetFactory(), makeConfig(store_dir, bs_dir)); + for (int i = 0; i < 6; ++i) { + regen.addBatch(surveyCell(43.0 + 0.02 * i, -70.0, 10.0f + i, 30.0f + i)); + } + const std::string scratch = regen.scratchDir(); + ASSERT_FALSE(scratch.empty()) << "scatter should have created a scratch dir"; + EXPECT_TRUE(std::filesystem::exists(scratch)) << "scratch dir should exist mid-run"; + + regen.finalize(); + + EXPECT_FALSE(std::filesystem::exists(scratch)) + << "finalize must delete the scatter scratch dir"; + EXPECT_TRUE(regen.scratchDir().empty()) << "finalize must clear the scratch path"; + EXPECT_GT(regen.bathyTilesPersisted(), 0u); + + std::filesystem::remove_all(root); +} + +} // namespace cube diff --git a/cube_bathymetry/test/test_import_eviction.cpp b/cube_bathymetry/test/test_import_eviction.cpp index 293794b..c0b2bef 100644 --- a/cube_bathymetry/test/test_import_eviction.cpp +++ b/cube_bathymetry/test/test_import_eviction.cpp @@ -58,7 +58,6 @@ namespace cube { namespace { -constexpr int64_t kStamp = 1234567890123456789LL; constexpr float kCellSize = 1.0f; std::string makeTempDir(const std::string & tag) @@ -95,12 +94,8 @@ ImportAccumulatorConfig makeConfig( { 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; } @@ -115,15 +110,9 @@ void runImport( 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); + // No StoreMetadata needed for the equivalence assertions (uma#248 moved provenance + // to an optional store-level sidecar; the tile bands carry the data under test). + accumulator.finalize(); } // Load every finite bathy cell of a store layer into a comparable map. @@ -132,11 +121,10 @@ std::map> loadBathyCells( { marine_bathymetry_store::BathymetryStore store = marine_bathymetry_store::BathymetryStore::fromCellSize(kCellSize); - marine_bathymetry_store::SourceRegistry reg; - marine_bathymetry_store::load(store, store_dir, ®); + marine_bathymetry_store::load(store, store_dir); std::map> out; for (const auto & grid_tile : store.tiles( - marine_bathymetry_store::SourceLayer::Processed)) + marine_bathymetry_store::SourceLayer::Survey)) { const auto & depth = grid_tile.second.depthBand(); const auto & unc = grid_tile.second.uncertaintyBand(); @@ -151,24 +139,43 @@ std::map> loadBathyCells( return out; } -// Load every finite backscatter cell of a store into a comparable map. +// Count finite cells in a store's `reference` layer (what a measured seed WOULD +// have contributed, for the reference-seed enforcement test). +std::size_t countReferenceFiniteCells(const std::string & store_dir) +{ + marine_bathymetry_store::BathymetryStore store = + marine_bathymetry_store::BathymetryStore::fromCellSize(kCellSize); + marine_bathymetry_store::load(store, store_dir); + std::size_t n = 0; + for (const auto & grid_tile : store.tiles( + marine_bathymetry_store::SourceLayer::Reference)) + { + for (double d : grid_tile.second.depthBand()) { + if (std::isfinite(d)) {++n;} + } + } + return n; +} + +// Load every finite backscatter cell of a store into a comparable map. Compares +// the RAW stored 3-band statistic (mean, sample_sd) so bounded vs unbounded is a +// byte-level losslessness check (no reconstruction, so no n=1-sentinel ambiguity). 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, ®); + marine_mbes_backscatter_store::load(store, bs_dir); std::map> out; for (const auto & grid_tile : store.tiles( - marine_mbes_backscatter_store::SourceLayer::Processed)) + marine_mbes_backscatter_store::SourceLayer::Survey)) { 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)); + out.emplace(*it, std::make_pair(c.mean, c.sample_sd)); } } } @@ -213,8 +220,8 @@ void expectStoresEqual( 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"; + EXPECT_NEAR(it->second.first, iv.first, 1e-3) << "backscatter mean differs"; + EXPECT_NEAR(it->second.second, iv.second, 1e-3) << "backscatter sample_sd differs"; } } } // namespace @@ -326,10 +333,10 @@ TEST(ImportEviction, ResurveyedBackscatterBlendsLosslessly) 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. + // And the sample standard deviation (the stored dispersion band) matches too -- + // same sample population => same sample_sd. EXPECT_NEAR(bs_b.at(xcell).second, bs_u.at(xcell).second, 1e-3f) - << "bounded build must reproduce the unbounded intensity variance"; + << "bounded build must reproduce the unbounded backscatter sample_sd"; } // The scratch spill directory is created during eviction and DELETED by finalize. @@ -348,9 +355,7 @@ TEST(ImportEviction, ScratchDirCleanedUpAfterFinalize) 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); + accumulator.finalize(); EXPECT_FALSE(std::filesystem::exists(scratch)) << "finalize must delete the scratch spill dir"; @@ -368,13 +373,13 @@ 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 + // layerDirName(Survey) == "survey": 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)); + marine_bathymetry_store::SourceLayer::Survey)); blocker << "not a directory"; } @@ -400,4 +405,74 @@ TEST(ImportEviction, NeverDropsTileWhenPersistFails) std::filesystem::remove_all(root); } +// Reference-layer seed precedence (#96): a tile seeded from the `reference` prior +// must NOT be counted as measured data. The coarse prior gates blunder rejection +// (seed_settled=false) but never settles a cell, so a sparse survey over a +// reference-seeded tile writes ONLY the resurveyed cells -- identical in count to +// the same survey with NO reference at all, and far fewer than the reference +// tile's own finite cells. +TEST(ImportEviction, ReferenceSeedDoesNotAddMeasuredData) +{ + const std::string root = makeTempDir("refseed"); + const std::string ref_dir = root + "/reference_store"; + + // Build a DENSE reference tile (many finite cells) and write it to the store's + // reference/ layer. Depth ~ 20 m so a same-depth survey later is accepted (not + // blunder-rejected) by the seeded predicted surface. + { + GeoMapSheet ref_sheet(kCellSize); + for (int i = 0; i < 20; ++i) { + // ~1.1 m steps (1e-5 deg) keep the cells within one ~960 m tile. + ref_sheet.addSoundings( + surveyCell(43.0 + i * 1e-5, -70.0 + i * 1e-5, 20.0f, 30.0f)); + } + auto tiles = mapSheetToTiles(ref_sheet); + ASSERT_FALSE(tiles.empty()); + marine_bathymetry_store::BathymetryStore ref_store = + marine_bathymetry_store::BathymetryStore::fromCellSize( + kCellSize, /*reference_writable=*/true); + ref_store.importTiles( + marine_bathymetry_store::SourceLayer::Reference, std::move(tiles)); + marine_bathymetry_store::save(ref_store, ref_dir); + } + const std::size_t reference_finite = countReferenceFiniteCells(ref_dir); + ASSERT_GT(reference_finite, 1u) << "reference tile should cover many cells"; + + // A sparse survey: soundings on ONE cell of that tile, at the same ~20 m depth. + std::vector> batches; + batches.push_back(surveyCell(43.00000, -70.00000, 20.0f, 40.0f)); + + // Run WITH reference seeding. + const std::string with_ref = root + "/with_ref"; + { + GeoMapSheet sheet(kCellSize); + ImportAccumulatorConfig cfg = makeConfig(with_ref, "", /*budget=*/0); + cfg.reference_store_dir = ref_dir; + ImportAccumulator acc(sheet, cfg); + for (const auto & b : batches) { + acc.addBatch(b); + } + acc.finalize(); + } + + // Run WITHOUT reference seeding (baseline). + const std::string no_ref = root + "/no_ref"; + runImport(batches, no_ref, "", /*budget=*/0); + + const auto with = loadBathyCells(with_ref); + const auto without = loadBathyCells(no_ref); + + ASSERT_FALSE(with.empty()) << "the sparse survey should settle at least one cell"; + // The reference seed contributed ZERO measured cells: the two survey stores hold + // the same cell COUNT (only the sparse survey's cells)... + EXPECT_EQ(with.size(), without.size()) + << "reference seeding must not add measured cells"; + // ...and that is far fewer than the reference tile's own finite cells (the prior + // fill was gated-only, never settled). + EXPECT_LT(with.size(), reference_finite) + << "the reference prior's cells must not be settled as survey data"; + + std::filesystem::remove_all(root); +} + } // namespace cube diff --git a/cube_bathymetry/test/test_persistence.cpp b/cube_bathymetry/test/test_persistence.cpp index 4b6c1d5..c4d83c9 100644 --- a/cube_bathymetry/test/test_persistence.cpp +++ b/cube_bathymetry/test/test_persistence.cpp @@ -68,18 +68,15 @@ std::vector makeSoundings() return soundings; } -constexpr int64_t kStamp = 1234567890123456789LL; - // Mirror cube_bathymetry_node::saveDirtyTiles() at library level: write each -// dirty grid as a draft tile under /draft/ (flat, no epoch segment). -std::size_t saveDirty( - GeoMapSheet & sheet, const std::string & dir, int64_t ts_ns) +// dirty grid as a survey tile under /survey/ (flat, no epoch segment). +std::size_t saveDirty(GeoMapSheet & sheet, const std::string & dir) { const std::set dirty = sheet.dirtyGrids(); const std::string out = dir + "/" + marine_bathymetry_store::layerDirName( - marine_bathymetry_store::SourceLayer::Draft); + marine_bathymetry_store::SourceLayer::Survey); std::filesystem::create_directories(out); std::size_t written = 0; for (const auto & index : dirty) { @@ -87,8 +84,7 @@ std::size_t saveDirty( if (!grid) { continue; } - marine_bathymetry_store::BathymetryTile tile = - geoGridToTile(*grid, ts_ns, 0); + marine_bathymetry_store::BathymetryTile tile = geoGridToTile(*grid); if (!tile.dirty()) { continue; } @@ -110,14 +106,14 @@ std::string makeTempDir(const std::string & tag) } } // namespace -// confirm layerDirName(Draft) == "draft" so the hand-built save path matches the -// directory load() scans. -TEST(Persistence, DraftLayerDirNameIsDraft) +// confirm layerDirName(Survey) == "survey" so the hand-built save path matches the +// directory load() scans (uma#248 renamed the live layer draft -> survey). +TEST(Persistence, SurveyLayerDirNameIsSurvey) { EXPECT_EQ( marine_bathymetry_store::layerDirName( - marine_bathymetry_store::SourceLayer::Draft), - "draft"); + marine_bathymetry_store::SourceLayer::Survey), + "survey"); } // Save dirty tiles, then load the whole store back and verify the depth/ @@ -130,12 +126,12 @@ TEST(Persistence, SaveDirtyThenStoreLoadRoundTrips) sheet.addSoundings(makeSoundings()); ASSERT_FALSE(sheet.dirtyGrids().empty()); - const std::size_t written = saveDirty(sheet, dir, kStamp); + const std::size_t written = saveDirty(sheet, dir); ASSERT_GT(written, 0u); EXPECT_TRUE(sheet.dirtyGrids().empty()) << "save must clear the dirty set"; // Reference: convert the same sheet to tiles directly (end-of-session export). - auto reference = mapSheetToTiles(sheet, kStamp, 0); + auto reference = mapSheetToTiles(sheet); ASSERT_FALSE(reference.empty()); // Load the store back through the public store-level load(). @@ -145,7 +141,7 @@ TEST(Persistence, SaveDirtyThenStoreLoadRoundTrips) EXPECT_EQ(loaded_count, written); const auto & draft_tiles = - loaded.tiles(marine_bathymetry_store::SourceLayer::Draft); + loaded.tiles(marine_bathymetry_store::SourceLayer::Survey); ASSERT_FALSE(draft_tiles.empty()) << "draft tiles must be present after load"; // Every finite reference cell must match a loaded cell (depth + uncertainty). @@ -192,10 +188,10 @@ TEST(Persistence, PeriodicSaveEqualsEndOfSessionExport) // and saving wholesale). GeoMapSheet sheet_ref(1.0f); sheet_ref.addSoundings(makeSoundings()); - auto ref_tiles = mapSheetToTiles(sheet_ref, kStamp, 0); + auto ref_tiles = mapSheetToTiles(sheet_ref); // Periodic-style save of every dirty grid. - const std::size_t written = saveDirty(sheet, dir_periodic, kStamp); + const std::size_t written = saveDirty(sheet, dir_periodic); ASSERT_EQ(written, ref_tiles.size()); // Load the periodic store and compare to the reference tiles. @@ -203,7 +199,7 @@ TEST(Persistence, PeriodicSaveEqualsEndOfSessionExport) marine_bathymetry_store::BathymetryStore::fromCellSize(1.0f); marine_bathymetry_store::load(loaded, dir_periodic); const auto & draft_tiles = - loaded.tiles(marine_bathymetry_store::SourceLayer::Draft); + loaded.tiles(marine_bathymetry_store::SourceLayer::Survey); ASSERT_FALSE(draft_tiles.empty()); EXPECT_EQ(draft_tiles.size(), ref_tiles.size()); @@ -232,7 +228,7 @@ std::size_t finiteCellCount( const marine_bathymetry_store::BathymetryStore & store, const gggs::GridIndex & index) { - const auto & tiles = store.tiles(marine_bathymetry_store::SourceLayer::Draft); + const auto & tiles = store.tiles(marine_bathymetry_store::SourceLayer::Survey); auto it = tiles.find(index); if (it == tiles.end()) { return 0; @@ -286,7 +282,7 @@ TEST(Persistence, RevisitAfterEvictPreservesData) // 1. Survey the tile fully and save it (the pre-eviction durable state). GeoMapSheet surveyed(1.0f); surveyed.addSoundings(makeSoundings()); - ASSERT_GT(saveDirty(surveyed, dir, kStamp), 0u); + ASSERT_GT(saveDirty(surveyed, dir), 0u); const std::size_t original_finite = finiteCellCount(loadDraft(dir), index); ASSERT_GT(original_finite, 1u) << "the survey must populate more than one cell"; @@ -295,10 +291,10 @@ TEST(Persistence, RevisitAfterEvictPreservesData) GeoMapSheet reloaded(1.0f); { marine_bathymetry_store::BathymetryStore store = loadDraft(dir); - loadIntoSheet(store, marine_bathymetry_store::SourceLayer::Draft, reloaded); + loadIntoSheet(store, marine_bathymetry_store::SourceLayer::Survey, reloaded); } reloaded.addSoundings(makeSparseResurvey()); // touches ~one cell - ASSERT_GT(saveDirty(reloaded, dir, kStamp + 1), 0u); + ASSERT_GT(saveDirty(reloaded, dir), 0u); const std::size_t after_reload_finite = finiteCellCount(loadDraft(dir), index); EXPECT_GE(after_reload_finite, original_finite) << "reload+resurvey+save must NOT lose any previously-surveyed cell"; @@ -308,10 +304,10 @@ TEST(Persistence, RevisitAfterEvictPreservesData) const std::string dir_ctrl = makeTempDir("lossless_ctrl"); GeoMapSheet seeded(1.0f); seeded.addSoundings(makeSoundings()); - ASSERT_GT(saveDirty(seeded, dir_ctrl, kStamp), 0u); + ASSERT_GT(saveDirty(seeded, dir_ctrl), 0u); GeoMapSheet no_reload(1.0f); no_reload.addSoundings(makeSparseResurvey()); - ASSERT_GT(saveDirty(no_reload, dir_ctrl, kStamp + 1), 0u); + ASSERT_GT(saveDirty(no_reload, dir_ctrl), 0u); const std::size_t no_reload_finite = finiteCellCount(loadDraft(dir_ctrl), index); EXPECT_LT(no_reload_finite, original_finite) << "without reload the overwrite-on-save DOES lose cells -- " diff --git a/cube_bathymetry/test/test_store_import.cpp b/cube_bathymetry/test/test_store_import.cpp index 248e770..eb0998e 100644 --- a/cube_bathymetry/test/test_store_import.cpp +++ b/cube_bathymetry/test/test_store_import.cpp @@ -95,8 +95,6 @@ std::vector makeDeepSoundings(float depth) return soundings; } -constexpr int64_t kStamp = 1234567890123456789LL; -constexpr uint16_t kSource = 7; constexpr float kIntensity = 42.5f; } // namespace @@ -118,7 +116,7 @@ TEST(StoreImport, TileCellsMatchGridValues) // the same way geoGridToTile does, then compare against the produced tile. const std::vector values = grid->values(); const marine_bathymetry_store::BathymetryTile tile = - geoGridToTile(*grid, kStamp, kSource); + geoGridToTile(*grid); gggs::CellAreaIterator it(grid->index()); std::size_t k = 0; @@ -134,8 +132,6 @@ TEST(StoreImport, TileCellsMatchGridValues) ASSERT_TRUE(cell.hasData()); EXPECT_DOUBLE_EQ(cell.depth, static_cast(values[k].depth)); EXPECT_DOUBLE_EQ(cell.uncertainty, static_cast(values[k].uncertainty)); - EXPECT_EQ(cell.timestamp, kStamp); - EXPECT_EQ(cell.source_index, kSource); } } finite_total += finite_in_grid; @@ -144,14 +140,15 @@ TEST(StoreImport, TileCellsMatchGridValues) } // Converting the same map sheet twice must yield an identical tile set: same -// grids, and byte-identical depth/uncertainty/timestamp/source bands. +// grids, and byte-identical depth/uncertainty bands (BathyCell is 2-band since +// uma#248 — the timestamp/source bands were dropped). TEST(StoreImport, ConversionIsDeterministic) { GeoMapSheet ms(1.0f); ms.addSoundings(makeSoundings()); - auto tiles_a = mapSheetToTiles(ms, kStamp, kSource); - auto tiles_b = mapSheetToTiles(ms, kStamp, kSource); + auto tiles_a = mapSheetToTiles(ms); + auto tiles_b = mapSheetToTiles(ms); ASSERT_FALSE(tiles_a.empty()); ASSERT_EQ(tiles_a.size(), tiles_b.size()); @@ -179,8 +176,6 @@ TEST(StoreImport, ConversionIsDeterministic) EXPECT_DOUBLE_EQ(ua[i], ub[i]); } } - EXPECT_EQ(ta.timestampBand(), tb.timestampBand()); - EXPECT_EQ(ta.sourceBand(), tb.sourceBand()); } } @@ -188,7 +183,7 @@ TEST(StoreImport, ConversionIsDeterministic) TEST(StoreImport, EmptyMapSheetYieldsNoTiles) { GeoMapSheet ms(1.0f); - auto tiles = mapSheetToTiles(ms, kStamp, kSource); + auto tiles = mapSheetToTiles(ms); EXPECT_TRUE(tiles.empty()); } @@ -204,7 +199,7 @@ TEST(StoreImport, PrimeFromTileSeedsFiniteCells) const std::vector values = grids.front()->values(); const marine_bathymetry_store::BathymetryTile tile = - geoGridToTile(*grids.front(), kStamp, kSource); + geoGridToTile(*grids.front()); // Prime a fresh sheet from that tile. GeoMapSheet primed(1.0f); @@ -238,7 +233,7 @@ TEST(StoreImport, LoadIntoSheetRoundTrip) GeoMapSheet source(1.0f); source.addSoundings(makeSoundings()); - auto tiles = mapSheetToTiles(source, kStamp, kSource); + auto tiles = mapSheetToTiles(source); ASSERT_FALSE(tiles.empty()); marine_bathymetry_store::BathymetryStore store = @@ -247,11 +242,11 @@ TEST(StoreImport, LoadIntoSheetRoundTrip) // Copy the tile map (importTiles consumes it) but keep a reference set. std::map tiles_copy = tiles; store.importTiles( - marine_bathymetry_store::SourceLayer::Draft, std::move(tiles)); + marine_bathymetry_store::SourceLayer::Survey, std::move(tiles)); GeoMapSheet loaded(1.0f); loadIntoSheet( - store, marine_bathymetry_store::SourceLayer::Draft, loaded); + store, marine_bathymetry_store::SourceLayer::Survey, loaded); std::size_t checked = 0; for (const auto & grid_tile : tiles_copy) { @@ -281,15 +276,15 @@ TEST(StoreImport, LoadIntoSheetEmptyLayerIsNoOp) marine_bathymetry_store::BathymetryStore::fromCellSize(1.0f); GeoMapSheet loaded(1.0f); loadIntoSheet( - store, marine_bathymetry_store::SourceLayer::Draft, loaded); + store, marine_bathymetry_store::SourceLayer::Survey, loaded); EXPECT_TRUE(loaded.grids().empty()); } // Intensity-bearing soundings must produce a non-empty backscatter cell map whose // cells match the grid's finite-intensity nodeRecords() entries cell-for-cell, -// carry the import timestamp/source, and (with a constant input intensity) surface -// that same constant value (the surfaced backscatter is uncorrected by default, -// BackscatterAngleCorrection::None, #80). +// encode the 3-band MbesCell sufficient statistic (mean + n=1 sentinel here), and +// (with a constant input intensity) surface that same constant mean (the surfaced +// backscatter is uncorrected by default, BackscatterAngleCorrection::None, #80). TEST(StoreImport, BackscatterCellsMatchGridRecords) { GeoMapSheet ms(1.0f); @@ -298,7 +293,7 @@ TEST(StoreImport, BackscatterCellsMatchGridRecords) auto grids = ms.grids(); ASSERT_FALSE(grids.empty()); - const auto cells = mapSheetToBackscatterCells(ms, kStamp, kSource); + const auto cells = mapSheetToBackscatterCells(ms); EXPECT_FALSE(cells.empty()) << "intensity-bearing soundings should co-estimate some backscatter cells"; @@ -320,18 +315,28 @@ TEST(StoreImport, BackscatterCellsMatchGridRecords) } else { ++finite_total; ASSERT_NE(found, cells.end()); - EXPECT_FLOAT_EQ(found->second.intensity, records[k].intensity); + // mean is the running-mean band, unchanged from the NodeRecord intensity. + EXPECT_FLOAT_EQ(found->second.mean, records[k].intensity); // Constant input intensity, uncorrected surface by default (None) -> - // emitted value is that constant (mean of equal per-beam intensities). - EXPECT_FLOAT_EQ(found->second.intensity, kIntensity); - EXPECT_EQ(found->second.timestamp, kStamp); - EXPECT_EQ(found->second.source_index, kSource); - // intensity_variance mirrors intensity_var (NaN with < 2 samples); a - // finite value must round-trip exactly. + // emitted mean is that constant (mean of equal per-beam intensities). + EXPECT_FLOAT_EQ(found->second.mean, kIntensity); + // Each cell here has a single intensity-bearing beam (the soundings spread + // one-per-cell), so intensity_var is NaN and the n=1 sentinel is written: + // sample_sd == 0, standard_error == 0. welfordFromCell reconstructs n = 1. if (std::isnan(records[k].intensity_var)) { - EXPECT_TRUE(std::isnan(found->second.intensity_variance)); + EXPECT_FLOAT_EQ(found->second.sample_sd, 0.0f); + EXPECT_FLOAT_EQ(found->second.standard_error, 0.0f); + const IntensityWelford w = welfordFromCell(found->second); + EXPECT_EQ(w.n, 1u); + EXPECT_DOUBLE_EQ(w.mean, static_cast(records[k].intensity)); } else { - EXPECT_FLOAT_EQ(found->second.intensity_variance, records[k].intensity_var); + // n >= 2: sample_sd = sqrt(var*n), standard_error = scale*sqrt(var). + const float n = static_cast(records[k].n_samples); + EXPECT_FLOAT_EQ( + found->second.sample_sd, std::sqrt(records[k].intensity_var * n)); + EXPECT_FLOAT_EQ( + found->second.standard_error, + kBackscatterConfidenceScale * std::sqrt(records[k].intensity_var)); } } } @@ -348,7 +353,7 @@ TEST(StoreImport, BackscatterCellsMatchGridRecords) std::size_t bathy_finite = 0; for (const auto & grid : grids) { const marine_bathymetry_store::BathymetryTile tile = - geoGridToTile(*grid, kStamp, kSource); + geoGridToTile(*grid); gggs::CellAreaIterator it(grid->index()); for (; it.valid(); it.next()) { const marine_bathymetry_store::BathyCell bcell = @@ -379,12 +384,12 @@ TEST(StoreImport, BackscatterNaNPropagation) // Sanity: the same soundings DO yield finite bathy tiles, so an empty // backscatter map is about missing intensity, not missing data. - auto bathy = mapSheetToTiles(ms, kStamp, kSource); + auto bathy = mapSheetToTiles(ms); ASSERT_FALSE(bathy.empty()); GeoMapSheet ms_bs(1.0f); ms_bs.addSoundings(makeSoundings()); - const auto cells = mapSheetToBackscatterCells(ms_bs, kStamp, kSource); + const auto cells = mapSheetToBackscatterCells(ms_bs); EXPECT_TRUE(cells.empty()) << "soundings without intensity must surface no backscatter cells"; } @@ -403,7 +408,7 @@ TEST(StoreImport, PredictedOnlyPrimeSeedsNoSettledHypothesis) const std::vector values = grids.front()->values(); const marine_bathymetry_store::BathymetryTile tile = - geoGridToTile(*grids.front(), kStamp, kSource); + geoGridToTile(*grids.front()); // Predicted-only prime: seeds the predicted surface but creates no hypothesis. GeoMapSheet predicted_only(1.0f); @@ -457,7 +462,7 @@ TEST(StoreImport, SeededPredictedSurfaceRejectsDeepBlunder) // like prior tile we seed from. GeoMapSheet shallow_src(1.0f); shallow_src.addSoundings(makeSoundings()); - auto shallow_tiles = mapSheetToTiles(shallow_src, kStamp, kSource); + auto shallow_tiles = mapSheetToTiles(shallow_src); ASSERT_FALSE(shallow_tiles.empty()); // A clearly-too-deep blunder at the SAME locations (identical touchdown cells), @@ -514,4 +519,92 @@ TEST(StoreImport, SeededPredictedSurfaceRejectsDeepBlunder) EXPECT_GT(gated_cells, 0u) << "the shallow tile should prime some cells to gate"; } +// Backscatter Welford round-trip incl. n>=2 real dispersion (#96): several beams +// with DIFFERENT intensities on one cell produce a finite intensity_var; the +// 3-band encode + welfordFromCell reconstruct the sample count, mean, and estimate +// variance -- the seed-from-store path an off-boat re-run relies on. +TEST(StoreImport, BackscatterWelfordRoundTripMultiSample) +{ + GeoMapSheet ms(1.0f); + // Several soundings at the SAME position (one cell) with distinct intensities so + // the node co-estimates n>=2 with real dispersion (intensity_var finite). + std::vector stacked; + const float intensities[] = {30.0f, 32.0f, 28.0f, 31.0f, 29.0f}; + for (float bs : intensities) { + gz4d::GeoPointLatLongDegrees point(43.07, -70.76, -12.0); + GeoSounding s(point); + s.sounding.vertical_error = 0.5f; + s.sounding.horizontal_error = 0.1f; + s.sounding.intensity = bs; + s.sounding.beam_angle = 0.0f; + stacked.push_back(s); + } + ms.addSoundings(stacked); + + const auto cells = mapSheetToBackscatterCells(ms); + ASSERT_FALSE(cells.empty()); + + std::size_t checked_multi = 0; + for (const auto & grid : ms.grids()) { + const std::vector records = grid->nodeRecords(); + gggs::CellAreaIterator it(grid->index()); + std::size_t k = 0; + for (; it.valid() && k < records.size(); it.next(), ++k) { + if (std::isnan(records[k].intensity)) {continue;} + const auto found = cells.find(*it); + ASSERT_NE(found, cells.end()); + if (!std::isnan(records[k].intensity_var)) { + ++checked_multi; + const IntensityWelford w = welfordFromCell(found->second); + EXPECT_EQ(w.n, records[k].n_samples); + EXPECT_NEAR(w.mean, static_cast(records[k].intensity), 1e-4); + // Estimate variance ((m2/(n-1))/n) round-trips to intensity_var. + const double recon_var = (w.m2 / (w.n - 1)) / w.n; + EXPECT_NEAR( + recon_var, static_cast(records[k].intensity_var), 1e-4); + } + } + } + EXPECT_GT(checked_multi, 0u) + << "stacked soundings should co-estimate at least one n>=2 cell"; +} + +// welfordFromCell unit round-trip: the n=1 sentinel, the no-data cell, and the +// n>=2 path each invert the 3-band encode exactly (#96). +TEST(StoreImport, WelfordFromCellInvertsEncode) +{ + namespace mbs = marine_mbes_backscatter_store; + // n = 1 sentinel: sample_sd == 0 with a finite mean -> n = 1, M2 = 0. + { + mbs::MbesCell cell; + cell.mean = 25.0f; + cell.standard_error = 0.0f; + cell.sample_sd = 0.0f; + const IntensityWelford w = welfordFromCell(cell); + EXPECT_EQ(w.n, 1u); + EXPECT_DOUBLE_EQ(w.mean, 25.0); + EXPECT_DOUBLE_EQ(w.m2, 0.0); + } + // no-data cell (mean NaN) -> empty Welford {n = 0}. + { + const mbs::MbesCell cell; // defaults: mean NaN + const IntensityWelford w = welfordFromCell(cell); + EXPECT_EQ(w.n, 0u); + } + // n >= 2: build the 3-band from a known (n, var), invert, confirm n recovered. + { + const uint32_t n = 7; + const float var = 0.5f; // variance of the mean + mbs::MbesCell cell; + cell.mean = 40.0f; + cell.sample_sd = std::sqrt(var * n); + cell.standard_error = kBackscatterConfidenceScale * std::sqrt(var); + const IntensityWelford w = welfordFromCell(cell); + EXPECT_EQ(w.n, n); + EXPECT_DOUBLE_EQ(w.mean, 40.0); + const double recon_var = (w.m2 / (w.n - 1)) / w.n; + EXPECT_NEAR(recon_var, static_cast(var), 1e-5); + } +} + } // namespace cube diff --git a/cube_bathymetry/test/test_tile_eviction_rss.cpp b/cube_bathymetry/test/test_tile_eviction_rss.cpp index a46e179..c9b1349 100644 --- a/cube_bathymetry/test/test_tile_eviction_rss.cpp +++ b/cube_bathymetry/test/test_tile_eviction_rss.cpp @@ -52,8 +52,6 @@ namespace cube namespace { -constexpr int64_t kStamp = 1234567890123456789LL; - // Add several soundings at one location so CUBE settles a finite estimate (a // single sample may not resolve) -- one tile's worth of survey. void surveyTile(GeoMapSheet & sheet, double lat, double lon) @@ -79,7 +77,7 @@ bool persistTile( if (!grid) { return false; } - marine_bathymetry_store::BathymetryTile tile = geoGridToTile(*grid, kStamp, 0); + marine_bathymetry_store::BathymetryTile tile = geoGridToTile(*grid); if (!tile.dirty()) { return false; } @@ -125,7 +123,7 @@ TEST(TileEvictionRss, ResidentCountBoundedWhileDiskGrows) const std::string draft_dir = dir + "/" + marine_bathymetry_store::layerDirName( - marine_bathymetry_store::SourceLayer::Draft); + marine_bathymetry_store::SourceLayer::Survey); std::filesystem::create_directories(draft_dir); GeoMapSheet sheet(1.0f);