ConformalLabpp/doc/tutorials/add-inversive-distance.md

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:

$ 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 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

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).

struct InversiveDistanceMaps {
    ConformalMesh::Property_map<Vertex_index, int>     v_idx;    // DOF index (−1 = pinned)
    ConformalMesh::Property_map<Vertex_index, double>  theta_v;  // target cone angle
    ConformalMesh::Property_map<Vertex_index, double>  r0;       // initial radius r_i^(0)
    ConformalMesh::Property_map<Edge_index,   double>  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:

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:

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:

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 ℓ.

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:

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:

    # ── 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:

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:

inline NewtonResult newton_inversive_distance(
    ConformalMesh&                  mesh,
    std::vector<double>             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/<name>_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:

find /Users/tarikmoussa/Desktop/conformallab -iname "*<feature>*"
grep -r "<ClassName>" /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.