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

# Tutorial: Porting the Inversive-Distance Functional (Phase 9a.2)

This is a complete, step-by-step example of how to add a new functional
to conformallab++.

> **Note:** This tutorial implements the **vertex-based** inversive-distance
> circle packing of Luo (2004) and Glickenstein (2011).  Despite an earlier
> draft of this document, the Java reference repository `varylab/conformallab`
> does **not** contain an `InversiveDistanceFunctional.java`.  This is a
> from-the-literature implementation, validated by cross-checking against
> the **face-based** circle-packing functional `CPEuclideanFunctional.java`
> (Bobenko-Pinkall-Springborn 2010, see Phase 9a.1) in the tangential
> limit `I_ij = 1` ⇔ `θ_e = 0`.

**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 circle packing parametrises each vertex by a circle of
radius `r_i = exp(u_i)`.  For each edge, the *inversive distance* `I_ij`
between the two adjacent circles is a fixed constant of the edge, computed
once from the initial geometry:

```
I_ij = (ℓ_ij² − r_i² − r_j²) / (2 r_i r_j)                       (Bowers-Stephenson)
```

| Range of `I_ij` | Geometric meaning |
|---|---|
| `I_ij = 1`            | Tangential circles (Koebe-style) |
| `I_ij = cos φ ∈ (0,1)` | Overlapping circles with intersection angle φ |
| `I_ij = 0`            | Orthogonal circles |
| `I_ij ∈ (−1, 0)`      | Disjoint circles (inversive-distance > 1) |

The edge length is then determined by the radii via:

```
ℓ_ij(u)² = exp(2 u_i) + exp(2 u_j) + 2 I_ij exp(u_i + u_j)        (Luo 2004 §3)
```

The angle formula is the same law of cosines (numerically stable
half-tangent form) as the Euclidean functional — only the edge-length
mapping changes.  The gradient

```
∂E/∂u_i = Θ_i − Σ_{T ∋ i} α_i(T)
```

is the standard Yamabe-flow gradient (Luo 2004 Lemma 3.1).

**Primary references:**
- Luo, F. (2004). *Combinatorial Yamabe Flow on Surfaces*.
  Communications in Contemporary Mathematics, 6(5), 765–780.
- Bowers, P. L. & Stephenson, K. (2004). *Uniformizing dessins and
  Belyĭ maps via circle packing*.  Memoirs of the AMS 170(805).
- Glickenstein, D. (2011). *Discrete conformal variations and scalar
  curvature on piecewise flat manifolds*.  J. Differential Geometry
  87(2), 201–238  (analytic Hessian, eq. 4.6).

**Java source:** none — this functional is implemented from the
above papers.  See [doc/roadmap/phases.md](../roadmap/phases.md) §9a for
the historical clarification.

---

## 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 <Eigen/Dense>
#include <cmath>

namespace conformallab {

struct InversiveDistanceMaps {
    CGAL::Surface_mesh<CGAL::Simple_cartesian<double>::Point_3>
        ::Property_map<Edge_index, double>   inv_dist0; ///< I_{ij} per edge
    CGAL::Surface_mesh<CGAL::Simple_cartesian<double>::Point_3>
        ::Property_map<Vertex_index, double> theta_v;   ///< target corner angles Θ_v
    CGAL::Surface_mesh<CGAL::Simple_cartesian<double>::Point_3>
        ::Property_map<Vertex_index, int>    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 <gtest/gtest.h>

// 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<double> 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<double>             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 → ✅