# 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.