ConformalLabpp/doc/reviewer/api-performance-audit-2026-05-31.md

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

# 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

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.

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

// 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

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

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.