diff --git a/doc/reviewer/dependency-license-audit-2026-05-31.md b/doc/reviewer/dependency-license-audit-2026-05-31.md new file mode 100644 index 0000000..42a2458 --- /dev/null +++ b/doc/reviewer/dependency-license-audit-2026-05-31.md @@ -0,0 +1,167 @@ +# Dependency & License-Compatibility Audit β€” ConformalLabpp + +**Date:** 2026-05-31 +**Auditor:** External reviewer (Claude Opus 4.8) +**Scope:** Vendored/auto-fetched dependencies, their license compatibility, and the +provenance/license of bundled **data assets** (test meshes). +**Focus:** Is the dependency + data licensing coherent β€” especially in light of the +unresolved provenance/license blocker (CGAL audit **G0/G1**)? + +Status legend: πŸ”΄ Critical Β· 🟑 Important Β· πŸ”΅ Polish + +> **Good news up front:** a `code/deps/THIRD-PARTY-LICENSES.md` already exists and is +> a genuinely good per-dependency compatibility matrix (CGAL, Eigen, libigl, GLFW, +> nlohmann/json). This audit does **not** redo that work. It flags two things that +> matrix cannot see: (D1) it is built on the MIT premise that G0 undermines, and +> (D2) the bundled **test meshes have no provenance or license at all**. + +--- + +## Summary table + +| ID | Sev | Title | Location | +|----|-----|-------|----------| +| D1 | πŸ”΄ | The dependency-license matrix is predicated on "conformallab++ is MIT" β€” a premise CGAL-audit G0 shows is unfounded; the whole compatibility conclusion is conditional on the unresolved license decision | `code/deps/THIRD-PARTY-LICENSES.md` | +| D2 | πŸ”΄ | Bundled test meshes (`cathead.obj`, `brezel*.obj`, torus OFFs) have **no provenance and no license** β€” likely third-party / Varylab-origin, redistributed without attribution | `code/data/` | +| D3 | 🟑 | No per-dependency license *files* under most vendored trees (only GLFW ships one); the matrix is the only record | `code/deps/*/` | +| D4 | πŸ”΅ | Vendored-tarball integrity is asserted ("bit-for-bit identical") but not checksum-verified in-repo | `code/deps/tarballs/` | + +--- + +## D1 β€” πŸ”΄ License matrix is conditional on the unresolved MIT premise + +### Evidence +`code/deps/THIRD-PARTY-LICENSES.md` opens its analysis with: +> *"## conformallab++ itself … **MIT** … Every C++ source file carries +> `SPDX-License-Identifier: MIT`"* + +and every dependency row evaluates "**Compatibility with MIT distribution**". + +### Problem +The matrix's central column β€” and therefore its conclusion that everything is +distributable β€” assumes conformallab++ *is* MIT-licensed and *may* be distributed as +MIT. CGAL-audit **G0** establishes that this premise is unfounded: the code is a +documented port (derivative work) of the unlicensed Varylab/TU-Berlin project, so the +MIT grant itself is on shaky ground. Two consequences cascade: + +1. If porting rights are **not** granted, there is no valid MIT (or any) license to + be "compatible with" β€” the matrix evaluates compatibility against a license the + project may not be entitled to use. +2. If the project relicenses to **GPLv3+** for CGAL submission (G1), the entire + "compatible with MIT distribution" analysis must be **redone against GPLv3+** + (e.g. the CGAL LGPL/GPL split, Eigen's MPL-2.0, libigl's MPL-2.0 are all + GPL-compatible β€” but that needs stating, not assuming). + +### Fix +- Resolve G0/G1 first (owner-level). +- Then re-run the matrix against the *actual* outbound license. Add a column + "compatibility with GPLv3+ distribution" if the CGAL path is taken. +- Until resolved, annotate `THIRD-PARTY-LICENSES.md` with a banner that its + conclusions are conditional on the pending G0/G1 decision. + +### Acceptance criteria +- The matrix's outbound-license assumption matches the resolved G0/G1 outcome and is + re-evaluated against it. + +--- + +## D2 β€” πŸ”΄ Bundled test meshes have no provenance or license + +### Evidence +``` +code/data/obj: brezel.obj brezel2.obj cathead.obj tetraflat.obj +code/data/off: simple_cupe.off torus_4x4.off torus_8x8.off torus_hex_6x6.off torus_skewed_4x4.off +``` +`find code/data -iname '*readme*' -o -iname '*licen*' -o -iname '*source*'` β†’ **nothing.** + +### Problem +These meshes are committed and redistributed with the project, used by the smoke and +parity tests (`test_scalability_smoke.cpp`, `test_geometry_utils.cpp` references +`brezel2.obj`, etc.). But: + +- `cathead.obj` is a **well-known third-party mesh** that circulates in graphics + courses/datasets β€” it is almost certainly not original to this project and carries + whoever's terms. +- `brezel.obj` / `brezel2.obj` ("Brezel" = pretzel; genus-1/genus-2) are the exact + example surfaces used in the **Varylab/TU-Berlin** discrete-conformal work β€” i.e. + likely the *same* upstream as the G0 code-provenance problem. + +Redistributing third-party data assets without attribution or a license is the same +class of issue as the code port (G0), and CGAL submission would flag it: CGAL example +data must have clear redistribution rights. + +### Fix +- Establish each mesh's origin and license. For `cathead.obj`, identify the standard + source and its terms. For the `brezel*` meshes, this is part of the **same + conversation with the Varylab/TU-Berlin authors** as G0. +- Add `code/data/PROVENANCE.md` listing each file: source, author, license, URL. +- Replace any mesh whose redistribution rights cannot be confirmed with a + cleanly-licensed or self-generated equivalent (the torus OFFs look procedurally + generated β€” confirm and document that, which would make them safe). + +### Acceptance criteria +- Every file under `code/data/` has documented provenance + redistribution rights; + unconfirmed assets are replaced. + +--- + +## D3 β€” 🟑 Most vendored trees ship no license file + +### Evidence +``` +eigen-3.4.0: NONE FOUND CGAL-6.1.1: NONE FOUND +libigl-2.6.0: NONE FOUND glfw-3.4: LICENSE.md βœ“ +``` +Only GLFW retains its upstream license file; the others rely solely on the central +`THIRD-PARTY-LICENSES.md`. + +### Problem +When vendoring source trees, the upstream LICENSE/COPYING file should travel with the +code (it is usually a license *requirement* β€” MPL-2.0 and LGPL both require preserving +license notices). A central summary is good practice but does not substitute for the +upstream notice files inside each tree. + +### Fix +Restore each upstream license file into its vendored tree (`eigen-3.4.0/COPYING.*`, +`CGAL-6.1.1/LICENSE*`, `libigl-2.6.0/LICENSE*`). They were likely stripped during +trimming of the vendored copies. + +### Acceptance criteria +- Each vendored dependency tree contains its upstream license notice file. + +--- + +## D4 β€” πŸ”΅ Vendored tarball integrity not checksum-verified + +### Evidence +`THIRD-PARTY-LICENSES.md` states the cached tarballs are "bit-for-bit identical to the +upstream releases", but no checksum manifest is committed. + +### Problem +The reproducibility claim (the stated rationale for vendoring) rests on an unverified +assertion. A corrupted or substituted tarball would not be detected. + +### Fix +Commit a `code/deps/tarballs/SHA256SUMS` and verify it in the extraction step of the +CMake/deps logic. + +--- + +## What is already good + +- `THIRD-PARTY-LICENSES.md` is a thorough, honest per-dependency matrix with real + compatibility reasoning (the CGAL LGPL/GPL split, Eigen's NonMPL2 gating, the + viewer-only scope of libigl/GLFW) β€” excellent practice, just built on a premise + (D1) that G0 destabilizes. +- The dependency choices themselves are clean: Eigen (MPL-2.0), libigl (MPL-2.0), + GLFW (zlib), nlohmann/json (MIT) are all permissive; only CGAL carries copyleft, and + its header-only LGPL consumption is correctly reasoned. +- Vendoring at fixed versions for build reproducibility is a defensible choice and is + documented. + +## Suggested order +1. **D2** (data provenance) β€” couple with the G0 conversation; CGAL-blocking and the + most clearly missing piece. +2. **D1** (re-base the matrix on the resolved license) β€” after G0/G1. +3. **D3** (restore upstream license files) β€” quick, and likely a license requirement. +4. **D4** (checksum manifest) β€” polish. diff --git a/doc/reviewer/input-validation-audit-2026-05-31.md b/doc/reviewer/input-validation-audit-2026-05-31.md new file mode 100644 index 0000000..e3d63a7 --- /dev/null +++ b/doc/reviewer/input-validation-audit-2026-05-31.md @@ -0,0 +1,197 @@ +# Input-Validation & Malformed-Input Audit β€” ConformalLabpp + +**Date:** 2026-05-31 +**Auditor:** External reviewer (Claude Opus 4.8) +**Scope:** Robustness of the library's input boundaries β€” mesh I/O (OFF/OBJ/PLY) and +result (de)serialization (JSON/XML) β€” against malformed, adversarial, or +out-of-domain input. +**Focus:** What happens with *bad* input. Distinct from the error-handling audit, +which covered *internal* error propagation and the deliberate `throw` paths. + +Status legend: πŸ”΄ Critical Β· 🟑 Important Β· πŸ”΅ Polish + +> **Threat model:** this is a scientific library, not a network service, so the bar +> is "fail cleanly and diagnosably", not "resist attackers". But meshes and result +> files routinely come from other tools, other people, and older versions of this +> code β€” so malformed input is a *when*, not an *if*. + +--- + +## Summary table + +| ID | Sev | Title | Location | +|----|-----|-------|----------| +| V1 | 🟑 | JSON deserialization has no error handling β€” nlohmann exceptions leak the abstraction, not the library's `runtime_error` | `serialization.hpp:96-110` | +| V2 | 🟑 | XML parser uses unchecked `std::stoi`/`std::stod` on attributes β†’ uncaught `std::invalid_argument` on garbage | `serialization.hpp:258-261` | +| V3 | 🟑 | No validation of vertex coordinate values on load β€” NaN/Inf coordinates flow silently into the solver | `mesh_io.hpp:56-66` | +| V4 | 🟑 | JSON loader accesses nested keys (`j["solver"]["converged"]`) without `contains` checks β†’ throws `type_error`/creates nulls on schema drift | `serialization.hpp:104-108` | +| V5 | πŸ”΅ | Hand-rolled XML parser assumes one element per line and no escaping β€” silently mis-reads valid-but-reformatted XML | `serialization.hpp:248-270` | +| V6 | πŸ”΅ | No upper bound / sanity check on DOF-vector length vs mesh size when loading | `serialization.hpp` | + +--- + +## V1 β€” 🟑 JSON deserialization leaks nlohmann exceptions + +### Evidence +`serialization.hpp:96-108`: +```cpp +std::ifstream ifs(path); +if (!ifs) throw std::runtime_error("Cannot open: " + path); +json j; ifs >> j; // ← malformed JSON β†’ nlohmann::parse_error +... +res->converged = j["solver"]["converged"].get(); // ← wrong type β†’ type_error +``` + +### Problem +The "file missing" case is handled with the library's own `std::runtime_error` +(consistent with `mesh_io.hpp`). But the moment the file *exists and is malformed*, +`ifs >> j` throws `nlohmann::json::parse_error` and the typed accessors throw +`nlohmann::json::type_error` β€” neither is the library's documented error type, and +neither carries the file path. A caller writing `catch (const std::runtime_error&)` +(the pattern the examples use) will **miss** these. (`nlohmann` exceptions derive from +`std::exception`, not `std::runtime_error`.) + +### Fix +Wrap the parse + extraction in `try { ... } catch (const json::exception& e)` and +rethrow as `std::runtime_error("conformallab: malformed JSON in " + path + ": " + e.what())`, +matching the `mesh_io.hpp` boundary convention. + +### Acceptance criteria +- Loading a truncated/garbage JSON throws `std::runtime_error` with the path. +- A negative test (`EXPECT_THROW(..., std::runtime_error)`) covers it. + +--- + +## V2 β€” 🟑 XML parser: unchecked stoi/stod on attributes + +### Evidence +`serialization.hpp:258-261`: +```cpp +res->iterations = std::stoi(detail_xml::xml_get_attr(line, "iterations")); +res->grad_inf_norm = std::stod(detail_xml::xml_get_attr(line, "grad_inf_norm")); +``` + +### Problem +If the attribute is missing or non-numeric (`xml_get_attr` returns `""` or garbage), +`std::stoi("")` throws `std::invalid_argument` and out-of-range values throw +`std::out_of_range` β€” both uncaught, both not the library's `runtime_error`, neither +carrying the path. Same abstraction leak as V1, different exception family. + +### Fix +Parse defensively: check the attribute is non-empty and numeric, or wrap in +try/catch and rethrow as `std::runtime_error` with context (path + attribute name). + +### Acceptance criteria +- An XML file with a non-numeric `iterations` attribute throws `std::runtime_error`, + not `std::invalid_argument`. + +--- + +## V3 β€” 🟑 No validation of vertex coordinates on load + +### Evidence +`mesh_io.hpp:56-66` β€” `load_mesh` checks read success and `is_triangle_mesh`, but +does **not** inspect the vertex coordinate values. + +### Problem +A mesh file containing `nan`/`inf` coordinates (common from a crashed upstream tool, +or `1e308`-style overflow) reads in cleanly. The NaN then propagates: edge lengths +β†’ `lambda0` β†’ gradient β†’ Hessian β†’ the Newton solver, where it silently poisons the +result (`converged` may even read `true` if `NaN < tol` evaluates falsely in a way +that exits the loop). This is the worst kind of failure: silent, late, and +hard to trace back to the input. + +### Fix +After `read_mesh`, scan vertex coordinates for finiteness +(`std::isfinite`) and throw `std::runtime_error("conformallab: non-finite vertex +coordinate in " + filename)` on the first offender. O(V), negligible cost at the I/O +boundary β€” the same philosophy as the existing `is_triangle_mesh` guard. + +### Acceptance criteria +- Loading a mesh with a NaN/Inf coordinate throws `std::runtime_error`. +- Covered by a negative test. + +--- + +## V4 β€” 🟑 JSON loader assumes schema without checking + +### Evidence +`serialization.hpp:104-108` reads `j["solver"]["converged"]`, `["iterations"]`, +`["grad_inf_norm"]` and `j.at("dof_vector")` β€” only `geometry` is guarded with +`j.contains(...)` (`:99`). + +### Problem +`operator[]` on a missing key in nlohmann **inserts a null** (for non-const json) or +throws; `.get()` on null throws `type_error`. So a result file written by a +*future or older* schema version fails with an opaque nlohmann error rather than a +clear "missing field X" message. (`dof_vector` correctly uses `.at()`, which throws a +clear `out_of_range` β€” but still not the library's type.) + +### Fix +Validate the expected keys up front (or use `.at()` uniformly + the V1 try/catch), +producing a single `std::runtime_error("conformallab: result JSON missing field +'solver.converged'")`-style message. + +### Acceptance criteria +- A JSON missing `solver.iterations` yields a `runtime_error` naming the field. + +--- + +## V5 β€” πŸ”΅ Hand-rolled XML parser is format-fragile + +### Evidence +`serialization.hpp:248-270` parses line-by-line with `line.find("…` text) fits patterns on one +line. + +### Problem +This is not an XML parser β€” it is a line-pattern matcher. A semantically identical +file that pretty-prints elements across multiple lines, or adds an XML declaration / +namespace / attribute reordering, will be **silently mis-read** (fields default to +zero/empty) rather than rejected. Round-trips written by this same code work; nothing +else is guaranteed. + +### Fix +Either (a) document the XML format as a strict internal-only subset (and reject input +that doesn't match, rather than mis-reading it), or (b) if interop matters, use a real +XML parser. Given JSON is the primary format, (a) is the pragmatic choice β€” but the +"silently mis-read" behavior must become "explicitly reject". + +### Acceptance criteria +- Reformatted-but-valid XML is either parsed correctly or rejected with an error β€” + never silently mis-read into zeros. + +--- + +## V6 β€” πŸ”΅ No DOF-vector vs mesh sanity check + +### Evidence +`load_result_json` returns `x = j.at("dof_vector")` with no check that `x.size()` +matches the DOF count of the mesh it will be applied to. + +### Problem +A result file from a *different* mesh loads happily; the size mismatch only surfaces +later (out-of-bounds or wrong-answer) when `x` is indexed against the new mesh. + +### Fix +Where the loaded `x` is paired with a mesh, assert/throw on a size mismatch with a +clear message. (May belong at the call site rather than the loader.) + +--- + +## What is already good + +- The **file-missing** path is handled correctly and consistently + (`runtime_error` + path) in both `mesh_io.hpp` and `serialization.hpp`. +- `load_mesh` already enforces the triangle-mesh precondition at the boundary β€” the + right place and the right philosophy; V3 just extends it to coordinate finiteness. +- `dof_vector` uses `.at()` (throwing) rather than `operator[]` β€” partial good + practice that V4 asks to make uniform. +- JSON (nlohmann) is a robust, well-tested parser β€” the gap is only the missing + try/catch *around* it (V1), not the parser itself. + +## Suggested order +1. **V3** (NaN/Inf coordinate guard) β€” highest impact, prevents silent solver poisoning. +2. **V1 + V4** (JSON error handling + schema checks) β€” do together. +3. **V2** (XML stoi/stod) β€” same pattern as V1. +4. **V5 + V6** (XML strictness, size check) β€” polish. diff --git a/doc/reviewer/numerical-stability-audit-2026-05-31.md b/doc/reviewer/numerical-stability-audit-2026-05-31.md new file mode 100644 index 0000000..71d957c --- /dev/null +++ b/doc/reviewer/numerical-stability-audit-2026-05-31.md @@ -0,0 +1,234 @@ +# Numerical-Stability Audit β€” ConformalLabpp + +**Date:** 2026-05-31 +**Auditor:** External reviewer (Claude Opus 4.8) +**Scope:** Tolerance hierarchy, conditioning, catastrophic cancellation, and the +hard-coded "magic" constants across `code/include/`. +**Focus:** Floating-point robustness of a numerical library β€” distinct from +port-faithfulness (`java-port-audit.md`) and from internal error propagation +(`test-coverage-error-handling-audit-2026-05-31.md`). + +Status legend: πŸ”΄ Critical Β· 🟑 Important Β· πŸ”΅ Polish + +> **Relationship to existing audits:** `java-port-audit.md` already verified several +> per-formula numerical edge cases (degenerate triangles #1, overflow centring #10). +> This audit takes the *cross-cutting* view those miss: how the tolerances relate to +> each other, where constants are hard-coded vs. centralized, and where the formulas +> can lose precision. No finding here duplicates java-port-audit. + +--- + +## Summary table + +| ID | Sev | Title | Location | +|----|-----|-------|----------| +| N1 | πŸ”΄ | FD-Hessian step (`1e-5`) sets an accuracy floor (~`1e-10`) too close to the Newton tol (`1e-8`) β€” can stall quadratic convergence | `newton_solver.hpp:408,569` | +| N2 | 🟑 | Tolerance hierarchy is undocumented and uncoordinated β€” 7+ independent literals with no stated relationship | many headers | +| N3 | 🟑 | Discontinuous clamp `b<0 β†’ 0.01` breaks gradient/Hessian continuity at the domain boundary | `hyper_ideal_functional.hpp:179,274` | +| N4 | 🟑 | Magic constants (`-30.0`, `0.01`, `1e-15`, `1.0-1e-15`) are unnamed, undocumented, and scattered | `spherical_functional.hpp`, `spherical_geometry.hpp`, `hyper_ideal_functional.hpp` | +| N5 | 🟑 | Cotangent/area formula is cancellation-prone for needle/cap triangles | `euclidean_hessian.hpp:81-93` | +| N6 | πŸ”΅ | `constants.hpp` holds only Ο€/2Ο€ β€” all tolerances live as bare literals | `constants.hpp` | +| N7 | πŸ”΅ | No conditioning diagnostic on the Hessian before the linear solve | `newton_solver.hpp` | + +--- + +## N1 β€” πŸ”΄ FD-Hessian step size fights the convergence tolerance + +### Evidence +- `newton_hyper_ideal` and `newton_inversive_distance` build the Hessian by **finite + differences** with `hess_eps = 1e-5` (`newton_solver.hpp:408`, `:569`). +- Both solvers converge against `tol = 1e-8` on `β€–Gβ€–βˆž` (`:406`, `:567`). + +### Problem +A central-difference Hessian with step `h = 1e-5` has error `O(hΒ²) + O(Ξ΅_machine/h)`. +With `h = 1e-5` the truncation term is ~`1e-10` and the round-off term is +~`1e-16/1e-5 = 1e-11`; the achievable Hessian accuracy floor is therefore +**~`1e-10`–`1e-11`**. That is only ~2–3 decades below the convergence target +`tol = 1e-8`. Near the solution, where Newton should converge quadratically, an +inaccurate Hessian degrades the step to linear convergence and can **stall** the +iteration just above `tol` β€” exactly the regime `test_newton_phase9a.cpp` exercises. + +This is also why the line-search "stall" exit (FINDING-H1 / error-handling audit) +can trigger on otherwise well-posed problems: the FD Hessian, not the geometry, is +the limiting factor. + +### Fix +- Primary: replace the FD Hessian with the analytic / block-FD paths (cross-ref + `api-performance-audit` B1 β€” the HyperIdeal block-FD already exists and is more + accurate as well as faster). +- Until then: document the coupling, and consider a richer FD stencil or + Richardson extrapolation if FD must stay. A relative step `h = 1e-6Β·(1+|x_i|)` + is more robust than an absolute `1e-5`. + +### Acceptance criteria +- The FD-step/tol relationship is documented, or the analytic/block-FD Hessian + replaces FD so the floor drops well below `tol`. +- A test demonstrates quadratic (not stalled) convergence near the solution. + +--- + +## N2 β€” 🟑 Tolerance hierarchy is undocumented and uncoordinated + +### Evidence (the full set of independent literals) +| Constant | Value | Location | Meaning | +|---|---|---|---| +| Newton `tol` | `1e-8` | `newton_solver.hpp` (all 5 solvers) | convergence on `β€–Gβ€–βˆž` | +| FD `hess_eps` | `1e-5` | `newton_solver.hpp:408,569`, `hyper_ideal_hessian.hpp:67` | Hessian FD step | +| `gradient_check` `eps`/`tol` | `1e-5` / `1e-4` | `hyper_ideal_functional.hpp:474-475` | test-only FD check | +| Armijo `c1` | `1e-4` | `newton_solver.hpp:149` | sufficient-decrease | +| sparsity drop | `1e-15` | `newton_solver.hpp:599` | triplet prune threshold | +| Gauss-Bonnet `tol` | `1e-8` | `gauss_bonnet.hpp:134,154` | Ο‡ residual | +| spherical clamp | `1e-15` | `spherical_geometry.hpp:34` | `asin` domain guard | + +### Problem +These seven thresholds were each chosen locally (several "to match the Java +`FunctionalTest`"), with **no documented relationship**. Yet they interact: N1 shows +`hess_eps` vs `tol`; the sparsity drop `1e-15` interacts with `hess_eps` (FD values +below `1e-15` are pruned, but FD noise is ~`1e-11`, so the prune threshold is +effectively never the limiting factor β€” possibly dead). A reviewer cannot tell +which are load-bearing and which are arbitrary. + +### Fix +Add a short `doc/math/tolerances.md` (or a section in `constants.hpp`) that states +each tolerance, its role, and the constraints between them (e.g. "FD floor must be +β‰ͺ Newton tol", "Armijo c1 ∈ [1e-4, 1e-1]"). Centralize them as named `constexpr`. + +### Acceptance criteria +- Every numerical threshold is named, documented, and its relationship to the + others stated. + +--- + +## N3 β€” 🟑 Discontinuous clamp breaks gradient/Hessian continuity + +### Evidence +`hyper_ideal_functional.hpp:179-181` and `:274-276`: +```cpp +if (v1b && b1 < 0.0) b1 = 0.01; // negative scale β†’ snap to 0.01 +``` + +### Problem +This snaps any negative `b` to a constant `0.01`. The map `b ↦ max(b, 0.01)` is +continuous but **not differentiable** at `b = 0`, and the snap-to-constant makes the +*returned angles* locally independent of `b` for `b < 0` β€” so the analytic gradient +and the FD Hessian disagree across the boundary, and Newton steps that cross `b = 0` +get an inconsistent local model. The comment says it "mirrors the Java original", +but the Java code is a reference, not a correctness proof: a hard clamp inside the +function being differentiated is a known source of Newton stalls. + +### Fix +- Prefer keeping iterates in the feasible region via the line search (reject steps + that drive `b < b_min`) rather than clamping inside the energy evaluation. +- If clamping must stay, use a smooth barrier / softplus so the derivative is + continuous, and make `b_min` a named constant (N4). + +### Acceptance criteria +- The energy/gradient used by Newton is CΒΉ across the feasibility boundary, or the + line search keeps iterates strictly feasible so the clamp never fires. + +--- + +## N4 β€” 🟑 Unnamed magic constants + +### Evidence +- `spherical_functional.hpp`: `m.lambda0[e] = -30.0;` (very-short-edge floor) +- `spherical_geometry.hpp:34`: `if (half >= 1.0) half = 1.0 - 1e-15;` +- `hyper_ideal_functional.hpp`: `0.01` clamp (N3), `b < 0.0` thresholds +- `newton_solver.hpp:599`: `if (std::abs(val) > 1e-15)` + +### Problem +Each is a bare literal with at most an inline comment. `-30.0` (β‰ˆ `exp(-15)`, i.e. +edge length ~`3e-7`) and `1.0 - 1e-15` (the `asin`/`acos` domain guard) encode real +numerical reasoning that is invisible to the next maintainer and impossible to tune +consistently. + +### Fix +Promote to named `constexpr` in `constants.hpp` with a one-line rationale each, e.g. +`constexpr double LOG_LENGTH_FLOOR = -30.0; // exp(-15): edge below ~3e-7 treated as 0`. + +### Acceptance criteria +- No bare numeric tolerance/floor literals remain in the math headers; each is a + documented named constant. + +--- + +## N5 β€” 🟑 Cotangent/area formula is cancellation-prone + +### Evidence +`euclidean_hessian.hpp:81-93`: +```cpp +const double t12 = -l12 + l23 + l31; +const double t23 = +l12 - l23 + l31; +const double t31 = +l12 + l23 - l31; +... +const double denom2 = 2.0 * std::sqrt(t12 * t23 * t31 * l123); // = 8Β·Area +``` + +### Problem +For a needle triangle (one edge β‰ˆ sum of the other two), one of `t12/t23/t31` is the +difference of nearly-equal lengths β†’ **catastrophic cancellation**, and the area via +`sqrt` of the product loses precision. The cotangent weights then carry large relative +error, which feeds directly into the Hessian and the linear solve. The `<= 0` guard +catches the fully-degenerate case but not the *ill-conditioned-but-valid* one. + +### Fix +Use a numerically stable area formula (Kahan's formula for triangle area from sorted +side lengths) instead of the raw product-under-sqrt. This is a well-known, drop-in +improvement for cotangent-Laplacian assembly. + +### Acceptance criteria +- A sliver-triangle test shows the cotangent weights match a high-precision + (e.g. `long double` / Kahan) reference to acceptable relative error. + +--- + +## N6 β€” πŸ”΅ constants.hpp holds only Ο€ + +### Evidence +`constants.hpp` defines exactly `PI` and `TWO_PI`. Every tolerance from N2 lives +elsewhere as a literal. + +### Fix +Extend `constants.hpp` (or a sibling `tolerances.hpp`) to be the single source of +truth for the N2/N4 constants. Tie to N2's documentation. + +--- + +## N7 β€” πŸ”΅ No conditioning diagnostic before the linear solve + +### Evidence +`solve_with_fallback` (`newton_solver.hpp:65`) tries LDLT then SparseQR but never +reports the conditioning of `H`. + +### Problem +On an ill-conditioned-but-non-singular Hessian, LDLT "succeeds" and returns a +low-accuracy step with no signal. The `sparse_qr_fallback_used` flag only fires on +outright LDLT failure, not on near-singularity. + +### Fix +Optionally expose an estimate (e.g. smallest pivot magnitude from LDLT, or a cheap +condition estimate) in the diagnostics, so callers/tests can flag near-singular +solves. Low priority; pairs naturally with the FINDING-I1 status enum +(error-handling audit). + +--- + +## What is already good + +- `constants.hpp` centralizes Ο€ at 30 digits β€” the right instinct, just not applied + to tolerances (N6). +- The spherical `asin` domain guard (`spherical_geometry.hpp:34`) exists at all β€” + many implementations forget it (it just needs a named constant, N4). +- Central differences (not forward) are used for FD gradients/Hessians β€” the correct + O(hΒ²) choice. +- `java-port-audit.md` #10 already flagged the inversive-distance edge-length overflow + centring β€” consistent with this audit's spirit; coordinate the two. + +## Suggested order + +1. **N1** (FD step vs tol) β€” highest correctness impact; largely subsumed by the + `api-performance-audit` B1 analytic/block-FD work. +2. **N3** (discontinuous clamp) β€” fixes a real Newton-stall source. +3. **N5** (stable area) β€” drop-in robustness for the Euclidean Hessian. +4. **N2 + N4 + N6** (document + centralize tolerances/constants) β€” do together. +5. **N7** (conditioning diagnostic) β€” polish, pairs with I1. diff --git a/doc/reviewer/thread-safety-audit-2026-05-31.md b/doc/reviewer/thread-safety-audit-2026-05-31.md new file mode 100644 index 0000000..1d14a92 --- /dev/null +++ b/doc/reviewer/thread-safety-audit-2026-05-31.md @@ -0,0 +1,139 @@ +# Thread-Safety & Reentrancy Audit β€” ConformalLabpp + +**Date:** 2026-05-31 +**Auditor:** External reviewer (Claude Opus 4.8) +**Scope:** Concurrency contract of the header-only library β€” reentrancy, hidden +mutable global state, and safety of parallel solves. +**Focus:** Can two threads use this library at the same time, and is that contract +documented? Currently undocumented everywhere. + +Status legend: πŸ”΄ Critical Β· 🟑 Important Β· πŸ”΅ Polish + +> **Good news up front:** a scan for mutable function-local `static` variables (the +> classic hidden-shared-state hazard) found **none**. All `static` in the headers are +> stateless factory/accessor functions (`MobiusMap::identity`, +> `Default_conformal_map_traits::theta_map`, etc.). So the library is *close* to +> trivially reentrant β€” the findings below are about the mesh-mutation contract and +> the missing documentation, not about data races in hidden globals. + +--- + +## Summary table + +| ID | Sev | Title | Location | +|----|-----|-------|----------| +| T1 | 🟑 | No documented threading contract anywhere β€” users cannot know what is safe | all public headers | +| T2 | 🟑 | The solve pipeline mutates the mesh via property maps β†’ concurrent solves on the *same* mesh are unsafe and silently so | `*_functional.hpp`, `newton_solver.hpp` | +| T3 | πŸ”΅ | Property-map *creation* (`add_property_map`) on a shared mesh from two threads is a data race on CGAL's internal map registry | `setup_*_maps`, traits accessors | +| T4 | πŸ”΅ | Eigen threading not pinned β€” relies on Eigen's defaults; undocumented | `newton_solver.hpp` | + +--- + +## T1 β€” 🟑 No documented threading contract + +### Evidence +No header, no `doc/` file, and no README section states whether the API is +thread-safe, reentrant, or single-threaded-only. `briefing.md` positions the library +as potential "infrastructure for your work" β€” which invites multi-threaded use it +never characterizes. + +### Problem +Without a stated contract, a user must read the implementation to discover (T2) that +distinct meshes are safe but a shared mesh is not. The safe pattern is real and +simple β€” it just needs to be written down. + +### Fix +Add a short "Thread safety" section to the main API doc (and a sentence in the key +headers) stating the contract, e.g.: +> *The library performs no synchronization and uses no shared mutable global state. +> Operations on **distinct** `ConformalMesh` objects are independent and may run +> concurrently. Operations on a **single** mesh (including any `setup_*_maps`, +> `assign_*`, gradient, Hessian, or `newton_*` call) mutate that mesh's property maps +> and must not run concurrently on the same mesh.* + +### Acceptance criteria +- A "Thread safety" subsection exists in the public API documentation and matches the + actual behavior described in T2/T3. + +--- + +## T2 β€” 🟑 Solve pipeline mutates the mesh; same-mesh concurrency is unsafe + +### Evidence +- `setup_*_maps(mesh)` adds property maps to `mesh`. +- `assign_*_dof_indices(mesh, m)` writes `m.v_idx` / `m.e_idx` into the mesh. +- `compute_*_lambda0_from_mesh(mesh, m)` writes `m.lambda0`. +- `newton_*` reads those maps and (via the gradient/Hessian) reads the mesh + repeatedly; the functions take `ConformalMesh& mesh` (non-const). + +### Problem +Two threads calling `newton_euclidean(mesh, …)` on the **same** `mesh` race on the +shared property maps β€” a silent data race, not a clean error. Because the signature +is `ConformalMesh&` (mutable) this is technically "expected", but nothing warns the +user, and the natural-looking "solve the same input mesh under two geometries in +parallel" is exactly the unsafe case. + +### Fix +- Document the contract (T1). +- Optionally, offer a value-semantics convenience overload that takes the mesh by + `const&` and operates on an internal copy, for users who want fire-and-forget + parallel solves. (Larger change; document first.) + +### Acceptance criteria +- The same-mesh hazard is documented; if a const/copying overload is added, a test + runs two concurrent solves on copies without a race (clean under ThreadSanitizer). + +--- + +## T3 β€” πŸ”΅ Concurrent property-map creation races CGAL's registry + +### Evidence +`setup_*_maps` and the traits accessors call `mesh.add_property_map<…>(name, default)`. + +### Problem +CGAL's `Surface_mesh::add_property_map` mutates the mesh's internal property +container. Two threads creating maps on the same mesh concurrently is undefined +behavior at the CGAL level β€” a sub-case of T2 but worth calling out because the +traits accessors (`Conformal_map_traits.h:158-179`) create-on-demand, so even a +"read-only-looking" first access mutates. + +### Fix +Covered by the T1 contract ("any `setup_*` … on the same mesh must not run +concurrently"). No code change beyond documentation unless the const overload (T2) is +pursued. + +--- + +## T4 β€” πŸ”΅ Eigen threading is unpinned and undocumented + +### Evidence +`newton_solver.hpp` uses `SimplicialLDLT` / `SparseQR` without calling +`Eigen::setNbThreads(…)` or documenting OpenMP behavior. + +### Problem +Eigen may use OpenMP threads internally depending on build flags. For a library that +might itself be called from parallel user code, the nesting behavior (and whether +`EIGEN_DONT_PARALLELIZE` is wanted) is unspecified. Distinct Eigen objects are +thread-safe, so this is low risk β€” but it should be a documented choice. + +### Fix +State the intended Eigen threading mode in the build/threading docs; consider pinning +it explicitly (`EIGEN_DONT_PARALLELIZE` for predictability, or document that callers +control it). + +--- + +## What is already good + +- **No mutable static locals** anywhere in the headers β€” the single most important + reentrancy property holds. +- All `static` members are pure functions (factories/accessors) β€” no shared state. +- The header-only, value-oriented design (results returned by value as + `NewtonResult` / `Conformal_map_result`) means outputs are not shared. +- Distinct-mesh concurrency genuinely *is* safe today β€” the contract just needs + documenting (T1), not implementing. + +## Suggested order +1. **T1** (write the contract) β€” trivial, unblocks safe parallel use immediately. +2. **T2/T3** (covered by the contract; optional const overload later). +3. **T4** (pin/document Eigen threading) β€” polish.