diff --git a/.gitea/workflows/cpp-tests.yml b/.gitea/workflows/cpp-tests.yml index 1c3db65..84a4396 100644 --- a/.gitea/workflows/cpp-tests.yml +++ b/.gitea/workflows/cpp-tests.yml @@ -41,7 +41,11 @@ jobs: run: cmake -S code -B build -DCMAKE_BUILD_TYPE=Release - name: Build - run: nice -n 19 cmake --build build --target conformallab_tests -j$(nproc) + # Serial build (-j1): the 800m-limited container OOM-kills cc1plus when + # several Eigen+GoogleTest translation units compile in parallel + # (`-j$(nproc)`), failing test-fast non-deterministically regardless of + # branch content. Mirrors the Pi-protection already used by test-cgal. + run: nice -n 19 cmake --build build --target conformallab_tests -j1 - name: Run tests run: > @@ -160,9 +164,18 @@ jobs: - name: shellcheck (scripts/**/*.sh, severity=warning, strict) run: bash scripts/quality/shellcheck.sh --strict + - name: Install lcov (coverage gate) + run: apt-get install -y --no-install-recommends lcov + + - name: Coverage gate (fast-test suite, ramp-up mode) + # SKIP_COVERAGE_GATE=1: reports numbers but does not fail on threshold. + # Remove once I5 is resolved (coverage is wired to the full CGAL suite, + # not just the 26 fast tests). Threshold: 80% line / 70% branch / 90% func. + run: SKIP_COVERAGE_GATE=1 bash scripts/quality/coverage.sh + - name: Summary if: always() run: | - echo "QUALITY ▸ all four gates passed." + echo "QUALITY ▸ all gates passed." echo " see scripts/quality/README.md for the full catalogue" - echo " (sanitizers, clang-tidy, coverage, etc. are local-only)" + echo " (sanitizers, clang-tidy, coverage report in build-coverage/)" diff --git a/code/include/CGAL/Discrete_circle_packing.h b/code/include/CGAL/Discrete_circle_packing.h index 3942e76..3882eee 100644 --- a/code/include/CGAL/Discrete_circle_packing.h +++ b/code/include/CGAL/Discrete_circle_packing.h @@ -8,7 +8,7 @@ \ingroup PkgConformalMapRef User-facing entry for the **face-based** circle-packing functional of -Bobenko-Pinkall-Springborn 2010. See `cp_euclidean_functional.hpp` +Bobenko-Pinkall-Springborn 2015. See `cp_euclidean_functional.hpp` for the underlying algorithm and `doc/architecture/phase-9a-validation.md` for the line-by-line mapping to the Java original `CPEuclideanFunctional.java`. @@ -44,7 +44,7 @@ namespace CGAL { \ingroup PkgConformalMapConcepts \brief Traits class for `discrete_circle_packing_euclidean()` — declares the kernel, mesh and property-map types used by the -BPS-2010 face-based circle-packing functional. +BPS-2015 face-based circle-packing functional. Primary template; specialise it for non-`Surface_mesh` triangle meshes. */ @@ -96,6 +96,9 @@ struct Default_cp_euclidean_traits, K> Result of `discrete_circle_packing_euclidean`. Carries face DOFs `ρ_f = log R_f` rather than the vertex DOFs of the classical modes. */ +/// Why the Newton solve terminated (re-exported from the solver; I1 audit). +using Newton_status = ::conformallab::NewtonStatus; + template struct Circle_packing_result { @@ -108,6 +111,14 @@ struct Circle_packing_result FT gradient_norm = FT(0); /// `true` iff `gradient_norm < gradient_tolerance` at exit. bool converged = false; + + /// Why the solve stopped: `Converged` / `MaxIterations` / + /// `LinearSolverFailed` / `LineSearchStalled`. + Newton_status status = Newton_status::MaxIterations; + /// `true` iff any Newton step fell back to SparseQR (gauge singularity). + bool sparse_qr_fallback_used = false; + /// Smallest `|Dᵢᵢ|` of the last LDLT (near-singularity proxy; 0 if none). + FT min_ldlt_pivot = FT(0); }; // ── Entry function ──────────────────────────────────────────────────────────── @@ -115,7 +126,7 @@ struct Circle_packing_result /*! \ingroup PkgConformalMapRef -Compute the BPS-2010 face-based circle-packing of `mesh`. +Compute the BPS-2015 face-based circle-packing of `mesh`. \tparam TriangleMesh A `CGAL::Surface_mesh

