Merge pull request #50 from SpM-lab/terasaki/fix-poly-roots

Resolve roots
SpM-lab · Dec 13, 2024 · f6a8854 · f6a8854
2 parents fe97bd9 + 1c695f6
commit f6a8854
Show file tree

Hide file tree

Showing 4 changed files with 31 additions and 85 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -25,6 +25,7 @@ project(
 
 # Eigen3
 find_package (Eigen3 3.4 REQUIRED NO_MODULE)
+find_package(Threads REQUIRED)
 
 # libxprec
 include(FetchContent)
@@ -52,7 +53,7 @@ set_target_properties(sparseir PROPERTIES
     CXX_STANDARD_REQUIRED ON
     )
 
-target_link_libraries(sparseir Eigen3::Eigen XPrec::xprec)
+target_link_libraries(sparseir Eigen3::Eigen XPrec::xprec Threads::Threads)
 
 # Use library convention.
 add_library(SparseIR::sparseir ALIAS sparseir)

diff --git a/include/sparseir/_root.hpp b/include/sparseir/_root.hpp
@@ -235,7 +235,7 @@ std::vector<double> discrete_extrema(F f, const std::vector<double> &xgrid)
     std::vector<double> res;
     for (size_t i = 0; i < a.size(); ++i) {
         // abs ∘ f
-        auto abf = [f](double x) { return std::abs(f(x)); };
+        auto abf = [f](double x) { return std::fabs(f(x)); };
         res.push_back(
             bisect_discr_extremum(abf, a[i], b[i], absf_a[i], absf_b[i]));
     }

diff --git a/include/sparseir/poly.hpp b/include/sparseir/poly.hpp
@@ -146,7 +146,9 @@ class PiecewiseLegendrePoly {
         double x_tilde;
         std::tie(i, x_tilde) = split(x);
         Eigen::VectorXd coeffs = data.col(i);
-        double value = legval(x_tilde, coeffs) * norms[i];
+        // convert coeffs to std::vector<double>
+        std::vector<double> coeffs_vec(coeffs.data(), coeffs.data() + coeffs.size());
+        double value = legval<double>(x_tilde, coeffs_vec) * norms[i];
         return value;
     }
 
@@ -215,7 +217,7 @@ class PiecewiseLegendrePoly {
 
 
     Eigen::VectorXd refine_grid(const Eigen::VectorXd& grid, int alpha) const {
-        Eigen::VectorXd refined(grid.size() * alpha);
+        Eigen::VectorXd refined((grid.size() - 1) * alpha + 1);
 
         for (size_t i = 0; i < grid.size() - 1; ++i) {
             double start = grid[i];
@@ -245,80 +247,19 @@ class PiecewiseLegendrePoly {
         }
     }
 
-    // Equivalent to find_all function in julia above
-    Eigen::VectorXd find_all(const Eigen::VectorXd& xgrid) const {
-        Eigen::VectorXd fx(xgrid.size());
-        for (size_t i = 0; i < xgrid.size(); ++i) {
-            fx(i) = static_cast<double>((*this)(xgrid[i]));
-        }
-        // Find direct hits (zeros)
-        std::vector<bool> hit(fx.size(), false);
-        std::vector<double> x_hit;
-        for (Eigen::Index i = 0; i < fx.size(); ++i) {
-            hit[i] = (fx(i) == 0.0);
-            if (hit[i]) {
-                x_hit.push_back(xgrid(i));
-            }
-        }
-
-        // Check for sign changes
-        std::vector<bool> sign_change(fx.size() - 1, false);
-        bool found_sign_change = false;
-
-        for (Eigen::Index i = 0; i < fx.size() - 1; ++i) {
-            bool different_signs = std::signbit(fx(i)) != std::signbit(fx(i + 1));
-            bool neither_zero = fx(i) != 0.0 && fx(i + 1) != 0.0;
-            sign_change[i] = different_signs && neither_zero;
-            if (sign_change[i]) {
-                found_sign_change = true;
-            }
-        }
-
-        if (!found_sign_change) {
-            // convert to Eigen::VectorXd
-            return Eigen::Map<Eigen::VectorXd>(x_hit.data(), x_hit.size());
-        }
-
-        // Collect points for bisection
-        std::vector<double> a, b;
-        std::vector<double> fa;
-
-        for (size_t i = 0; i < sign_change.size(); ++i) {
-            if (sign_change[i]) {
-                a.push_back(xgrid(i));
-                b.push_back(xgrid(i + 1));
-                fa.push_back(fx(i));
-            }
-        }
-
-        // Calculate epsilon for floating point types
-        double eps_x;
-        if (std::is_floating_point<double>::value) {
-            eps_x = std::numeric_limits<double>::epsilon() * xgrid.cwiseAbs().maxCoeff();
-        } else {
-            eps_x = 0;
-        }
-
-        // Perform bisection for each interval
-        std::vector<double> x_bisect;
-        for (size_t i = 0; i < a.size(); ++i) {
-            x_bisect.push_back(bisect(a[i], b[i], fa[i], eps_x));
-        }
-
-        // Combine and sort results
-        x_hit.insert(x_hit.end(), x_bisect.begin(), x_bisect.end());
-        std::sort(x_hit.begin(), x_hit.end());
-        // convert to Eigen::VectorXd
-        return Eigen::Map<Eigen::VectorXd>(x_hit.data(), x_hit.size());
-    }
-
     // Roots function
     Eigen::VectorXd roots(double tol = 1e-10) const
     {
         Eigen::VectorXd grid = this->knots;
+
         Eigen::VectorXd refined_grid = refine_grid(grid, 2);
-        std::cout << "refined_grid: " << refined_grid.size() << "\n";
-        return find_all(refined_grid);
+        auto f = [this](double x) { return this->operator()(x); };
+        // convert to std::vector<double>
+        std::vector<double> refined_grid_vec(refined_grid.data(),
+                                             refined_grid.data() +
+                                                 refined_grid.size());
+        std::vector<double> roots = find_all(f, refined_grid_vec);
+        return Eigen::Map<Eigen::VectorXd>(roots.data(), roots.size());
     }
 
     // Overloaded operators
@@ -366,6 +307,7 @@ class PiecewiseLegendrePoly {
     int get_polyorder() const { return polyorder; }
 
 private:
+    /*
     // Helper function to compute legval
     static double legval(double x, const Eigen::VectorXd &coeffs)
     {
@@ -385,6 +327,7 @@ class PiecewiseLegendrePoly {
         }
         return result;
     }
+    */
 
     // Helper function to split x into segment index i and x_tilde
     std::pair<int, double> split(double x) const

diff --git a/test/poly.cxx b/test/poly.cxx
@@ -299,8 +299,8 @@ TEST_CASE("Overlap")
 TEST_CASE("Roots")
 {
     // Initialize data and knots (from Julia code)
-    Eigen::MatrixXd data(16, 2);
-    data << 0.16774734206553019, 0.49223680914312595, -0.8276728567928646,
+    Eigen::VectorXd v(32);
+    v << 0.16774734206553019, 0.49223680914312595, -0.8276728567928646,
         0.16912891046582143, -0.0016231275318572044, 0.00018381683946452256,
         -9.699355027805034e-7, 7.60144228530804e-8, -2.8518324490258146e-10,
         1.7090590205708293e-11, -5.0081401126025e-14, 2.1244236198427895e-15,
@@ -311,32 +311,34 @@ TEST_CASE("Roots")
         2.8518324490258146e-10, 1.7090590205708293e-11, 5.0081401126025e-14,
         2.1244236198427895e-15, -2.0478095258000225e-16, -2.676573801530628e-16,
         -2.338165820094204e-16, -1.2050663212312096e-16;
-    //data.resize(16, 2);
+
+    Eigen::Map<Eigen::MatrixXd> data(v.data(), 16, 2);
 
     Eigen::VectorXd knots(3);
     knots << 0.0, 0.5, 1.0;
     int l = 3;
 
     sparseir::PiecewiseLegendrePoly pwlp(data, knots, l);
-    /*
+
     // Find roots
     Eigen::VectorXd roots = pwlp.roots();
 
-    // Print roots for debugging
-    std::cout << "Found roots: " << roots.transpose() << "\n";
-
+    /*
     // Expected roots from Julia
     Eigen::VectorXd expected_roots(3);
     expected_roots << 0.1118633448586015, 0.4999999999999998, 0.8881366551413985;
 
-    // REQUIRE(roots.size() == expected_roots.size());
+    REQUIRE(roots.size() == expected_roots.size());
     for(Eigen::Index i = 0; i < roots.size(); ++i) {
-        //REQUIRE(std::abs(roots[i] - expected_roots[i]) < 1e-10);
+        REQUIRE(std::fabs(roots[i] - expected_roots[i]) < 1e-10);
+        // Verifying the polynomial evaluates to zero
+        //at the roots with tight tolerance
+        REQUIRE(std::fabs(pwlp(roots[i])) < 1e-10);
         // Verify roots are in domain
-        //REQUIRE(roots[i] >= knots[0]);
-        //REQUIRE(roots[i] <= knots[knots.size()-1]);
+        REQUIRE(roots[i] >= knots[0]);
+        REQUIRE(roots[i] <= knots[knots.size() - 1]);
         // Verify these are actually roots
-        //REQUIRE(std::abs(pwlp(roots[i])) < 1e-10);
+        REQUIRE(std::abs(pwlp(roots[i])) < 1e-10);
     }
     */
 }