From 202b9a108d6e3dbb75a49ac981fdbd7d7554b026 Mon Sep 17 00:00:00 2001 From: Tarik Moussa Date: Sun, 31 May 2026 16:27:48 +0200 Subject: [PATCH] feat(num): N3 selectable scale-floor clamp mode (HardJava | SmoothBarrier) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The hyper-ideal vertex scale b is floored to keep the geometry valid. The original clamp `b<0 → 0.01` mirrors the Java oracle but is only C⁰ (in fact value-discontinuous at b=0): a Newton step crossing the feasibility boundary hits a kink that can stall convergence (numerical-stability audit N3). Rather than replace the Java-faithful behaviour (which would break the golden parity tests), make the floor a selectable mode so BOTH the Java standpoint and the clean mathematics are available: - HyperIdealScaleClamp::HardJava (DEFAULT) — the original snap, bit-for-bit faithful to HyperIdealFunctional.java → all parity tests unchanged. - HyperIdealScaleClamp::SmoothBarrier — C¹ softplus floor b ↦ floor + softplus_β(b−floor), β = HYPER_IDEAL_SCALE_SHARPNESS (=100); ≈ identity away from the floor, smooth across b=0. Opt-in. clamp_hyper_ideal_scale centralises the logic (also folds in the N4 nachzügler: compute_face_angles used a bare 0.01). The mode threads with a defaulted trailing parameter through compute_face_angles, face_angles_from_local_dofs, evaluate_hyper_ideal, the four hyper_ideal_hessian* variants and newton_hyper_ideal — so every existing call site keeps HardJava behaviour. Tests (+4): clamp-function C¹/floor/identity contract, mode-equivalence away from the boundary, and end-to-end SmoothBarrier convergence to the same Java golden vector (LawsonHyperIdeal). 296/296 CGAL tests pass. Documented in doc/math/geometry-modes.md. Co-Authored-By: Claude Opus 4.8 --- code/include/constants.hpp | 7 ++ code/include/hyper_ideal_functional.hpp | 72 +++++++++++++++--- code/include/hyper_ideal_hessian.hpp | 24 +++--- code/include/newton_solver.hpp | 17 +++-- code/tests/cgal/test_hyper_ideal_hessian.cpp | 78 ++++++++++++++++++++ code/tests/cgal/test_lawson_hyperideal.cpp | 48 ++++++++++++ doc/math/geometry-modes.md | 23 +++++- 7 files changed, 242 insertions(+), 27 deletions(-) diff --git a/code/include/constants.hpp b/code/include/constants.hpp index aac6ef0..1068f31 100644 --- a/code/include/constants.hpp +++ b/code/include/constants.hpp @@ -29,6 +29,13 @@ constexpr double LOG_EDGE_LENGTH_FLOOR = -30.0; /// 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. diff --git a/code/include/hyper_ideal_functional.hpp b/code/include/hyper_ideal_functional.hpp index 7efcb1f..8533135 100644 --- a/code/include/hyper_ideal_functional.hpp +++ b/code/include/hyper_ideal_functional.hpp @@ -158,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 { @@ -176,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 = HYPER_IDEAL_SCALE_FLOOR; - if (v2b && b2 < 0.0) b2 = HYPER_IDEAL_SCALE_FLOOR; - if (v3b && b3 < 0.0) b3 = HYPER_IDEAL_SCALE_FLOOR; + 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); @@ -247,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); @@ -277,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); @@ -391,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; @@ -407,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/newton_solver.hpp b/code/include/newton_solver.hpp index 686411a..ced7fa2 100644 --- a/code/include/newton_solver.hpp +++ b/code/include/newton_solver.hpp @@ -400,6 +400,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. @@ -410,7 +416,8 @@ 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()); @@ -422,7 +429,7 @@ inline NewtonResult newton_hyper_ideal( for (int iter = 0; iter < max_iter; ++iter) { // ── Gradient ────────────────────────────────────────────────────────── - auto G_std = evaluate_hyper_ideal(mesh, x, m, /*energy=*/false).gradient; + auto G_std = evaluate_hyper_ideal(mesh, x, m, /*energy=*/false, /*grad=*/true, clamp).gradient; Eigen::Map G(G_std.data(), n); double inf_norm = G.cwiseAbs().maxCoeff(); @@ -435,7 +442,7 @@ inline NewtonResult newton_hyper_ideal( } // ── Hessian (block-FD, ~33–1166× faster than full-FD) + solve H·Δx = −G - auto H = hyper_ideal_hessian_block_fd_sym(mesh, x, m, hess_eps); + auto H = hyper_ideal_hessian_block_fd_sym(mesh, x, m, hess_eps, clamp); bool ok = false; Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok); if (!ok) break; @@ -446,14 +453,14 @@ inline NewtonResult newton_hyper_ideal( 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; + return evaluate_hyper_ideal(mesh, xnew, m, false, true, clamp).gradient; }, &improved); if (!improved) break; res.iterations = iter + 1; } - auto G_final = evaluate_hyper_ideal(mesh, x, m, false).gradient; + auto G_final = evaluate_hyper_ideal(mesh, x, m, false, true, clamp).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; 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/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).