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

// newton_solver.hpp
//
// Phase 4a — Newton solver for all three discrete conformal functionals.
//
// Solves G(x) = 0 where G is the gradient of the discrete conformal energy.
//
// ┌──────────────────────────────────────────────────────────────────────────┐
// │  Gradient sign conventions                                               │
// │    Euclidean / Spherical:  G_v = Θ_v − Σ α_v    (target − actual)      │
// │    HyperIdeal:             G_v = Σ β_v − Θ_v    (actual − target)      │
// │                                                                          │
// │  All solvers use the same Newton step  Δx = −H⁻¹·G                     │
// │                                                                          │
// │  Hessian sign at equilibrium                                             │
// │    Euclidean:   H PSD  → SimplicialLDLT on  H                          │
// │    Spherical:   H NSD  → SimplicialLDLT on −H  (solve (−H)Δx = G)      │
// │    HyperIdeal:  H PSD  → SimplicialLDLT on  H  (analytical H: future)  │
// └──────────────────────────────────────────────────────────────────────────┘
//
// SparseQR fallback:
//   When SimplicialLDLT reports a failure (e.g. singular H on a closed mesh
//   without a pinned vertex), the solver automatically retries with
//   Eigen::SparseQR, which finds the minimum-norm Newton step orthogonal to
//   the null space.  This handles the gauge mode on closed surfaces without
//   requiring the caller to pin a vertex explicitly.
//
// Requires:
//   Eigen::SimplicialLDLT, Eigen::SparseQR (Eigen sparse module)

#include "euclidean_hessian.hpp"
#include "spherical_hessian.hpp"
#include "hyper_ideal_hessian.hpp"
#include "cp_euclidean_functional.hpp"
#include "inversive_distance_functional.hpp"
#include "inversive_distance_hessian.hpp"
#include <Eigen/SparseCholesky>
#include <Eigen/SparseQR>
#include <Eigen/OrderingMethods>
#include <Eigen/Dense>
#include <algorithm>
#include <cmath>

