ConformalLabpp/code/include/inversive_distance_functional.hpp

#pragma once
// Copyright (c) 2024-2026 Tarik Moussa.
// SPDX-License-Identifier: MIT

// inversive_distance_functional.hpp
//
// Phase 9a.2 — Inversive-distance circle-packing functional (Luo 2004).
//
// VERTEX-based circle packing.  Each vertex carries a circle of radius
// r_i = exp(u_i).  The inversive distance I_ij between two adjacent
// circles is a constant of the edge, derived once from the initial
// geometry via Bowers-Stephenson 2004.
//
// This is the FACE-DUAL of CPEuclideanFunctional (Phase 9a.1).  The
// correspondence is I_ij = cos θ_e (Glickenstein 2011 §5).
//
// ┌──────────────────────────────────────────────────────────────────────────┐
// │  Mathematical model                                                       │
// │  ──────────────────                                                       │
// │                                                                           │
// │  Variables:    u_i = log r_i           (per vertex; r_i is the radius)   │
// │  Constants:    I_ij                    (per edge; inversive distance)    │
// │                Θ_v                     (per vertex; target cone angle)   │
// │                                                                           │
// │  Bowers-Stephenson (init from initial geometry):                          │
// │    I_ij = ( ℓ_ij² − r_i² − r_j² ) / ( 2 r_i r_j )                        │
// │                                                                           │
// │  Edge length (Luo 2004 §3, Glickenstein 2011 eq. 2.1):                   │
// │    ℓ_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                              │
// │                                                                           │
// │  Triangle angles:    same half-tangent law of cosines as the              │
// │                      Euclidean functional (numerically stable).           │
// │                                                                           │
// │  Gradient  (Luo 2004 Lemma 3.1):                                          │
// │    ∂E/∂u_v  =  Θ_v  −  Σ_{T ∋ v}  α_v(T)                                 │
// │                                                                           │
// │  Energy:     path integral E(u) = ∫₀¹ ⟨G(tu), u⟩ dt                      │
// │              (Luo's 1-form is closed; we use 10-point Gauss-Legendre      │
// │               quadrature, identical to euclidean_functional.hpp)          │
// │                                                                           │
// │  Hessian:    finite-difference for the MVP port; an analytic form is      │
// │              given in Glickenstein 2011 eq. (4.6) and may be added        │
// │              later for performance.                                       │
// └──────────────────────────────────────────────────────────────────────────┘
//
// Relation to euclidean_functional.hpp
// ────────────────────────────────────
// The two are structurally identical in:
//   • DOF layout (per vertex), DOF index sentinel (−1 = pinned)
//   • Gradient pattern (Θ − Σ α)
//   • Energy via path integral (same Gauss-Legendre constants)
//   • Halfedge convention (h0/h1/h2, source pattern, α opposite-edge)
//
// They differ ONLY in:
//   • Per-edge constant: λ°_ij (log²-length) vs I_ij (inversive distance)
//   • Edge-length formula:
//       Euclidean:           ℓ_ij  = exp((λ°_ij + u_i + u_j) / 2)
//       Inversive distance:  ℓ_ij² = exp(2u_i) + exp(2u_j)
//                                  + 2 I_ij exp(u_i + u_j)
//
// In particular at the tangential limit I_ij = 1 the inversive-distance length
// reduces to (exp(u_i) + exp(u_j))² ⇒ ℓ_ij = r_i + r_j (tangential circles),
// which is *different* from the Euclidean-conformal length even at the same
// initial geometry.  The two functionals describe distinct geometric objects.
//
// Property-map name prefix: "iv:" (vertex) and "ie:" (edge).

#include "conformal_mesh.hpp"
#include "constants.hpp"
#include "gauss_legendre.hpp"
#include "euclidean_geometry.hpp"   // euclidean_angles(λ12, λ23, λ31)
#include <CGAL/boost/graph/iterator.h>
#include <vector>
#include <cmath>
#include <cstdint>
#include <iostream>

