From 52fcaa530de936f2c9dfc076bd662c72300f796e Mon Sep 17 00:00:00 2001 From: SatoshiTerasaki Date: Wed, 11 Dec 2024 00:30:23 +0900 Subject: [PATCH 1/5] Define default constructor --- include/sparseir/gauss.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/include/sparseir/gauss.hpp b/include/sparseir/gauss.hpp index aa55025..f8a26de 100644 --- a/include/sparseir/gauss.hpp +++ b/include/sparseir/gauss.hpp @@ -46,6 +46,8 @@ class Rule { // TODO: Define x and w as Eigen::VectorXd std::vector x, w, x_forward, x_backward; T a, b; + // Default constructor + Rule() {}; Rule(const std::vector& x, const std::vector& w, std::vector x_forward, std::vector x_backward, T a = -1, T b = 1) : x(x), w(w), x_forward(x_forward), x_backward(x_backward), a(a), b(b) {} Rule(const std::vector& x, const std::vector& w, T a = -1, T b = 1) From dc396bdc9f286ec97d4db383375eca06329d10bf Mon Sep 17 00:00:00 2001 From: SatoshiTerasaki Date: Wed, 11 Dec 2024 00:30:37 +0900 Subject: [PATCH 2/5] Add comments --- include/sparseir/gauss.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/include/sparseir/gauss.hpp b/include/sparseir/gauss.hpp index f8a26de..d42ae8a 100644 --- a/include/sparseir/gauss.hpp +++ b/include/sparseir/gauss.hpp @@ -48,8 +48,10 @@ class Rule { T a, b; // Default constructor Rule() {}; + // Constructor with x, w, x_forward, x_backward, a, b Rule(const std::vector& x, const std::vector& w, std::vector x_forward, std::vector 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& x, const std::vector& w, T a = -1, T b = 1) : x(x), w(w), a(a), b(b) { this->x_forward = x_forward.empty() ? std::vector(x.size(), 0) : x_forward; From 0bd46882beb728e0755536a5ff6390348ce97b8b Mon Sep 17 00:00:00 2001 From: SatoshiTerasaki Date: Wed, 11 Dec 2024 00:31:05 +0900 Subject: [PATCH 3/5] don't use shared_ptr --- include/sparseir/sve.hpp | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/include/sparseir/sve.hpp b/include/sparseir/sve.hpp index 2dd7469..69980a3 100644 --- a/include/sparseir/sve.hpp +++ b/include/sparseir/sve.hpp @@ -130,7 +130,7 @@ class AbstractSVE { template class SamplingSVE : public AbstractSVE { public: - std::shared_ptr kernel; + K kernel; double epsilon; int n_gauss; int nsvals_hint; @@ -142,25 +142,25 @@ class SamplingSVE : public AbstractSVE { Rule gauss_y; // Constructor - SamplingSVE(std::shared_ptr 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(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); + n_gauss = (n_gauss_ > 0) ? n_gauss_ : sve_hints(kernel, epsilon).ngauss(); + rule = convert(legendre(n_gauss)); + auto hints = sve_hints(kernel, epsilon); // TODO: This is not correct, we need to use the segments from the hints nsvals_hint = hints.nsvals(); segs_x = hints.template segments_x(); segs_y = hints.template segments_y(); + gauss_x = rule.piecewise(segs_x); + gauss_y = rule.piecewise(segs_y); } // Compute matrices for SVD std::vector> matrices() const override { std::vector> mats; - Eigen::MatrixX A = matrix_from_gauss(*kernel, gauss_x, gauss_y); + Eigen::MatrixX 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 { PiecewiseLegendrePolyVector ulx(polyvec_u); PiecewiseLegendrePolyVector vly(polyvec_v); canonicalize(ulx, vly); - return SVEResult(ulx, s, vly, *kernel, epsilon); + return SVEResult(ulx, s, vly, kernel, epsilon); } }; From 9afd981f580b6747fa730ceb1d937deeecd176b1 Mon Sep 17 00:00:00 2001 From: SatoshiTerasaki Date: Wed, 11 Dec 2024 00:31:17 +0900 Subject: [PATCH 4/5] create an instance of SamplingSVE --- test/sve.cxx | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/test/sve.cxx b/test/sve.cxx index 12893f6..0788565 100644 --- a/test/sve.cxx +++ b/test/sve.cxx @@ -58,7 +58,18 @@ void check_smooth(const std::function& u, const std::vector segs_x = hints.template segments_x(); + std::vector segs_y = hints.template segments_y(); + + // Ensure `convert` is declared before this line + sparseir::Rule rule = sparseir::convert(sparseir::legendre(n_gauss)); + + sparseir::Rule gauss_x = rule.piecewise(segs_x); + sparseir::Rule gauss_y = rule.piecewise(segs_y); + sparseir::SamplingSVE sve(lk, 1e-6); //REQUIRE(true); // auto sve = sparseir::SamplingSVE(std::make_shared(sparseir::LogisticKernel(10.0)), 1e-6); // REQUIRE(sve.nsvals_hint == 10); From 5b914f85b498c268b4129c289fa51a7009d79e1c Mon Sep 17 00:00:00 2001 From: SatoshiTerasaki Date: Wed, 11 Dec 2024 00:38:42 +0900 Subject: [PATCH 5/5] update --- include/sparseir/sve.hpp | 2 +- test/sve.cxx | 25 ++++++++++++++----------- 2 files changed, 15 insertions(+), 12 deletions(-) diff --git a/include/sparseir/sve.hpp b/include/sparseir/sve.hpp index 69980a3..3ac162c 100644 --- a/include/sparseir/sve.hpp +++ b/include/sparseir/sve.hpp @@ -147,9 +147,9 @@ class SamplingSVE : public AbstractSVE { epsilon(epsilon_) { n_gauss = (n_gauss_ > 0) ? n_gauss_ : sve_hints(kernel, epsilon).ngauss(); + // TODO: Implement Rule(n_gauss) rule = convert(legendre(n_gauss)); auto hints = sve_hints(kernel, epsilon); - // TODO: This is not correct, we need to use the segments from the hints nsvals_hint = hints.nsvals(); segs_x = hints.template segments_x(); segs_y = hints.template segments_y(); diff --git a/test/sve.cxx b/test/sve.cxx index 0788565..6f83c10 100644 --- 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& u, const std::vector& s, double uscale, double fudge_factor) { /* @@ -69,15 +67,20 @@ TEST_CASE("AbstractSVE", "[sve]"){ sparseir::Rule gauss_x = rule.piecewise(segs_x); sparseir::Rule gauss_y = rule.piecewise(segs_y); - sparseir::SamplingSVE sve(lk, 1e-6); - //REQUIRE(true); - // auto sve = sparseir::SamplingSVE(std::make_shared(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 ssve1 = sparseir::SamplingSVE(lk, 1e-6); + REQUIRE(ssve1.n_gauss == n_gauss); + auto ssve2 = sparseir::SamplingSVE(lk, 1e-6, 12); + REQUIRE(ssve2.n_gauss == 12); + + auto ssve1_double = sparseir::SamplingSVE(lk, 1e-6); + REQUIRE(ssve1_double.n_gauss == n_gauss); + auto ssve2_double = sparseir::SamplingSVE(lk, 1e-6, 12); + REQUIRE(ssve2_double.n_gauss == 12); + + //auto ssve1_ddouble = sparseir::SamplingSVE(lk, 1e-6); + //REQUIRE(ssve1_ddouble.n_gauss == n_gauss); + //auto ssve2_ddouble = sparseir::SamplingSVE(lk, 1e-6, 12); + //REQUIRE(ssve2_ddouble.n_gauss == 12); } TEST_CASE("sve.cpp", "[sve]")