From b235666725b8a63e5fff1c645f623618e5ebfc61 Mon Sep 17 00:00:00 2001 From: Tarik Moussa Date: Mon, 18 May 2026 01:19:44 +0200 Subject: [PATCH] docs: Doxygen-API + Validierungsprotokoll + Porting-Tutorial MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Für einen Mathematiker der unabhängig validieren und eigene Forschung einbringen möchte. Doxygen-Kommentare (code/include/): newton_solver.hpp — newton_euclidean(), newton_spherical(), newton_hyper_ideal() je mit \param, \return, \note, \see inkl. mathematischer Begründung (Konvexität, Vorzeichenkonvention, SparseQR-Fallback-Erklärung) layout.hpp — euclidean_layout(), spherical_layout(), hyper_ideal_layout() mit vollständiger Parameter-Doku, halfedge_uv-Semantik, Poincaré-Disk-Note Neues Dokument: doc/math/validation-protocol.md 7 reproduzierbare Checks mit konkreten Befehlen und erwartetem Output: 0. 170 Tests, 1 Skip 1. Gauss–Bonnet exakt (1e-10) 2. FD-Gradientencheck < 1e-6 für alle 3 Geometrien 3. Newton-Konvergenz < 50 Iterationen 4. τ ∈ SL(2,ℤ)-Fundamentaldomäne (3 Invarianten) 5. Möbius-Arithmetik (Inverse, Compose, from_three) 6. End-to-End-Pipeline 7. Manueller τ-Check für torus_4x4.off (Codebeispiel) Neues Tutorial: doc/tutorials/add-inversive-distance.md Vollständiger Step-by-Step-Port von Phase 9a (Luo 2004): Header anlegen, Energie/Gradient implementieren, FD-Check, Newton-Wrapper, CMakeLists, Java-Referenzvergleich, Checkliste. doc/getting-started.md: Abschnitt "Known issues": macOS-Finder-Duplikate (rm-Befehl), Warnung "First build 30–90s" (Tarball-Extraktion) README.md: Zwei neue Links in der Dokumentationstabelle (validation-protocol, tutorial) Co-Authored-By: Claude Sonnet 4.6 --- README.md | 2 + code/include/layout.hpp | 60 +++++++ code/include/newton_solver.hpp | 88 ++++++++-- doc/getting-started.md | 24 +++ doc/math/validation-protocol.md | 213 ++++++++++++++++++++++++ doc/tutorials/add-inversive-distance.md | 202 ++++++++++++++++++++++ 6 files changed, 571 insertions(+), 18 deletions(-) create mode 100644 doc/math/validation-protocol.md create mode 100644 doc/tutorials/add-inversive-distance.md diff --git a/README.md b/README.md index 0e3be55..d2d3292 100644 --- a/README.md +++ b/README.md @@ -85,6 +85,8 @@ Layout2D layout = euclidean_layout(mesh, res.x, maps); | **Project structure** — directory tree + build targets | [doc/architecture/project-structure.md](doc/architecture/project-structure.md) | | **Discrete conformal theory** — mathematical background for collaborators | [doc/math/discrete-conformal-theory.md](doc/math/discrete-conformal-theory.md) | | **Validation** — known analytic results + how to verify them | [doc/math/validation.md](doc/math/validation.md) | +| **Validation protocol** — concrete commands with expected outputs | [doc/math/validation-protocol.md](doc/math/validation-protocol.md) | +| **Tutorial: add a new functional** — step-by-step Inversive-Distance port | [doc/tutorials/add-inversive-distance.md](doc/tutorials/add-inversive-distance.md) | | **Geometry modes** — Euclidean / Spherical / HyperIdeal comparison | [doc/math/geometry-modes.md](doc/math/geometry-modes.md) | | **References** — all papers by module | [doc/math/references.md](doc/math/references.md) | | **Roadmap** — Phases 1–10 | [doc/roadmap/phases.md](doc/roadmap/phases.md) | diff --git a/code/include/layout.hpp b/code/include/layout.hpp index 8b5f1ed..70e7fe1 100644 --- a/code/include/layout.hpp +++ b/code/include/layout.hpp @@ -452,6 +452,29 @@ inline void set_root_huv_2d( } // namespace detail // ── Euclidean layout ────────────────────────────────────────────────────────── + +/// Embed the mesh in ℝ² using the Euclidean metric encoded in x. +/// +/// Runs priority-BFS trilateration: places vertices in order of BFS depth +/// from the root face (largest 3D area), so errors accumulate last. +/// For closed genus-g surfaces a CutGraph must be supplied — otherwise +/// the layout will have a seam discontinuity (Layout2D::has_seam = true). +/// +/// \param mesh Input surface mesh. Must have lambda0 and v_idx set in maps. +/// \param x DOF vector returned by newton_euclidean(). +/// \param maps EuclideanMaps (lambda0, v_idx, e_idx). +/// \param cut Optional cut graph (compute_cut_graph()). Pass nullptr for +/// open meshes or if seams are acceptable. +/// \param holonomy If non-null and cut != nullptr, receives the lattice +/// translations ω_i ∈ ℂ per cut edge. +/// Pass to compute_period_matrix() for the conformal modulus τ. +/// \param normalise If true, calls normalise_euclidean() on the result: +/// centroid → origin, major axis → x-axis (PCA). +/// \return Layout2D with .uv[v] (per-vertex UV) and +/// .halfedge_uv[h] (per-halfedge UV for texture atlasing). +/// +/// \note halfedge_uv differs from uv at seam edges: the two sides of a cut +/// carry different UV coordinates for proper GPU texture atlasing. inline Layout2D euclidean_layout( ConformalMesh& mesh, const std::vector& x, @@ -571,6 +594,20 @@ inline Layout2D euclidean_layout( } // ── Spherical layout ────────────────────────────────────────────────────────── + +/// Embed the mesh on the unit sphere S² using the spherical metric encoded in x. +/// +/// Runs priority-BFS trilateration using the spherical law of cosines. +/// Typical use: genus-0 (sphere-like) surfaces after newton_spherical(). +/// +/// \param mesh Input genus-0 surface mesh. +/// \param x DOF vector returned by newton_spherical(). +/// \param maps SphericalMaps. +/// \param cut Optional cut graph (rarely needed for genus-0). +/// \param holonomy If non-null, receives rotational holonomies (spherical). +/// \param normalise If true, calls normalise_spherical(): rotates the centroid +/// to the north pole (Rodrigues rotation formula). +/// \return Layout3D with .xyz[v] ∈ S² ⊂ ℝ³ for each vertex. inline Layout3D spherical_layout( ConformalMesh& mesh, const std::vector& x, @@ -670,6 +707,29 @@ inline Layout3D spherical_layout( } // ── HyperIdeal layout (Poincaré disk) — exact trilateration ────────────────── + +/// Embed the mesh in the Poincaré disk (H²) using the hyperbolic metric encoded in x. +/// +/// Runs priority-BFS trilateration using exact Möbius-isometric placement: +/// each new vertex is located by solving the hyperbolic law of cosines and +/// applying a Möbius map to position it in the disk. +/// For closed genus-g surfaces (g ≥ 1) a CutGraph is required. +/// +/// \param mesh Input genus-g surface mesh (g ≥ 1). +/// \param x DOF vector returned by newton_hyper_ideal() +/// (vertex b_v and edge a_e variables). +/// \param maps HyperIdealMaps (lambda0, v_idx, e_idx). +/// \param cut CutGraph from compute_cut_graph(). Required for closed surfaces. +/// \param holonomy If non-null, receives the Möbius maps T_i ∈ SU(1,1) per cut edge. +/// Pass to compute_period_matrix() for holonomy analysis. +/// \param normalise If true, calls normalise_hyperbolic(): iterative face-area-weighted +/// Möbius centring (Fréchet mean, 30 iterations) → disk origin. +/// \return Layout2D with .uv[v] ∈ Poincaré disk (|uv| < 1) for each vertex, +/// and .halfedge_uv[h] for seam-aware texture atlasing. +/// +/// \note All vertex positions satisfy |uv[v]| < 1 (inside the Poincaré disk) +/// if the metric is hyperbolic. Points on or outside the boundary indicate +/// a non-hyperbolic metric (Gauss–Bonnet violation). inline Layout2D hyper_ideal_layout( ConformalMesh& mesh, const std::vector& x, diff --git a/code/include/newton_solver.hpp b/code/include/newton_solver.hpp index dcb222c..ae26a6e 100644 --- a/code/include/newton_solver.hpp +++ b/code/include/newton_solver.hpp @@ -139,9 +139,29 @@ inline std::vector line_search( } // namespace detail // ── Euclidean Newton solver ──────────────────────────────────────────────────── -// -// Minimises the Euclidean discrete conformal energy by solving G(x) = 0. -// The Hessian H is PSD; Eigen::SimplicialLDLT is used directly. + +/// Solve the Euclidean discrete conformal problem: find u ∈ ℝ^V such that +/// Σ_{faces adj v} α_v(u) = Θ_v for all vertices v. +/// +/// Starting from x0, Newton's method minimises E(u) (the Euclidean DCE energy, +/// which is convex) by iterating u ← u − H⁻¹·G with backtracking line search. +/// The Hessian H is the cotangent Laplacian — PSD with one zero eigenvalue on +/// closed surfaces (gauge mode). A SparseQR fallback handles this automatically. +/// +/// \param mesh Input triangulated surface (edges must carry lambda0 + theta_v). +/// \param x0 Initial DOF vector (length = number of free vertices). +/// Pass all-zeros for a flat start (typical). +/// \param m EuclideanMaps: lambda0[e], theta_v[v], v_idx[v] must be set. +/// Call setup_euclidean_maps() + compute_euclidean_lambda0_from_mesh() +/// + enforce_gauss_bonnet() before passing here. +/// \param tol Convergence threshold on max |G_i|. Default: 1e-8. +/// \param max_iter Maximum Newton iterations. Default: 200. +/// \return NewtonResult{x*, iterations, grad_inf_norm, converged}. +/// +/// \note On closed meshes without a pinned vertex, SimplicialLDLT detects the +/// gauge singularity and falls back to SparseQR automatically. +/// +/// \see doc/math/discrete-conformal-theory.md §3 for the mathematical background. inline NewtonResult newton_euclidean( ConformalMesh& mesh, std::vector x0, @@ -199,10 +219,28 @@ inline NewtonResult newton_euclidean( } // ── Spherical Newton solver ─────────────────────────────────────────────────── -// -// Solves G(x) = 0 for the spherical discrete conformal functional. -// The Hessian H is NSD at the solution; −H is PSD, so we factorise −H and -// solve (−H)·Δx = G ⟺ H·Δx = −G. + +/// Solve the spherical discrete conformal problem: find u ∈ ℝ^V such that +/// Σ_{faces adj v} α_v(u) = Θ_v for all vertices v (genus-0 / sphere-like surfaces). +/// +/// The spherical DCE energy is *concave*, so the Hessian H is NSD at the solution. +/// The solver factorises −H (which is PSD) and solves (−H)·Δx = G. +/// A gauge vertex must be pinned (set v_idx = -1) to remove the rotational mode. +/// +/// \param mesh Input triangulated surface, genus 0. +/// \param x0 Initial DOF vector (length = free vertices, excluding gauge_vertex). +/// All-zeros is a good start. +/// \param m SphericalMaps: lambda0[e], theta_v[v], v_idx[v], gauge_vertex set. +/// Call setup_spherical_maps() + compute_spherical_lambda0_from_mesh() +/// + enforce_gauss_bonnet() (checks Σ(2π-Θ) > 0) before passing here. +/// \param tol Convergence threshold on max |G_i|. Default: 1e-8. +/// \param max_iter Maximum Newton iterations. Default: 200. +/// \return NewtonResult{x*, iterations, grad_inf_norm, converged}. +/// +/// \note Unlike the Euclidean solver, the spherical solver does NOT need a SparseQR +/// fallback — the gauge vertex pins the null mode directly. +/// +/// \see doc/math/geometry-modes.md §Spherical for sign-convention details. inline NewtonResult newton_spherical( ConformalMesh& mesh, std::vector x0, @@ -258,17 +296,31 @@ inline NewtonResult newton_spherical( } // ── HyperIdeal Newton solver ────────────────────────────────────────────────── -// -// Solves G(x) = 0 for the hyper-ideal discrete conformal functional. -// -// Gradient sign convention (opposite to Euclidean/Spherical): -// G_v = Σ β_v − Θ_v, G_e = Σ α_e − θ_e (actual − target) -// -// The hyper-ideal energy is strictly convex (Springborn 2020), so H is PSD -// and SimplicialLDLT (with SparseQR fallback) applies directly. -// -// The Hessian is computed by symmetric finite differences of G (see -// hyper_ideal_hessian.hpp); replace with an analytical Hessian in Phase 5. + +/// Solve the hyper-ideal discrete conformal problem: find (b, a) ∈ ℝ^{V+E} such that +/// Σ β_v(b,a) = Θ_v and Σ α_e(b,a) = θ_e for all vertices v and edges e. +/// +/// Used for genus-g surfaces (g ≥ 1) under hyperbolic cone metrics. +/// The energy is *strictly convex* (Springborn 2020, Theorem 1.3), so Newton +/// converges globally from any starting point. +/// +/// DOF layout: first V_free entries are vertex variables b_v (hyper-ideal radii), +/// followed by E entries for edge variables a_e (intersection angles). +/// Use assign_all_dof_indices(mesh, maps) to set v_idx and e_idx automatically — +/// no vertex needs to be pinned. +/// +/// \param mesh Input triangulated surface, genus g ≥ 1. +/// \param x0 Initial DOF vector (length = V + E). All-zeros typical. +/// \param m HyperIdealMaps: lambda0[e], theta_v[v], v_idx[v], e_idx[e] set. +/// Call setup_hyper_ideal_maps() + compute_hyper_ideal_lambda0_from_mesh(). +/// \param tol Convergence threshold on max |G_i|. Default: 1e-8. +/// \param max_iter Maximum Newton iterations. Default: 200. +/// \param hess_eps Finite-difference step for Hessian approximation. Default: 1e-5. +/// (Phase 9b will replace this with an analytic Hessian.) +/// \return NewtonResult{x*, iterations, grad_inf_norm, converged}. +/// +/// \see Springborn (2020), Theorem 1.3 for the strict convexity proof. +/// \see doc/math/geometry-modes.md §Hyper-ideal for DOF layout details. inline NewtonResult newton_hyper_ideal( ConformalMesh& mesh, std::vector x0, diff --git a/doc/getting-started.md b/doc/getting-started.md index a86372e..a5a70f4 100644 --- a/doc/getting-started.md +++ b/doc/getting-started.md @@ -127,6 +127,30 @@ After a full build (`-DWITH_CGAL=ON`): --- +## Known issues + +### macOS Finder duplicates + +macOS Finder sometimes creates duplicate files named `foo 2.hpp` when copying +the repository. These cause confusing compile errors ("redefinition of …"). + +**Fix (run once from the repo root):** +```bash +find code/include -name "* 2.*" -delete +find code/include -name "*\ 2.*" -delete +``` + +Files without a ` 2` suffix are always canonical — the duplicates are safe to delete. + +### First build is slow + +The first CMake configure extracts four tarballs (Eigen 3.4, CGAL 6.1.1, +libigl 2.6, GLFW 3.4) and downloads GTest via `FetchContent`. +Allow **30–90 seconds** for the first configure. Subsequent builds are fast +(< 10 s incremental). + +--- + ## Rebuilding the CI Docker image The CI runner is a self-hosted Raspberry Pi (ARM64). After changes to diff --git a/doc/math/validation-protocol.md b/doc/math/validation-protocol.md new file mode 100644 index 0000000..3fb60a0 --- /dev/null +++ b/doc/math/validation-protocol.md @@ -0,0 +1,213 @@ +# Validation Protocol + +Concrete, reproducible steps to verify the mathematical correctness of +conformallab++. Every check below has a deterministic expected outcome. + +--- + +## Prerequisites + +```bash +cmake -S code -B build -DWITH_CGAL=ON -DCMAKE_BUILD_TYPE=Release +cmake --build build --target conformallab_cgal_tests -j$(nproc) +``` + +--- + +## Check 0 — All tests pass + +```bash +ctest --test-dir build -R "^cgal\." --output-on-failure +``` + +**Expected output (last lines):** +``` +100% tests passed, 0 tests failed out of 170 + +The following tests did not run: + 206 - cgal.HomologyGenerators.Genus2_FourGeneratorPaths_BLOCKED (Skipped) +``` + +If any test fails, stop — the implementation is broken. + +--- + +## Check 1 — Gauss–Bonnet (topological identity, error < 1e-10) + +```bash +./build/conformallab_cgal_tests --gtest_filter="GaussBonnet.*" -v +``` + +**Expected:** all 8 tests `[ PASSED ]` + +What is verified: for each test mesh (tetrahedron χ=2, torus χ=0, open mesh χ=1): + +``` +Σᵥ (2π − Θᵥ) = 2π · χ(M) ± 1e-10 +``` + +This is a **pure topology check** — it fails only if vertex/face counts or +property-map assignments are wrong. It does not depend on the Newton solver. + +--- + +## Check 2 — Euclidean gradient consistency (FD vs. analytic, error < 1e-6) + +```bash +./build/conformallab_cgal_tests --gtest_filter="EuclideanFunctional.GradientCheck*" -v +``` + +**Expected:** all `GradientCheck_*` tests `[ PASSED ]` + +What is verified: for ε = 1e-5, + +``` +|G(u)ᵢ − (E(u + εeᵢ) − E(u − εeᵢ)) / (2ε)| < 1e-6 +``` + +This check **proves that the energy and its gradient are mathematically consistent**. +A failing FD check means Newton will converge to the wrong point — it is the most +important correctness check for any new functional. + +Also run for Spherical and HyperIdeal: +```bash +./build/conformallab_cgal_tests --gtest_filter="SphericalFunctional.GradientCheck*" -v +./build/conformallab_cgal_tests --gtest_filter="HyperIdealFunctional.GradientCheck*" -v +``` + +--- + +## Check 3 — Newton convergence on canonical test meshes + +```bash +./build/conformallab_cgal_tests --gtest_filter="NewtonSolver.*" -v +``` + +**Expected:** all 11 tests `[ PASSED ]` + +Each test verifies: +- `res.converged == true` +- `res.grad_inf_norm < 1e-8` +- `res.iterations < 50` (typically 5–20 for the small test meshes) + +--- + +## Check 4 — Period matrix: SL(2,ℤ)-reduction invariants + +```bash +./build/conformallab_cgal_tests --gtest_filter="PeriodMatrix.*" -v +``` + +**Expected:** all 7 tests `[ PASSED ]` + +The three mathematical invariants checked for **any** genus-1 output τ: + +| Property | Condition | Why | +|---|---|---| +| Upper half-plane | `Im(τ) > 0` | τ encodes a positive-area lattice | +| Outside unit disk | `|τ| ≥ 1 − 1e-10` | SL(2,ℤ) reduction step | +| Vertical strip | `|Re(τ)| ≤ 0.5 + 1e-10` | SL(2,ℤ) reduction step | + +These hold for **any** well-formed genus-1 mesh — they are topology, not geometry. + +--- + +## Check 5 — Möbius arithmetic (complex analysis correctness) + +```bash +./build/conformallab_cgal_tests --gtest_filter="MobiusMap.*" -v +``` + +**Expected:** all 8 tests `[ PASSED ]` + +What is verified: +- `T ∘ T⁻¹ = Id` (inverse is correct) +- `(T₁ ∘ T₂)(z) = T₁(T₂(z))` (composition is associative) +- `from_three(z₁, z₂, z₃)` maps z₁→0, z₂→1, z₃→∞ (unique Möbius transformation) + +A bug here would corrupt **all** hyperbolic holonomy computation. + +--- + +## Check 6 — End-to-end pipeline (build + solve + layout) + +```bash +./build/conformallab_cgal_tests --gtest_filter="Pipeline.*" -v +``` + +**Expected:** all 5 tests `[ PASSED ]` + +What is verified: starting from a mesh file, the full pipeline +(setup → Gauss-Bonnet → Newton → layout → serialise → reload) +produces a consistent result. + +--- + +## Check 7 — Manual torus τ verification + +This check requires adding a small program (or modifying an existing test). +It validates that `torus_4x4.off` produces τ in the fundamental domain: + +```cpp +#include "conformal_mesh.hpp" +#include "euclidean_functional.hpp" +#include "newton_solver.hpp" +#include "cut_graph.hpp" +#include "layout.hpp" +#include "period_matrix.hpp" +#include "mesh_io.hpp" +#include + +int main() { + conformallab::ConformalMesh mesh; + conformallab::load_mesh(mesh, "code/data/off/torus_4x4.off"); + + auto maps = conformallab::setup_euclidean_maps(mesh); + conformallab::compute_euclidean_lambda0_from_mesh(mesh, maps); + conformallab::enforce_gauss_bonnet(mesh, maps); + + auto res = conformallab::newton_euclidean(mesh, std::vector(maps.n_dof, 0.0), maps); + std::cout << "Converged: " << res.converged + << " iterations: " << res.iterations + << " |G|∞: " << res.grad_inf_norm << "\n"; + + auto cg = conformallab::compute_cut_graph(mesh); + conformallab::HolonomyData hol; + conformallab::euclidean_layout(mesh, res.x, maps, &cg, &hol, true); + auto pd = conformallab::compute_period_matrix(hol); + + std::cout << "τ = " << pd.tau_reduced.real() + << " + " << pd.tau_reduced.imag() << "i\n"; + std::cout << "|τ| = " << std::abs(pd.tau_reduced) << "\n"; + std::cout << "|Re(τ)| = " << std::abs(pd.tau_reduced.real()) << "\n"; +} +``` + +**Expected output (torus_4x4.off, R=2, r=1 torus of revolution):** +``` +Converged: 1 iterations: <30 |G|∞: <1e-8 +τ = [small] + [positive]i (Re close to 0 by 4-fold symmetry) +|τ| ≥ 1.0 (fundamental domain) +|Re(τ)| ≤ 0.5 (fundamental domain) +``` + +The exact value of Im(τ) depends on the 3D embedding (R=2, r=1 gives unequal +inner/outer edge lengths). Use `torus_8x8.off` for a finer approximation. + +--- + +## Summary checklist + +``` +[ ] Check 0: 170 tests pass, 1 skip +[ ] Check 1: Gauss–Bonnet exact (1e-10) +[ ] Check 2: FD gradient < 1e-6 for all 3 geometries +[ ] Check 3: Newton convergence < 50 iterations +[ ] Check 4: τ in SL(2,ℤ) fundamental domain +[ ] Check 5: Möbius arithmetic (inverse, compose, from_three) +[ ] Check 6: End-to-end pipeline +[ ] Check 7: Torus τ in upper half-plane with correct symmetry +``` + +All checks are deterministic and do not depend on random initialization or +floating-point non-determinism beyond standard IEEE-754. diff --git a/doc/tutorials/add-inversive-distance.md b/doc/tutorials/add-inversive-distance.md new file mode 100644 index 0000000..96bc2b1 --- /dev/null +++ b/doc/tutorials/add-inversive-distance.md @@ -0,0 +1,202 @@ +# Tutorial: Porting the Inversive-Distance Functional (Phase 9a) + +This is a complete, step-by-step example of how to add a new functional +to conformallab++. It ports `InversiveDistanceFunctional.java` from the +Java reference implementation (Luo 2004). + +**Prerequisite:** Read [doc/api/extending.md](../api/extending.md) first +for the general pattern. This tutorial fills in every detail for one +specific case. + +--- + +## Mathematical background + +Inversive distance (Luo 2004) uses a different edge-length update formula: + +``` +Given inversive distances I_{ij} ∈ ℝ for each edge {i,j}: + + cosh(l̃_{ij}) = I_{ij} · cosh((u_i + u_j) / 2) + + (cosh²((u_i - u_j) / 2) - 1) · ... +``` + +For the Euclidean version the angle formula is the same law of cosines, +but the log-lengths Λ̃ are computed differently from the u-vector. + +**Reference:** Luo, F. (2004). *Combinatorial Yamabe Flow on Surfaces*. +Communications in Contemporary Mathematics, 6(5), 765–780. + +**Java source:** `de.varylab.discreteconformal.functional.InversiveDistanceFunctional` + +--- + +## Step 1 — Create the header + +```bash +cp code/include/euclidean_functional.hpp \ + code/include/inversive_distance_functional.hpp +``` + +Edit the new file: + +```cpp +#pragma once +// inversive_distance_functional.hpp +// +// Phase 9a — Inversive-distance discrete conformal functional (Luo 2004). +// Ported from de.varylab.discreteconformal.functional.InversiveDistanceFunctional. +// +// Usage: identical to euclidean_functional.hpp — replace lambda0 with +// inversive_distance0 (the initial inversive distances per edge). + +#include "conformal_mesh.hpp" +#include +#include + +namespace conformallab { + +struct InversiveDistanceMaps { + CGAL::Surface_mesh::Point_3> + ::Property_map inv_dist0; ///< I_{ij} per edge + CGAL::Surface_mesh::Point_3> + ::Property_map theta_v; ///< target corner angles Θ_v + CGAL::Surface_mesh::Point_3> + ::Property_map v_idx; ///< DOF index (-1 = pinned) +}; + +// ... (follow euclidean_functional.hpp exactly, replacing lambda0 with inv_dist0 +// and the angle formula with the Luo (2004) version) +``` + +--- + +## Step 2 — Implement the energy, gradient, and angle formula + +The Euclidean angle formula uses the law of cosines on edge lengths +derived from log-lengths. For inversive distance, the edge lengths +are derived from inversive distances I_{ij} and the u-vector: + +```cpp +// Euclidean (reference): +// lambda_ij = lambda0_ij + u_i + u_j +// l_ij = exp(lambda_ij / 2) + +// Inversive distance (Luo 2004): +// cosh(l̃_ij) = I_ij * cosh((u_i + u_j) / 2) [simplified form] +// l̃_ij = acosh(I_ij * cosh((u_i + u_j) / 2)) +``` + +Then the **corner angle** at vertex k opposite edge {i,j} in triangle {i,j,k} +is computed by the standard law of cosines from l̃_{ij}, l̃_{jk}, l̃_{ki}. + +The **gradient** is the same as Euclidean: + +```cpp +G_v = Σ_{faces containing v} alpha_v(face, u) − Theta_v +``` + +--- + +## Step 3 — Write the gradient-check test + +Create `code/tests/cgal/test_inversive_distance.cpp`: + +```cpp +#include "inversive_distance_functional.hpp" +#include "mesh_factory.hpp" +#include + +// Finite-difference gradient check: copy the pattern from +// test_euclidean_functional.cpp :: EuclideanFunctional.GradientCheck_Triangle +TEST(InversiveDistance, GradientCheck_Triangle) +{ + // Build a single equilateral triangle (open mesh, no boundary issues). + ConformalMesh mesh = MeshFactory::make_open_mesh(MeshFactory::Kind::triangle); + InversiveDistanceMaps maps = setup_inversive_distance_maps(mesh); + compute_inversive_distance0_from_mesh(mesh, maps); // I_ij from 3D edge lengths + + const int n = /* count free DOFs */; + std::vector x0(n, 0.0); + constexpr double eps = 1e-5; + + auto G = inversive_distance_gradient(mesh, x0, maps); + + for (int i = 0; i < n; ++i) { + auto xp = x0; xp[i] += eps; + auto xm = x0; xm[i] -= eps; + double Ep = inversive_distance_energy(mesh, xp, maps); + double Em = inversive_distance_energy(mesh, xm, maps); + double fd = (Ep - Em) / (2.0 * eps); + EXPECT_NEAR(G[i], fd, 1e-6) << "gradient mismatch at DOF " << i; + } +} +``` + +Run with: + +```bash +./build/conformallab_cgal_tests --gtest_filter="InversiveDistance.*" +``` + +--- + +## Step 4 — Register in CMakeLists + +Add to `code/tests/cgal/CMakeLists.txt` inside the `add_executable` block: + +```cmake + # ── Phase 9a: Inversive-distance functional ──────────────────────────────── + test_inversive_distance.cpp +``` + +--- + +## Step 5 — Add a Newton wrapper + +In `code/include/newton_solver.hpp` add: + +```cpp +/// Solve the inversive-distance conformal problem. +/// \see newton_euclidean() — identical structure. +inline NewtonResult newton_inversive_distance( + ConformalMesh& mesh, + std::vector x0, + const InversiveDistanceMaps& m, + double tol = 1e-8, + int max_iter = 200) +{ + // Copy newton_euclidean() exactly, replacing euclidean_* with + // inversive_distance_*. The Hessian is PSD (Luo 2004, Thm. 1), + // so SimplicialLDLT with SparseQR fallback applies directly. +} +``` + +--- + +## Step 6 — Verify against Java output + +The Java library outputs text results for a triangulated torus. To compare: + +1. Run Java ConformalLab on the same OFF mesh with Inversive-Distance mode. +2. Record the final gradient norm and angle sums. +3. Run the C++ equivalent and check: + +```cpp +EXPECT_LT(res.grad_inf_norm, 1e-8); +EXPECT_LT(res.iterations, 50); // inversive-distance converges faster than hyper-ideal +``` + +**Java reference:** `InversiveDistanceFunctionalTest.java` in `de.varylab.discreteconformal.test` + +--- + +## Checklist + +- [ ] `code/include/inversive_distance_functional.hpp` compiles +- [ ] `GradientCheck_Triangle` passes +- [ ] `GradientCheck_OpenMesh` passes (copy from Euclidean tests) +- [ ] Newton converges on `torus_4x4.off` +- [ ] Result matches Java output (gradient norm < 1e-8 on same mesh) +- [ ] Registered in `CMakeLists.txt` +- [ ] `doc/roadmap/java-parity.md` updated: Phase 9a → ✅