namespace conformallab {

// ── Property-map type aliases ────────────────────────────────────────────────

/// Property map vertex → `int`    for the Inversive-Distance functional.
using IDVMapI = ConformalMesh::Property_map<Vertex_index, int>;
/// Property map vertex → `double` for the Inversive-Distance functional.
using IDVMapD = ConformalMesh::Property_map<Vertex_index, double>;
/// Property map edge   → `double` for the Inversive-Distance functional.
using IDEMapD = ConformalMesh::Property_map<Edge_index, double>;

// ── Persistent map bundle ─────────────────────────────────────────────────────

/// Bundle of the four property maps consumed by the Inversive-Distance
/// circle-packing functional (Luo 2004 / Bowers-Stephenson 2004).
struct InversiveDistanceMaps {
    IDVMapI v_idx;       ///< DOF index per vertex (−1 = pinned / u_v = 0)
    IDVMapD theta_v;     ///< target cone angle Θ_v   (default 2π)
    IDVMapD r0;          ///< initial radius r_i^(0)  (default 1)
    IDEMapD I_e;         ///< inversive distance I_ij (per edge, constant)
};

/// Attach the four inversive-distance property maps to `mesh` and
/// return their handles.
///
/// Defaults are intentionally trivial — every real use of this
/// functional must call `compute_inversive_distance_init_from_mesh()`
/// next to populate `r0` and `I_e` from the input geometry.
/// * `v_idx[v] = -1` (all vertices pinned initially)
/// * `theta_v[v] = 2π` (regular interior vertex)
/// * `r0[v] = 1.0` (placeholder)
/// * `I_e[e] = 1.0` (tangential default — overwritten by init step)
///
/// The maps use the `"iv:"` / `"ie:"` prefix so they do not collide
/// with the Euclidean / Spherical / HyperIdeal / CP-Euclidean maps.
inline InversiveDistanceMaps setup_inversive_distance_maps(ConformalMesh& mesh)
{
    InversiveDistanceMaps m;
    m.v_idx   = mesh.add_property_map<Vertex_index, int>   ("iv:idx",   -1     ).first;
    m.theta_v = mesh.add_property_map<Vertex_index, double>("iv:theta", TWO_PI ).first;
    m.r0      = mesh.add_property_map<Vertex_index, double>("iv:r0",    1.0    ).first;
    m.I_e     = mesh.add_property_map<Edge_index,   double>("ie:I",     1.0    ).first;
    return m;
}

/// Assign sequential DOF indices `0..n-1` to every vertex.
///
/// **Note:** this overload assigns indices to ALL vertices unconditionally.
/// Any `v_idx` set before the call is overwritten.  To pin a gauge vertex,
/// either use the two-argument overload below, or set `m.v_idx[v] = -1`
/// **after** this call.  Pinning before this call has no effect.
///
/// For a closed mesh, exactly one pin is required to remove the
/// global rotational mode.
inline int assign_inversive_distance_vertex_dof_indices(ConformalMesh&         mesh,
                                                          InversiveDistanceMaps& m)
{
    int idx = 0;
    for (auto v : mesh.vertices()) m.v_idx[v] = idx++;
    return idx;
}

/// Assign sequential DOF indices to all vertices, pinning `gauge`
/// (`m.v_idx[gauge] = -1`).  Use this overload on closed meshes to fix
/// the rotational gauge mode in a single call.
///
/// \returns  The number of free DOFs assigned (`num_vertices − 1`).
inline int assign_inversive_distance_vertex_dof_indices(ConformalMesh&         mesh,
                                                          InversiveDistanceMaps& m,
                                                          Vertex_index           gauge)
{
    int idx = 0;
    for (auto v : mesh.vertices())
        m.v_idx[v] = (v == gauge) ? -1 : idx++;
    return idx;
}

/// Count the free DOFs (vertices with `v_idx >= 0`).
inline int inversive_distance_dimension(const ConformalMesh&         mesh,
                                          const InversiveDistanceMaps& m)
{
    int dim = 0;
    for (auto v : mesh.vertices()) if (m.v_idx[v] >= 0) ++dim;
    return dim;
}

/// Two-phase initialisation from initial mesh geometry.  Mirrors the
/// role of `compute_lambda0_from_mesh` in the Euclidean functional, but
/// adapted to Luo's vertex-based radius parametrisation.
///
/// **Phase 1.** Pick a positive radius per vertex:
/// \code
///   r_i^(0) = (1/3) · min{ℓ_e : e adjacent to v_i}
/// \endcode
/// This is a heuristic — the user may override `m.r0[v]` for any
/// vertex between `setup_inversive_distance_maps()` and this call.
///
/// **Phase 2.** Compute the per-edge inversive distance via the
/// Bowers-Stephenson 2004 identity:
/// \code
///   I_ij = ( ℓ_ij² − r_i² − r_j² ) / ( 2 r_i r_j )
/// \endcode
///
/// \pre  Every edge has positive 3-D length.
/// \pre  Radii produced in Phase 1 are positive (degenerate isolated
///       vertices fall back to `r_i = 1`).
/// \post Every `I_e[e] > -1` for a valid packing.  The chosen
///       Phase-1 heuristic keeps `I_e > 0` for most real meshes.
inline void compute_inversive_distance_init_from_mesh(ConformalMesh&         mesh,
                                                        InversiveDistanceMaps& m)
{
    // Phase 1: r_i = (1/3) · min adjacent edge length.
    for (auto v : mesh.vertices()) {
        double min_len = std::numeric_limits<double>::infinity();
        for (auto h : CGAL::halfedges_around_target(v, mesh)) {
            auto p1 = mesh.point(mesh.source(h));
            auto p2 = mesh.point(mesh.target(h));
            double dx = p1.x() - p2.x();
            double dy = p1.y() - p2.y();
            double dz = p1.z() - p2.z();
            double len = std::sqrt(dx*dx + dy*dy + dz*dz);
            if (len < min_len) min_len = len;
        }
        m.r0[v] = (std::isfinite(min_len) && min_len > 1e-15)
                ? min_len / 3.0
                : 1.0;
    }

    // Phase 2: I_ij from initial geometry.
    for (auto e : mesh.edges()) {
        auto h  = mesh.halfedge(e);
        auto vi = mesh.source(h);
        auto vj = mesh.target(h);
        auto p1 = mesh.point(vi);
        auto p2 = mesh.point(vj);
        double dx = p1.x() - p2.x();
        double dy = p1.y() - p2.y();
        double dz = p1.z() - p2.z();
        double l2 = dx*dx + dy*dy + dz*dz;
        double ri = m.r0[vi];
        double rj = m.r0[vj];
        m.I_e[e] = (l2 - ri*ri - rj*rj) / (2.0 * ri * rj);
    }
}

// ── Internal helpers ──────────────────────────────────────────────────────────

namespace id_detail {

inline double dof_val(int idx, const std::vector<double>& x) noexcept
{
    return idx >= 0 ? x[static_cast<std::size_t>(idx)] : 0.0;
}

// halfedge_to_index is defined in conformal_mesh.hpp.
inline std::size_t hidx(Halfedge_index h) noexcept { return halfedge_to_index(h); }

// Inversive-distance edge length squared: ℓ² = exp(2u_i) + exp(2u_j) + 2 I r_i r_j
// where r_i = exp(u_i), so:  ℓ² = r_i² + r_j² + 2 I r_i r_j.
// Returns -1 if the result is non-positive (degenerate; the caller skips the face).
inline double edge_length_squared(double u_i, double u_j, double I_ij) noexcept
{
    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

/// Inversive-Distance gradient `G_v = Θ_v − Σ_faces α_v(face)`.  Same
/// half-edge corner-angle storage convention as `euclidean_gradient`.
inline std::vector<double> inversive_distance_gradient(
    const ConformalMesh&             mesh,
    const std::vector<double>&       x,
    const InversiveDistanceMaps&     m)
{
    const int n = inversive_distance_dimension(mesh, m);
    std::vector<double> G(static_cast<std::size_t>(n), 0.0);

    const std::size_t nh = mesh.number_of_halfedges();
    std::vector<double> h_alpha(nh, 0.0);

    // Pass 1 — per face, compute corner angles via the law of cosines.
    // We reuse euclidean_angles(λ12, λ23, λ31) which takes 2·log(ℓ) per edge.
    for (auto f : mesh.faces()) {
        Halfedge_index h0 = mesh.halfedge(f);
        Halfedge_index h1 = mesh.next(h0);
        Halfedge_index h2 = mesh.next(h1);

        Vertex_index v1 = mesh.source(h0);
        Vertex_index v2 = mesh.source(h1);
        Vertex_index v3 = mesh.source(h2);

        Edge_index e12 = mesh.edge(h0);
        Edge_index e23 = mesh.edge(h1);
        Edge_index e31 = mesh.edge(h2);

        double u1 = id_detail::dof_val(m.v_idx[v1], x);
        double u2 = id_detail::dof_val(m.v_idx[v2], x);
        double u3 = id_detail::dof_val(m.v_idx[v3], x);

        double l12sq = id_detail::edge_length_squared(u1, u2, m.I_e[e12]);
        double l23sq = id_detail::edge_length_squared(u2, u3, m.I_e[e23]);
        double l31sq = id_detail::edge_length_squared(u3, u1, m.I_e[e31]);

        // A non-real circle configuration (ℓ² ≤ 0) has no limiting angle — skip it.
        if (l12sq <= 0 || l23sq <= 0 || l31sq <= 0) continue;

        // euclidean_angles expects 2·log(ℓ) per edge — feed log(ℓ²).
        // For a triangle-inequality-violating face euclidean_angles returns the
        // *limiting* angles (π opposite the over-long edge, 0/0 otherwise) with
        // valid=false.  We deliberately do NOT skip on !fa.valid: using those
        // limiting angles is the convex C¹ extension onto the infeasible region,
        // so Newton can pass through a flip instead of stalling (Finding 9 —
        // mirrors the Euclidean/Spherical fix in Finding 1).  The angles come
        // from genuine Euclidean side lengths, so the extension is geometric.
        auto fa = euclidean_angles(std::log(l12sq), std::log(l23sq), std::log(l31sq));

        h_alpha[id_detail::hidx(h0)] = fa.alpha3;
        h_alpha[id_detail::hidx(h1)] = fa.alpha1;
        h_alpha[id_detail::hidx(h2)] = fa.alpha2;
    }

    // Pass 2 — accumulate vertex gradient.
    for (auto v : mesh.vertices()) {
        int iv = m.v_idx[v];
        if (iv < 0) continue;
        double sum_alpha = 0.0;
        for (auto h : CGAL::halfedges_around_target(v, mesh)) {
            if (mesh.is_border(h)) continue;
            sum_alpha += h_alpha[id_detail::hidx(mesh.prev(h))];
        }
        G[static_cast<std::size_t>(iv)] = m.theta_v[v] - sum_alpha;
    }

    return G;
}

/// Per-face contribution to the Inversive-Distance gradient, as a pure
/// 3→3 kernel of the three local vertex DOFs `(u1,u2,u3)` and the three
/// edge inversive distances `(I12,I23,I31)`.  No mesh, no property maps —
/// used by the per-face block-FD Hessian in `inversive_distance_hessian.hpp`.
///
/// Returns the three values this face *adds* to the global gradient at
/// `(v1,v2,v3)`.  Since `G_v = Θ_v − Σ_faces α_v`, the per-face contribution
/// is the **negative** corner angles `(−α₁,−α₂,−α₃)`.  A non-real circle
/// configuration (any `ℓ² ≤ 0`) contributes nothing — exactly mirroring the
/// face-skip (`continue`) in `inversive_distance_gradient`, so the block-FD
/// Hessian and the full-FD Hessian see the same per-face support.
struct IDFaceGradContribs {
    double g1; ///< contribution to G at v₁  (= −α₁)
    double g2; ///< contribution to G at v₂  (= −α₂)
    double g3; ///< contribution to G at v₃  (= −α₃)
};

inline IDFaceGradContribs inversive_distance_face_grad_contribs(
    double u1, double u2, double u3,
    double I12, double I23, double I31)
{
    double l12sq = id_detail::edge_length_squared(u1, u2, I12);
    double l23sq = id_detail::edge_length_squared(u2, u3, I23);
    double l31sq = id_detail::edge_length_squared(u3, u1, I31);

    // Same face-skip as inversive_distance_gradient: a non-real circle
    // configuration has no limiting angle, so the face contributes 0.
    if (l12sq <= 0.0 || l23sq <= 0.0 || l31sq <= 0.0)
        return {0.0, 0.0, 0.0};

    // euclidean_angles deliberately keeps its limiting (valid=false) angles
    // here — the convex C¹ extension — matching the gradient's choice not to
    // skip on !fa.valid.  Corner at v_k = fa.alpha_k (see gradient Pass-2 trace).
    auto fa = euclidean_angles(std::log(l12sq), std::log(l23sq), std::log(l31sq));
    return {-fa.alpha1, -fa.alpha2, -fa.alpha3};
}

/// Inversive-Distance energy `E(u) = ∫₀¹ ⟨G(t·u), u⟩ dt`, evaluated
/// with 10-point Gauss-Legendre (constants shared with `euclidean_energy`).
inline double inversive_distance_energy(
    const ConformalMesh&             mesh,
    const std::vector<double>&       x,
    const InversiveDistanceMaps&     m)
{
    const double* gl_s = gl10_nodes();
    const double* gl_w = gl10_weights();

    const std::size_t n = x.size();
    double E = 0.0;
    std::vector<double> tx(n);
    for (int k = 0; k < 10; ++k) {
        double t  = (1.0 + gl_s[k]) * 0.5;
        double wt = gl_w[k] * 0.5;
        for (std::size_t i = 0; i < n; ++i) tx[i] = t * x[i];
        auto G = inversive_distance_gradient(mesh, tx, m);
        double dot = 0.0;
        for (std::size_t i = 0; i < n; ++i) dot += G[i] * x[i];
        E += wt * dot;
    }
    return E;
}

/// FD gradient check for the Inversive-Distance functional (central diff).
/// Uses the same **relative** error criterion as every other gradient check:
/// `|analytic − fd| / max(1, |analytic|) < tol`.
inline bool gradient_check_inversive_distance(
    const ConformalMesh&             mesh,
    const std::vector<double>&       x,
    const InversiveDistanceMaps&     m,
    double eps = 1e-5,
    double tol = 1e-4)
{
    auto G = inversive_distance_gradient(mesh, x, m);
    const std::size_t n = G.size();
    bool ok = true;

    for (std::size_t i = 0; i < n; ++i) {
        std::vector<double> xp = x, xm = x;
        xp[i] += eps;
        xm[i] -= eps;
        double Ep    = inversive_distance_energy(mesh, xp, m);
        double Em    = inversive_distance_energy(mesh, xm, m);
        double fd    = (Ep - Em) / (2.0 * eps);
        double err   = std::abs(G[i] - fd);
        double scale = std::max(1.0, std::abs(G[i]));
        if (err / scale > tol) {
            std::cerr << "[inversive-distance] FD gradient mismatch at DOF " << i
                      << ": analytic=" << G[i]
                      << "  FD=" << fd
                      << "  rel-err=" << (err / scale) << "\n";
            ok = false;
        }
    }
    return ok;
}

/// Newton equilibrium check: returns `true` iff the gradient at `x`
/// is below `tol` in infinity norm (Σ adj-face angles equal Θ_v).
inline bool is_inversive_distance_equilibrium(
    const ConformalMesh&             mesh,
    const std::vector<double>&       x,
    const InversiveDistanceMaps&     m,
    double tol = 1e-8)
{
    auto G = inversive_distance_gradient(mesh, x, m);
    for (double g : G)
        if (std::abs(g) > tol) return false;
    return true;
}

} // namespace conformallab