namespace conformallab {

// ── Result ────────────────────────────────────────────────────────────────────

/// Result of `newton_solve(...)` — converged DOF vector + diagnostics.
struct NewtonResult {
    std::vector<double> x;            ///< DOF vector at termination.
    int                 iterations;   ///< Newton steps taken.
    double              grad_inf_norm;///< max |Gᵢ| at termination.
    bool                converged;    ///< `true` iff `grad_inf_norm < tol`.
};

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

namespace detail {

// Solve A·Δx = rhs with SimplicialLDLT; on failure fall back to SparseQR.
// Returns Δx.  ok is set to false only if both solvers fail.
// If fallback_used is non-null, it is set to true iff SparseQR was needed.
inline Eigen::VectorXd solve_with_fallback(
    const Eigen::SparseMatrix<double>& A,
    const Eigen::VectorXd&             rhs,
    bool&                              ok,
    bool*                              fallback_used = nullptr)
{
    if (fallback_used) *fallback_used = false;

    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt(A);
    if (ldlt.info() == Eigen::Success) {
        Eigen::VectorXd dx = ldlt.solve(rhs);
        if (ldlt.info() == Eigen::Success) { ok = true; return dx; }
    }
    // Fallback: SparseQR — handles singular/rank-deficient H (gauge modes).
    if (fallback_used) *fallback_used = true;
    Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> qr(A);
    if (qr.info() == Eigen::Success) {
        Eigen::VectorXd dx = qr.solve(rhs);
        if (qr.info() == Eigen::Success) { ok = true; return dx; }
    }
    ok = false;
    return Eigen::VectorXd::Zero(rhs.size());
}

} // namespace detail

// ── Public linear-system solver (SparseQR fallback) ──────────────────────────
//
// Solve A·x = rhs with Eigen::SimplicialLDLT; if that fails (singular or
// rank-deficient A), retry with Eigen::SparseQR which finds the minimum-norm
// solution orthogonal to the null space.
//
// This is the same primitive used internally by all three Newton solvers.
// Exposing it publicly lets callers (tests, downstream code) reuse the logic
// and — via the optional fallback_used pointer — verify which code path ran.
//
//   fallback_used  – if non-null, set to true iff SparseQR was invoked
//   Returns Eigen::VectorXd::Zero(rhs.size()) if both solvers fail.
/// Solve `A·x = rhs` with the same SimplicialLDLT → SparseQR fallback
/// strategy used inside all three Newton solvers.  If `fallback_used`
/// is non-null, it is set to `true` iff the SparseQR fallback ran.
/// If `ok` is non-null, it is set to `true` iff at least one solver succeeded.
/// Returns `Eigen::VectorXd::Zero(rhs.size())` if both solvers fail.
inline Eigen::VectorXd solve_linear_system(
    const Eigen::SparseMatrix<double>& A,
    const Eigen::VectorXd&             rhs,
    bool*                              fallback_used = nullptr,
    bool*                              ok            = nullptr)
{
    bool ok_local = false;
    auto x = detail::solve_with_fallback(A, rhs, ok_local, fallback_used);
    if (ok) *ok = ok_local;
    return x;
}

namespace detail {   // re-open for the remaining helpers

// Globalised line search for the Newton system G(x) = 0.
//
// Merit function  f(x) = ½‖G(x)‖²₂.  Driving f to its minimum drives the
// residual G to zero; because the merit only depends on ‖G‖ it is sign-agnostic
// and works identically for the convex Euclidean / HyperIdeal / CP / InvDist
// energies and the concave Spherical energy.
//
//   Phase 1 — Newton direction `dx` (satisfies H·dx = −G, so the merit slope
//             ∇f·dx = GᵀH·dx = −‖G‖² < 0 — always a descent direction).
//             Backtrack α ∈ {1, ½, ¼, …} until the Armijo sufficient-decrease
//             condition holds:
//                   ‖G(x + α·dx)‖²  ≤  (1 − 2·c1·α)·‖G(x)‖²
//
//   Phase 2 — if Phase 1 exhausts its halvings, fall back to the steepest-
//             descent direction of the merit, `d_sd = −H·G` (slope
//             ∇f·d_sd = −‖H·G‖² ≤ 0 regardless of the definiteness of H), with
//             the analogous Armijo test:
//                   ‖G(x + α·d_sd)‖²  ≤  ‖G(x)‖² − 2·c1·α·‖d_sd‖²
//
// If neither phase satisfies Armijo, return the best (smallest-residual) point
// visited.  If nothing beat ‖G(x)‖, return x unchanged and set *improved=false,
// so the caller can stop cleanly instead of taking the old divergent full step.
template <typename GradFn>
inline std::vector<double> line_search(
    const std::vector<double>& x,
    const Eigen::VectorXd&     dx,
    const Eigen::VectorXd&     d_sd,
    double                     norm0,
    GradFn&&                   grad_fn,
    bool*                      improved     = nullptr,
    int                        max_halvings = 20,
    double                     c1           = 1e-4)
{
    const int    n        = static_cast<int>(x.size());
    const double norm0_sq = norm0 * norm0;

    std::vector<double> xnew(static_cast<std::size_t>(n));
    std::vector<double> best_x    = x;
    double              best_norm = norm0;

    // Evaluate ‖G(x + α·dir)‖₂, leaving the trial point in `xnew`.
    auto eval = [&](const Eigen::VectorXd& dir, double alpha) -> double {
        for (int i = 0; i < n; ++i)
            xnew[static_cast<std::size_t>(i)] =
                x[static_cast<std::size_t>(i)] + alpha * dir[i];
        auto Gnew = grad_fn(xnew);
        double s = 0.0;
        for (double v : Gnew) s += v * v;
        return std::sqrt(s);
    };

    // ── Phase 1: Newton direction, Armijo backtracking ────────────────────────
    double alpha = 1.0;
    for (int ls = 0; ls < max_halvings; ++ls) {
        double norm_new = eval(dx, alpha);
        if (norm_new < best_norm) { best_norm = norm_new; best_x = xnew; }
        if (norm_new * norm_new <= (1.0 - 2.0 * c1 * alpha) * norm0_sq) {
            if (improved) *improved = true;
            return xnew;
        }
        alpha *= 0.5;
    }

    // ── Phase 2: steepest-descent fallback (−H·G), Armijo backtracking ────────
    const double dsd_sq = d_sd.squaredNorm();
    if (dsd_sq > 0.0) {
        alpha = 1.0;
        for (int ls = 0; ls < max_halvings; ++ls) {
            double thresh_sq = norm0_sq - 2.0 * c1 * alpha * dsd_sq;
            double norm_new  = eval(d_sd, alpha);
            if (norm_new < best_norm) { best_norm = norm_new; best_x = xnew; }
            if (thresh_sq >= 0.0 && norm_new * norm_new <= thresh_sq) {
                if (improved) *improved = true;
                return xnew;
            }
            alpha *= 0.5;
        }
    }

    // ── Both phases failed Armijo — never take the divergent full step. ───────
    //    Return the best point seen; if none improved, stay put and signal stall.
    if (improved) *improved = (best_norm < norm0);
    return best_x;
}

} // namespace detail

// ── Euclidean Newton solver ────────────────────────────────────────────────────

/// Solve the Euclidean discrete conformal problem: find u ∈ ℝ^V such that
/// Σ_{faces adj v} α_v(u) = Θ_v for all vertices v.
///
/// Starting from x0, Newton's method minimises E(u) (the Euclidean DCE energy,
/// which is convex) by iterating  u ← u − H⁻¹·G  with backtracking line search.
/// The Hessian H is the cotangent Laplacian — PSD with one zero eigenvalue on
/// closed surfaces (gauge mode).  A SparseQR fallback handles this automatically.
///
/// \param mesh      Input triangulated surface (edges must carry lambda0 + theta_v).
/// \param x0        Initial DOF vector (length = number of free vertices).
///                  Pass all-zeros for a flat start (typical).
/// \param m         EuclideanMaps: lambda0[e], theta_v[v], v_idx[v] must be set.
///                  Call setup_euclidean_maps() + compute_euclidean_lambda0_from_mesh()
///                  + enforce_gauss_bonnet() before passing here.
/// \param tol       Convergence threshold on max |G_i|.  Default: 1e-8.
/// \param max_iter  Maximum Newton iterations.  Default: 200.
/// \return          NewtonResult{x*, iterations, grad_inf_norm, converged}.
///
/// \note On closed meshes without a pinned vertex, SimplicialLDLT detects the
///       gauge singularity and falls back to SparseQR automatically.
///
/// \see doc/math/discrete-conformal-theory.md §3 for the mathematical background.
inline NewtonResult newton_euclidean(
    ConformalMesh&             mesh,
    std::vector<double>        x0,
    const EuclideanMaps&       m,
    double                     tol      = 1e-8,
    int                        max_iter = 200)
{
    std::vector<double> x = x0;
    const int n = static_cast<int>(x.size());

    NewtonResult res;
    res.converged      = false;
    res.iterations     = 0;
    res.grad_inf_norm  = 0.0;

    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;

    for (int iter = 0; iter < max_iter; ++iter) {
        // ── Gradient ──────────────────────────────────────────────────────────
        auto G_std = euclidean_gradient(mesh, x, m);
        Eigen::Map<const Eigen::VectorXd> G(G_std.data(), n);

        double inf_norm = G.cwiseAbs().maxCoeff();
        if (inf_norm < tol) {
            res.converged     = true;
            res.grad_inf_norm = inf_norm;
            res.iterations    = iter;
            res.x             = x;
            return res;
        }

        // ── Hessian + solve H·Δx = −G  (SparseQR fallback for singular H) ──
        // Cyclic layout (edge DOFs present) → analytic cyclic Hessian, covering
        // the vertex-edge / edge-edge blocks the vertex-only cotangent Laplacian
        // omits.  Vertex-only layout → analytic cotangent Laplacian.
        bool has_edge_dof = false;
        for (auto e : mesh.edges())
            if (m.e_idx[e] >= 0) { has_edge_dof = true; break; }
        auto H = has_edge_dof ? euclidean_hessian_analytic(mesh, x, m)
                              : euclidean_hessian(mesh, x, m);
        bool ok = false;
        Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok);
        if (!ok) break;

        // ── Globalised line search (Armijo + steepest-descent fallback) ───────
        double norm0 = G.norm();
        Eigen::VectorXd d_sd = -(H * G);   // merit-function steepest descent
        bool improved = true;
        x = detail::line_search(x, dx, d_sd, norm0,
            [&](const std::vector<double>& xnew) {
                return euclidean_gradient(mesh, xnew, m);
            }, &improved);
        if (!improved) break;   // line search stalled — stop cleanly

        res.iterations = iter + 1;
    }

