202b9a108d
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>
236 lines
11 KiB
C++
236 lines
11 KiB
C++
#pragma once
|
||
// Copyright (c) 2024-2026 Tarik Moussa.
|
||
// SPDX-License-Identifier: MIT
|
||
|
||
// hyper_ideal_hessian.hpp
|
||
//
|
||
// Phase 4a — Hessian of the hyper-ideal discrete conformal functional.
|
||
// Phase 9b — Block-finite-difference Hessian (intermediate optimisation).
|
||
//
|
||
// ┌──────────────────────────────────────────────────────────────────────────┐
|
||
// │ Implementation strategy │
|
||
// │ │
|
||
// │ The hyper-ideal functional involves angle functions (ζ, σ, α, β) │
|
||
// │ composed through several nested layers (lij → ζ13/14/15 → β/α). │
|
||
// │ Deriving closed-form Hessian entries analytically through all these │
|
||
// │ layers is feasible but lengthy and is deferred to a future PR. │
|
||
// │ │
|
||
// │ TWO Hessian implementations are provided here: │
|
||
// │ │
|
||
// │ 1. `hyper_ideal_hessian` — full finite-difference baseline. │
|
||
// │ Cost ≈ n × (cost of full gradient evaluation) │
|
||
// │ = O(n · F) where n = #DOFs and F = #faces. │
|
||
// │ Used for correctness reference and small meshes. │
|
||
// │ │
|
||
// │ 2. `hyper_ideal_hessian_block_fd` — block-local finite-difference, │
|
||
// │ Phase 9b. Exploits the fact that each face contributes to the │
|
||
// │ gradient through exactly 6 DOFs (3 vertex b_i + 3 edge a_e). │
|
||
// │ Cost ≈ F × 6 × (cost of a single face-angle evaluation) │
|
||
// │ = O(36 · F). │
|
||
// │ Speed-up factor ≈ n / 36, i.e. typically 10–50× on V > 200. │
|
||
// │ │
|
||
// │ Both produce the same Hessian to O(ε²) and pass identical PSD checks. │
|
||
// │ The block-FD variant is the production default; the full-FD variant is │
|
||
// │ kept for cross-validation tests. │
|
||
// │ │
|
||
// │ An analytic Hessian via Schläfli-type differentiation through the chain │
|
||
// │ (bᵢ, aₑ) → lᵢⱼ → ζ₁₃/ζ₁₄/ζ₁₅ → αᵢⱼ / βᵢ │
|
||
// │ is deferred to a future PR (Phase 9b-analytic). Speed-up would be │
|
||
// │ another ~6×, taking the cost to O(F). │
|
||
// │ │
|
||
// │ The hyper-ideal energy is strictly convex (Springborn 2020), so H is │
|
||
// │ positive semi-definite everywhere and Eigen::SimplicialLDLT applies │
|
||
// │ directly to either Hessian variant. │
|
||
// └──────────────────────────────────────────────────────────────────────────┘
|
||
//
|
||
// Note on the Java reference: HyperIdealFunctional.java line 295-298 declares
|
||
// public boolean hasHessian() { return false; }
|
||
// — i.e. the upstream Java implementation supplies NO Hessian, analytic or
|
||
// numerical. Both `hyper_ideal_hessian` and `hyper_ideal_hessian_block_fd`
|
||
// are conformallab++ additions beyond Java parity.
|
||
|
||
#include "hyper_ideal_functional.hpp"
|
||
#include <Eigen/Sparse>
|
||
#include <vector>
|
||
#include <cmath>
|
||
#include <cstdint>
|
||
|
||
namespace conformallab {
|
||
|
||
/// Full finite-difference HyperIdeal Hessian (baseline, Phase 4a).
|
||
/// 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<double> hyper_ideal_hessian(
|
||
ConformalMesh& mesh,
|
||
const std::vector<double>& x,
|
||
const HyperIdealMaps& m,
|
||
double eps = 1e-5,
|
||
HyperIdealScaleClamp clamp = HyperIdealScaleClamp::HardJava)
|
||
{
|
||
const int n = hyper_ideal_dimension(mesh, m);
|
||
std::vector<Eigen::Triplet<double>> trips;
|
||
trips.reserve(static_cast<std::size_t>(n * n));
|
||
|
||
std::vector<double> xp = x, xm = x;
|
||
|
||
for (int j = 0; j < n; ++j) {
|
||
const std::size_t sj = static_cast<std::size_t>(j);
|
||
xp[sj] = x[sj] + eps;
|
||
xm[sj] = x[sj] - eps;
|
||
|
||
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
|
||
|
||
for (int i = 0; i < n; ++i) {
|
||
double val = (Gp[static_cast<std::size_t>(i)]
|
||
- Gm[static_cast<std::size_t>(i)]) / (2.0 * eps);
|
||
if (std::abs(val) > 1e-15)
|
||
trips.emplace_back(i, j, val);
|
||
}
|
||
}
|
||
|
||
Eigen::SparseMatrix<double> H(n, n);
|
||
H.setFromTriplets(trips.begin(), trips.end());
|
||
return H;
|
||
}
|
||
|
||
/// Symmetrised full-FD HyperIdeal Hessian: returns `(H + Hᵀ) / 2` to
|
||
/// scrub the tiny asymmetries introduced by floating-point rounding.
|
||
inline Eigen::SparseMatrix<double> hyper_ideal_hessian_sym(
|
||
ConformalMesh& mesh,
|
||
const std::vector<double>& x,
|
||
const HyperIdealMaps& m,
|
||
double eps = 1e-5,
|
||
HyperIdealScaleClamp clamp = HyperIdealScaleClamp::HardJava)
|
||
{
|
||
auto H = hyper_ideal_hessian(mesh, x, m, eps, clamp);
|
||
Eigen::SparseMatrix<double> Ht = H.transpose();
|
||
return (H + Ht) * 0.5;
|
||
}
|
||
|
||
// ── Block-FD Hessian (Phase 9b) ──────────────────────────────────────────────
|
||
//
|
||
// Computes the Hessian by FD on each face's 6×6 local block. The 6 local
|
||
// DOFs of a face f are:
|
||
// (b_{v1}, b_{v2}, b_{v3}, a_{e12}, a_{e23}, a_{e31}).
|
||
// For each face we recompute the 6 output angles (β₁,β₂,β₃,α₁₂,α₂₃,α₃₁)
|
||
// at x ± ε along each local axis and read off the 6×6 Jacobian. The result
|
||
// scatters into the global Hessian via the DOF-index lookup.
|
||
//
|
||
// Why this is correct:
|
||
// ─────────────────────
|
||
// The global gradient decomposes by face:
|
||
// G_b_v = Σ_{f ∋ v} β_v(f) − Θ_v
|
||
// G_a_e = Σ_{f ∋ e} α_e(f) − θ_e
|
||
// Since β and α at face f depend ONLY on the 6 local DOFs of f, the
|
||
// Hessian also decomposes:
|
||
// ∂G_x/∂y = Σ_{f: x,y ∈ local(f)} ∂(β or α)/∂y at f.
|
||
// So accumulating per-face 6×6 blocks reproduces the full Hessian.
|
||
//
|
||
// Cost: F × 12 face-angle evaluations (6 DOFs × 2 directions).
|
||
// On a tetrahedron (F=4, n≈10): 48 face evaluations
|
||
// vs full-FD ≈ 80 → ~1.7× speed-up.
|
||
// On cathead.obj (F=248, n≈400): 2976 face evaluations
|
||
// vs full-FD ≈ 99,200 → ~33× speed-up.
|
||
// On brezel.obj (F=13824, n≈14000): 165 888 face evaluations
|
||
// vs full-FD ≈ 193 M → ~1166× speed-up.
|
||
/// Per-face block-FD HyperIdeal Hessian (Phase 9b). Uses the locality
|
||
/// lemma `∂G_x/∂y = Σ_{f: x,y ∈ local(f)} ∂(β or α)/∂y` to perturb only
|
||
/// the 6 face-local DOFs at a time, giving an `F·12` face-evaluation
|
||
/// budget vs `n·F` for full-FD (~96× speed-up on brezel.obj).
|
||
inline Eigen::SparseMatrix<double> hyper_ideal_hessian_block_fd(
|
||
ConformalMesh& mesh,
|
||
const std::vector<double>& x,
|
||
const HyperIdealMaps& m,
|
||
double eps = 1e-5,
|
||
HyperIdealScaleClamp clamp = HyperIdealScaleClamp::HardJava)
|
||
{
|
||
const int n = hyper_ideal_dimension(mesh, m);
|
||
std::vector<Eigen::Triplet<double>> trips;
|
||
trips.reserve(36 * 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: (b1, b2, b3, a12, a23, a31). Pinned slots = -1.
|
||
const int idx[6] = {
|
||
m.v_idx[v1], m.v_idx[v2], m.v_idx[v3],
|
||
m.e_idx[e12], m.e_idx[e23], m.e_idx[e31]
|
||
};
|
||
const bool v1b = idx[0] >= 0;
|
||
const bool v2b = idx[1] >= 0;
|
||
const bool v3b = idx[2] >= 0;
|
||
|
||
// Local DOF values (0 for pinned).
|
||
const double vals[6] = {
|
||
dof_val(idx[0], x), dof_val(idx[1], x), dof_val(idx[2], x),
|
||
dof_val(idx[3], x), dof_val(idx[4], x), dof_val(idx[5], x)
|
||
};
|
||
|
||
// For each free local DOF, evaluate the 6 outputs at ±ε.
|
||
// We never perturb a pinned DOF (its column would be physically zero
|
||
// because it is not part of the DOF vector at all).
|
||
for (int j = 0; j < 6; ++j) {
|
||
if (idx[j] < 0) continue;
|
||
|
||
double vp[6], vm[6];
|
||
for (int k = 0; k < 6; ++k) { vp[k] = vm[k] = vals[k]; }
|
||
vp[j] += eps;
|
||
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, clamp);
|
||
auto Om = face_angles_from_local_dofs(
|
||
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,
|
||
Op.alpha12, Op.alpha23, Op.alpha31
|
||
};
|
||
const double Gm[6] = {
|
||
Om.beta1, Om.beta2, Om.beta3,
|
||
Om.alpha12, Om.alpha23, Om.alpha31
|
||
};
|
||
|
||
for (int i = 0; i < 6; ++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<double> H(n, n);
|
||
H.setFromTriplets(trips.begin(), trips.end());
|
||
return H;
|
||
}
|
||
|
||
/// Symmetrised block-FD HyperIdeal Hessian: returns `(H + Hᵀ) / 2` of
|
||
/// `hyper_ideal_hessian_block_fd(...)` for downstream solvers that
|
||
/// require strict symmetry.
|
||
inline Eigen::SparseMatrix<double> hyper_ideal_hessian_block_fd_sym(
|
||
ConformalMesh& mesh,
|
||
const std::vector<double>& x,
|
||
const HyperIdealMaps& m,
|
||
double eps = 1e-5,
|
||
HyperIdealScaleClamp clamp = HyperIdealScaleClamp::HardJava)
|
||
{
|
||
auto H = hyper_ideal_hessian_block_fd(mesh, x, m, eps, clamp);
|
||
Eigen::SparseMatrix<double> Ht = H.transpose();
|
||
return (H + Ht) * 0.5;
|
||
}
|
||
|
||
} // namespace conformallab
|