From 822a27da692967b051c0182c50251dcc481aa921 Mon Sep 17 00:00:00 2001 From: Tarik Moussa Date: Sat, 30 May 2026 00:18:41 +0200 Subject: [PATCH] feat(euclidean): full analytic edge-DOF (cyclic) Hessian; Tier-2 Wente finding MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Upgrades the cyclic Euclidean Hessian from block-FD to closed-form, satisfying the novelty-statement §3.2 "analytic Hessians, not finite difference" claim for the Euclidean path. - euclidean_hessian.hpp: `euclidean_hessian_analytic` — closed-form cyclic Hessian from the law-of-cosines angle derivatives ∂α_i/∂s_i = ℓ_i²/4A, ∂α_i/∂s_j = ½cot α_i − ℓ_j²/4A (Σ_j = 0), chained to (u, λ_e) and sign-mapped to the gradient outputs (−α vertex, +α_opp edge). Reuses euclidean_cot_weights. Block-FD kept as cross-check. - newton_solver.hpp: newton_euclidean cyclic path now uses the analytic Hessian. - tests: CyclicHessian_Analytic_MatchesBlockFD_Tetrahedron — analytic == block-FD (1e-6), == gradient FD (1e-5), symmetric (1e-9). Existing cyclic convergence oracle still GREEN with the analytic Hessian routed in. Tier-2 (Wente) finding: wente_torus02.obj is a QUAD mesh (1240 quads) and the Java golden comes from cyclic (quad-net) uniformization; the C++ period-matrix pipeline is triangle-based, so a faithful bit-vs-Java τ comparison needs a quad/cyclic pipeline (Phase 9f). Deferred and documented; golden τ = ½+i√3/2. 244/244 cgal tests pass. Co-Authored-By: Claude Opus 4.8 --- code/include/euclidean_hessian.hpp | 116 ++++++++++++++++++ code/include/newton_solver.hpp | 8 +- code/tests/cgal/test_euclidean_functional.cpp | 58 +++++++++ doc/reviewer/java-ignore-crossvalidation.md | 23 +++- 4 files changed, 198 insertions(+), 7 deletions(-) diff --git a/code/include/euclidean_hessian.hpp b/code/include/euclidean_hessian.hpp index bd39924..5918c79 100644 --- a/code/include/euclidean_hessian.hpp +++ b/code/include/euclidean_hessian.hpp @@ -298,6 +298,122 @@ inline Eigen::SparseMatrix euclidean_hessian_block_fd_sym( return (H + Ht) * 0.5; } +// ── Analytic cyclic Hessian (vertex + edge DOFs) ───────────────────────────── +// +// Closed-form Jacobian of the cyclic gradient. With s_i = effective log-length +// (Λ̃) of the edge OPPOSITE vertex i (s = 2·ln ℓ), the Euclidean corner-angle +// derivatives are (derived from the law of cosines; Σ_j ∂α_i/∂s_j = 0): +// +// ∂α_i/∂s_i = ℓ_i² / (4A) +// ∂α_i/∂s_j = ½ cot α_i − ℓ_j² / (4A) (j ≠ i) +// +// where A is the triangle area (4A = denom2/2, denom2 = 8A from the cot helper). +// Chaining through s₁ = Λ̃(e₂₃), s₂ = Λ̃(e₃₁), s₃ = Λ̃(e₁₂) and +// Λ̃_e = λ⁰_e + u_i + u_j + λ_e gives the per-face 6×6 block over the local +// DOFs (u₁,u₂,u₃,λ₁₂,λ₂₃,λ₃₁); per-face outputs carry the gradient signs +// (−α vertex, +α_opp edge), so the assembled matrix equals ∂G/∂x exactly. +// This is the analytic counterpart of `euclidean_hessian_block_fd` (no FD). + +/// Analytic Euclidean cyclic Hessian (vertex + edge DOFs), sparse. +/// Closed-form (no finite differences); matches `euclidean_hessian_block_fd` +/// to round-off and the analytic vertex-only Laplacian on vertex-only layouts. +inline Eigen::SparseMatrix euclidean_hessian_analytic( + ConformalMesh& mesh, + const std::vector& x, + const EuclideanMaps& m) +{ + const int n = euclidean_dimension(mesh, m); + std::vector> trips; + trips.reserve(static_cast(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); + + const double u1 = eucl_dof_val(m.v_idx[v1], x); + const double u2 = eucl_dof_val(m.v_idx[v2], x); + const double u3 = eucl_dof_val(m.v_idx[v3], x); + + const double lam12 = m.lambda0[e12] + u1 + u2 + eucl_dof_val(m.e_idx[e12], x); // s3 + const double lam23 = m.lambda0[e23] + u2 + u3 + eucl_dof_val(m.e_idx[e23], x); // s1 + const double lam31 = m.lambda0[e31] + u3 + u1 + eucl_dof_val(m.e_idx[e31], x); // s2 + + // Centered side lengths (scale-invariant for ℓ²/4A and cot). + const double mu = (lam12 + lam23 + lam31) / 6.0; + const double l12 = std::exp((lam12 - 2.0 * mu) * 0.5); + const double l23 = std::exp((lam23 - 2.0 * mu) * 0.5); + const double l31 = std::exp((lam31 - 2.0 * mu) * 0.5); + + const double t12 = -l12 + l23 + l31; + const double t23 = l12 - l23 + l31; + const double t31 = l12 + l23 - l31; + if (t12 <= 0.0 || t23 <= 0.0 || t31 <= 0.0) continue; + const double l123 = l12 + l23 + l31; + const double denom2_sq = t12 * t23 * t31 * l123; + if (denom2_sq <= 0.0) continue; + const double denom2 = 2.0 * std::sqrt(denom2_sq); // = 8A + const double fourA = denom2 * 0.5; // = 4A + + auto cw = euclidean_cot_weights(l12, l23, l31); + if (!cw.valid) continue; + const double cot[3] = { cw.cot1, cw.cot2, cw.cot3 }; // at v1, v2, v3 + + // ℓ_i² = squared length of edge opposite vertex i: + // opp(v1)=e23=l23, opp(v2)=e31=l31, opp(v3)=e12=l12 + const double Lsq[3] = { l23 * l23, l31 * l31, l12 * l12 }; + + // dA[i][j] = ∂α_i/∂s_j (i,j ∈ {0,1,2}; s_j opposite v_{j+1}) + double dA[3][3]; + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + dA[i][j] = (i == j) ? (Lsq[i] / fourA) + : (0.5 * cot[i] - Lsq[j] / fourA); + + // ∂α_i/∂(local dof), dof order (u1,u2,u3, λ12,λ23,λ31). + // u1 ↔ s2,s3 ; u2 ↔ s1,s3 ; u3 ↔ s1,s2 ; λ12 ↔ s3 ; λ23 ↔ s1 ; λ31 ↔ s2. + double g[3][6]; + for (int i = 0; i < 3; ++i) { + g[i][0] = dA[i][1] + dA[i][2]; // u1 + g[i][1] = dA[i][0] + dA[i][2]; // u2 + g[i][2] = dA[i][0] + dA[i][1]; // u3 + g[i][3] = dA[i][2]; // λ12 (s3) + g[i][4] = dA[i][0]; // λ23 (s1) + g[i][5] = dA[i][1]; // λ31 (s2) + } + + // Output slot → (angle index, sign): (−α1,−α2,−α3,+α3,+α1,+α2). + const int out_a[6] = { 0, 1, 2, 2, 0, 1 }; + const double out_sg[6] = { -1.0, -1.0, -1.0, +1.0, +1.0, +1.0 }; + + 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] + }; + + for (int oi = 0; oi < 6; ++oi) { + if (idx[oi] < 0) continue; + for (int dj = 0; dj < 6; ++dj) { + if (idx[dj] < 0) continue; + const double val = out_sg[oi] * g[out_a[oi]][dj]; + if (std::abs(val) > 1e-15) + trips.emplace_back(idx[oi], idx[dj], val); + } + } + } + + Eigen::SparseMatrix H(n, n); + H.setFromTriplets(trips.begin(), trips.end()); + return H; +} + // ── Finite-difference Hessian check ────────────────────────────────────────── /// FD Hessian check for the Euclidean functional. Compares analytic /// `H` column-by-column to `(G(x+εeⱼ) − G(x−εeⱼ)) / (2ε)`; returns diff --git a/code/include/newton_solver.hpp b/code/include/newton_solver.hpp index fb77062..b882e46 100644 --- a/code/include/newton_solver.hpp +++ b/code/include/newton_solver.hpp @@ -258,13 +258,13 @@ inline NewtonResult newton_euclidean( } // ── Hessian + solve H·Δx = −G (SparseQR fallback for singular H) ── - // Cyclic layout (edge DOFs present) → block-FD Hessian, which covers the - // vertex-edge / edge-edge blocks the analytic cotangent Laplacian omits. - // Vertex-only layout → analytic cotangent Laplacian (cheaper, exact). + // 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_block_fd_sym(mesh, x, m) + 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); diff --git a/code/tests/cgal/test_euclidean_functional.cpp b/code/tests/cgal/test_euclidean_functional.cpp index 1c9f9b2..9857ae5 100644 --- a/code/tests/cgal/test_euclidean_functional.cpp +++ b/code/tests/cgal/test_euclidean_functional.cpp @@ -587,3 +587,61 @@ TEST(EuclideanFunctional, CyclicHessian_BlockFD_MatchesGradientFD_Tetrahedron) EXPECT_LT(max_err, 1e-5) << "block-FD cyclic Hessian disagrees with the gradient finite difference"; } + +// ════════════════════════════════════════════════════════════════════════════ +// Analytic cyclic Hessian == block-FD (and == gradient FD) +// +// `euclidean_hessian_analytic` is the closed-form counterpart of +// `euclidean_hessian_block_fd`. On the full cyclic layout they must agree to +// round-off, and both must match the gradient finite difference. This validates +// the analytic derivation (∂α_i/∂s_i = ℓ_i²/4A, ∂α_i/∂s_j = ½cot α_i − ℓ_j²/4A). +// ════════════════════════════════════════════════════════════════════════════ +TEST(EuclideanFunctional, CyclicHessian_Analytic_MatchesBlockFD_Tetrahedron) +{ + auto mesh = make_tetrahedron(); + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + const int n = assign_euclidean_all_dof_indices(mesh, maps); // 4 + 6 = 10 + + std::vector x(static_cast(n), 0.0); + for (int i = 0; i < 4; ++i) x[static_cast(i)] = -0.15; + for (int i = 4; i < n; ++i) x[static_cast(i)] = 0.05; + + auto Ha = euclidean_hessian_analytic(mesh, x, maps); + auto Hb = euclidean_hessian_block_fd(mesh, x, maps); + + // Analytic vs block-FD. + double max_ab = 0.0; + for (int i = 0; i < n; ++i) + for (int j = 0; j < n; ++j) + max_ab = std::max(max_ab, std::abs(Ha.coeff(i, j) - Hb.coeff(i, j))); + EXPECT_LT(max_ab, 1e-6) + << "analytic cyclic Hessian disagrees with block-FD"; + + // Analytic vs gradient FD (independent ground truth). + const double eps = 1e-6; + std::vector xp = x, xm = x; + double max_af = 0.0; + 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 = euclidean_gradient(mesh, xp, maps); + auto Gm = euclidean_gradient(mesh, xm, maps); + xp[sj] = xm[sj] = x[sj]; + for (int i = 0; i < n; ++i) { + const double fd = (Gp[static_cast(i)] + - Gm[static_cast(i)]) / (2.0 * eps); + max_af = std::max(max_af, std::abs(Ha.coeff(i, j) - fd)); + } + } + EXPECT_LT(max_af, 1e-5) + << "analytic cyclic Hessian disagrees with the gradient finite difference"; + + // Symmetry (Hessian of a scalar energy). + double max_sym = 0.0; + for (int i = 0; i < n; ++i) + for (int j = 0; j < n; ++j) + max_sym = std::max(max_sym, std::abs(Ha.coeff(i, j) - Ha.coeff(j, i))); + EXPECT_LT(max_sym, 1e-9) << "analytic cyclic Hessian is not symmetric"; +} diff --git a/doc/reviewer/java-ignore-crossvalidation.md b/doc/reviewer/java-ignore-crossvalidation.md index 65176f8..f821b81 100644 --- a/doc/reviewer/java-ignore-crossvalidation.md +++ b/doc/reviewer/java-ignore-crossvalidation.md @@ -151,9 +151,26 @@ cotangent Laplacian. Spherical convergence (Tier-1 #2) was found **already covered** by `test_newton_solver` (`Spherical_ConvergesFromPerturbation` et al.), and `make_octahedron_face()` is a single face, not a closed octahedron — so it adds -little. → next genuinely-new target: **Tier 2 (Wente uniformization XML)**. A -full *analytic* edge-DOF Hessian (vs the current block-FD) remains a future -optimisation. +little. + +### Tier-2 (Wente) finding (2026-05-30) — needs a quad/cyclic pipeline + +`wente_torus02.obj` is a **quad mesh** (1240 quads, genus 1), and the Java +`wente_uniformization.xml` golden (uniformizing group → τ = (√2+i√6)/(2√2) = +½ + i·√3/2, the hexagonal torus) comes from Java's **cyclic (quad) net** +uniformization. The C++ period-matrix pipeline (`run_torus_pipeline`) is +**triangle-based**, so a faithful bit-vs-Java comparison would require triangulating +(which changes the discrete conformal structure → no exact match) or a +**quad/cyclic period-matrix pipeline** — that belongs to **Phase 9f** (general +polygon meshes), not the triangle path. → Tier 2 deferred to Phase 9f; the +hexagonal golden τ is recorded here for when that lands. + +### Done instead (2026-05-30): full analytic edge-DOF Hessian +With Tier-2 gated on a quad pipeline, the higher-value clean win was upgrading +the cyclic Hessian from block-FD to **fully analytic** (closed-form +`∂α_i/∂s_i = ℓ_i²/4A`, `∂α_i/∂s_j = ½cot α_i − ℓ_j²/4A`), satisfying the +`novelty-statement.md §3.2` "analytic Hessians, not finite difference" claim for +the Euclidean path. Verified against block-FD and gradient-FD. ### Other generators in the Java tree `HyperellipticCurveGenerator` (→ Phase 13), `SchottkyGenerator` (→ Phase 11a),