135bcf0bba
Implement Phase-Session P1 quick wins (4 independent additions):
9h.1: Add --tol and --max-iter CLI options to conformallab_core
- Newton solver tolerance [default 1e-8]
- Newton iteration limit [default 200]
- Thread both through run_euclidean / run_spherical / run_hyper_ideal
- Update CLI parameter table in documentation
9h.2: Add -g cp_euclidean and -g inversive_distance geometry routes
- run_cp_euclidean() & run_inversive_distance() pipelines (~60 lines each)
- Face-based DOF assignment for CP-Euclidean
- Vertex-based DOF assignment for Inversive-Distance
- Both integrated into CLI geometry validator (IsMember)
9g.1: Create conformal_quality.hpp with validation measures
- IsothermicityMeasure: metric anisotropy (conformality deviation)
- DiscreteConformalEquivalenceMeasure: length-cross-ratio residuals
- FlippedTriangles: detects inverted/degenerate triangles
- LengthCrossRatio: discrete conformal invariant computation
- ConvergenceUtility: aggregated convergence statistics (max/mean/sum)
- Ported from Java: plugin/visualizer + convergence utilities
- Includes sanity tests validating finite outputs on valid layouts
9d.3: Create stereographic_layout.hpp for S² → ℂ projection
- Stereographic projection from north pole: S² → ℂ ∪ {∞}
- Inverse projection: ℂ → S² for round-trip validation
- Möbius centring: centres the 2-D point cloud at origin
- stereographic_layout(Layout3D) -> Layout2D conversion
- Round-trip tests: south pole, equator, random sphere points
- Tests: projection/inverse consistency, north pole handling
Test results: 336/336 CGAL tests pass (272 pre-existing + 64 new from all phases)
- conformal_quality.cpp: 13 new tests (measures, isothermic, dce, convergence)
- stereographic_layout.cpp: 10 new tests (projection, inverse, round-trip, layout)
Co-Authored-By: Claude Haiku 4.5 <noreply@anthropic.com>
844 lines
39 KiB
C++
844 lines
39 KiB
C++
// Copyright (c) 2024-2026 Tarik Moussa.
|
||
// SPDX-License-Identifier: MIT
|
||
|
||
// test_newton_solver.cpp
|
||
//
|
||
// Phase 4 — Newton solver tests.
|
||
//
|
||
// Design principle:
|
||
// We test convergence to a KNOWN equilibrium. For the spherical tetrahedron
|
||
// x* = 0 is built-in (G(0) ≈ 0 by construction). For Euclidean meshes we
|
||
// use "natural theta": set theta_v[v] = actual angle sum at x=0, which makes
|
||
// x* = 0 the exact equilibrium by definition.
|
||
// For HyperIdeal we use the same "natural target" trick at a valid base point
|
||
// (b=1.0, a=0.5), since x=0 is degenerate for the HyperIdeal functional.
|
||
//
|
||
// Tests:
|
||
// Spherical:
|
||
// 1. Converges from x=[-0.2,...] to x*=0 (spherical tetrahedron).
|
||
// 2. Converges in few iterations (quadratic convergence near equilibrium).
|
||
// 3. Converges from a large perturbation x=[-0.5,...].
|
||
// 4. Result fields (x.size, grad_inf_norm) are self-consistent.
|
||
//
|
||
// Euclidean:
|
||
// 5. Converges (triangle, 1 pinned vertex, natural theta).
|
||
// 6. Converges (quad strip, 1 pinned vertex, natural theta).
|
||
// 7. Converges with explicitly chosen mixed pinned/variable layout.
|
||
//
|
||
// HyperIdeal:
|
||
// 8. Converges on triangle (all variable, natural targets).
|
||
// 9. Result fields self-consistent.
|
||
// 10. Converges on tetrahedron (10 DOFs, larger mesh).
|
||
// 11. SparseQR: result consistent (valid starting region).
|
||
//
|
||
// SparseQR fallback (direct unit tests):
|
||
// 12. solve_linear_system recovers correct solution on rank-deficient matrix.
|
||
// 13. solve_linear_system sets fallback_used=true on a singular matrix.
|
||
// 14. Euclidean Newton on closed tetrahedron (no pinned vertex) converges
|
||
// via SparseQR gauge-mode handling.
|
||
|
||
#include "conformal_mesh.hpp"
|
||
#include "mesh_builder.hpp"
|
||
#include "euclidean_functional.hpp"
|
||
#include "spherical_functional.hpp"
|
||
#include "hyper_ideal_functional.hpp"
|
||
#include "newton_solver.hpp"
|
||
#include <Eigen/Dense>
|
||
#include <gtest/gtest.h>
|
||
#include <cmath>
|
||
#include <vector>
|
||
|
||
using namespace conformallab;
|
||
|
||
// ────────────────────────────────────────────────────────────────────────────
|
||
// Helper: set theta_v[v] = actual angle sum at x=0 for each variable vertex.
|
||
//
|
||
// Requires that DOF indices have already been assigned (v_idx populated).
|
||
// Uses euclidean_gradient directly so the formula is exact and consistent.
|
||
// ────────────────────────────────────────────────────────────────────────────
|
||
static void set_natural_euclidean_theta(ConformalMesh& mesh, EuclideanMaps& maps, int n)
|
||
{
|
||
std::vector<double> x0(static_cast<std::size_t>(n), 0.0);
|
||
// G[iv] = theta_v[v] - sum_alpha(x=0)
|
||
// so sum_alpha(x=0) = theta_v[v] - G[iv]
|
||
auto G = euclidean_gradient(mesh, x0, maps);
|
||
for (auto v : mesh.vertices()) {
|
||
int iv = maps.v_idx[v];
|
||
if (iv < 0) continue;
|
||
maps.theta_v[v] -= G[static_cast<std::size_t>(iv)];
|
||
}
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Spherical 1 — Converges from moderate perturbation
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(NewtonSolver, Spherical_ConvergesFromPerturbation)
|
||
{
|
||
auto mesh = make_spherical_tetrahedron();
|
||
auto maps = setup_spherical_maps(mesh);
|
||
compute_spherical_lambda0_from_mesh(mesh, maps);
|
||
int n = assign_spherical_vertex_dof_indices(mesh, maps);
|
||
|
||
std::vector<double> x0(static_cast<std::size_t>(n), -0.2);
|
||
auto res = newton_spherical(mesh, x0, maps, /*tol=*/1e-8, /*max_iter=*/50);
|
||
|
||
EXPECT_TRUE(res.converged)
|
||
<< "Newton (spherical) should converge; grad_inf_norm = " << res.grad_inf_norm;
|
||
EXPECT_LT(res.grad_inf_norm, 1e-8);
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Spherical 2 — Quadratic convergence: few iterations suffice
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(NewtonSolver, Spherical_FewIterations)
|
||
{
|
||
auto mesh = make_spherical_tetrahedron();
|
||
auto maps = setup_spherical_maps(mesh);
|
||
compute_spherical_lambda0_from_mesh(mesh, maps);
|
||
int n = assign_spherical_vertex_dof_indices(mesh, maps);
|
||
|
||
std::vector<double> x0(static_cast<std::size_t>(n), -0.2);
|
||
auto res = newton_spherical(mesh, x0, maps, /*tol=*/1e-8, /*max_iter=*/50);
|
||
|
||
EXPECT_LE(res.iterations, 20)
|
||
<< "Newton should converge in ≤ 20 iterations; took " << res.iterations;
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Spherical 3 — Large perturbation: global convergence via line search
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(NewtonSolver, Spherical_ConvergesFromLargePerturbation)
|
||
{
|
||
auto mesh = make_spherical_tetrahedron();
|
||
auto maps = setup_spherical_maps(mesh);
|
||
compute_spherical_lambda0_from_mesh(mesh, maps);
|
||
int n = assign_spherical_vertex_dof_indices(mesh, maps);
|
||
|
||
std::vector<double> x0(static_cast<std::size_t>(n), -0.5);
|
||
auto res = newton_spherical(mesh, x0, maps, /*tol=*/1e-8, /*max_iter=*/100);
|
||
|
||
EXPECT_TRUE(res.converged)
|
||
<< "Newton (spherical, large perturbation) should converge; "
|
||
"grad_inf_norm = " << res.grad_inf_norm;
|
||
EXPECT_LT(res.grad_inf_norm, 1e-8);
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Spherical 4 — Result fields are self-consistent
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(NewtonSolver, Spherical_ResultFieldsConsistent)
|
||
{
|
||
auto mesh = make_spherical_tetrahedron();
|
||
auto maps = setup_spherical_maps(mesh);
|
||
compute_spherical_lambda0_from_mesh(mesh, maps);
|
||
int n = assign_spherical_vertex_dof_indices(mesh, maps);
|
||
|
||
std::vector<double> x0(static_cast<std::size_t>(n), -0.1);
|
||
auto res = newton_spherical(mesh, x0, maps, /*tol=*/1e-8);
|
||
|
||
EXPECT_EQ(static_cast<int>(res.x.size()), n);
|
||
|
||
// Reported grad_inf_norm must match re-computed gradient at res.x
|
||
auto G = spherical_gradient(mesh, res.x, maps);
|
||
double actual_inf = 0.0;
|
||
for (double v : G) actual_inf = std::max(actual_inf, std::abs(v));
|
||
EXPECT_NEAR(actual_inf, res.grad_inf_norm, 1e-10);
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Euclidean 1 — Triangle with 1 pinned vertex + natural theta
|
||
//
|
||
// make_triangle(): 3 vertices. Pin v0 → 2 free DOFs.
|
||
// With natural theta: x* = 0 (G(0) = 0 by construction).
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(NewtonSolver, Euclidean_ConvergesTrianglePinned)
|
||
{
|
||
auto mesh = make_triangle();
|
||
auto maps = setup_euclidean_maps(mesh);
|
||
compute_euclidean_lambda0_from_mesh(mesh, maps);
|
||
|
||
// Pin first vertex, assign sequential DOFs to the other two
|
||
auto vit = mesh.vertices().begin();
|
||
Vertex_index v0 = *vit++;
|
||
maps.v_idx[v0] = -1; // pinned
|
||
int idx = 0;
|
||
for (; vit != mesh.vertices().end(); ++vit)
|
||
maps.v_idx[*vit] = idx++;
|
||
int n = idx; // = 2
|
||
|
||
set_natural_euclidean_theta(mesh, maps, n);
|
||
|
||
std::vector<double> x0(static_cast<std::size_t>(n), -0.1);
|
||
auto res = newton_euclidean(mesh, x0, maps, /*tol=*/1e-8, /*max_iter=*/50);
|
||
|
||
EXPECT_TRUE(res.converged)
|
||
<< "Newton (Euclidean, triangle, pinned) should converge; "
|
||
"grad_inf_norm = " << res.grad_inf_norm;
|
||
EXPECT_LT(res.grad_inf_norm, 1e-8);
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Euclidean 2 — Quad strip with 1 pinned vertex + natural theta
|
||
//
|
||
// make_quad_strip(): 4 vertices, 2 faces. Pin v0 → 3 free DOFs.
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(NewtonSolver, Euclidean_ConvergesQuadStripPinned)
|
||
{
|
||
auto mesh = make_quad_strip();
|
||
auto maps = setup_euclidean_maps(mesh);
|
||
compute_euclidean_lambda0_from_mesh(mesh, maps);
|
||
|
||
// Pin first vertex
|
||
auto vit = mesh.vertices().begin();
|
||
Vertex_index v0 = *vit++;
|
||
maps.v_idx[v0] = -1;
|
||
int idx = 0;
|
||
for (; vit != mesh.vertices().end(); ++vit)
|
||
maps.v_idx[*vit] = idx++;
|
||
int n = idx; // = 3
|
||
|
||
set_natural_euclidean_theta(mesh, maps, n);
|
||
|
||
std::vector<double> x0(static_cast<std::size_t>(n), -0.15);
|
||
auto res = newton_euclidean(mesh, x0, maps, /*tol=*/1e-8, /*max_iter=*/50);
|
||
|
||
EXPECT_TRUE(res.converged)
|
||
<< "Newton (Euclidean, quad strip, pinned) should converge; "
|
||
"grad_inf_norm = " << res.grad_inf_norm;
|
||
EXPECT_LT(res.grad_inf_norm, 1e-8);
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Euclidean 3 — Mixed pinned layout: explicit vertex assignment
|
||
//
|
||
// Quad strip: v0 pinned, v1/v2/v3 free.
|
||
// Natural theta set AFTER DOF assignment so that x* = 0 is the equilibrium
|
||
// for the free vertices (with v0 fixed at u0=0).
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(NewtonSolver, Euclidean_ConvergesMixedPinned)
|
||
{
|
||
auto mesh = make_quad_strip();
|
||
auto maps = setup_euclidean_maps(mesh);
|
||
compute_euclidean_lambda0_from_mesh(mesh, maps);
|
||
|
||
// Explicitly assign DOF indices
|
||
auto vit = mesh.vertices().begin();
|
||
Vertex_index v0 = *vit++;
|
||
Vertex_index v1 = *vit++;
|
||
Vertex_index v2 = *vit++;
|
||
Vertex_index v3 = *vit;
|
||
|
||
maps.v_idx[v0] = -1; // pinned at u0 = 0
|
||
maps.v_idx[v1] = 0;
|
||
maps.v_idx[v2] = 1;
|
||
maps.v_idx[v3] = 2;
|
||
const int n = 3;
|
||
|
||
// Set natural theta AFTER pinning so that x* = [0,0,0] is the equilibrium
|
||
set_natural_euclidean_theta(mesh, maps, n);
|
||
|
||
std::vector<double> x0 = {-0.1, -0.15, -0.05};
|
||
auto res = newton_euclidean(mesh, x0, maps, /*tol=*/1e-8, /*max_iter=*/50);
|
||
|
||
EXPECT_TRUE(res.converged)
|
||
<< "Newton (Euclidean, mixed pinned) should converge; "
|
||
"grad_inf_norm = " << res.grad_inf_norm;
|
||
EXPECT_LT(res.grad_inf_norm, 1e-8);
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// Helper: set HyperIdeal target angles to actual sums at a non-degenerate
|
||
// base point (b_base, a_base), making that point the equilibrium x*.
|
||
//
|
||
// Note: x = 0 is degenerate for the HyperIdeal functional (log-space; the
|
||
// functional requires b_i > 0 / a_e > 0). We therefore choose a valid base
|
||
// point, evaluate G there, and absorb G into the targets so that G(xbase) = 0.
|
||
// Newton tests then start from a perturbation of xbase.
|
||
//
|
||
// Returns xbase so callers can construct a perturbed starting point.
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
static std::vector<double> set_natural_hyper_ideal_targets(
|
||
ConformalMesh& mesh, HyperIdealMaps& maps, int n,
|
||
double b_base = 1.0, double a_base = 0.5)
|
||
{
|
||
const auto sz = static_cast<std::size_t>(n);
|
||
std::vector<double> xbase(sz, 0.0);
|
||
|
||
for (auto v : mesh.vertices()) {
|
||
int iv = maps.v_idx[v];
|
||
if (iv >= 0) xbase[static_cast<std::size_t>(iv)] = b_base;
|
||
}
|
||
for (auto e : mesh.edges()) {
|
||
int ie = maps.e_idx[e];
|
||
if (ie >= 0) xbase[static_cast<std::size_t>(ie)] = a_base;
|
||
}
|
||
|
||
// G = Σβ - theta_target (initial target = 0 → G = Σβ = "actual" angles)
|
||
auto G = evaluate_hyper_ideal(mesh, xbase, maps, /*energy=*/false).gradient;
|
||
|
||
// Set target := actual so that G(xbase) = actual - target = 0
|
||
for (auto v : mesh.vertices()) {
|
||
int iv = maps.v_idx[v];
|
||
if (iv < 0) continue;
|
||
maps.theta_v[v] += G[static_cast<std::size_t>(iv)];
|
||
}
|
||
for (auto e : mesh.edges()) {
|
||
int ie = maps.e_idx[e];
|
||
if (ie < 0) continue;
|
||
maps.theta_e[e] += G[static_cast<std::size_t>(ie)];
|
||
}
|
||
|
||
return xbase;
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// HyperIdeal 1 — Triangle, all DOFs variable, converges from perturbation
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(NewtonSolver, HyperIdeal_ConvergesTriangleAllVariable)
|
||
{
|
||
auto mesh = make_triangle();
|
||
auto maps = setup_hyper_ideal_maps(mesh);
|
||
int n = assign_hyper_ideal_all_dof_indices(mesh, maps);
|
||
|
||
// xbase = (b=1.0, a=0.5) is the equilibrium after natural-target setup.
|
||
auto xbase = set_natural_hyper_ideal_targets(mesh, maps, n);
|
||
|
||
// Perturb by +0.2 uniformly
|
||
std::vector<double> x0 = xbase;
|
||
for (auto& v : x0) v += 0.2;
|
||
|
||
auto res = newton_hyper_ideal(mesh, x0, maps, /*tol=*/1e-7, /*max_iter=*/100);
|
||
|
||
EXPECT_TRUE(res.converged)
|
||
<< "Newton (HyperIdeal, triangle) should converge; "
|
||
"grad_inf_norm = " << res.grad_inf_norm;
|
||
EXPECT_LT(res.grad_inf_norm, 1e-7);
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// HyperIdeal 2 — Result fields self-consistent
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(NewtonSolver, HyperIdeal_ResultFieldsConsistent)
|
||
{
|
||
auto mesh = make_triangle();
|
||
auto maps = setup_hyper_ideal_maps(mesh);
|
||
int n = assign_hyper_ideal_all_dof_indices(mesh, maps);
|
||
|
||
auto xbase = set_natural_hyper_ideal_targets(mesh, maps, n);
|
||
|
||
std::vector<double> x0 = xbase;
|
||
for (auto& v : x0) v += 0.1;
|
||
|
||
auto res = newton_hyper_ideal(mesh, x0, maps, /*tol=*/1e-7, /*max_iter=*/100);
|
||
|
||
EXPECT_EQ(static_cast<int>(res.x.size()), n);
|
||
|
||
// Reported grad_inf_norm must match re-computed gradient at res.x
|
||
auto G = evaluate_hyper_ideal(mesh, res.x, maps, false).gradient;
|
||
double actual_inf = 0.0;
|
||
for (double v : G) actual_inf = std::max(actual_inf, std::abs(v));
|
||
EXPECT_NEAR(actual_inf, res.grad_inf_norm, 1e-9);
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// HyperIdeal 3 — Tetrahedron (10 DOFs): 4 vertex b-vals + 6 edge a-vals
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(NewtonSolver, HyperIdeal_ConvergesTetrahedron)
|
||
{
|
||
auto mesh = make_tetrahedron();
|
||
auto maps = setup_hyper_ideal_maps(mesh);
|
||
int n = assign_hyper_ideal_all_dof_indices(mesh, maps);
|
||
|
||
auto xbase = set_natural_hyper_ideal_targets(mesh, maps, n);
|
||
|
||
// Perturb by +0.15
|
||
std::vector<double> x0 = xbase;
|
||
for (auto& v : x0) v += 0.15;
|
||
|
||
auto res = newton_hyper_ideal(mesh, x0, maps, /*tol=*/1e-7, /*max_iter=*/200);
|
||
|
||
EXPECT_TRUE(res.converged)
|
||
<< "Newton (HyperIdeal, tetrahedron) should converge; "
|
||
"grad_inf_norm = " << res.grad_inf_norm;
|
||
EXPECT_LT(res.grad_inf_norm, 1e-7);
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// HyperIdeal 4 — SparseQR fallback: solver returns a result (no crash)
|
||
//
|
||
// With all targets = 0 the equilibrium is not at x=0 but the solver should
|
||
// at minimum not crash and return a consistent result struct.
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(NewtonSolver, HyperIdeal_SparseQRFallbackNoCrash)
|
||
{
|
||
auto mesh = make_triangle();
|
||
auto maps = setup_hyper_ideal_maps(mesh);
|
||
int n = assign_hyper_ideal_all_dof_indices(mesh, maps);
|
||
|
||
// Leave targets at their default (0): solver tries to solve but the
|
||
// "equilibrium" is at some unknown x*. With valid starting point the
|
||
// Hessian is positive-definite and the solver should not crash.
|
||
// We don't assert convergence — just that the result struct is consistent.
|
||
std::vector<double> x0(static_cast<std::size_t>(n), 1.0);
|
||
// Mix vertex / edge DOFs: b=1.0, a=0.5 (valid region of the functional)
|
||
for (auto e : mesh.edges()) {
|
||
int ie = maps.e_idx[e];
|
||
if (ie >= 0) x0[static_cast<std::size_t>(ie)] = 0.5;
|
||
}
|
||
auto res = newton_hyper_ideal(mesh, x0, maps, /*tol=*/1e-7, /*max_iter=*/50);
|
||
|
||
// Struct fields must always be populated
|
||
EXPECT_EQ(static_cast<int>(res.x.size()), n);
|
||
EXPECT_GE(res.iterations, 0);
|
||
EXPECT_FALSE(std::isnan(res.grad_inf_norm));
|
||
EXPECT_FALSE(std::isinf(res.grad_inf_norm));
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// SparseQR fallback — Test 12: solve_linear_system recovers correct solution
|
||
//
|
||
// The public API solve_linear_system(A, rhs) must return the correct answer
|
||
// for a well-conditioned full-rank system (LDLT path taken).
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SparseQRFallback, FullRankSystem_CorrectSolution)
|
||
{
|
||
// Build a simple 3×3 diagonal PD matrix: A = diag(1, 2, 3)
|
||
Eigen::SparseMatrix<double> A(3, 3);
|
||
A.insert(0, 0) = 1.0;
|
||
A.insert(1, 1) = 2.0;
|
||
A.insert(2, 2) = 3.0;
|
||
A.makeCompressed();
|
||
|
||
Eigen::VectorXd rhs(3);
|
||
rhs << 1.0, 4.0, 9.0; // solution = [1, 2, 3]
|
||
|
||
bool fallback = true; // expect it to be set to false (LDLT succeeds)
|
||
Eigen::VectorXd x = conformallab::solve_linear_system(A, rhs, &fallback);
|
||
|
||
EXPECT_FALSE(fallback) << "Full-rank system: LDLT should succeed (no SparseQR needed)";
|
||
EXPECT_NEAR(x[0], 1.0, 1e-12);
|
||
EXPECT_NEAR(x[1], 2.0, 1e-12);
|
||
EXPECT_NEAR(x[2], 3.0, 1e-12);
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// SparseQR fallback — Test 13: fallback_used=true on a singular matrix
|
||
//
|
||
// Construct a symmetric 3×3 matrix of rank 1 where LDLT fails (the (2,2)
|
||
// pivot is zero). SparseQR finds the minimum-norm least-squares solution.
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SparseQRFallback, SingularMatrix_FallbackActivated)
|
||
{
|
||
// A = [[2, 0, 0],
|
||
// [0, 0, 0], ← zero pivot → LDLT failure
|
||
// [0, 0, 3]]
|
||
// rhs compatible with the row space: [2, 0, 3] → solution [1, 0, 1]
|
||
Eigen::SparseMatrix<double> A(3, 3);
|
||
A.insert(0, 0) = 2.0;
|
||
// row/col 1 deliberately all-zero
|
||
A.insert(2, 2) = 3.0;
|
||
A.makeCompressed();
|
||
|
||
Eigen::VectorXd rhs(3);
|
||
rhs << 2.0, 0.0, 3.0;
|
||
|
||
bool fallback = false;
|
||
Eigen::VectorXd x = conformallab::solve_linear_system(A, rhs, &fallback);
|
||
|
||
EXPECT_TRUE(fallback) << "Singular matrix: SparseQR fallback must be triggered";
|
||
// SparseQR min-norm solution: x[0]=1, x[1]=0, x[2]=1
|
||
EXPECT_NEAR(x[0], 1.0, 1e-10);
|
||
EXPECT_NEAR(x[1], 0.0, 1e-10);
|
||
EXPECT_NEAR(x[2], 1.0, 1e-10);
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// SparseQR fallback — Test 14: Euclidean Newton on a closed mesh, no pinning
|
||
//
|
||
// make_tetrahedron() is a closed surface (4 vertices, 4 faces). Without a
|
||
// pinned vertex the Euclidean Hessian has a 1-D null space (uniform scale
|
||
// gauge mode): H·1 = 0. SimplicialLDLT fails on this rank-deficient H;
|
||
// SparseQR finds the min-norm Newton step orthogonal to the null space.
|
||
//
|
||
// The gradient always lives in the row space of H (Σ G_v = 0 by angle-sum
|
||
// invariance), so the SparseQR step is also the Newton step and the solver
|
||
// converges to the natural equilibrium.
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SparseQRFallback, Euclidean_ClosedMeshNoPinConverges)
|
||
{
|
||
auto mesh = make_tetrahedron();
|
||
auto maps = setup_euclidean_maps(mesh);
|
||
compute_euclidean_lambda0_from_mesh(mesh, maps);
|
||
|
||
// Assign all 4 vertices as free DOFs (no pinning).
|
||
int idx = 0;
|
||
for (auto v : mesh.vertices())
|
||
maps.v_idx[v] = idx++;
|
||
const int n = idx; // = 4
|
||
|
||
// Natural theta: equilibrium at x* = 0.
|
||
set_natural_euclidean_theta(mesh, maps, n);
|
||
|
||
std::vector<double> x0(static_cast<std::size_t>(n), -0.1);
|
||
auto res = newton_euclidean(mesh, x0, maps, /*tol=*/1e-8, /*max_iter=*/100);
|
||
|
||
EXPECT_TRUE(res.converged)
|
||
<< "Euclidean Newton on closed tetrahedron (no pin) must converge via SparseQR; "
|
||
"grad_inf_norm = " << res.grad_inf_norm;
|
||
EXPECT_LT(res.grad_inf_norm, 1e-8);
|
||
EXPECT_EQ(res.status, NewtonStatus::Converged);
|
||
// N7: the closed mesh with no pinned vertex has a gauge-singular Hessian.
|
||
// That near-singularity must be observable in the diagnostics — either the
|
||
// SparseQR fallback fired, or (when LDLT "succeeds" with a ~0 pivot, which is
|
||
// exactly the silent case N7 targets) the smallest LDLT pivot is tiny.
|
||
EXPECT_TRUE(res.sparse_qr_fallback_used || res.min_ldlt_pivot < 1e-6)
|
||
<< "gauge singularity should surface via fallback flag or min_ldlt_pivot; "
|
||
<< "fallback=" << res.sparse_qr_fallback_used
|
||
<< " min_pivot=" << res.min_ldlt_pivot;
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// C1: solve_linear_system exposes the ok flag via the new out-parameter
|
||
//
|
||
// The C1 audit finding was that solve_linear_system silently dropped the
|
||
// internal ok flag, making double-solver failure invisible to callers.
|
||
// These tests verify that the new optional bool* ok parameter is correctly
|
||
// wired: ok=true for successful solves, and the parameter is optional (null
|
||
// is safe).
|
||
//
|
||
// Note: constructing a matrix where Eigen's SparseQR hard-fails is not
|
||
// practical — SparseQR reports Eigen::Success even for rank-0 inputs,
|
||
// returning a zero min-norm solution. The observable contract for callers is
|
||
// therefore: ok=true whenever a solver reports success; ok=false only when
|
||
// both report failure (an edge case Eigen essentially never produces). The
|
||
// fix in C1 ensures that if that edge case ever fires it propagates to the
|
||
// caller — previously it was silently swallowed.
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(SparseQRFallback, OkFlag_TrueOnFullRankSystem)
|
||
{
|
||
// 3×3 diagonal PD matrix: LDLT succeeds → ok must be true.
|
||
Eigen::SparseMatrix<double> A(3, 3);
|
||
A.insert(0, 0) = 1.0;
|
||
A.insert(1, 1) = 2.0;
|
||
A.insert(2, 2) = 3.0;
|
||
A.makeCompressed();
|
||
|
||
Eigen::VectorXd rhs(3);
|
||
rhs << 1.0, 4.0, 9.0;
|
||
|
||
bool fallback = true;
|
||
bool ok = false;
|
||
Eigen::VectorXd x = conformallab::solve_linear_system(A, rhs, &fallback, &ok);
|
||
|
||
EXPECT_TRUE(ok) << "Full-rank system: ok must be true";
|
||
EXPECT_FALSE(fallback) << "Full-rank system: LDLT path, no fallback expected";
|
||
EXPECT_NEAR(x[0], 1.0, 1e-12);
|
||
EXPECT_NEAR(x[1], 2.0, 1e-12);
|
||
EXPECT_NEAR(x[2], 3.0, 1e-12);
|
||
}
|
||
|
||
TEST(SparseQRFallback, OkFlag_TrueOnSingularSystemViaFallback)
|
||
{
|
||
// Singular matrix (zero row/col 1) — LDLT fails, SparseQR succeeds → ok=true.
|
||
Eigen::SparseMatrix<double> A(3, 3);
|
||
A.insert(0, 0) = 2.0;
|
||
A.insert(2, 2) = 3.0;
|
||
A.makeCompressed();
|
||
|
||
Eigen::VectorXd rhs(3);
|
||
rhs << 2.0, 0.0, 3.0;
|
||
|
||
bool fallback = false;
|
||
bool ok = false;
|
||
Eigen::VectorXd x = conformallab::solve_linear_system(A, rhs, &fallback, &ok);
|
||
|
||
EXPECT_TRUE(ok) << "Singular-but-consistent system: SparseQR succeeds → ok must be true";
|
||
EXPECT_TRUE(fallback) << "Singular system: fallback must have been triggered";
|
||
}
|
||
|
||
TEST(SparseQRFallback, OkFlag_NullPointerIsSafe)
|
||
{
|
||
// Calling with ok=nullptr must not crash (backward-compatibility).
|
||
Eigen::SparseMatrix<double> A(2, 2);
|
||
A.insert(0, 0) = 1.0;
|
||
A.insert(1, 1) = 1.0;
|
||
A.makeCompressed();
|
||
|
||
Eigen::VectorXd rhs(2);
|
||
rhs << 3.0, 5.0;
|
||
|
||
EXPECT_NO_THROW({
|
||
Eigen::VectorXd x = conformallab::solve_linear_system(A, rhs, nullptr, nullptr);
|
||
EXPECT_NEAR(x[0], 3.0, 1e-12);
|
||
EXPECT_NEAR(x[1], 5.0, 1e-12);
|
||
});
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// H2 / B2 — newton_core reuses the line-search gradient
|
||
//
|
||
// All five solvers now share detail::newton_core. B2: the gradient computed at
|
||
// the accepted line-search point is carried into the next iteration's
|
||
// convergence check instead of being recomputed at the loop top. On a
|
||
// unit-Hessian quadratic, exact Newton reaches the minimum in one step, so the
|
||
// total gradient-eval count is exactly 2 (1 initial + 1 accepted trial); the
|
||
// pre-B2 loop would have recomputed G at the top of iteration 1 → 3.
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(NewtonCore, ReusesLineSearchGradient_B2_Convex)
|
||
{
|
||
const std::vector<double> xstar = {1.0, -2.0, 3.5};
|
||
int grad_calls = 0;
|
||
|
||
auto grad = [&](const std::vector<double>& x) {
|
||
++grad_calls;
|
||
std::vector<double> g(x.size());
|
||
for (std::size_t i = 0; i < x.size(); ++i) g[i] = x[i] - xstar[i];
|
||
return g; // G(x) = x − x*, H = +I
|
||
};
|
||
auto hess = [&](const std::vector<double>&) {
|
||
Eigen::SparseMatrix<double> H(3, 3);
|
||
for (int i = 0; i < 3; ++i) H.insert(i, i) = 1.0;
|
||
H.makeCompressed();
|
||
return H;
|
||
};
|
||
|
||
auto res = conformallab::detail::newton_core(
|
||
std::vector<double>{0.0, 0.0, 0.0}, grad, hess,
|
||
/*concave=*/false, /*tol=*/1e-9, /*max_iter=*/50);
|
||
|
||
ASSERT_TRUE(res.converged);
|
||
EXPECT_EQ(res.iterations, 1);
|
||
EXPECT_NEAR(res.x[0], 1.0, 1e-9);
|
||
EXPECT_NEAR(res.x[1], -2.0, 1e-9);
|
||
EXPECT_NEAR(res.x[2], 3.5, 1e-9);
|
||
EXPECT_EQ(grad_calls, 2)
|
||
<< "B2: the loop-top gradient must be reused from the line search";
|
||
}
|
||
|
||
// Concave path (spherical-style): H is NSD, newton_core factors −H and solves
|
||
// (−H)·Δx = G. G(x) = x* − x with H = −I makes x* a maximiser; the core must
|
||
// still converge in one step.
|
||
TEST(NewtonCore, ConcavePathConverges)
|
||
{
|
||
const std::vector<double> xstar = {0.5, -1.5, 2.0};
|
||
int grad_calls = 0;
|
||
|
||
auto grad = [&](const std::vector<double>& x) {
|
||
++grad_calls;
|
||
std::vector<double> g(x.size());
|
||
for (std::size_t i = 0; i < x.size(); ++i) g[i] = xstar[i] - x[i];
|
||
return g; // G(x) = x* − x, H = −I
|
||
};
|
||
auto hess = [&](const std::vector<double>&) {
|
||
Eigen::SparseMatrix<double> H(3, 3);
|
||
for (int i = 0; i < 3; ++i) H.insert(i, i) = -1.0;
|
||
H.makeCompressed();
|
||
return H;
|
||
};
|
||
|
||
auto res = conformallab::detail::newton_core(
|
||
std::vector<double>{0.0, 0.0, 0.0}, grad, hess,
|
||
/*concave=*/true, /*tol=*/1e-9, /*max_iter=*/50);
|
||
|
||
ASSERT_TRUE(res.converged);
|
||
EXPECT_EQ(res.iterations, 1);
|
||
EXPECT_NEAR(res.x[0], 0.5, 1e-9);
|
||
EXPECT_NEAR(res.x[1], -1.5, 1e-9);
|
||
EXPECT_NEAR(res.x[2], 2.0, 1e-9);
|
||
EXPECT_EQ(grad_calls, 2);
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// I1 / H1 / N7 — newton_core termination status, iteration count, diagnostics
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
namespace {
|
||
// 1×1 sparse identity / scalar helpers for the synthetic newton_core tests.
|
||
inline Eigen::SparseMatrix<double> diag_sparse(std::initializer_list<double> d)
|
||
{
|
||
const int n = static_cast<int>(d.size());
|
||
Eigen::SparseMatrix<double> H(n, n);
|
||
int i = 0;
|
||
for (double v : d) { H.insert(i, i) = v; ++i; }
|
||
H.makeCompressed();
|
||
return H;
|
||
}
|
||
} // namespace
|
||
|
||
// I1: a converging solve reports status == Converged and the N7 pivot is the
|
||
// (well-conditioned) unit diagonal; no SparseQR fallback.
|
||
TEST(NewtonCore, Status_Converged_AndDiagnostics_N7)
|
||
{
|
||
const std::vector<double> xstar = {2.0, -1.0};
|
||
auto grad = [&](const std::vector<double>& x) {
|
||
return std::vector<double>{ x[0] - xstar[0], x[1] - xstar[1] };
|
||
};
|
||
auto hess = [&](const std::vector<double>&) { return diag_sparse({1.0, 1.0}); };
|
||
|
||
auto res = conformallab::detail::newton_core(
|
||
std::vector<double>{0.0, 0.0}, grad, hess, /*concave=*/false, 1e-9, 50);
|
||
|
||
EXPECT_TRUE(res.converged);
|
||
EXPECT_EQ(res.status, conformallab::NewtonStatus::Converged);
|
||
EXPECT_EQ(res.iterations, 1);
|
||
EXPECT_FALSE(res.sparse_qr_fallback_used); // N7
|
||
EXPECT_NEAR(res.min_ldlt_pivot, 1.0, 1e-12); // N7: D = I
|
||
}
|
||
|
||
// I1: a slow (linearly-convergent cubic) problem that does not reach tol within
|
||
// max_iter reports MaxIterations, with iterations == max_iter (H1).
|
||
TEST(NewtonCore, Status_MaxIterations)
|
||
{
|
||
// G(x) = x³, H = 3x² → Newton map x ← (2/3)x (linear convergence).
|
||
auto grad = [&](const std::vector<double>& x) {
|
||
return std::vector<double>{ x[0]*x[0]*x[0] };
|
||
};
|
||
auto hess = [&](const std::vector<double>& x) {
|
||
return diag_sparse({ 3.0 * x[0] * x[0] });
|
||
};
|
||
|
||
auto res = conformallab::detail::newton_core(
|
||
std::vector<double>{1.0}, grad, hess, /*concave=*/false, /*tol=*/1e-8, /*max_iter=*/3);
|
||
|
||
EXPECT_FALSE(res.converged);
|
||
EXPECT_EQ(res.status, conformallab::NewtonStatus::MaxIterations);
|
||
EXPECT_EQ(res.iterations, 3);
|
||
}
|
||
|
||
// I1/H1: a constant non-zero gradient cannot be reduced by any step, so the
|
||
// line search stalls on the first iteration → LineSearchStalled, iterations == 0.
|
||
TEST(NewtonCore, Status_LineSearchStalled)
|
||
{
|
||
auto grad = [&](const std::vector<double>&) {
|
||
return std::vector<double>{ 1.0, 1.0 }; // constant, independent of x
|
||
};
|
||
auto hess = [&](const std::vector<double>&) { return diag_sparse({1.0, 1.0}); };
|
||
|
||
auto res = conformallab::detail::newton_core(
|
||
std::vector<double>{0.0, 0.0}, grad, hess, /*concave=*/false, 1e-9, 50);
|
||
|
||
EXPECT_FALSE(res.converged);
|
||
EXPECT_EQ(res.status, conformallab::NewtonStatus::LineSearchStalled);
|
||
EXPECT_EQ(res.iterations, 0); // H1: no step completed
|
||
}
|
||
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
// H5 (test-coverage audit, 2026-06-01): degenerate-triangle integration test
|
||
//
|
||
// Finding H5: euclidean_hessian.hpp:85-90 returns {0,0,0,false} for degenerate
|
||
// triangles (triangle inequality violated or area = 0), making the assembled
|
||
// Hessian singular. This path had no integration test: the behavior on a
|
||
// near-degenerate mesh was undefined.
|
||
//
|
||
// Test strategy: build a mesh with a very thin/sliver triangle (aspect ratio
|
||
// ~1000:1) so that euclidean_cot_weights returns valid=true but the Hessian
|
||
// is severely ill-conditioned (the cotangent weights blow up for a near-zero
|
||
// area). Then feed this through newton_euclidean and characterize the result:
|
||
// either converges (the SparseQR fallback handles the ill-conditioned H) or
|
||
// reports a non-Converged status. In either case the solver must not crash,
|
||
// must not produce NaN in the result, and the behavior is documented.
|
||
//
|
||
// We also test the exact-degenerate case (zero-area triangle), where
|
||
// euclidean_cot_weights explicitly returns valid=false and the Hessian row/col
|
||
// for those DOFs is zero → the SparseQR fallback must handle it without crash.
|
||
// ════════════════════════════════════════════════════════════════════════════
|
||
|
||
TEST(NewtonSolver, Euclidean_SliverTriangle_CharacterizedBehavior)
|
||
{
|
||
// Build a very thin sliver triangle: v0=(0,0), v1=(1,0), v2=(0,1e-4).
|
||
// Area ≈ 5e-5, aspect ratio ≈ 10000. The cot weights are valid (triangle
|
||
// inequality holds) but the cotangent at v2 is huge (≈ l01/Area).
|
||
ConformalMesh mesh;
|
||
auto v0 = mesh.add_vertex(Point3(0.0, 0.0, 0.0));
|
||
auto v1 = mesh.add_vertex(Point3(1.0, 0.0, 0.0));
|
||
auto v2 = mesh.add_vertex(Point3(0.0, 1e-4, 0.0));
|
||
mesh.add_face(v0, v1, v2);
|
||
|
||
auto maps = setup_euclidean_maps(mesh);
|
||
compute_euclidean_lambda0_from_mesh(mesh, maps);
|
||
|
||
// Pin v0; assign DOF indices to v1 and v2.
|
||
maps.v_idx[v0] = -1;
|
||
maps.v_idx[v1] = 0;
|
||
maps.v_idx[v2] = 1;
|
||
const int n = 2;
|
||
|
||
// Natural theta: equilibrium at x* = 0 by construction.
|
||
set_natural_euclidean_theta(mesh, maps, n);
|
||
|
||
std::vector<double> x0(n, 0.0);
|
||
auto res = newton_euclidean(mesh, x0, maps, /*tol=*/1e-8, /*max_iter=*/100);
|
||
|
||
// H5 acceptance criterion: behavior is characterized, not undefined.
|
||
// The solver must not crash or produce NaN.
|
||
EXPECT_EQ(static_cast<int>(res.x.size()), n)
|
||
<< "Result vector must always be populated";
|
||
for (double xi : res.x)
|
||
EXPECT_FALSE(std::isnan(xi)) << "NaN in result x — degenerate-triangle path";
|
||
EXPECT_FALSE(std::isnan(res.grad_inf_norm))
|
||
<< "NaN in grad_inf_norm — degenerate-triangle path";
|
||
|
||
// Document the outcome: the sliver has valid cotangent weights (they are
|
||
// large but finite), so the Hessian is positive-definite; Newton converges
|
||
// (possibly via SparseQR for numerical stability).
|
||
// We tolerate both converged and non-converged outcomes; what matters is
|
||
// that the result is finite and the status is meaningful.
|
||
EXPECT_NE(res.status, NewtonStatus::LinearSolverFailed)
|
||
<< "A sliver triangle should not cause both LDLT and SparseQR to fail;"
|
||
" the system is still consistent (just ill-conditioned).";
|
||
}
|
||
|
||
TEST(NewtonSolver, Euclidean_ExactDegenerateTriangle_NoCrash)
|
||
{
|
||
// Build a degenerate triangle: all three vertices collinear → area = 0.
|
||
// v0=(0,0), v1=(1,0), v2=(2,0). This forces kahan <= 0 in
|
||
// euclidean_cot_weights → {0,0,0,false}. The assembled Hessian is the
|
||
// zero matrix → both LDLT and SparseQR fall through gracefully.
|
||
ConformalMesh mesh;
|
||
auto v0 = mesh.add_vertex(Point3(0.0, 0.0, 0.0));
|
||
auto v1 = mesh.add_vertex(Point3(1.0, 0.0, 0.0));
|
||
auto v2 = mesh.add_vertex(Point3(2.0, 0.0, 0.0));
|
||
mesh.add_face(v0, v1, v2);
|
||
|
||
auto maps = setup_euclidean_maps(mesh);
|
||
compute_euclidean_lambda0_from_mesh(mesh, maps);
|
||
|
||
maps.v_idx[v0] = -1;
|
||
maps.v_idx[v1] = 0;
|
||
maps.v_idx[v2] = 1;
|
||
const int n = 2;
|
||
|
||
// Use zero theta (not natural theta) — we just want to verify no crash.
|
||
std::vector<double> x0(n, 0.0);
|
||
|
||
// H5 acceptance criterion: no crash, no UB, result struct populated.
|
||
NewtonResult res;
|
||
ASSERT_NO_THROW(res = newton_euclidean(mesh, x0, maps, /*tol=*/1e-8, /*max_iter=*/5));
|
||
|
||
EXPECT_EQ(static_cast<int>(res.x.size()), n);
|
||
// A zero Hessian cannot be solved → either solver fails → LinearSolverFailed,
|
||
// OR SparseQR finds a trivially-zero step and the loop exits via MaxIterations.
|
||
// Either is an acceptable documented outcome; what must NOT happen is a crash.
|
||
EXPECT_TRUE(res.status == NewtonStatus::LinearSolverFailed
|
||
|| res.status == NewtonStatus::MaxIterations
|
||
|| res.status == NewtonStatus::LineSearchStalled)
|
||
<< "Exact-degenerate triangle: expected documented failure status, got "
|
||
<< to_string(res.status);
|
||
}
|