fix(solver+io): B1 block-FD Hessian, V3 NaN/Inf guard, C1 ok-flag (audit quick-wins)

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 <noreply@anthropic.com>

code/include/mesh_io.hpp

@@ -29,6 +29,7 @@
 #include <CGAL/boost/graph/helpers.h>
 #include <string>
 #include <stdexcept>
+#include <cmath>
 
 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;
 }
 

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<double>& 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;

code/tests/cgal/test_mesh_io.cpp

@@ -18,6 +18,7 @@
 #include <cmath>
 #include <cstdio>
 #include <filesystem>
+#include <fstream>
 #include <stdexcept>
 
 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);
+}

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<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);
+    });
+}