conformallab/ConformalLabpp
3
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 1 Wiki Activity

feat(num): N3 selectable scale-floor clamp mode (HardJava | SmoothBarrier)

Browse Source 
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 <noreply@anthropic.com>
This commit is contained in:
Tarik Moussa
2026-05-31 16:27:48 +02:00
parent b7c7d108b5
commit 186f11b11a
7 changed files with 242 additions and 27 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
code/include/constants.hpp
View File

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

72
code/include/hyper_ideal_functional.hpp
View File

@@ -153,6 +153,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 {
@@ -171,15 +220,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);
@@ -242,7 +292,8 @@ static FaceAngles compute_face_angles(
const ConformalMesh& mesh,
Face_index f,
const std::vector<double>& x,
const HyperIdealMaps& m)
const HyperIdealMaps& m,
HyperIdealScaleClamp clamp = HyperIdealScaleClamp::HardJava)
{
Halfedge_index h0 = mesh.halfedge(f);
Halfedge_index h1 = mesh.next(h0);
@@ -272,9 +323,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);
@@ -386,7 +437,8 @@ inline HyperIdealResult evaluate_hyper_ideal(
const std::vector<double>& x,
const HyperIdealMaps& m,
bool need_energy = true,
bool need_gradient = true)
bool need_gradient = true,
HyperIdealScaleClamp clamp = HyperIdealScaleClamp::HardJava)
{
HyperIdealResult res;
@@ -402,7 +454,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.

24
code/include/hyper_ideal_hessian.hpp
View File

@@ -64,7 +64,8 @@ inline Eigen::SparseMatrix<double> hyper_ideal_hessian(
ConformalMesh& mesh,
const std::vector<double>& 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<Eigen::Triplet<double>> trips;
@@ -77,8 +78,8 @@ inline Eigen::SparseMatrix<double> 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<double> hyper_ideal_hessian_sym(
ConformalMesh& mesh,
const std::vector<double>& 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<double> Ht = H.transpose();
return (H + Ht) * 0.5;
}
@@ -142,7 +144,8 @@ inline Eigen::SparseMatrix<double> hyper_ideal_hessian_block_fd(
ConformalMesh& mesh,
const std::vector<double>& 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<Eigen::Triplet<double>> trips;
@@ -187,9 +190,9 @@ inline Eigen::SparseMatrix<double> 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<double> hyper_ideal_hessian_block_fd_sym(
ConformalMesh& mesh,
const std::vector<double>& 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<double> Ht = H.transpose();
return (H + Ht) * 0.5;
}

17
code/include/newton_solver.hpp
View File

@@ -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<double> x = x0;
const int n = static_cast<int>(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<const Eigen::VectorXd> 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<double>& 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;