Unvalid lambdas when matrix.colPivHouseholderQr().solve(z) is not a solution

2 files changed, 23 insertions, 6 deletions

diff --git a/src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h b/src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h
index 25f688de..dca97331 100644
--- a/src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h
+++ b/src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h
@@ -17,8 +17,11 @@
 #include <gudhi/Functions/Constant_function.h>
 #include <gudhi/Functions/PL_approximation.h>
 #include <gudhi/Query_result.h>
+#include <gudhi/Debug_utils.h>  // for GUDHI_CHECK
 
 #include <vector>
+#include <limits>  // for std::numeric_limits<>
+#include <cmath>   // for std::fabs
 
 namespace Gudhi {
 
@@ -53,6 +56,11 @@ class Implicit_manifold_intersection_oracle {
     z(0) = 1;
     for (std::size_t i = 1; i < cod_d + 1; ++i) z(i) = 0;
     Eigen::VectorXd lambda = matrix.colPivHouseholderQr().solve(z);
+    if (!z.isApprox(matrix*lambda)) {
+      // NaN non valid results
+      for (std::size_t i = 0; i < (std::size_t)lambda.size(); ++i) lambda(i) =
+        std::numeric_limits<double>::quiet_NaN();
+    }
     return lambda;
   }
 
@@ -75,6 +83,11 @@ class Implicit_manifold_intersection_oracle {
     z(0) = 1;
     for (std::size_t i = 1; i < cod_d + 2; ++i) z(i) = 0;
     Eigen::VectorXd lambda = matrix.colPivHouseholderQr().solve(z);
+    if (!z.isApprox(matrix*lambda)) {
+      // NaN non valid results
+      for (std::size_t i = 0; i < (std::size_t)lambda.size(); ++i) lambda(i) =
+        std::numeric_limits<double>::quiet_NaN();
+    }
     return lambda;
   }
 
@@ -85,10 +98,13 @@ class Implicit_manifold_intersection_oracle {
     using QR = Query_result<Simplex_handle>;
     std::size_t amb_d = triangulation.dimension();
     std::size_t cod_d = simplex.dimension();
-
-    for (std::size_t i = 0; i < (std::size_t)lambda.size(); ++i)
-      if (lambda(i) < 0 || lambda(i) > 1) return QR({Eigen::VectorXd(), false});
-
+    for (std::size_t i = 0; i < (std::size_t)lambda.size(); ++i) {
+      if (std::isnan(lambda(i))) return QR({Eigen::VectorXd(), false});
+      GUDHI_CHECK((std::fabs(lambda(i) - 1.) > std::numeric_limits<double>::epsilon() &&
+                   std::fabs(lambda(i) - 0.) > std::numeric_limits<double>::epsilon()),
+                  std::invalid_argument("A vertex of the triangulation lies exactly on the manifold"));
+      if (lambda(i) < 0. || lambda(i) > 1.) return QR({Eigen::VectorXd(), false});
+    }
     Eigen::MatrixXd vertex_matrix(cod_d + 1, amb_d);
     auto v_range = simplex.vertex_range();
     auto v_it = v_range.begin();
@@ -159,6 +175,7 @@ class Implicit_manifold_intersection_oracle {
   template <class Simplex_handle, class Triangulation>
   Query_result<Simplex_handle> intersects_boundary(const Simplex_handle& simplex,
                                                    const Triangulation& triangulation) const {
+    //std::cout << "intersects_boundary" << std::endl;
     Eigen::VectorXd lambda = compute_boundary_lambda(simplex, triangulation);
     return intersection_result(lambda, simplex, triangulation);
   }
diff --git a/src/Coxeter_triangulation/test/manifold_tracing_test.cpp b/src/Coxeter_triangulation/test/manifold_tracing_test.cpp
index 68b79487..2d7d6087 100644
--- a/src/Coxeter_triangulation/test/manifold_tracing_test.cpp
+++ b/src/Coxeter_triangulation/test/manifold_tracing_test.cpp
@@ -40,7 +40,7 @@ BOOST_AUTO_TEST_CASE(manifold_tracing) {
     BOOST_CHECK(si_pair.second.size() == (long int)oracle.function().amb_d());
   }
   std::cout << "out_simplex_map.size() = " << out_simplex_map.size() << "\n";
-  BOOST_CHECK(out_simplex_map.size() == 1140);
+  BOOST_CHECK(out_simplex_map.size() == 1118);
 
   // manifold with boundary
   Function_Sm_in_Rd fun_boundary(3.0, 2, fun_sph.seed());
@@ -58,5 +58,5 @@ BOOST_AUTO_TEST_CASE(manifold_tracing) {
     BOOST_CHECK(si_pair.second.size() == (long int)oracle.function().amb_d());
   }
   std::cout << "boundary_simplex_map.size() = " << boundary_simplex_map.size() << "\n";
-  BOOST_CHECK(boundary_simplex_map.size() == 55);
+  BOOST_CHECK(boundary_simplex_map.size() == 54);
 }