`. \tparam NamedParameters Optional CGAL named-parameter pack. @@ -132,7 +143,7 @@ Compute the BPS-2010 face-based circle-packing of `mesh`. \returns A `Circle_packing_result` with `ρ_f` per face. \pre `mesh` is a triangle mesh. -\pre `φ_f` and `θ_e` satisfy the BPS-2010 admissibility conditions +\pre `φ_f` and `θ_e` satisfy the BPS-2015 admissibility conditions (Σ_f φ_f = 2π·χ + Σ_e (π − θ_e), see paper §6). */ template (nr.min_ldlt_pivot); // ── output_uv_map (Phase 8b-Lite extension) ──────────────────────────── // @@ -201,7 +215,7 @@ auto discrete_circle_packing_euclidean( // // For Phase 8b-Lite we deliberately don't fake it. If the caller // supplies `output_uv_map(pmap)` we throw `std::runtime_error` with a - // clear pointer to Phase 9c (BPS-2010 §6 face-based circle-packing + // clear pointer to Phase 9c (BPS-2015 §6 face-based circle-packing // layout, ~150 lines, on the porting roadmap). Failing loudly is // better than silently writing zeros. // diff --git a/code/include/CGAL/Discrete_conformal_map.h b/code/include/CGAL/Discrete_conformal_map.h index 4021796..9e81359 100644 --- a/code/include/CGAL/Discrete_conformal_map.h +++ b/code/include/CGAL/Discrete_conformal_map.h @@ -78,6 +78,10 @@ namespace CGAL { // Result type // ════════════════════════════════════════════════════════════════════════════ +/// Why the Newton solve terminated (re-exported from the solver; I1 audit). +/// `Converged` / `MaxIterations` / `LinearSolverFailed` / `LineSearchStalled`. +using Newton_status = ::conformallab::NewtonStatus; + /*! \ingroup PkgConformalMapRef @@ -100,9 +104,18 @@ struct Conformal_map_result /// `true` iff `gradient_norm < gradient_tolerance`. bool converged = false; + /// Why the solve stopped (more informative than `converged` alone): one of + /// `Converged` / `MaxIterations` / `LinearSolverFailed` / `LineSearchStalled`. + Newton_status status = Newton_status::MaxIterations; + /// `true` iff the linear solver used the SparseQR fallback at any /// Newton step (gauge mode on closed mesh without pinned vertex). bool sparse_qr_fallback_used = false; + + /// Smallest `|Dᵢᵢ|` of the last LDLT factorisation — a cheap + /// near-singularity proxy (small ⇒ ill-conditioned solve). `0` if LDLT + /// never succeeded. + FT min_ldlt_pivot = FT(0); }; @@ -272,6 +285,9 @@ auto discrete_conformal_map_euclidean( result.iterations = nr.iterations; result.gradient_norm = nr.grad_inf_norm; result.converged = nr.converged; + result.status = nr.status; + result.sparse_qr_fallback_used = nr.sparse_qr_fallback_used; + result.min_ldlt_pivot = static_cast(nr.min_ldlt_pivot); // ── 8. Optional layout step (Phase 8b-Lite extension) ────────────────── // @@ -410,6 +426,9 @@ auto discrete_conformal_map_spherical( result.iterations = nr.iterations; result.gradient_norm = nr.grad_inf_norm; result.converged = nr.converged; + result.status = nr.status; + result.sparse_qr_fallback_used = nr.sparse_qr_fallback_used; + result.min_ldlt_pivot = static_cast(nr.min_ldlt_pivot); // Optional 3-D layout step (point on S²) auto uv_param = parameters::get_parameter( @@ -459,6 +478,14 @@ struct Hyper_ideal_map_result FT gradient_norm = FT(0); /// `true` iff `gradient_norm < gradient_tolerance` at exit. bool converged = false; + + /// Why the solve stopped: `Converged` / `MaxIterations` / + /// `LinearSolverFailed` / `LineSearchStalled`. + Newton_status status = Newton_status::MaxIterations; + /// `true` iff any Newton step fell back to SparseQR (gauge singularity). + bool sparse_qr_fallback_used = false; + /// Smallest `|Dᵢᵢ|` of the last LDLT (near-singularity proxy; 0 if none). + FT min_ldlt_pivot = FT(0); }; /*! @@ -538,6 +565,9 @@ auto discrete_conformal_map_hyper_ideal( result.iterations = nr.iterations; result.gradient_norm = nr.grad_inf_norm; result.converged = nr.converged; + result.status = nr.status; + result.sparse_qr_fallback_used = nr.sparse_qr_fallback_used; + result.min_ldlt_pivot = static_cast(nr.min_ldlt_pivot); // Optional Poincaré-disk layout (2-D in the unit disk). auto uv_param = parameters::get_parameter( diff --git a/code/include/CGAL/Discrete_inversive_distance.h b/code/include/CGAL/Discrete_inversive_distance.h index 50f43b5..dcaebf1 100644 --- a/code/include/CGAL/Discrete_inversive_distance.h +++ b/code/include/CGAL/Discrete_inversive_distance.h @@ -8,7 +8,7 @@ \ingroup PkgConformalMapRef User-facing entry for the **vertex-based** inversive-distance circle- -packing functional of Luo (2004), with the Bowers-Stephenson (2004) +packing functional of Luo 2004, with the Bowers-Stephenson (2004) initialisation. See `inversive_distance_functional.hpp` for the underlying algorithm and `doc/roadmap/research-track.md` (item 9a.2) for the research-track classification — this functional has **no Java @@ -98,7 +98,7 @@ struct Default_inversive_distance_traits /*! \ingroup PkgConformalMapRef -Compute the Luo-2004 vertex-based inversive-distance circle packing of `mesh`. +Compute the Luo 2004 vertex-based inversive-distance circle packing of `mesh`. The per-edge constant `I_ij` is computed once at the start from the input 3-D geometry via the Bowers-Stephenson identity @@ -207,6 +207,9 @@ auto discrete_inversive_distance_map( result.iterations = nr.iterations; result.gradient_norm = nr.grad_inf_norm; result.converged = nr.converged; + result.status = nr.status; + result.sparse_qr_fallback_used = nr.sparse_qr_fallback_used; + result.min_ldlt_pivot = static_cast(nr.min_ldlt_pivot); // ── Optional layout step (Phase 8b-Lite extension) ───────────────────── // diff --git a/code/include/constants.hpp b/code/include/constants.hpp index a4d5afa..1068f31 100644 --- a/code/include/constants.hpp +++ b/code/include/constants.hpp @@ -17,4 +17,28 @@ constexpr double PI = 3.14159265358979323846264338328; /// 2π (full turn). constexpr double TWO_PI = 2.0 * PI; +// ── Numerical thresholds and bounds (N4/N6 audit) ────────────────────────── + +/// Edge length floor for log-scale functionals (euclidean_functional, spherical_functional). +/// exp(-30.0) ≈ 3e-7: edges below this length are treated as degenerate/zero. +/// Rationale: avoids log(tiny) overflow and enforces a minimum edge scale. +constexpr double LOG_EDGE_LENGTH_FLOOR = -30.0; + +/// Hyper-ideal scale factor lower bound (hyper_ideal_functional:179-181). +/// If b (scale) becomes negative (infeasible), clamp to 0.01 to keep geometry valid. +/// Rationale: ensures b > 0 maintains the Lorentzian model's causal structure. +constexpr double HYPER_IDEAL_SCALE_FLOOR = 0.01; + +/// Sharpness β of the optional C¹ smooth-barrier scale floor (N3 audit). +/// Only used when the hyper-ideal clamp mode is `SmoothBarrier`; larger β +/// makes the softplus transition tighter around `HYPER_IDEAL_SCALE_FLOOR` +/// (β·floor ≈ 1, so β ≈ 100 keeps the barrier ≈ identity for b ≳ 0.05 while +/// staying C¹ across b = 0). See `clamp_hyper_ideal_scale`. +constexpr double HYPER_IDEAL_SCALE_SHARPNESS = 100.0; + +/// Domain guard for asin in spherical_l (spherical_geometry:34). +/// Clamps exp(λ/2) to slightly below 1.0 so asin stays in [−π/2, π/2]. +/// Rationale: asin(x) is undefined for |x| > 1; this prevents IEEE Inf/NaN. +constexpr double ASIN_DOMAIN_GUARD = 1.0 - 1e-15; + } // namespace conformallab diff --git a/code/include/cp_euclidean_functional.hpp b/code/include/cp_euclidean_functional.hpp index 8267591..f39998f 100644 --- a/code/include/cp_euclidean_functional.hpp +++ b/code/include/cp_euclidean_functional.hpp @@ -19,7 +19,7 @@ // │ The variable is ρ_f = log R_f. │ // │ Adjacent face-circles intersect at a prescribed angle θ_e per edge. │ // │ │ -// │ This is the FACE-DUAL of the classical vertex-based Luo (2004) │ +// │ This is the FACE-DUAL of the classical vertex-based Luo 2004 │ // │ inversive-distance circle packing implemented in │ // │ inversive_distance_functional.hpp (Phase 9a.2). The relation │ // │ I_ij = cos θ_e │ @@ -33,7 +33,7 @@ // │ θ_e per edge intersection angle of the two face-circles │ // │ φ_f per face target sum of corner-angles inside the face │ // │ │ -// │ Energy (BPS-2010 §6) │ +// │ Energy (BPS-2015 §6) │ // │ E(ρ) = Σ_f φ_f · ρ_f │ // │ + Σ_{(h,f=face(h)): │ // │ [ if opposite face exists ] │ @@ -52,7 +52,7 @@ // │ Per interior hf: −(p + θ*) added to G[face(h)] │ // │ Per boundary hf: −2 θ* added to G[face(h)] │ // │ │ -// │ Hessian (analytic, BPS-2010 eq. 6.8; Java getHessian lines 127-166) │ +// │ Hessian (analytic, BPS-2015 eq. 6.8; Java getHessian lines 127-166) │ // │ Per interior undirected edge e (connecting faces j and k): │ // │ h_jk = sin θ / (cosh Δρ − cos θ) │ // │ H[j,j] += h_jk, H[k,k] += h_jk, H[j,k] −= h_jk, H[k,j] −= h_jk │ @@ -90,7 +90,7 @@ using CPEMapD = ConformalMesh::Property_map; // ── Persistent map bundle ───────────────────────────────────────────────────── /// Bundle of the three property maps consumed by the CP-Euclidean -/// (Bobenko-Pinkall-Springborn 2010) circle-packing functional. +/// (Bobenko-Pinkall-Springborn 2015) circle-packing functional. struct CPEuclideanMaps { CPFMapI f_idx; ///< DOF index per face (−1 = pinned) CPEMapD theta_e; ///< intersection angle per edge (default π/2 = orthogonal) diff --git a/code/include/euclidean_functional.hpp b/code/include/euclidean_functional.hpp index f8cac0c..5a36966 100644 --- a/code/include/euclidean_functional.hpp +++ b/code/include/euclidean_functional.hpp @@ -162,7 +162,7 @@ inline void compute_euclidean_lambda0_from_mesh(ConformalMesh& mesh, EuclideanMa if (len > 1e-15) m.lambda0[e] = 2.0 * std::log(len); else - m.lambda0[e] = -30.0; // degenerate edge + m.lambda0[e] = LOG_EDGE_LENGTH_FLOOR; // degenerate edge } } diff --git a/code/include/euclidean_hessian.hpp b/code/include/euclidean_hessian.hpp index 83f473d..b4706ac 100644 --- a/code/include/euclidean_hessian.hpp +++ b/code/include/euclidean_hessian.hpp @@ -17,7 +17,8 @@ // │ side lengths lij = exp(Λ̃ij/2): │ // │ │ // │ t12 = −l12+l23+l31, t23 = l12−l23+l31, t31 = l12+l23−l31 │ -// │ l123 = l12+l23+l31, denom2 = 2·sqrt(t12·t23·t31·l123) = 8·Area │ +// │ l123 = l12+l23+l31, denom2 = 8·Area (Area via Kahan's stable │ +// │ side-length formula — see euclidean_cot_weights)│ // │ │ // │ Cotangent at vertex k (opposite t_opp, adjacent t_a and t_b): │ // │ cot_k = (t_opp · l123 − t_a · t_b) / denom2 │ @@ -48,6 +49,7 @@ #include #include #include +#include #include namespace conformallab { @@ -85,12 +87,26 @@ inline EuclCotWeights euclidean_cot_weights(double l12, double l23, double l31) if (t12 <= 0.0 || t23 <= 0.0 || t31 <= 0.0) return {0.0, 0.0, 0.0, false}; - const double l123 = l12 + l23 + l31; - const double denom2_sq = t12 * t23 * t31 * l123; - if (denom2_sq <= 0.0) return {0.0, 0.0, 0.0, false}; + const double l123 = l12 + l23 + l31; - // denom2 = 2·sqrt(t12·t23·t31·l123) = 8·Area - const double denom2 = 2.0 * std::sqrt(denom2_sq); + // Triangle area via Kahan's numerically stable formula. Sort the side + // lengths a ≥ b ≥ c and evaluate with the cancellation-avoiding grouping + // Area = ¼·√[ (a+(b+c))·(c−(a−b))·(c+(a−b))·(a+(b−c)) ]. + // This replaces the naive 2·√(t12·t23·t31·l123): that product of + // near-equal-length differences loses precision for needle/cap triangles, + // where it would feed large relative error straight into the cotangent + // weights and the linear solve. The result denom2 = 8·Area is identical + // up to rounding for well-shaped triangles, but accurate for slivers. + double a = l12, b = l23, c = l31; + if (a < b) std::swap(a, b); + if (a < c) std::swap(a, c); + if (b < c) std::swap(b, c); // now a ≥ b ≥ c + + const double kahan = (a + (b + c)) * (c - (a - b)) + * (c + (a - b)) * (a + (b - c)); + if (kahan <= 0.0) return {0.0, 0.0, 0.0, false}; + + const double denom2 = 8.0 * (0.25 * std::sqrt(kahan)); // = 8·Area // Formula: cot_k = (t_opp · l123 − t_a · t_b) / denom2 // cot1: t_opp=t23, t_a=t12, t_b=t31 (v1 opposite l23) diff --git a/code/include/hyper_ideal_functional.hpp b/code/include/hyper_ideal_functional.hpp index d5e028c..8533135 100644 --- a/code/include/hyper_ideal_functional.hpp +++ b/code/include/hyper_ideal_functional.hpp @@ -31,6 +31,7 @@ // // res.energy, res.gradient #include "conformal_mesh.hpp" +#include "constants.hpp" #include "hyper_ideal_geometry.hpp" #include "hyper_ideal_utility.hpp" #include @@ -157,6 +158,55 @@ static inline std::size_t hidx(Halfedge_index h) { return halfedge_to_index(h); // HyperIdealFunctional.java's defensive behaviour (lines 122-127 of the // Java original); this keeps the FD perturbation regime well-defined. +// ── Scale-floor clamp mode (N3 audit) ───────────────────────────────────────── +// +// The vertex scale b must stay positive for the hyper-ideal geometry to be +// valid. Two ways to enforce that are offered, selectable per call: +// +// • HardJava (DEFAULT) — the original `b < 0 → HYPER_IDEAL_SCALE_FLOOR` snap. +// Bit-for-bit faithful to HyperIdealFunctional.java, so the Java +// golden-oracle parity tests hold exactly. It is only C⁰ at b = 0: a +// Newton step that crosses the feasibility boundary sees a kink, which can +// stall convergence (numerical-stability audit N3). This is the +// production default precisely to preserve parity. +// +// • SmoothBarrier — a C¹ softplus floor +// b ↦ floor + softplus_β(b − floor), softplus_β(x) = log(1+e^{βx})/β +// which is ≈ b for b well above the floor (to machine precision once +// β·(b−floor) ≳ 35) and decays smoothly to `floor` as b → −∞, with a +// continuous derivative everywhere. This removes the N3 kink and is the +// mathematically clean choice, but it perturbs values near the boundary +// and therefore deviates from the Java oracle — opt in when robustness of +// a boundary-crossing solve matters more than strict parity. +// +// Both modes share the same floor (`HYPER_IDEAL_SCALE_FLOOR`); SmoothBarrier +// additionally uses `HYPER_IDEAL_SCALE_SHARPNESS` (β). The `a < 0 → 0` edge +// clamp is unaffected (a = 0 is a genuine geometric floor, not flagged by N3). +enum class HyperIdealScaleClamp { + HardJava, ///< b<0 → floor. C⁰, Java-parity-faithful (default). + SmoothBarrier ///< C¹ softplus floor; clean but deviates from Java near b=0. +}; + +/// Apply the vertex-scale floor to `b` under the chosen clamp `mode`. +/// `HardJava` reproduces the original snap; `SmoothBarrier` is the C¹ +/// softplus floor (see `HyperIdealScaleClamp`). `variable` mirrors the +/// call-site guard (only clamp DOFs that are actually free / variable). +inline double clamp_hyper_ideal_scale(double b, bool variable, + HyperIdealScaleClamp mode) noexcept +{ + if (!variable) return b; + if (mode == HyperIdealScaleClamp::HardJava) + return b < 0.0 ? HYPER_IDEAL_SCALE_FLOOR : b; + + // SmoothBarrier: floor + softplus_β(b − floor), evaluated stably. + const double beta = HYPER_IDEAL_SCALE_SHARPNESS; + const double x = beta * (b - HYPER_IDEAL_SCALE_FLOOR); + // softplus_β(x) = log1p(e^{βx})/β, with the standard large-x guard + // (for x ≳ 35, log1p(e^x) == x to double precision) to avoid overflow. + const double softplus = (x > 35.0) ? x : std::log1p(std::exp(x)); + return HYPER_IDEAL_SCALE_FLOOR + softplus / beta; +} + /// Six per-face angle outputs computed from local DOFs (see /// `face_angles_from_local_dofs`). Used by the block-FD Hessian. struct FaceAngleOutputs { @@ -175,15 +225,16 @@ struct FaceAngleOutputs { inline FaceAngleOutputs face_angles_from_local_dofs( double b1, double b2, double b3, double a12, double a23, double a31, - bool v1b, bool v2b, bool v3b) + bool v1b, bool v2b, bool v3b, + HyperIdealScaleClamp clamp = HyperIdealScaleClamp::HardJava) { // Same defensive clamps as compute_face_angles. if (v1b && v2b && a12 < 0.0) a12 = 0.0; if (v2b && v3b && a23 < 0.0) a23 = 0.0; if (v3b && v1b && a31 < 0.0) a31 = 0.0; - if (v1b && b1 < 0.0) b1 = 0.01; - if (v2b && b2 < 0.0) b2 = 0.01; - if (v3b && b3 < 0.0) b3 = 0.01; + b1 = clamp_hyper_ideal_scale(b1, v1b, clamp); + b2 = clamp_hyper_ideal_scale(b2, v2b, clamp); + b3 = clamp_hyper_ideal_scale(b3, v3b, clamp); double l12 = lij(b1, b2, a12, v1b, v2b); double l23 = lij(b2, b3, a23, v2b, v3b); @@ -246,7 +297,8 @@ static FaceAngles compute_face_angles( const ConformalMesh& mesh, Face_index f, const std::vector& x, - const HyperIdealMaps& m) + const HyperIdealMaps& m, + HyperIdealScaleClamp clamp = HyperIdealScaleClamp::HardJava) { Halfedge_index h0 = mesh.halfedge(f); Halfedge_index h1 = mesh.next(h0); @@ -276,9 +328,9 @@ static FaceAngles compute_face_angles( if (fa.v1b && fa.v2b && fa.a12 < 0.0) fa.a12 = 0.0; if (fa.v2b && fa.v3b && fa.a23 < 0.0) fa.a23 = 0.0; if (fa.v3b && fa.v1b && fa.a31 < 0.0) fa.a31 = 0.0; - if (fa.v1b && fa.b1 < 0.0) fa.b1 = 0.01; - if (fa.v2b && fa.b2 < 0.0) fa.b2 = 0.01; - if (fa.v3b && fa.b3 < 0.0) fa.b3 = 0.01; + fa.b1 = clamp_hyper_ideal_scale(fa.b1, fa.v1b, clamp); + fa.b2 = clamp_hyper_ideal_scale(fa.b2, fa.v2b, clamp); + fa.b3 = clamp_hyper_ideal_scale(fa.b3, fa.v3b, clamp); double l12 = lij(fa.b1, fa.b2, fa.a12, fa.v1b, fa.v2b); double l23 = lij(fa.b2, fa.b3, fa.a23, fa.v2b, fa.v3b); @@ -390,7 +442,8 @@ inline HyperIdealResult evaluate_hyper_ideal( const std::vector& x, const HyperIdealMaps& m, bool need_energy = true, - bool need_gradient = true) + bool need_gradient = true, + HyperIdealScaleClamp clamp = HyperIdealScaleClamp::HardJava) { HyperIdealResult res; @@ -406,7 +459,7 @@ inline HyperIdealResult evaluate_hyper_ideal( Halfedge_index h1 = mesh.next(h0); Halfedge_index h2 = mesh.next(h1); - FaceAngles fa = compute_face_angles(mesh, f, x, m); + FaceAngles fa = compute_face_angles(mesh, f, x, m, clamp); // Store computed angles into temporary arrays. // h_alpha[h] = α for the edge of h in this face. diff --git a/code/include/hyper_ideal_hessian.hpp b/code/include/hyper_ideal_hessian.hpp index 110062f..cacad36 100644 --- a/code/include/hyper_ideal_hessian.hpp +++ b/code/include/hyper_ideal_hessian.hpp @@ -64,7 +64,8 @@ inline Eigen::SparseMatrix hyper_ideal_hessian( ConformalMesh& mesh, const std::vector& x, const HyperIdealMaps& m, - double eps = 1e-5) + double eps = 1e-5, + HyperIdealScaleClamp clamp = HyperIdealScaleClamp::HardJava) { const int n = hyper_ideal_dimension(mesh, m); std::vector> trips; @@ -77,8 +78,8 @@ inline Eigen::SparseMatrix hyper_ideal_hessian( xp[sj] = x[sj] + eps; xm[sj] = x[sj] - eps; - auto Gp = evaluate_hyper_ideal(mesh, xp, m, /*energy=*/false).gradient; - auto Gm = evaluate_hyper_ideal(mesh, xm, m, /*energy=*/false).gradient; + auto Gp = evaluate_hyper_ideal(mesh, xp, m, /*energy=*/false, /*grad=*/true, clamp).gradient; + auto Gm = evaluate_hyper_ideal(mesh, xm, m, /*energy=*/false, /*grad=*/true, clamp).gradient; xp[sj] = xm[sj] = x[sj]; // restore @@ -101,9 +102,10 @@ inline Eigen::SparseMatrix hyper_ideal_hessian_sym( ConformalMesh& mesh, const std::vector& x, const HyperIdealMaps& m, - double eps = 1e-5) + double eps = 1e-5, + HyperIdealScaleClamp clamp = HyperIdealScaleClamp::HardJava) { - auto H = hyper_ideal_hessian(mesh, x, m, eps); + auto H = hyper_ideal_hessian(mesh, x, m, eps, clamp); Eigen::SparseMatrix Ht = H.transpose(); return (H + Ht) * 0.5; } @@ -142,7 +144,8 @@ inline Eigen::SparseMatrix hyper_ideal_hessian_block_fd( ConformalMesh& mesh, const std::vector& x, const HyperIdealMaps& m, - double eps = 1e-5) + double eps = 1e-5, + HyperIdealScaleClamp clamp = HyperIdealScaleClamp::HardJava) { const int n = hyper_ideal_dimension(mesh, m); std::vector> trips; @@ -187,9 +190,9 @@ inline Eigen::SparseMatrix hyper_ideal_hessian_block_fd( vm[j] -= eps; auto Op = face_angles_from_local_dofs( - vp[0], vp[1], vp[2], vp[3], vp[4], vp[5], v1b, v2b, v3b); + vp[0], vp[1], vp[2], vp[3], vp[4], vp[5], v1b, v2b, v3b, clamp); auto Om = face_angles_from_local_dofs( - vm[0], vm[1], vm[2], vm[3], vm[4], vm[5], v1b, v2b, v3b); + vm[0], vm[1], vm[2], vm[3], vm[4], vm[5], v1b, v2b, v3b, clamp); const double Gp[6] = { Op.beta1, Op.beta2, Op.beta3, @@ -221,9 +224,10 @@ inline Eigen::SparseMatrix hyper_ideal_hessian_block_fd_sym( ConformalMesh& mesh, const std::vector& x, const HyperIdealMaps& m, - double eps = 1e-5) + double eps = 1e-5, + HyperIdealScaleClamp clamp = HyperIdealScaleClamp::HardJava) { - auto H = hyper_ideal_hessian_block_fd(mesh, x, m, eps); + auto H = hyper_ideal_hessian_block_fd(mesh, x, m, eps, clamp); Eigen::SparseMatrix Ht = H.transpose(); return (H + Ht) * 0.5; } diff --git a/code/include/inversive_distance_functional.hpp b/code/include/inversive_distance_functional.hpp index 8d67a32..c973cf5 100644 --- a/code/include/inversive_distance_functional.hpp +++ b/code/include/inversive_distance_functional.hpp @@ -315,6 +315,43 @@ inline std::vector inversive_distance_gradient( return G; } +/// Per-face contribution to the Inversive-Distance gradient, as a pure +/// 3→3 kernel of the three local vertex DOFs `(u1,u2,u3)` and the three +/// edge inversive distances `(I12,I23,I31)`. No mesh, no property maps — +/// used by the per-face block-FD Hessian in `inversive_distance_hessian.hpp`. +/// +/// Returns the three values this face *adds* to the global gradient at +/// `(v1,v2,v3)`. Since `G_v = Θ_v − Σ_faces α_v`, the per-face contribution +/// is the **negative** corner angles `(−α₁,−α₂,−α₃)`. A non-real circle +/// configuration (any `ℓ² ≤ 0`) contributes nothing — exactly mirroring the +/// face-skip (`continue`) in `inversive_distance_gradient`, so the block-FD +/// Hessian and the full-FD Hessian see the same per-face support. +struct IDFaceGradContribs { + double g1; ///< contribution to G at v₁ (= −α₁) + double g2; ///< contribution to G at v₂ (= −α₂) + double g3; ///< contribution to G at v₃ (= −α₃) +}; + +inline IDFaceGradContribs inversive_distance_face_grad_contribs( + double u1, double u2, double u3, + double I12, double I23, double I31) +{ + double l12sq = id_detail::edge_length_squared(u1, u2, I12); + double l23sq = id_detail::edge_length_squared(u2, u3, I23); + double l31sq = id_detail::edge_length_squared(u3, u1, I31); + + // Same face-skip as inversive_distance_gradient: a non-real circle + // configuration has no limiting angle, so the face contributes 0. + if (l12sq <= 0.0 || l23sq <= 0.0 || l31sq <= 0.0) + return {0.0, 0.0, 0.0}; + + // euclidean_angles deliberately keeps its limiting (valid=false) angles + // here — the convex C¹ extension — matching the gradient's choice not to + // skip on !fa.valid. Corner at v_k = fa.alpha_k (see gradient Pass-2 trace). + auto fa = euclidean_angles(std::log(l12sq), std::log(l23sq), std::log(l31sq)); + return {-fa.alpha1, -fa.alpha2, -fa.alpha3}; +} + /// Inversive-Distance energy `E(u) = ∫₀¹ ⟨G(t·u), u⟩ dt`, evaluated /// with 10-point Gauss-Legendre (constants shared with `euclidean_energy`). inline double inversive_distance_energy( diff --git a/code/include/inversive_distance_hessian.hpp b/code/include/inversive_distance_hessian.hpp new file mode 100644 index 0000000..80a6d5e --- /dev/null +++ b/code/include/inversive_distance_hessian.hpp @@ -0,0 +1,187 @@ +#pragma once +// Copyright (c) 2024-2026 Tarik Moussa. +// SPDX-License-Identifier: MIT + +// inversive_distance_hessian.hpp +// +// Phase 9a.2 — Hessian of the inversive-distance circle-packing functional +// (Luo 2004 / Bowers-Stephenson 2004). +// +// ┌──────────────────────────────────────────────────────────────────────────┐ +// │ Implementation strategy │ +// │ │ +// │ TWO finite-difference Hessian implementations are provided here, │ +// │ mirroring the HyperIdeal pair in `hyper_ideal_hessian.hpp`: │ +// │ │ +// │ 1. `inversive_distance_hessian` — full finite-difference baseline. │ +// │ Cost ≈ n × (cost of a full gradient evaluation) = O(n · F). │ +// │ Kept as the cross-validation reference for the block-FD variant. │ +// │ │ +// │ 2. `inversive_distance_hessian_block_fd` — per-face block-FD. │ +// │ Each face contributes to the gradient through exactly 3 vertex │ +// │ DOFs (u₁,u₂,u₃); we FD the 3×3 local Jacobian of that face's │ +// │ gradient contribution and scatter it. Cost ≈ F × 6 face-angle │ +// │ evaluations = O(F). Speed-up factor ≈ n/6 over full-FD. │ +// │ │ +// │ Why the block-FD is correct (locality lemma): │ +// │ G_v = Θ_v − Σ_{f ∋ v} α_v(f), and α_v(f) depends ONLY on the 3 │ +// │ vertex DOFs of face f. Hence ∂G_x/∂y = Σ_{f: x,y ∈ {v1,v2,v3}(f)} │ +// │ ∂(−α_x)/∂y at f, so accumulating per-face 3×3 blocks reproduces the │ +// │ full Hessian (identical to O(ε²)). │ +// │ │ +// │ An analytic Hessian via Glickenstein 2011 eq. (4.6) is tracked in │ +// │ `doc/roadmap/research-track.md` as Phase 9a.2-analytic; it would take │ +// │ the cost from O(F)·(FD constant) to a single O(F) analytic pass. │ +// └──────────────────────────────────────────────────────────────────────────┘ + +#include "inversive_distance_functional.hpp" +#include +#include +#include + +namespace conformallab { + +/// Full finite-difference Inversive-Distance Hessian (baseline). +/// Cost: `n` full-gradient evaluations ≈ `O(n·F)`. Use for small meshes +/// or as a correctness reference for the block-FD variant. +inline Eigen::SparseMatrix inversive_distance_hessian( + const ConformalMesh& mesh, + const std::vector& x, + const InversiveDistanceMaps& m, + double eps = 1e-5) +{ + const int n = inversive_distance_dimension(mesh, m); + std::vector> trips; + trips.reserve(static_cast(n) * 16); + + std::vector xp = x, xm = x; + for (int j = 0; j < n; ++j) { + const std::size_t sj = static_cast(j); + xp[sj] = x[sj] + eps; + xm[sj] = x[sj] - eps; + + auto Gp = inversive_distance_gradient(mesh, xp, m); + auto Gm = inversive_distance_gradient(mesh, xm, m); + + xp[sj] = xm[sj] = x[sj]; // restore + + for (int i = 0; i < n; ++i) { + double val = (Gp[static_cast(i)] + - Gm[static_cast(i)]) / (2.0 * eps); + if (std::abs(val) > 1e-15) + trips.emplace_back(i, j, val); + } + } + + Eigen::SparseMatrix H(n, n); + H.setFromTriplets(trips.begin(), trips.end()); + return H; +} + +/// Symmetrised full-FD Inversive-Distance Hessian: `(H + Hᵀ)/2`. +inline Eigen::SparseMatrix inversive_distance_hessian_sym( + const ConformalMesh& mesh, + const std::vector& x, + const InversiveDistanceMaps& m, + double eps = 1e-5) +{ + auto H = inversive_distance_hessian(mesh, x, m, eps); + Eigen::SparseMatrix Ht = H.transpose(); + return (H + Ht) * 0.5; +} + +// ── Block-FD Hessian ────────────────────────────────────────────────────────── +// +// The 3 local DOFs of a face f are (u_{v1}, u_{v2}, u_{v3}). For each free +// local DOF we recompute the face's gradient contribution (−α₁,−α₂,−α₃) at +// x ± ε along that axis and read off the 3×3 Jacobian. The result scatters +// into the global Hessian via the DOF-index lookup. +// +// Cost: F × 6 face-angle evaluations (3 DOFs × 2 directions) vs n×F for +// full-FD — a speed-up of ≈ n/6, i.e. ~hundreds× on large closed meshes. +/// Per-face block-FD Inversive-Distance Hessian. Uses the locality lemma +/// `∂G_x/∂y = Σ_{f: x,y ∈ {v1,v2,v3}(f)} ∂(−α_x)/∂y` to perturb only the 3 +/// face-local DOFs at a time, giving an `F·6` face-evaluation budget vs `n·F` +/// for full-FD. Mathematically equivalent to `inversive_distance_hessian` +/// up to O(ε²) FD rounding. +inline Eigen::SparseMatrix inversive_distance_hessian_block_fd( + const ConformalMesh& mesh, + const std::vector& x, + const InversiveDistanceMaps& m, + double eps = 1e-5) +{ + const int n = inversive_distance_dimension(mesh, m); + std::vector> trips; + trips.reserve(9 * mesh.number_of_faces()); + + for (auto f : mesh.faces()) { + Halfedge_index h0 = mesh.halfedge(f); + Halfedge_index h1 = mesh.next(h0); + Halfedge_index h2 = mesh.next(h1); + + Vertex_index v1 = mesh.source(h0); + Vertex_index v2 = mesh.source(h1); + Vertex_index v3 = mesh.source(h2); + Edge_index e12 = mesh.edge(h0); + Edge_index e23 = mesh.edge(h1); + Edge_index e31 = mesh.edge(h2); + + // Local DOF indices: (u1, u2, u3). Pinned slots = -1. + const int idx[3] = { m.v_idx[v1], m.v_idx[v2], m.v_idx[v3] }; + + const double I12 = m.I_e[e12]; + const double I23 = m.I_e[e23]; + const double I31 = m.I_e[e31]; + + // Local DOF values (0 for pinned). + const double vals[3] = { + id_detail::dof_val(idx[0], x), + id_detail::dof_val(idx[1], x), + id_detail::dof_val(idx[2], x) + }; + + for (int j = 0; j < 3; ++j) { + if (idx[j] < 0) continue; // never perturb a pinned DOF + + double vp[3], vm[3]; + for (int k = 0; k < 3; ++k) { vp[k] = vm[k] = vals[k]; } + vp[j] += eps; + vm[j] -= eps; + + auto Cp = inversive_distance_face_grad_contribs( + vp[0], vp[1], vp[2], I12, I23, I31); + auto Cm = inversive_distance_face_grad_contribs( + vm[0], vm[1], vm[2], I12, I23, I31); + + const double Gp[3] = { Cp.g1, Cp.g2, Cp.g3 }; + const double Gm[3] = { Cm.g1, Cm.g2, Cm.g3 }; + + for (int i = 0; i < 3; ++i) { + if (idx[i] < 0) continue; // pinned: contributes nothing + const double val = (Gp[i] - Gm[i]) / (2.0 * eps); + if (std::abs(val) > 1e-15) + trips.emplace_back(idx[i], idx[j], val); + } + } + } + + Eigen::SparseMatrix H(n, n); + H.setFromTriplets(trips.begin(), trips.end()); + return H; +} + +/// Symmetrised block-FD Inversive-Distance Hessian: `(H + Hᵀ)/2` of +/// `inversive_distance_hessian_block_fd(...)` for solvers requiring strict +/// symmetry. +inline Eigen::SparseMatrix inversive_distance_hessian_block_fd_sym( + const ConformalMesh& mesh, + const std::vector& x, + const InversiveDistanceMaps& m, + double eps = 1e-5) +{ + auto H = inversive_distance_hessian_block_fd(mesh, x, m, eps); + Eigen::SparseMatrix Ht = H.transpose(); + return (H + Ht) * 0.5; +} + +} // namespace conformallab diff --git a/code/include/mesh_io.hpp b/code/include/mesh_io.hpp index 48bf145..84a4e40 100644 --- a/code/include/mesh_io.hpp +++ b/code/include/mesh_io.hpp @@ -29,6 +29,7 @@ #include #include #include +#include namespace conformallab { @@ -62,6 +63,12 @@ inline ConformalMesh load_mesh(const std::string& filename) throw std::runtime_error( "conformallab: mesh from " + filename + " is not triangulated (functionals require triangle faces)"); + for (auto v : mesh.vertices()) { + const auto& p = mesh.point(v); + if (!std::isfinite(p.x()) || !std::isfinite(p.y()) || !std::isfinite(p.z())) + throw std::runtime_error( + "conformallab: non-finite vertex coordinate in " + filename); + } return mesh; } diff --git a/code/include/newton_solver.hpp b/code/include/newton_solver.hpp index 25944bc..1be5041 100644 --- a/code/include/newton_solver.hpp +++ b/code/include/newton_solver.hpp @@ -36,6 +36,7 @@ #include "hyper_ideal_hessian.hpp" #include "cp_euclidean_functional.hpp" #include "inversive_distance_functional.hpp" +#include "inversive_distance_hessian.hpp" #include #include #include @@ -47,12 +48,42 @@ namespace conformallab { // ── Result ──────────────────────────────────────────────────────────────────── -/// Result of `newton_solve(...)` — converged DOF vector + diagnostics. +/// Why a Newton solve terminated (I1 audit). Distinguishes the three +/// non-convergent exits that `converged == false` previously conflated. +enum class NewtonStatus { + Converged, ///< `grad_inf_norm < tol` — success. + MaxIterations, ///< hit `max_iter` without reaching `tol`. + LinearSolverFailed, ///< both LDLT and SparseQR failed on H·Δx = −G. + LineSearchStalled ///< the line search could not reduce the residual. +}; + +/// Human-readable name for a `NewtonStatus` (for logs / test messages). +inline const char* to_string(NewtonStatus s) noexcept +{ + switch (s) { + case NewtonStatus::Converged: return "Converged"; + case NewtonStatus::MaxIterations: return "MaxIterations"; + case NewtonStatus::LinearSolverFailed: return "LinearSolverFailed"; + case NewtonStatus::LineSearchStalled: return "LineSearchStalled"; + } + return "Unknown"; +} + +/// Result of a Newton solve — converged DOF vector + diagnostics. struct NewtonResult { std::vector x; ///< DOF vector at termination. - int iterations; ///< Newton steps taken. + int iterations; ///< Newton steps actually completed. double grad_inf_norm;///< max |Gᵢ| at termination. - bool converged; ///< `true` iff `grad_inf_norm < tol`. + bool converged; ///< `true` iff `status == Converged`. + + // I1: why it stopped (more informative than the `converged` bool alone). + NewtonStatus status = NewtonStatus::MaxIterations; + + // N7: linear-algebra conditioning diagnostics. + bool sparse_qr_fallback_used = false; ///< any iteration fell back to SparseQR. + double min_ldlt_pivot = 0.0; ///< smallest |Dᵢᵢ| of the last LDLT + ///< factorisation (≈ near-singularity + ///< proxy); 0 if LDLT never succeeded. }; // ── Internal helpers ────────────────────────────────────────────────────────── @@ -103,14 +134,18 @@ inline Eigen::VectorXd solve_with_fallback( /// Solve `A·x = rhs` with the same SimplicialLDLT → SparseQR fallback /// strategy used inside all three Newton solvers. If `fallback_used` /// is non-null, it is set to `true` iff the SparseQR fallback ran. +/// If `ok` is non-null, it is set to `true` iff at least one solver succeeded. /// Returns `Eigen::VectorXd::Zero(rhs.size())` if both solvers fail. inline Eigen::VectorXd solve_linear_system( const Eigen::SparseMatrix& A, const Eigen::VectorXd& rhs, - bool* fallback_used = nullptr) + bool* fallback_used = nullptr, + bool* ok = nullptr) { - bool ok = false; - return detail::solve_with_fallback(A, rhs, ok, fallback_used); + bool ok_local = false; + auto x = detail::solve_with_fallback(A, rhs, ok_local, fallback_used); + if (ok) *ok = ok_local; + return x; } namespace detail { // re-open for the remaining helpers @@ -137,6 +172,14 @@ namespace detail { // re-open for the remaining helpers // If neither phase satisfies Armijo, return the best (smallest-residual) point // visited. If nothing beat ‖G(x)‖, return x unchanged and set *improved=false, // so the caller can stop cleanly instead of taking the old divergent full step. +// +// B2 (api-performance audit): the gradient computed at the accepted trial point +// is normally discarded (only its norm is kept), forcing the next Newton +// iteration to recompute G(x) from scratch. When `accepted_grad` is non-null it +// receives the gradient vector at the returned point, so the caller can reuse it +// as the next iteration's convergence gradient. It is only written when the +// returned point differs from `x` (i.e. `*improved == true`); on a pure stall +// (`*improved == false`) the caller must recompute G(x) itself. template inline std::vector line_search( const std::vector& x, @@ -144,9 +187,10 @@ inline std::vector line_search( const Eigen::VectorXd& d_sd, double norm0, GradFn&& grad_fn, - bool* improved = nullptr, - int max_halvings = 20, - double c1 = 1e-4) + bool* improved = nullptr, + std::vector* accepted_grad = nullptr, + int max_halvings = 20, + double c1 = 1e-4) { const int n = static_cast(x.size()); const double norm0_sq = norm0 * norm0; @@ -155,14 +199,18 @@ inline std::vector line_search( std::vector best_x = x; double best_norm = norm0; - // Evaluate ‖G(x + α·dir)‖₂, leaving the trial point in `xnew`. + std::vector trial_grad; // gradient at the most recent trial point + std::vector best_grad; // gradient at best_x (B2) + + // Evaluate ‖G(x + α·dir)‖₂, leaving the trial point in `xnew` and the + // gradient there in `trial_grad`. auto eval = [&](const Eigen::VectorXd& dir, double alpha) -> double { for (int i = 0; i < n; ++i) xnew[static_cast(i)] = x[static_cast(i)] + alpha * dir[i]; - auto Gnew = grad_fn(xnew); + trial_grad = grad_fn(xnew); double s = 0.0; - for (double v : Gnew) s += v * v; + for (double v : trial_grad) s += v * v; return std::sqrt(s); }; @@ -170,9 +218,10 @@ inline std::vector line_search( double alpha = 1.0; for (int ls = 0; ls < max_halvings; ++ls) { double norm_new = eval(dx, alpha); - if (norm_new < best_norm) { best_norm = norm_new; best_x = xnew; } + if (norm_new < best_norm) { best_norm = norm_new; best_x = xnew; best_grad = trial_grad; } if (norm_new * norm_new <= (1.0 - 2.0 * c1 * alpha) * norm0_sq) { if (improved) *improved = true; + if (accepted_grad) *accepted_grad = trial_grad; return xnew; } alpha *= 0.5; @@ -185,9 +234,10 @@ inline std::vector line_search( for (int ls = 0; ls < max_halvings; ++ls) { double thresh_sq = norm0_sq - 2.0 * c1 * alpha * dsd_sq; double norm_new = eval(d_sd, alpha); - if (norm_new < best_norm) { best_norm = norm_new; best_x = xnew; } + if (norm_new < best_norm) { best_norm = norm_new; best_x = xnew; best_grad = trial_grad; } if (thresh_sq >= 0.0 && norm_new * norm_new <= thresh_sq) { if (improved) *improved = true; + if (accepted_grad) *accepted_grad = trial_grad; return xnew; } alpha *= 0.5; @@ -196,10 +246,169 @@ inline std::vector line_search( // ── Both phases failed Armijo — never take the divergent full step. ─────── // Return the best point seen; if none improved, stay put and signal stall. - if (improved) *improved = (best_norm < norm0); + const bool any_improvement = (best_norm < norm0); + if (improved) *improved = any_improvement; + if (accepted_grad && any_improvement) *accepted_grad = best_grad; return best_x; } +// ── Cached linear solver (B4) ───────────────────────────────────────────────── +// +// solve_with_fallback rebuilds a fresh SimplicialLDLT every call, re-running the +// symbolic analysis (fill-reducing reordering / COLAMD) each Newton iteration — +// even though the Hessian sparsity pattern is iteration-invariant for a fixed +// mesh + DOF layout. This object keeps one persistent LDLT and runs +// `analyzePattern` once, then only `factorize` per iteration. +// +// FD-built Hessians (hyper-ideal, inversive-distance) prune near-zero triplets, +// so their pattern can shift slightly between iterations; we detect that via a +// change in `nonZeros()` and re-analyze, which keeps the cache correct and +// self-healing. The SparseQR fallback for singular/rank-deficient systems is +// preserved verbatim (built fresh on the rare path). The fallback semantics are +// identical to `solve_with_fallback`, so behaviour is unchanged. +struct NewtonLinearSolver { + Eigen::SimplicialLDLT> ldlt; + bool analyzed = false; + Eigen::Index last_nnz = -1; + bool fallback_used = false; ///< true iff the last solve used SparseQR + double last_min_pivot = 0.0; ///< N7: smallest |Dᵢᵢ| of the last + ///< successful LDLT (0 if it fell back) + + Eigen::VectorXd solve(const Eigen::SparseMatrix& A, + const Eigen::VectorXd& rhs, + bool& ok) + { + fallback_used = false; + last_min_pivot = 0.0; + if (!analyzed || A.nonZeros() != last_nnz) { + ldlt.analyzePattern(A); + analyzed = true; + last_nnz = A.nonZeros(); + } + ldlt.factorize(A); + if (ldlt.info() == Eigen::Success) { + Eigen::VectorXd x = ldlt.solve(rhs); + if (ldlt.info() == Eigen::Success) { + ok = true; + // N7: smallest |Dᵢᵢ| in the LDLᵀ diagonal — a cheap + // near-singularity proxy (small ⇒ ill-conditioned step). + const auto D = ldlt.vectorD(); + last_min_pivot = (D.size() > 0) ? D.cwiseAbs().minCoeff() : 0.0; + return x; + } + } + // Fallback: SparseQR — handles singular/rank-deficient A (gauge modes). + fallback_used = true; + Eigen::SparseQR, Eigen::COLAMDOrdering> qr(A); + if (qr.info() == Eigen::Success) { + Eigen::VectorXd x = qr.solve(rhs); + if (qr.info() == Eigen::Success) { ok = true; return x; } + } + ok = false; + return Eigen::VectorXd::Zero(rhs.size()); + } +}; + +// ── Unified Newton core (H2) ────────────────────────────────────────────────── +// +// All five DCE Newton solvers (Euclidean, Spherical, HyperIdeal, CP-Euclidean, +// Inversive-Distance) ran a near-identical loop that differed only in (a) the +// gradient function, (b) the Hessian function, and (c) the spherical sign quirk +// (the concave spherical energy factorises −H). This template captures that one +// loop so a fix lands once instead of five times. +// +// `grad(x) -> std::vector` : the energy gradient G(x). +// `hess(x) -> Eigen::SparseMatrix`: the TRUE Hessian H(x) (un-negated). +// `concave` : if true, solve (−H)·Δx = G (spherical); +// otherwise H·Δx = −G. Both are the same +// Newton system; `concave` only selects +// the PSD matrix actually factorised. +// +// Folds in B2 (reuse the line-search gradient as the next convergence gradient), +// B4 (cached symbolic factorization), and B5 (no dead solver variable). +template +inline NewtonResult newton_core( + std::vector x, + GradFn&& grad, + HessFn&& hess, + bool concave, + double tol, + int max_iter) +{ + const int n = static_cast(x.size()); + + NewtonResult res; + res.converged = false; + res.iterations = 0; + res.grad_inf_norm = 0.0; + + NewtonLinearSolver solver; + + // B2: gradient carried across the loop; the first one is the only "extra" + // evaluation — every later iteration reuses the line-search gradient. + std::vector G_std = grad(x); + + for (int iter = 0; iter < max_iter; ++iter) { + Eigen::Map G(G_std.data(), n); + + const double inf_norm = G.cwiseAbs().maxCoeff(); + if (inf_norm < tol) { + res.converged = true; + res.status = NewtonStatus::Converged; // I1 + res.grad_inf_norm = inf_norm; + res.iterations = iter; + res.x = x; + return res; + } + + Eigen::SparseMatrix H = hess(x); + + // Newton system: factor the PSD matrix (H for convex, −H for concave). + Eigen::SparseMatrix A = concave ? Eigen::SparseMatrix(-H) : H; + Eigen::VectorXd rhs = concave ? Eigen::VectorXd(G) : Eigen::VectorXd(-G); + + bool ok = false; + Eigen::VectorXd dx = solver.solve(A, rhs, ok); + // N7: accumulate conditioning diagnostics from this solve. + res.sparse_qr_fallback_used = res.sparse_qr_fallback_used || solver.fallback_used; + res.min_ldlt_pivot = solver.last_min_pivot; + if (!ok) { // I1/H1: linear solver failed this step + res.iterations = iter; // (iter steps actually completed) + res.status = NewtonStatus::LinearSolverFailed; + break; + } + + // Globalised line search (Armijo + steepest-descent fallback). + // d_sd = −(H·G) is the merit-function steepest descent for f = ½‖G‖², + // using the TRUE (un-negated) Hessian for both signs. + const double norm0 = G.norm(); + Eigen::VectorXd d_sd = -(H * G); + bool improved = true; + std::vector G_next; + x = detail::line_search(x, dx, d_sd, norm0, grad, &improved, &G_next); + if (!improved) { // I1/H1: line search stalled — stop cleanly + res.iterations = iter; + res.status = NewtonStatus::LineSearchStalled; + break; + } + + G_std = std::move(G_next); // B2: reuse for the next convergence check + res.iterations = iter + 1; + } + + // If we fell out of the loop without converging or breaking, it was max_iter + // (the default `res.status` is already MaxIterations; a break has set its own). + + // Final gradient norm (reached only on max-iter / break — recomputed to be + // correct after a stalled or failed step). + auto G_final = grad(x); + double inf_final = 0.0; + for (double v : G_final) inf_final = std::max(inf_final, std::abs(v)); + res.grad_inf_norm = inf_final; + res.x = x; + return res; +} + } // namespace detail // ── Euclidean Newton solver ──────────────────────────────────────────────────── @@ -233,63 +442,21 @@ inline NewtonResult newton_euclidean( double tol = 1e-8, int max_iter = 200) { - std::vector x = x0; - const int n = static_cast(x.size()); + // Layout is loop-invariant — decide the Hessian variant once (B3): cyclic + // layout (edge DOFs present) → analytic cyclic Hessian, which covers the + // vertex-edge / edge-edge blocks the vertex-only cotangent Laplacian omits. + bool has_edge_dof = false; + for (auto e : mesh.edges()) + if (m.e_idx[e] >= 0) { has_edge_dof = true; break; } - NewtonResult res; - res.converged = false; - res.iterations = 0; - res.grad_inf_norm = 0.0; - - Eigen::SimplicialLDLT> solver; - - for (int iter = 0; iter < max_iter; ++iter) { - // ── Gradient ────────────────────────────────────────────────────────── - auto G_std = euclidean_gradient(mesh, x, m); - Eigen::Map G(G_std.data(), n); - - double inf_norm = G.cwiseAbs().maxCoeff(); - if (inf_norm < tol) { - res.converged = true; - res.grad_inf_norm = inf_norm; - res.iterations = iter; - res.x = x; - return res; - } - - // ── Hessian + solve H·Δx = −G (SparseQR fallback for singular H) ── - // Cyclic layout (edge DOFs present) → analytic cyclic Hessian, covering - // the vertex-edge / edge-edge blocks the vertex-only cotangent Laplacian - // omits. Vertex-only layout → analytic cotangent Laplacian. - bool has_edge_dof = false; - for (auto e : mesh.edges()) - if (m.e_idx[e] >= 0) { has_edge_dof = true; break; } - auto H = has_edge_dof ? euclidean_hessian_analytic(mesh, x, m) - : euclidean_hessian(mesh, x, m); - bool ok = false; - Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok); - if (!ok) break; - - // ── Globalised line search (Armijo + steepest-descent fallback) ─────── - double norm0 = G.norm(); - Eigen::VectorXd d_sd = -(H * G); // merit-function steepest descent - bool improved = true; - x = detail::line_search(x, dx, d_sd, norm0, - [&](const std::vector& xnew) { - return euclidean_gradient(mesh, xnew, m); - }, &improved); - if (!improved) break; // line search stalled — stop cleanly - - res.iterations = iter + 1; - } - - // Report final gradient norm - auto G_final = euclidean_gradient(mesh, x, m); - double inf_final = 0.0; - for (double v : G_final) inf_final = std::max(inf_final, std::abs(v)); - res.grad_inf_norm = inf_final; - res.x = x; - return res; + return detail::newton_core( + std::move(x0), + [&](const std::vector& xc) { return euclidean_gradient(mesh, xc, m); }, + [&](const std::vector& xc) { + return has_edge_dof ? euclidean_hessian_analytic(mesh, xc, m) + : euclidean_hessian(mesh, xc, m); + }, + /*concave=*/false, tol, max_iter); } // ── Spherical Newton solver ─────────────────────────────────────────────────── @@ -322,55 +489,14 @@ inline NewtonResult newton_spherical( double tol = 1e-8, int max_iter = 200) { - std::vector x = x0; - const int n = static_cast(x.size()); - - NewtonResult res; - res.converged = false; - res.iterations = 0; - res.grad_inf_norm = 0.0; - - for (int iter = 0; iter < max_iter; ++iter) { - // ── Gradient ────────────────────────────────────────────────────────── - auto G_std = spherical_gradient(mesh, x, m); - Eigen::Map G(G_std.data(), n); - - double inf_norm = G.cwiseAbs().maxCoeff(); - if (inf_norm < tol) { - res.converged = true; - res.grad_inf_norm = inf_norm; - res.iterations = iter; - res.x = x; - return res; - } - - // ── Hessian: negate to get PSD; solve (−H)·Δx = G ────────────────── - auto H = spherical_hessian(mesh, x, m); - auto negH = Eigen::SparseMatrix(-H); - bool ok = false; - Eigen::VectorXd dx = detail::solve_with_fallback(negH, G, ok); - if (!ok) break; - - // ── Globalised line search (Armijo + steepest-descent fallback) ─────── - // d_sd uses the actual (un-negated) Hessian: ∇f = H·G for f = ½‖G‖². - double norm0 = G.norm(); - Eigen::VectorXd d_sd = -(H * G); - bool improved = true; - x = detail::line_search(x, dx, d_sd, norm0, - [&](const std::vector& xnew) { - return spherical_gradient(mesh, xnew, m); - }, &improved); - if (!improved) break; - - res.iterations = iter + 1; - } - - auto G_final = spherical_gradient(mesh, x, m); - double inf_final = 0.0; - for (double v : G_final) inf_final = std::max(inf_final, std::abs(v)); - res.grad_inf_norm = inf_final; - res.x = x; - return res; + // Concave energy: the Hessian H is NSD, so newton_core factors −H (PSD) and + // solves (−H)·Δx = G — the same Newton system, just the sign that keeps LDLT + // well-defined. The merit steepest-descent inside still uses the true H. + return detail::newton_core( + std::move(x0), + [&](const std::vector& xc) { return spherical_gradient(mesh, xc, m); }, + [&](const std::vector& xc) { return spherical_hessian(mesh, xc, m); }, + /*concave=*/true, tol, max_iter); } // ── HyperIdeal Newton solver ────────────────────────────────────────────────── @@ -395,6 +521,12 @@ inline NewtonResult newton_spherical( /// \param max_iter Maximum Newton iterations. Default: 200. /// \param hess_eps Finite-difference step for Hessian approximation. Default: 1e-5. /// (Phase 9b will replace this with an analytic Hessian.) +/// \param clamp Vertex-scale floor mode (N3 audit). Default `HardJava` +/// reproduces the Java oracle's hard `b<0 → floor` snap +/// (C⁰, parity-faithful). `SmoothBarrier` uses the C¹ +/// softplus floor — smoother near the feasibility boundary +/// but deviates from the Java golden values. See +/// `HyperIdealScaleClamp` in hyper_ideal_functional.hpp. /// \return NewtonResult{x*, iterations, grad_inf_norm, converged}. /// /// \see Springborn (2020), Theorem 1.3 for the strict convexity proof. @@ -405,55 +537,20 @@ inline NewtonResult newton_hyper_ideal( const HyperIdealMaps& m, double tol = 1e-8, int max_iter = 200, - double hess_eps = 1e-5) + double hess_eps = 1e-5, + HyperIdealScaleClamp clamp = HyperIdealScaleClamp::HardJava) { - std::vector x = x0; - const int n = static_cast(x.size()); - - NewtonResult res; - res.converged = false; - res.iterations = 0; - res.grad_inf_norm = 0.0; - - for (int iter = 0; iter < max_iter; ++iter) { - // ── Gradient ────────────────────────────────────────────────────────── - auto G_std = evaluate_hyper_ideal(mesh, x, m, /*energy=*/false).gradient; - Eigen::Map G(G_std.data(), n); - - double inf_norm = G.cwiseAbs().maxCoeff(); - if (inf_norm < tol) { - res.converged = true; - res.grad_inf_norm = inf_norm; - res.iterations = iter; - res.x = x; - return res; - } - - // ── Hessian (numerical FD) + solve H·Δx = −G ───────────────────────── - auto H = hyper_ideal_hessian_sym(mesh, x, m, hess_eps); - bool ok = false; - Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok); - if (!ok) break; - - // ── Globalised line search (Armijo + steepest-descent fallback) ─────── - double norm0 = G.norm(); - Eigen::VectorXd d_sd = -(H * G); - bool improved = true; - x = detail::line_search(x, dx, d_sd, norm0, - [&](const std::vector& xnew) { - return evaluate_hyper_ideal(mesh, xnew, m, false).gradient; - }, &improved); - if (!improved) break; - - res.iterations = iter + 1; - } - - auto G_final = evaluate_hyper_ideal(mesh, x, m, false).gradient; - double inf_final = 0.0; - for (double v : G_final) inf_final = std::max(inf_final, std::abs(v)); - res.grad_inf_norm = inf_final; - res.x = x; - return res; + // block-FD Hessian (~33–1166× faster than full-FD); `clamp` selects the + // vertex-scale floor mode (N3). Convex energy → factor H directly. + return detail::newton_core( + std::move(x0), + [&](const std::vector& xc) { + return evaluate_hyper_ideal(mesh, xc, m, /*energy=*/false, /*grad=*/true, clamp).gradient; + }, + [&](const std::vector& xc) { + return hyper_ideal_hessian_block_fd_sym(mesh, xc, m, hess_eps, clamp); + }, + /*concave=*/false, tol, max_iter); } // ── CP-Euclidean Newton solver (Phase 9a.1) ─────────────────────────────────── @@ -461,7 +558,7 @@ inline NewtonResult newton_hyper_ideal( /// Solve the CP-Euclidean circle-packing problem: find ρ ∈ ℝ^F such that the /// per-face angle sums match φ_f at every free face. /// -/// The CP-Euclidean energy (Bobenko-Pinkall-Springborn 2010 §6) is strictly +/// The CP-Euclidean energy (Bobenko-Pinkall-Springborn 2015 §6) is strictly /// convex on its open domain of validity, so the Hessian H is PSD and the /// solution is unique up to the gauge mode pinned by `f_idx == −1`. /// `cp_euclidean_hessian` provides the analytic 2×2-per-edge formula @@ -479,7 +576,7 @@ inline NewtonResult newton_hyper_ideal( /// \note Unlike the Euclidean solver, the CP-Euclidean Hessian is exact /// (analytic), so the SparseQR fallback only triggers in genuine /// gauge-singular situations (no pinned face). -/// \see doc/architecture/phase-9a-validation.md §1 for the BPS-2010 mapping. +/// \see doc/architecture/phase-9a-validation.md §1 for the BPS-2015 mapping. inline NewtonResult newton_cp_euclidean( ConformalMesh& mesh, std::vector x0, @@ -487,50 +584,12 @@ inline NewtonResult newton_cp_euclidean( double tol = 1e-8, int max_iter = 200) { - std::vector x = x0; - const int n = static_cast(x.size()); - - NewtonResult res; - res.converged = false; - res.iterations = 0; - res.grad_inf_norm = 0.0; - - for (int iter = 0; iter < max_iter; ++iter) { - auto G_std = cp_euclidean_gradient(mesh, x, m); - Eigen::Map G(G_std.data(), n); - - double inf_norm = G.cwiseAbs().maxCoeff(); - if (inf_norm < tol) { - res.converged = true; - res.grad_inf_norm = inf_norm; - res.iterations = iter; - res.x = x; - return res; - } - - auto H = cp_euclidean_hessian(mesh, x, m); - bool ok = false; - Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok); - if (!ok) break; - - double norm0 = G.norm(); - Eigen::VectorXd d_sd = -(H * G); - bool improved = true; - x = detail::line_search(x, dx, d_sd, norm0, - [&](const std::vector& xnew) { - return cp_euclidean_gradient(mesh, xnew, m); - }, &improved); - if (!improved) break; - - res.iterations = iter + 1; - } - - auto G_final = cp_euclidean_gradient(mesh, x, m); - double inf_final = 0.0; - for (double v : G_final) inf_final = std::max(inf_final, std::abs(v)); - res.grad_inf_norm = inf_final; - res.x = x; - return res; + // Convex energy with an exact analytic Hessian (BPS-2015 2×2-per-edge). + return detail::newton_core( + std::move(x0), + [&](const std::vector& xc) { return cp_euclidean_gradient(mesh, xc, m); }, + [&](const std::vector& xc) { return cp_euclidean_hessian(mesh, xc, m); }, + /*concave=*/false, tol, max_iter); } // ── Inversive-Distance Newton solver (Phase 9a.2) ───────────────────────────── @@ -542,10 +601,11 @@ inline NewtonResult newton_cp_euclidean( /// domain where every triangle satisfies the inequalities. Luo's 1-form is /// closed there, so the path-integral energy is well-defined. /// -/// MVP implementation: the Hessian is computed by **finite differences** of -/// the analytic gradient (same pattern as the Phase 4a HyperIdeal solver). -/// An analytic Hessian via Glickenstein 2011 eq. (4.6) is tracked in -/// `doc/roadmap/research-track.md` as Phase 9a.2-analytic. +/// The Hessian is computed by **per-face block finite differences** +/// (`inversive_distance_hessian_block_fd_sym`, same pattern as the HyperIdeal +/// solver after Phase 9b) — O(F) face evaluations instead of the O(n·F) of the +/// full-FD baseline. An analytic Hessian via Glickenstein 2011 eq. (4.6) is +/// tracked in `doc/roadmap/research-track.md` as Phase 9a.2-analytic. /// /// \param mesh Input triangle mesh. /// \param x0 Initial DOF vector (length = number of free vertices). @@ -568,81 +628,14 @@ inline NewtonResult newton_inversive_distance( int max_iter = 200, double hess_eps = 1e-5) { - std::vector x = x0; - const int n = static_cast(x.size()); - - NewtonResult res; - res.converged = false; - res.iterations = 0; - res.grad_inf_norm = 0.0; - - // Local FD Hessian builder — n × (cost of gradient eval). - auto build_hessian = [&](const std::vector& xc) -> Eigen::SparseMatrix { - std::vector> trips; - trips.reserve(static_cast(n) * 16); // sparse heuristic - - std::vector xp = xc, xm = xc; - for (int j = 0; j < n; ++j) { - const std::size_t sj = static_cast(j); - xp[sj] = xc[sj] + hess_eps; - xm[sj] = xc[sj] - hess_eps; - - auto Gp = inversive_distance_gradient(mesh, xp, m); - auto Gm = inversive_distance_gradient(mesh, xm, m); - - xp[sj] = xm[sj] = xc[sj]; // restore - - for (int i = 0; i < n; ++i) { - double val = (Gp[static_cast(i)] - - Gm[static_cast(i)]) - / (2.0 * hess_eps); - if (std::abs(val) > 1e-15) - trips.emplace_back(i, j, val); - } - } - Eigen::SparseMatrix H(n, n); - H.setFromTriplets(trips.begin(), trips.end()); - // Symmetrise — FD rounding may introduce tiny asymmetries. - Eigen::SparseMatrix Ht = H.transpose(); - return (H + Ht) * 0.5; - }; - - for (int iter = 0; iter < max_iter; ++iter) { - auto G_std = inversive_distance_gradient(mesh, x, m); - Eigen::Map G(G_std.data(), n); - - double inf_norm = G.cwiseAbs().maxCoeff(); - if (inf_norm < tol) { - res.converged = true; - res.grad_inf_norm = inf_norm; - res.iterations = iter; - res.x = x; - return res; - } - - auto H = build_hessian(x); - bool ok = false; - Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok); - if (!ok) break; - - double norm0 = G.norm(); - Eigen::VectorXd d_sd = -(H * G); - bool improved = true; - x = detail::line_search(x, dx, d_sd, norm0, - [&](const std::vector& xnew) { - return inversive_distance_gradient(mesh, xnew, m); - }, &improved); - if (!improved) break; - - res.iterations = iter + 1; - } - - auto G_final = inversive_distance_gradient(mesh, x, m); - double inf_final = 0.0; - for (double v : G_final) inf_final = std::max(inf_final, std::abs(v)); - res.grad_inf_norm = inf_final; - res.x = x; - return res; + // Convex energy; block-FD Hessian (~n/6× faster than full-FD). + return detail::newton_core( + std::move(x0), + [&](const std::vector& xc) { return inversive_distance_gradient(mesh, xc, m); }, + [&](const std::vector& xc) { + return inversive_distance_hessian_block_fd_sym(mesh, xc, m, hess_eps); + }, + /*concave=*/false, tol, max_iter); } } // namespace conformallab diff --git a/code/include/serialization.hpp b/code/include/serialization.hpp index 6e87376..4265deb 100644 --- a/code/include/serialization.hpp +++ b/code/include/serialization.hpp @@ -93,31 +93,60 @@ inline std::vector load_result_json( { using json = nlohmann::json; std::ifstream ifs(path); - if (!ifs) throw std::runtime_error("Cannot open: " + path); - json j; ifs >> j; + if (!ifs) throw std::runtime_error("conformallab: cannot open: " + path); - if (geom && j.contains("geometry")) - *geom = j["geometry"].get(); - - std::vector x = j.at("dof_vector").get>(); - - if (res) { - res->x = x; - res->converged = j["solver"]["converged"].get(); - res->iterations = j["solver"]["iterations"].get(); - res->grad_inf_norm = j["solver"]["grad_inf_norm"].get(); + // V1: wrap parse + field extraction so nlohmann exceptions (parse_error, + // type_error, out_of_range) surface as std::runtime_error with the path. + json j; + try { + ifs >> j; + } catch (const json::exception& e) { + throw std::runtime_error( + "conformallab: malformed JSON in " + path + ": " + e.what()); } - if (layout2d && j.contains("layout") && j["layout"]["dim"] == 2) { - auto uv = j["layout"]["uv"]; - layout2d->uv.resize(uv.size()); - for (std::size_t i = 0; i < uv.size(); ++i) - layout2d->uv[i] = { uv[i][0].get(), uv[i][1].get() }; - layout2d->has_seam = j["layout"].value("has_seam", false); - layout2d->success = true; - } + try { + if (geom && j.contains("geometry")) + *geom = j["geometry"].get(); - return x; + // V4: validate required top-level key before accessing it. + if (!j.contains("dof_vector")) + throw std::runtime_error( + "conformallab: result JSON missing field 'dof_vector' in " + path); + std::vector x = j.at("dof_vector").get>(); + + if (res) { + // V4: validate nested solver keys before accessing. + if (!j.contains("solver")) + throw std::runtime_error( + "conformallab: result JSON missing field 'solver' in " + path); + const auto& s = j.at("solver"); + for (const char* key : {"converged", "iterations", "grad_inf_norm"}) { + if (!s.contains(key)) + throw std::runtime_error( + std::string("conformallab: result JSON missing field 'solver.") + + key + "' in " + path); + } + res->x = x; + res->converged = s.at("converged").get(); + res->iterations = s.at("iterations").get(); + res->grad_inf_norm = s.at("grad_inf_norm").get(); + } + + if (layout2d && j.contains("layout") && j["layout"]["dim"] == 2) { + auto uv = j["layout"]["uv"]; + layout2d->uv.resize(uv.size()); + for (std::size_t i = 0; i < uv.size(); ++i) + layout2d->uv[i] = { uv[i][0].get(), uv[i][1].get() }; + layout2d->has_seam = j["layout"].value("has_seam", false); + layout2d->success = true; + } + + return x; + } catch (const json::exception& e) { + throw std::runtime_error( + "conformallab: malformed JSON in " + path + ": " + e.what()); + } } // ════════════════════════════════════════════════════════════════════════════ @@ -255,9 +284,23 @@ inline std::vector load_result_xml( // Solver metadata else if (line.find("converged = (detail_xml::xml_get_attr(line, "converged") == "true"); - 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")); + res->converged = (detail_xml::xml_get_attr(line, "converged") == "true"); + // V2: stoi/stod throw std::invalid_argument on empty or non-numeric + // attribute values; wrap and rethrow as runtime_error with context. + try { + auto iter_str = detail_xml::xml_get_attr(line, "iterations"); + auto grad_str = detail_xml::xml_get_attr(line, "grad_inf_norm"); + if (iter_str.empty()) + throw std::runtime_error("missing attribute 'iterations'"); + if (grad_str.empty()) + throw std::runtime_error("missing attribute 'grad_inf_norm'"); + res->iterations = std::stoi(iter_str); + res->grad_inf_norm = std::stod(grad_str); + } catch (const std::exception& e) { + throw std::runtime_error( + "conformallab: malformed XML Solver element in " + + path + ": " + e.what()); + } } } // DOF vector diff --git a/code/include/spherical_functional.hpp b/code/include/spherical_functional.hpp index adbe623..3cb69be 100644 --- a/code/include/spherical_functional.hpp +++ b/code/include/spherical_functional.hpp @@ -35,6 +35,7 @@ // └──────────────────────────────────────────────────────────────────────────┘ #include "conformal_mesh.hpp" +#include "constants.hpp" #include "spherical_geometry.hpp" #include "gauss_legendre.hpp" #include @@ -176,7 +177,7 @@ inline void compute_spherical_lambda0_from_mesh(ConformalMesh& mesh, SphericalMa if (half_sin > 1e-15) m.lambda0[e] = 2.0 * std::log(half_sin); else - m.lambda0[e] = -30.0; // very short edge: essentially 0 + m.lambda0[e] = LOG_EDGE_LENGTH_FLOOR; // very short edge: essentially 0 } } diff --git a/code/include/spherical_geometry.hpp b/code/include/spherical_geometry.hpp index 991ebb7..1702f09 100644 --- a/code/include/spherical_geometry.hpp +++ b/code/include/spherical_geometry.hpp @@ -31,7 +31,7 @@ constexpr double PI_SPHER = PI; inline double spherical_l(double lambda) { double half = std::exp(lambda * 0.5); - if (half >= 1.0) half = 1.0 - 1e-15; + if (half >= 1.0) half = ASIN_DOMAIN_GUARD; if (half <= 0.0) return 0.0; return 2.0 * std::asin(half); } diff --git a/code/tests/cgal/test_cgal_traits_mvp.cpp b/code/tests/cgal/test_cgal_traits_mvp.cpp index 5767502..36e1ac7 100644 --- a/code/tests/cgal/test_cgal_traits_mvp.cpp +++ b/code/tests/cgal/test_cgal_traits_mvp.cpp @@ -111,6 +111,11 @@ TEST(CGALDiscreteConformalMap, SingleTriangleConverges) EXPECT_GE(result.iterations, 0); EXPECT_EQ(result.u_per_vertex.size(), num_vertices(mesh)); + // I1/N7 propagated through the CGAL layer: status enum + diagnostics. + EXPECT_EQ(result.status, CGAL::Newton_status::Converged); + EXPECT_FALSE(result.sparse_qr_fallback_used); // pinned vertex → no gauge mode + EXPECT_GE(result.min_ldlt_pivot, 0.0); + // With the default flat-disc target curvature and the first vertex pinned, // the natural-theta equilibrium is at u = 0 — Newton should accept x0=0. for (double u : result.u_per_vertex) diff --git a/code/tests/cgal/test_euclidean_hessian.cpp b/code/tests/cgal/test_euclidean_hessian.cpp index 80305c6..72e7ca6 100644 --- a/code/tests/cgal/test_euclidean_hessian.cpp +++ b/code/tests/cgal/test_euclidean_hessian.cpp @@ -61,6 +61,47 @@ TEST(EuclideanHessian, CotWeights_RightIsoscelesTriangle) EXPECT_NEAR(cw.cot3, 1.0, 1e-12); // 45° at v3 } +// ════════════════════════════════════════════════════════════════════════════ +// N5: cotangent weights on a SLIVER triangle match a high-precision reference. +// +// The area is now computed by Kahan's stable side-length formula instead of +// the naive 2·√(t12·t23·t31·l123). On a thin (sliver) triangle the naive +// product of near-equal-length differences loses precision, biasing every +// cotangent weight. Here we cross-check against an independent law-of-cosines +// reference in long double (80-bit on x86 CI — a genuine high-precision oracle; +// equal to double on ARM64, where it still serves as a regression cross-check). +// ════════════════════════════════════════════════════════════════════════════ + +TEST(EuclideanHessian, CotWeights_SliverMatchesHighPrecisionReference) +{ + // Thin isosceles: short base l12, two near-unit legs (apex angle ≈ 0.01 rad). + const double l12 = 0.01, l23 = 1.0, l31 = 1.0; + auto cw = euclidean_cot_weights(l12, l23, l31); + ASSERT_TRUE(cw.valid); + + // cot at the vertex opposite side `opp`, with adjacent sides s1, s2: + // cos = (s1² + s2² − opp²) / (2·s1·s2), sin = √((1−cos)(1+cos)). + auto cot_ref = [](long double opp, long double s1, long double s2) { + long double cosA = (s1 * s1 + s2 * s2 - opp * opp) / (2.0L * s1 * s2); + long double sinA = std::sqrt((1.0L - cosA) * (1.0L + cosA)); + return static_cast(cosA / sinA); + }; + const double c1 = cot_ref(l23, l12, l31); // v1 opposite l23 + const double c2 = cot_ref(l31, l12, l23); // v2 opposite l31 + const double c3 = cot_ref(l12, l23, l31); // v3 opposite l12 (needle angle) + + auto rel = [](double a, double b) { + return std::abs(a - b) / std::max(1.0, std::abs(b)); + }; + EXPECT_LT(rel(cw.cot1, c1), 1e-9); + EXPECT_LT(rel(cw.cot2, c2), 1e-9); + EXPECT_LT(rel(cw.cot3, c3), 1e-9); + + EXPECT_TRUE(std::isfinite(cw.cot1) && + std::isfinite(cw.cot2) && + std::isfinite(cw.cot3)); +} + // ════════════════════════════════════════════════════════════════════════════ // Hessian is symmetric: H[i,j] == H[j,i] // ════════════════════════════════════════════════════════════════════════════ diff --git a/code/tests/cgal/test_hyper_ideal_hessian.cpp b/code/tests/cgal/test_hyper_ideal_hessian.cpp index 8a9043f..3722134 100644 --- a/code/tests/cgal/test_hyper_ideal_hessian.cpp +++ b/code/tests/cgal/test_hyper_ideal_hessian.cpp @@ -338,3 +338,81 @@ TEST(HyperIdealHessian, BlockFD_FasterThanFullFD) EXPECT_GE(ms_full, 3 * ms_block) << "Block-FD should be at least 3× faster than full-FD on this mesh"; } + +// ════════════════════════════════════════════════════════════════════════════ +// N3 — Scale-floor clamp modes (HardJava vs SmoothBarrier) +// +// The vertex-scale floor can be applied two ways (HyperIdealScaleClamp): +// • HardJava (default): b<0 → floor. Faithful to Java, but only C⁰ — and +// in fact value-discontinuous at b=0 (jumps from `floor` to 0). +// • SmoothBarrier: floor + softplus_β(b−floor). C¹ everywhere, ≈ b away +// from the floor, and keeps b ≥ floor > 0 smoothly. +// These tests pin the contract of both modes and the default-preservation. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(HyperIdealClamp, SmoothBarrierIsContinuousAndC1AcrossZero) +{ + const auto Hard = HyperIdealScaleClamp::HardJava; + const auto Smooth = HyperIdealScaleClamp::SmoothBarrier; + auto f = [](double b, HyperIdealScaleClamp mode) { + return clamp_hyper_ideal_scale(b, /*variable=*/true, mode); + }; + + const double eps = 1e-6; + + // HardJava jumps in VALUE at b=0 (floor on the left, 0 on the right). + EXPECT_GT(std::abs(f(-eps, Hard) - f(+eps, Hard)), 0.5 * HYPER_IDEAL_SCALE_FLOOR); + + // SmoothBarrier is continuous in value across b=0 … + EXPECT_NEAR(f(-eps, Smooth), f(+eps, Smooth), 1e-4); + + // … and C¹: the one-sided slopes across b=0 agree (the HardJava slopes, + // 0 on the left and 1 on the right, would not). + auto slope = [&](double b0, HyperIdealScaleClamp mode) { + return (f(b0 + eps, mode) - f(b0 - eps, mode)) / (2.0 * eps); + }; + const double sL = slope(-1e-3, Smooth); + const double sR = slope(+1e-3, Smooth); + EXPECT_NEAR(sL, sR, 5e-2) << "SmoothBarrier derivative should be continuous"; +} + +TEST(HyperIdealClamp, SmoothBarrierStaysAboveFloorAndApproachesIdentity) +{ + const auto Smooth = HyperIdealScaleClamp::SmoothBarrier; + auto f = [&](double b) { return clamp_hyper_ideal_scale(b, true, Smooth); }; + + // At or above the floor for any input (mathematically > floor; for very + // negative b the softplus underflows to 0 in double, giving exactly floor). + EXPECT_GE(f(-1e6), HYPER_IDEAL_SCALE_FLOOR); + EXPECT_GE(f(-1.0), HYPER_IDEAL_SCALE_FLOOR); + EXPECT_GT(f(0.0), 0.0); + EXPECT_LT(f(-1e6) - HYPER_IDEAL_SCALE_FLOOR, 1e-9); // converges down to floor + + // ≈ identity well above the floor (the normal operating regime). + EXPECT_NEAR(f(1.0), 1.0, 1e-9); + EXPECT_NEAR(f(5.0), 5.0, 1e-12); + + // Pinned DOFs are never clamped, regardless of mode. + EXPECT_EQ(clamp_hyper_ideal_scale(-3.0, /*variable=*/false, Smooth), -3.0); +} + +// Away from the feasibility boundary (all b = 1 ≫ floor), the two modes must +// produce the same Hessian — SmoothBarrier only differs near b ≈ floor, so the +// default-mode parity results are not perturbed in the normal regime. +TEST(HyperIdealClamp, ModesAgreeAwayFromBoundary) +{ + auto mesh = make_tetrahedron(); + auto m = setup_hyper_ideal_maps(mesh); + const int n = assign_all_dof_indices(mesh, m); + auto x = natural_x(mesh, m); // b = 1, a = 0.5 — all well above the floor + + auto H_hard = hyper_ideal_hessian_block_fd_sym( + mesh, x, m, 1e-5, HyperIdealScaleClamp::HardJava); + auto H_soft = hyper_ideal_hessian_block_fd_sym( + mesh, x, m, 1e-5, HyperIdealScaleClamp::SmoothBarrier); + + Eigen::MatrixXd Dh(H_hard), Ds(H_soft); + EXPECT_LT((Dh - Ds).cwiseAbs().maxCoeff(), 1e-9) + << "clamp modes must agree when every scale is well above the floor"; + (void)n; +} diff --git a/code/tests/cgal/test_lawson_hyperideal.cpp b/code/tests/cgal/test_lawson_hyperideal.cpp index 74b8290..539f386 100644 --- a/code/tests/cgal/test_lawson_hyperideal.cpp +++ b/code/tests/cgal/test_lawson_hyperideal.cpp @@ -194,6 +194,54 @@ TEST(LawsonHyperIdeal, ConvergenceGoldenVector_JavaXVal) } } +// ════════════════════════════════════════════════════════════════════════════ +// N3 — SmoothBarrier clamp mode converges to the same golden solution +// +// The C¹ smooth-barrier scale floor is an opt-in alternative to the Java hard +// clamp. On this well-posed problem the scales stay well above the floor +// throughout the solve, so the barrier is ≈ identity and the converged vector +// must match the same golden classes as the HardJava run above — demonstrating +// the smooth mode is a drop-in that does not move the solution when the clamp +// region is never entered. +// ════════════════════════════════════════════════════════════════════════════ +TEST(LawsonHyperIdeal, SmoothBarrierConvergesToGoldenVector) +{ + std::set original; + ConformalMesh m = make_lawson_square_tiled(&original); + + HyperIdealMaps maps = setup_hyper_ideal_maps(m); + const int n = assign_all_dof_indices(m, maps); + ASSERT_EQ(n, 4 + 18); + + for (auto e : m.edges()) + if (original.count(e)) maps.theta_e[e] = PI / 2.0; + + std::vector x0(static_cast(n), 1.0); + auto res = newton_hyper_ideal(m, x0, maps, /*tol=*/1e-10, /*max_iter=*/200, + /*hess_eps=*/1e-5, + HyperIdealScaleClamp::SmoothBarrier); + ASSERT_TRUE(res.converged) + << "SmoothBarrier HyperIdeal Newton did not converge; ||G||=" + << res.grad_inf_norm; + + constexpr double b_gold = 1.1462158341786262; + constexpr double a_orig_gold = 1.7627471737467797; + constexpr double a_aux_gold = 2.633915794495759; + const double tol = 1e-5; + + for (auto v : m.vertices()) { + const int iv = maps.v_idx[v]; + ASSERT_GE(iv, 0); + EXPECT_NEAR(res.x[static_cast(iv)], b_gold, tol); + } + for (auto e : m.edges()) { + const int ie = maps.e_idx[e]; + ASSERT_GE(ie, 0); + const double gold = original.count(e) ? a_orig_gold : a_aux_gold; + EXPECT_NEAR(res.x[static_cast(ie)], gold, tol); + } +} + // ════════════════════════════════════════════════════════════════════════════ // Branch-points variant — Java HyperIdealConvergenceTest...WithBranchPoints // diff --git a/code/tests/cgal/test_layout.cpp b/code/tests/cgal/test_layout.cpp index 82b8692..bce1528 100644 --- a/code/tests/cgal/test_layout.cpp +++ b/code/tests/cgal/test_layout.cpp @@ -24,6 +24,7 @@ #include "serialization.hpp" #include #include +#include #include #include @@ -370,3 +371,67 @@ TEST(Serialization, XML_RoundTrip) std::filesystem::remove(path); } + +// ════════════════════════════════════════════════════════════════════════════ +// I2: serialization throws on malformed / schema-violating input +// +// These tests cover the throw paths in load_result_json / load_result_xml +// that were previously untested (I2 in the test-coverage audit). +// ════════════════════════════════════════════════════════════════════════════ + +TEST(Serialization, LoadResultJson_ThrowsOnMissingFile) +{ + EXPECT_THROW(load_result_json("/tmp/does_not_exist_conflab.json"), + std::runtime_error); +} + +TEST(Serialization, LoadResultJson_ThrowsOnMalformedJson) +{ + const std::string path = "/tmp/conflab_bad.json"; + { std::ofstream ofs(path); ofs << "{not valid json!!!"; } + EXPECT_THROW(load_result_json(path), std::runtime_error); + std::filesystem::remove(path); +} + +TEST(Serialization, LoadResultJson_ThrowsOnMissingDofVector) +{ + // V4: a valid JSON object but without the required "dof_vector" key. + const std::string path = "/tmp/conflab_nodof.json"; + { std::ofstream ofs(path); ofs << R"({"geometry":"euclidean"})"; } + EXPECT_THROW(load_result_json(path), std::runtime_error); + std::filesystem::remove(path); +} + +TEST(Serialization, LoadResultJson_ThrowsOnMissingSolverField) +{ + // V4: has dof_vector but solver block is absent when res != nullptr. + const std::string path = "/tmp/conflab_nosolver.json"; + { std::ofstream ofs(path); ofs << R"({"dof_vector":[0.1,0.2]})"; } + NewtonResult res; + EXPECT_THROW(load_result_json(path, &res), std::runtime_error); + std::filesystem::remove(path); +} + +TEST(Serialization, LoadResultXml_ThrowsOnMissingFile) +{ + EXPECT_THROW(load_result_xml("/tmp/does_not_exist_conflab.xml"), + std::runtime_error); +} + +TEST(Serialization, LoadResultXml_ThrowsOnMalformedSolverAttribute) +{ + // V2: non-numeric attribute causes stoi/stod error that must surface + // as runtime_error, not std::invalid_argument. + const std::string path = "/tmp/conflab_badxml.xml"; + { + std::ofstream ofs(path); + ofs << R"( + + + 0.0 +)"; + } + NewtonResult res; + EXPECT_THROW(load_result_xml(path, &res), std::runtime_error); + std::filesystem::remove(path); +} diff --git a/code/tests/cgal/test_mesh_io.cpp b/code/tests/cgal/test_mesh_io.cpp index f767211..2602339 100644 --- a/code/tests/cgal/test_mesh_io.cpp +++ b/code/tests/cgal/test_mesh_io.cpp @@ -18,6 +18,7 @@ #include #include #include +#include #include using namespace conformallab; @@ -111,6 +112,30 @@ TEST(MeshIO, LoadMeshThrowsOnMissingFile) EXPECT_THROW(load_mesh("/tmp/does_not_exist_conflab.off"), std::runtime_error); } +// ════════════════════════════════════════════════════════════════════════════ +// I3: load_mesh throws on non-triangulated mesh +// +// The quad-strip mesh has 4-vertex faces; load_mesh must reject it at the +// I/O boundary rather than letting the quad faces flow silently into the +// triangle-only functionals. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(MeshIO, LoadMeshThrowsOnNonTriangulatedMesh) +{ + // Write a minimal OFF quad-face mesh (2 quads, not triangles). + auto path = tmp_path("quad.off"); + { + std::ofstream ofs(path); + ofs << "OFF\n6 2 0\n" + << "0 0 0\n1 0 0\n1 1 0\n0 1 0\n" + << "2 0 0\n2 1 0\n" + << "4 0 1 2 3\n" + << "4 1 4 5 2\n"; + } + EXPECT_THROW(load_mesh(path), std::runtime_error); + rm(path); +} + // ════════════════════════════════════════════════════════════════════════════ // Vertex positions survive a round-trip (OFF) // ════════════════════════════════════════════════════════════════════════════ @@ -171,3 +196,42 @@ TEST(MeshIO, SaveLoadConvenienceWrappers) rm(path); } + +// ════════════════════════════════════════════════════════════════════════════ +// V3: load_mesh throws on non-finite vertex coordinate (NaN / Inf) +// +// A mesh file containing NaN or Inf vertex coordinates must be rejected at +// the I/O boundary — before the coordinates can silently poison the solver. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(MeshIO, LoadMeshThrowsOnNaNCoordinate) +{ + // Write a minimal OFF triangle with a NaN in the first vertex. + auto path = tmp_path("nan.off"); + { + std::ofstream ofs(path); + ofs << "OFF\n3 1 0\n" + << "nan 0.0 0.0\n" + << "1.0 0.0 0.0\n" + << "0.0 1.0 0.0\n" + << "3 0 1 2\n"; + } + EXPECT_THROW(load_mesh(path), std::runtime_error); + rm(path); +} + +TEST(MeshIO, LoadMeshThrowsOnInfCoordinate) +{ + // Write a minimal OFF triangle with Inf in the y-coordinate. + auto path = tmp_path("inf.off"); + { + std::ofstream ofs(path); + ofs << "OFF\n3 1 0\n" + << "0.0 0.0 0.0\n" + << "1.0 inf 0.0\n" + << "0.0 1.0 0.0\n" + << "3 0 1 2\n"; + } + EXPECT_THROW(load_mesh(path), std::runtime_error); + rm(path); +} diff --git a/code/tests/cgal/test_newton_phase9a.cpp b/code/tests/cgal/test_newton_phase9a.cpp index 03e5a1e..689e9d1 100644 --- a/code/tests/cgal/test_newton_phase9a.cpp +++ b/code/tests/cgal/test_newton_phase9a.cpp @@ -16,6 +16,7 @@ #include "newton_solver.hpp" #include "cp_euclidean_functional.hpp" #include "inversive_distance_functional.hpp" +#include "inversive_distance_hessian.hpp" #include "mesh_builder.hpp" #include "conformal_mesh.hpp" @@ -232,6 +233,49 @@ TEST(NewtonPhase9a, InversiveDistance_PerturbedTetrahedron_Converges) EXPECT_LT(res.grad_inf_norm, 1e-8); } +// ════════════════════════════════════════════════════════════════════════════ +// 6b. Inversive-Distance block-FD Hessian == full-FD Hessian (B1 port) +// +// The solver was switched from the O(n·F) full-FD Hessian to the O(F) +// per-face block-FD Hessian. The two must agree to FD rounding on any +// configuration — including a perturbed (non-equilibrium) one, where the +// off-diagonal coupling is non-trivial. This is the locality-lemma +// cross-validation that licenses the faster path. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(NewtonPhase9a, InversiveDistance_BlockFDHessianMatchesFullFD) +{ + auto mesh = make_tetrahedron(); + auto m = setup_inversive_distance_maps(mesh); + compute_inversive_distance_init_from_mesh(mesh, m); + + // Pin one vertex; index the rest. + auto vit = mesh.vertices().begin(); + m.v_idx[*vit++] = -1; + int n = 0; + for (; vit != mesh.vertices().end(); ++vit) m.v_idx[*vit] = n++; + + // Evaluate the two Hessians away from equilibrium so off-diagonals are live. + std::vector x(static_cast(n), 0.0); + for (int i = 0; i < n; ++i) x[static_cast(i)] = 0.07 * (i + 1); + + auto H_full = inversive_distance_hessian_sym(mesh, x, m); + auto H_block = inversive_distance_hessian_block_fd_sym(mesh, x, m); + + ASSERT_EQ(H_full.rows(), H_block.rows()); + ASSERT_EQ(H_full.cols(), H_block.cols()); + + Eigen::MatrixXd Df = Eigen::MatrixXd(H_full); + Eigen::MatrixXd Db = Eigen::MatrixXd(H_block); + double max_abs_diff = (Df - Db).cwiseAbs().maxCoeff(); + EXPECT_LT(max_abs_diff, 1e-7) + << "block-FD and full-FD Hessians must agree to FD rounding;\n" + << "max |Δ| = " << max_abs_diff; + + // Sanity: the matrices are non-trivial (not both accidentally zero). + EXPECT_GT(Df.cwiseAbs().maxCoeff(), 1e-3); +} + // ════════════════════════════════════════════════════════════════════════════ // 7. CP-Euclidean Newton — uses analytic Hessian (NOT FD) // diff --git a/code/tests/cgal/test_newton_solver.cpp b/code/tests/cgal/test_newton_solver.cpp index d9ce3c3..b60da71 100644 --- a/code/tests/cgal/test_newton_solver.cpp +++ b/code/tests/cgal/test_newton_solver.cpp @@ -500,4 +500,240 @@ TEST(SparseQRFallback, Euclidean_ClosedMeshNoPinConverges) << "Euclidean Newton on closed tetrahedron (no pin) must converge via SparseQR; " "grad_inf_norm = " << res.grad_inf_norm; EXPECT_LT(res.grad_inf_norm, 1e-8); + EXPECT_EQ(res.status, NewtonStatus::Converged); + // N7: the closed mesh with no pinned vertex has a gauge-singular Hessian. + // That near-singularity must be observable in the diagnostics — either the + // SparseQR fallback fired, or (when LDLT "succeeds" with a ~0 pivot, which is + // exactly the silent case N7 targets) the smallest LDLT pivot is tiny. + EXPECT_TRUE(res.sparse_qr_fallback_used || res.min_ldlt_pivot < 1e-6) + << "gauge singularity should surface via fallback flag or min_ldlt_pivot; " + << "fallback=" << res.sparse_qr_fallback_used + << " min_pivot=" << res.min_ldlt_pivot; +} + +// ════════════════════════════════════════════════════════════════════════════ +// C1: solve_linear_system exposes the ok flag via the new out-parameter +// +// The C1 audit finding was that solve_linear_system silently dropped the +// internal ok flag, making double-solver failure invisible to callers. +// These tests verify that the new optional bool* ok parameter is correctly +// wired: ok=true for successful solves, and the parameter is optional (null +// is safe). +// +// Note: constructing a matrix where Eigen's SparseQR hard-fails is not +// practical — SparseQR reports Eigen::Success even for rank-0 inputs, +// returning a zero min-norm solution. The observable contract for callers is +// therefore: ok=true whenever a solver reports success; ok=false only when +// both report failure (an edge case Eigen essentially never produces). The +// fix in C1 ensures that if that edge case ever fires it propagates to the +// caller — previously it was silently swallowed. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SparseQRFallback, OkFlag_TrueOnFullRankSystem) +{ + // 3×3 diagonal PD matrix: LDLT succeeds → ok must be true. + Eigen::SparseMatrix A(3, 3); + A.insert(0, 0) = 1.0; + A.insert(1, 1) = 2.0; + A.insert(2, 2) = 3.0; + A.makeCompressed(); + + Eigen::VectorXd rhs(3); + rhs << 1.0, 4.0, 9.0; + + bool fallback = true; + bool ok = false; + Eigen::VectorXd x = conformallab::solve_linear_system(A, rhs, &fallback, &ok); + + EXPECT_TRUE(ok) << "Full-rank system: ok must be true"; + EXPECT_FALSE(fallback) << "Full-rank system: LDLT path, no fallback expected"; + EXPECT_NEAR(x[0], 1.0, 1e-12); + EXPECT_NEAR(x[1], 2.0, 1e-12); + EXPECT_NEAR(x[2], 3.0, 1e-12); +} + +TEST(SparseQRFallback, OkFlag_TrueOnSingularSystemViaFallback) +{ + // Singular matrix (zero row/col 1) — LDLT fails, SparseQR succeeds → ok=true. + Eigen::SparseMatrix A(3, 3); + A.insert(0, 0) = 2.0; + A.insert(2, 2) = 3.0; + A.makeCompressed(); + + Eigen::VectorXd rhs(3); + rhs << 2.0, 0.0, 3.0; + + bool fallback = false; + bool ok = false; + Eigen::VectorXd x = conformallab::solve_linear_system(A, rhs, &fallback, &ok); + + EXPECT_TRUE(ok) << "Singular-but-consistent system: SparseQR succeeds → ok must be true"; + EXPECT_TRUE(fallback) << "Singular system: fallback must have been triggered"; +} + +TEST(SparseQRFallback, OkFlag_NullPointerIsSafe) +{ + // Calling with ok=nullptr must not crash (backward-compatibility). + Eigen::SparseMatrix A(2, 2); + A.insert(0, 0) = 1.0; + A.insert(1, 1) = 1.0; + A.makeCompressed(); + + Eigen::VectorXd rhs(2); + rhs << 3.0, 5.0; + + EXPECT_NO_THROW({ + Eigen::VectorXd x = conformallab::solve_linear_system(A, rhs, nullptr, nullptr); + EXPECT_NEAR(x[0], 3.0, 1e-12); + EXPECT_NEAR(x[1], 5.0, 1e-12); + }); +} + +// ════════════════════════════════════════════════════════════════════════════ +// H2 / B2 — newton_core reuses the line-search gradient +// +// All five solvers now share detail::newton_core. B2: the gradient computed at +// the accepted line-search point is carried into the next iteration's +// convergence check instead of being recomputed at the loop top. On a +// unit-Hessian quadratic, exact Newton reaches the minimum in one step, so the +// total gradient-eval count is exactly 2 (1 initial + 1 accepted trial); the +// pre-B2 loop would have recomputed G at the top of iteration 1 → 3. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(NewtonCore, ReusesLineSearchGradient_B2_Convex) +{ + const std::vector xstar = {1.0, -2.0, 3.5}; + int grad_calls = 0; + + auto grad = [&](const std::vector& x) { + ++grad_calls; + std::vector g(x.size()); + for (std::size_t i = 0; i < x.size(); ++i) g[i] = x[i] - xstar[i]; + return g; // G(x) = x − x*, H = +I + }; + auto hess = [&](const std::vector&) { + Eigen::SparseMatrix H(3, 3); + for (int i = 0; i < 3; ++i) H.insert(i, i) = 1.0; + H.makeCompressed(); + return H; + }; + + auto res = conformallab::detail::newton_core( + std::vector{0.0, 0.0, 0.0}, grad, hess, + /*concave=*/false, /*tol=*/1e-9, /*max_iter=*/50); + + ASSERT_TRUE(res.converged); + EXPECT_EQ(res.iterations, 1); + EXPECT_NEAR(res.x[0], 1.0, 1e-9); + EXPECT_NEAR(res.x[1], -2.0, 1e-9); + EXPECT_NEAR(res.x[2], 3.5, 1e-9); + EXPECT_EQ(grad_calls, 2) + << "B2: the loop-top gradient must be reused from the line search"; +} + +// Concave path (spherical-style): H is NSD, newton_core factors −H and solves +// (−H)·Δx = G. G(x) = x* − x with H = −I makes x* a maximiser; the core must +// still converge in one step. +TEST(NewtonCore, ConcavePathConverges) +{ + const std::vector xstar = {0.5, -1.5, 2.0}; + int grad_calls = 0; + + auto grad = [&](const std::vector& x) { + ++grad_calls; + std::vector g(x.size()); + for (std::size_t i = 0; i < x.size(); ++i) g[i] = xstar[i] - x[i]; + return g; // G(x) = x* − x, H = −I + }; + auto hess = [&](const std::vector&) { + Eigen::SparseMatrix H(3, 3); + for (int i = 0; i < 3; ++i) H.insert(i, i) = -1.0; + H.makeCompressed(); + return H; + }; + + auto res = conformallab::detail::newton_core( + std::vector{0.0, 0.0, 0.0}, grad, hess, + /*concave=*/true, /*tol=*/1e-9, /*max_iter=*/50); + + ASSERT_TRUE(res.converged); + EXPECT_EQ(res.iterations, 1); + EXPECT_NEAR(res.x[0], 0.5, 1e-9); + EXPECT_NEAR(res.x[1], -1.5, 1e-9); + EXPECT_NEAR(res.x[2], 2.0, 1e-9); + EXPECT_EQ(grad_calls, 2); +} + +// ════════════════════════════════════════════════════════════════════════════ +// I1 / H1 / N7 — newton_core termination status, iteration count, diagnostics +// ════════════════════════════════════════════════════════════════════════════ + +namespace { +// 1×1 sparse identity / scalar helpers for the synthetic newton_core tests. +inline Eigen::SparseMatrix diag_sparse(std::initializer_list d) +{ + const int n = static_cast(d.size()); + Eigen::SparseMatrix H(n, n); + int i = 0; + for (double v : d) { H.insert(i, i) = v; ++i; } + H.makeCompressed(); + return H; +} +} // namespace + +// I1: a converging solve reports status == Converged and the N7 pivot is the +// (well-conditioned) unit diagonal; no SparseQR fallback. +TEST(NewtonCore, Status_Converged_AndDiagnostics_N7) +{ + const std::vector xstar = {2.0, -1.0}; + auto grad = [&](const std::vector& x) { + return std::vector{ x[0] - xstar[0], x[1] - xstar[1] }; + }; + auto hess = [&](const std::vector&) { return diag_sparse({1.0, 1.0}); }; + + auto res = conformallab::detail::newton_core( + std::vector{0.0, 0.0}, grad, hess, /*concave=*/false, 1e-9, 50); + + EXPECT_TRUE(res.converged); + EXPECT_EQ(res.status, conformallab::NewtonStatus::Converged); + EXPECT_EQ(res.iterations, 1); + EXPECT_FALSE(res.sparse_qr_fallback_used); // N7 + EXPECT_NEAR(res.min_ldlt_pivot, 1.0, 1e-12); // N7: D = I +} + +// I1: a slow (linearly-convergent cubic) problem that does not reach tol within +// max_iter reports MaxIterations, with iterations == max_iter (H1). +TEST(NewtonCore, Status_MaxIterations) +{ + // G(x) = x³, H = 3x² → Newton map x ← (2/3)x (linear convergence). + auto grad = [&](const std::vector& x) { + return std::vector{ x[0]*x[0]*x[0] }; + }; + auto hess = [&](const std::vector& x) { + return diag_sparse({ 3.0 * x[0] * x[0] }); + }; + + auto res = conformallab::detail::newton_core( + std::vector{1.0}, grad, hess, /*concave=*/false, /*tol=*/1e-8, /*max_iter=*/3); + + EXPECT_FALSE(res.converged); + EXPECT_EQ(res.status, conformallab::NewtonStatus::MaxIterations); + EXPECT_EQ(res.iterations, 3); +} + +// I1/H1: a constant non-zero gradient cannot be reduced by any step, so the +// line search stalls on the first iteration → LineSearchStalled, iterations == 0. +TEST(NewtonCore, Status_LineSearchStalled) +{ + auto grad = [&](const std::vector&) { + return std::vector{ 1.0, 1.0 }; // constant, independent of x + }; + auto hess = [&](const std::vector&) { return diag_sparse({1.0, 1.0}); }; + + auto res = conformallab::detail::newton_core( + std::vector{0.0, 0.0}, grad, hess, /*concave=*/false, 1e-9, 50); + + EXPECT_FALSE(res.converged); + EXPECT_EQ(res.status, conformallab::NewtonStatus::LineSearchStalled); + EXPECT_EQ(res.iterations, 0); // H1: no step completed } diff --git a/code/tests/cgal/test_spherical_hessian.cpp b/code/tests/cgal/test_spherical_hessian.cpp index 495b2dd..b888bd8 100644 --- a/code/tests/cgal/test_spherical_hessian.cpp +++ b/code/tests/cgal/test_spherical_hessian.cpp @@ -206,3 +206,31 @@ TEST(SphericalHessian, FDCheck_MixedPinnedVertices) EXPECT_TRUE(hessian_check_spherical(mesh, x, maps)) << "FD Hessian check failed for mixed pinned/variable vertices"; } + +// ════════════════════════════════════════════════════════════════════════════ +// I4: spherical_hessian throws on edge DOFs +// +// The spherical Hessian only implements the vertex-block cotangent Laplacian; +// if any edge has a free DOF index (e_idx >= 0) it must fail loudly rather +// than returning a silently rank-deficient matrix. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SphericalHessian, ThrowsOnEdgeDOF) +{ + // Build a tetrahedron and assign one edge as a free DOF. + auto mesh = make_tetrahedron(); + auto maps = setup_spherical_maps(mesh); + compute_lambda0_from_mesh(mesh, maps); // spherical variant (A2 rename pending) + + int idx = 0; + for (auto v : mesh.vertices()) + maps.v_idx[v] = idx++; + const int n_v = idx; + + // Give one edge a free DOF index to trigger the guard. + auto e0 = *mesh.edges().begin(); + maps.e_idx[e0] = n_v; // first free edge DOF + + std::vector x(static_cast(n_v + 1), 0.0); + EXPECT_THROW(spherical_hessian(mesh, x, maps), std::logic_error); +} diff --git a/doc/math/geometry-modes.md b/doc/math/geometry-modes.md index baa67f6..bfb3b1f 100644 --- a/doc/math/geometry-modes.md +++ b/doc/math/geometry-modes.md @@ -80,8 +80,27 @@ lᵢⱼ(ζ) — edge length from ζ value βᵢ(l₁, l₂, l₃) — vertex angle sum contribution ``` -**Hessian:** currently a symmetric finite-difference approximation (see `hyper_ideal_hessian.hpp`). -The analytic Hessian through the `(bᵢ, aₑ) → lᵢⱼ → ζ → αᵢⱼ/βᵢ` chain is Phase 9b. +**Hessian:** a symmetric finite-difference approximation (see `hyper_ideal_hessian.hpp`). +The solver uses the **per-face block-FD** variant by default (`O(F)` face evaluations, +~33–1166× faster than full-FD); the full-FD baseline is retained for cross-validation. +The analytic Hessian through the `(bᵢ, aₑ) → lᵢⱼ → ζ → αᵢⱼ/βᵢ` chain is Phase 9b-analytic. + +**Scale-floor clamp mode (`HyperIdealScaleClamp`).** The vertex scale `bᵢ` must stay +positive for the geometry to be valid. When a Newton iterate drives `bᵢ` toward or below +zero, one of two floors is applied — selectable per solver call (last argument of +`newton_hyper_ideal`, default `HardJava`): + +| Mode | Behaviour | When to use | +|---|---|---| +| `HardJava` *(default)* | `bᵢ < 0 → HYPER_IDEAL_SCALE_FLOOR` (= 0.01). Bit-for-bit faithful to `HyperIdealFunctional.java`, so the Java golden-oracle parity tests hold **exactly**. It is only C⁰ (in fact value-discontinuous at `bᵢ = 0`), so a step crossing the feasibility boundary sees a kink that can stall Newton. | Default — preserves parity; correct whenever scales stay positive (the normal case). | +| `SmoothBarrier` | `bᵢ ↦ floor + softplusβ(bᵢ − floor)` with `β = HYPER_IDEAL_SCALE_SHARPNESS` (= 100). C¹ everywhere, ≈ `bᵢ` to machine precision once `β·(bᵢ−floor) ≳ 35`, and decays smoothly to `floor` as `bᵢ → −∞`. Removes the N3 kink, but perturbs values near the boundary and therefore **deviates from the Java oracle** there. | Opt in when a solve is expected to cross the feasibility boundary and convergence robustness matters more than strict Java parity. | + +Both modes share `HYPER_IDEAL_SCALE_FLOOR`; the edge clamp (`aₑ < 0 → 0`) is unchanged. +Away from the floor (the normal `bᵢ ≫ 0.01` regime) the two modes are numerically +identical — the `LawsonHyperIdeal.SmoothBarrierConvergesToGoldenVector` test confirms the +smooth mode reaches the same golden solution as `HardJava`. Implemented in +`clamp_hyper_ideal_scale` (`hyper_ideal_functional.hpp`); origin: numerical-stability +audit finding **N3**. **Holonomy:** Möbius isometries Tᵢ ∈ SU(1,1) (orientation-preserving isometries of the Poincaré disk). `MobiusMap` in `layout.hpp`: T(z) = (az+b)/(cz+d). diff --git a/doc/math/references.md b/doc/math/references.md index 8ca8e0b..ce90980 100644 --- a/doc/math/references.md +++ b/doc/math/references.md @@ -29,12 +29,14 @@ Java reference implementation: [github.com/varylab/conformallab](https://github. | Reference | Used in | |---|---| | ✅ **Springborn** — *Ideal Hyperbolic Polyhedra and Discrete Uniformization*, Discrete & Computational Geometry **64** (2020), pp. 63–108. DOI: [10.1007/s00454-019-00132-8](https://doi.org/10.1007/s00454-019-00132-8) | `hyper_ideal_geometry.hpp` — ζ₁₃/ζ₁₄/ζ₁₅ functions; `hyper_ideal_functional.hpp` | +| ✅ **Kolpakov, Mednykh** — *A Formula for the Volume of a Hyperbolic Tetrahedron*, arXiv: [math/0603097](https://arxiv.org/abs/math/0603097) (2006) | Tetrahedron volume with one ideal vertex: `calculateTetrahedronVolumeWithIdealVertexAtGamma` in `hyper_ideal_utility.hpp` (Phase 9b analytic Hessian) | +| ✅ **Meyerhoff, Ushijima** — *A Note on the Dirichlet Domain*, in: The Epstein Birthday Schrift (2006) | Tetrahedron volume with three ideal vertices: `calculateTetrahedronVolumeFullyIdeal` in `hyper_ideal_utility.hpp` | | **Pinkall, Polthier** — *Computing Discrete Minimal Surfaces and Their Conjugates*, Experimental Mathematics (1993) | `euclidean_hessian.hpp` — cotangent Laplacian | | **Bobenko, Springborn** — *Variational Principles for Circle Patterns and Koebe's Theorem*, Transactions AMS (2004) | Variational angle-sum framework underlying all three functionals | | **Luo** — *Combinatorial Yamabe Flow on Surfaces*, Communications in Contemporary Mathematics (2004) | Inversive-distance functional — **new research** in Phase 9a.2 (no Java original; implemented from this paper + Glickenstein 2011 + Bowers-Stephenson 2004) | | **Bowers, Stephenson** — *Uniformizing dessins and Belyĭ maps via circle packing*, Memoirs of the AMS 170(805) (2004) | Introduces **inversive-distance circle packings** (used in Phase 9a.2). *Hinweis:* die zur Initialisierung benutzte Formel I_ij = (ℓ²−r_i²−r_j²)/(2 r_i r_j) ist die **klassische** inversive Distanz (vgl. Glickenstein §5.2: ℓ²=r_i²+r_j²+2r_ir_jη), nicht eine eigene „Bowers-Stephenson-Identität" — B–S liefern die Packungstheorie, nicht diese Formel. | | **Glickenstein** — *Discrete conformal variations and scalar curvature on piecewise flat two- and three-dimensional manifolds*, J. Differential Geometry **87**(2) (2011), pp. 201–238 | Analytic Hessian of the inversive-distance functional. ⚠️ *Korrektur:* die Arbeit nummeriert Gleichungen **nicht** im Format „(4.6)" — der Verweis ist durch die **§5.2**-Parametrisierung ℓ²_ij = r²_i + r²_j + 2 r_i r_j η_ij zu ersetzen. Cross-correspondence: η_ij ist die inversive Distanz und entspricht dem Kosinus des **Supplements** des Schnittwinkels (Schnitt bei arccos(−η_ij)) — also I_ij = cos θ_e **nur bis aufs Vorzeichen/Supplement**, nicht wörtlich. | -| **Bobenko, Pinkall, Springborn** — *Discrete conformal maps and ideal hyperbolic polyhedra*, Geometry & Topology **19**(4) (2015), pp. 2155–2215. arXiv: [1005.2698](https://arxiv.org/abs/1005.2698) | Face-based circle-packing functional (`CPEuclideanFunctional.java` → `cp_euclidean_functional.hpp`, Phase 9a.1) | +| ✅ **Bobenko, Pinkall, Springborn** — *Discrete conformal maps and ideal hyperbolic polyhedra*, Geometry & Topology **19**(4) (2015), pp. 2155–2215. arXiv: [1005.2698](https://arxiv.org/abs/1005.2698) (first posted 2010) | Face-based circle-packing functional (`CPEuclideanFunctional.java` → `cp_euclidean_functional.hpp`, Phase 9a.1) | | **Schläfli** — *On the multiple integral ∫dx dy …*, Quarterly Journal of Pure and Applied Mathematics (1858/60) | Klassische Schläfli-Differentialformel (dV = −½ Σ_e ℓ_e dθ_e). ⚠️ *Hinweis:* die in Phase 9b-analytic benutzte **Randterm-Form** `2 dV = Σ_e aₑ dαₑ + Σ_v bᵥ dβᵥ` steht **nicht** bei Schläfli 1858, sondern ist die verallgemeinerte Fassung für Mannigfaltigkeiten mit Rand → korrekter Beleg: **Rivin–Schlenker 1999** (Phase-10-Liste). Schläfli 1858 nur als historischer Ursprung zitieren. | | **Erickson, Whittlesey** — *Greedy Optimal Homotopy and Homology Generators*, SODA (2005) | `cut_graph.hpp` — tree-cotree algorithm | | **Bobenko, Springborn** — *A Discrete Laplace–Beltrami Operator for Simplicial Surfaces*, Discrete & Computational Geometry (2007) | Background for cotangent weights | diff --git a/doc/reviewer/README.md b/doc/reviewer/README.md index 2a339f6..d867ca5 100644 --- a/doc/reviewer/README.md +++ b/doc/reviewer/README.md @@ -37,6 +37,21 @@ a concrete fix, and acceptance criteria. > submission. Referenced from three audits: `cgal-submission-readiness` (G0/G1), > `dependency-license` (D1/D2 — incl. the unlicensed test meshes), and indirectly > `math-derivation-citation` (provenance of the Sechelmann-2016 primary source). +> +> **🔄 Status 2026-05-31:** the original authors (Sechelmann / Springborn, +> TU Berlin / Varylab) have been **contacted by email** to clarify +> porting/relicensing rights. **Awaiting reply.** Until a written answer exists, +> all G0-gated work (G1–G12, D1/D2, A4/A5, public release, CGAL submission) +> stays paused. + +## 🗺️ Working plan across sessions & models + +See [`finding-orchestration.md`](finding-orchestration.md) — the meta-plan that +maps every finding to a session, an implementing model (Haiku/Sonnet/Opus), and +an Opus review gate. **Session 1 (2026-05-31) is complete:** 26 findings +resolved across Sonnet/Haiku/Opus, 298/298 CGAL tests green (branch +`fix/b1-v3-c1-quick-wins`). The orchestration doc lists the remaining sessions +S2–S6. ## Correctness & port fidelity @@ -45,33 +60,48 @@ a concrete fix, and acceptance criteria. | [`external-audit-2026-05-30.md`](external-audit-2026-05-30.md) | Correctness bugs | 3 findings (face_energy, Gauss–Bonnet, cotangent doc) — **all resolved** | | [`java-port-audit.md`](java-port-audit.md) | Java→C++ math fidelity | line-by-line port anomalies — mostly resolved | | [`java-ignore-crossvalidation.md`](java-ignore-crossvalidation.md) | `@Ignore` cross-validation | which Java tests were intentionally skipped and why | -| [`numerical-stability-audit-2026-05-31.md`](numerical-stability-audit-2026-05-31.md) | FP robustness | N1: FD-Hessian step (`1e-5`) fights Newton tol (`1e-8`); clamps, cancellation, magic constants | -| [`math-derivation-citation-audit-2026-05-31.md`](math-derivation-citation-audit-2026-05-31.md) | Derivations & citations | M1 Kolpakov–Mednykh missing from refs; M2 BPS cited "2010" vs "2015" | +| [`numerical-stability-audit-2026-05-31.md`](numerical-stability-audit-2026-05-31.md) | FP robustness | **N3/N4/N5/N6 ✅** (clamp mode, named constants, Kahan area), **N7 ✅** (conditioning diag, S2); **N1 ✅\*** subsumed by B1; **N2 open** | +| [`math-derivation-citation-audit-2026-05-31.md`](math-derivation-citation-audit-2026-05-31.md) | Derivations & citations | **M1/M2/M4 ✅** (refs + year + citation format); **M3 open**, **M5 needs expert** | ## API, performance & packaging | Document | Focus | Headline | |---|---|---| -| [`api-performance-audit-2026-05-31.md`](api-performance-audit-2026-05-31.md) | API naming + Newton perf | A1–A5 naming (**decided**, see below); B1: block-FD Hessian exists but solver uses slow full-FD | -| [`cgal-submission-readiness-audit-2026-05-31.md`](cgal-submission-readiness-audit-2026-05-31.md) | CGAL submission | G0 (rights), G1 (MIT≠GPL), layout, concept, user manual, test harness | -| [`dependency-license-audit-2026-05-31.md`](dependency-license-audit-2026-05-31.md) | Deps + data licensing | D1 matrix predicated on the G0-undermined MIT premise; D2 test meshes have no provenance | +| [`api-performance-audit-2026-05-31.md`](api-performance-audit-2026-05-31.md) | API naming + Newton perf | **A1–A3 ✅**, **B1 ✅** (both halves), **B2–B5 ✅** (newton_core refactor); **A4/A5 deferred** (G0/G1) | +| [`cgal-submission-readiness-audit-2026-05-31.md`](cgal-submission-readiness-audit-2026-05-31.md) | CGAL submission | **⛔ blocked by G0** (author contacted); G1 (MIT≠GPL), layout, concept, manual, test harness all pending | +| [`dependency-license-audit-2026-05-31.md`](dependency-license-audit-2026-05-31.md) | Deps + data licensing | **⛔ blocked by G0**; D1 matrix predicated on the G0-undermined MIT premise; D2 test meshes have no provenance | | [`usability-audit-2026-05-31.md`](usability-audit-2026-05-31.md) | Examples & docs | U1–U11 — **all resolved** (PR #34) | ## Robustness & operational | Document | Focus | Headline | |---|---|---| -| [`test-coverage-error-handling-audit-2026-05-31.md`](test-coverage-error-handling-audit-2026-05-31.md) | Coverage + error handling | C1 solver `ok` flag dropped; coverage gate **decided: 80 % line / 70 % branch / 90 % func** | -| [`input-validation-audit-2026-05-31.md`](input-validation-audit-2026-05-31.md) | Malformed input | V3: NaN/Inf vertex coords flow silently into the solver | -| [`thread-safety-audit-2026-05-31.md`](thread-safety-audit-2026-05-31.md) | Concurrency | no mutable static state (good); same-mesh concurrency contract undocumented | +| [`test-coverage-error-handling-audit-2026-05-31.md`](test-coverage-error-handling-audit-2026-05-31.md) | Coverage + error handling | **C1/C2/C3 ✅**, **I2/I3/I4 ✅**, **H2 ✅** (newton_core), **I1/H1 ✅** (status enum, S2); **I5/H3/H4/H5 open**; gate in ramp-up until I5 | +| [`input-validation-audit-2026-05-31.md`](input-validation-audit-2026-05-31.md) | Malformed input | **V1/V2/V3/V4 ✅** (NaN/Inf guard + JSON/XML error handling); **V5/V6 open** (polish) | +| [`thread-safety-audit-2026-05-31.md`](thread-safety-audit-2026-05-31.md) | Concurrency | no mutable static state (good); same-mesh concurrency contract **still undocumented (open, S4)** | -## Decisions taken (2026-05-31) +## Session 1 — implemented (2026-05-31) -- **Coverage gate:** 80 % line / 70 % branch / 90 % function (after fixing the - coverage script so it measures honestly — `test-coverage` C2/I5). +Branch `fix/b1-v3-c1-quick-wins`, 9 commits, **298/298 CGAL tests green**. +See [`finding-orchestration.md`](finding-orchestration.md) for the full plan. + +- **Sonnet:** B1(HyperIdeal)/V3/C1; C2/C3 (coverage gate, ramp-up); V1/V2/V4; I2/I3/I4. +- **Haiku:** N4/N6 (named constants); M1/M2/M4 (citations); A1–A3 (naming + `[[deprecated]]` aliases). +- **Opus:** B1 inversive-distance block-FD port; N5 (Kahan area); N3 (selectable + clamp mode `HardJava`|`SmoothBarrier`); H2 + B2/B3/B4/B5 (newton_core refactor). +- **🔍 Opus review:** all Sonnet/Haiku commits validated (value-identity, error + handling, coverage-gate logic, alias forwarding) — all correct. + +**Session 2 (Opus, same day):** I1 (NewtonStatus enum) + H1 (iteration count) + +N7 (`sparse_qr_fallback_used` / `min_ldlt_pivot` diagnostics) — done directly in +the warm post-refactor session, commit `90e966a`, 301/301 tests green. + +### Decisions standing +- **Coverage gate:** 80 % line / 70 % branch / 90 % function. Script now measures + honestly (C2 ✅); gate is in **ramp-up** (`SKIP_COVERAGE_GATE=1`) until **I5** + wires it to the full CGAL suite. - **API naming A1–A3 (internal):** standardized on `__` with - `[[deprecated]]` aliases — **implemented**, open as **PR #36** - (`refactor/api-naming-a1-a3`), 277/277 tests green. + `[[deprecated]]` aliases — **implemented**. - **API naming A4–A5 (public CGAL surface):** decided (`discrete_inversive_distance_euclidean`; both packers return `Circle_packing_result`) but **deferred** until the G0/G1 outcome — no point stabilizing a public API on a package whose release status is open. @@ -79,9 +109,10 @@ a concrete fix, and acceptance criteria. ## How to act on a finding in a new session Open the relevant audit, pick a finding, follow its "Fix" + "Acceptance criteria". -Build/verify with the commands in that document's header. Suggested first batch -(low-risk, high-value quick wins): **B1** (HyperIdeal block-FD one-liner), **B3**, -**C1**, **C2**, **V3**. +Build/verify with the commands in that document's header. **The Session-1 quick-win +batch (B1, B3, C1, C2, V3, …) is done** — for what to pick up next, follow the +session plan in [`finding-orchestration.md`](finding-orchestration.md): **S1+S2 +done**; the next pending batch is **S3** (H3/H4/H5 + V5/V6, Sonnet→Opus review). ## When to send to the reviewer diff --git a/doc/reviewer/api-performance-audit-2026-05-31.md b/doc/reviewer/api-performance-audit-2026-05-31.md index b788d37..c641691 100644 --- a/doc/reviewer/api-performance-audit-2026-05-31.md +++ b/doc/reviewer/api-performance-audit-2026-05-31.md @@ -14,6 +14,12 @@ act on it without prior context. Each finding includes: Status legend: 🔴 Critical · 🟡 Important · 🔵 Hint / nice-to-have +> **✅ Resolution status (2026-05-31, Session 1):** **A1–A3 ✅** (renamed with +> `[[deprecated]]` aliases), **B1 ✅** both halves (HyperIdeal + Inversive-Distance +> block-FD), **B2/B3/B4/B5 ✅** (folded into the `newton_core` refactor — see +> test-coverage **H2**). **A4/A5 deferred** until the G0/G1 licensing outcome +> (public CGAL surface). 298/298 tests green. Plan: [`finding-orchestration.md`](finding-orchestration.md). + > **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 diff --git a/doc/reviewer/cgal-submission-readiness-audit-2026-05-31.md b/doc/reviewer/cgal-submission-readiness-audit-2026-05-31.md index f011aa7..0b79ae9 100644 --- a/doc/reviewer/cgal-submission-readiness-audit-2026-05-31.md +++ b/doc/reviewer/cgal-submission-readiness-audit-2026-05-31.md @@ -12,6 +12,13 @@ act on it without prior context. Status legend: ⛔ Blocker (submission rejected on sight) · 🔴 Critical · 🟡 Important · 🔵 Polish +> **🔄 Status (2026-05-31):** the whole package — **G0** (porting rights), and +> everything it gates (**G1–G12**) — is **blocked pending the original authors' +> reply**. The authors have been **contacted by email** asking to clarify +> porting/relicensing rights. No layout/license/header work proceeds until a +> written answer exists. Scheduled as **S6** (Opus, after G0/G1) in +> [`finding-orchestration.md`](finding-orchestration.md). + > **Companion documents:** > - `external-audit-2026-05-30.md` — port-faithfulness bugs > - `test-coverage-error-handling-audit-2026-05-31.md` — test gaps + error handling diff --git a/doc/reviewer/dependency-license-audit-2026-05-31.md b/doc/reviewer/dependency-license-audit-2026-05-31.md index 42a2458..6343822 100644 --- a/doc/reviewer/dependency-license-audit-2026-05-31.md +++ b/doc/reviewer/dependency-license-audit-2026-05-31.md @@ -9,6 +9,12 @@ unresolved provenance/license blocker (CGAL audit **G0/G1**)? Status legend: 🔴 Critical · 🟡 Important · 🔵 Polish +> **🔄 Status (2026-05-31):** **D1/D2 blocked by G0** — the dependency/licensing +> matrix rests on the MIT premise that G0 (porting rights) undermines. The +> original authors have been **contacted by email**; no licensing action until +> their reply. Scheduled **S6** (after G0/G1) in +> [`finding-orchestration.md`](finding-orchestration.md). + > **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 diff --git a/doc/reviewer/finding-orchestration.md b/doc/reviewer/finding-orchestration.md new file mode 100644 index 0000000..be1de3f --- /dev/null +++ b/doc/reviewer/finding-orchestration.md @@ -0,0 +1,176 @@ +# Finding Orchestration — Sessions × Models × Review Gates + +**Updated:** 2026-05-31 +**Purpose:** a structured plan to work through every audit finding across +**multiple short sessions**, each run by the **model best suited** to the work, +with an **external Opus reviewer session after every implementation session**. +This lets you pick up any "Session" below cold, hand it to the named model, and +know exactly what it covers and how its output gets validated. + +> Companion: the per-finding detail lives in the 11 audit documents in this +> folder (see [`README.md`](README.md)). This file is the *meta-plan* — the +> order, the model assignment, and the review cadence. + +--- + +> **Ready-to-paste prompts** for every pending session (S3–S6 + the review gate) +> live in [`session-prompts.md`](session-prompts.md) — copy one block, set the +> named model, go. + +## How to use this + +1. Pick the next **⬜ pending** session from the [sequence](#session-sequence). +2. Start a session with the **assigned model**, point it at the listed findings + in their audit docs (each finding has `file:line`, fix, acceptance criteria). +3. When it finishes, start the paired **🔍 Opus review session** — it validates + the diff (correctness, value-identity, tests, no parity regression) and only + then is the batch "done". +4. Mark the session ✅ here and move on. + +Rationale for the cadence: implementation models (Sonnet/Haiku) are cheaper and +fast for well-specified work; Opus is reserved for (a) findings that need real +numerical/architectural judgement and (b) the **review gate**, because an +independent high-capability pass catches subtle parity/numerics regressions that +a self-review from the implementing model tends to miss. + +--- + +## Model-assignment heuristic + +| Model | Take findings that are… | Examples | +|---|---|---| +| **Haiku** | mechanical, local, low-risk; docs, citations, renames, named-constant extraction | M-series, N4/N6, A1–A3, doc files | +| **Sonnet** | normal coding with a clear acceptance criterion; tests, error-handling, CI, small API changes | V-series, C-series, I-series, H-series | +| **Opus** | numerics, math correctness, architecture, public-API/irreversible decisions — **and every review gate** | B1-port, N1/N3/N5/N7, H2 refactor, A4/A5, all 🔍 reviews | + +> Decision rule isn't "how severe" but "how much must be *understood* to get it +> right." A 🔴 one-line switch can be Sonnet; a 🟡 numerical reformulation is Opus. + +--- + +## Master finding table + +Status: ✅ done · ⬜ open (actionable) · ⏸ deferred (intentional) · ⛔ blocked (G0) · 👤 needs human expert + +| ID | Audit | Sev | Status | Model | Session | +|----|-------|-----|--------|-------|---------| +| B1 (HyperIdeal) | api-perf | 🔴 | ✅ | Sonnet | S1 | +| B1 (Inv-Dist port) | api-perf | 🔴 | ✅ | Opus | S1 | +| V3 | input-val | 🟡 | ✅ | Sonnet | S1 | +| C1 | test-cov | 🔴 | ✅ | Sonnet | S1 | +| N4, N6 | numerics | 🟡/🔵 | ✅ | Haiku | S1 | +| M1, M2, M4 | math-cite | 🟡/🔵 | ✅ | Haiku | S1 | +| C2, C3 | test-cov | 🔴 | ✅ | Sonnet | S1 | +| V1, V2, V4 | input-val | 🟡 | ✅ | Sonnet | S1 | +| I2, I3, I4 | test-cov | 🟡 | ✅ | Sonnet | S1 | +| A1, A2, A3 | api-perf | 🔴/🟡 | ✅ | Haiku | S1 | +| N5 | numerics | 🟡 | ✅ | Opus | S1 | +| N3 | numerics | 🟡 | ✅ | Opus | S1 | +| H2, B2, B3, B4, B5 | test-cov/api-perf | 🔵/🟡 | ✅ | Opus | S1 | +| I1 | test-cov | 🟡 | ✅ | Opus | S2 | +| H1 | test-cov | 🔵 | ✅ | Opus | S2 | +| N7 | numerics | 🔵 | ✅ | Opus | S2 | +| **H3** | test-cov | 🔵 | ⬜ | Sonnet→🔍Opus | **S3** | +| **H4** | test-cov | 🔵 | ⬜ | Sonnet→🔍Opus | **S3** | +| **H5** | test-cov | 🔵 | ⬜ | Sonnet→🔍Opus | **S3** | +| **V5, V6** | input-val | 🔵 | ⬜ | Sonnet→🔍Opus | **S3** | +| **N2** | numerics | 🟡 | ⬜ | Haiku→🔍Opus | **S4** | +| **thread-safety doc** | thread-safety | 🟡 | ⬜ | Haiku→🔍Opus | **S4** | +| **M3** | math-cite | 🟡 | ⬜ | Haiku→🔍Opus | **S4** | +| **I5** | test-cov | 🟡 | ⬜ | Sonnet→🔍Opus | **S5** | +| A4, A5 | api-perf | 🟡 | ⏸ | Opus | S6 (after G0/G1) | +| N1 | numerics | 🔴 | ✅* | — | (subsumed by B1; doc note folds into N2/S4) | +| G1–G12 | cgal | ⛔/🔴 | ⛔ | Opus | (after G0) | +| D1, D2 | dep-license | 🔴 | ⛔ | Opus | (after G0) | +| M5 | math-cite | 🔵 | 👤 | — | human domain expert | +| G0 | cgal | ⛔ | 🔄 in progress | owner | **author contacted by email — awaiting reply** | + +\* N1 (FD-step vs Newton tol) is largely subsumed by the B1 block-FD work; the +remaining piece is the tolerance-coupling note, folded into N2 (Session S4). + +--- + +## Session sequence + +### ✅ S1 — Audit quick-wins + numerics + refactor (DONE, 2026-05-31) +Branch `fix/b1-v3-c1-quick-wins`, 9 commits, 298/298 CGAL tests green. +- **Sonnet:** B1(HyperIdeal)/V3/C1; C2/C3, V1/V2/V4, I2/I3/I4. +- **Haiku:** N4/N6; M1/M2/M4; A1–A3. +- **Opus:** B1 inv-dist port; N5; N3; H2+B2/B3/B4/B5. +- **🔍 Opus review:** validated all Sonnet/Haiku commits (value-identity, error + handling, coverage gate, alias forwarding) — all correct. + +### ✅ S2 — Solver result diagnostics (DONE, 2026-05-31, Opus impl + self-review) +Done directly in the warm post-refactor Opus session (all three findings live in +the just-written `newton_core`/`NewtonResult`/`NewtonLinearSolver`). Commit +`90e966a`; 301/301 CGAL tests green. **I1** status enum + **H1** iteration fix + +**N7** `sparse_qr_fallback_used` / `min_ldlt_pivot` diagnostics. Follow-up **done** +(commit `957f506`): `NewtonStatus` + both diagnostics propagated into the public +CGAL result types (`Conformal_map_result`, `Hyper_ideal_map_result`, +`Circle_packing_result`) via the `CGAL::Newton_status` alias. + +

original S2 plan + +#### S2 — Solver result diagnostics (Sonnet → 🔍 Opus) +- **I1** — give `NewtonResult` a status enum + (`Converged`/`MaxIterations`/`LinearSolverFailed`/`LineSearchStalled`). The + `newton_core` refactor centralised every exit point, so this lands in one + place now. Update the 5 wrappers + any callers reading `converged`. +- **H1** — `res.iterations` is stale when the solver breaks on `!ok`/stall; + set it correctly at each `break` in `newton_core`. +- **N7** *(Opus within this session)* — optional conditioning diagnostic + (smallest LDLT pivot / cheap cond estimate) surfaced in the result; pairs + naturally with the I1 status enum (near-singular ⇒ a distinct status/flag). +- **🔍 Opus review:** confirm the enum doesn't break the CGAL-layer result + mapping and that `converged` semantics are unchanged for existing callers. + +
+ +### ⬜ S3 — Robustness & test-gap closure (Sonnet → 🔍 Opus) +- **H3** — `enforce_gauss_bonnet` reports the magnitude of the applied correction. +- **H4** — test `reduce_to_fundamental_domain` boundary `Im(τ)==0.0`. +- **H5** — integration test for the degenerate-triangle path in the Newton solve. +- **V5** — make the hand-rolled XML reader *reject* reformatted-but-valid XML + instead of silently mis-reading it into zeros. +- **V6** — DOF-vector-vs-mesh size check on load. +- **🔍 Opus review:** verify the new rejections don't break valid round-trips. + +### ⬜ S4 — Documentation & citations (Haiku → 🔍 Opus) +- **N2** — `doc/math/tolerances.md`: every numerical threshold, its role, and + the constraints between them (incl. the N1 FD-step ≪ Newton-tol coupling). +- **thread-safety** — document the same-mesh concurrency contract (the code has + no mutable static state; this is purely the written contract). +- **M3** — re-verify the future-dated / arXiv-only citations against their + published versions; update `references.md`. +- **🔍 Opus review:** sanity-check the tolerance relationships are stated correctly. + +### ⬜ S5 — Coverage truthfulness (Sonnet → 🔍 Opus) +- **I5** — wire `coverage.sh` to the full CGAL suite (not just the 26 fast + tests, ~9.6%), then remove `SKIP_COVERAGE_GATE=1` so the agreed + 80%/70%/90% gate goes live (C3 ramp-up → enforced). +- **🔍 Opus review:** confirm the measured numbers are honest and the gate fails + correctly when coverage drops. + +### ⏸ S6 — Public CGAL surface + packaging (Opus, **only after G0/G1**) +Gated by the author's reply on porting/relicensing rights. +- **A4, A5** — public entry-point renames + unified `Circle_packing_result`. +- **G1–G12, D1, D2** — license headers, CGAL package layout, concept doc, user + manual, CGAL test harness, dependency/data provenance. +- **G0 status:** *author contacted by email (2026-05-31), awaiting reply.* No + public release / submission / relicensing until a written answer exists. + +### 👤 Out of model scope +- **M5** — the large hand-derivations need a domain-expert prose review + (the numerical results are already validated; the prose is not). + +--- + +## Review-gate checklist (what every 🔍 Opus session checks) + +- [ ] Builds clean; full CGAL suite green (no count regression). +- [ ] No Java golden-vector / parity test perturbed (HardJava defaults intact). +- [ ] Numeric changes are value-identical where claimed, or justified + tested. +- [ ] New public surface (result types, enums, API) is intentional and documented. +- [ ] Commit message attributes the implementing model + (`Co-Authored-By: Claude `). +- [ ] Finding marked ✅ in the master table above with the commit ref. diff --git a/doc/reviewer/input-validation-audit-2026-05-31.md b/doc/reviewer/input-validation-audit-2026-05-31.md index e3d63a7..da378f1 100644 --- a/doc/reviewer/input-validation-audit-2026-05-31.md +++ b/doc/reviewer/input-validation-audit-2026-05-31.md @@ -10,6 +10,12 @@ which covered *internal* error propagation and the deliberate `throw` paths. Status legend: 🔴 Critical · 🟡 Important · 🔵 Polish +> **✅ Resolution status (2026-05-31, Session 1):** **V1/V2/V4 ✅** (JSON parse + +> field checks, XML stoi/stod guarded — all surface as `std::runtime_error` with +> path), **V3 ✅** (NaN/Inf vertex-coordinate guard in `load_mesh`). **V5** (XML +> strict-reject) and **V6** (DOF-size check) **open** — scheduled S3 in +> [`finding-orchestration.md`](finding-orchestration.md). 298/298 tests green. + > **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 diff --git a/doc/reviewer/math-derivation-citation-audit-2026-05-31.md b/doc/reviewer/math-derivation-citation-audit-2026-05-31.md index 2db1ba5..9368fae 100644 --- a/doc/reviewer/math-derivation-citation-audit-2026-05-31.md +++ b/doc/reviewer/math-derivation-citation-audit-2026-05-31.md @@ -10,6 +10,13 @@ checks code/docs vs. *the published mathematics*. Status legend: 🔴 Critical · 🟡 Important · 🔵 Polish +> **✅ Resolution status (2026-05-31, Session 1):** **M1 ✅** (Kolpakov–Mednykh + +> Meyerhoff/Ushijima added to `references.md`), **M2 ✅** (BPS year unified to the +> published 2015), **M4 ✅** (citation format normalized). **M3** (re-verify +> future-dated/arXiv citations) **open** → S4; **M5** (prose derivation review) +> needs a **domain expert** (out of model scope). See +> [`finding-orchestration.md`](finding-orchestration.md). + > **Good news up front:** `doc/math/references.md` is unusually scholarly and already > contains several *self-corrections* (the Glickenstein "eq. (4.6)" numbering fix → > §5.2 parametrization; the Schläfli boundary-term re-attribution to Rivin–Schlenker diff --git a/doc/reviewer/numerical-stability-audit-2026-05-31.md b/doc/reviewer/numerical-stability-audit-2026-05-31.md index 71d957c..f885bdf 100644 --- a/doc/reviewer/numerical-stability-audit-2026-05-31.md +++ b/doc/reviewer/numerical-stability-audit-2026-05-31.md @@ -10,6 +10,14 @@ port-faithfulness (`java-port-audit.md`) and from internal error propagation Status legend: 🔴 Critical · 🟡 Important · 🔵 Polish +> **✅ Resolution status (2026-05-31, Session 1):** **N3 ✅** (selectable clamp mode +> `HardJava`|`SmoothBarrier`), **N4/N6 ✅** (named `constexpr` constants), **N5 ✅** +> (Kahan stable triangle area), **N7 ✅** (S2 — `min_ldlt_pivot` / +> `sparse_qr_fallback_used` in `NewtonResult`). **N1 ✅\*** largely subsumed by the +> B1 block-FD Hessian; the remaining FD-step/tol coupling note folds into **N2**. +> **N2** (tolerances.md) **open** — scheduled S4 in +> [`finding-orchestration.md`](finding-orchestration.md). 301/301 tests green. + > **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 diff --git a/doc/reviewer/session-prompts.md b/doc/reviewer/session-prompts.md new file mode 100644 index 0000000..e8cedf9 --- /dev/null +++ b/doc/reviewer/session-prompts.md @@ -0,0 +1,146 @@ +# Ready-to-paste session prompts (S3–S6 + review gate) + +Copy one block into a fresh session, set the **model named in the prompt**, and go. +Each prompt is self-contained: it names the findings, the audit doc to follow, the +build/test commands, and the branch/push/PR + review-gate workflow. + +Shared conventions (already baked into each prompt): +- Repo: `/Users/tarikmoussa/Desktop/ConformalLabpp`, base branch `main`. +- Branch + push to the **eulernest fork** = remote `origin`; open the PR via the + Gitea API (`https://git.eulernest.eu/api/v1/repos/conformallab/ConformalLabpp/pulls`, + basic-auth with the token embedded in the `origin` URL), base `main`. +- Build/test: `cmake -S code -B build-cgal -DWITH_CGAL_TESTS=ON && cmake --build + build-cgal --target conformallab_cgal_tests -j8 && ctest --test-dir build-cgal -R '^cgal\.'` +- Plan + finding details: `doc/reviewer/finding-orchestration.md` and the per-finding + audit docs in `doc/reviewer/`. +- After each implementation session, run the **Review gate** prompt (Opus). + +--- + +## S3 — Robustness & test-gap closure (model: **Sonnet**) + +``` +Use the Sonnet model. Work in /Users/tarikmoussa/Desktop/ConformalLabpp on a new +branch off main called `fix/s3-robustness-gaps`. + +Implement audit findings H3, H4, H5, V5, V6 — details (file:line, fix, acceptance +criteria) are in doc/reviewer/test-coverage-error-handling-audit-2026-05-31.md +(H3/H4/H5) and doc/reviewer/input-validation-audit-2026-05-31.md (V5/V6): +- H3: enforce_gauss_bonnet reports the magnitude of the applied correction. +- H4: test reduce_to_fundamental_domain at the boundary Im(τ)==0.0. +- H5: integration test for the degenerate-triangle path in the Newton solver. +- V5: make the hand-rolled XML reader REJECT reformatted-but-valid XML instead of + silently mis-reading it into zeros (document the strict internal-only subset). +- V6: DOF-vector-vs-mesh size check on load, with a clear error. + +Add a negative/boundary test for each. Build and run the full CGAL suite +(cmake -S code -B build-cgal -DWITH_CGAL_TESTS=ON && cmake --build build-cgal +--target conformallab_cgal_tests -j8 && ctest --test-dir build-cgal -R '^cgal\.') — +it must stay green with your new tests added. Do NOT change the default behaviour +of any existing valid round-trip. + +Commit per finding (or in two logical commits) with trailer +`Co-Authored-By: Claude Sonnet 4.6 `. Push to origin and +open a PR (base main) via the Gitea API using the token in the origin remote URL. +Then update the status in doc/reviewer/finding-orchestration.md (H3/H4/H5/V5/V6 → ✅, S3) +and the audit banners. Report the PR URL and test count. +``` + +--- + +## S4 — Documentation & citations (model: **Haiku**) + +``` +Use the Haiku model. Work in /Users/tarikmoussa/Desktop/ConformalLabpp on a new +branch off main called `docs/s4-tolerances-threadsafety-citations`. + +Three documentation findings — no code logic changes: +- N2 (numerical-stability audit): create doc/math/tolerances.md listing EVERY + numerical threshold (Newton tol, FD hess_eps, Armijo c1, sparsity drop, + Gauss-Bonnet tol, spherical asin guard, the named constants in constants.hpp), + each with its role and the constraints between them — include the N1 coupling + "FD floor must be ≪ Newton tol". Cross-link from constants.hpp. +- thread-safety audit: document the same-mesh concurrency contract (the code has + NO mutable static state; this is purely the written contract — state what is and + isn't safe to call concurrently on the same vs different meshes). +- M3 (math-citation audit): re-verify the future-dated / arXiv-only citations in + doc/math/references.md against their published versions; correct any that have + since appeared in a journal/proceedings. + +Build still succeeds; run the full CGAL suite to confirm no accidental code change +(should stay green). Commit with trailer +`Co-Authored-By: Claude Haiku 4.5 `. Push to origin, open a +PR (base main) via the Gitea API. Update doc/reviewer/finding-orchestration.md +(N2/thread-safety/M3 → ✅, S4) and the audit banners. Report the PR URL. +``` + +--- + +## S5 — Coverage truthfulness (model: **Sonnet**) + +``` +Use the Sonnet model. Work in /Users/tarikmoussa/Desktop/ConformalLabpp on a new +branch off main called `ci/s5-coverage-full-suite`. + +Implement audit finding I5 (test-coverage-error-handling audit): the coverage +script currently measures only the 26 fast tests (~9.6%), not the CGAL suite. +- Wire scripts/quality/coverage.sh to build + run the full CGAL test suite + (WITH_CGAL_TESTS=ON) under coverage instrumentation, so code/include/ is measured + honestly. +- Once the measured numbers are real, remove the SKIP_COVERAGE_GATE=1 ramp-up flag + from .gitea/workflows/cpp-tests.yml so the agreed 80% line / 70% branch / 90% + function gate goes live. If real coverage is below threshold, either raise it + with targeted tests or document the gap and keep ramp-up — but state which. +- Verify locally: bash scripts/quality/coverage.sh should report the full-suite + numbers and pass/fail the gate correctly. + +Commit with trailer `Co-Authored-By: Claude Sonnet 4.6 `. +Push to origin, open a PR (base main) via the Gitea API. Update +doc/reviewer/finding-orchestration.md (I5 → ✅, S5; flip the gate note in the +test-coverage banner). Report the PR URL and the measured coverage numbers. +``` + +--- + +## S6 — Public CGAL surface + packaging (model: **Opus**) — ⛔ BLOCKED until G0 + +``` +PRECONDITION: do NOT start this session until the original Varylab/TU-Berlin +authors have replied in writing granting porting/relicensing rights (finding G0). +If there is no written reply yet, stop and report that S6 is still blocked. + +Use the Opus model. Work in /Users/tarikmoussa/Desktop/ConformalLabpp. +Once G0 is cleared, implement, in dependency order, from the CGAL submission audit +(doc/reviewer/cgal-submission-readiness-audit-2026-05-31.md) and api-performance +audit (A4/A5): +- A4/A5: public entry-point renames (discrete_inversive_distance_euclidean; both + packers return Circle_packing_result), with [[deprecated]] aliases + CHANGELOG. +- G1 (license headers per the chosen license), G2/G3/G8 (CGAL package layout, + drop vendored CGAL on the submission branch), G9 (public/internal header split), + G4/G5 (concept doc + user manual + PackageDescription), G6 (CGAL test harness), + G7 (deliver or narrow genericity), G10/G12 (benchmark + headers). + +This is a large, mostly-mechanical effort — split into several PRs. Keep the +standalone MIT build working in parallel. Report progress against the G-series list. +``` + +--- + +## Review gate (run after S3 / S4 / S5) (model: **Opus**) + +``` +Use the Opus model. Work in /Users/tarikmoussa/Desktop/ConformalLabpp. Review the +open PR as an independent reviewer. Check, and report a +pass/fail per item: +- Builds clean; full CGAL suite green with no test-count regression + (ctest --test-dir build-cgal -R '^cgal\.'). +- No Java golden-vector / parity test perturbed; HardJava clamp default intact. +- Numeric changes are value-identical where claimed, or justified + covered by a test. +- Any new public surface (result types, enums, API, file formats) is intentional + and documented; existing valid round-trips still work. +- Commit messages attribute the implementing model. +- The finding is marked ✅ in doc/reviewer/finding-orchestration.md with the commit ref. +Read the actual diff (git diff main...HEAD) and the touched audit doc. If you find +a real problem, fix it directly (small) or list precise required changes (larger), +then re-run the suite. Conclude with an explicit APPROVE / CHANGES-REQUESTED. +``` diff --git a/doc/reviewer/test-coverage-error-handling-audit-2026-05-31.md b/doc/reviewer/test-coverage-error-handling-audit-2026-05-31.md index 05743ea..eb6ab7e 100644 --- a/doc/reviewer/test-coverage-error-handling-audit-2026-05-31.md +++ b/doc/reviewer/test-coverage-error-handling-audit-2026-05-31.md @@ -14,6 +14,15 @@ act on it without prior context. Each finding includes: Status legend: 🔴 Critical · 🟡 Important · 🔵 Hint / nice-to-have +> **✅ Resolution status (2026-05-31, Session 1):** **C1 ✅** (solver `ok` flag +> exposed), **C2/C3 ✅** (coverage script honest + CI gate in ramp-up), +> **I2/I3/I4 ✅** (serialization / mesh / spherical-hessian negative tests), +> **H2 ✅** (five Newton loops unified into `newton_core`, folding in B2/B3/B4/B5). +> **I1 ✅** (`NewtonStatus` enum) + **H1 ✅** (iteration count) done in S2. +> **Open:** **H3/H4/H5** → S3; **I5** (coverage on full suite, un-gates the +> 80/70/90 threshold) → S5. +> See [`finding-orchestration.md`](finding-orchestration.md). 298/298 tests green. + > **Note:** This audit is complementary to `external-audit-2026-05-30.md` (which > covers port-faithfulness bugs). There is no overlap in findings — this one is > strictly about *test coverage* and *error-handling robustness*. diff --git a/doc/reviewer/thread-safety-audit-2026-05-31.md b/doc/reviewer/thread-safety-audit-2026-05-31.md index 1d14a92..588b9f8 100644 --- a/doc/reviewer/thread-safety-audit-2026-05-31.md +++ b/doc/reviewer/thread-safety-audit-2026-05-31.md @@ -9,6 +9,11 @@ documented? Currently undocumented everywhere. Status legend: 🔴 Critical · 🟡 Important · 🔵 Polish +> **Status (2026-05-31):** no code change needed (no mutable static state — good). +> The one open item is **documenting the same-mesh concurrency contract** — +> scheduled **S4** (Haiku→Opus review) in +> [`finding-orchestration.md`](finding-orchestration.md). + > **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`, diff --git a/scripts/quality/coverage.sh b/scripts/quality/coverage.sh index 8b9e050..f512f98 100755 --- a/scripts/quality/coverage.sh +++ b/scripts/quality/coverage.sh @@ -9,9 +9,9 @@ # * build-coverage/lcov-html/index.html — browseable HTML report # * stdout: per-file summary + grand total # -# Local-only. Not gated in CI yet; once a coverage threshold is agreed -# with the reviewer (e.g. 80 %), the gate can be a single line in -# cpp-tests.yml. +# Local + CI. Threshold gate: 80 % line / 70 % branch / 90 % function. +# Set SKIP_COVERAGE_GATE=1 to measure without failing (ramp-up mode). +# CI: invoked from cpp-tests.yml test-cgal job (after the full suite runs). # # Usage: # bash scripts/quality/coverage.sh # gcc default @@ -20,9 +20,11 @@ # Prerequisites: gcov + lcov (apt install lcov / brew install lcov) # # Exit codes: -# 0 coverage report generated; prints % +# 0 coverage report generated and all thresholds met # 1 tests failed (no usable trace) -# 2 prerequisite missing +# 2 prerequisite missing (lcov / compiler not found) +# 3 coverage.info is empty after lcov run (version mismatch) +# 4 coverage below threshold set -euo pipefail ROOT="$(cd "$(dirname "${BASH_SOURCE[0]}")/../.." && pwd)" @@ -92,20 +94,24 @@ LCOV_TOLERANT=( --rc lcov_branch_coverage=1 ) +# C2 fix: capture and extract must succeed (or we have no valid trace). +# The --ignore-errors flags above suppress the gcov-version noise; any +# remaining error is real (e.g. no .gcda files, compiler mismatch) and +# should be visible as a non-zero exit. lcov --capture --directory "$BUILD_DIR" \ --output-file "$BUILD_DIR/coverage.raw.info" \ --no-external \ "${LCOV_TOLERANT[@]}" \ - >/dev/null 2>&1 || true + >/dev/null 2>&1 # Restrict to code/include/ (our public API surface; ignore deps/tests). lcov --extract "$BUILD_DIR/coverage.raw.info" \ "*/code/include/*" \ --output-file "$BUILD_DIR/coverage.info" \ "${LCOV_TOLERANT[@]}" \ - >/dev/null 2>&1 || true + >/dev/null 2>&1 -# ── HTML report ────────────────────────────────────────────────────────────── +# ── HTML report (best-effort; failure here does not fail the script) ────────── genhtml --branch-coverage --legend \ --output-directory "$BUILD_DIR/lcov-html" \ "${LCOV_TOLERANT[@]}" \ @@ -114,15 +120,69 @@ genhtml --branch-coverage --legend \ # ── Summary to stdout ──────────────────────────────────────────────────────── echo echo "── Coverage summary (code/include/) ─────────────────────────" -if [ -s "$BUILD_DIR/coverage.info" ]; then - lcov --summary "$BUILD_DIR/coverage.info" "${LCOV_TOLERANT[@]}" 2>/dev/null \ - | grep -E "lines\.\.\.\.|functions|branches" \ - | sed 's/^/ /' - echo - echo "HTML report: $BUILD_DIR/lcov-html/index.html" - echo " open $BUILD_DIR/lcov-html/index.html" -else - echo " WARN: coverage.info is empty — likely an lcov/gcov version" - echo " mismatch. Tests passed; raw .gcda files are in $BUILD_DIR." - echo " Inspect with: find $BUILD_DIR -name '*.gcda' | head" +if [ ! -s "$BUILD_DIR/coverage.info" ]; then + echo " FAIL: coverage.info is empty — likely an lcov/gcov version" >&2 + echo " mismatch. Tests passed; raw .gcda files are in $BUILD_DIR." >&2 + echo " Inspect with: find $BUILD_DIR -name '*.gcda' | head" >&2 + exit 3 fi + +lcov --summary "$BUILD_DIR/coverage.info" "${LCOV_TOLERANT[@]}" 2>/dev/null \ + | grep -E "lines\.\.\.\.|functions|branches" \ + | sed 's/^/ /' +echo +echo "HTML report: $BUILD_DIR/lcov-html/index.html" +echo " open $BUILD_DIR/lcov-html/index.html" +echo + +# ── Coverage threshold gate (C3) ───────────────────────────────────────────── +# Agreed thresholds (2026-05-31): 80 % line / 70 % branch / 90 % function. +# Set SKIP_COVERAGE_GATE=1 to print without failing (useful during ramp-up). +THRESHOLD_LINE=80 +THRESHOLD_BRANCH=70 +THRESHOLD_FUNC=90 + +if [ "${SKIP_COVERAGE_GATE:-0}" = "1" ]; then + echo "NOTE: SKIP_COVERAGE_GATE=1 — thresholds not enforced." + exit 0 +fi + +SUMMARY=$(lcov --summary "$BUILD_DIR/coverage.info" "${LCOV_TOLERANT[@]}" 2>/dev/null) + +extract_pct() { + echo "$SUMMARY" | grep -i "$1" | grep -oE '[0-9]+\.[0-9]+' | head -1 +} + +LINE_PCT=$(extract_pct "lines") +BRANCH_PCT=$(extract_pct "branches") +FUNC_PCT=$(extract_pct "functions") + +FAIL=0 +check_threshold() { + local label="$1" actual="$2" threshold="$3" + if [ -z "$actual" ]; then + echo " WARN: could not parse $label coverage — skipping gate" >&2 + return + fi + # Use awk for float comparison (bash cannot do floats) + if awk "BEGIN { exit ($actual >= $threshold) ? 0 : 1 }"; then + printf " ✓ %-10s %s%% >= %s%%\n" "$label" "$actual" "$threshold" + else + printf " ✗ %-10s %s%% < %s%% (threshold: %s%%)\n" \ + "$label" "$actual" "$threshold" "$threshold" >&2 + FAIL=1 + fi +} + +echo "── Coverage gate ─────────────────────────────────────────────" +check_threshold "lines" "$LINE_PCT" "$THRESHOLD_LINE" +check_threshold "branches" "$BRANCH_PCT" "$THRESHOLD_BRANCH" +check_threshold "functions" "$FUNC_PCT" "$THRESHOLD_FUNC" +echo + +if [ "$FAIL" -eq 1 ]; then + echo "FAIL: coverage below threshold — see above." >&2 + echo " To suppress (ramp-up): SKIP_COVERAGE_GATE=1 bash $0" >&2 + exit 4 +fi +echo "PASS: all coverage thresholds met."