    // Report final gradient norm
    auto G_final = euclidean_gradient(mesh, x, m);
    double inf_final = 0.0;
    for (double v : G_final) inf_final = std::max(inf_final, std::abs(v));
    res.grad_inf_norm = inf_final;
    res.x             = x;
    return res;
}

// ── Spherical Newton solver ───────────────────────────────────────────────────

/// Solve the spherical discrete conformal problem: find u ∈ ℝ^V such that
/// Σ_{faces adj v} α_v(u) = Θ_v for all vertices v (genus-0 / sphere-like surfaces).
///
/// The spherical DCE energy is *concave*, so the Hessian H is NSD at the solution.
/// The solver factorises −H (which is PSD) and solves (−H)·Δx = G.
/// A gauge vertex must be pinned (set v_idx = -1) to remove the rotational mode.
///
/// \param mesh      Input triangulated surface, genus 0.
/// \param x0        Initial DOF vector (length = free vertices, excluding gauge_vertex).
///                  All-zeros is a good start.
/// \param m         SphericalMaps: lambda0[e], theta_v[v], v_idx[v], gauge_vertex set.
///                  Call setup_spherical_maps() + compute_spherical_lambda0_from_mesh()
///                  + enforce_gauss_bonnet() (checks Σ(2π-Θ) > 0) before passing here.
/// \param tol       Convergence threshold on max |G_i|.  Default: 1e-8.
/// \param max_iter  Maximum Newton iterations.  Default: 200.
/// \return          NewtonResult{x*, iterations, grad_inf_norm, converged}.
///
/// \note Unlike the Euclidean solver, the spherical solver does NOT need a SparseQR
///       fallback — the gauge vertex pins the null mode directly.
///
/// \see doc/math/geometry-modes.md §Spherical for sign-convention details.
inline NewtonResult newton_spherical(
    ConformalMesh&             mesh,
    std::vector<double>        x0,
    const SphericalMaps&       m,
    double                     tol      = 1e-8,
    int                        max_iter = 200)
{
    std::vector<double> x = x0;
    const int n = static_cast<int>(x.size());

    NewtonResult res;
    res.converged      = false;
    res.iterations     = 0;
    res.grad_inf_norm  = 0.0;

    for (int iter = 0; iter < max_iter; ++iter) {
        // ── Gradient ──────────────────────────────────────────────────────────
        auto G_std = spherical_gradient(mesh, x, m);
        Eigen::Map<const Eigen::VectorXd> G(G_std.data(), n);

        double inf_norm = G.cwiseAbs().maxCoeff();
        if (inf_norm < tol) {
            res.converged     = true;
            res.grad_inf_norm = inf_norm;
            res.iterations    = iter;
            res.x             = x;
            return res;
        }

        // ── Hessian: negate to get PSD; solve (−H)·Δx = G  ──────────────────
        auto H    = spherical_hessian(mesh, x, m);
        auto negH = Eigen::SparseMatrix<double>(-H);
        bool ok   = false;
        Eigen::VectorXd dx = detail::solve_with_fallback(negH, G, ok);
        if (!ok) break;

        // ── Globalised line search (Armijo + steepest-descent fallback) ───────
        // d_sd uses the actual (un-negated) Hessian: ∇f = H·G for f = ½‖G‖².
        double norm0 = G.norm();
        Eigen::VectorXd d_sd = -(H * G);
        bool improved = true;
        x = detail::line_search(x, dx, d_sd, norm0,
            [&](const std::vector<double>& xnew) {
                return spherical_gradient(mesh, xnew, m);
            }, &improved);
        if (!improved) break;

        res.iterations = iter + 1;
    }

    auto G_final = spherical_gradient(mesh, x, m);
    double inf_final = 0.0;
    for (double v : G_final) inf_final = std::max(inf_final, std::abs(v));
    res.grad_inf_norm = inf_final;
    res.x             = x;
    return res;
}

