docs(reviewer): add external audits — test/error-handling, API/perf, CGAL readiness

Three complementary, self-contained audit documents (each actionable by a fresh session, with file:line, code snippets, fixes, acceptance criteria): - test-coverage-error-handling-audit: coverage gaps + error-handling robustness; agreed gate 80% line / 70% branch / 90% function. - api-performance-audit: API-naming consistency (A1–A5, decided) + Newton-solver performance (B1 block-FD, redundant gradient evals, factorization reuse). - cgal-submission-readiness-audit: G0 porting-rights blocker (original is unlicensed Varylab/TU-Berlin code; this is a derivative work — clarify with authors before any license/release), G1 license, layout, concept, manual, tests. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
2026-05-31 09:15:39 +02:00
parent d725166c6a
commit 7902ae8e2c
3 changed files with 1382 additions and 0 deletions
--- a/doc/reviewer/api-performance-audit-2026-05-31.md
+++ b/doc/reviewer/api-performance-audit-2026-05-31.md
@@ -0,0 +1,420 @@
+# API-Consistency & Performance Audit — ConformalLabpp
+
+**Date:** 2026-05-31
+**Auditor:** External reviewer (Claude Opus 4.8)
+**Scope:** Low-level API (`code/include/*.hpp`), high-level CGAL API (`code/include/CGAL/`), Newton-solver performance
+**Focus:** Public-API naming consistency and runtime performance — NOT port-faithfulness or test coverage.
+
+This document is self-contained. A new session can pick up any finding below and
+act on it without prior context. Each finding includes:
+- exact file path + line numbers (verified by direct file read)
+- a minimal reproduction of the problematic code
+- the recommended fix
+- acceptance criteria for "done"
+
+Status legend: 🔴 Critical · 🟡 Important · 🔵 Hint / nice-to-have
+
+> **Companion documents (no overlap in findings):**
+> - `external-audit-2026-05-30.md` — port-faithfulness bugs
+> - `test-coverage-error-handling-audit-2026-05-31.md` — test gaps + error handling
+> - **this file** — API consistency + performance
+>
+> Cross-reference: performance findings B2/B3/B5 below share root cause with
+> `FINDING-H2` ("five near-identical Newton loops") in the test-coverage audit.
+> A single `newton_core` refactor can resolve all of them together.
+
+---
+
+## How to read this document in a new session
+
+```bash
+# Build the CGAL suite (the API + solvers live behind WITH_CGAL_TESTS)
+cmake -S code -B build-cgal -DWITH_CGAL_TESTS=ON
+cmake --build build-cgal --target conformallab_cgal_tests -j$(nproc)
+
+# Baseline before any change
+ctest --test-dir build-cgal -R '^cgal\.' --output-on-failure
+# Expected: 246 passed, 0 failed
+```
+
+For performance findings, the smoke meshes used by the project are:
+- `cathead.obj` (open, F≈248)
+- `brezel.obj` (genus-1, F≈13824)
+- `brezel2.obj` (genus-2)
+exercised by `code/tests/cgal/test_scalability_smoke.cpp`.
+
+---
+
+## Summary table
+
+| ID | Sev | Title | Primary location |
+|----|-----|-------|------------------|
+| B1 | 🔴 | `newton_hyper_ideal` uses full-FD Hessian; the 96–1166× faster block-FD exists but is unused. Inversive-Distance has no fast path at all | `newton_solver.hpp:433`, `:580` |
+| A1 | 🔴 | `assign_*_dof_indices`: three different naming conventions across the 5 functionals | 5 functional headers |
+| B2 | 🟡 | `line_search` discards the gradient at the accepted point → ~1 redundant full gradient eval per Newton iteration (×5 solvers) | `newton_solver.hpp` |
+| B3 | 🟡 | `has_edge_dof` O(E) scan runs every Newton iteration (layout is loop-invariant) | `newton_solver.hpp:264` |
+| A4 | 🟡 | High-level CGAL entry points have inconsistent suffix (`_map` only on one) | `code/include/CGAL/*.h` |
+| A5 | 🟡 | The two circle-packing models return different result types | `Discrete_circle_packing.h`, `Discrete_inversive_distance.h` |
+| A2 | 🟡 | `compute_*_from_mesh`: inconsistent prefix + verb (`lambda0` vs `init`) | 3 functional headers |
+| A3 | 🟡 | `gradient_check` (HyperIdeal) lacks the `_<geom>` suffix the other 4 have | `hyper_ideal_functional.hpp:470` |
+| B4 | 🟡 | SimplicialLDLT symbolic factorization rebuilt every iteration (sparsity pattern is constant) | `newton_solver.hpp:73` |
+| B5 | 🔵 | Dead `SimplicialLDLT solver;` variable + per-call gradient heap allocation | `newton_solver.hpp:244`, `euclidean_functional.hpp:193` |
+
+---
+
+# Part A — API consistency
+
+> Context: the project has two API layers. The **high-level CGAL layer**
+> (`CGAL::discrete_*` + named parameters) intentionally follows CGAL conventions —
+> that is correct and not a finding by itself. The **low-level free-function layer**
+> (`setup_*`, `assign_*`, `*_gradient`) is the one used directly by tests and
+> internal callers, and it is internally inconsistent. `setup_<geom>_maps` IS
+> uniform across all five functionals, which proves the convention is achievable —
+> it simply was not held for the other verbs.
+
+## FINDING-A1 — 🔴 `assign_*_dof_indices`: three conventions for one operation
+
+### Location
+- `euclidean_functional.hpp:107,132`
+- `spherical_functional.hpp:99,123`
+- `hyper_ideal_functional.hpp:107`
+- `cp_euclidean_functional.hpp:140`
+- `inversive_distance_functional.hpp:133`
+
+### Problem
+| Functional | Vertex variant | All variant |
+|---|---|---|
+| Euclidean | `assign_euclidean_vertex_dof_indices` | `assign_euclidean_all_dof_indices` |
+| Spherical | `assign_vertex_dof_indices` ⚠️ no prefix | `assign_all_spherical_dof_indices` ⚠️ infix |
+| HyperIdeal | — | `assign_all_dof_indices` ⚠️ no prefix |
+| CP-Euclidean | `assign_cp_euclidean_face_dof_indices` | — |
+| InversiveDist | `assign_inversive_distance_vertex_dof_indices` | — |
+
+Three patterns: `assign_<geom>_<scope>_dof_indices`, `assign_<scope>_dof_indices`
+(no geometry), and `assign_all_<geom>_dof_indices` (geometry in the **middle**).
+The prefix-less `assign_vertex_dof_indices` (spherical) and `assign_all_dof_indices`
+(hyper-ideal) read as generic but are geometry-specific, which is misleading at
+the call site (e.g. `test_newton_solver.cpp:81` calls `assign_vertex_dof_indices`).
+
+### Fix — DECIDED (2026-05-31)
+Standardize on `assign_<geom>_<scope>_dof_indices` (matches the already-consistent
+`setup_<geom>_maps`). Confirmed renames:
+
+| Old | New |
+|---|---|
+| `assign_vertex_dof_indices` (spherical) | `assign_spherical_vertex_dof_indices` |
+| `assign_all_spherical_dof_indices` | `assign_spherical_all_dof_indices` |
+| `assign_all_dof_indices` (hyper-ideal) | `assign_hyper_ideal_all_dof_indices` |
+| `assign_euclidean_vertex_dof_indices` | unchanged ✓ |
+| `assign_euclidean_all_dof_indices` | unchanged ✓ |
+| `assign_cp_euclidean_face_dof_indices` | unchanged ✓ |
+| `assign_inversive_distance_vertex_dof_indices` | unchanged ✓ |
+
+Keep the old names as `[[deprecated]]` inline wrappers for one release so existing
+call sites (tests, examples) keep compiling, then update all call sites.
+
+### Acceptance criteria
+- All five functionals follow `assign_<geom>_<scope>_dof_indices`.
+- Old names exist as `[[deprecated]]` aliases (or all call sites updated in the same PR).
+- Full CGAL suite still 246/246.
+
+---
+
+## FINDING-A2 — 🟡 `compute_*_from_mesh`: inconsistent prefix and verb
+
+### Location
+- `euclidean_functional.hpp:152` — `compute_euclidean_lambda0_from_mesh`
+- `spherical_functional.hpp:148` — `compute_lambda0_from_mesh` ⚠️ no prefix
+- `inversive_distance_functional.hpp:187` — `compute_inversive_distance_init_from_mesh` ⚠️ `init` not `lambda0`
+
+### Problem
+`compute_lambda0_from_mesh` (spherical) has no geometry prefix and is used directly
+in `test_newton_solver.cpp:80,99`. Inversive-Distance uses the verb `init` instead
+of `lambda0` for the conceptually parallel step.
+
+### Fix — DECIDED (2026-05-31)
+- `compute_lambda0_from_mesh` → `compute_spherical_lambda0_from_mesh` (add prefix).
+- `compute_inversive_distance_init_from_mesh` → **keep `init`**. Rationale: this step
+  initializes `I_e` and `r0`, i.e. genuinely more than just λ₀, so `init` is the
+  honest verb. The distinction (`lambda0` vs `init`) will be documented in
+  `doc/api/extending.md`.
+- `compute_euclidean_lambda0_from_mesh` → unchanged ✓.
+- `[[deprecated]]` alias for the spherical rename.
+
+### Acceptance criteria
+- Spherical compute fn carries the `spherical` prefix.
+- Naming rationale for `init` vs `lambda0` documented in `doc/api/extending.md`.
+
+---
+
+## FINDING-A3 — 🟡 `gradient_check` (HyperIdeal) breaks the `_<geom>` suffix pattern
+
+### Location
+- `euclidean_functional.hpp:324` — `gradient_check_euclidean`
+- `spherical_functional.hpp:367` — `gradient_check_spherical`
+- `cp_euclidean_functional.hpp:331` — `gradient_check_cp_euclidean`
+- `inversive_distance_functional.hpp:346` — `gradient_check_inversive_distance`
+- `hyper_ideal_functional.hpp:470` — `gradient_check` ⚠️ no suffix
+
+### Fix — DECIDED (2026-05-31)
+Rename `gradient_check` → `gradient_check_hyper_ideal` with a `[[deprecated]]` alias.
+The other four (`gradient_check_euclidean/_spherical/_cp_euclidean/_inversive_distance`)
+are unchanged ✓.
+
+### Acceptance criteria
+- All five follow `gradient_check_<geom>`.
+
+---
+
+## FINDING-A4 — 🟡 High-level CGAL entry points: inconsistent suffix
+
+### Location
+`code/include/CGAL/Discrete_conformal_map.h`, `Discrete_circle_packing.h`,
+`Discrete_inversive_distance.h`
+
+### Problem
+```cpp
+CGAL::discrete_conformal_map_euclidean       // no trailing word
+CGAL::discrete_conformal_map_spherical
+CGAL::discrete_conformal_map_hyper_ideal
+CGAL::discrete_circle_packing_euclidean      // no _map
+CGAL::discrete_inversive_distance_map        // ⚠️ trailing _map
+```
+`discrete_inversive_distance_map` ends in `_map`; `discrete_circle_packing_euclidean`
+does not — yet both are circle-packing models.
+
+### Fix — DECIDED (2026-05-31)
+`discrete_inversive_distance_map` → **`discrete_inversive_distance_euclidean`** (drop
+the trailing `_map`, add the geometry word), so it matches
+`discrete_circle_packing_euclidean`. Because this is the **public CGAL API**, treat
+it as a breaking change: add the new name, mark the old `[[deprecated]]`, announce in
+`CHANGELOG.md`.
+
+> ⚠️ Sequencing: A4 touches the **public CGAL surface**. Do this rename only after the
+> G0/G1 license/provenance outcome is known (see CGAL audit) — there is no point
+> stabilizing a public API on a package whose release status is unresolved. A1–A3 are
+> internal and can proceed now.
+
+### Acceptance criteria
+- The two circle-packing entry points share a naming shape
+  (`discrete_*_euclidean`).
+- CHANGELOG documents the rename; deprecated alias provided.
+
+---
+
+## FINDING-A5 — 🟡 The two circle-packing models return different result types
+
+### Location
+- `Discrete_circle_packing.h:100` → returns `Circle_packing_result<FT>`
+- `Discrete_inversive_distance.h:133,146` → returns `Conformal_map_result<FT>`
+
+### Problem
+CP-Euclidean and Inversive-Distance are both circle-packing models, but
+`discrete_circle_packing_euclidean` returns `Circle_packing_result` while
+`discrete_inversive_distance_map` returns `Conformal_map_result`. A user switching
+between the two packers must rewrite their result-handling code, even though the
+semantics (per-element radii + Newton diagnostics) are the same.
+
+### Fix — DECIDED (2026-05-31)
+**Option (a):** both circle-packing entry points return `Circle_packing_result`.
+`discrete_inversive_distance_*` switches from `Conformal_map_result` to
+`Circle_packing_result` (semantically truthful — both produce circle radii). The
+unified-result option (b) is rejected as over-engineered.
+
+> ⚠️ Same sequencing as A4 — public-API change, do after the G0/G1 outcome.
+
+### Acceptance criteria
+- The two packing entry points return `Circle_packing_result`.
+- If a field is renamed (e.g. `u_per_vertex` → `radius_per_*`), document it in CHANGELOG.
+
+---
+
+# Part B — Performance
+
+## FINDING-B1 — 🔴 Fast Hessian paths exist but the solvers don't use them
+
+### Location
+- `newton_solver.hpp:433` — `newton_hyper_ideal` calls `hyper_ideal_hessian_sym(...)`
+  (the **full-FD** O(n·F) version)
+- `hyper_ideal_hessian.hpp:220` — `hyper_ideal_hessian_block_fd_sym(...)` (the fast
+  block-FD version) **exists and is tested** but is never called by the solver
+- `newton_solver.hpp:580-608` — `newton_inversive_distance` builds a **full-FD**
+  Hessian inline (2 gradient evals per DOF), with no fast path at all
+
+### Problem
+**HyperIdeal:** the solver uses the slow full-FD Hessian even though a block-FD
+version with a documented **33× speedup on cathead, ~1166× on brezel** already
+exists (`hyper_ideal_hessian.hpp:130-136`). The fast code is implemented, tested
+(`test_hyper_ideal_hessian.cpp`), and simply not wired into the solver.
+
+```cpp
+// newton_solver.hpp:433 — current (slow):
+auto H  = hyper_ideal_hessian_sym(mesh, x, m, hess_eps);          // full-FD, O(n·F)
+// available (fast, same result up to FD noise):
+auto H  = hyper_ideal_hessian_block_fd_sym(mesh, x, m, hess_eps); // block-FD
+```
+
+**Inversive-Distance:** the inline `build_hessian` lambda runs `2n` full gradient
+evaluations per Hessian (each O(F)), i.e. **O(n·F) per Newton iteration** — quadratic
+in mesh size. No block-FD or analytic path exists yet.
+
+```cpp
+// newton_solver.hpp:585-601 — current:
+for (int j = 0; j < n; ++j) {
+    auto Gp = inversive_distance_gradient(mesh, xp, m);   // O(F)
+    auto Gm = inversive_distance_gradient(mesh, xm, m);   // O(F)  → 2n total per Hessian
+    ...
+}
+```
+
+### Fix
+1. **HyperIdeal (cheap win):** switch line 433 to `hyper_ideal_hessian_block_fd_sym`.
+   Verify the existing convergence tests (`test_newton_solver.cpp` HyperIdeal cases,
+   `test_newton_phase9a.cpp`) still pass within tolerance — block-FD is mathematically
+   equivalent up to FD rounding (locality lemma, documented at `hyper_ideal_hessian.hpp:120-128`).
+2. **Inversive-Distance:** port a block-FD Hessian mirroring the HyperIdeal one
+   (the gradient is also face-decomposable). Track the analytic Glickenstein-2011
+   Hessian separately as already noted in `doc/roadmap/research-track.md`.
+
+### Acceptance criteria
+- `newton_hyper_ideal` uses block-FD; HyperIdeal convergence tests still pass.
+- A timing assertion or smoke test demonstrates the speedup on `brezel.obj`
+  (or at minimum, the change is benchmarked and recorded in `doc/architecture/`).
+- Inversive-Distance has a block-FD path; `test_newton_phase9a.cpp` still passes.
+
+---
+
+## FINDING-B2 — 🟡 line_search discards the gradient at the accepted point
+
+### Location
+`newton_solver.hpp` — `detail::line_search` (lines 140-201) and all 5 solver loops
+(e.g. Euclidean lines 248, 277-280, 287)
+
+### Problem
+Per Newton iteration the gradient over the whole mesh is computed redundantly:
+1. line 248 — gradient for the convergence check
+2. inside `line_search` — gradient at each trial point; the **vector is discarded**,
+   only its norm is kept
+3. next iteration's line 248 — recomputes the gradient at the just-accepted point
+   (already computed in step 2)
+4. line 287 (`G_final`) — recomputes once more after the loop
+
+That is roughly **one redundant full-mesh gradient per iteration** plus one at the end.
+
+### Fix
+Have `line_search` return (via out-param) the gradient vector at the accepted point,
+and feed it into the next iteration's convergence check instead of recomputing. The
+end-of-loop `G_final` can reuse the last accepted gradient too.
+
+> Best done as part of the `newton_core` template refactor (test-coverage audit
+> FINDING-H2) so the fix lands once instead of five times.
+
+### Acceptance criteria
+- Gradient evaluations per iteration drop by ~1 (instrument with a counter in a test).
+- All convergence tests pass unchanged.
+
+---
+
+## FINDING-B3 — 🟡 has_edge_dof scan runs every Newton iteration
+
+### Location
+`newton_solver.hpp:264-268` (inside the `newton_euclidean` iteration loop)
+
+### Problem
+```cpp
+bool has_edge_dof = false;
+for (auto e : mesh.edges())          // O(E), runs EVERY iteration
+    if (m.e_idx[e] >= 0) { has_edge_dof = true; break; }
+```
+The DOF layout (`m.e_idx`) does not change during the solve, so this scan is
+loop-invariant. Over up to 200 iterations the mesh edge set is scanned 200× for no
+reason.
+
+### Fix
+Hoist the scan above the `for (iter...)` loop; compute `has_edge_dof` once.
+
+### Acceptance criteria
+- `has_edge_dof` computed exactly once per `newton_euclidean` call.
+- Euclidean tests pass unchanged.
+
+---
+
+## FINDING-B4 — 🟡 Symbolic factorization rebuilt every iteration
+
+### Location
+`newton_solver.hpp:73` (`detail::solve_with_fallback`)
+
+### Problem
+```cpp
+Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt(A);  // analyze + factorize
+```
+A fresh `SimplicialLDLT` is constructed each iteration, which re-runs the symbolic
+analysis (fill-reducing reordering / COLAMD). But the **sparsity pattern of the
+Hessian is constant across iterations** (only the numeric values change). Eigen
+supports `analyzePattern()` once + `factorize()` per iteration, skipping the
+repeated symbolic step.
+
+### Fix
+This needs a structural change: `solve_with_fallback` is stateless, so it cannot
+cache the analysis. Either:
+- thread a persistent solver object through the Newton loop (cleanest with the
+  `newton_core` refactor), calling `analyzePattern` on the first iteration and
+  `factorize` thereafter, or
+- key a cached symbolic factorization on the matrix sparsity pattern.
+
+Keep the SparseQR fallback for the singular-pattern case (the pattern there is
+also constant, so the same caching applies).
+
+### Acceptance criteria
+- Symbolic analysis runs once per solve, not once per iteration.
+- All Newton + fallback tests (incl. the 3 SparseQR tests) pass.
+- Benchmarked on `brezel.obj`; improvement recorded.
+
+---
+
+## FINDING-B5 — 🔵 Dead variable + per-call gradient allocation
+
+### Location
+- `newton_solver.hpp:244` — `Eigen::SimplicialLDLT<...> solver;` declared in
+  `newton_euclidean`, **never used** (the solve goes through `solve_with_fallback`)
+- `euclidean_functional.hpp:193,197` — each `euclidean_gradient` call freshly
+  allocates `std::vector<double> G(n)` and `std::vector<double> h_alpha(nh)`
+
+### Problem
+- The dead `solver` variable is harmless but misleading (suggests reuse that
+  doesn't happen — directly relevant once B4 is addressed).
+- In the FD Hessian (B1) the per-call allocation means `2n × 2` fresh heap
+  allocations per Hessian build.
+
+### Fix
+- Remove the dead `solver` declaration.
+- When addressing B1/B2, pass reusable scratch buffers (`G`, `h_alpha`) into the
+  gradient functions, or provide an overload that writes into a caller-owned buffer.
+
+### Acceptance criteria
+- Dead variable removed (compiles clean with `-Wunused`).
+- (Optional, with B1/B2) gradient functions support buffer reuse; allocation count
+  in the FD Hessian path drops measurably.
+
+---
+
+## Suggested order of work
+
+1. **B1 (HyperIdeal half)** — one-line switch to `hyper_ideal_hessian_block_fd_sym`;
+   biggest win for the smallest change. Verify tests, benchmark.
+2. **B3** — trivial loop hoist.
+3. **A1, A2, A3** — low-level renames with `[[deprecated]]` aliases; do together.
+4. **A4, A5** — high-level CGAL API politur; coordinate with CHANGELOG (breaking).
+5. **newton_core refactor** (test-coverage audit FINDING-H2) — then fold in **B2,
+   B4, B5** in the same pass.
+6. **B1 (Inversive-Distance half)** — port a block-FD Hessian; larger effort.
+
+## What is already good (do not "fix")
+
+- `setup_<geom>_maps` is uniform across all five functionals — the consistency
+  target is proven achievable.
+- The block-FD HyperIdeal Hessian is correct, tested, and well-documented — it
+  just needs to be wired into the solver (B1).
+- The high-level CGAL API correctly uses CGAL named parameters and `Kernel_traits` —
+  idiomatic for a CGAL package.
+- `Conformal_map_result` carries `sparse_qr_fallback_used` — good observability.