Fix tests related to _roots.cxx

.github/workflows/CI_cmake.yml

@@ -37,5 +37,7 @@ jobs:
       working-directory: ${{github.workspace}}/build
       # Execute tests defined by the CMake configuration.
       # See https://cmake.org/cmake/help/latest/manual/ctest.1.html for more detail
-      run: ctest -C ${{env.BUILD_TYPE}}
+      run: ./test/libsparseirtests
+      #run: ctest -C ${{env.BUILD_TYPE}}
+
 

include/sparseir/_root.hpp

@@ -1,26 +1,80 @@
 #pragma once
 
 #include <algorithm>
-#include <bitset>
 #include <cmath>
 #include <functional>
+#include <limits>
+#include <type_traits>
 #include <vector>
 
+#include <iostream>
+
+#include <Eigen/Dense>
+
 namespace sparseir {
 
-template <typename T>
-T midpoint(T lo, T hi)
-{
-    if (std::is_integral<T>::value) {
-        return lo + ((hi - lo) >> 1);
-    } else {
-        return lo + ((hi - lo) * static_cast<T>(0.5));
+// Midpoint function for floating-point types
+template<typename T1, typename T2>
+typename std::enable_if<
+    std::is_floating_point<T1>::value && std::is_floating_point<T2>::value,
+    typename std::common_type<T1, T2>::type
+>::type
+inline midpoint(T1 a, T2 b) {
+    typedef typename std::common_type<T1, T2>::type CommonType;
+    return static_cast<CommonType>(a) +
+           (static_cast<CommonType>(b) - static_cast<CommonType>(a)) *
+           static_cast<CommonType>(0.5);
+}
+
+// Midpoint function for integral types
+template<typename T>
+typename std::enable_if<std::is_integral<T>::value, T>::type
+inline midpoint(T a, T b) {
+    return a + ((b - a) / 2);
+}
+
+// Close enough function for floating-point types
+template<typename T>
+inline bool closeenough(T a, T b, T epsilon) {
+    return std::abs(a - b) <= epsilon;
+}
+
+template<typename T>
+inline bool closeenough(int a, int b, T _dummyepsilon) {
+    return a == b;
+}
+
+// Signbit function (handles both floating-point and integral types)
+template<typename T>
+inline bool signbit(T x) {
+    return x < static_cast<T>(0);
+}
+
+// Bisection method to find a root of function f in [a, b]
+template<typename F, typename T>
+T bisect(F f, T a, T b, T fa, T epsilon_x) {
+    while (true) {
+        T mid = midpoint(a, b);
+        assert(epsilon_x > 0);
+        if (closeenough(a, mid, epsilon_x)) {
+            return mid;
+        }
+        T fmid = static_cast<T>(f(mid));
+        if (signbit(fa) != signbit(fmid)) {
+            b = mid;
+        } else {
+            a = mid;
+            fa = fmid;
+        }
     }
 }
 
 template <typename F, typename T>
 std::vector<T> find_all(F f, const std::vector<T> &xgrid)
 {
+    if (xgrid.empty()) {
+        return {};
+    }
     std::vector<double> fx;
     std::transform(xgrid.begin(), xgrid.end(), std::back_inserter(fx),
                    [&](T x) { return static_cast<double>(f(x)); });
@@ -68,48 +122,23 @@ std::vector<T> find_all(F f, const std::vector<T> &xgrid)
         }
     }
 
-    double epsilon_x =
-        std::numeric_limits<T>::epsilon() *
-        *std::max_element(xgrid.begin(), xgrid.end(),
-                          [](T a, T b) { return std::abs(a) < std::abs(b); });
+    T max_elm = std::abs(xgrid[0]);
+    for (size_t i = 1; i < xgrid.size(); ++i) {
+        max_elm = std::max(max_elm, std::abs(xgrid[i]));
+    }
+    T epsilon_x =std::numeric_limits<T>::epsilon() * max_elm;
 
     std::vector<T> x_bisect;
     for (size_t i = 0; i < a.size(); ++i) {
-        x_bisect.push_back(bisect(f, a[i], b[i], fa[i], epsilon_x));
+        double root = bisect(f, a[i], b[i], fa[i], epsilon_x);
+        x_bisect.push_back(static_cast<T>(root));
     }
 
     x_hit.insert(x_hit.end(), x_bisect.begin(), x_bisect.end());
     std::sort(x_hit.begin(), x_hit.end());
     return x_hit;
 }
 
-template <typename F, typename T>
-T bisect(F f, T a, T b, double fa, double epsilon_x)
-{
-    while (true) {
-        T mid = midpoint(a, b);
-        if (closeenough(a, mid, epsilon_x))
-            return mid;
-        double fmid = f(mid);
-        if (std::signbit(fa) != std::signbit(fmid)) {
-            b = mid;
-        } else {
-            a = mid;
-            fa = fmid;
-        }
-    }
-}
-
-template <typename T>
-bool closeenough(T a, T b, double epsilon)
-{
-    if (std::is_floating_point<T>::value) {
-        return std::abs(a - b) <= epsilon;
-    } else {
-        return a == b;
-    }
-}
-
 template <typename T>
 std::vector<T> refine_grid(const std::vector<T> &grid, int alpha)
 {
@@ -130,18 +159,18 @@ std::vector<T> refine_grid(const std::vector<T> &grid, int alpha)
     return newgrid;
 }
 
-template <typename F, typename T>
-T bisect_discr_extremum(F absf, T a, T b, double absf_a, double absf_b)
+template <typename F>
+double bisect_discr_extremum(F absf, double a, double b, double absf_a, double absf_b)
 {
-    T d = b - a;
+    double d = b - a;
 
     if (d <= 1)
         return absf_a > absf_b ? a : b;
     if (d == 2)
         return a + 1;
 
-    T m = midpoint(a, b);
-    T n = m + 1;
+    double m = midpoint(a, b);
+    double n = m + 1;
     double absf_m = absf(m);
     double absf_n = absf(n);
 
@@ -152,36 +181,45 @@ T bisect_discr_extremum(F absf, T a, T b, double absf_a, double absf_b)
     }
 }
 