// ── HyperIdeal Newton solver ──────────────────────────────────────────────────

/// Solve the hyper-ideal discrete conformal problem: find (b, a) ∈ ℝ^{V+E} such that
/// Σ β_v(b,a) = Θ_v  and  Σ α_e(b,a) = θ_e  for all vertices v and edges e.
///
/// Used for genus-g surfaces (g ≥ 1) under hyperbolic cone metrics.
/// The energy is *strictly convex* (Springborn 2020, Theorem 1.3), so Newton
/// converges globally from any starting point.
///
/// DOF layout: first V_free entries are vertex variables b_v (hyper-ideal radii),
/// followed by E entries for edge variables a_e (intersection angles).
/// Use assign_hyper_ideal_all_dof_indices(mesh, maps) to set v_idx and e_idx automatically —
/// no vertex needs to be pinned.
///
/// \param mesh      Input triangulated surface, genus g ≥ 1.
/// \param x0        Initial DOF vector (length = V + E).  All-zeros typical.
/// \param m         HyperIdealMaps: lambda0[e], theta_v[v], v_idx[v], e_idx[e] set.
///                  Call setup_hyper_ideal_maps() + compute_hyper_ideal_lambda0_from_mesh().
/// \param tol       Convergence threshold on max |G_i|.  Default: 1e-8.
/// \param max_iter  Maximum Newton iterations.  Default: 200.
/// \param hess_eps  Finite-difference step for Hessian approximation.  Default: 1e-5.
///                  (Phase 9b will replace this with an analytic Hessian.)
/// \return          NewtonResult{x*, iterations, grad_inf_norm, converged}.
///
/// \see Springborn (2020), Theorem 1.3 for the strict convexity proof.
/// \see doc/math/geometry-modes.md §Hyper-ideal for DOF layout details.
inline NewtonResult newton_hyper_ideal(
    ConformalMesh&             mesh,
    std::vector<double>        x0,
    const HyperIdealMaps&      m,
    double                     tol      = 1e-8,
    int                        max_iter = 200,
    double                     hess_eps = 1e-5)
{
    std::vector<double> x = x0;
    const int n = static_cast<int>(x.size());

    NewtonResult res;
    res.converged      = false;
    res.iterations     = 0;
    res.grad_inf_norm  = 0.0;

    for (int iter = 0; iter < max_iter; ++iter) {
        // ── Gradient ──────────────────────────────────────────────────────────
        auto G_std = evaluate_hyper_ideal(mesh, x, m, /*energy=*/false).gradient;
        Eigen::Map<const Eigen::VectorXd> G(G_std.data(), n);

        double inf_norm = G.cwiseAbs().maxCoeff();
        if (inf_norm < tol) {
            res.converged     = true;
            res.grad_inf_norm = inf_norm;
            res.iterations    = iter;
            res.x             = x;
            return res;
        }

        // ── Hessian (block-FD, ~33–1166× faster than full-FD) + solve H·Δx = −G
        auto H  = hyper_ideal_hessian_block_fd_sym(mesh, x, m, hess_eps);
        bool ok = false;
        Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok);
        if (!ok) break;

        // ── Globalised line search (Armijo + steepest-descent fallback) ───────
        double norm0 = G.norm();
        Eigen::VectorXd d_sd = -(H * G);
        bool improved = true;
        x = detail::line_search(x, dx, d_sd, norm0,
            [&](const std::vector<double>& xnew) {
                return evaluate_hyper_ideal(mesh, xnew, m, false).gradient;
            }, &improved);
        if (!improved) break;

        res.iterations = iter + 1;
    }

    auto G_final = evaluate_hyper_ideal(mesh, x, m, false).gradient;
    double inf_final = 0.0;
    for (double v : G_final) inf_final = std::max(inf_final, std::abs(v));
    res.grad_inf_norm = inf_final;
    res.x             = x;
    return res;
}

