From ba83974525fbc7baa04bd0a2a81a2a797558a1e7 Mon Sep 17 00:00:00 2001 From: Tarik Moussa Date: Fri, 29 May 2026 19:08:37 +0200 Subject: [PATCH] test: Java golden-value oracles for the five DCE math cores + P1-2/P1-3 fixes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add bit-for-bit (1e-12) golden-value oracle tests pinning the C++ pure-math and functional cores against the compiled upstream Java library (openjdk 17): - HyperIdealGoldenJava: Clausen/Л/ImLi2, ζ13/14/15/ζ, both tetrahedron-volume formulas (real de.varylab…Clausen / HyperIdealUtility). - EuclideanGoldenJava / SphericalGoldenJava: angle formulas + β relations + Л energy terms, plus FULL-MESH oracles driving the real EuclideanCyclicFunctional / SphericalFunctional on a shared tetrahedron — per-vertex gradient (Θ−Σα) and ΔE = E(x)−E(0) (C++ Gauss-Legendre path integral vs Java closed form). - SphericalGoldenJava.FullMeshEdgeDofGradient: edge-DOF gradient (vertex + edge components, α_opp⁺+α_opp⁻−θ_e) vs raw conformalEnergyAndGradient — locks Finding 3 at the solution level (audit items 4 & 5). - PeriodMatrix.NormalizeModulus_GoldenJava: τ-reduction fold convention vs the real DiscreteEllipticUtility.normalizeModulus (audit items 7 & 8). Subtlety documented: the spherical oracles call Java's raw conformalEnergyAndGradient, not evaluate() (which pre-runs a Brent gauge maximization that C++ factors into the Newton solver's spherical_gauge_shift). Also: - P1-2 (layout.hpp): Euclidean holonomy now uses a per-cut-edge rigid-motion fit g(z)=a·z+b, exposing residual_rotation = |arg(a)| as a diagnostic; non- regressive (flat case a=1 reduces to the old midpoint formula). - P1-3 (period_matrix.hpp): is_in_fundamental_domain fixed to the correct half-open SL(2,ℤ) domain (−½ ≤ Re < ½). Updated the now-exposed ComputePeriodMatrix_ReducedTau_InFD to assert the normalizeModulus domain (closed +½ edge) instead. Test counts (single source of truth = doc/api/tests.md): 272/272 pass, 0 skipped (26 non-CGAL + 246 CGAL). Co-Authored-By: Claude Opus 4.8 --- CLAUDE.md | 12 +- code/include/layout.hpp | 55 ++-- code/include/period_matrix.hpp | 11 +- code/tests/cgal/test_euclidean_functional.cpp | 115 +++++++++ code/tests/cgal/test_phase7.cpp | 44 +++- code/tests/cgal/test_spherical_functional.cpp | 236 ++++++++++++++++++ code/tests/test_hyper_ideal_utility.cpp | 42 ++++ doc/api/tests.md | 10 +- doc/reviewer/java-port-audit.md | 116 ++++++--- 9 files changed, 580 insertions(+), 61 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 6fe13ae..b792932 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -282,7 +282,7 @@ Expected results: `test-fast` + `quality-gates` green on every push, 0 skipped, ## Release state Current release: **v0.10.0** (tag on `main`, released 2026-05-26 — the "reviewer-ready" release). -Phases 1–9a complete, Phase 8b-Lite CGAL API surface complete (all 5 DCE models reachable via ``), Phase 9b block-FD HyperIdeal Hessian shipped (~96× speed-up). v0.10.0 added: 100 % Doxygen public-API coverage (396/396 symbols, 0 warnings), a 14-gate structural quality suite (4 required in CI), the reviewer materials package (`doc/reviewer/`), `output_uv_map` for 4 of 5 DCE entries, six new roadmap phases + three RESEARCH phases, and a six-mode compile-time workflow matrix. Numbers at release: **259/259 tests pass, 0 skipped** (23 non-CGAL + 236 CGAL). Next planned milestones: Phase 9c (4g-polygon, genus g > 1) and Phase 9b-analytic (Schläfli identity). See `doc/release-policy.md` for the version-tag policy, `CHANGELOG.md` for full release notes, and `doc/roadmap/phases.md` for the phase plan. +Phases 1–9a complete, Phase 8b-Lite CGAL API surface complete (all 5 DCE models reachable via ``), Phase 9b block-FD HyperIdeal Hessian shipped (~96× speed-up). v0.10.0 added: 100 % Doxygen public-API coverage (396/396 symbols, 0 warnings), a 14-gate structural quality suite (4 required in CI), the reviewer materials package (`doc/reviewer/`), `output_uv_map` for 4 of 5 DCE entries, six new roadmap phases + three RESEARCH phases, and a six-mode compile-time workflow matrix. Numbers (single source of truth = `doc/api/tests.md`): **272/272 tests pass, 0 skipped** (26 non-CGAL + 246 CGAL). Next planned milestones: Phase 9c (4g-polygon, genus g > 1) and Phase 9b-analytic (Schläfli identity). See `doc/release-policy.md` for the version-tag policy, `CHANGELOG.md` for full release notes, and `doc/roadmap/phases.md` for the phase plan. ## Phase 8 CGAL-package decisions (frozen 2026-05-19) @@ -302,14 +302,18 @@ The 2026-05-21 audit corrected four mis-labels (now fixed): `InversiveDistanceFu ### Java↔C++ math-correctness audit (2026-05-29) -Full line-by-line audit of all math-critical headers against `de.varylab.discreteconformal.*` lives in **`doc/reviewer/java-port-audit.md`** (read it before re-investigating any of these). All 11 findings are resolved or noted; the four that needed code changes are now ✅ FIXED and **all CGAL tests pass** (240 after the Euclidean holonomy/τ end-to-end tests + the spherical edge-DOF closed-form oracle landed on top — see `doc/api/tests.md`): +Full line-by-line audit of all math-critical headers against `de.varylab.discreteconformal.*` lives in **`doc/reviewer/java-port-audit.md`** (read it before re-investigating any of these). All 11 findings are resolved or noted; the four that needed code changes are now ✅ FIXED and **all CGAL tests pass** (246 after the Euclidean holonomy/τ end-to-end tests, the spherical edge-DOF closed-form oracle, and the Java golden-value oracle suites — incl. three *full-mesh* oracles driving the real `EuclideanCyclicFunctional`/`SphericalFunctional` on a shared tetrahedron, one of them an edge-DOF gradient oracle — landed on top; see `doc/api/tests.md`): - **Finding 3** (`spherical_functional.hpp`) — spherical edge-DOF now uses Java's *replacement* parameterization (`Λ = λ_e` when the edge is a DOF, via helper `spher_eff_lambda`); edge gradient is `α_opp⁺ + α_opp⁻ − θ_e` (dropped the extra `−(S_f⁺+S_f⁻)/2` term). Vertex-only path is bit-for-bit unchanged. - **Finding 4** (`spherical_hessian.hpp`) — added an always-compiled `throw std::logic_error` edge-DOF guard (mirrors Finding 2 for Euclidean). - **Finding 6** (`period_matrix.hpp`) — `compute_period_matrix` now calls the faithful `normalizeModulus` (matches Java oracle: `0 ≤ Re ≤ ½`, `Im ≥ 0`, `|τ| ≥ 1`); `reduce_to_fundamental_domain` retained for the canonical SL(2,ℤ) domain. - **Finding 9** (`inversive_distance_functional.hpp`) — degenerate-face gradient now uses limiting angles instead of skipping (mirrors Finding 1); the genuinely-non-real `l*sq <= 0` skip is kept. -**Open follow-ups (tests only, not bugs):** the edge-DOF paths for Findings 3 & 4 are verified for compilation + the unchanged vertex-only case but **not yet against a Java oracle at the solution level** — missing-test items 4, 7, 10 in the audit doc (golden-value oracle tests, period-matrix sign-convention test, inversive-distance degenerate test). Build/test reminder: the CGAL build dir used for this audit is **`build-cgal`** at the *repo root* (not under `code/`), target `conformallab_cgal_tests`, filter `ctest -R '^cgal\.'`. +**Java golden-value oracles (2026-05-29, P0):** five suites now pin the C++ pure-math core bit-for-bit (1e-12) against the compiled Java library (openjdk 17, real `Clausen.Л` / `HyperIdealUtility` / `DiscreteEllipticUtility.normalizeModulus`): `HyperIdealGoldenJava` (Clausen/Л/ImLi₂, ζ₁₃/₁₄/₁₅/ζ, both tetrahedron-volume formulas), `EuclideanGoldenJava` (angle formula + 2·Л energy), `SphericalGoldenJava` (law-of-cosines angles + β relations + Л energy), and `PeriodMatrix.NormalizeModulus_GoldenJava` (audit missing-test item 7, ✅ done). Oracle harnesses live in `/tmp/oracle/*.java` (recipe in the audit doc). + +**Full-mesh oracles (2026-05-29):** `EuclideanGoldenJava.FullMeshGradientAndEnergy_Tetrahedron` and `SphericalGoldenJava.FullMeshGradientAndEnergy_Tetrahedron` drive the *real* `EuclideanCyclicFunctional` / `SphericalFunctional` on a shared tetrahedron (`/tmp/oracle/{tet.obj,EucMeshOracle.java,SphereMeshOracle.java}`) and pin both the per-vertex gradient (`Θ−Σα`) and `ΔE = E(x)−E(0)` to 1e-12 — the energy cross-checks C++'s Gauss-Legendre *path integral* against Java's *closed-form* functional (audit missing-test item 5, ✅ done). **Finding surfaced:** the spherical oracle must call Java's raw `conformalEnergyAndGradient`, NOT `evaluate()` — the latter pre-runs a 1-D Brent maximization over the global-scale gauge (`maximizeInNegativeDirection`), which C++ deliberately factors into the Newton solver's `spherical_gauge_shift` instead (both correct; different factoring). + +**Open follow-ups (tests only, not bugs):** (a) the spherical edge-DOF *gradient* is now Java-oracle'd at the solution level (`SphericalGoldenJava.FullMeshEdgeDofGradient_Tetrahedron` — vertex + edge components vs raw `conformalEnergyAndGradient`, locking Finding 3); the only remaining edge-DOF gap is *Newton-to-convergence* with edge DOFs, blocked by the Finding-4 spherical-Hessian guard (needs an FD-Hessian or guard relaxation first) — a solver feature, not a correctness gap; (b) inversive-distance degenerate-face robustness (item 10) has no Java oracle because the functional is research with no Java parent — only a NaN-free limiting-angle regression test applies. Build/test reminder: the CGAL build dir used for this audit is **`build-cgal`** at the *repo root* (not under `code/`), target `conformallab_cgal_tests`, filter `ctest -R '^cgal\.'`. ## Documentation map @@ -422,7 +426,7 @@ Every turn re-sends the whole conversation, so context accumulation is the large ### Token hygiene — command discipline (Tier 2) -- **Scope + filter together.** Path-scope searches (`grep -rn "newton_" code/include/`, not repo-wide). Filter test output (`ctest -R "cgal.NewtonSolver" --output-on-failure`, or `--gtest_filter` for one suite) instead of dumping all 259 results. +- **Scope + filter together.** Path-scope searches (`grep -rn "newton_" code/include/`, not repo-wide). Filter test output (`ctest -R "cgal.NewtonSolver" --output-on-failure`, or `--gtest_filter` for one suite) instead of dumping all 272 results. - **Never read raw logs into context.** Redirect to a file, then `grep`/`tail` it. With `run_in_background`, read the output file selectively rather than pulling it whole. - **Don't re-read a file you just edited.** `Edit` errors on stale state — the harness tracks it; a `Read`-back after a successful edit is wasted context. - **Reference by line number** (`newton_solver.hpp:147`) instead of re-pasting code blocks. diff --git a/code/include/layout.hpp b/code/include/layout.hpp index 2d8f9ed..6a52f08 100644 --- a/code/include/layout.hpp +++ b/code/include/layout.hpp @@ -176,6 +176,14 @@ struct HolonomyData { std::vector translations; ///< Euclidean / spherical translation per cut edge. std::vector mobius_maps; ///< Hyperbolic Möbius isometry per cut edge (Phase 7). std::vector cut_edge_indices; ///< Index (in the cut-graph edge list) of each holonomy entry. + + /// Residual rotation |arg(a)| (radians) of the Euclidean deck isometry + /// z ↦ a·z + b per cut edge. Zero for a perfectly flat cone metric; + /// a non-zero value signals under-convergence or a genuine cone-angle + /// defect, in which case `translations[i]` (the isometry's translation + /// part b) is only an approximate lattice generator. Empty for the + /// hyperbolic (Möbius) path. + std::vector residual_rotation; }; // ── Internal helpers ────────────────────────────────────────────────────────── @@ -490,11 +498,23 @@ inline void set_root_huv_2d( // different locations. The rigid motion identifying those two copies is the // deck transformation of the generator loop (cut edge + tree path) — i.e. the // holonomy. For a flat cone metric (Θ ≡ 2π) the linear part is trivial, so the -// holonomy is the pure translation given by the displacement of the shared -// edge's midpoint between the two developments. That translation is exactly -// the lattice generator ω consumed by compute_period_matrix(). +// holonomy is the pure translation that identifies the two developments of the +// shared edge. We recover it as a full rigid motion z ↦ a·z + b fitted to the +// correspondence (Bs,Bt) ↦ (As,At): the shared edge endpoints S,T as developed +// in face B map to the same endpoints as developed in face A. Because both +// faces use the same edge length, |a| = 1 and the fit is an exact +// orientation-preserving isometry. Its translation part b is the lattice +// generator ω consumed by compute_period_matrix(); the rotation magnitude +// |arg(a)| is returned as a convergence diagnostic. For a perfectly flat cone +// metric (Θ ≡ 2π) a = 1 and b reduces to the midpoint displacement midA − midB +// (the previous formula), so converged flat tori are bit-for-bit unchanged. +struct EuclideanHolonomyResult { + std::vector translations; ///< ω_i = translation part b + std::vector residual_rotation; ///< |arg(a_i)| per cut edge +}; + template -inline std::vector euclidean_holonomy( +inline EuclideanHolonomyResult euclidean_holonomy( const ConformalMesh& mesh, const CutGraph& cut, EdgeLenFn&& edge_len) @@ -565,9 +585,10 @@ inline std::vector euclidean_holonomy( for (auto f : mesh.faces()) if (!face_done[static_cast(f.idx())]) develop(f); - // ── Holonomy translation per cut edge ───────────────────────────────────── - std::vector omega; - omega.reserve(cut.cut_edge_indices.size()); + // ── Holonomy isometry per cut edge ───────────────────────────────────────── + EuclideanHolonomyResult out; + out.translations.reserve(cut.cut_edge_indices.size()); + out.residual_rotation.reserve(cut.cut_edge_indices.size()); for (std::size_t ce : cut.cut_edge_indices) { Edge_index e = *std::next(mesh.edges().begin(), static_cast(ce)); Halfedge_index h = mesh.halfedge(e); @@ -575,7 +596,8 @@ inline std::vector euclidean_holonomy( if (mesh.is_border(h) || mesh.is_border(ho) || !face_done[static_cast(mesh.face(h).idx())] || !face_done[static_cast(mesh.face(ho).idx())]) { - omega.push_back(Eigen::Vector2d::Zero()); + out.translations.push_back(Eigen::Vector2d::Zero()); + out.residual_rotation.push_back(0.0); continue; } // Face A = face(h) places edge e as halfedge h: S=source(h), T=target(h) @@ -585,12 +607,15 @@ inline std::vector euclidean_holonomy( // target(ho)=S → S=pos of source(next(ho)), T=pos of source(ho) C Bt = hpos[static_cast(ho.idx())]; C Bs = hpos[static_cast(mesh.next(ho).idx())]; - C midA = 0.5 * (As + At); - C midB = 0.5 * (Bs + Bt); - C w = midA - midB; // pure translation for a flat cone metric - omega.push_back(Eigen::Vector2d(w.real(), w.imag())); + // Fit the deck isometry g(z) = a·z + b with (Bs,Bt) ↦ (As,At). + // |a| = 1 exactly (shared edge length); b is the translation part. + C edgeB = Bt - Bs; + C a = (std::abs(edgeB) > 1e-300) ? (At - As) / edgeB : C(1.0, 0.0); + C b = As - a * Bs; + out.translations.push_back(Eigen::Vector2d(b.real(), b.imag())); + out.residual_rotation.push_back(std::abs(std::arg(a))); } - return omega; + return out; } } // namespace detail @@ -719,7 +744,9 @@ inline Layout2D euclidean_layout( if (cut && holonomy) { holonomy->cut_edge_indices = cut->cut_edge_indices; holonomy->mobius_maps.clear(); - holonomy->translations = detail::euclidean_holonomy(mesh, *cut, edge_len); + auto eh = detail::euclidean_holonomy(mesh, *cut, edge_len); + holonomy->translations = std::move(eh.translations); + holonomy->residual_rotation = std::move(eh.residual_rotation); // Preserve the per-cut-edge seam UV in halfedge_uv (texture atlas), as // before, so HalfedgeUV-based tests still see the seam-crossing layout. diff --git a/code/include/period_matrix.hpp b/code/include/period_matrix.hpp index 0b081dd..921e5ba 100644 --- a/code/include/period_matrix.hpp +++ b/code/include/period_matrix.hpp @@ -114,11 +114,18 @@ inline std::complex reduce_to_fundamental_domain(std::complex ta // is_in_fundamental_domain — check membership in F with tolerance tol. // ───────────────────────────────────────────────────────────────────────────── /// `true` iff `τ` lies inside the standard SL(2,ℤ) fundamental domain -/// with tolerance `tol`. +/// `F = { Im τ > 0, −½ ≤ Re τ < ½, |τ| ≥ 1 }` with tolerance `tol`. +/// +/// This is the half-open domain produced by `reduce_to_fundamental_domain` +/// (the right boundary `Re τ = +½ ≡ −½` is excluded via `T`). It is NOT the +/// mirror-folded `0 ≤ Re τ ≤ ½` domain produced by `normalizeModulus` (used +/// inside `compute_period_matrix`); a τ with `Re τ < 0` is a legitimate member +/// of `F` here but would be folded to `Re ≥ 0` by `normalizeModulus`. inline bool is_in_fundamental_domain(std::complex tau, double tol = 1e-9) { if (tau.imag() <= 0.0) return false; - if (std::abs(tau.real()) > 0.5 + tol) return false; + if (tau.real() < -0.5 - tol) return false; // left boundary closed + if (tau.real() > 0.5 - tol) return false; // right boundary Re = +½ excluded (≡ −½) if (std::abs(tau) < 1.0 - tol) return false; return true; } diff --git a/code/tests/cgal/test_euclidean_functional.cpp b/code/tests/cgal/test_euclidean_functional.cpp index 2046262..5f9e249 100644 --- a/code/tests/cgal/test_euclidean_functional.cpp +++ b/code/tests/cgal/test_euclidean_functional.cpp @@ -26,6 +26,7 @@ #include "euclidean_geometry.hpp" #include "euclidean_functional.hpp" #include "euclidean_hessian.hpp" +#include "clausen.hpp" #include #include #include @@ -290,3 +291,117 @@ TEST(EuclideanFunctional, GradientCheck_MixedPinnedVertices) EXPECT_TRUE(gradient_check_euclidean(mesh, x, maps)) << "Gradient check failed for mixed pinned/variable vertices"; } + +// ───────────────────────────────────────────────────────────────────────────── +// Golden-value oracle — pin the Euclidean angle formula and the 2·Л(α) energy +// term bit-for-bit against the upstream Java reference (EuclideanCyclicFunctional +// .triangleEnergyAndAlphas, lines 341-361), captured by running the compiled Java +// library (openjdk 17) with the real de.varylab…Clausen.Л on these exact edge +// lengths. Companion to HyperIdealGoldenJava and SphericalGoldenJava: locks the +// absolute angle/energy values against an independent implementation, catching +// silent index/sign drift the curl-free path-integral gradient check cannot see. +// +// To regenerate: /tmp/oracle/EucOracle.java. Values are Java printf %.17g. +// ───────────────────────────────────────────────────────────────────────────── +TEST(EuclideanGoldenJava, AngleAndLobachevskyEnergyFromLengths) +{ + auto check = [](double l12, double l23, double l31, + double a1_g, double a2_g, double a3_g, double L_g) { + auto fa = euclidean_angles_from_lengths(l12, l23, l31); + EXPECT_TRUE(fa.valid); + EXPECT_NEAR(fa.alpha1, a1_g, 1e-12); + EXPECT_NEAR(fa.alpha2, a2_g, 1e-12); + EXPECT_NEAR(fa.alpha3, a3_g, 1e-12); + const double Lterm = 2.0 * Lobachevsky(fa.alpha1) + + 2.0 * Lobachevsky(fa.alpha2) + + 2.0 * Lobachevsky(fa.alpha3); + EXPECT_NEAR(Lterm, L_g, 1e-12); + }; + check(1.0, 1.2, 0.9, + 1.3637649752769678, 0.82416964552030680, 0.95365803279251860, + 1.9456273836230942); + check(1.0, 1.0, 1.0, + 1.0471975511965979, 1.0471975511965979, 1.0471975511965979, + 2.0298832128193070); +} + +// ───────────────────────────────────────────────────────────────────────────── +// FULL-MESH golden oracle — the strongest cross-check: drives the REAL upstream +// EuclideanCyclicFunctional (openjdk 17) on a tetrahedron loaded from a shared +// OBJ (identical topology + geometry to make_tetrahedron()), and pins BOTH the +// per-vertex gradient G_v = Θ_v − Σα AND the energy difference ΔE = E(x) − E(0) +// bit-for-bit against it. +// +// This closes audit missing-test item 5 (full-mesh energy + gradient at a known +// x). It is genuinely independent of the C++ implementation in two ways: +// • the gradient is the upstream library's own analytic gradient (not an FD +// check, which only proves curl-freeness of the C++ self-consistent energy); +// • the energy is Java's CLOSED-FORM functional value, whereas C++ computes it +// as a Gauss-Legendre PATH INTEGRAL of its gradient — two different methods +// that must agree on ΔE (the initialEnergy constant and the φ·λ⁰ term cancel +// in the difference, and there are no edge DOFs here). +// +// Setup parity (verified against UnwrapUtility.prepareInvariantDataEuclidean): +// closed mesh, ALL 4 vertices variable (no pin), Θ_v = 2π, no edge DOFs, +// λ°_e = 2·log(|p_i − p_j|), per-vertex u(P) = 0.10·X − 0.07·Y + 0.13·Z. +// +// To regenerate: /tmp/oracle/{tet.obj,EucMeshOracle.java}. Values are Java %.17g. +// ───────────────────────────────────────────────────────────────────────────── +TEST(EuclideanGoldenJava, FullMeshGradientAndEnergy_Tetrahedron) +{ + constexpr double TWO_PI = 2.0 * 3.14159265358979323846264338328; + + auto mesh = make_tetrahedron(); // same 4 vertices as /tmp/oracle/tet.obj + auto maps = setup_euclidean_maps(mesh); + compute_euclidean_lambda0_from_mesh(mesh, maps); + + // All four vertices are DOFs (Java prepareInvariantData makes every interior + // vertex variable on a closed mesh); Θ_v = 2π; no edge DOFs. + int idx = 0; + for (auto v : mesh.vertices()) { + maps.v_idx[v] = idx++; + maps.theta_v[v] = TWO_PI; + } + + auto u_of = [](const Point3& p) { + return 0.10 * p.x() - 0.07 * p.y() + 0.13 * p.z(); + }; + + std::vector x(static_cast(idx), 0.0); + for (auto v : mesh.vertices()) + x[static_cast(maps.v_idx[v])] = u_of(mesh.point(v)); + + auto G = euclidean_gradient(mesh, x, maps); + + // Java golden gradients keyed by vertex position (order-independent lookup). + struct GoldRow { double X, Y, Z, G; }; + const GoldRow gold[4] = { + { 1, 1, 1, 3.5277511803984396}, + { 1, -1, -1, 3.2837611905358440}, + {-1, 1, -1, 2.3437485291237100}, + {-1, -1, 1, 3.4111097143011780}, + }; + for (auto v : mesh.vertices()) { + const auto& p = mesh.point(v); + const double g = G[static_cast(maps.v_idx[v])]; + bool matched = false; + for (const auto& row : gold) { + if (std::abs(p.x() - row.X) < 1e-9 && + std::abs(p.y() - row.Y) < 1e-9 && + std::abs(p.z() - row.Z) < 1e-9) { + EXPECT_NEAR(g, row.G, 1e-12) + << "gradient mismatch at (" << p.x() << "," << p.y() + << "," << p.z() << ")"; + matched = true; + break; + } + } + EXPECT_TRUE(matched) << "unexpected vertex position"; + } + + // ΔE = E(x) − E(0). C++ path integral vs Java closed-form functional value. + std::vector x0(static_cast(idx), 0.0); + const double dE = euclidean_energy(mesh, x, maps) + - euclidean_energy(mesh, x0, maps); + EXPECT_NEAR(dE, 0.15962619236187336, 1e-12); +} diff --git a/code/tests/cgal/test_phase7.cpp b/code/tests/cgal/test_phase7.cpp index a870799..06c9647 100644 --- a/code/tests/cgal/test_phase7.cpp +++ b/code/tests/cgal/test_phase7.cpp @@ -338,12 +338,52 @@ TEST(PeriodMatrix, ComputePeriodMatrix_UnitSquare) TEST(PeriodMatrix, ComputePeriodMatrix_ReducedTau_InFD) { // ω_1 = (1, 0), ω_2 = (0.5, 0.25) → τ = 0.5 + 0.25i - // |τ| = sqrt(0.25 + 0.0625) ≈ 0.559 < 1 → needs S step + // |τ| = sqrt(0.25 + 0.0625) ≈ 0.559 < 1 → needs S step. + // compute_period_matrix reduces with normalizeModulus (Finding 6), whose + // mirror-folded target domain is { 0 ≤ Re ≤ ½, Im > 0, |τ| ≥ 1 } — note the + // RIGHT boundary Re = +½ is CLOSED here. This is NOT the half-open SL(2,ℤ) + // domain of is_in_fundamental_domain (−½ ≤ Re < ½), which excludes Re = +½ + // because +½ ≡ −½ under T. For this input normalizeModulus lands exactly on + // τ = ½ + i, so we must check the normalizeModulus domain, not the SL(2,ℤ) + // one (asserting is_in_fundamental_domain here would wrongly fail on +½). HolonomyData hol; hol.translations = { Eigen::Vector2d(1.0, 0.0), Eigen::Vector2d(0.5, 0.25) }; PeriodData pd = compute_period_matrix(hol, /*reduce=*/true); EXPECT_TRUE(pd.in_fundamental_domain); - EXPECT_TRUE(is_in_fundamental_domain(pd.tau, 1e-9)); + const double tol = 1e-9; + EXPECT_GT(pd.tau.imag(), 0.0); + EXPECT_GE(pd.tau.real(), 0.0 - tol); // 0 ≤ Re (mirror fold) + EXPECT_LE(pd.tau.real(), 0.5 + tol); // Re ≤ ½ (closed right edge) + EXPECT_GE(std::abs(pd.tau), 1.0 - tol); // |τ| ≥ 1 + // Concretely: τ = ½ + i. + EXPECT_NEAR(pd.tau.real(), 0.5, 1e-12); + EXPECT_NEAR(pd.tau.imag(), 1.0, 1e-12); +} + +// ───────────────────────────────────────────────────────────────────────────── +// Golden-value oracle — pin normalizeModulus (the τ reduction used by +// compute_period_matrix, Finding 6) bit-for-bit against the upstream Java +// reference (de.varylab.discreteconformal.util.DiscreteEllipticUtility. +// normalizeModulus), captured by calling the compiled Java method (openjdk 17) +// on these exact τ. This locks the sign/fold conventions of the SL(2,ℤ)+mirror +// reduction (0 ≤ Re ≤ ½, Im ≥ 0, |τ| ≥ 1) against an independent implementation, +// catching drift the existing in-FD membership checks cannot (they only assert +// the result lies in F, not that it is the SAME representative Java picks). +// +// To regenerate: /tmp/oracle/TauOracle.java. Values are Java printf %.17g. +// ───────────────────────────────────────────────────────────────────────────── +TEST(PeriodMatrix, NormalizeModulus_GoldenJava) +{ + auto chk = [](double re, double im, double re_g, double im_g) { + C n = conformallab::normalizeModulus(C(re, im)); + EXPECT_NEAR(n.real(), re_g, 1e-12); + EXPECT_NEAR(n.imag(), im_g, 1e-12); + }; + chk(0.3, 0.5, 0.11764705882352933, 1.4705882352941178); // |τ|<1 → S + folds + chk(-0.4, 1.3, 0.40000000000000000, 1.3000000000000000); // Re<0 mirror fold + chk(2.7, 0.8, 0.41095890410958880, 1.0958904109589043); // large Re → T + chk(0.1, 2.0, 0.10000000000000000, 2.0000000000000000); // already in F + chk(-1.6, 0.9, 0.41237113402061850, 0.92783505154639180); // T + S + mirror } // ════════════════════════════════════════════════════════════════════════════ diff --git a/code/tests/cgal/test_spherical_functional.cpp b/code/tests/cgal/test_spherical_functional.cpp index dc3edc3..a8d745f 100644 --- a/code/tests/cgal/test_spherical_functional.cpp +++ b/code/tests/cgal/test_spherical_functional.cpp @@ -27,6 +27,8 @@ #include "mesh_builder.hpp" #include "spherical_functional.hpp" #include "spherical_hessian.hpp" +#include "spherical_geometry.hpp" +#include "clausen.hpp" #include #include #include @@ -416,3 +418,237 @@ TEST(SphericalFunctional, GaugeFix_AlreadyAtGaugeReturnsTNearZero) EXPECT_NEAR(t, 0.0, 1e-5) << "Gauge shift from the symmetric point should be ~0; got " << t; } + +// ───────────────────────────────────────────────────────────────────────────── +// Golden-value oracle — pin the spherical law-of-cosines angle formula, the β +// half-angle relations, and the Lobachevsky energy term bit-for-bit against the +// upstream Java reference (SphericalFunctional.triangleEnergyAndAlphas, lines +// 401-453), captured by running the compiled Java library (openjdk 17) with the +// real de.varylab…Clausen.Л on these exact arc lengths. Unlike the Schläfli +// path-integral gradient check (which only verifies curl-freeness), this locks +// the absolute angle/β/energy values against an independent implementation, +// catching any silent index/sign/convention drift in the spherical port. +// +// To regenerate: /tmp/oracle/SphereOracle.java (recipe in doc/reviewer/ +// java-port-audit.md). Values are Java printf %.17g. Tolerance 1e-12. +// +// Index map (C++ spherical_angles ↔ Java): with l12=lij, l23=ljk, l31=lki the +// returned alpha1/alpha2/alpha3 are Java αi/αj/αk (angle opposite ljk/lki/lij). +// ───────────────────────────────────────────────────────────────────────────── +TEST(SphericalGoldenJava, AngleBetaEnergyFromLengths) +{ + auto check = [](double lij, double ljk, double lki, + double ai_g, double aj_g, double ak_g, + double bi_g, double bj_g, double bk_g, double L_g) { + auto fa = spherical_angles(lij, ljk, lki); + EXPECT_TRUE(fa.valid); + EXPECT_NEAR(fa.alpha1, ai_g, 1e-12); + EXPECT_NEAR(fa.alpha2, aj_g, 1e-12); + EXPECT_NEAR(fa.alpha3, ak_g, 1e-12); + + const double ai = fa.alpha1, aj = fa.alpha2, ak = fa.alpha3; + const double bi = 0.5 * (PI + ai - aj - ak); + const double bj = 0.5 * (PI - ai + aj - ak); + const double bk = 0.5 * (PI - ai - aj + ak); + EXPECT_NEAR(bi, bi_g, 1e-12); + EXPECT_NEAR(bj, bj_g, 1e-12); + EXPECT_NEAR(bk, bk_g, 1e-12); + + const double Lterm = + Lobachevsky(ai) + Lobachevsky(aj) + Lobachevsky(ak) + + Lobachevsky(bi) + Lobachevsky(bj) + Lobachevsky(bk) + + Lobachevsky(0.5 * (PI - ai - aj - ak)); + EXPECT_NEAR(Lterm, L_g, 1e-12); + }; + + // Scalene spherical triangle (all Δ > 0, Δ_ijk < 2π). + check(1.0, 1.2, 0.9, + 1.5305813141072122, 0.99684981272794600, 1.1246069519101605, + 1.2753586015294494, 0.74162710015018330, 0.86938423933239760, + 1.3576100185550408); + // Equilateral spherical triangle. + check(0.7, 0.7, 0.7, + 1.1225596283199812, 1.1225596283199812, 1.1225596283199812, + 1.0095165126349062, 1.0095165126349060, 1.0095165126349060, + 1.6806828584976297); +} + +// ───────────────────────────────────────────────────────────────────────────── +// FULL-MESH golden oracle (spherical) — drives the REAL upstream SphericalFunctional +// on the shared tetrahedron (scaled so every arc length stays < π) and pins BOTH +// the per-vertex gradient G_v = Θ_v − Σα AND ΔE = E(x) − E(0) bit-for-bit. +// +// IMPORTANT — the oracle calls Java's `conformalEnergyAndGradient` (the RAW +// θ−Σα gradient), NOT `evaluate()`: `evaluate()` first runs a 1-D Brent +// maximization over the global-scale gauge direction (`maximizeInNegativeDirection`), +// which the C++ `spherical_gradient` deliberately does NOT — C++ factors that +// gauge into the Newton solver's `spherical_gauge_shift` instead. Comparing +// against `evaluate()` would (wrongly) drive every component to ~0. The raw +// gradient is the piece that corresponds 1:1 to the C++ functional. +// +// Setup parity (UnwrapUtility.prepareInvariantDataHyperbolicAndSpherical, scale): +// closed mesh, ALL 4 vertices variable, Θ_v = 2π, no edge DOFs, +// λ°_e = 2·log(SCALE·|p_i − p_j|) (Java uses chord length × scale, whereas the +// C++ compute_lambda0_from_mesh helper assumes unit-sphere vertices and uses +// ARC length — so we set λ° directly here to match Java exactly), +// per-vertex u(P) = 0.10·X − 0.07·Y + 0.13·Z, SCALE = 0.2. +// +// To regenerate: /tmp/oracle/{tet.obj,SphereMeshOracle.java}. Values are Java %.17g. +// ───────────────────────────────────────────────────────────────────────────── +TEST(SphericalGoldenJava, FullMeshGradientAndEnergy_Tetrahedron) +{ + constexpr double TWO_PI = 2.0 * 3.14159265358979323846264338328; + constexpr double SCALE = 0.2; + + auto mesh = make_tetrahedron(); // same 4 vertices as /tmp/oracle/tet.obj + auto maps = setup_spherical_maps(mesh); + + // λ°_e = 2·log(SCALE · chord), matching Java's prepareInvariantData(scale). + for (auto e : mesh.edges()) { + auto h = mesh.halfedge(e); + const auto& p1 = mesh.point(mesh.source(h)); + const auto& p2 = mesh.point(mesh.target(h)); + const double dx = p1.x() - p2.x(), dy = p1.y() - p2.y(), dz = p1.z() - p2.z(); + const double chord = std::sqrt(dx*dx + dy*dy + dz*dz); + maps.lambda0[e] = 2.0 * std::log(SCALE * chord); + } + + int idx = 0; + for (auto v : mesh.vertices()) { + maps.v_idx[v] = idx++; + maps.theta_v[v] = TWO_PI; + } + + auto u_of = [](const Point3& p) { + return 0.10 * p.x() - 0.07 * p.y() + 0.13 * p.z(); + }; + + std::vector x(static_cast(idx), 0.0); + for (auto v : mesh.vertices()) + x[static_cast(maps.v_idx[v])] = u_of(mesh.point(v)); + + auto G = spherical_gradient(mesh, x, maps); + + struct GoldRow { double X, Y, Z, G; }; + const GoldRow gold[4] = { + { 1, 1, 1, 2.7671034786104927}, + { 1, -1, -1, 2.5216834054857546}, + {-1, 1, -1, 1.5401761310866633}, + {-1, -1, 1, 2.6518569531861100}, + }; + for (auto v : mesh.vertices()) { + const auto& p = mesh.point(v); + const double g = G[static_cast(maps.v_idx[v])]; + bool matched = false; + for (const auto& row : gold) { + if (std::abs(p.x() - row.X) < 1e-9 && + std::abs(p.y() - row.Y) < 1e-9 && + std::abs(p.z() - row.Z) < 1e-9) { + EXPECT_NEAR(g, row.G, 1e-12) + << "gradient mismatch at (" << p.x() << "," << p.y() + << "," << p.z() << ")"; + matched = true; + break; + } + } + EXPECT_TRUE(matched) << "unexpected vertex position"; + } + + std::vector x0(static_cast(idx), 0.0); + const double dE = spherical_energy(mesh, x, maps) + - spherical_energy(mesh, x0, maps); + EXPECT_NEAR(dE, 0.16409141487397116, 1e-12); +} + +// ───────────────────────────────────────────────────────────────────────────── +// FULL-MESH EDGE-DOF golden oracle (spherical, Finding 3) — the missing solution- +// level Java oracle for the edge-DOF gradient path. Makes two opposite edges of +// the shared tetrahedron variable (replacement parameterization: λ_e = x[e_idx], +// no u-shift) and pins BOTH the vertex gradient (Θ−Σα) AND the edge gradient +// (α_opp⁺ + α_opp⁻ − θ_e, θ_e = π) bit-for-bit against the REAL upstream +// SphericalFunctional.conformalEnergyAndGradient (raw gradient, not evaluate()). +// +// This is a pure GRADIENT oracle (no ΔE): with an edge DOF, x = 0 means λ_e = 0 +// ⇒ spherical arc length = π (degenerate), so the path-integral-from-origin +// energy reference is ill-defined — the gradient at a fixed non-degenerate x is +// the unambiguous Finding-3 quantity. It does NOT touch the spherical Hessian, +// so it is unaffected by the Finding-4 edge-DOF Hessian guard. +// +// To regenerate: /tmp/oracle/SphereEdgeOracle.java. Values are Java %.17g. +// ───────────────────────────────────────────────────────────────────────────── +TEST(SphericalGoldenJava, FullMeshEdgeDofGradient_Tetrahedron) +{ + constexpr double TWO_PI = 2.0 * 3.14159265358979323846264338328; + constexpr double PI_C = 3.14159265358979323846264338328; + constexpr double SCALE = 0.2; + + auto mesh = make_tetrahedron(); + auto maps = setup_spherical_maps(mesh); + + for (auto e : mesh.edges()) { + auto h = mesh.halfedge(e); + const auto& p1 = mesh.point(mesh.source(h)); + const auto& p2 = mesh.point(mesh.target(h)); + const double dx = p1.x()-p2.x(), dy = p1.y()-p2.y(), dz = p1.z()-p2.z(); + maps.lambda0[e] = 2.0 * std::log(SCALE * std::sqrt(dx*dx + dy*dy + dz*dz)); + maps.theta_e[e] = PI_C; + maps.e_idx[e] = -1; + } + + int idx = 0; + for (auto v : mesh.vertices()) { maps.v_idx[v] = idx++; maps.theta_v[v] = TWO_PI; } + + auto pos = [&](Vertex_index v){ return mesh.point(v); }; + auto same = [](const Point3& p, double X, double Y, double Z) { + return std::abs(p.x()-X)<1e-9 && std::abs(p.y()-Y)<1e-9 && std::abs(p.z()-Z)<1e-9; + }; + // Edge connects positions (aX,aY,aZ)–(bX,bY,bZ) in either order? + auto connects = [&](Edge_index e, const double a[3], const double b[3]) { + auto h = mesh.halfedge(e); + const Point3& s = pos(mesh.source(h)); + const Point3& t = pos(mesh.target(h)); + return (same(s,a[0],a[1],a[2]) && same(t,b[0],b[1],b[2])) || + (same(s,b[0],b[1],b[2]) && same(t,a[0],a[1],a[2])); + }; + const double A[3]={1,1,1}, B[3]={1,-1,-1}, C[3]={-1,1,-1}, D[3]={-1,-1,1}; + + // Make edges A–B and C–D variable (replacement parameterization). + int edof = idx; // edge DOFs start after the 4 vertex DOFs + Edge_index eAB, eCD; + for (auto e : mesh.edges()) { + if (connects(e, A, B)) { maps.e_idx[e] = edof++; eAB = e; } + else if (connects(e, C, D)) { maps.e_idx[e] = edof++; eCD = e; } + } + ASSERT_EQ(edof, 6); // 4 vertex + 2 edge DOFs + + auto u_of = [](const Point3& p){ return 0.10*p.x() - 0.07*p.y() + 0.13*p.z(); }; + + std::vector x(6, 0.0); + for (auto v : mesh.vertices()) + x[static_cast(maps.v_idx[v])] = u_of(mesh.point(v)); + // Edge DOF value = λ⁰_e + 0.1 (same perturbation as the Java oracle). + x[static_cast(maps.e_idx[eAB])] = maps.lambda0[eAB] + 0.1; + x[static_cast(maps.e_idx[eCD])] = maps.lambda0[eCD] + 0.1; + + auto G = spherical_gradient(mesh, x, maps); + + // Vertex gradients keyed by position. + struct VRow { double X, Y, Z, G; }; + const VRow vg[4] = { + { 1, 1, 1, 2.5036751008175546}, + { 1, -1, -1, 2.2303362601050205}, + {-1, 1, -1, 1.8463710314817632}, + {-1, -1, 1, 2.7499719487341570}, + }; + for (auto v : mesh.vertices()) { + const auto& p = mesh.point(v); + const double g = G[static_cast(maps.v_idx[v])]; + bool m = false; + for (auto& r : vg) + if (same(p, r.X, r.Y, r.Z)) { EXPECT_NEAR(g, r.G, 1e-12); m = true; break; } + EXPECT_TRUE(m); + } + // Edge gradients: A–B and C–D (Java golden values). + EXPECT_NEAR(G[static_cast(maps.e_idx[eAB])], -0.35189517043413690, 1e-12); + EXPECT_NEAR(G[static_cast(maps.e_idx[eCD])], -0.44101986058895950, 1e-12); +} diff --git a/code/tests/test_hyper_ideal_utility.cpp b/code/tests/test_hyper_ideal_utility.cpp index 8c68577..7688f5c 100644 --- a/code/tests/test_hyper_ideal_utility.cpp +++ b/code/tests/test_hyper_ideal_utility.cpp @@ -4,10 +4,12 @@ // Port of de.varylab.discreteconformal.functional.HyperIdealUtilityTest (Java/JUnit). #include "hyper_ideal_utility.hpp" +#include "hyper_ideal_geometry.hpp" #include "clausen.hpp" #include #include +#include using conformallab::calculateTetrahedronVolume; using conformallab::calculateTetrahedronVolumeWithIdealVertexAtGamma; @@ -93,3 +95,43 @@ TEST(HyperIdealUtilityTest, CompareGeneralAndIdealFormulaCase2) { double V = calculateTetrahedronVolume(bi, bj, bk, ai, aj, ak); EXPECT_NEAR(Ve, V, 1e-12); } + +// ───────────────────────────────────────────────────────────────────────────── +// Golden-value oracle tests — pin the C++ pure-math core bit-for-bit against the +// upstream Java reference (de.varylab.discreteconformal.functional.{Clausen, +// HyperIdealUtility}), captured by running the compiled Java library (openjdk 17) +// on these exact inputs. Unlike the FD gradient checks (which only verify +// curl-freeness, since the C++ energy is the path-integral of its own gradient), +// these lock the absolute values of the math-critical helpers against an +// independent implementation, catching any silent convention/formula drift. +// +// To regenerate: see doc/reviewer/java-port-audit.md (oracle harness recipe). +// Values are Java's Double.toString output (shortest round-trip). Tolerance is +// 1e-12 (well above the ~1e-15 inter-platform libm divergence for these ranges). +// ───────────────────────────────────────────────────────────────────────────── +TEST(HyperIdealGoldenJava, ClausenLobachevskyImLi2) { + EXPECT_NEAR(conformallab::clausen2(0.7), 0.954448086482735, 1e-12); + EXPECT_NEAR(conformallab::clausen2(2.5), 0.4335982032355327, 1e-12); + EXPECT_NEAR(conformallab::clausen2(-1.3), -0.9897032532295984, 1e-12); + EXPECT_NEAR(Lobachevsky(0.9), 0.4121734067662043, 1e-12); + EXPECT_NEAR(conformallab::ImLi2(std::complex(0.3, 0.4)), + 0.46136289181910894, 1e-12); +} + +TEST(HyperIdealGoldenJava, ZetaFamily) { + EXPECT_NEAR(conformallab::zeta13(0.5, 0.7, 0.9), 2.663195966482385, 1e-12); + EXPECT_NEAR(conformallab::zeta14(0.4, 0.8), 1.8262295633065202, 1e-12); + EXPECT_NEAR(conformallab::zeta15(0.6), 2.216976794676588, 1e-12); + EXPECT_NEAR(conformallab::zeta(0.5, 0.7, 0.9), 1.6156519307269948, 1e-12); +} + +TEST(HyperIdealGoldenJava, TetrahedronVolumes) { + // Generic non-degenerate generalized hyperbolic tetrahedron. + EXPECT_NEAR(calculateTetrahedronVolume(0.6, 0.7, 0.8, 0.5, 0.9, 0.4), + 2.9074633382516435, 1e-12); + // One ideal vertex at gamma, non-symmetric angles (Java test2 inputs). + EXPECT_NEAR(calculateTetrahedronVolumeWithIdealVertexAtGamma( + 0.6623267054958116, 1.437248992086214, 1.0420169560077686, + 0.6896178197389236, 0.5195634857410114, 0.6304500578493993), + 2.0459750326926214, 1e-12); +} diff --git a/doc/api/tests.md b/doc/api/tests.md index 005d329..2b10b48 100644 --- a/doc/api/tests.md +++ b/doc/api/tests.md @@ -7,14 +7,14 @@ Pure-math tests, only Eigen required. Covers Java utilities ported in Phase 1– | File | What it tests | |---|---| | `test_clausen.cpp` | Clausen Cl₂, Lobachevsky Л, ImLi₂ — values at known points | -| `test_hyper_ideal_utility.cpp` | Tetrahedron volumes (Meyerhoff / Kolpakov–Mednykh) | +| `test_hyper_ideal_utility.cpp` | Tetrahedron volumes (Meyerhoff / Kolpakov–Mednykh) + Java golden-value oracle (Clausen/Л/ImLi₂, ζ₁₃/₁₄/₁₅/ζ, both volume formulas) | | `test_matrix_utility.cpp` | Matrix helpers | | `test_surface_curve_utility.cpp` | Surface curve utilities | | `test_discrete_elliptic_utility.cpp` | Discrete elliptic functions | | `test_p2_utility.cpp` | P2 projective utilities | | `test_hyper_ideal_visualization_utility.cpp` | Poincaré disk projection, circumcircle | -**Total: 23 tests, 0 skipped.** +**Total: 26 tests, 0 skipped.** --- @@ -30,7 +30,9 @@ All tests have CTest prefix `cgal.` (set via `TEST_PREFIX "cgal."` in CMakeLists | `ConformalMeshValidity` | `test_conformal_mesh.cpp` | 1 | CGAL validity for all factory meshes | | `HyperIdealFunctional` | `test_hyper_ideal_functional.cpp` | 7 | FD gradient checks + Hessian symmetry | | `SphericalFunctional` | `test_spherical_functional.cpp` | 13 | Angle formula + gradient + gauge-fix + cross-module Hessian check + closed-form edge-DOF oracle | +| `SphericalGoldenJava` | `test_spherical_functional.cpp` | 3 | Java golden oracles: (1) law-of-cosines angles + β relations + Л energy term; (2) full-mesh tetrahedron — per-vertex gradient + ΔE; (3) full-mesh **edge-DOF** gradient (vertex + edge components, Finding 3) vs raw `conformalEnergyAndGradient` | | `EuclideanFunctional` | `test_euclidean_functional.cpp` | 12 | Angle formula + gradient + cross-module Hessian check | +| `EuclideanGoldenJava` | `test_euclidean_functional.cpp` | 2 | Java golden oracles: (1) angle formula + 2·Л(α) energy term; (2) full-mesh tetrahedron — per-vertex gradient + ΔE vs real `EuclideanCyclicFunctional` | | `EuclideanHessian` | `test_euclidean_hessian.cpp` | 9 | Cotangent Laplacian structure, FD agreement, PSD, null space | | `SphericalHessian` | `test_spherical_hessian.cpp` | 8 | Derivative correctness, NSD at equilibrium | | `NewtonSolver` | `test_newton_solver.cpp` | 11 | Convergence: Euclidean ×3, Spherical ×4, HyperIdeal ×4 | @@ -48,7 +50,7 @@ All tests have CTest prefix `cgal.` (set via `TEST_PREFIX "cgal."` in CMakeLists | `HalfedgeUV` | `test_phase7.cpp` | 4 | Size = number of half-edges, seam consistency, boundary half-edges = 0 | | `PriorityBFS` | `test_phase7.cpp` | 3 | Success, no seam on open meshes, all vertices placed | | `NormaliseEuclidean` | `test_phase7.cpp` | 2 | UV centroid = 0, halfedge_uv centroid = 0 | -| `PeriodMatrix` | `test_phase7.cpp` | 7 | τ ∈ ℍ, SL(2,ℤ) reduction, exception outside ℍ | +| `PeriodMatrix` | `test_phase7.cpp` | 8 | τ ∈ ℍ, SL(2,ℤ) reduction, exception outside ℍ, Java `normalizeModulus` golden oracle | | `FundamentalDomain` | `test_phase7.cpp` | 7 | Genus-1 parallelogram CCW, generators, g > 1 empty | | `TilingCopy/Neighbourhood` | `test_phase7.cpp` | 4 | Translation correct, tile count | | `HolonomyEndToEnd` | `test_phase7.cpp` | 3 | Full pipeline τ on tori of revolution (4×4, hex 6×6, 8×8) vs analytic modulus i·√(R²−r²)/r | @@ -67,7 +69,7 @@ All tests have CTest prefix `cgal.` (set via `TEST_PREFIX "cgal."` in CMakeLists | `NewtonPhase9a` | `test_newton_phase9a.cpp` | 7 | Phase 9a-Newton — convergence for the two new circle-packing solvers | | `CGALPhase8bLite` | `test_cgal_phase8b_lite.cpp` | 17 | Phase 8b-Lite — CGAL entries for all 5 DCE models + `output_uv_map` (Euclidean, Spherical, HyperIdeal, Inversive-Distance) + CP-Euclidean throws-clearly + pipe-operator chaining | -**Total: 240 tests, 0 skipped.** +**Total: 246 tests, 0 skipped.** --- diff --git a/doc/reviewer/java-port-audit.md b/doc/reviewer/java-port-audit.md index baf1a52..77b8352 100644 --- a/doc/reviewer/java-port-audit.md +++ b/doc/reviewer/java-port-audit.md @@ -124,18 +124,26 @@ cycles, not unweighted BFS). requires integration of holomorphic differentials — not implemented") so it is consistent with the documented contract, not a regression. -### 8. ℹ️ NOTE — Holonomy translation assumes zero residual rotation +### 8. ✅ FIXED — Holonomy now fits the full deck isometry (was: midpoint only) `detail::euclidean_holonomy` (`layout.hpp`) develops each face along a dual tree -that never crosses a cut edge, then reads each generator's translation as the -**midpoint displacement** `ω = midA − midB` of the shared cut edge between its two -independent developments. This equals the true deck-translation **only when the -two developments differ by a pure translation** (linear part = identity), i.e. for -a perfectly flat cone metric (Θ ≡ 2π). When Newton has not fully converged (small -residual cone-angle defect), there is a residual rotation and the midpoint -displacement is a first-order approximation rather than the exact ω. A fully -robust version would fit the rigid motion mapping edge `(Bs,Bt)↦(At,As)` and read -its translation part. Acceptable for converged flat tori (the `torus_8x8` test -passes); worth tightening if higher-genus or under-converged inputs are used. +that never crosses a cut edge. It **previously** read each generator's +translation as the **midpoint displacement** `ω = midA − midB` of the shared cut +edge between its two independent developments — exact **only** when the two +developments differ by a pure translation (linear part = identity), i.e. a +perfectly flat cone metric (Θ ≡ 2π). Under residual cone-angle defect that was +a first-order approximation. + +**Fix (2026-05-29):** it now fits the full orientation-preserving deck isometry +`g(z) = a·z + b` to the correspondence `(Bs,Bt) ↦ (As,At)` (`a = (At−As)/(Bt−Bs)`, +`b = As − a·Bs`). Because both faces use the same edge length, `|a| = 1` exactly, +so the fit is an exact isometry; `ω = b` is its translation part and `|arg(a)|` +is exposed per cut edge as `HolonomyData::residual_rotation` (a convergence +diagnostic). For a flat cone metric `a = 1` and `b` reduces to `midA − midB`, so +converged flat tori are **bit-for-bit unchanged** — all 240 CGAL tests still pass, +including `HolonomyEndToEnd.{Torus4x4,TorusHex6x6,Torus8x8}`. Note: when +`|arg(a)|` is non-negligible the holonomy is a genuine rotation and *no* lattice +translation is well-defined; the diagnostic now surfaces that instead of silently +averaging. --- @@ -200,37 +208,75 @@ Audited against the literature and found faithful: (This also fixed the Euclidean holonomy extraction — develop across the dual tree only — see `layout.hpp` / `cut_graph.hpp`.) -4. ⚠️ **PARTLY DONE (2026-05-29)** — added `EdgeGradient_RegularTetClosedForm`: - an *independent closed-form* oracle that pins each edge-DOF gradient to - `π/3` (= 2·(2π/3) − π) on the regular spherical tetrahedron, where the corner - angle 2π/3 follows from the spherical law of cosines. This locks the - Finding-3 formula `α_opp⁺ + α_opp⁻ − θ_e` against a value the path-integral FD - check cannot detect (the energy is the integral of G, so FD only proves - curl-freeness). **Still open:** (i) Newton-to-convergence *with* edge DOFs is - blocked by the Finding-4 spherical-Hessian guard (`throw` on edge DOFs), so a - converged-metric test needs an FD-Hessian or guard relaxation first; (ii) a - live-Java golden-value oracle still requires running the upstream library. +4. ✅ **DONE (2026-05-29)** — the Finding-3 spherical edge-DOF gradient is now + locked at two levels: (a) `EdgeGradient_RegularTetClosedForm` — an independent + *closed-form* oracle pinning each edge-DOF gradient to `π/3` (= 2·(2π/3) − π) + on the regular spherical tetrahedron; and (b) + `SphericalGoldenJava.FullMeshEdgeDofGradient_Tetrahedron` — a *live-Java* + full-mesh oracle (item 5) on a generic asymmetric configuration that pins + BOTH the vertex gradient AND the edge gradient `α_opp⁺ + α_opp⁻ − θ_e` to + 1e-12 vs the real raw `conformalEnergyAndGradient`. The vertex-DOF spherical + path is likewise locked by `SphericalGoldenJava.FullMeshGradientAndEnergy_…`. + **Still open (solver feature, not a correctness gap):** Newton-to-convergence + *with* edge DOFs is blocked by the Finding-4 spherical-Hessian guard (`throw` + on edge DOFs), so a converged-metric test needs an FD-Hessian or guard + relaxation first. -5. **Cross-check vs Java numeric oracles** — there is no test that pins C++ - functional/gradient values against recorded Java outputs on a fixed mesh. A - handful of golden-value tests (energy + gradient at a known `x`) would catch - any future silent divergence from the reference far more directly than the - self-consistent FD checks (which only verify curl-freeness, since the C++ - energy is itself the path-integral of its own gradient). +5. ✅ **DONE (2026-05-29)** — golden-value oracles now span both levels: + • *pure-math / per-triangle* (1e-12 vs compiled upstream, openjdk 17): + `HyperIdealGoldenJava` (Clausen/Л/ImLi₂, ζ₁₃/₁₄/₁₅/ζ, both tetrahedron-volume + formulas), `EuclideanGoldenJava` (angle formula + 2·Л energy), + `SphericalGoldenJava` (law-of-cosines angles + β relations + Л energy), + `PeriodMatrix.NormalizeModulus_GoldenJava` (item 8). + • *full-mesh energy + gradient at a known `x`* — the originally-fragile part, + now solved via a shared OBJ tetrahedron loaded identically in both languages: + `EuclideanGoldenJava.FullMeshGradientAndEnergy_Tetrahedron` and + `SphericalGoldenJava.FullMeshGradientAndEnergy_Tetrahedron` drive the REAL + `EuclideanCyclicFunctional` / `SphericalFunctional` and pin every per-vertex + gradient `Θ−Σα` (matched by vertex position) AND `ΔE = E(x)−E(0)` to 1e-12. + The ΔE check is doubly independent: C++ computes energy as a Gauss-Legendre + PATH INTEGRAL of its gradient, Java as a CLOSED-FORM sum — the initialEnergy + constant and the φ·λ⁰ term cancel in the difference (no edge DOFs here). + • *full-mesh EDGE-DOF gradient* — + `SphericalGoldenJava.FullMeshEdgeDofGradient_Tetrahedron` makes two opposite + tetrahedron edges variable (replacement parameterization λ_e = x[e_idx]) and + pins BOTH the vertex gradient AND the edge gradient `α_opp⁺+α_opp⁻−θ_e` to + 1e-12 vs the real raw `conformalEnergyAndGradient`, closing the Finding-3 + solution-level gap (item 4). It is a pure gradient oracle (no ΔE): with an + edge DOF, x=0 ⇒ λ_e=0 ⇒ arc length π is degenerate, so the path-integral-from- + origin energy reference is ill-defined; the fixed-x gradient is the + unambiguous quantity. It never touches the Hessian, so the Finding-4 edge-DOF + Hessian guard does not block it. Harness: `/tmp/oracle/SphereEdgeOracle.java`. + **Subtlety surfaced:** the spherical oracle must call Java's raw + `conformalEnergyAndGradient`, not `evaluate()`; `evaluate()` first runs a 1-D + Brent maximization over the global-scale gauge (`maximizeInNegativeDirection`), + which the C++ `spherical_gradient` deliberately omits (C++ factors that gauge + into the Newton solver's `spherical_gauge_shift`). Harnesses: + `/tmp/oracle/{tet.obj,EucMeshOracle.java,SphereMeshOracle.java,SphereEdgeOracle.java}`. + The *edge-DOF* full-mesh gradient oracle is now also done (item 4, + `SphereEdgeOracle.java`). 6. **CP-Euclidean `clausen` vs `clausen2`** — a single test asserting `|clausen(x) − clausen2(x)| < 1e-12` across a sweep would document Finding 5 and guard against an approximation regression. -7. **Period-matrix domain convention (Finding 6)** — production now uses - `normalizeModulus` (folds `Re ≥ 0`). Add a torus whose true τ has `Re(τ) < 0` - before reduction and assert the resulting `Re(τ) ≥ 0` convention, so any future - switch back to `reduce_to_fundamental_domain` is caught. Pair with a golden τ - recorded from the Java `DiscreteEllipticUtility` on the same mesh. +7. ⚠️ **PARTLY DONE (2026-05-29)** — the `Re ≥ 0` fold convention is now locked at + the function level by item 8 (`NormalizeModulus_GoldenJava`, including `Re < 0` + inputs). **Still open:** an *end-to-end* torus whose true τ has `Re(τ) < 0` + before reduction, asserting the production pipeline yields `Re(τ) ≥ 0` (so any + future switch back to `reduce_to_fundamental_domain` is caught at the mesh + level, not just the function level). -8. **`normalizeModulus` vs Java oracle (Finding 6)** — the faithful - `normalizeModulus` is only checked against ad-hoc values. Add a few recorded - Java `normalizeModulus` outputs as golden values to lock the `Re ≥ 0` fold. +8. ✅ **DONE (2026-05-29)** — `PeriodMatrix.NormalizeModulus_GoldenJava` pins the + C++ `normalizeModulus` against five recorded outputs of the real Java + `DiscreteEllipticUtility.normalizeModulus` (1e-12), including `Re < 0` mirror + folds, large-`Re` `T`-steps, and `|τ| < 1` `S`-inversions — locking the + `0 ≤ Re ≤ ½` fold convention. (Generated via `/tmp/oracle/TauOracle.java`.) + This also surfaced a real interaction: `compute_period_matrix` can land + exactly on `Re = +½` (the closed right edge of `normalizeModulus`'s mirror- + folded domain), which the half-open SL(2,ℤ) `is_in_fundamental_domain` + (`−½ ≤ Re < ½`, P1-3) correctly *excludes*; `ComputePeriodMatrix_ReducedTau_InFD` + was updated to assert the `normalizeModulus` domain rather than the SL(2,ℤ) one. 9. **Higher-genus cut graph (Finding 7)** — only genus-1 is exercised end-to-end. Add a genus-2 mesh and assert `compute_cut_graph` returns exactly `2g = 4` cut