From ea5f01d73d2746133901f7c8cd31601ca2144a4f Mon Sep 17 00:00:00 2001 From: Tarik Moussa Date: Sat, 30 May 2026 00:09:37 +0200 Subject: [PATCH] =?UTF-8?q?feat(euclidean):=20block-FD=20edge-DOF=20Hessia?= =?UTF-8?q?n=20=E2=86=92=20cyclic=20Newton=20+=20Java=20convergence=20orac?= =?UTF-8?q?le?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Implements the edge-DOF (cyclic) Euclidean Hessian, unblocking the full cyclic Newton solve, and enables the Java EuclideanCyclicConvergenceTest cross-validation. - euclidean_hessian.hpp: `euclidean_hessian_block_fd` / `_sym` — per-face 6×6 block FD over (u1,u2,u3,λ12,λ23,λ31), mirroring hyper_ideal_hessian_block_fd. Per-face outputs carry the gradient signs (−α vertex, +α_opp edge), so the result equals ∂G/∂x by construction (locality lemma). Analytic vertex-only cotangent Hessian unchanged (still used for vertex-only layouts). - newton_solver.hpp: newton_euclidean routes cyclic layouts (edge DOFs present) through the block-FD Hessian; vertex-only path unchanged. - tests: * CyclicCircularEdge_CatHead_JavaXVal (now GREEN) — prescribe φ=π−0.1 on one interior edge, solve, assert realised α_opp+α_opp = π−0.1 @1e-9. * CyclicCircularEdge_PhiEntersGradient_CatHead — solver-free φ-wiring check. * CyclicHessian_BlockFD_MatchesGradientFD_Tetrahedron — Hessian correctness. 243/243 cgal tests pass; vertex-only Euclidean Newton unaffected. Co-Authored-By: Claude Opus 4.8 --- code/include/euclidean_hessian.hpp | 111 ++++++++++++++++++ code/include/newton_solver.hpp | 9 +- code/tests/cgal/test_euclidean_functional.cpp | 90 +++++++++----- doc/reviewer/java-ignore-crossvalidation.md | 52 ++++---- 4 files changed, 210 insertions(+), 52 deletions(-) diff --git a/code/include/euclidean_hessian.hpp b/code/include/euclidean_hessian.hpp index e9a61c3..bd39924 100644 --- a/code/include/euclidean_hessian.hpp +++ b/code/include/euclidean_hessian.hpp @@ -43,6 +43,7 @@ #include "euclidean_functional.hpp" #include #include +#include #include #include @@ -187,6 +188,116 @@ inline Eigen::SparseMatrix euclidean_hessian( return H; } +// ── Block-FD Hessian (cyclic: vertex + edge DOFs) ──────────────────────────── +// +// The analytic `euclidean_hessian` above covers only the vertex-block cotangent +// Laplacian. The *cyclic* functional adds edge DOFs λ_e, giving vertex-edge and +// edge-edge Hessian blocks. We obtain the full Hessian by per-face block FD, +// mirroring `hyper_ideal_hessian_block_fd`. +// +// Locality lemma: the gradient decomposes by face, +// G_v = Θ_v − Σ_{f∋v} α_v(f) (vertex output = −α_v) +// G_e = Σ_{f∋e} α_opp(f) − φ_e (edge output = +α_opp) +// and α at face f depends ONLY on f's 6 local DOFs (u₁,u₂,u₃,λ₁₂,λ₂₃,λ₃₁), so +// ∂G_x/∂y = Σ_{f: x,y ∈ local(f)} ∂(output)/∂y at f. +// Accumulating per-face 6×6 blocks therefore reproduces the full Hessian and is +// consistent with `euclidean_gradient` by construction (it is its FD). + +/// Per-face contributions to the cyclic gradient, in slot order +/// (v₁, v₂, v₃, e₁₂, e₂₃, e₃₁). Vertex slots carry −α (since G_v = Θ−Σα), +/// edge slots carry +α_opp (since G_e = Σα_opp − φ). +inline std::array eucl_face_local_outputs( + double u1, double u2, double u3, + double d12, double d23, double d31, + double lam0_12, double lam0_23, double lam0_31) +{ + const double lam12 = lam0_12 + u1 + u2 + d12; + const double lam23 = lam0_23 + u2 + u3 + d23; + const double lam31 = lam0_31 + u3 + u1 + d31; + const auto fa = euclidean_angles(lam12, lam23, lam31); + // edge slot e12 ↔ opposite corner α3, e23 ↔ α1, e31 ↔ α2 (see gradient Pass 3). + return { -fa.alpha1, -fa.alpha2, -fa.alpha3, + +fa.alpha3, +fa.alpha1, +fa.alpha2 }; +} + +/// Block-FD Euclidean cyclic Hessian (vertex + edge DOFs), sparse. +/// Supports both vertex-only and cyclic DOF layouts; cost `F·12` face-angle +/// evaluations. Consistent with `euclidean_gradient` to `O(ε²)`. +inline Eigen::SparseMatrix euclidean_hessian_block_fd( + ConformalMesh& mesh, + const std::vector& x, + const EuclideanMaps& m, + double eps = 1e-5) +{ + 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 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 double lam0_12 = m.lambda0[e12]; + const double lam0_23 = m.lambda0[e23]; + const double lam0_31 = m.lambda0[e31]; + + const double vals[6] = { + eucl_dof_val(idx[0], x), eucl_dof_val(idx[1], x), eucl_dof_val(idx[2], x), + eucl_dof_val(idx[3], x), eucl_dof_val(idx[4], x), eucl_dof_val(idx[5], x) + }; + + for (int j = 0; j < 6; ++j) { + if (idx[j] < 0) continue; // pinned: no column + + 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 = eucl_face_local_outputs( + vp[0], vp[1], vp[2], vp[3], vp[4], vp[5], lam0_12, lam0_23, lam0_31); + auto Om = eucl_face_local_outputs( + vm[0], vm[1], vm[2], vm[3], vm[4], vm[5], lam0_12, lam0_23, lam0_31); + + for (int i = 0; i < 6; ++i) { + if (idx[i] < 0) continue; + const double val = (Op[static_cast(i)] + - Om[static_cast(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 Euclidean cyclic Hessian: `(H + Hᵀ)/2`. +inline Eigen::SparseMatrix euclidean_hessian_block_fd_sym( + ConformalMesh& mesh, + const std::vector& x, + const EuclideanMaps& m, + double eps = 1e-5) +{ + auto H = euclidean_hessian_block_fd(mesh, x, m, eps); + Eigen::SparseMatrix Ht = H.transpose(); + return (H + Ht) * 0.5; +} + // ── 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 c8802f5..fb77062 100644 --- a/code/include/newton_solver.hpp +++ b/code/include/newton_solver.hpp @@ -258,7 +258,14 @@ inline NewtonResult newton_euclidean( } // ── Hessian + solve H·Δx = −G (SparseQR fallback for singular H) ── - auto H = euclidean_hessian(mesh, x, m); + // 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). + 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) + : euclidean_hessian(mesh, x, m); bool ok = false; Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok); if (!ok) break; diff --git a/code/tests/cgal/test_euclidean_functional.cpp b/code/tests/cgal/test_euclidean_functional.cpp index 2f1fc87..1c9f9b2 100644 --- a/code/tests/cgal/test_euclidean_functional.cpp +++ b/code/tests/cgal/test_euclidean_functional.cpp @@ -485,18 +485,10 @@ TEST(EuclideanFunctional, CyclicCircularEdge_PhiEntersGradient_CatHead) // is the prescribed φ); the test therefore FAILS if the solver ignores a // non-default φ (the sum would stay at its natural value, not π − 0.1). // -// ⚠️ DISABLED (build-verified 2026-05-29): blocked by a missing feature, not a -// bug. `newton_euclidean` uses `euclidean_hessian`, which throws -// "euclidean_hessian: edge DOFs are not supported -// (only the vertex-block cotangent Laplacian is implemented)". -// The cyclic functional needs vertex+edge DOFs, so the full Newton solve is not -// yet possible in C++ (the gradient supports edge DOFs; the analytic Hessian -// does not). PREREQUISITE: edge-DOF Euclidean Hessian — see -// `doc/roadmap/research-track.md` and `doc/reviewer/java-ignore-crossvalidation.md`. -// The golden semantics (φ = π−0.1 ⇒ realised α_opp+α_opp = π−0.1) are kept here -// so this test auto-activates once that Hessian lands; just drop the DISABLED_. +// Enabled 2026-05-30: `newton_euclidean` now uses the block-FD edge-DOF Hessian +// (`euclidean_hessian_block_fd_sym`) for cyclic layouts, so the full solve runs. // ════════════════════════════════════════════════════════════════════════════ -TEST(EuclideanFunctional, DISABLED_CyclicCircularEdge_CatHead_JavaXVal) +TEST(EuclideanFunctional, CyclicCircularEdge_CatHead_JavaXVal) { const std::string path = std::string(CONFORMALLAB_DATA_DIR) + "/obj/cathead.obj"; ConformalMesh mesh; @@ -505,26 +497,13 @@ TEST(EuclideanFunctional, DISABLED_CyclicCircularEdge_CatHead_JavaXVal) auto maps = setup_euclidean_maps(mesh); compute_euclidean_lambda0_from_mesh(mesh, maps); - // Cyclic DOFs: interior vertices (border pinned) + all edges. + // DOFs: interior vertices (border pinned) + exactly ONE edge DOF on the + // "circular" edge (Java marks a single circularHoleEdge). Giving every edge + // a DOF would make λ_e redundant with u_i+u_j (a V-dim null space) and stall + // Newton; the single extra edge variable keeps the system well-posed. int idx = 0; for (auto v : mesh.vertices()) maps.v_idx[v] = mesh.is_border(v) ? -1 : idx++; - for (auto e : mesh.edges()) - maps.e_idx[e] = idx++; - const int n = idx; - ASSERT_GT(n, 0); - - // Natural targets: set Θ_v / φ_e so that x = 0 is the equilibrium (G(0)=0). - std::vector x0(static_cast(n), 0.0); - auto G0 = euclidean_gradient(mesh, x0, maps); - for (auto v : mesh.vertices()) { - int iv = maps.v_idx[v]; - if (iv >= 0) maps.theta_v[v] -= G0[static_cast(iv)]; - } - for (auto e : mesh.edges()) { - int ie = maps.e_idx[e]; - if (ie >= 0) maps.phi_e[e] += G0[static_cast(ie)]; - } // Pick one interior edge: both incident faces present, both endpoints interior. Edge_index circular{}; @@ -537,8 +516,21 @@ TEST(EuclideanFunctional, DISABLED_CyclicCircularEdge_CatHead_JavaXVal) circular = e; found = true; break; } ASSERT_TRUE(found) << "no interior edge found on cathead"; + maps.e_idx[circular] = idx++; // the single circular-edge DOF + const int n = idx; + ASSERT_GT(n, 0); const std::size_t ie = static_cast(maps.e_idx[circular]); + // Natural targets: set Θ_v / φ_e so that x = 0 is the equilibrium (G(0)=0), + // so the only deviation is the prescribed circular φ below. + std::vector x0(static_cast(n), 0.0); + auto G0 = euclidean_gradient(mesh, x0, maps); + for (auto v : mesh.vertices()) { + int iv = maps.v_idx[v]; + if (iv >= 0) maps.theta_v[v] -= G0[static_cast(iv)]; + } + maps.phi_e[circular] += G0[ie]; + // Prescribe the circular edge turn angle φ = π − 0.1 (Java CustomEdgeInfo.phi). const double phi_target = PI - 0.1; maps.phi_e[circular] = phi_target; @@ -553,3 +545,45 @@ TEST(EuclideanFunctional, DISABLED_CyclicCircularEdge_CatHead_JavaXVal) EXPECT_NEAR(phi_target, realised, 1e-9) << "prescribed circular edge turn angle π−0.1 not realised at the solution"; } + +// ════════════════════════════════════════════════════════════════════════════ +// Correctness of the block-FD edge-DOF (cyclic) Hessian +// +// Validates `euclidean_hessian_block_fd` directly: on a tetrahedron with the +// full cyclic DOF layout (4 vertex + 6 edge), every entry must match the +// column-wise finite difference of `euclidean_gradient` (the true Jacobian of +// G). This pins the per-face output sign mapping (−α vertex / +α_opp edge) and +// the scatter independently of the convergence test above. +// ════════════════════════════════════════════════════════════════════════════ +TEST(EuclideanFunctional, CyclicHessian_BlockFD_MatchesGradientFD_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; // vertices + for (int i = 4; i < n; ++i) x[static_cast(i)] = 0.05; // edges + + auto H = euclidean_hessian_block_fd(mesh, x, maps); + + const double eps = 1e-6; + std::vector xp = x, xm = x; + double max_err = 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_err = std::max(max_err, std::abs(H.coeff(i, j) - fd)); + } + } + EXPECT_LT(max_err, 1e-5) + << "block-FD cyclic Hessian disagrees with the gradient finite difference"; +} diff --git a/doc/reviewer/java-ignore-crossvalidation.md b/doc/reviewer/java-ignore-crossvalidation.md index c70c0db..65176f8 100644 --- a/doc/reviewer/java-ignore-crossvalidation.md +++ b/doc/reviewer/java-ignore-crossvalidation.md @@ -121,33 +121,39 @@ Assertion strategy: because the values are symmetric per DOF-class, assert XML reader; strongest deterministic pipeline oracle that maps to landed C++). 4. **Tier 3 — Lawson generator** (only after Tiers 1–2). -### Build-verification finding (2026-05-29) — Tier 1 cyclic is blocked +### Build-verification finding (2026-05-29) → resolved (2026-05-30) Attempting the Tier-1 `EuclideanCyclicConvergenceTest` port **build-verified** a -blocker: `newton_euclidean` calls `euclidean_hessian`, which throws -*"edge DOFs are not supported (only the vertex-block cotangent Laplacian is -implemented)"*. The cyclic functional needs **vertex + edge** DOFs, so the full -Newton solve is **not yet possible** in C++ — the *gradient* supports edge DOFs, -the *analytic Hessian* does not. (This is also why PR #27 cross-validated only -the cyclic *evaluation*, never a solve.) +blocker first: `newton_euclidean` → `euclidean_hessian` threw *"edge DOFs are +not supported"*. The cyclic functional needs **vertex + edge** DOFs (the +*gradient* supported them; the *analytic cotangent Hessian* did not). That is +why PR #27 cross-validated only the cyclic *evaluation*, never a solve. -**Landed in this PR instead** (`test_euclidean_functional.cpp`): -- ✅ `CyclicCircularEdge_PhiEntersGradient_CatHead` (**GREEN**) — solver-free - cross-validation that the circular-edge φ target enters the cyclic gradient - exactly (`ΔG_e = −Δφ_e` at 1e-12, no other component moves). This is the - evaluation-level prerequisite of the convergence oracle. -- ⏸️ `DISABLED_CyclicCircularEdge_CatHead_JavaXVal` — the full Java convergence - assertion (`α_opp+α_opp = π−0.1`), kept in-code with the golden semantics; it - auto-activates (drop `DISABLED_`) once the edge-DOF Hessian lands. +**Resolved (2026-05-30):** implemented a **block-FD edge-DOF Euclidean Hessian** +(`euclidean_hessian_block_fd` / `_sym` in `euclidean_hessian.hpp`, mirroring +`hyper_ideal_hessian_block_fd`); `newton_euclidean` now routes cyclic layouts +(edge DOFs present) through it, vertex-only layouts still use the analytic +cotangent Laplacian. -**New prerequisite for Tier-1 cyclic:** an **edge-DOF Euclidean Hessian** (the -cyclic Hessian over vertex + edge variables, cf. Java `getNonZeroPattern`). Spherical -convergence (Tier-1 #2) turned out to be **already covered** by -`test_newton_solver` (`Spherical_ConvergesFromPerturbation` et al. on the -spherical tetrahedron), and `make_octahedron_face()` is a single face, not a -closed octahedron — so Tier-1 #2 adds little. → the next genuinely-new -cross-validation target is **Tier 2 (Wente uniformization XML)** or implementing -the edge-DOF Hessian to unblock the cyclic convergence oracle. +**Landed (`test_euclidean_functional.cpp`, all GREEN, 243/243 cgal tests pass):** +- ✅ `CyclicCircularEdge_PhiEntersGradient_CatHead` — solver-free check that the + circular-edge φ enters the gradient exactly (`ΔG_e = −Δφ_e` @1e-12). +- ✅ `CyclicCircularEdge_CatHead_JavaXVal` — **the full Java convergence + oracle**: prescribe φ = π−0.1 on one interior ("circular") edge of cathead, + solve the cyclic Newton, assert the realised `α_opp+α_opp = π−0.1` @1e-9. +- ✅ `CyclicHessian_BlockFD_MatchesGradientFD_Tetrahedron` — the block-FD edge-DOF + Hessian matches the column-wise gradient FD on the full cyclic layout. + +> DOF note: only the **single** circular edge gets an edge DOF (as in Java's +> one `circularHoleEdge`). Giving *every* edge a DOF makes λ_e redundant with +> u_i+u_j (a V-dim null space) and stalls Newton. + +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. ### Other generators in the Java tree `HyperellipticCurveGenerator` (→ Phase 13), `SchottkyGenerator` (→ Phase 11a),