// ── CP-Euclidean Newton solver (Phase 9a.1) ───────────────────────────────────

/// Solve the CP-Euclidean circle-packing problem: find ρ ∈ ℝ^F such that the
/// per-face angle sums match φ_f at every free face.
///
/// The CP-Euclidean energy (Bobenko-Pinkall-Springborn 2015 §6) is strictly
/// convex on its open domain of validity, so the Hessian H is PSD and the
/// solution is unique up to the gauge mode pinned by `f_idx == −1`.
/// `cp_euclidean_hessian` provides the analytic 2×2-per-edge formula
/// `h_jk = sin θ / (cosh Δρ − cos θ)`; no FD machinery is required.
///
/// \param mesh      Input triangle mesh (closed or with boundary).
/// \param x0        Initial DOF vector (length = number of free faces).
///                  All-zeros is a valid start.
/// \param m         CPEuclideanMaps: f_idx must have one pinned face
///                  (`f_idx[f0] == −1`); theta_e and phi_f set by the caller.
/// \param tol       Convergence threshold on `‖G‖∞`.  Default: 1e-8.
/// \param max_iter  Newton iteration limit.  Default: 200.
/// \return          NewtonResult{x*, iterations, grad_inf_norm, converged}.
///
/// \note  Unlike the Euclidean solver, the CP-Euclidean Hessian is exact
///        (analytic), so the SparseQR fallback only triggers in genuine
///        gauge-singular situations (no pinned face).
/// \see   doc/architecture/phase-9a-validation.md §1 for the BPS-2015 mapping.
inline NewtonResult newton_cp_euclidean(
    ConformalMesh&             mesh,
    std::vector<double>        x0,
    const CPEuclideanMaps&     m,
    double                     tol      = 1e-8,
    int                        max_iter = 200)
{
    std::vector<double> x = x0;
    const int n = static_cast<int>(x.size());

    NewtonResult res;
    res.converged      = false;
    res.iterations     = 0;
    res.grad_inf_norm  = 0.0;

    for (int iter = 0; iter < max_iter; ++iter) {
        auto G_std = cp_euclidean_gradient(mesh, x, m);
        Eigen::Map<const Eigen::VectorXd> G(G_std.data(), n);

        double inf_norm = G.cwiseAbs().maxCoeff();
        if (inf_norm < tol) {
            res.converged     = true;
            res.grad_inf_norm = inf_norm;
            res.iterations    = iter;
            res.x             = x;
            return res;
        }

        auto H = cp_euclidean_hessian(mesh, x, m);
        bool ok = false;
        Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok);
        if (!ok) break;

        double norm0 = G.norm();
        Eigen::VectorXd d_sd = -(H * G);
        bool improved = true;
        x = detail::line_search(x, dx, d_sd, norm0,
            [&](const std::vector<double>& xnew) {
                return cp_euclidean_gradient(mesh, xnew, m);
            }, &improved);
        if (!improved) break;

        res.iterations = iter + 1;
    }

    auto G_final = cp_euclidean_gradient(mesh, x, m);
    double inf_final = 0.0;
    for (double v : G_final) inf_final = std::max(inf_final, std::abs(v));
    res.grad_inf_norm = inf_final;
    res.x             = x;
    return res;
}