-template <typename F, typename T>
-std::vector<T> discrete_extrema(F f, const std::vector<T> &xgrid)
+template <typename F>
+std::vector<double> discrete_extrema(F f, const std::vector<double> &xgrid)
 {
     std::vector<double> fx(xgrid.size());
-    std::transform(xgrid.begin(), xgrid.end(), fx.begin(), f);
+    for (size_t i = 0; i < xgrid.size(); ++i) {
+        fx[i] = f(xgrid[i]);
+    }
 
     std::vector<double> absfx(fx.size());
-    std::transform(fx.begin(), fx.end(), absfx.begin(),
-                   [](double val) { return std::abs(val); });
+    for (size_t i = 0; i < fx.size(); ++i) {
+        absfx[i] = std::abs(fx[i]);
+    }
 
     std::vector<bool> signdfdx(fx.size() - 1);
     for (size_t i = 0; i < fx.size() - 1; ++i) {
-        signdfdx[i] = std::signbit(fx[i]) != std::signbit(fx[i + 1]);
+        signdfdx[i] = std::signbit(fx[i + 1] - fx[i]);
     }
 
     std::vector<bool> derivativesignchange(signdfdx.size() - 1);
     for (size_t i = 0; i < signdfdx.size() - 1; ++i) {
         derivativesignchange[i] = signdfdx[i] != signdfdx[i + 1];
     }
 
-    std::vector<bool> derivativesignchange_a(derivativesignchange.size() + 2,
-                                             false);
-    std::vector<bool> derivativesignchange_b(derivativesignchange.size() + 2,
-                                             false);
-    for (size_t i = 0; i < derivativesignchange.size(); ++i) {
-        derivativesignchange_a[i] = derivativesignchange[i];
-        derivativesignchange_b[i + 2] = derivativesignchange[i];
-    }
-
-    std::vector<T> a, b;
+    // create copy of derivativesignchange and add two false at the end
+    std::vector<bool> derivativesignchange_a(derivativesignchange);
+    derivativesignchange_a.push_back(false);
+    derivativesignchange_a.push_back(false);
+
+    std::vector<bool> derivativesignchange_b;
+    derivativesignchange_b.reserve(derivativesignchange.size() + 2);
+    derivativesignchange_b.push_back(false);
+    derivativesignchange_b.push_back(false);
+    derivativesignchange_b.insert(
+        derivativesignchange_b.end(),
+        derivativesignchange.begin(),
+        derivativesignchange.end()
+    );
+
+    std::vector<double> a, b;
     std::vector<double> absf_a, absf_b;
     for (size_t i = 0; i < derivativesignchange_a.size(); ++i) {
         if (derivativesignchange_a[i]) {
@@ -194,15 +232,21 @@ std::vector<T> discrete_extrema(F f, const std::vector<T> &xgrid)
         }
     }
 
-    std::vector<T> res;
+    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)); };
         res.push_back(
-            bisect_discr_extremum(f, a[i], b[i], absf_a[i], absf_b[i]));
+            bisect_discr_extremum(abf, a[i], b[i], absf_a[i], absf_b[i]));
     }
 
+    // We consider the outer points to be extrema if there is a decrease
+    // in magnitude or a sign change inwards
+
     std::vector<bool> sfx(fx.size());
-    std::transform(fx.begin(), fx.end(), sfx.begin(),
-                   [](double val) { return std::signbit(val); });
+    for (size_t i = 0; i < fx.size(); ++i) {
+        sfx[i] = std::signbit(fx[i]);
+    }
 
     if (absfx.front() > absfx[1] || sfx.front() != sfx[1]) {
         res.insert(res.begin(), xgrid.front());