Merge pull request #39 from SpM-lab/terasaki/write-test-samplingsve

Write tests related to SamplingSVE
SpM-lab · Dec 10, 2024 · 6c55333 · 6c55333
2 parents 859a674 + 5b914f8
commit 6c55333
Show file tree

Hide file tree

Showing 3 changed files with 40 additions and 22 deletions.
diff --git a/include/sparseir/gauss.hpp b/include/sparseir/gauss.hpp
@@ -46,8 +46,12 @@ class Rule {
     // TODO: Define x and w as Eigen::VectorXd
     std::vector<T> x, w, x_forward, x_backward;
     T a, b;
+    // Default constructor
+    Rule() {};
+    // Constructor with x, w, x_forward, x_backward, a, b
     Rule(const std::vector<T>& x, const std::vector<T>& w, std::vector<T> x_forward, std::vector<T> x_backward, T a = -1, T b = 1)
         : x(x), w(w), x_forward(x_forward), x_backward(x_backward), a(a), b(b) {}
+    // Constructor with x, w, a, b
     Rule(const std::vector<T>& x, const std::vector<T>& w, T a = -1, T b = 1)
         : x(x), w(w), a(a), b(b) {
         this->x_forward = x_forward.empty() ? std::vector<T>(x.size(), 0) : x_forward;

diff --git a/include/sparseir/sve.hpp b/include/sparseir/sve.hpp
@@ -130,7 +130,7 @@ class AbstractSVE {
 template <typename K, typename T = double>
 class SamplingSVE : public AbstractSVE<K, T> {
 public:
-    std::shared_ptr<const K> kernel;
+    K kernel;
     double epsilon;
     int n_gauss;
     int nsvals_hint;
@@ -142,25 +142,25 @@ class SamplingSVE : public AbstractSVE<K, T> {
     Rule<T> gauss_y;
 
     // Constructor
-    SamplingSVE(std::shared_ptr<const K> kernel_, double epsilon_, int n_gauss_ = -1)
+    SamplingSVE(K kernel_, double epsilon_, int n_gauss_ = -1)
         : kernel(kernel_),
-          epsilon(epsilon_),
-          n_gauss((n_gauss_ > 0) ? n_gauss_ : sve_hints(*kernel, epsilon).ngauss()),
-          rule(convert<T>(legendre(n_gauss))),
-          gauss_x(rule.piecewise(segs_x)), // Initialize gauss_x
-          gauss_y(rule.piecewise(segs_y))  // Initialize gauss_y
+          epsilon(epsilon_)
     {
-        auto hints = sve_hints(*kernel, epsilon);
-        // TODO: This is not correct, we need to use the segments from the hints
+        n_gauss = (n_gauss_ > 0) ? n_gauss_ : sve_hints(kernel, epsilon).ngauss();
+        // TODO: Implement Rule<T>(n_gauss)
+        rule = convert<T>(legendre(n_gauss));
+        auto hints = sve_hints(kernel, epsilon);
         nsvals_hint = hints.nsvals();
         segs_x = hints.template segments_x<T>();
         segs_y = hints.template segments_y<T>();
+        gauss_x = rule.piecewise(segs_x);
+        gauss_y = rule.piecewise(segs_y);
     }
 
     // Compute matrices for SVD
     std::vector<Eigen::MatrixX<T>> matrices() const override {
         std::vector<Eigen::MatrixX<T>> mats;
-        Eigen::MatrixX<T> A = matrix_from_gauss(*kernel, gauss_x, gauss_y);
+        Eigen::MatrixX<T> A = matrix_from_gauss(kernel, gauss_x, gauss_y);
         // Element-wise multiplication with square roots of weights
         for (int i = 0; i < gauss_x.w.size(); ++i)
         {
@@ -304,7 +304,7 @@ class SamplingSVE : public AbstractSVE<K, T> {
         PiecewiseLegendrePolyVector ulx(polyvec_u);
         PiecewiseLegendrePolyVector vly(polyvec_v);
         canonicalize(ulx, vly);
-        return SVEResult<K>(ulx, s, vly, *kernel, epsilon);
+        return SVEResult<K>(ulx, s, vly, kernel, epsilon);
     }
 };
 

diff --git a/test/sve.cxx b/test/sve.cxx
@@ -13,8 +13,6 @@
 
 using std::invalid_argument;
 
-using xprec::DDouble;
-
 // Function to check smoothness
 void check_smooth(const std::function<double(double)>& u, const std::vector<double>& s, double uscale, double fudge_factor) {
     /*
@@ -58,15 +56,31 @@ void check_smooth(const std::function<double(double)>& u, const std::vector<doub
 
 TEST_CASE("AbstractSVE", "[sve]"){
     sparseir::LogisticKernel lk(10.0);
-    //auto sve = sparseir::compute_sve(lk);
-    //REQUIRE(true);
-    // auto sve = sparseir::SamplingSVE<sparseir::LogisticKernel, double>(std::make_shared<sparseir::LogisticKernel>(sparseir::LogisticKernel(10.0)), 1e-6);
-    // REQUIRE(sve.nsvals_hint == 10);
-    // REQUIRE(sve.n_gauss == 10);
-    // REQUIRE(sve.segs_x.size() == 10);
-    // REQUIRE(sve.segs_y.size() == 10);
-    // REQUIRE(sve.gauss_x.size() == 10);
-    // REQUIRE(sve.gauss_y.size() == 10);
+    auto hints = sparseir::sve_hints(lk, 1e-6);
+    int nsvals_hint = hints.nsvals();
+    int n_gauss = hints.ngauss();
+    std::vector<double> segs_x = hints.template segments_x<double>();
+    std::vector<double> segs_y = hints.template segments_y<double>();
+
+    // Ensure `convert` is declared before this line
+    sparseir::Rule<double> rule = sparseir::convert<double>(sparseir::legendre(n_gauss));
+
+    sparseir::Rule<double> gauss_x = rule.piecewise(segs_x);
+    sparseir::Rule<double> gauss_y = rule.piecewise(segs_y);
+    auto ssve1 = sparseir::SamplingSVE<sparseir::LogisticKernel>(lk, 1e-6);
+    REQUIRE(ssve1.n_gauss == n_gauss);
+    auto ssve2 = sparseir::SamplingSVE<sparseir::LogisticKernel>(lk, 1e-6, 12);
+    REQUIRE(ssve2.n_gauss == 12);
+
+    auto ssve1_double = sparseir::SamplingSVE<sparseir::LogisticKernel, double>(lk, 1e-6);
+    REQUIRE(ssve1_double.n_gauss == n_gauss);
+    auto ssve2_double = sparseir::SamplingSVE<sparseir::LogisticKernel, double>(lk, 1e-6, 12);
+    REQUIRE(ssve2_double.n_gauss == 12);
+
+    //auto ssve1_ddouble = sparseir::SamplingSVE<sparseir::LogisticKernel, xprec::DDouble>(lk, 1e-6);
+    //REQUIRE(ssve1_ddouble.n_gauss == n_gauss);
+    //auto ssve2_ddouble = sparseir::SamplingSVE<sparseir::LogisticKernel, xprec::DDouble>(lk, 1e-6, 12);
+    //REQUIRE(ssve2_ddouble.n_gauss == 12);
 }
 
 TEST_CASE("sve.cpp", "[sve]")