1bf1defc70
Builds on the newton_core refactor; all three findings live in exactly that code.
I1 (test-coverage): add NewtonStatus { Converged, MaxIterations,
LinearSolverFailed, LineSearchStalled } and a `status` field to NewtonResult, so
the three non-convergent exits that `converged == false` previously conflated are
now distinguishable. `converged` is kept (== status==Converged) for back-compat;
+ to_string(NewtonStatus) for logs/tests. newton_core sets the status at each
exit point (centralised by the H2 refactor).
H1 (test-coverage): set res.iterations explicitly at the LinearSolverFailed and
LineSearchStalled breaks (= completed steps), instead of relying on the last
successful iteration's stale value.
N7 (numerical-stability): surface linear-algebra conditioning in the result —
`sparse_qr_fallback_used` (any iteration fell back to SparseQR) and
`min_ldlt_pivot` (smallest |Dᵢᵢ| of the last LDLT, a cheap near-singularity
proxy). This catches the silent case the audit flagged: on a gauge-singular
Hessian SimplicialLDLT "succeeds" with a ~0 pivot and no fallback fires — now
min_ldlt_pivot is tiny and observable.
Tests (+3 synthetic newton_core status/diagnostic tests; +1 assertion on the
closed-mesh-no-pin fallback test). All purely additive — no control-flow or
convergence behaviour changes; 301/301 CGAL tests pass incl. Java parity.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
642 lines
32 KiB
C++
642 lines
32 KiB
C++
#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 ────────────────────────────────────────────────────────────────────
|
||
|
||
/// Why a Newton solve terminated (I1 audit). Distinguishes the three
|
||
/// non-convergent exits that `converged == false` previously conflated.
|
||
enum class NewtonStatus {
|
||
Converged, ///< `grad_inf_norm < tol` — success.
|
||
MaxIterations, ///< hit `max_iter` without reaching `tol`.
|
||
LinearSolverFailed, ///< both LDLT and SparseQR failed on H·Δx = −G.
|
||
LineSearchStalled ///< the line search could not reduce the residual.
|
||
};
|
||
|
||
/// Human-readable name for a `NewtonStatus` (for logs / test messages).
|
||
inline const char* to_string(NewtonStatus s) noexcept
|
||
{
|
||
switch (s) {
|
||
case NewtonStatus::Converged: return "Converged";
|
||
case NewtonStatus::MaxIterations: return "MaxIterations";
|
||
case NewtonStatus::LinearSolverFailed: return "LinearSolverFailed";
|
||
case NewtonStatus::LineSearchStalled: return "LineSearchStalled";
|
||
}
|
||
return "Unknown";
|
||
}
|
||
|
||
/// Result of a Newton solve — converged DOF vector + diagnostics.
|
||
struct NewtonResult {
|
||
std::vector<double> x; ///< DOF vector at termination.
|
||
int iterations; ///< Newton steps actually completed.
|
||
double grad_inf_norm;///< max |Gᵢ| at termination.
|
||
bool converged; ///< `true` iff `status == Converged`.
|
||
|
||
// I1: why it stopped (more informative than the `converged` bool alone).
|
||
NewtonStatus status = NewtonStatus::MaxIterations;
|
||
|
||
// N7: linear-algebra conditioning diagnostics.
|
||
bool sparse_qr_fallback_used = false; ///< any iteration fell back to SparseQR.
|
||
double min_ldlt_pivot = 0.0; ///< smallest |Dᵢᵢ| of the last LDLT
|
||
///< factorisation (≈ near-singularity
|
||
///< proxy); 0 if LDLT never succeeded.
|
||
};
|
||
|
||
// ── 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.
|
||
//
|
||
// B2 (api-performance audit): the gradient computed at the accepted trial point
|
||
// is normally discarded (only its norm is kept), forcing the next Newton
|
||
// iteration to recompute G(x) from scratch. When `accepted_grad` is non-null it
|
||
// receives the gradient vector at the returned point, so the caller can reuse it
|
||
// as the next iteration's convergence gradient. It is only written when the
|
||
// returned point differs from `x` (i.e. `*improved == true`); on a pure stall
|
||
// (`*improved == false`) the caller must recompute G(x) itself.
|
||
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,
|
||
std::vector<double>* accepted_grad = 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;
|
||
|
||
std::vector<double> trial_grad; // gradient at the most recent trial point
|
||
std::vector<double> best_grad; // gradient at best_x (B2)
|
||
|
||
// Evaluate ‖G(x + α·dir)‖₂, leaving the trial point in `xnew` and the
|
||
// gradient there in `trial_grad`.
|
||
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];
|
||
trial_grad = grad_fn(xnew);
|
||
double s = 0.0;
|
||
for (double v : trial_grad) 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; best_grad = trial_grad; }
|
||
if (norm_new * norm_new <= (1.0 - 2.0 * c1 * alpha) * norm0_sq) {
|
||
if (improved) *improved = true;
|
||
if (accepted_grad) *accepted_grad = trial_grad;
|
||
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; best_grad = trial_grad; }
|
||
if (thresh_sq >= 0.0 && norm_new * norm_new <= thresh_sq) {
|
||
if (improved) *improved = true;
|
||
if (accepted_grad) *accepted_grad = trial_grad;
|
||
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.
|
||
const bool any_improvement = (best_norm < norm0);
|
||
if (improved) *improved = any_improvement;
|
||
if (accepted_grad && any_improvement) *accepted_grad = best_grad;
|
||
return best_x;
|
||
}
|
||
|
||
// ── Cached linear solver (B4) ─────────────────────────────────────────────────
|
||
//
|
||
// solve_with_fallback rebuilds a fresh SimplicialLDLT every call, re-running the
|
||
// symbolic analysis (fill-reducing reordering / COLAMD) each Newton iteration —
|
||
// even though the Hessian sparsity pattern is iteration-invariant for a fixed
|
||
// mesh + DOF layout. This object keeps one persistent LDLT and runs
|
||
// `analyzePattern` once, then only `factorize` per iteration.
|
||
//
|
||
// FD-built Hessians (hyper-ideal, inversive-distance) prune near-zero triplets,
|
||
// so their pattern can shift slightly between iterations; we detect that via a
|
||
// change in `nonZeros()` and re-analyze, which keeps the cache correct and
|
||
// self-healing. The SparseQR fallback for singular/rank-deficient systems is
|
||
// preserved verbatim (built fresh on the rare path). The fallback semantics are
|
||
// identical to `solve_with_fallback`, so behaviour is unchanged.
|
||
struct NewtonLinearSolver {
|
||
Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt;
|
||
bool analyzed = false;
|
||
Eigen::Index last_nnz = -1;
|
||
bool fallback_used = false; ///< true iff the last solve used SparseQR
|
||
double last_min_pivot = 0.0; ///< N7: smallest |Dᵢᵢ| of the last
|
||
///< successful LDLT (0 if it fell back)
|
||
|
||
Eigen::VectorXd solve(const Eigen::SparseMatrix<double>& A,
|
||
const Eigen::VectorXd& rhs,
|
||
bool& ok)
|
||
{
|
||
fallback_used = false;
|
||
last_min_pivot = 0.0;
|
||
if (!analyzed || A.nonZeros() != last_nnz) {
|
||
ldlt.analyzePattern(A);
|
||
analyzed = true;
|
||
last_nnz = A.nonZeros();
|
||
}
|
||
ldlt.factorize(A);
|
||
if (ldlt.info() == Eigen::Success) {
|
||
Eigen::VectorXd x = ldlt.solve(rhs);
|
||
if (ldlt.info() == Eigen::Success) {
|
||
ok = true;
|
||
// N7: smallest |Dᵢᵢ| in the LDLᵀ diagonal — a cheap
|
||
// near-singularity proxy (small ⇒ ill-conditioned step).
|
||
const auto D = ldlt.vectorD();
|
||
last_min_pivot = (D.size() > 0) ? D.cwiseAbs().minCoeff() : 0.0;
|
||
return x;
|
||
}
|
||
}
|
||
// Fallback: SparseQR — handles singular/rank-deficient A (gauge modes).
|
||
fallback_used = true;
|
||
Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> qr(A);
|
||
if (qr.info() == Eigen::Success) {
|
||
Eigen::VectorXd x = qr.solve(rhs);
|
||
if (qr.info() == Eigen::Success) { ok = true; return x; }
|
||
}
|
||
ok = false;
|
||
return Eigen::VectorXd::Zero(rhs.size());
|
||
}
|
||
};
|
||
|
||
// ── Unified Newton core (H2) ──────────────────────────────────────────────────
|
||
//
|
||
// All five DCE Newton solvers (Euclidean, Spherical, HyperIdeal, CP-Euclidean,
|
||
// Inversive-Distance) ran a near-identical loop that differed only in (a) the
|
||
// gradient function, (b) the Hessian function, and (c) the spherical sign quirk
|
||
// (the concave spherical energy factorises −H). This template captures that one
|
||
// loop so a fix lands once instead of five times.
|
||
//
|
||
// `grad(x) -> std::vector<double>` : the energy gradient G(x).
|
||
// `hess(x) -> Eigen::SparseMatrix<double>`: the TRUE Hessian H(x) (un-negated).
|
||
// `concave` : if true, solve (−H)·Δx = G (spherical);
|
||
// otherwise H·Δx = −G. Both are the same
|
||
// Newton system; `concave` only selects
|
||
// the PSD matrix actually factorised.
|
||
//
|
||
// Folds in B2 (reuse the line-search gradient as the next convergence gradient),
|
||
// B4 (cached symbolic factorization), and B5 (no dead solver variable).
|
||
template <typename GradFn, typename HessFn>
|
||
inline NewtonResult newton_core(
|
||
std::vector<double> x,
|
||
GradFn&& grad,
|
||
HessFn&& hess,
|
||
bool concave,
|
||
double tol,
|
||
int max_iter)
|
||
{
|
||
const int n = static_cast<int>(x.size());
|
||
|
||
NewtonResult res;
|
||
res.converged = false;
|
||
res.iterations = 0;
|
||
res.grad_inf_norm = 0.0;
|
||
|
||
NewtonLinearSolver solver;
|
||
|
||
// B2: gradient carried across the loop; the first one is the only "extra"
|
||
// evaluation — every later iteration reuses the line-search gradient.
|
||
std::vector<double> G_std = grad(x);
|
||
|
||
for (int iter = 0; iter < max_iter; ++iter) {
|
||
Eigen::Map<const Eigen::VectorXd> G(G_std.data(), n);
|
||
|
||
const double inf_norm = G.cwiseAbs().maxCoeff();
|
||
if (inf_norm < tol) {
|
||
res.converged = true;
|
||
res.status = NewtonStatus::Converged; // I1
|
||
res.grad_inf_norm = inf_norm;
|
||
res.iterations = iter;
|
||
res.x = x;
|
||
return res;
|
||
}
|
||
|
||
Eigen::SparseMatrix<double> H = hess(x);
|
||
|
||
// Newton system: factor the PSD matrix (H for convex, −H for concave).
|
||
Eigen::SparseMatrix<double> A = concave ? Eigen::SparseMatrix<double>(-H) : H;
|
||
Eigen::VectorXd rhs = concave ? Eigen::VectorXd(G) : Eigen::VectorXd(-G);
|
||
|
||
bool ok = false;
|
||
Eigen::VectorXd dx = solver.solve(A, rhs, ok);
|
||
// N7: accumulate conditioning diagnostics from this solve.
|
||
res.sparse_qr_fallback_used = res.sparse_qr_fallback_used || solver.fallback_used;
|
||
res.min_ldlt_pivot = solver.last_min_pivot;
|
||
if (!ok) { // I1/H1: linear solver failed this step
|
||
res.iterations = iter; // (iter steps actually completed)
|
||
res.status = NewtonStatus::LinearSolverFailed;
|
||
break;
|
||
}
|
||
|
||
// Globalised line search (Armijo + steepest-descent fallback).
|
||
// d_sd = −(H·G) is the merit-function steepest descent for f = ½‖G‖²,
|
||
// using the TRUE (un-negated) Hessian for both signs.
|
||
const double norm0 = G.norm();
|
||
Eigen::VectorXd d_sd = -(H * G);
|
||
bool improved = true;
|
||
std::vector<double> G_next;
|
||
x = detail::line_search(x, dx, d_sd, norm0, grad, &improved, &G_next);
|
||
if (!improved) { // I1/H1: line search stalled — stop cleanly
|
||
res.iterations = iter;
|
||
res.status = NewtonStatus::LineSearchStalled;
|
||
break;
|
||
}
|
||
|
||
G_std = std::move(G_next); // B2: reuse for the next convergence check
|
||
res.iterations = iter + 1;
|
||
}
|
||
|
||
// If we fell out of the loop without converging or breaking, it was max_iter
|
||
// (the default `res.status` is already MaxIterations; a break has set its own).
|
||
|
||
// Final gradient norm (reached only on max-iter / break — recomputed to be
|
||
// correct after a stalled or failed step).
|
||
auto G_final = grad(x);
|
||
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 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)
|
||
{
|
||
// Layout is loop-invariant — decide the Hessian variant once (B3): cyclic
|
||
// layout (edge DOFs present) → analytic cyclic Hessian, which covers the
|
||
// vertex-edge / edge-edge blocks the vertex-only cotangent Laplacian omits.
|
||
bool has_edge_dof = false;
|
||
for (auto e : mesh.edges())
|
||
if (m.e_idx[e] >= 0) { has_edge_dof = true; break; }
|
||
|
||
return detail::newton_core(
|
||
std::move(x0),
|
||
[&](const std::vector<double>& xc) { return euclidean_gradient(mesh, xc, m); },
|
||
[&](const std::vector<double>& xc) {
|
||
return has_edge_dof ? euclidean_hessian_analytic(mesh, xc, m)
|
||
: euclidean_hessian(mesh, xc, m);
|
||
},
|
||
/*concave=*/false, tol, max_iter);
|
||
}
|
||
|
||
// ── 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)
|
||
{
|
||
// Concave energy: the Hessian H is NSD, so newton_core factors −H (PSD) and
|
||
// solves (−H)·Δx = G — the same Newton system, just the sign that keeps LDLT
|
||
// well-defined. The merit steepest-descent inside still uses the true H.
|
||
return detail::newton_core(
|
||
std::move(x0),
|
||
[&](const std::vector<double>& xc) { return spherical_gradient(mesh, xc, m); },
|
||
[&](const std::vector<double>& xc) { return spherical_hessian(mesh, xc, m); },
|
||
/*concave=*/true, tol, max_iter);
|
||
}
|
||
|
||
// ── 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.)
|
||
/// \param clamp Vertex-scale floor mode (N3 audit). Default `HardJava`
|
||
/// reproduces the Java oracle's hard `b<0 → floor` snap
|
||
/// (C⁰, parity-faithful). `SmoothBarrier` uses the C¹
|
||
/// softplus floor — smoother near the feasibility boundary
|
||
/// but deviates from the Java golden values. See
|
||
/// `HyperIdealScaleClamp` in hyper_ideal_functional.hpp.
|
||
/// \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,
|
||
HyperIdealScaleClamp clamp = HyperIdealScaleClamp::HardJava)
|
||
{
|
||
// block-FD Hessian (~33–1166× faster than full-FD); `clamp` selects the
|
||
// vertex-scale floor mode (N3). Convex energy → factor H directly.
|
||
return detail::newton_core(
|
||
std::move(x0),
|
||
[&](const std::vector<double>& xc) {
|
||
return evaluate_hyper_ideal(mesh, xc, m, /*energy=*/false, /*grad=*/true, clamp).gradient;
|
||
},
|
||
[&](const std::vector<double>& xc) {
|
||
return hyper_ideal_hessian_block_fd_sym(mesh, xc, m, hess_eps, clamp);
|
||
},
|
||
/*concave=*/false, tol, max_iter);
|
||
}
|
||
|
||
// ── 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)
|
||
{
|
||
// Convex energy with an exact analytic Hessian (BPS-2015 2×2-per-edge).
|
||
return detail::newton_core(
|
||
std::move(x0),
|
||
[&](const std::vector<double>& xc) { return cp_euclidean_gradient(mesh, xc, m); },
|
||
[&](const std::vector<double>& xc) { return cp_euclidean_hessian(mesh, xc, m); },
|
||
/*concave=*/false, tol, max_iter);
|
||
}
|
||
|
||
// ── 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)
|
||
{
|
||
// Convex energy; block-FD Hessian (~n/6× faster than full-FD).
|
||
return detail::newton_core(
|
||
std::move(x0),
|
||
[&](const std::vector<double>& xc) { return inversive_distance_gradient(mesh, xc, m); },
|
||
[&](const std::vector<double>& xc) {
|
||
return inversive_distance_hessian_block_fd_sym(mesh, xc, m, hess_eps);
|
||
},
|
||
/*concave=*/false, tol, max_iter);
|
||
}
|
||
|
||
} // namespace conformallab
|