From 1e629f108401f7395385f4bf58b108cf15e25509 Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Wed, 20 Nov 2024 18:30:53 +0900 Subject: [PATCH 1/5] [WIP] Implement PiecewiseLegendrePolyVector --- include/sparseir/poly.hpp | 60 +++++++++++++++++++++++++-------------- 1 file changed, 38 insertions(+), 22 deletions(-) diff --git a/include/sparseir/poly.hpp b/include/sparseir/poly.hpp index dda8ec0..d3cdb09 100644 --- a/include/sparseir/poly.hpp +++ b/include/sparseir/poly.hpp @@ -315,15 +315,15 @@ class PiecewiseLegendrePoly { } // namespace sparseir -/* namespace sparseir { - class PiecewiseLegendrePolyVector { + +class PiecewiseLegendrePolyVector { public: // Member variable std::vector polyvec; // Constructors - PiecewiseLegendrePolyVector() {} + // PiecewiseLegendrePolyVector() {} // Constructor with polyvec PiecewiseLegendrePolyVector(const std::vector& polyvec) @@ -338,27 +338,27 @@ namespace sparseir { if (!symm.empty() && symm.size() != npolys) { throw std::runtime_error("Sizes of data and symm don't match"); } + + + int nrows = data.dimension(0); + int ncols = data.dimension(1); + polyvec.reserve(npolys); for (int i = 0; i < npolys; ++i) { - Eigen::MatrixXd data_i = data.chip(i, 2); - int sym = symm.empty() ? 0 : symm[i]; - polyvec.emplace_back(data_i, knots, i, Eigen::VectorXd(), sym); - } - } + // Evaluate the tensor slice + Eigen::Tensor data_i_tensor = data.chip(i, 2).eval(); - // Constructor with polys, new knots, delta_x, and symm - PiecewiseLegendrePolyVector(const PiecewiseLegendrePolyVector& polys, - const Eigen::VectorXd& knots, - const Eigen::VectorXd& delta_x = Eigen::VectorXd(), - const std::vector& symm = std::vector()) - { - if (!symm.empty() && symm.size() != polys.size()) { - throw std::runtime_error("Sizes of polys and symm don't match"); + // Create an Eigen::MatrixXd + Eigen::MatrixXd data_i(nrows, ncols); + + // Copy data from the tensor to the matrix + for (int r = 0; r < nrows; ++r) { + for (int c = 0; c < ncols; ++c) { + data_i(r, c) = data_i_tensor(r, c); + } } - polyvec.reserve(polys.size()); - for (size_t i = 0; i < polys.size(); ++i) { - int sym = symm.empty() ? 0 : symm[i]; - polyvec.emplace_back(polys.polyvec[i].data, knots, polys.polyvec[i].l, delta_x, sym); + int sym = symm.empty() ? 0 : symm[i]; + polyvec.emplace_back(data_i, knots, i, Eigen::VectorXd(), sym); } } @@ -371,8 +371,24 @@ namespace sparseir { throw std::runtime_error("Sizes of data and polys don't match"); } polyvec.reserve(npolys); + + int nrows = data.dimension(0); + int ncols = data.dimension(1); + for (int i = 0; i < npolys; ++i) { - Eigen::MatrixXd data_i = data.chip(i, 2); + // Evaluate the tensor slice + Eigen::Tensor data_i_tensor = data.chip(i, 2).eval(); + + // Create an Eigen::MatrixXd + Eigen::MatrixXd data_i(nrows, ncols); + + // Copy data from the tensor to the matrix + for (int r = 0; r < nrows; ++r) { + for (int c = 0; c < ncols; ++c) { + data_i(r, c) = data_i_tensor(r, c); + } + } + // TODO: fix me. polyvec.emplace_back(data_i, polys.polyvec[i]); } } @@ -460,4 +476,4 @@ namespace sparseir { } }; } // namespace sparseir -*/ + From c647ee9425acaf08b587c83f5afd07c6c266fc2a Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Fri, 22 Nov 2024 11:36:08 +0900 Subject: [PATCH 2/5] Resolve Test properties --- include/sparseir/poly.hpp | 6 ++++-- test/poly.cxx | 2 +- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/include/sparseir/poly.hpp b/include/sparseir/poly.hpp index d3cdb09..dcd3745 100644 --- a/include/sparseir/poly.hpp +++ b/include/sparseir/poly.hpp @@ -314,7 +314,6 @@ class PiecewiseLegendrePoly { } // namespace sparseir - namespace sparseir { class PiecewiseLegendrePolyVector { @@ -329,6 +328,7 @@ class PiecewiseLegendrePolyVector { PiecewiseLegendrePolyVector(const std::vector& polyvec) : polyvec(polyvec) {} + /* // Constructor with data tensor, knots, and optional symm vector PiecewiseLegendrePolyVector(const Eigen::Tensor& data, const Eigen::VectorXd& knots, @@ -361,7 +361,6 @@ class PiecewiseLegendrePolyVector { polyvec.emplace_back(data_i, knots, i, Eigen::VectorXd(), sym); } } - // Constructor with data tensor and existing polys PiecewiseLegendrePolyVector(const Eigen::Tensor& data, const PiecewiseLegendrePolyVector& polys) @@ -392,6 +391,7 @@ class PiecewiseLegendrePolyVector { polyvec.emplace_back(data_i, polys.polyvec[i]); } } + */ // Accessors size_t size() const { return polyvec.size(); } @@ -419,6 +419,7 @@ class PiecewiseLegendrePolyVector { } return symms; } + /* // Function to retrieve data as Eigen Tensor Eigen::Tensor get_data() const { @@ -474,6 +475,7 @@ class PiecewiseLegendrePolyVector { os << "on [" << polys.xmin() << ", " << polys.xmax() << "]"; return os; } + */ }; } // namespace sparseir diff --git a/test/poly.cxx b/test/poly.cxx index afbd2f5..3807e95 100644 --- a/test/poly.cxx +++ b/test/poly.cxx @@ -205,7 +205,6 @@ TEST_CASE("sparseir::PiecewiseLegendrePolyVector") { // Create sparseir::PiecewiseLegendrePolyVector std::vector polyvec = {pwlp1, pwlp2, pwlp3}; - /* sparseir::PiecewiseLegendrePolyVector polys(polyvec); // Test length @@ -223,6 +222,7 @@ TEST_CASE("sparseir::PiecewiseLegendrePolyVector") { std::vector expected_symm = {pwlp1.symm, pwlp2.symm, pwlp3.symm}; std::vector polys_symm = polys.get_symm(); REQUIRE(polys_symm == expected_symm); + /* // Test evaluation at a random point x double x = 0.5; // Example point From 84250550689a494c6de613d2291e11c1ed0cd817 Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Fri, 22 Nov 2024 11:38:33 +0900 Subject: [PATCH 3/5] Resolve tests for polys(x) --- include/sparseir/poly.hpp | 2 +- test/poly.cxx | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/include/sparseir/poly.hpp b/include/sparseir/poly.hpp index dcd3745..c89bf1b 100644 --- a/include/sparseir/poly.hpp +++ b/include/sparseir/poly.hpp @@ -419,7 +419,6 @@ class PiecewiseLegendrePolyVector { } return symms; } - /* // Function to retrieve data as Eigen Tensor Eigen::Tensor get_data() const { @@ -448,6 +447,7 @@ class PiecewiseLegendrePolyVector { } return results; } + /* // Evaluate the vector of polynomials at multiple x Eigen::MatrixXd operator()(const Eigen::VectorXd& xs) const { diff --git a/test/poly.cxx b/test/poly.cxx index 3807e95..4df00c0 100644 --- a/test/poly.cxx +++ b/test/poly.cxx @@ -222,7 +222,6 @@ TEST_CASE("sparseir::PiecewiseLegendrePolyVector") { std::vector expected_symm = {pwlp1.symm, pwlp2.symm, pwlp3.symm}; std::vector polys_symm = polys.get_symm(); REQUIRE(polys_symm == expected_symm); - /* // Test evaluation at a random point x double x = 0.5; // Example point @@ -230,6 +229,7 @@ TEST_CASE("sparseir::PiecewiseLegendrePolyVector") { Eigen::VectorXd expected_x(3); expected_x << pwlp1(x), pwlp2(x), pwlp3(x); REQUIRE(polys_x.isApprox(expected_x)); + /* // Test data Eigen::Tensor data_tensor = polys.get_data(); From 98c881b92cdb60d838377afba23338b107824142 Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Fri, 22 Nov 2024 11:42:26 +0900 Subject: [PATCH 4/5] Resolve tests in "sparseir::PiecewiseLegendrePolyVector" --- include/sparseir/poly.hpp | 2 +- test/poly.cxx | 3 +-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/include/sparseir/poly.hpp b/include/sparseir/poly.hpp index c89bf1b..9cb0ea1 100644 --- a/include/sparseir/poly.hpp +++ b/include/sparseir/poly.hpp @@ -447,7 +447,6 @@ class PiecewiseLegendrePolyVector { } return results; } - /* // Evaluate the vector of polynomials at multiple x Eigen::MatrixXd operator()(const Eigen::VectorXd& xs) const { @@ -457,6 +456,7 @@ class PiecewiseLegendrePolyVector { } return results; } + /* // Overlap function std::vector overlap(std::function f, double rtol = std::numeric_limits::epsilon(), diff --git a/test/poly.cxx b/test/poly.cxx index 4df00c0..bf876ea 100644 --- a/test/poly.cxx +++ b/test/poly.cxx @@ -229,7 +229,6 @@ TEST_CASE("sparseir::PiecewiseLegendrePolyVector") { Eigen::VectorXd expected_x(3); expected_x << pwlp1(x), pwlp2(x), pwlp3(x); REQUIRE(polys_x.isApprox(expected_x)); - /* // Test data Eigen::Tensor data_tensor = polys.get_data(); @@ -240,12 +239,12 @@ TEST_CASE("sparseir::PiecewiseLegendrePolyVector") { Eigen::VectorXd xs(4); xs << -0.8, -0.2, 0.2, 0.8; Eigen::MatrixXd polys_xs = polys(xs); // Should return a matrix of size (3, 4) + Eigen::MatrixXd expected_xs(3, 4); for (int i = 0; i < 4; ++i) { expected_xs.col(i) << pwlp1(xs[i]), pwlp2(xs[i]), pwlp3(xs[i]); } REQUIRE(polys_xs.isApprox(expected_xs)); - */ } TEST_CASE("Deriv") { From ea9995e04510dedd6b515ddfc79f164dd8634321 Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Fri, 22 Nov 2024 11:44:05 +0900 Subject: [PATCH 5/5] arrange headers --- include/sparseir/poly.hpp | 22 +++++----------------- 1 file changed, 5 insertions(+), 17 deletions(-) diff --git a/include/sparseir/poly.hpp b/include/sparseir/poly.hpp index 9cb0ea1..1ee0faf 100644 --- a/include/sparseir/poly.hpp +++ b/include/sparseir/poly.hpp @@ -1,27 +1,15 @@ -// Include necessary headers -#include -#include -#include +// C++ Standard Library headers #include -#include -#include #include +#include #include - -#include -#include #include -#include -#include #include -#include -#include -#include -#include -#include #include -#include +#include +// Eigen headers +#include #include