From ca936b7652fc40e8b7a824bf9b7ac9516ddd04ad Mon Sep 17 00:00:00 2001 From: Tarik Moussa Date: Thu, 28 May 2026 07:41:23 +0200 Subject: [PATCH 1/2] docs: note high-precision requirement for Phase 9c/10 uniformization MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The Java uniformization classes (FundamentalPolygonUtility, CanonicalFormUtility) rely on the *Big arbitrary-precision geometry (MathContext(50)) because products of hyperbolic isometry generators grow exponentially and double fails to verify the group relation ∏gᵢ = Id. Record this planning-relevant prerequisite and resolve the contradiction in java-parity.md, which previously listed *Big as permanently out of scope. - CLAUDE.md: † note on the Phase-9 not-yet-ported table - java-parity.md: *Big exception (localized high-precision substrate for 9c/10) - phases.md: precision prerequisite as 9c sub-task + effort estimate - design-decisions.md: new "Scalar type: double, with one localized exception" Core flattening (Newton/energy/Eigen solver) stays double; the substrate (cpp_dec_float_50 / mpreal) is localized to the uniformization module only. Co-Authored-By: Claude Opus 4.7 --- CLAUDE.md | 6 ++++-- doc/architecture/design-decisions.md | 20 ++++++++++++++++++++ doc/roadmap/java-parity.md | 11 ++++++++++- doc/roadmap/phases.md | 9 ++++++++- 4 files changed, 42 insertions(+), 4 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 59458b4..3b35b4c 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -191,13 +191,15 @@ The Java library under `de.varylab.discreteconformal` contains these items not y |---|---|---| | `InversiveDistanceFunctional` | `inversive_distance_functional.hpp` | 9a | | Analytic HyperIdeal Hessian | `hyper_ideal_hessian.hpp` (replace FD) | 9b | -| 4g-polygon boundary walk in `FundamentalDomainUtility` | `fundamental_domain.hpp` (extend) | 9c | +| 4g-polygon boundary walk in `FundamentalDomainUtility` | `fundamental_domain.hpp` (extend) | 9c † | | `DiscreteHarmonicFormUtility` | Phase 10a prerequisite | 10 | | `DiscreteHolomorphicFormUtility` | Phase 10a | 10 | -| `HomologyUtility`, `CanonicalBasisUtility` | Phase 10 | 10 | +| `HomologyUtility`, `CanonicalBasisUtility` | Phase 10 † | 10 | When porting a Java class, locate the original in `de.varylab.discreteconformal.*` at [github.com/varylab/conformallab](https://github.com/varylab/conformallab) and use it as the reference implementation. +**† High-precision requirement (Phase 9c / 10):** The Java uniformization classes (`FundamentalPolygon`, `CanonicalFormUtility`) use `RnBig`/`PnBig`/`P2Big` with `MathContext(50)` — 50 significant decimal digits. Reason: products of hyperbolic isometry generators grow exponentially, so `double` fails when verifying the group relation ∏gᵢ = Id. When porting, replicate this **locally** with `boost::multiprecision::cpp_dec_float_50` (or MPFR `mpreal`) — only inside the uniformization module, NOT globally and NOT in the Eigen solver. The core flattening (Newton/energy) stays `double` (see `conformal_mesh.hpp:45`). + ## Test design patterns ### "Natural theta" — constructing a known equilibrium at x* = 0 diff --git a/doc/architecture/design-decisions.md b/doc/architecture/design-decisions.md index 9e40bb2..6416071 100644 --- a/doc/architecture/design-decisions.md +++ b/doc/architecture/design-decisions.md @@ -46,6 +46,26 @@ and linear system infrastructure to serve all three without branching. --- +## Scalar type: `double`, with one localized exception + +The kernel is `CGAL::Simple_cartesian` and the whole numerical core +(functionals, Hessians, Newton, Eigen sparse solvers) is `double`. Conformal +flattening minimises a smooth energy over transcendental (floating-point) lengths +and angles, so exact/extended arithmetic buys nothing there — the Java original is +`double` here too. The CGAL wrapper (`include/CGAL/*.h`) is templated on `FT` for +API genericity, but delegates to the `double` core. + +**The one exception (deferred, Phase 9c/10):** hyperbolic uniformization +(`FundamentalPolygonUtility`, `CanonicalFormUtility`) composes products of isometry +generators whose entries grow exponentially; `double` then fails to verify the group +relation ∏gᵢ = Id. The Java original handles this with `RnBig`/`PnBig`/`P2Big` at +`MathContext(50)`. When ported, this needs a **localized** high-precision substrate +(`boost::multiprecision::cpp_dec_float_50` or MPFR `mpreal`) **inside the +uniformization module only** — never the core or the Eigen solver. See `CLAUDE.md` +Phase-9 note and `doc/roadmap/java-parity.md` (`*Big` exception). + +--- + ## Priority-BFS layout A naive BFS layout places faces in arbitrary order; trilateration errors accumulate diff --git a/doc/roadmap/java-parity.md b/doc/roadmap/java-parity.md index dd357f7..e2cc0d1 100644 --- a/doc/roadmap/java-parity.md +++ b/doc/roadmap/java-parity.md @@ -89,7 +89,16 @@ items were unrecorded and have now been added above: Everything else unmentioned is intentionally out of scope: the `plugin/*` jReality Swing GUI, the `MTJ*`/`Tao*`/`*PETSc` solver bindings (replaced by Eigen), the convergence/`*Series` test harness, and the `*Big` arbitrary-precision geometry -(`P2Big`/`PnBig`/`RnBig` — deliberately dropped in favour of `Simple_cartesian`). +(`P2Big`/`PnBig`/`RnBig` — deliberately dropped from the core in favour of `Simple_cartesian`). + +> **Exception — `*Big` is NOT permanently out of scope.** The Phase 9c/10 uniformization +> classes (`FundamentalPolygonUtility`, `CanonicalFormUtility`, rows 55–56 above) build on +> the `*Big` classes in Java for a reason: products of hyperbolic isometry generators grow +> exponentially, so `double` fails when verifying the group relation ∏gᵢ = Id (`MathContext(50)` +> = 50 significant digits in Java). When those classes are ported, a localized high-precision +> substrate (`boost::multiprecision::cpp_dec_float_50` or MPFR `mpreal`) must be reintroduced — +> **only inside the uniformization module**, not in the core flattening or the Eigen solver. +> See the Phase-9 note in `CLAUDE.md`. --- diff --git a/doc/roadmap/phases.md b/doc/roadmap/phases.md index ae311f1..84c3c8b 100644 --- a/doc/roadmap/phases.md +++ b/doc/roadmap/phases.md @@ -131,8 +131,15 @@ mesh type. Mathematical source: Poincaré 1882 + Sechelmann 2016 §5 Research component: bridging to conformallab++ cut_graph.hpp + holonomy infrastructure. + † PRECISION PREREQUISITE: the canonicalisation composes products of + hyperbolic isometry generators, whose entries grow exponentially. + double fails to verify the group relation ∏gᵢ = Id; Java uses + MathContext(50) (RnBig/PnBig/P2Big). Port a LOCALIZED high-precision + substrate (boost::multiprecision::cpp_dec_float_50 or MPFR mpreal) + inside the uniformization module only — NOT the core or the Eigen + solver. See CLAUDE.md Phase-9 note + java-parity.md exception. Effort: ~2 weeks for fundamental polygon, +2 weeks for surgery - layer, +1 week integration. + layer, +1 week integration, +~3 days high-precision substrate. ``` 9d — Cone metrics + sphere utilities (Java port + research extension) From a3ee9576d4df7564ca62b8a3a50054a847802a8c Mon Sep 17 00:00:00 2001 From: Tarik Moussa Date: Fri, 29 May 2026 12:50:16 +0200 Subject: [PATCH 2/2] =?UTF-8?q?fix+test:=20Euclidean=20holonomy/=CF=84=20e?= =?UTF-8?q?nd-to-end=20+=20spherical=20edge-DOF=20oracle=20(2026-05-29=20a?= =?UTF-8?q?udit)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Bundles the 2026-05-29 Java↔C++ math-correctness audit (doc/reviewer/ java-port-audit.md, 11 findings) with two follow-up fixes. Audit code changes: - Finding 3 (spherical_functional): edge-DOF replacement parameterization via spher_eff_lambda; edge gradient α_opp⁺+α_opp⁻−θ_e (drops additive −(S⁺+S⁻)/2) - Finding 4 (spherical_hessian): always-compiled edge-DOF throw guard - Finding 6 (period_matrix): faithful normalizeModulus (0≤Re≤½, Im≥0, |τ|≥1) - Finding 9 (inversive_distance): degenerate-face limiting angles, no skip - Findings 1/2 (euclidean): degenerate gradient limiting angles + Hessian guard Euclidean holonomy/τ fix: develop the cut surface across the dual spanning tree only (cut_graph now exposes is_dual_tree), so genus-1 cut edges yield non-degenerate lattice generators. Previously τ came out 0 / NaN / 1e13 on the bundled tori; now matches the analytic revolution modulus i·√(R²−r²)/r. Re-enabled τ reporting in the Euclidean CLI; rewrote validation.md §3/§4 accordingly. Tests (240 CGAL, 0 skipped): - HolonomyEndToEnd ×3 — tori of revolution (4×4, hex 6×6, 8×8) vs analytic modulus - SphericalFunctional.EdgeGradient_RegularTetClosedForm — independent closed-form π/3 oracle locking the Finding-3 edge formula (the path-integral FD check cannot detect a wrong-but-conservative gradient) Also documents the latent spherical/hyperbolic holonomy-extraction bug (same single-development pattern, dead code today) in research-track.md (Phase 9c/10), and adds favour/normalisations to the codespell ignore list. Co-Authored-By: Claude Opus 4.7 --- .codespellrc | 3 +- CLAUDE.md | 11 + README.md | 11 +- code/include/cut_graph.hpp | 20 ++ code/include/euclidean_functional.hpp | 6 +- code/include/euclidean_geometry.hpp | 13 +- code/include/euclidean_hessian.hpp | 13 + code/include/gauss_bonnet.hpp | 13 +- .../include/inversive_distance_functional.hpp | 9 +- code/include/layout.hpp | 155 ++++++++++- code/include/mesh_io.hpp | 10 + code/include/newton_solver.hpp | 130 +++++++--- code/include/period_matrix.hpp | 15 +- code/include/spherical_functional.hpp | 72 +++--- code/include/spherical_geometry.hpp | 16 +- code/include/spherical_hessian.hpp | 13 + code/src/apps/v0/conformallab_cli.cpp | 133 +++++++--- code/tests/cgal/test_phase7.cpp | 120 +++++++++ code/tests/cgal/test_spherical_functional.cpp | 61 ++++- doc/api/tests.md | 5 +- doc/math/validation.md | 143 ++++++---- doc/reviewer/java-port-audit.md | 244 ++++++++++++++++++ doc/roadmap/research-track.md | 14 + 23 files changed, 1065 insertions(+), 165 deletions(-) create mode 100644 doc/reviewer/java-port-audit.md diff --git a/.codespellrc b/.codespellrc index 6e7b5b3..b6a5455 100644 --- a/.codespellrc +++ b/.codespellrc @@ -27,7 +27,8 @@ ignore-words-list = bessel,ist,sinces,nd,te,inout,nin,numer,neet,anc,sinks,doubl behaviour,behaviours,behavioural, analogue,analogues, initialise,initialised,initialises,initialising,initialisation, - normalise,normalised,normalises,normalising,normalisation, + normalise,normalised,normalises,normalising,normalisation,normalisations, + favour,favours,favoured,favouring,favourite,favourites, centralise,centralised,centralises,centralising, serialise,serialised,serialises,serialising,serialisation, parameterise,parameterised,parameterises,parameterising, diff --git a/CLAUDE.md b/CLAUDE.md index 3b35b4c..6fe13ae 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -300,6 +300,17 @@ grep -r "ClassName" /Users/tarikmoussa/Desktop/conformallab/src The 2026-05-21 audit corrected four mis-labels (now fixed): `InversiveDistanceFunctional`, both HyperIdeal Hessians (FD + analytic), and the inversive-distance tutorial were all wrongly tagged "Java port" — they are research (no Java parent; Java declares `hasHessian()==false`). Details in `research-track.md`. +### 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`): + +- **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\.'`. + ## Documentation map 38 documents across 7 categories. Read the relevant one before reasoning from scratch diff --git a/README.md b/README.md index b00318c..80473e0 100644 --- a/README.md +++ b/README.md @@ -34,7 +34,16 @@ ctest --test-dir build -R "^cgal\." --output-on-failure # Full build with CLI + viewer (requires Wayland/X11 dev headers) cmake -S code -B build -DWITH_CGAL=ON && cmake --build build -j$(nproc) -./bin/conformallab_core -i input.off -g euclidean -o layout.off -j result.json + +# Conformal flattening (Θ_v = 2π target). Closed meshes pin one vertex + +# enforce Gauss-Bonnet; open meshes pin the boundary and flatten the interior. +./bin/conformallab_core -i code/data/off/torus_8x8.off -g euclidean -v \ + -o layout.off -j result.json +# topology: closed, free DOFs=63, genus=1 +# Euclidean: converged=yes iter=3 |grad|_inf≈5e-15 +./bin/conformallab_core -i code/data/obj/cathead.obj -g euclidean -v -o cat.off +# topology: open (boundary pinned), free DOFs=119 +# Euclidean: converged=yes iter=4 # API documentation (requires doxygen: brew/apt install doxygen) cmake --build build --target doc diff --git a/code/include/cut_graph.hpp b/code/include/cut_graph.hpp index febb41f..d2731a6 100644 --- a/code/include/cut_graph.hpp +++ b/code/include/cut_graph.hpp @@ -49,6 +49,13 @@ struct CutGraph { /// Indices of the 2g cut edges in order (size = 2g). std::vector cut_edge_indices; + /// dual_tree_edge_flags[e.idx()] = true ↔ edge `e` is a dual-spanning-tree + /// edge (its dual is in T*). Size = mesh.number_of_edges(). + /// Crossing only these edges develops the surface onto a topological disk + /// (the fundamental polygon), so that the `2g` cut edges become genuine + /// boundary identifications carrying the holonomy generators. + std::vector dual_tree_edge_flags; + /// Genus of the surface (0 for topological spheres and open patches). int genus = 0; @@ -58,6 +65,14 @@ struct CutGraph { return static_cast(e.idx()) < cut_edge_flags.size() && cut_edge_flags[static_cast(e.idx())]; } + + /// `true` iff edge `e` is a dual-spanning-tree edge (crossable when + /// developing the surface onto a disk). + bool is_dual_tree(Edge_index e) const + { + return static_cast(e.idx()) < dual_tree_edge_flags.size() + && dual_tree_edge_flags[static_cast(e.idx())]; + } }; /// Compute the cut graph of `mesh` via the standard tree-cotree @@ -145,6 +160,11 @@ inline CutGraph compute_cut_graph(const ConformalMesh& mesh) cg.cut_edge_indices.push_back(idx); } + // Expose the dual spanning tree T*: developing across only these edges + // unfolds the surface onto a disk, making the cut edges the boundary + // identifications that carry the holonomy generators. + cg.dual_tree_edge_flags = std::move(dual_tree_edge); + return cg; } diff --git a/code/include/euclidean_functional.hpp b/code/include/euclidean_functional.hpp index a20c195..37e766f 100644 --- a/code/include/euclidean_functional.hpp +++ b/code/include/euclidean_functional.hpp @@ -202,7 +202,11 @@ inline std::vector euclidean_gradient( double lam31 = m.lambda0[e31] + u3 + u1 + eucl_dof_val(m.e_idx[e31], x); auto fa = euclidean_angles(lam12, lam23, lam31); - if (!fa.valid) continue; // degenerate triangle: contributes 0 + // NOTE: do NOT skip degenerate faces. euclidean_angles() returns the + // limiting angles (one corner = π, others = 0) when the triangle + // inequality is violated; using them is exactly what the Java reference + // does and is required for the BPS energy to be the convex C¹ extension + // onto the infeasible region (otherwise Newton can stall at a flip). // h_alpha[h] = corner angle OPPOSITE to h's edge: // h0 (edge v1v2) → opposite corner at v3 → α3 diff --git a/code/include/euclidean_geometry.hpp b/code/include/euclidean_geometry.hpp index f78a0f8..d9f8f69 100644 --- a/code/include/euclidean_geometry.hpp +++ b/code/include/euclidean_geometry.hpp @@ -28,6 +28,7 @@ // rescales all three sides by the same factor, leaving angles unchanged but // keeping the arguments of exp in a safe numerical range. +#include "constants.hpp" #include namespace conformallab { @@ -50,8 +51,16 @@ inline EuclideanFaceAngles euclidean_angles_from_lengths( const double t23 = +l12 - l23 + l31; // 2*(s − l23) const double t31 = +l12 + l23 - l31; // 2*(s − l31) - if (t12 <= 0.0 || t23 <= 0.0 || t31 <= 0.0) - return {0.0, 0.0, 0.0, false}; + // Degenerate (triangle inequality violated): return the *limiting* angles + // of the flat-out triangle — the corner opposite the over-long edge is π, + // the other two are 0. This is the convex C¹ extension of the BPS energy + // onto the infeasible region and is what the Java reference + // (EuclideanCyclicFunctional.triangleEnergyAndAlphas) does. `valid` stays + // false so the cotangent Hessian still skips this face. At most one t can + // be ≤ 0 (t12+t23 = 2·l31 > 0, etc.), so the order of these checks is moot. + if (t23 <= 0.0) return {PI, 0.0, 0.0, false}; // l23 too long → α₁ = π + if (t31 <= 0.0) return {0.0, PI, 0.0, false}; // l31 too long → α₂ = π + if (t12 <= 0.0) return {0.0, 0.0, PI, false}; // l12 too long → α₃ = π const double l123 = l12 + l23 + l31; const double denom2 = t12 * t23 * t31 * l123; // = (4·Area)² diff --git a/code/include/euclidean_hessian.hpp b/code/include/euclidean_hessian.hpp index ed87cdd..e9a61c3 100644 --- a/code/include/euclidean_hessian.hpp +++ b/code/include/euclidean_hessian.hpp @@ -44,6 +44,7 @@ #include #include #include +#include namespace conformallab { @@ -104,6 +105,18 @@ inline Eigen::SparseMatrix euclidean_hessian( { const int n = euclidean_dimension(mesh, m); + // Only the vertex block of the cyclic Hessian is implemented here. If any + // edge DOF is variable (assign_euclidean_all_dof_indices), the edge-edge and + // vertex-edge blocks present in the Java reference (conformalHessian) are + // missing, which would leave singular zero rows/cols. Fail loudly instead + // of silently returning a rank-deficient matrix (matches the doc contract). + for (auto e : mesh.edges()) { + if (m.e_idx[e] >= 0) + throw std::logic_error( + "euclidean_hessian: edge DOFs are not supported " + "(only the vertex-block cotangent Laplacian is implemented)"); + } + // Collect triplets (row, col, value) — setFromTriplets sums duplicates. std::vector> trips; trips.reserve(static_cast(n) * 7); // rough estimate diff --git a/code/include/gauss_bonnet.hpp b/code/include/gauss_bonnet.hpp index 15825f0..1d92c68 100644 --- a/code/include/gauss_bonnet.hpp +++ b/code/include/gauss_bonnet.hpp @@ -16,6 +16,13 @@ // If this fails, no conformal factor can realise the target angles and // Newton will silently fail to converge. // +// PRECONDITION — closed meshes only. Every function here sums (2π − Θ_v) +// over ALL vertices. On a mesh with boundary the boundary vertices carry a +// (π − Θ_v) term instead, so the identity Σ(2π−Θ_v) = 2π·χ does NOT hold and +// `check_gauss_bonnet` will (correctly) throw. For open meshes pin the +// boundary directly and skip the Gauss–Bonnet check (see the CLI's +// flattening path). +// // API: // int euler_characteristic(mesh) // int genus(mesh) @@ -128,10 +135,10 @@ inline void check_gauss_bonnet(const ConformalMesh& mesh, // ── enforce_gauss_bonnet — adjust θ_v by uniform Δ ─────────────────────────── // -// Adds δ = (rhs − lhs) / V to every θ_v so that Gauss–Bonnet holds exactly. +// Adds δ = (lhs − rhs) / V to every θ_v so that Gauss–Bonnet holds exactly. // After this call, check_gauss_bonnet() will not throw (up to floating-point). -// Only modifies free vertices (v_idx[v] >= 0 for EuclideanMaps / SphericalMaps; -// always all vertices for the raw property-map overload). +// Modifies ALL vertices' θ_v (no v_idx filtering) — the shift is a property +// of the target angles, independent of which vertices are free DOFs. /// Distribute the Gauss-Bonnet deficit uniformly across all `Θ_v`: /// add `δ = (lhs − rhs) / V` to every entry so that the identity holds diff --git a/code/include/inversive_distance_functional.hpp b/code/include/inversive_distance_functional.hpp index 2f709f6..2d48e3e 100644 --- a/code/include/inversive_distance_functional.hpp +++ b/code/include/inversive_distance_functional.hpp @@ -269,11 +269,18 @@ inline std::vector inversive_distance_gradient( double l23sq = id_detail::edge_length_squared(u2, u3, m.I_e[e23]); double l31sq = id_detail::edge_length_squared(u3, u1, m.I_e[e31]); + // A non-real circle configuration (ℓ² ≤ 0) has no limiting angle — skip it. if (l12sq <= 0 || l23sq <= 0 || l31sq <= 0) continue; // euclidean_angles expects 2·log(ℓ) per edge — feed log(ℓ²). + // For a triangle-inequality-violating face euclidean_angles returns the + // *limiting* angles (π opposite the over-long edge, 0/0 otherwise) with + // valid=false. We deliberately do NOT skip on !fa.valid: using those + // limiting angles is the convex C¹ extension onto the infeasible region, + // so Newton can pass through a flip instead of stalling (Finding 9 — + // mirrors the Euclidean/Spherical fix in Finding 1). The angles come + // from genuine Euclidean side lengths, so the extension is geometric. auto fa = euclidean_angles(std::log(l12sq), std::log(l23sq), std::log(l31sq)); - if (!fa.valid) continue; h_alpha[id_detail::hidx(h0)] = fa.alpha3; h_alpha[id_detail::hidx(h1)] = fa.alpha1; diff --git a/code/include/layout.hpp b/code/include/layout.hpp index 790a59c..2d8f9ed 100644 --- a/code/include/layout.hpp +++ b/code/include/layout.hpp @@ -474,6 +474,125 @@ inline void set_root_huv_2d( huv[static_cast(hf.idx())] = uv[static_cast(mesh.source(hf).idx())]; } +// ───────────────────────────────────────────────────────────────────────────── +// Euclidean holonomy via the developing map (genus-g closed surfaces). +// +// The per-vertex layout produced by euclidean_layout() places every face in a +// SINGLE consistent global frame (each vertex is placed once). In that frame +// the holonomy is identically trivial — it is exactly the obstruction to such a +// single frame existing on the uncut surface. To recover it we develop every +// face INDEPENDENTLY along a spanning tree of the dual graph that does not cross +// any cut edge, storing each face's own copy of its three corner positions +// (`hpos[h]` = global position of source(h) as developed inside face(h)). +// +// Two faces adjacent across a cut edge are NOT tree-adjacent, so they are +// developed via different tree branches and place the shared edge at two +// 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(). +template +inline std::vector euclidean_holonomy( + const ConformalMesh& mesh, + const CutGraph& cut, + EdgeLenFn&& edge_len) +{ + using C = std::complex; + const std::size_t nh = mesh.number_of_halfedges(); + const std::size_t nf = mesh.number_of_faces(); + + std::vector hpos(nh, C(0.0, 0.0)); // pos of source(h) inside face(h) + std::vector face_done(nf, false); + + auto place_root = [&](Face_index f) { + Halfedge_index h0 = mesh.halfedge(f); + Halfedge_index h1 = mesh.next(h0), h2 = mesh.next(h1); + double lAB = edge_len(h0), lBC = edge_len(h1), lCA = edge_len(h2); + Eigen::Vector2d A(0.0, 0.0), B(lAB, 0.0); + Eigen::Vector2d Cc = trilaterate_2d(A, B, lCA, lBC); // apex = source(h2) + hpos[static_cast(h0.idx())] = C(A.x(), A.y()); + hpos[static_cast(h1.idx())] = C(B.x(), B.y()); + hpos[static_cast(h2.idx())] = C(Cc.x(), Cc.y()); + face_done[static_cast(f.idx())] = true; + }; + + auto develop = [&](Face_index root) { + place_root(root); + std::queue q; // halfedges pointing INTO unplaced faces + auto enqueue = [&](Face_index f) { + for (auto hf : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) { + Halfedge_index ho = mesh.opposite(hf); + if (mesh.is_border(ho)) continue; + // Develop across the dual spanning tree T* ONLY. Crossing any + // other edge (a cut/generator edge OR a primal-tree edge) would + // over-connect the development: the surface minus the 2g cut + // edges is still non-simply-connected, so the immersion would + // wrap around and place the two copies of a generator edge on + // top of each other (zero/garbage holonomy). Crossing only T* + // unfolds the surface onto a disk (the fundamental polygon). + if (!cut.is_dual_tree(mesh.edge(hf))) continue; + Face_index fa = mesh.face(ho); + if (!face_done[static_cast(fa.idx())]) q.push(ho); + } + }; + enqueue(root); + while (!q.empty()) { + Halfedge_index h = q.front(); q.pop(); + Face_index f = mesh.face(h); + if (face_done[static_cast(f.idx())]) continue; + + Halfedge_index ho = mesh.opposite(h); + // Shared edge endpoints, taken from the PARENT face's development: + // source(h) = target(ho) → parent pos hpos[next(ho)] + // target(h) = source(ho) → parent pos hpos[ho] + C ps = hpos[static_cast(mesh.next(ho).idx())]; + C pt = hpos[static_cast(ho.idx())]; + Eigen::Vector2d A(ps.real(), ps.imag()), B(pt.real(), pt.imag()); + Eigen::Vector2d apex = trilaterate_2d( + A, B, edge_len(mesh.prev(h)), edge_len(mesh.next(h))); + + hpos[static_cast(h.idx())] = ps; + hpos[static_cast(mesh.next(h).idx())] = pt; + hpos[static_cast(mesh.prev(h).idx())] = C(apex.x(), apex.y()); + face_done[static_cast(f.idx())] = true; + enqueue(f); + } + }; + + develop(best_root_face(mesh)); + 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()); + 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); + Halfedge_index ho = mesh.opposite(h); + 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()); + continue; + } + // Face A = face(h) places edge e as halfedge h: S=source(h), T=target(h) + C As = hpos[static_cast(h.idx())]; + C At = hpos[static_cast(mesh.next(h).idx())]; + // Face B = face(ho) places the same edge as halfedge ho: source(ho)=T, + // 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())); + } + return omega; +} + } // namespace detail // ── Euclidean layout ────────────────────────────────────────────────────────── @@ -592,25 +711,30 @@ inline Layout2D euclidean_layout( result.success = true; // ── Holonomy ────────────────────────────────────────────────────────────── + // The single-frame BFS layout above puts every face in one consistent frame, + // so the holonomy read off it is identically trivial. Instead develop each + // face independently along a dual spanning tree that never crosses a cut edge + // (detail::euclidean_holonomy): the displacement of each cut edge between the + // two faces that share it is the lattice generator ω_i. if (cut && holonomy) { holonomy->cut_edge_indices = cut->cut_edge_indices; - holonomy->translations.clear(); - holonomy->translations.reserve(cut->cut_edge_indices.size()); holonomy->mobius_maps.clear(); + holonomy->translations = detail::euclidean_holonomy(mesh, *cut, edge_len); + + // Preserve the per-cut-edge seam UV in halfedge_uv (texture atlas), as + // before, so HalfedgeUV-based tests still see the seam-crossing layout. for (std::size_t ce_idx : cut->cut_edge_indices) { Edge_index e = *std::next(mesh.edges().begin(), static_cast(ce_idx)); Halfedge_index h = mesh.halfedge(e); Halfedge_index ho = mesh.opposite(h); Halfedge_index hx = mesh.is_border(ho) ? h : ho; - if (mesh.is_border(hx)) { holonomy->translations.push_back(Eigen::Vector2d::Zero()); continue; } + if (mesh.is_border(hx)) continue; Vertex_index vs = mesh.source(hx), vt = mesh.target(hx), vn = mesh.target(mesh.next(hx)); - if (!vertex_placed[vn.idx()]) { holonomy->translations.push_back(Eigen::Vector2d::Zero()); continue; } + if (!vertex_placed[vn.idx()]) continue; Eigen::Vector2d p_tri = detail::trilaterate_2d( result.uv[vs.idx()], result.uv[vt.idx()], edge_len(mesh.prev(hx)), edge_len(mesh.next(hx))); - // Store seam UV for the cut-crossing halfedges detail::set_face_huv_2d(result.halfedge_uv, mesh, hx, result.uv, p_tri); - holonomy->translations.push_back(p_tri - result.uv[vn.idx()]); } } @@ -707,6 +831,15 @@ inline Layout3D spherical_layout( for (auto f : mesh.faces()) if (!face_placed[f.idx()]) place_component(f); result.success = true; + // KNOWN LIMITATION (latent — every current caller passes holonomy == nullptr). + // This block extracts holonomy from a single full-surface development (the BFS + // above crosses every non-cut edge), then reads off an apex trilateration on one + // side of each seam. That is the same flawed pattern that produced garbage τ for + // the Euclidean path; the correct approach is detail::euclidean_holonomy, which + // develops across only the dual spanning tree (is_dual_tree) and measures the + // shared-edge displacement between two independent developments. Until a + // detail::spherical_holonomy mirror exists (Phase 9c/10, see research-track.md), + // these spherical translations are not trustworthy for genus g ≥ 1. if (cut && holonomy) { holonomy->cut_edge_indices = cut->cut_edge_indices; holonomy->translations.clear(); @@ -838,6 +971,16 @@ inline Layout2D hyper_ideal_layout( result.success = true; // ── Möbius-map holonomy ─────────────────────────────────────────────────── + // KNOWN LIMITATION (latent — every current caller passes holonomy == nullptr). + // Like the spherical block above, this reads the Möbius deck transformation from + // a single full-surface development (BFS crosses all non-cut edges) and one-sided + // apex trilateration — the same flawed pattern fixed for the Euclidean path by + // detail::euclidean_holonomy (develop across the dual tree only, measure seam + // displacement between two independent developments). A faithful + // detail::hyperbolic_holonomy is Phase 9c/10 work and additionally requires + // cpp_dec_float_50: products of these generators grow exponentially, so verifying + // the group relation ∏gᵢ = Id overflows double (see CLAUDE.md, research-track.md). + // Until then these mobius_maps are NOT correct for genus g ≥ 2 uniformization. if (cut && holonomy) { holonomy->cut_edge_indices = cut->cut_edge_indices; holonomy->translations.clear(); diff --git a/code/include/mesh_io.hpp b/code/include/mesh_io.hpp index 78d9948..48bf145 100644 --- a/code/include/mesh_io.hpp +++ b/code/include/mesh_io.hpp @@ -26,6 +26,7 @@ #include "conformal_mesh.hpp" #include +#include #include #include @@ -47,11 +48,20 @@ inline bool write_mesh(const std::string& filename, const ConformalMesh& mesh) /// Throwing wrapper around `read_mesh`: returns the mesh by value /// or throws `std::runtime_error` on read failure. +/// +/// Also enforces that the mesh is triangulated, because every functional +/// in this library assumes triangle faces — a quad/polygon mesh would be +/// read in silently and then mis-handled by the angle/length formulas. +/// Fail loudly here at the I/O boundary instead. inline ConformalMesh load_mesh(const std::string& filename) { ConformalMesh mesh; if (!read_mesh(filename, mesh)) throw std::runtime_error("conformallab: failed to read mesh from " + filename); + if (!CGAL::is_triangle_mesh(mesh)) + throw std::runtime_error( + "conformallab: mesh from " + filename + + " is not triangulated (functionals require triangle faces)"); return mesh; } diff --git a/code/include/newton_solver.hpp b/code/include/newton_solver.hpp index 2bba207..c8802f5 100644 --- a/code/include/newton_solver.hpp +++ b/code/include/newton_solver.hpp @@ -115,35 +115,89 @@ inline Eigen::VectorXd solve_linear_system( namespace detail { // re-open for the remaining helpers -// Backtracking line search: find the largest α in {1, 0.5, 0.25, …} such that -// ||G(x + α·Δx)||₂ < ||G(x)||₂. Returns the accepted step (α may stay 1). +// Globalised line search for the Newton system G(x) = 0. +// +// Merit function f(x) = ½‖G(x)‖²₂. Driving f to its minimum drives the +// residual G to zero; because the merit only depends on ‖G‖ it is sign-agnostic +// and works identically for the convex Euclidean / HyperIdeal / CP / InvDist +// energies and the concave Spherical energy. +// +// Phase 1 — Newton direction `dx` (satisfies H·dx = −G, so the merit slope +// ∇f·dx = GᵀH·dx = −‖G‖² < 0 — always a descent direction). +// Backtrack α ∈ {1, ½, ¼, …} until the Armijo sufficient-decrease +// condition holds: +// ‖G(x + α·dx)‖² ≤ (1 − 2·c1·α)·‖G(x)‖² +// +// Phase 2 — if Phase 1 exhausts its halvings, fall back to the steepest- +// descent direction of the merit, `d_sd = −H·G` (slope +// ∇f·d_sd = −‖H·G‖² ≤ 0 regardless of the definiteness of H), with +// the analogous Armijo test: +// ‖G(x + α·d_sd)‖² ≤ ‖G(x)‖² − 2·c1·α·‖d_sd‖² +// +// If neither phase satisfies Armijo, return the best (smallest-residual) point +// visited. If nothing beat ‖G(x)‖, return x unchanged and set *improved=false, +// so the caller can stop cleanly instead of taking the old divergent full step. template inline std::vector line_search( const std::vector& x, const Eigen::VectorXd& dx, + const Eigen::VectorXd& d_sd, double norm0, GradFn&& grad_fn, - int max_halvings = 20) + bool* improved = nullptr, + int max_halvings = 20, + double c1 = 1e-4) { - const int n = static_cast(x.size()); - double alpha = 1.0; - std::vector xnew(static_cast(n)); + const int n = static_cast(x.size()); + const double norm0_sq = norm0 * norm0; - for (int ls = 0; ls < max_halvings; ++ls) { + std::vector xnew(static_cast(n)); + std::vector best_x = x; + double best_norm = norm0; + + // Evaluate ‖G(x + α·dir)‖₂, leaving the trial point in `xnew`. + auto eval = [&](const Eigen::VectorXd& dir, double alpha) -> double { for (int i = 0; i < n; ++i) - xnew[static_cast(i)] = x[static_cast(i)] - + alpha * dx[i]; + xnew[static_cast(i)] = + x[static_cast(i)] + alpha * dir[i]; auto Gnew = grad_fn(xnew); - double norm_new = 0.0; - for (double v : Gnew) norm_new += v * v; - norm_new = std::sqrt(norm_new); - if (norm_new < norm0) return xnew; + double s = 0.0; + for (double v : Gnew) s += v * v; + return std::sqrt(s); + }; + + // ── Phase 1: Newton direction, Armijo backtracking ──────────────────────── + double alpha = 1.0; + for (int ls = 0; ls < max_halvings; ++ls) { + double norm_new = eval(dx, alpha); + if (norm_new < best_norm) { best_norm = norm_new; best_x = xnew; } + if (norm_new * norm_new <= (1.0 - 2.0 * c1 * alpha) * norm0_sq) { + if (improved) *improved = true; + return xnew; + } alpha *= 0.5; } - // No improvement found — return best attempt (full step) - for (int i = 0; i < n; ++i) - xnew[static_cast(i)] = x[static_cast(i)] + dx[i]; - return xnew; + + // ── Phase 2: steepest-descent fallback (−H·G), Armijo backtracking ──────── + const double dsd_sq = d_sd.squaredNorm(); + if (dsd_sq > 0.0) { + alpha = 1.0; + for (int ls = 0; ls < max_halvings; ++ls) { + double thresh_sq = norm0_sq - 2.0 * c1 * alpha * dsd_sq; + double norm_new = eval(d_sd, alpha); + if (norm_new < best_norm) { best_norm = norm_new; best_x = xnew; } + if (thresh_sq >= 0.0 && norm_new * norm_new <= thresh_sq) { + if (improved) *improved = true; + return xnew; + } + alpha *= 0.5; + } + } + + // ── Both phases failed Armijo — never take the divergent full step. ─────── + // Return the best point seen; if none improved, stay put and signal stall. + if (improved) *improved = (best_norm < norm0); + return best_x; } } // namespace detail @@ -209,12 +263,15 @@ inline NewtonResult newton_euclidean( Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok); if (!ok) break; - // ── Backtracking line search ────────────────────────────────────────── + // ── Globalised line search (Armijo + steepest-descent fallback) ─────── double norm0 = G.norm(); - x = detail::line_search(x, dx, norm0, + Eigen::VectorXd d_sd = -(H * G); // merit-function steepest descent + bool improved = true; + x = detail::line_search(x, dx, d_sd, norm0, [&](const std::vector& xnew) { return euclidean_gradient(mesh, xnew, m); - }); + }, &improved); + if (!improved) break; // line search stalled — stop cleanly res.iterations = iter + 1; } @@ -287,12 +344,16 @@ inline NewtonResult newton_spherical( Eigen::VectorXd dx = detail::solve_with_fallback(negH, G, ok); if (!ok) break; - // ── Backtracking line search ────────────────────────────────────────── + // ── Globalised line search (Armijo + steepest-descent fallback) ─────── + // d_sd uses the actual (un-negated) Hessian: ∇f = H·G for f = ½‖G‖². double norm0 = G.norm(); - x = detail::line_search(x, dx, norm0, + Eigen::VectorXd d_sd = -(H * G); + bool improved = true; + x = detail::line_search(x, dx, d_sd, norm0, [&](const std::vector& xnew) { return spherical_gradient(mesh, xnew, m); - }); + }, &improved); + if (!improved) break; res.iterations = iter + 1; } @@ -367,12 +428,15 @@ inline NewtonResult newton_hyper_ideal( Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok); if (!ok) break; - // ── Backtracking line search ────────────────────────────────────────── + // ── Globalised line search (Armijo + steepest-descent fallback) ─────── double norm0 = G.norm(); - x = detail::line_search(x, dx, norm0, + Eigen::VectorXd d_sd = -(H * G); + bool improved = true; + x = detail::line_search(x, dx, d_sd, norm0, [&](const std::vector& xnew) { return evaluate_hyper_ideal(mesh, xnew, m, false).gradient; - }); + }, &improved); + if (!improved) break; res.iterations = iter + 1; } @@ -443,10 +507,13 @@ inline NewtonResult newton_cp_euclidean( if (!ok) break; double norm0 = G.norm(); - x = detail::line_search(x, dx, norm0, + Eigen::VectorXd d_sd = -(H * G); + bool improved = true; + x = detail::line_search(x, dx, d_sd, norm0, [&](const std::vector& xnew) { return cp_euclidean_gradient(mesh, xnew, m); - }); + }, &improved); + if (!improved) break; res.iterations = iter + 1; } @@ -552,10 +619,13 @@ inline NewtonResult newton_inversive_distance( if (!ok) break; double norm0 = G.norm(); - x = detail::line_search(x, dx, norm0, + Eigen::VectorXd d_sd = -(H * G); + bool improved = true; + x = detail::line_search(x, dx, d_sd, norm0, [&](const std::vector& xnew) { return inversive_distance_gradient(mesh, xnew, m); - }); + }, &improved); + if (!improved) break; res.iterations = iter + 1; } diff --git a/code/include/period_matrix.hpp b/code/include/period_matrix.hpp index 3d8956c..0b081dd 100644 --- a/code/include/period_matrix.hpp +++ b/code/include/period_matrix.hpp @@ -41,6 +41,7 @@ // std::complex reduce_to_fundamental_domain(τ) — apply SL(2,ℤ) #include "layout.hpp" +#include "discrete_elliptic_utility.hpp" // normalizeModulus (Java-faithful) #include #include #include @@ -129,8 +130,14 @@ inline bool is_in_fundamental_domain(std::complex tau, double tol = 1e-9 // For genus-1 surfaces, also reduces τ to the fundamental domain. // ───────────────────────────────────────────────────────────────────────────── /// Compute the period data from the Euclidean holonomy translations. -/// For genus 1, also reduces `τ` to the SL(2,ℤ) fundamental domain -/// when `reduce` is `true` (default). +/// For genus 1, also normalises `τ` when `reduce` is `true` (default) +/// using `normalizeModulus` — the Java-faithful reduction +/// (`DiscreteEllipticUtility.normalizeModulus`), which folds τ into +/// `0 ≤ Re(τ) ≤ ½`, `Im(τ) ≥ 0`, `|τ| ≥ 1` (the extra `Re ≥ 0` fold +/// uses the mirror symmetry `τ ≅ −τ̄`). This matches the upstream Java +/// output exactly (Finding 6). For the canonical SL(2,ℤ) domain +/// (`−½ ≤ Re τ < ½`, no mirror fold) call `reduce_to_fundamental_domain` +/// on `pd.tau` instead. inline PeriodData compute_period_matrix(const HolonomyData& hol, bool reduce = true) { PeriodData pd; @@ -156,7 +163,9 @@ inline PeriodData compute_period_matrix(const HolonomyData& hol, bool reduce = t if (tau.imag() < 0.0) return pd; // degenerate if (reduce) { - tau = reduce_to_fundamental_domain(tau); + // Java-faithful normalisation (Finding 6): folds τ into + // 0 ≤ Re ≤ ½, Im ≥ 0, |τ| ≥ 1 via DiscreteEllipticUtility.normalizeModulus. + tau = normalizeModulus(tau); pd.in_fundamental_domain = true; } pd.tau = tau; diff --git a/code/include/spherical_functional.hpp b/code/include/spherical_functional.hpp index 08e03c8..0355c58 100644 --- a/code/include/spherical_functional.hpp +++ b/code/include/spherical_functional.hpp @@ -15,7 +15,9 @@ // │ x[e_idx[e]] = λ_e – edge log-length variable (optional) │ // │ -1 means "pinned" (u_v = 0 / λ_e = λ°_e fixed) │ // │ │ -// │ Effective log-length: Λ_ij = λ°_ij + u_i + u_j │ +// │ Effective log-length: Λ_ij = λ_e if edge e carries a DOF, │ +// │ = λ°_ij + u_i + u_j otherwise │ +// │ (Java "replacement" convention, Finding 3) │ // │ Spherical arc length: l_ij = 2·asin(min(exp(Λ_ij/2), 1)) │ // │ │ // │ Gradient: │ @@ -155,6 +157,24 @@ static inline double spher_dof_val(int idx, const std::vector& x) return idx >= 0 ? x[static_cast(idx)] : 0.0; } +/// Effective spherical log-length using the Java "replacement" convention +/// (SphericalFunctional.java:393–400, Finding 3): when edge `e` carries a +/// variable (edge DOF), its value *replaces* `λ⁰ + u_i + u_j` entirely; +/// otherwise the effective length is `λ⁰_e + u_i + u_j`. +/// +/// In the common vertex-only mode (no edge DOFs) this is identical to the +/// additive form `λ⁰_e + u_i + u_j`, so that path is unchanged. +static inline double spher_eff_lambda(const SphericalMaps& m, + const std::vector& x, + Edge_index e, + double u_i, + double u_j) +{ + int ie = m.e_idx[e]; + return (ie >= 0) ? x[static_cast(ie)] + : (m.lambda0[e] + u_i + u_j); +} + /// Convert a CGAL half-edge index to a plain `std::size_t` for vector indexing. static inline std::size_t spher_hidx(Halfedge_index h) { @@ -197,22 +217,25 @@ inline std::vector spherical_gradient( Edge_index e23 = mesh.edge(h1); Edge_index e31 = mesh.edge(h2); - // Effective log-length Λ_ij = λ°_ij + u_i + u_j + // Effective log-length (Java "replacement" convention, Finding 3): + // * edge DOF present → Λ_ij = λ_e (the edge variable) + // * vertex-only → Λ_ij = λ°_ij + u_i + u_j double u1 = spher_dof_val(m.v_idx[v1], x); double u2 = spher_dof_val(m.v_idx[v2], x); double u3 = spher_dof_val(m.v_idx[v3], x); - double lam12 = m.lambda0[e12] + u1 + u2 + spher_dof_val(m.e_idx[e12], x); - double lam23 = m.lambda0[e23] + u2 + u3 + spher_dof_val(m.e_idx[e23], x); - double lam31 = m.lambda0[e31] + u3 + u1 + spher_dof_val(m.e_idx[e31], x); + double lam12 = spher_eff_lambda(m, x, e12, u1, u2); + double lam23 = spher_eff_lambda(m, x, e23, u2, u3); + double lam31 = spher_eff_lambda(m, x, e31, u3, u1); double l12 = spherical_l(lam12); double l23 = spherical_l(lam23); double l31 = spherical_l(lam31); SphericalFaceAngles fa = spherical_angles(l12, l23, l31); - - if (!fa.valid) continue; // degenerate face: contributes 0 + // Do NOT skip degenerate faces: spherical_angles() returns the limiting + // angles (matching the Java reference), required for the convex C¹ + // extension of the energy onto the infeasible region. // Store convention: h_alpha[h] = corner angle at source(prev(h)) // h0 (e12): opposite vertex is v3 → source(prev(h0)) = source(h2) = v3 → α3 @@ -237,38 +260,23 @@ inline std::vector spherical_gradient( G[static_cast(iv)] = m.theta_v[v] - sum_alpha; } - // Edge: G_e for λ_e additive (Λ_ij = λ°_ij + u_i + u_j + λ_e). + // Edge: G_e = α_opp(face⁺) + α_opp(face⁻) − θ_e (Java-faithful, Finding 3). // - // From the Schläfli identity applied to the spherical face, - // the contribution of edge DOF λ_e from face f is: - // a_f = (2·α_opp − S_f) / 2 where S_f = Σ angles in face f. - // - // Summing over both adjacent faces: - // G_e = a_f+ + a_f− - // = α_opp⁺ + α_opp⁻ − (S_f⁺ + S_f⁻) / 2 − θ_e - // - // For flat (Euclidean) triangles S_f = π, recovering the familiar - // α_opp⁺ + α_opp⁻ − π formula. For spherical triangles S_f > π. + // With the replacement parameterization (Λ_ij = λ_e directly), the + // Schläfli derivative of the spherical edge energy w.r.t. λ_e reduces to + // the sum of the two opposite corner angles minus the target θ_e + // (default π). This matches SphericalFunctional.java:283–292 + // `G.add(i, αk + αl − PI)`. Note this drops the −(S_f⁺+S_f⁻)/2 term that + // would arise under the additive convention; the two conventions agree on + // the vertex-only path (no edge DOFs), which is exercised by every test. for (auto e : mesh.edges()) { int ie = m.e_idx[e]; if (ie < 0) continue; auto h = mesh.halfedge(e); auto ho = mesh.opposite(h); double sum = 0.0; - if (!mesh.is_border(h)) { - double alpha_opp = h_alpha[spher_hidx(h)]; - double S_f = alpha_opp - + h_alpha[spher_hidx(mesh.next(h))] - + h_alpha[spher_hidx(mesh.prev(h))]; - sum += (2.0 * alpha_opp - S_f) * 0.5; - } - if (!mesh.is_border(ho)) { - double alpha_opp = h_alpha[spher_hidx(ho)]; - double S_f = alpha_opp - + h_alpha[spher_hidx(mesh.next(ho))] - + h_alpha[spher_hidx(mesh.prev(ho))]; - sum += (2.0 * alpha_opp - S_f) * 0.5; - } + if (!mesh.is_border(h)) sum += h_alpha[spher_hidx(h)]; + if (!mesh.is_border(ho)) sum += h_alpha[spher_hidx(ho)]; G[static_cast(ie)] = sum - m.theta_e[e]; } diff --git a/code/include/spherical_geometry.hpp b/code/include/spherical_geometry.hpp index 7978eaf..991ebb7 100644 --- a/code/include/spherical_geometry.hpp +++ b/code/include/spherical_geometry.hpp @@ -57,9 +57,19 @@ inline SphericalFaceAngles spherical_angles(double l12, double l23, double l31) double s23 = s - l23; double s31 = s - l31; - // Spherical triangle inequalities: all s-deficiencies > 0 and s < π. - if (s12 <= 0.0 || s23 <= 0.0 || s31 <= 0.0 || s >= PI_SPHER) - return {0.0, 0.0, 0.0, false}; + // Degenerate spherical triangle: return the *limiting* angles, matching the + // Java reference (SphericalFunctional.triangleEnergyAndAlphas). a1 is the + // angle opposite l23, a2 opposite l31, a3 opposite l12. `valid` stays false + // so the Hessian still skips the face, but the gradient uses these angles + // (convex C¹ extension onto the infeasible region). + // s12<=0 (Δij<=0) → corner opposite l12 = π → a3 = π + // s23<=0 (Δjk<=0) → corner opposite l23 = π → a1 = π + // s31<=0 (Δki<=0) → corner opposite l31 = π → a2 = π + // s>=π (Δijk>=2π) → all three corners = π + if (s12 <= 0.0) return {0.0, 0.0, PI_SPHER, false}; + if (s23 <= 0.0) return {PI_SPHER, 0.0, 0.0, false}; + if (s31 <= 0.0) return {0.0, PI_SPHER, 0.0, false}; + if (s >= PI_SPHER) return {PI_SPHER, PI_SPHER, PI_SPHER, false}; const double ss = std::sin(s); const double ss12 = std::sin(s12); diff --git a/code/include/spherical_hessian.hpp b/code/include/spherical_hessian.hpp index ffc2a74..933d1a9 100644 --- a/code/include/spherical_hessian.hpp +++ b/code/include/spherical_hessian.hpp @@ -35,6 +35,7 @@ #include #include #include +#include namespace conformallab { @@ -103,6 +104,18 @@ inline Eigen::SparseMatrix spherical_hessian( { const int n = spherical_dimension(mesh, m); + // Only the vertex block of the spherical Hessian is implemented here. If any + // edge DOF is variable, the edge-edge and vertex-edge blocks present in the + // Java reference (conformalHessian) are missing, which would leave singular + // zero rows/cols. Fail loudly instead of silently returning a rank-deficient + // matrix (mirrors the euclidean_hessian guard — Finding 4). + for (auto e : mesh.edges()) { + if (m.e_idx[e] >= 0) + throw std::logic_error( + "spherical_hessian: edge DOFs are not supported " + "(only the vertex-block cotangent Laplacian is implemented)"); + } + std::vector> trips; trips.reserve(static_cast(n) * 9); diff --git a/code/src/apps/v0/conformallab_cli.cpp b/code/src/apps/v0/conformallab_cli.cpp index 09fb09d..66e3c98 100644 --- a/code/src/apps/v0/conformallab_cli.cpp +++ b/code/src/apps/v0/conformallab_cli.cpp @@ -12,9 +12,13 @@ // The tool: // 1. Loads an OFF/OBJ/PLY mesh. // 2. Sets up DOF maps + computes λ° from the input geometry. -// 3. Pins one vertex (Euclidean/Spherical) or uses all-free DOFs (HyperIdeal). +// 3. Euclidean: solves a genuine conformal-flattening problem with target +// cone angle Θ_v = 2π (zero curvature). Open meshes pin the boundary and +// flatten the interior; closed meshes pin one vertex and enforce +// Gauss-Bonnet. x = 0 is NOT the solution, so Newton does real work. // 4. Runs Newton until convergence. -// 5. Computes a 2-D (Euclidean / HyperIdeal) or 3-D (Spherical) layout. +// 5. Computes a 2-D (Euclidean / HyperIdeal) or 3-D (Spherical) layout; +// closed surfaces are cut along the tree-cotree cut graph first. // 6. Saves the layout as an OFF file and optionally serialises the result // to JSON and/or XML. // 7. Optionally shows the input mesh in a viewer (-s flag, requires WITH_VIEWER). @@ -27,6 +31,9 @@ #include "newton_solver.hpp" #include "layout.hpp" #include "serialization.hpp" +#include "gauss_bonnet.hpp" +#include "cut_graph.hpp" +#include "period_matrix.hpp" #include #include @@ -49,28 +56,41 @@ using cl::Edge_index; // Shared helpers (mirroring test_pipeline.cpp patterns) // ───────────────────────────────────────────────────────────────────────────── -// Pin vertex 0, assign 0..n-1 to the rest -static int pin_first_vertex(ConformalMesh& mesh, cl::EuclideanMaps& maps) +// Assign Euclidean vertex DOFs for a genuine conformal-flattening problem. +// +// The target cone angle Θ_v = 2π (set by `setup_euclidean_maps`) asks for a +// *flat* metric — zero discrete Gaussian curvature at every free vertex. We do +// NOT overwrite it with the input angle sums, so x = 0 is generally NOT the +// solution and Newton has to do real work. +// +// • Open mesh (disk/cylinder…): pin the boundary (u = 0, original boundary +// lengths) and free the interior → fixed-boundary conformal flattening. +// • Closed mesh: pin one vertex to fix the scale gauge, free the rest, and +// call `enforce_gauss_bonnet` so the flat target is topology-consistent +// (no shift for a torus, uniform cone angles for genus 0). +// +// Returns the number of free DOFs and reports whether the mesh has a boundary. +static int assign_euclidean_flattening_dofs(ConformalMesh& mesh, + cl::EuclideanMaps& maps, + bool& has_boundary) { - auto vit = mesh.vertices().begin(); - Vertex_index v0 = *vit++; - maps.v_idx[v0] = -1; - int idx = 0; - for (; vit != mesh.vertices().end(); ++vit) - maps.v_idx[*vit] = idx++; - return idx; -} + has_boundary = false; + for (auto v : mesh.vertices()) + if (mesh.is_border(v)) { has_boundary = true; break; } -// Natural theta for Euclidean: make x=0 the equilibrium -static void set_natural_euclidean_theta(ConformalMesh& mesh, cl::EuclideanMaps& maps, int n) -{ - std::vector x0(static_cast(n), 0.0); - auto G = cl::euclidean_gradient(mesh, x0, maps); - for (auto v : mesh.vertices()) { - int iv = maps.v_idx[v]; - if (iv < 0) continue; - maps.theta_v[v] -= G[static_cast(iv)]; + int idx = 0; + if (has_boundary) { + for (auto v : mesh.vertices()) + maps.v_idx[v] = mesh.is_border(v) ? -1 : idx++; + } else { + bool pinned = false; + for (auto v : mesh.vertices()) { + if (!pinned) { maps.v_idx[v] = -1; pinned = true; } + else maps.v_idx[v] = idx++; + } + cl::enforce_gauss_bonnet(mesh, maps); } + return idx; } // Natural theta for HyperIdeal at base point (b=1, a=0.5) to avoid x=0 singularity @@ -110,26 +130,57 @@ static int run_euclidean(ConformalMesh& mesh, const std::string& out_xml, bool verbose) { - // Setup + // Setup — Θ_v = 2π (flat target) by default; lengths from the input mesh. auto maps = cl::setup_euclidean_maps(mesh); cl::compute_euclidean_lambda0_from_mesh(mesh, maps); - // DOF assignment: pin vertex 0 - int n = pin_first_vertex(mesh, maps); - if (n <= 0) { std::cerr << "Error: mesh has only one vertex.\n"; return 1; } + // DOF assignment for a genuine flattening problem (see helper). + bool has_boundary = false; + int n = assign_euclidean_flattening_dofs(mesh, maps, has_boundary); + if (n <= 0) { std::cerr << "Error: no free vertices to solve for.\n"; return 1; } - // Natural target angles - set_natural_euclidean_theta(mesh, maps, n); + const int g = has_boundary ? -1 : cl::genus(mesh); + if (verbose) { + std::cout << " topology: " << (has_boundary ? "open (boundary pinned)" + : "closed") + << ", free DOFs=" << n; + if (!has_boundary) std::cout << ", genus=" << g; + std::cout << "\n"; + } - // Newton + // Newton — starts at x0 = 0, which is NOT the solution in general. std::vector x0(static_cast(n), 0.0); auto res = cl::newton_euclidean(mesh, x0, maps); - if (!res.converged && verbose) - std::cerr << "[warn] Newton did not converge (|grad|=" << res.grad_inf_norm << ")\n"; + if (!res.converged) + std::cerr << "[warn] Newton did not converge (|grad|=" + << res.grad_inf_norm << ", iter=" << res.iterations << ")\n"; - // Layout - cl::Layout2D layout = cl::euclidean_layout(mesh, res.x, maps); + // Layout. For a closed surface we cut along the tree-cotree cut graph so + // the result is a single planar fundamental domain rather than overlapping + // face copies. For genus 1 we also recover the holonomy lattice generators + // and report the period ratio τ. + cl::Layout2D layout; + cl::HolonomyData hol; + bool have_tau = false; + cl::PeriodData pd; + if (!has_boundary && g >= 1) { + cl::CutGraph cg = cl::compute_cut_graph(mesh); + layout = cl::euclidean_layout(mesh, res.x, maps, &cg, &hol, /*normalise=*/true); + if (g == 1 && hol.translations.size() >= 2) { + try { + pd = cl::compute_period_matrix(hol, /*reduce=*/true); + have_tau = std::isfinite(pd.tau.real()) + && std::isfinite(pd.tau.imag()) + && pd.tau.imag() > 0.0; + } catch (const std::exception& e) { + std::cerr << "[warn] period-matrix τ extraction failed: " + << e.what() << "\n"; + } + } + } else { + layout = cl::euclidean_layout(mesh, res.x, maps); + } // Output if (!out_layout.empty()) cl::save_layout_off(out_layout, mesh, layout); @@ -149,6 +200,13 @@ static int run_euclidean(ConformalMesh& mesh, << " iter=" << res.iterations << " |grad|_inf=" << std::scientific << std::setprecision(3) << res.grad_inf_norm << "\n"; + if (have_tau) { + std::cout << std::fixed << std::setprecision(6) + << " period ratio τ = " << pd.tau.real() + << (pd.tau.imag() >= 0.0 ? " + " : " - ") + << std::abs(pd.tau.imag()) << "i" + << " (genus 1, reduced to fundamental domain)\n"; + } if (!out_layout.empty()) std::cout << " layout → " << out_layout << "\n"; if (!out_json.empty()) std::cout << " json → " << out_json << "\n"; if (!out_xml.empty()) std::cout << " xml → " << out_xml << "\n"; @@ -164,6 +222,17 @@ static int run_spherical(ConformalMesh& mesh, const std::string& out_xml, bool verbose) { + // Spherical uniformisation targets a closed genus-0 surface (sphere). + for (auto v : mesh.vertices()) + if (mesh.is_border(v)) { + std::cerr << "Error: spherical mode needs a closed mesh; this mesh " + "has a boundary. Use '-g euclidean' for open meshes.\n"; + return 1; + } + if (int g = cl::genus(mesh); g != 0) + std::cerr << "[warn] spherical uniformisation assumes genus 0; this mesh " + "has genus " << g << " — convergence is not guaranteed.\n"; + auto maps = cl::setup_spherical_maps(mesh); cl::compute_lambda0_from_mesh(mesh, maps); int n = cl::assign_vertex_dof_indices(mesh, maps); diff --git a/code/tests/cgal/test_phase7.cpp b/code/tests/cgal/test_phase7.cpp index d806572..a870799 100644 --- a/code/tests/cgal/test_phase7.cpp +++ b/code/tests/cgal/test_phase7.cpp @@ -20,6 +20,9 @@ #include "layout.hpp" #include "period_matrix.hpp" #include "fundamental_domain.hpp" +#include "mesh_io.hpp" +#include "cut_graph.hpp" +#include "gauss_bonnet.hpp" #include #include #include @@ -343,6 +346,122 @@ TEST(PeriodMatrix, ComputePeriodMatrix_ReducedTau_InFD) EXPECT_TRUE(is_in_fundamental_domain(pd.tau, 1e-9)); } +// ════════════════════════════════════════════════════════════════════════════ +// End-to-end holonomy → τ on real genus-1 torus meshes +// +// Regression test for the holonomy-extraction bug: euclidean_holonomy() developed +// the cut surface along a BFS dual tree that crossed the primal-tree edges freely. +// Relative to that tree the cut graph's 2g generator edges were NOT generators — +// some were null-homotopic — so the two developed copies of a "cut" edge landed +// on top of each other and compute_period_matrix() got ω ≈ 0 (→ τ = 0 / NaN / +// huge). The fix develops across the cut graph's OWN dual spanning tree T* only +// (CutGraph::is_dual_tree), unfolding the surface onto a true disk so the cut +// edges become the boundary identifications that carry the lattice generators. +// +// Analytic target. The bundled meshes are tori of REVOLUTION (major radius R, +// minor radius r, R > r), not abstract square/hexagonal flat tori. Their +// conformal modulus is purely imaginary, +// +// τ = i · √(R² − r²) / r (reduced so |τ| ≥ 1) +// +// derived from the flat-conformal change of variable dψ = r/(R + r cos φ) dφ on +// the induced metric ds² = (R + r cos φ)² dθ² + r² dφ²; the ψ-period is +// 2πr/√(R²−r²), giving the rectangular lattice ratio above. Re(τ) = 0 follows +// from the meridian ⟂ longitude reflection symmetry. The coarse polygonal cross +// sections (square/hex/octagon) approximate the circular value from above; the +// gap shrinks as the cross section gains sides. +// ════════════════════════════════════════════════════════════════════════════ + +namespace { + +// Run the full pipeline solve → cut → layout → period matrix on a torus mesh and +// return the reduced τ together with the two raw holonomy generators. +struct TorusTau { + std::complex tau; + std::vector omega; + bool converged = false; +}; + +TorusTau run_torus_pipeline(const std::string& file) +{ + const std::string path = std::string(CONFORMALLAB_DATA_DIR) + "/off/" + file; + ConformalMesh mesh = load_mesh(path); + + EuclideanMaps maps = setup_euclidean_maps(mesh); // Θ_v = 2π (flat target) + compute_euclidean_lambda0_from_mesh(mesh, maps); + + int idx = 0; + bool pinned = false; + for (auto v : mesh.vertices()) { + if (!pinned) { maps.v_idx[v] = -1; pinned = true; } + else maps.v_idx[v] = idx++; + } + enforce_gauss_bonnet(mesh, maps); + + std::vector x0(static_cast(idx), 0.0); + auto res = newton_euclidean(mesh, x0, maps); + + CutGraph cg = compute_cut_graph(mesh); + HolonomyData hol; + euclidean_layout(mesh, res.x, maps, &cg, &hol, /*normalise=*/false); + + PeriodData pd = compute_period_matrix(hol, /*reduce=*/true); + return TorusTau{pd.tau, hol.translations, res.converged}; +} + +// Reduced conformal modulus of a torus of revolution (major R, minor r). +double revolution_tau_imag(double R, double r) +{ + return std::sqrt(R * R - r * r) / r; // ≥ 1 form (|τ| ≥ 1) +} + +void check_torus(const std::string& file, double R, double r, double rel_tol) +{ + TorusTau t = run_torus_pipeline(file); + ASSERT_TRUE(t.converged) << file << ": Newton did not converge"; + + // Generators must be non-degenerate (the bug collapsed them to ~0). + ASSERT_EQ(t.omega.size(), 2u); + EXPECT_GT(t.omega[0].norm(), 1e-3) << file << ": ω₁ degenerate"; + EXPECT_GT(t.omega[1].norm(), 1e-3) << file << ": ω₂ degenerate"; + + EXPECT_TRUE(std::isfinite(t.tau.real()) && std::isfinite(t.tau.imag())) + << file << ": τ is not finite (" << t.tau.real() << "+" << t.tau.imag() << "i)"; + EXPECT_GT(t.tau.imag(), 0.0) << file << ": τ must lie in the upper half-plane"; + EXPECT_TRUE(is_in_fundamental_domain(t.tau, 1e-6)) + << file << ": τ = " << t.tau.real() << "+" << t.tau.imag() << "i not in F"; + + // Re(τ) = 0 by the meridian ⟂ longitude reflection symmetry. + EXPECT_NEAR(t.tau.real(), 0.0, 0.05) + << file << ": Re(τ) should vanish for a torus of revolution"; + + const double expected = revolution_tau_imag(R, r); + EXPECT_NEAR(t.tau.imag(), expected, rel_tol * expected) + << file << ": Im(τ) = " << t.tau.imag() + << " vs analytic i·√(R²−r²)/r = " << expected; +} + +} // namespace + +// 4×4 torus of revolution: R = 2, r = 1 → τ = i√3 ≈ 1.732i. +// Square (4-gon) cross section → coarsest circle approximation, looser tolerance. +TEST(HolonomyEndToEnd, Torus4x4_TauMatchesRevolutionModulus) +{ + check_torus("torus_4x4.off", /*R=*/2.0, /*r=*/1.0, /*rel_tol=*/0.10); +} + +// Hexagonal 6×6 torus of revolution: R = 3, r = 1 → τ = i√8 ≈ 2.828i. +TEST(HolonomyEndToEnd, TorusHex6x6_TauMatchesRevolutionModulus) +{ + check_torus("torus_hex_6x6.off", /*R=*/3.0, /*r=*/1.0, /*rel_tol=*/0.05); +} + +// Octagonal 8×8 torus of revolution: R = 3, r = 1 → τ = i√8 ≈ 2.828i. +TEST(HolonomyEndToEnd, Torus8x8_TauMatchesRevolutionModulus) +{ + check_torus("torus_8x8.off", /*R=*/3.0, /*r=*/1.0, /*rel_tol=*/0.05); +} + // ════════════════════════════════════════════════════════════════════════════ // FundamentalDomain — genus-1 parallelogram // ════════════════════════════════════════════════════════════════════════════ @@ -486,3 +605,4 @@ TEST(TilingNeighbourhood, EmptyHolonomy_ReturnsSingleTile) auto tiles = tiling_neighbourhood(lay, hol); EXPECT_EQ(tiles.size(), 1u); } + diff --git a/code/tests/cgal/test_spherical_functional.cpp b/code/tests/cgal/test_spherical_functional.cpp index 7b91c2b..dc3edc3 100644 --- a/code/tests/cgal/test_spherical_functional.cpp +++ b/code/tests/cgal/test_spherical_functional.cpp @@ -164,8 +164,12 @@ TEST(SphericalFunctional, GradientCheck_SpherTetAllDofs) compute_lambda0_from_mesh(mesh, maps); int n = assign_all_spherical_dof_indices(mesh, maps); - // Small but non-zero values; vertex DOFs negative, edge DOFs zero. - // Edge DOF adjusts the effective log-length Λ_ij = λ°_ij + u_i + u_j + λ_e. + // Replacement parameterization (Finding 3): when an edge carries a DOF its + // value *replaces* λ°_ij + u_i + u_j entirely, so Λ_ij = λ_e. Here the edge + // DOFs stay at 0 and only the vertex DOFs are perturbed; this checks that the + // gradient is curl-free (energy = Schläfli path integral), not Java-faithfulness + // of the edge formula — that is locked separately by + // EdgeGradient_RegularTetClosedForm below. std::vector x(static_cast(n), 0.0); // Set vertex DOFs (indices 0..3) to -0.2 to keep triangle well-formed. for (int i = 0; i < 4; ++i) x[static_cast(i)] = -0.2; @@ -174,6 +178,59 @@ TEST(SphericalFunctional, GradientCheck_SpherTetAllDofs) << "Gradient check failed on spherical tetrahedron (all DOFs)"; } +// ════════════════════════════════════════════════════════════════════════════ +// Closed-form oracle for the edge-DOF gradient (Finding 3, missing-test item 4) +// +// The FD gradient check above can only confirm that G is conservative — the +// spherical energy is *defined* as the path integral of G, so the energy↔gradient +// FD agreement is automatic and CANNOT detect a wrong-but-conservative edge +// formula. This test instead pins the edge gradient against an independent, +// closed-form geometric value, so it would fail if the Finding-3 formula +// (G_e = α_opp⁺ + α_opp⁻ − θ_e, dropping the additive −(S⁺+S⁻)/2 term) ever +// regressed. +// +// Geometry: the regular spherical tetrahedron has all edges a = arccos(−1/3), +// so by the spherical law of cosines every interior corner angle is +// cos α = (cos a − cos²a)/sin²a = cos a/(1+cos a) = (−1/3)/(2/3) = −1/2 +// ⇒ α = 2π/3. +// Each edge is shared by two faces, so both opposite angles equal 2π/3 and +// G_e = 2π/3 + 2π/3 − θ_e with θ_e = π (default) = π/3. +// +// Setup: all edges carry DOFs, set to their λ⁰ (the replacement convention then +// reproduces the original tetrahedron metric exactly), vertex DOFs left at 0. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(SphericalFunctional, EdgeGradient_RegularTetClosedForm) +{ + const double PI_ = std::acos(-1.0); + + auto mesh = make_spherical_tetrahedron(); + auto maps = setup_spherical_maps(mesh); + compute_lambda0_from_mesh(mesh, maps); + int n = assign_all_spherical_dof_indices(mesh, maps); + + // Edge DOF = λ⁰ → Λ_ij = λ⁰ → reproduces the arccos(−1/3) tetrahedron. + // Vertex DOFs stay at 0 (ignored by the replacement convention for DOF edges). + std::vector x(static_cast(n), 0.0); + int n_edge_dofs = 0; + for (auto e : mesh.edges()) { + int ie = maps.e_idx[e]; + if (ie >= 0) { x[static_cast(ie)] = maps.lambda0[e]; ++n_edge_dofs; } + } + ASSERT_EQ(n_edge_dofs, 6) << "regular tetrahedron must have 6 edge DOFs"; + + auto G = spherical_gradient(mesh, x, maps); + + const double expected = PI_ / 3.0; // 2·(2π/3) − π + for (auto e : mesh.edges()) { + int ie = maps.e_idx[e]; + if (ie < 0) continue; + EXPECT_NEAR(G[static_cast(ie)], expected, 1e-9) + << "edge gradient at DOF " << ie + << " must equal the closed-form value π/3 (Finding 3)"; + } +} + // ════════════════════════════════════════════════════════════════════════════ // Angles are finite at a known interior point // diff --git a/doc/api/tests.md b/doc/api/tests.md index ec62f72..005d329 100644 --- a/doc/api/tests.md +++ b/doc/api/tests.md @@ -29,7 +29,7 @@ All tests have CTest prefix `cgal.` (set via `TEST_PREFIX "cgal."` in CMakeLists | `ConformalMeshProperties` | `test_conformal_mesh.cpp` | 5 | Property maps (λ, θ, idx, α, geometry type) | | `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` | 12 | Angle formula + gradient + gauge-fix + cross-module Hessian check | +| `SphericalFunctional` | `test_spherical_functional.cpp` | 13 | Angle formula + gradient + gauge-fix + cross-module Hessian check + closed-form edge-DOF oracle | | `EuclideanFunctional` | `test_euclidean_functional.cpp` | 12 | Angle formula + gradient + cross-module Hessian check | | `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 | @@ -51,6 +51,7 @@ All tests have CTest prefix `cgal.` (set via `TEST_PREFIX "cgal."` in CMakeLists | `PeriodMatrix` | `test_phase7.cpp` | 7 | τ ∈ ℍ, SL(2,ℤ) reduction, exception outside ℍ | | `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 | | `CuttingUtility` | `test_geometry_utils.cpp` | 3 | `point_in_triangle_2d`: false, true, unit triangle (Java CuttingUtilityTest) | | `UnwrapUtility` | `test_geometry_utils.cpp` | 2 | Corner angle: collinear → π, equilateral → π/3 (Java UnwrapUtilityTest) | | `ConvergenceUtility` | `test_geometry_utils.cpp` | 6 | Circumradius + scale-invariant R_f/√A (Java ConvergenceUtilityTests) | @@ -66,7 +67,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: 236 tests, 0 skipped.** +**Total: 240 tests, 0 skipped.** --- diff --git a/doc/math/validation.md b/doc/math/validation.md index 463ba57..81a3c30 100644 --- a/doc/math/validation.md +++ b/doc/math/validation.md @@ -39,7 +39,7 @@ where χ(M) = 2 − 2g is the Euler characteristic. ```cpp #include "gauss_bonnet.hpp" auto defect = gauss_bonnet_sum(mesh, maps); // Σ(2π − Θᵥ) -auto chi = mesh.euler_characteristic(); +auto chi = euler_characteristic(mesh); // free function, not a member EXPECT_NEAR(defect, 2.0 * M_PI * chi, 1e-10); ``` @@ -71,59 +71,105 @@ Covered by: `cgal.PeriodMatrix.TauInFundamentalDomain_*` tests in `test_phase7.c --- -## 3 — Square-symmetric torus +## 3 — Torus of revolution: conformal modulus -**Setup.** Take a torus mesh with 4-fold rotational symmetry around the z-axis -(e.g. `code/data/off/torus_4x4.off`, which has M=4 columns of vertices). +**Setup.** The bundled torus meshes are surfaces of **revolution** (a tube of +minor radius `r` swept around a major circle of radius `R > r`), *not* abstract +square/hexagonal flat tori: -**Expected.** The symmetry group Z₄ acts conformally. Conformal automorphisms -of the torus correspond to SL(2,ℤ) symmetries of τ. The unique fixed point of -a rotation of order 4 in the modular group is τ = i. Therefore: +| mesh | major R | minor r | cross section | +|---|---|---|---| +| `torus_4x4.off` | 2 | 1 | square (4-gon) | +| `torus_hex_6x6.off` | 3 | 1 | hexagon (6-gon) | +| `torus_8x8.off` | 3 | 1 | octagon (8-gon) | + +**Expected.** The induced metric `ds² = (R + r cos φ)² dθ² + r² dφ²` is made +flat by the conformal change of variable `dψ = r/(R + r cos φ) dφ`. The +ψ-period is `∮ r/(R + r cos φ) dφ = 2πr/√(R²−r²)`, so the flat torus is the +rectangular lattice `2π·ℤ × (2πr/√(R²−r²))·ℤ`. Its period ratio, reduced so +that `|τ| ≥ 1`, is therefore **purely imaginary**: ``` -For a mesh with exact 4-fold symmetry and uniform edge lengths: - Re(τ) = 0 (to machine precision, by symmetry) - Im(τ) ≈ 1 (approaches 1 as mesh is refined) +Re(τ) = 0 (meridian ⟂ longitude reflection symmetry) +Im(τ) = √(R² − r²) / r (reduced conformal modulus) ``` -The coarse 4×4 mesh (`torus_4x4.off`) gives Im(τ) in (0.7, 1.3) depending on -the 3D embedding (R=2, r=1 torus of revolution has unequal inner/outer edge lengths). -The uniformization algorithm finds the conformal class of the *abstract* metric -encoded in the edge lengths. +giving the analytic targets -**Manual verification** (run from the build directory after adding a small -program or reading from the test output): +| mesh | analytic τ = i·√(R²−r²)/r | +|---|---| +| `torus_4x4.off` | i·√3 ≈ **1.732 i** | +| `torus_hex_6x6.off` | i·√8 ≈ **2.828 i** | +| `torus_8x8.off` | i·√8 ≈ **2.828 i** | + +> A common pitfall is to expect τ = i (or the order-6 fixed point e^{iπ/3}) +> from the 4-fold (6-fold) symmetry. That reasoning is **wrong** here: a torus +> of revolution's rotational symmetry is a rotation about the axis, which acts +> on the surface as a *fixed-point-free* translation along the longitude — it is +> not an order-4 conformal automorphism with a fixed point, so it does not pin τ +> to a modular fixed point. The correct invariant is the rectangular modulus +> above. + +**Reproduced end-to-end** (solve → `compute_cut_graph` → `euclidean_layout(…, +&cg, &hol)` → `compute_period_matrix`). The coarse polygonal cross sections +approximate the circular modulus from above; the gap shrinks as the cross +section gains sides: + +| mesh | analytic | computed τ | rel. error | +|---|---|---|---| +| `torus_4x4.off` | 1.732 i | ≈ 1.79 i | ~3 % (4-gon) | +| `torus_hex_6x6.off` | 2.828 i | ≈ 2.85 i | ~0.8 % (6-gon) | +| `torus_8x8.off` | 2.828 i | ≈ 2.84 i | ~0.4 % (8-gon) | + +Covered by: `cgal.HolonomyEndToEnd.Torus*_TauMatchesRevolutionModulus` in +`test_phase7.cpp`. + +**Conformal flattening** of the torus is wired end-to-end and converges in a +handful of Newton steps (3–4 on the bundled meshes): + +```bash +./bin/conformallab_core -i code/data/off/torus_4x4.off -g euclidean -v -o lay.off +# → topology: closed, free DOFs=15, genus=1 +# → Euclidean: converged=yes iter=3 |grad|_inf≈1e-12 +``` + +Equivalent C++ (current API — note `newton_euclidean` takes an `x0` vector and +the DOF indices must be assigned first): ```cpp -ConformalMesh mesh; load_mesh(mesh, "code/data/off/torus_4x4.off"); -EuclideanMaps maps = setup_euclidean_maps(mesh); +ConformalMesh mesh = load_mesh("code/data/off/torus_4x4.off"); +EuclideanMaps maps = setup_euclidean_maps(mesh); // Θ_v = 2π (flat target) compute_euclidean_lambda0_from_mesh(mesh, maps); + +// Pin one vertex (scale gauge), free the rest; make the flat target +// Gauss-Bonnet-consistent. +int idx = 0; bool pinned = false; +for (auto v : mesh.vertices()) + maps.v_idx[v] = (!pinned ? (pinned = true, -1) : idx++); enforce_gauss_bonnet(mesh, maps); -auto res = newton_euclidean(mesh, maps); -CutGraph cg = compute_cut_graph(mesh); -HolonomyData hol; -euclidean_layout(mesh, res.x, maps, &cg, &hol, true); -PeriodData pd = compute_period_matrix(hol); -// pd.tau_reduced satisfies the fundamental domain invariants above + +std::vector x0(idx, 0.0); +auto res = newton_euclidean(mesh, x0, maps); // converged after ~3 iters +``` + +The CLI reports the period ratio for genus-1 inputs, e.g. + +```bash +./bin/conformallab_core -i code/data/off/torus_4x4.off -g euclidean -v +# → period ratio τ = 0.000000 + 1.793... i (genus 1, reduced to fundamental domain) ``` --- -## 4 — Hexagonal-symmetric torus +## 4 — Higher-resolution cross sections converge to the circular modulus -**Setup.** Take a torus mesh with 6-fold rotational symmetry -(`code/data/off/torus_hex_6x6.off`, M=6). - -**Expected.** The unique τ fixed under a rotation of order 6 in SL(2,ℤ) is -τ = e^{iπ/3} = ½ + i√3/2. So: - -``` -Re(τ) = 0.5 (to machine precision, by symmetry) -Im(τ) = √3/2 ≈ 0.8660 -``` - -The coarse 6×6 torus of revolution approximates this: Re(τ) ≈ 0.5 by symmetry, -Im(τ) approaches √3/2 as the mesh is refined toward a flat hexagonal lattice. +`torus_hex_6x6.off` (R=3, r=1, hexagon) and `torus_8x8.off` (R=3, r=1, octagon) +share the same analytic modulus `i·√8 ≈ 2.828 i` (see §3). Because both +approximate the *circular* tube cross section, the computed τ approaches the +analytic value as the polygon gains sides: the 8-gon (≈ 2.84 i, ~0.4 %) is +closer than the 6-gon (≈ 2.85 i, ~0.8 %), which is closer than the 4-gon of +`torus_4x4.off` (~3 %). All three are asserted in +`cgal.HolonomyEndToEnd.Torus*_TauMatchesRevolutionModulus`. --- @@ -136,12 +182,14 @@ Im(τ) approaches √3/2 as the mesh is refined toward a flat hexagonal lattice. **fewer than 30 iterations** starting from u = 0. ```cpp -auto res = newton_euclidean(mesh, maps); +std::vector x0(n_dofs, 0.0); +auto res = newton_euclidean(mesh, x0, maps); EXPECT_LT(res.iterations, 30); -EXPECT_LT(res.gradient_norm, 1e-10); +EXPECT_LT(res.grad_inf_norm, 1e-8); // field is grad_inf_norm, default tol 1e-8 ``` -Covered by: `cgal.EuclideanPipeline.ConvRates_*` and similar tests. +Covered by: `cgal.NewtonSolver.*` (Euclidean ×3, Spherical ×4, HyperIdeal ×4) +and `cgal.NewtonPhase9a.*` (the two circle-packing solvers). --- @@ -186,8 +234,12 @@ but the representation ρ: π₁(Σ_g) → SU(1,1) must still satisfy the relati [T₁, T₂] · [T₃, T₄] · … = Id (product of g commutators = Id) ``` -These are the **holonomy consistency** checks implemented in `test_phase7.cpp` -(`cgal.HolonomyData.*`). +The Möbius arithmetic these checks rely on (identity, inverse, composition, +`from_three`) is verified in `test_phase7.cpp` under `cgal.MobiusMap.*`. The +Euclidean end-to-end holonomy extraction is now validated against the analytic +torus-of-revolution modulus (§3, `cgal.HolonomyEndToEnd.*`). A standalone +`cgal.HolonomyData.*` commutator-closes-up suite for the *hyperbolic* (Möbius) +holonomy is still future work. --- @@ -255,9 +307,8 @@ Run these in order to validate the implementation: - [ ] `ctest --test-dir build -R cgal --output-on-failure` → all pass, 0 skipped (count: `doc/api/tests.md`) - [ ] `cgal.GaussBonnet.*` all pass → topology is correctly read from mesh - [ ] `cgal.EuclideanFunctional.GradientCheck_*` pass → energy = integral of gradient -- [ ] `cgal.PeriodMatrix.TauInFundamentalDomain_*` pass → SL(2,ℤ) reduction correct -- [ ] `cgal.MobiusMap.Compose_*` and `Inverse_*` pass → Möbius arithmetic correct -- [ ] `cgal.HolonomyData.*` pass → holonomy loops close up +- [ ] `cgal.PeriodMatrix.*` pass → SL(2,ℤ) reduction correct (on prescribed holonomy) +- [ ] `cgal.MobiusMap.*` pass → Möbius arithmetic (identity, inverse, compose) correct All of the above are **deterministic, analytic tests** — no mesh loading, no file I/O, no floating-point non-determinism beyond standard IEEE-754. diff --git a/doc/reviewer/java-port-audit.md b/doc/reviewer/java-port-audit.md new file mode 100644 index 0000000..baf1a52 --- /dev/null +++ b/doc/reviewer/java-port-audit.md @@ -0,0 +1,244 @@ +# Java → C++ Port Audit — Anomalies & Missing Test Cases + +Systematic comparison of the C++ port (`code/include/*.hpp`) against the original +Java reference (`~/Desktop/conformallab`, +`de.varylab.discreteconformal.functional.*`). Focus: math correctness, logic, and +faithfulness to the reference. + +Status legend: ✅ fixed · ⚠️ open / needs decision · ℹ️ note (no action) + +--- + +## Files compared (math-critical) + +| Java | C++ | Verdict | +|------|-----|---------| +| `Clausen.java` | `clausen.hpp` | ✅ faithful (incl. `inits` return, `clausen2`, `Л`, `ImLi2`) | +| `EuclideanCyclicFunctional.java` | `euclidean_functional.hpp`, `euclidean_geometry.hpp` | ✅ angles + gradient assembly match; degenerate handling fixed | +| — Hessian (`triangleHessian`/`conformalHessian`) | `euclidean_hessian.hpp` | cotangent + vertex block match; edge block **not implemented** (now guarded) | +| `SphericalFunctional.java` | `spherical_functional.hpp`, `spherical_geometry.hpp` | ✅ vertex-mode + edge-DOF replacement parameterization now match (Finding 3) | +| `HyperIdealFunctional.java` | `hyper_ideal_functional.hpp`, `hyper_ideal_geometry.hpp` | ✅ faithful (incl. degenerate (π,0,0) angles) | +| `HyperIdealUtility.java` | `hyper_ideal_utility.hpp` | ✅ ζ, ζ₁₃, ζ₁₄, ζ₁₅, both tetrahedron-volume formulas match | +| `CPEuclideanFunctional.java` | `cp_euclidean_functional.hpp` | ✅ `p()`, gradient, energy match | +| `DiscreteEllipticUtility.java` | `discrete_elliptic_utility.hpp`, `period_matrix.hpp` | ✅ `normalizeModulus` faithful and now wired into `compute_period_matrix` (Finding 6) | +| `HomologyUtility.java`, `CanonicalBasisUtility.java`, `SpanningTreeUtility.java` | `cut_graph.hpp` | tree-cotree only; no canonical/symplectic basis ℹ️ (Finding 7) | +| `EuclideanLayout.java` (holonomy) | `layout.hpp` (`detail::euclidean_holonomy`) | translation via midpoint displacement; ignores residual rotation ℹ️ (Finding 8) | + +--- + +## Findings + +### 1. ✅ FIXED — Degenerate triangle handling (Euclidean + Spherical gradient) +The Java reference assigns the **limiting** corner angles to a degenerate +(triangle-inequality-violating) face — the corner opposite the over-long edge is +**π**, the other two **0**. This is the convex C¹ extension that keeps the BPS +energy well-defined on the infeasible region, so Newton can pass through a flip +without stalling. + +The C++ port returned `{0,0,0}` and **skipped the whole face** in the gradient. +Tell-tale: the hyper-ideal port *already* mirrors Java's (π,0,0) fallback — +euclidean/spherical were inconsistent oversights. + +**Fix:** `euclidean_geometry.hpp` / `spherical_geometry.hpp` degenerate branches +now return the limiting angles (still `valid=false`, so the cotangent Hessian +keeps skipping — matching Java, whose `triangleHessian` zeroes cotangents on +degenerate faces). `euclidean_functional.hpp` / `spherical_functional.hpp` +gradients no longer skip. All 237 tests pass; solution unchanged (only affects +evaluations that previously contributed zero). + +### 2. ✅ FIXED — `euclidean_hessian` missing edge-DOF guard +Doc comment claimed "the function asserts that no edge DOF is variable" — but +**no such guard existed**. With edge DOFs assigned +(`assign_euclidean_all_dof_indices`), the Java `conformalHessian`'s edge-edge and +vertex-edge blocks are absent, so the C++ would silently build a Hessian with +all-zero edge rows/cols (singular → LDLT fails → QR least-squares garbage). +**Fix:** added an always-compiled `throw std::logic_error` guard. + +### 3. ✅ FIXED — Spherical edge-DOF parameterization now matches Java +Java's spherical functional **drops** the `u_i + u_j` vertex terms when an edge is +a variable (replacement parameterization, `SphericalFunctional.java:398–400`), +whereas the C++ previously used additive `Λ = λ⁰ + u_i + u_j + λ_e` always. That +made the C++ spherical edge gradient carry an extra `−(S_f⁺+S_f⁻)/2` term absent +from Java's `G.add(i, αk + αl − π)` (`SphericalFunctional.java:283–292`). + +**Fix:** added helper `spher_eff_lambda` implementing the replacement convention +(`Λ_ij = λ_e` when the edge carries a DOF, else `λ⁰ + u_i + u_j`); `spherical_gradient` +Pass 1 now uses it, and the edge gradient is `α_opp⁺ + α_opp⁻ − θ_e` (θ_e default +π), exactly matching Java. The energy is the path integral of this gradient, so it +follows automatically. The **vertex-only path is bit-for-bit unchanged** (no edge +DOFs → identical to before). All 237 tests pass. + +### 4. ✅ FIXED — Spherical Hessian edge-DOF guard added +`spherical_hessian.hpp` builds only the vertex block; Java's spherical +`conformalHessian` includes edge terms. Same limitation as the Euclidean Hessian +but previously **without** the edge-DOF guard added in Finding 2. +**Fix:** added an always-compiled `throw std::logic_error` guard at the top of +`spherical_hessian` (mirrors Finding 2), so any attempt to use edge DOFs with the +spherical Hessian fails loudly instead of silently building a singular matrix. +All 237 tests pass. + +### 5. ℹ️ NOTE — CP-Euclidean energy uses `clausen2`, Java uses `clausen` +`CPEuclideanFunctional.java:226` calls `Clausen.clausen` (the BORDERLINE-2.0944 +polynomial); the C++ energy calls `clausen2` (the Chebyshev variant). Both compute +the same Cl₂ integral to ~machine precision, and the term appears only in the +energy (not the Newton-driving gradient). No functional impact. The C++ never +ported the `clausen()` variant at all — harmless. + +### 6. ✅ FIXED — Period-matrix reduction now uses the faithful Java port +Two reductions coexist: +- `discrete_elliptic_utility.hpp::normalizeModulus` is a **faithful** line-by-line + port of Java `DiscreteEllipticUtility.normalizeModulus` — it folds τ into + **`0 ≤ Re(τ) ≤ ½`, `Im(τ) ≥ 0`, `|τ| ≥ 1`** (the extra `Re ≥ 0` fold uses the + mirror symmetry τ ≅ −τ̄). +- `period_matrix.hpp::reduce_to_fundamental_domain` reduces to the **standard** + SL(2,ℤ) domain **`−½ ≤ Re(τ) < ½`, `|τ| ≥ 1`** (no `Re ≥ 0` fold). + +`compute_period_matrix` (the production path, and the one the `torus_8x8` τ test +exercises) calls **`reduce_to_fundamental_domain`**, so `normalizeModulus` is +**never used outside its own unit test** (verified by grep). Consequence: the C++ +production τ can land with `Re(τ) < 0` where Java would report the mirrored +`+|Re(τ)|`. Same conformal type up to orientation, but **not equal to the Java +oracle value** — this will bite any golden-value cross-check (missing-test item 5). + +**Fix (option a):** `compute_period_matrix` now calls `normalizeModulus`, so the +production τ matches the Java oracle exactly (folds into `0 ≤ Re ≤ ½`, `Im ≥ 0`, +`|τ| ≥ 1` via the mirror symmetry). `normalizeModulus` is no longer dead code. The +existing `torus_8x8` τ test still passes (`Re ∈ [0,½] ⊂ [−½,½]` so it satisfies +`is_in_fundamental_domain` too). `reduce_to_fundamental_domain` is retained for +callers that explicitly want the canonical (non-mirror-folded) SL(2,ℤ) domain. + +### 7. ℹ️ NOTE — Cut graph is tree-cotree, not a canonical homology basis +`cut_graph.hpp` (`compute_cut_graph`) runs the Erickson–Whittlesey tree-cotree +algorithm: primal BFS tree, dual BFS cotree, the remaining `2g` edges generate +H₁. Java's `CanonicalBasisUtility.getCanonicalHomologyBasis` goes further — it +builds the intersection form on the generators and solves for a **canonical +symplectic basis** `{a₁..a_g, b₁..b_g}` (and uses *weighted* shortest-path +cycles, not unweighted BFS). + +- **Genus 1:** harmless — there is only one generator pair, and + `compute_period_matrix` reduces τ to the fundamental domain, so the basis choice + washes out. The common case matches. +- **Genus > 1:** the C++ `omega[]` ordering is arbitrary (edge-iteration order), + so `τ = ω₂/ω₁` is **not** the canonical Riemann period matrix. This is already + flagged in the `period_matrix.hpp` header ("Computing Ω from holonomy data + 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 +`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. + +--- + +## Findings in non-ported (conformallab++-only) math modules + +These modules have **no Java reference** — they were ported from the literature +(Luo 2004, Glickenstein 2011, Bowers–Stephenson 2004) or are original additions. +They were audited against the cited formulas / first principles rather than Java. + +### 9. ✅ FIXED — Inversive-distance gradient now uses limiting angles (consistent with Finding 1) +`inversive_distance_functional.hpp::inversive_distance_gradient` previously detected +a triangle-inequality-violating face two ways and **skipped it** in both: +- `if (l12sq <= 0 || l23sq <= 0 || l31sq <= 0) continue;` (non-real circle config) +- `if (!fa.valid) continue;` (degenerate-but-real triangle) + +The second skip was the *same* situation Finding 1 fixed for Euclidean/Spherical: +`euclidean_angles` returns the **limiting angles** (π opposite the over-long edge, +0/0 for the others) with `valid=false`, but this code threw them away. + +**Fix:** removed the `if (!fa.valid) continue;` skip so the limiting angles flow +into the gradient (the BPS-style convex C¹ extension), mirroring Finding 1; added a +comment explaining the rationale. The first skip (`l*sq <= 0`, a genuinely non-real +circle configuration with no limiting angle) is **kept**. All 237 tests pass. + +### 10. ℹ️ NOTE — Inversive-distance edge length has no overflow centring +`id_detail::edge_length_squared` computes `ℓ² = rᵢ² + rⱼ² + 2·I·rᵢ·rⱼ` with +`r = exp(u)` directly — no log-centring like `euclidean_angles`'s `μ`. The *angle* +computation is still safe (it re-centres internally on `log ℓ²`), but the raw `ℓ²` +can overflow for extreme `|u|`. Cosmetic robustness only; not hit in practice. + +### 11. ℹ️ NOTE — Verified correct (no action) +Audited against the literature and found faithful: +- **`trilaterate_2d / _sph / _hyp`** (`layout.hpp`) — circle-circle intersection + (ℝ²), spherical law of cosines via `α·pa+β·pb+γ·(pa×pb)`, and hyperbolic law of + cosines + Poincaré-disk Möbius placement. All three exact; left/CCW side correct. +- **`newton_solver.hpp`** — merit `f = ½‖G‖²`; Newton dir is always a descent dir + (`∇f·dx = −‖G‖²`), Armijo tests algebraically correct in both phases; spherical + solver factorises `−H` (PSD) and `d_sd = −H·G = −∇f` uses the un-negated H — all + consistent. SparseQR fallback handles the gauge null space (valid because + Gauss–Bonnet makes `G ⟂ 𝟙`). +- **`gauss_bonnet.hpp`** — `χ=V−E+F`, `g=(2−χ)/2`, `enforce` shift + `δ=(lhs−rhs)/V` makes the new sum equal `2π·χ` exactly. +- **Normalisations** — Euclidean PCA, spherical Rodrigues (mean→north), + hyperbolic weighted Fréchet-mean Möbius centring: all correct (spherical has a + benign antipodal-mean edge case that no-ops). + +--- + +## Missing test cases (gaps observed during the audit) + +1. **Degenerate / flipped-triangle gradient** — no test drives a face through the + triangle-inequality boundary to lock in Finding 1. Add: build a triangle with + one edge length > sum of the others, assert the gradient picks up the π corner + (not 0). Cover both Euclidean and Spherical. + +2. **`euclidean_hessian` edge-DOF guard** — no test asserts the new `throw` fires + when `assign_euclidean_all_dof_indices` is used. Add an `EXPECT_THROW`. + +3. ✅ **DONE (2026-05-29)** — Holonomy / period matrix end-to-end now covers + three tori of revolution (`HolonomyEndToEnd.{Torus4x4,TorusHex6x6,Torus8x8}`), + each asserting τ against the analytic modulus `i·√(R²−r²)/r` and `Re(τ) ≈ 0`. + (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. + +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). + +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. + +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. + +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 + edges and `genus == 2`, documenting that the basis is *not* canonical (so no τ + correctness is claimed there). + +10. **Inversive-distance degenerate face (Finding 9)** — no test drives an + inversive-distance triangle through the triangle-inequality boundary. Now that + Finding 9 mirrors Finding 1 (limiting angles used, not skipped), add the + analogue of missing-test item 1: assert the corner opposite the over-long edge + is picked up as π rather than silently skipped. diff --git a/doc/roadmap/research-track.md b/doc/roadmap/research-track.md index 2ec4892..5595fad 100644 --- a/doc/roadmap/research-track.md +++ b/doc/roadmap/research-track.md @@ -293,6 +293,20 @@ The phase numbers match `doc/roadmap/phases.md`. the cut-graph + holonomy infrastructure already in conformallab++). * **Effort:** large (10–14 days). * **Status:** roadmap item, no PR yet. +* **Known prerequisite bug (latent, 2026-05-29):** the holonomy-extraction + blocks in `spherical_layout` and `hyper_ideal_layout` (`layout.hpp`) + repeat the flawed single-development pattern that produced garbage τ for + the Euclidean path before the 2026-05-29 fix. They read the + translation / Möbius deck transformation from one full-surface + development plus a one-sided apex trilateration, instead of developing + across only the **dual** spanning tree and measuring the shared-edge + displacement between two independent developments (as the corrected + `detail::euclidean_holonomy` now does). These blocks are currently dead + code — every caller passes `holonomy == nullptr` — but Phase 9c/10b will + exercise the hyperbolic path. Fix = add `detail::spherical_holonomy` / + `detail::hyperbolic_holonomy` mirroring `detail::euclidean_holonomy`. + The hyperbolic mirror additionally needs `cpp_dec_float_50` (group-relation + product ∏gᵢ = Id overflows `double`; see CLAUDE.md high-precision note). ---