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:
>
> ```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 eq. 4.6) 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<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:

```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<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:

```bash
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.