From 101f3138ce76ccc62059792f34e9e7a238e882a6 Mon Sep 17 00:00:00 2001 From: Tarik Moussa Date: Sun, 31 May 2026 14:03:45 +0200 Subject: [PATCH] fix(solver+io): B1 block-FD Hessian, V3 NaN/Inf guard, C1 ok-flag (audit quick-wins) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit B1 (api-performance): wire hyper_ideal_hessian_block_fd_sym into newton_hyper_ideal instead of the full-FD path — 33×/1166× faster on cathead/brezel, also fixes the FD-step vs Newton-tol accuracy floor (N1 cross-fix). V3 (input-validation): add isfinite check on all vertex coordinates in load_mesh; NaN/Inf now throws runtime_error at the I/O boundary before poisoning the solver. C1 (test-coverage): expose the internal ok flag via an optional bool* parameter on solve_linear_system so double-solver failure is no longer invisible to callers. +5 new tests (LoadMeshThrowsOnNaN/Inf, OkFlag_True*, OkFlag_NullPointerIsSafe). 282/282 tests pass. Co-Authored-By: Claude Sonnet 4.6 --- code/include/mesh_io.hpp | 7 +++ code/include/newton_solver.hpp | 14 +++-- code/tests/cgal/test_mesh_io.cpp | 40 +++++++++++++ code/tests/cgal/test_newton_solver.cpp | 78 ++++++++++++++++++++++++++ 4 files changed, 134 insertions(+), 5 deletions(-) diff --git a/code/include/mesh_io.hpp b/code/include/mesh_io.hpp index 48bf145..84a4e40 100644 --- a/code/include/mesh_io.hpp +++ b/code/include/mesh_io.hpp @@ -29,6 +29,7 @@ #include #include #include +#include namespace conformallab { @@ -62,6 +63,12 @@ inline ConformalMesh load_mesh(const std::string& filename) throw std::runtime_error( "conformallab: mesh from " + filename + " is not triangulated (functionals require triangle faces)"); + for (auto v : mesh.vertices()) { + const auto& p = mesh.point(v); + if (!std::isfinite(p.x()) || !std::isfinite(p.y()) || !std::isfinite(p.z())) + throw std::runtime_error( + "conformallab: non-finite vertex coordinate in " + filename); + } return mesh; } diff --git a/code/include/newton_solver.hpp b/code/include/newton_solver.hpp index 25944bc..193e801 100644 --- a/code/include/newton_solver.hpp +++ b/code/include/newton_solver.hpp @@ -103,14 +103,18 @@ inline Eigen::VectorXd solve_with_fallback( /// 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& A, const Eigen::VectorXd& rhs, - bool* fallback_used = nullptr) + bool* fallback_used = nullptr, + bool* ok = nullptr) { - bool ok = false; - return detail::solve_with_fallback(A, rhs, ok, fallback_used); + 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 @@ -429,8 +433,8 @@ inline NewtonResult newton_hyper_ideal( return res; } - // ── Hessian (numerical FD) + solve H·Δx = −G ───────────────────────── - auto H = hyper_ideal_hessian_sym(mesh, x, m, hess_eps); + // ── 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; diff --git a/code/tests/cgal/test_mesh_io.cpp b/code/tests/cgal/test_mesh_io.cpp index f767211..0ee1447 100644 --- a/code/tests/cgal/test_mesh_io.cpp +++ b/code/tests/cgal/test_mesh_io.cpp @@ -18,6 +18,7 @@ #include #include #include +#include #include using namespace conformallab; @@ -171,3 +172,42 @@ TEST(MeshIO, SaveLoadConvenienceWrappers) rm(path); } + +// ════════════════════════════════════════════════════════════════════════════ +// V3: load_mesh throws on non-finite vertex coordinate (NaN / Inf) +// +// A mesh file containing NaN or Inf vertex coordinates must be rejected at +// the I/O boundary — before the coordinates can silently poison the solver. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(MeshIO, LoadMeshThrowsOnNaNCoordinate) +{ + // Write a minimal OFF triangle with a NaN in the first vertex. + auto path = tmp_path("nan.off"); + { + std::ofstream ofs(path); + ofs << "OFF\n3 1 0\n" + << "nan 0.0 0.0\n" + << "1.0 0.0 0.0\n" + << "0.0 1.0 0.0\n" + << "3 0 1 2\n"; + } + EXPECT_THROW(load_mesh(path), std::runtime_error); + rm(path); +} + +TEST(MeshIO, LoadMeshThrowsOnInfCoordinate) +{ + // Write a minimal OFF triangle with Inf in the y-coordinate. + auto path = tmp_path("inf.off"); + { + std::ofstream ofs(path); + ofs << "OFF\n3 1 0\n" + << "0.0 0.0 0.0\n" + << "1.0 inf 0.0\n" + << "0.0 1.0 0.0\n" + << "3 0 1 2\n"; + } + EXPECT_THROW(load_mesh(path), std::runtime_error); + rm(path); +} diff --git a/code/tests/cgal/test_newton_solver.cpp b/code/tests/cgal/test_newton_solver.cpp index d9ce3c3..a63789c 100644 --- a/code/tests/cgal/test_newton_solver.cpp +++ b/code/tests/cgal/test_newton_solver.cpp @@ -501,3 +501,81 @@ TEST(SparseQRFallback, Euclidean_ClosedMeshNoPinConverges) "grad_inf_norm = " << res.grad_inf_norm; EXPECT_LT(res.grad_inf_norm, 1e-8); } + +// ════════════════════════════════════════════════════════════════════════════ +// 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 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 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 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); + }); +}