# Tutorial: Implementing the Inversive-Distance Functional (Phase 9a.2) This tutorial walks through adding a **new** discrete-conformal functional to conformallab++. The running example is the **vertex-based inversive- distance functional** of Luo (2004), used as Phase 9a.2 of the roadmap. > ## ⚠️ This is research, not a port > > An earlier draft of this document claimed this functional was a port of > `de.varylab.discreteconformal.functional.InversiveDistanceFunctional`. > **That Java class does not exist.** Verified empirically: > > ```bash > $ find /Users/tarikmoussa/Desktop/conformallab -iname "*nversive*" > (zero results) > $ grep -r "InversiveDistance" /Users/tarikmoussa/Desktop/conformallab/src > (zero matches) > ``` > > The closest Java cousin is `CPEuclideanFunctional.java`, which implements > the **face-based** circle-packing variant (Phase 9a.1). The vertex-based > inversive-distance functional (this tutorial) is built **from the > literature**, not from a Java reference, and the correctness validation > is cross-checked against three sources: > > 1. **Luo, F.** (2004). *Combinatorial Yamabe Flow on Surfaces.* Comm. Contemp. Math. 6(5), 765–780. > 2. **Bowers, P. L. & Stephenson, K.** (2004). *Uniformizing dessins and Belyĭ maps via circle packing.* Mem. AMS 170(805). > 3. **Glickenstein, D.** (2011). *Discrete conformal variations and scalar curvature on piecewise flat manifolds.* J. Differential Geometry 87(2), 201–238. > > The tutorial below has been re-written to match this reality. **Prerequisite:** Read [doc/api/extending.md](../api/extending.md) first for the general functional-porting pattern. This tutorial fills in the mathematical and code details for one specific case. --- ## Mathematical background Inversive-distance circle packing parametrises **each vertex** by a circle of radius `r_i = exp(u_i)`. Two adjacent circles have an *inversive distance* `I_ij` that is a fixed constant of the edge, derived once from the initial geometry via the Bowers–Stephenson identity: ``` I_ij = ( ℓ_ij² − r_i² − r_j² ) / ( 2 r_i r_j ) (Bowers-Stephenson 2004) ``` Geometric interpretation of `I_ij`: | Range | Configuration | |---|---| | `I_ij = +1` | tangent circles (Koebe-style) | | `I_ij ∈ (0, 1)` | overlapping with intersection angle `φ`, `I = cos φ` | | `I_ij = 0` | orthogonal circles | | `I_ij ∈ (−1, 0)` | disjoint circles, inversive distance > 1 | | `I_ij ≤ −1` | impossible packing | The edge length under a state `u` is then determined by Luo's formula: ``` ℓ_ij(u)² = exp(2 u_i) + exp(2 u_j) + 2 I_ij exp(u_i + u_j) = r_i² + r_j² + 2 I_ij r_i r_j (Luo 2004 §3) ``` The angle formula is the same numerically-stable half-tangent law of cosines used by `euclidean_functional.hpp`; only the way `ℓ_ij` is computed from `u` is different. The gradient is the standard Yamabe-flow gradient: ``` ∂E/∂u_v = Θ_v − Σ_{T ∋ v} α_v(T) (Luo 2004 Lemma 3.1) ``` The energy is a path integral of the gradient (Luo's 1-form is closed on the domain where every triangle is valid); we use the same 10-point Gauss-Legendre quadrature as `euclidean_functional.hpp`. The Hessian is finite-difference for the MVP; an analytic form (Glickenstein 2011 §5.2) is documented in the research-track roadmap. --- ## Step 1 — Create the header ```bash cp code/include/euclidean_functional.hpp \ code/include/inversive_distance_functional.hpp ``` Modify the maps struct: replace `lambda0` (Euclidean log-length) with the inversive-distance constant `I_e` and the initial radius `r0` (used for the Bowers-Stephenson init). ```cpp struct InversiveDistanceMaps { ConformalMesh::Property_map v_idx; // DOF index (−1 = pinned) ConformalMesh::Property_map theta_v; // target cone angle ConformalMesh::Property_map r0; // initial radius r_i^(0) ConformalMesh::Property_map I_e; // inversive distance per edge }; ``` Then implement the four entry points that any functional needs in conformallab++: - `setup_inversive_distance_maps(mesh)` — create maps with defaults. - `compute_inversive_distance_init_from_mesh(mesh, m)` — choose `r_i^(0)` from the input geometry, then compute `I_ij` via Bowers-Stephenson. - `inversive_distance_gradient(mesh, x, m)` — Luo's `Θ − Σ α`. - `inversive_distance_energy(mesh, x, m)` — 10-point Gauss-Legendre path integral. For the full implementation, see `code/include/inversive_distance_functional.hpp` (part of PR #8). --- ## Step 2 — Edge-length kernel The single new pure-math primitive is the Luo edge-length formula. Wrap it in a small detail helper so the gradient function reads cleanly: ```cpp namespace id_detail { // ℓ² = exp(2u_i) + exp(2u_j) + 2 I exp(u_i + u_j) // Returns -1 on degenerate input (no valid packing). inline double edge_length_squared(double u_i, double u_j, double I_ij) { double ri = std::exp(u_i); double rj = std::exp(u_j); double l2 = ri*ri + rj*rj + 2.0 * I_ij * ri * rj; return l2 > 0.0 ? l2 : -1.0; } } // namespace id_detail ``` This is the only place where the inversive-distance model differs from the Euclidean one. All downstream code (angle computation, gradient accumulation, energy integration) is structurally identical. --- ## Step 3 — Reuse `euclidean_angles()` The half-tangent law of cosines is independent of how lengths were obtained. Feed `log(ℓ²)` to the existing helper to compute the three corner angles per face: ```cpp auto fa = euclidean_angles(std::log(l12sq), std::log(l23sq), std::log(l31sq)); ``` This is the **non-trivial reuse** that justifies the structural similarity to the Euclidean functional — we get the law-of-cosines numerics for free, and only the edge-length input changes. --- ## Step 4 — Validation tests The acceptance criteria for this functional are stricter than for a Java port because there is no reference implementation to compare against. We need **three independent validations**: ### 4.1 Limit-case edge lengths Each of Luo's special cases (`I = 1` tangent, `I = 0` orthogonal, `I = −1` inside-tangent) gives a closed-form `ℓ` that must be reproduced to machine precision: ```cpp TEST(InversiveDistanceFunctional, EdgeLengthFormula_TangentialLimit) { // r_i=1, r_j=2, I=1: ℓ² = 1 + 4 + 2·1·1·2 = 9 ⇒ ℓ = 3 = r_i + r_j double l2 = id_detail::edge_length_squared(0.0, std::log(2.0), 1.0); EXPECT_NEAR(std::sqrt(l2), 3.0, 1e-12); } ``` Three such tests cover the diagnostic special cases. ### 4.2 Bowers-Stephenson round-trip The initialisation `compute_inversive_distance_init_from_mesh` must be self-consistent: starting from `(ℓ_3d, r_i, r_j)` and computing `I_ij`, the round-trip back through Luo's formula must give the original `ℓ`. ```cpp TEST(InversiveDistanceFunctional, BowersStephensonRoundTrip) { auto mesh = make_triangle(); auto m = setup_inversive_distance_maps(mesh); compute_inversive_distance_init_from_mesh(mesh, m); for (auto e : mesh.edges()) { double l_3d = /* mesh 3-D edge length */; double ri = m.r0[mesh.source(mesh.halfedge(e))]; double rj = m.r0[mesh.target(mesh.halfedge(e))]; double l_rec = std::sqrt(ri*ri + rj*rj + 2.0 * m.I_e[e] * ri * rj); EXPECT_NEAR(l_rec, l_3d, 1e-12); } } ``` ### 4.3 FD-vs-analytic gradient check Standard pattern from every functional in conformallab++ — see `test_euclidean_functional.cpp`. Compare the analytic gradient to a symmetric finite difference of the energy. ### 4.4 Cross-validation against `euclidean_functional.hpp` At `u = 0`, both functionals reconstruct the input 3-D edge length exactly (Euclidean via `compute_lambda0`, inversive distance via Bowers-Stephenson). Therefore the actual angle sums are identical, and the two gradients (with default `Θ_v = 2π`) must match component-wise: ```cpp TEST(InversiveDistanceFunctional, AngleDefectAtU0_AgreesWithEuclideanAtU0) { auto G_id = inversive_distance_gradient(mesh, /*x=0*/, m_id); auto G_eu = euclidean_gradient (mesh, /*x=0*/, m_eu); for (size_t i = 0; i < G_id.size(); ++i) EXPECT_NEAR(G_id[i], G_eu[i], 1e-10); } ``` This is the empirical statement of Glickenstein 2011 §5: different parametrisations of the same initial discrete metric produce the same Newton-time-zero gradient. --- ## Step 5 — Register the tests In `code/tests/cgal/CMakeLists.txt`: ```cmake # ── Phase 9a.2: InversiveDistance (Luo 2004 + Glickenstein 2011) ───────── # Vertex-based inversive-distance circle-packing functional. No Java # reference; implemented from the literature. Cross-validated against # EuclideanCyclicFunctional at the natural initial geometry (u = 0). test_inversive_distance_functional.cpp ``` Run: ```bash ctest --test-dir build -R "InversiveDistance" --output-on-failure ``` --- ## Step 6 — Newton solver Once the functional passes its tests, wire a Newton wrapper into `newton_solver.hpp`: ```cpp inline NewtonResult newton_inversive_distance( ConformalMesh& mesh, std::vector x0, const InversiveDistanceMaps& m, double tol = 1e-8, int max_iter = 200); ``` The body is structurally identical to `newton_euclidean()` — same SimplicialLDLT + SparseQR fallback, same termination test. Only the inner gradient / Hessian calls differ. --- ## Checklist for a new functional - [ ] `code/include/_functional.hpp` compiles - [ ] Limit-case edge-length tests pass at machine precision - [ ] Round-trip identity (init ⇄ length formula) verified - [ ] FD-vs-analytic gradient check passes on triangle, quad strip, tetra - [ ] Cross-validation test against an existing functional at `u = 0` - [ ] Newton wrapper added to `newton_solver.hpp` - [ ] Registered in `code/tests/cgal/CMakeLists.txt` - [ ] `doc/roadmap/java-parity.md` updated (port status or research note) - [ ] `doc/math/references.md` extended with the primary paper(s) - [ ] If this is *new research* beyond Java: add an entry in `doc/roadmap/research-track.md` with citations and acceptance criteria --- ## How to know if it's a port or research Run the local Java-repo check **first** before writing any tutorial doc: ```bash find /Users/tarikmoussa/Desktop/conformallab -iname "**" grep -r "" /Users/tarikmoussa/Desktop/conformallab/src ``` If both return zero matches, the feature is **not** in Java and any C++ implementation is **new research**, not a port. The tutorial framing and the `doc/roadmap/research-track.md` entry should reflect this from day one.