// ── Inversive-Distance Newton solver (Phase 9a.2) ─────────────────────────────

/// Solve the inversive-distance circle-packing problem: find u ∈ ℝ^V such that
/// Σ_{faces adj v} α_v(u) = Θ_v at every free vertex (Luo 2004 Lemma 3.1).
///
/// The inversive-distance energy is (locally) strictly convex on the open
/// domain where every triangle satisfies the inequalities.  Luo's 1-form is
/// closed there, so the path-integral energy is well-defined.
///
/// The Hessian is computed by **per-face block finite differences**
/// (`inversive_distance_hessian_block_fd_sym`, same pattern as the HyperIdeal
/// solver after Phase 9b) — O(F) face evaluations instead of the O(n·F) of the
/// full-FD baseline.  An analytic Hessian via Glickenstein 2011 eq. (4.6) is
/// tracked in `doc/roadmap/research-track.md` as Phase 9a.2-analytic.
///
/// \param mesh      Input triangle mesh.
/// \param x0        Initial DOF vector (length = number of free vertices).
/// \param m         InversiveDistanceMaps: v_idx has at least one pinned
///                  vertex; I_e and r0 set by compute_inversive_distance_init.
/// \param tol       Convergence threshold on `‖G‖∞`.  Default: 1e-8.
/// \param max_iter  Newton iteration limit.  Default: 200.
/// \param hess_eps  FD step size for the Hessian.  Default: 1e-5.
/// \return          NewtonResult{x*, iterations, grad_inf_norm, converged}.
///
/// \note  Convergence is sensitive to the initial point: u = 0 is the
///        natural choice when `compute_inversive_distance_init_from_mesh`
///        has been called, since the Bowers-Stephenson identity reconstructs
///        the input edge lengths at u = 0.
inline NewtonResult newton_inversive_distance(
    ConformalMesh&                mesh,
    std::vector<double>           x0,
    const InversiveDistanceMaps&  m,
    double                        tol      = 1e-8,
    int                           max_iter = 200,
    double                        hess_eps = 1e-5)
{
    std::vector<double> x = x0;
    const int n = static_cast<int>(x.size());

    NewtonResult res;
    res.converged      = false;
    res.iterations     = 0;
    res.grad_inf_norm  = 0.0;

    for (int iter = 0; iter < max_iter; ++iter) {
        auto G_std = inversive_distance_gradient(mesh, x, m);
        Eigen::Map<const Eigen::VectorXd> G(G_std.data(), n);

        double inf_norm = G.cwiseAbs().maxCoeff();
        if (inf_norm < tol) {
            res.converged     = true;
            res.grad_inf_norm = inf_norm;
            res.iterations    = iter;
            res.x             = x;
            return res;
        }

        // Hessian (block-FD, ~n/6× faster than full-FD) + solve H·Δx = −G.
        auto H = inversive_distance_hessian_block_fd_sym(mesh, x, m, hess_eps);
        bool ok = false;
        Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok);
        if (!ok) break;

        double norm0 = G.norm();
        Eigen::VectorXd d_sd = -(H * G);
        bool improved = true;
        x = detail::line_search(x, dx, d_sd, norm0,
            [&](const std::vector<double>& xnew) {
                return inversive_distance_gradient(mesh, xnew, m);
            }, &improved);
        if (!improved) break;

        res.iterations = iter + 1;
    }

    auto G_final = inversive_distance_gradient(mesh, x, m);
    double inf_final = 0.0;
    for (double v : G_final) inf_final = std::max(inf_final, std::abs(v));
    res.grad_inf_norm = inf_final;
    res.x             = x;
    return res;
}

} // namespace conformallab