diff --git a/include/sparseir/basis.hpp b/include/sparseir/basis.hpp index 677abb8..013a802 100644 --- a/include/sparseir/basis.hpp +++ b/include/sparseir/basis.hpp @@ -95,22 +95,33 @@ class FiniteTempBasis : public AbstractBasis { this->accuracy = sve_result.s(s_.size() - 1) / sve_result_s0; } - auto uk = u_[0].knots; - auto vk = v_[0].knots; - - Eigen::VectorXd u_knots = (beta / 2) * (uk.array() + 1); - Eigen::VectorXd v_knots = wmax * vk; - - Eigen::VectorXd deltax4u = (beta / 2) * u_[0].get_delta_x(); - Eigen::VectorXd deltax4v = wmax * v_[0].get_delta_x(); - + /* + Port the following Julia code to C++: + # The polynomials are scaled to the new variables by transforming + the # knots according to: tau = β/2 * (x + 1), w = ωmax * y. Scaling # the + data is not necessary as the normalization is inferred. ωmax = Λ(kernel) / β + u_knots = (β / 2) .* (knots(u_) .+ 1) + v_knots = ωmax .* knots(v_) + u = PiecewiseLegendrePolyVector(u_, u_knots; Δx=(β / 2) .* Δx(u_), + symm=symm(u_)) v = PiecewiseLegendrePolyVector(v_, v_knots; Δx=ωmax .* + Δx(v_), symm=symm(v_)) + */ + + auto u_knots_ = u_.polyvec[0].knots; + auto v_knots_ = v_.polyvec[0].knots; + + Eigen::VectorXd u_knots = (beta / 2) * (u_knots_.array() + 1); + Eigen::VectorXd v_knots = wmax * v_knots_; + + Eigen::VectorXd deltax4u = (beta / 2) * u_.get_delta_x(); + Eigen::VectorXd deltax4v = wmax * v_.get_delta_x(); std::vector u_symm_vec; for (int i = 0; i < u_.size(); ++i) { - u_symm_vec.push_back(u_[i].get_symm()); + u_symm_vec.push_back(u_.polyvec[i].get_symm()); } std::vector v_symm_vec; for (int i = 0; i < v_.size(); ++i) { - v_symm_vec.push_back(v_[i].get_symm()); + v_symm_vec.push_back(v_.polyvec[i].get_symm()); } Eigen::VectorXi u_symm = Eigen::Map(u_symm_vec.data(), u_symm_vec.size()); @@ -275,7 +286,7 @@ inline Eigen::VectorXd default_sampling_points(const PiecewiseLegendrePolyVector if (L < u.size()) { // TODO: Resolve this errors. - return u[L].roots(); + return u.polyvec[L].roots(); } else { // Approximate roots by extrema // TODO: resolve this error diff --git a/include/sparseir/poly.hpp b/include/sparseir/poly.hpp index cdbe6be..e7c3890 100644 --- a/include/sparseir/poly.hpp +++ b/include/sparseir/poly.hpp @@ -114,15 +114,18 @@ class PiecewiseLegendrePoly { knots(knots), l(l), symm(symm) { - polyorder = data.rows(); + this->polyorder = data.rows(); int nsegments = data.cols(); if (knots.size() != nsegments + 1) { throw std::invalid_argument("Invalid knots array"); } - this->delta_x = delta_x_.size() > 0 ? delta_x_ : knots.tail(knots.size() - 1) - knots.head(knots.size() - 1); - this->xm = 0.5 * (knots.head(knots.size() - 1) + knots.tail(knots.size() - 1)); + this->delta_x = delta_x_.size() > 0 ? delta_x_ : diff(knots); + this->xm = Eigen::VectorXd::Zero(nsegments); + for (int i = 0; i < nsegments; ++i) { + this->xm[i] = 0.5 * (knots[i] + knots[i + 1]); + } this->inv_xs = 2.0 / this->delta_x.array(); this->norms = this->inv_xs.array().sqrt(); @@ -318,7 +321,6 @@ class PiecewiseLegendrePoly { const Eigen::VectorXd &get_norms() const { return norms; } int get_polyorder() const { return polyorder; } -private: /* // Helper function to compute legval static double legval(double x, const Eigen::VectorXd &coeffs) @@ -357,6 +359,7 @@ class PiecewiseLegendrePoly { return std::make_pair(i, x_tilde); } +private: // Placeholder for Gauss-Legendre quadrature over [a, b] double gauss_legendre_quadrature(double a, double b, std::function f) const @@ -669,11 +672,8 @@ class PiecewiseLegendrePolyVector { "Sizes of polys and symm don't match " + std::to_string(polys.size()) + " " + std::to_string(symm.size())); } for (size_t i = 0; i < polys.size(); ++i) { - polyvec[i] = PiecewiseLegendrePoly(polys[i].get_data(), - knots, - polys[i].l, - polys[i].get_delta_x(), - symm(i)); + polyvec[i] = PiecewiseLegendrePoly(polys[i].get_data(), knots, + polys[i].l, Δx, symm(i)); } } diff --git a/include/sparseir/sve.hpp b/include/sparseir/sve.hpp index ad0e80c..cf9d345 100644 --- a/include/sparseir/sve.hpp +++ b/include/sparseir/sve.hpp @@ -325,16 +325,8 @@ class SamplingSVE : public AbstractSVE { } } - // Manually compute differences for dsegs_x and dsegs_y - Eigen::VectorX dsegs_x(segs_x.size() - 1); - for (int i = 0; i < segs_x.size() - 1; ++i) { - dsegs_x[i] = segs_x[i + 1] - segs_x[i]; - } - - Eigen::VectorX dsegs_y(segs_y.size() - 1); - for (int i = 0; i < segs_y.size() - 1; ++i) { - dsegs_y[i] = segs_y[i + 1] - segs_y[i]; - } + auto dsegs_x = diff(segs_x); + auto dsegs_y = diff(segs_y); // Using nested for loops to multiply u_data for (int j = 0; j < u_data.dimension(1); ++j) { @@ -356,6 +348,17 @@ class SamplingSVE : public AbstractSVE { std::vector polyvec_u; std::vector polyvec_v; + std::vector segs_x_double(segs_x.size()); + std::vector segs_y_double(segs_y.size()); + + for (int i = 0; i < segs_x.size(); ++i) { + segs_x_double[i] = static_cast(segs_x[i]); + } + for (int i = 0; i < segs_y.size(); ++i) { + segs_y_double[i] = static_cast(segs_y[i]); + } + Eigen::VectorXd knots_x = Eigen::Map(segs_x_double.data(), segs_x_double.size()); + Eigen::VectorXd knots_y = Eigen::Map(segs_y_double.data(), segs_y_double.size()); for (int i = 0; i < u_data.dimension(2); ++i) { Eigen::MatrixXd slice_double(u_data.dimension(0), @@ -365,16 +368,15 @@ class SamplingSVE : public AbstractSVE { slice_double(j, k) = static_cast(u_data(j, k, i)); } } - std::vector segs_x_double; - segs_x_double.reserve(segs_x.size()); - for (const auto &x : segs_x) { - segs_x_double.push_back(static_cast(x)); - } - polyvec_u.push_back(PiecewiseLegendrePoly( - slice_double, - Eigen::VectorXd::Map(segs_x_double.data(), - segs_x_double.size()), - i)); + + polyvec_u.push_back( + PiecewiseLegendrePoly( + slice_double, + knots_x, + i, + diff(knots_x) + ) + ); } // Repeat similar changes for v_data @@ -387,22 +389,16 @@ class SamplingSVE : public AbstractSVE { } } - std::vector segs_y_double; - segs_y_double.reserve(segs_y.size()); - for (const auto &y : segs_y) { - segs_y_double.push_back(static_cast(y)); - } - - polyvec_v.push_back(PiecewiseLegendrePoly( - slice_double, - Eigen::VectorXd::Map(segs_y_double.data(), - segs_y_double.size()), - i)); + polyvec_v.push_back( + PiecewiseLegendrePoly( + slice_double, + knots_y, + i, + diff(knots_y) + ) + ); } - - - PiecewiseLegendrePolyVector ulx(polyvec_u); PiecewiseLegendrePolyVector vly(polyvec_v); canonicalize(ulx, vly); diff --git a/include/sparseir/utils.hpp b/include/sparseir/utils.hpp index e79e429..129f6ad 100644 --- a/include/sparseir/utils.hpp +++ b/include/sparseir/utils.hpp @@ -61,4 +61,22 @@ bool tensorIsApprox(const Eigen::Tensor &a, const Eigen::Tensor &b, return true; } -} \ No newline at end of file +template +inline std::vector diff(const std::vector &xs) { + std::vector diff(xs.size() - 1); + for (size_t i = 0; i < xs.size() - 1; ++i) { + diff[i] = xs[i+1] - xs[i]; + } + return diff; +} + +template +inline Eigen::VectorX diff(const Eigen::VectorX &xs) { + Eigen::VectorX diff(xs.size() - 1); + for (Eigen::Index i = 0; i < xs.size() - 1; ++i) { + diff[i] = xs[i+1] - xs[i]; + } + return diff; +} + +} // namespace sparseir \ No newline at end of file diff --git a/test/basis.cxx b/test/basis.cxx index 8c25f00..c2e7b9a 100644 --- a/test/basis.cxx +++ b/test/basis.cxx @@ -17,6 +17,109 @@ using namespace std; using ComplexF64 = std::complex; +TEST_CASE("basis.u[0] test", "[basis]") +{ + /* + using SparseIR + begin + β = 1.0 + Λ = 10. + sve_result = SparseIR.SVEResult(LogisticKernel(Λ), ε=1e-15) + basis = FiniteTempBasis{Bosonic}(β, Λ; sve_result) + u = basis.u[begin] + end + */ + double beta = 1.0; + double Lambda = 10.0; + auto kernel = LogisticKernel(beta * Lambda); + auto sve_result = compute_sve(kernel, 1e-15); + auto basis = make_shared>(beta, Lambda, 1e-15, + kernel, sve_result); + + std::vector s_ref_vec = { + 1.2621489299919293, 0.8363588547029699, 0.3462622585830318, + 0.12082626967769121, 0.03387861935965415, 0.00796130085778543, + 0.0015966925515801561, 0.00027823051505153205, 4.276437930624593e-5, + 5.871774564004103e-6, 7.279433734090833e-7, 8.221881611234219e-8, + 8.525219704563244e-9, 8.168448057712933e-10, 7.273189647675939e-11, + 6.0477895415959716e-12, 4.716593209003674e-13, 3.4631886072022945e-14, + 2.4022217530858486e-15, + }; + Eigen::Map s(basis->s.data(), basis->s.size()); + REQUIRE(s.size() == s_ref_vec.size()); + Eigen::VectorXd s_ref = + Eigen::Map(s_ref_vec.data(), s_ref_vec.size()); + REQUIRE(s.isApprox(s_ref)); + + auto u = basis->u[0]; + + REQUIRE(u.xmin == 0.0); + REQUIRE(u.xmax == 1.0); + REQUIRE(u.get_polyorder() == 16); + REQUIRE(u.l == 0); + REQUIRE(u.symm == 1); + + std::vector u_knot_ref_vec = { + 0.0, + 0.013623446212733203, + 0.0292485379483044, + 0.04713527623719438, + 0.0675600047521634, + 0.09080712995504114, + 0.11715549463550101, + 0.14685782470149988, + 0.18011200903273533, + 0.21702410043575965, + 0.25756521020520734, + 0.3015280008924536, + 0.3484926000123403, + 0.3978146340870691, + 0.44864700611131486, + 0.5, + 0.5513529938886851, + 0.6021853659129309, + 0.6515073999876597, + 0.6984719991075464, + 0.7424347897947927, + 0.7829758995642404, + 0.8198879909672647, + 0.8531421752985001, + 0.882844505364499, + 0.9091928700449589, + 0.9324399952478366, + 0.9528647237628056, + 0.9707514620516956, + 0.9863765537872669, + 1.0, + }; + Eigen::Map u_knots(u.knots.data(), u.knots.size()); + REQUIRE(u_knots.isApprox(Eigen::Map( + u_knot_ref_vec.data(), u_knot_ref_vec.size()))); + + std::vector u_delta_x_ref_vec = { + 0.013623446212733203, 0.015625091735571195, 0.01788673828888998, + 0.020424728514969015, 0.02324712520287775, 0.026348364680459868, + 0.029702330065998872, 0.03325418433123545, 0.036912091403024316, + 0.04054110976944769, 0.043962790687246206, 0.04696459911988671, + 0.049322034074728835, 0.050832372024245794, 0.05135299388868512, + 0.05135299388868512, 0.050832372024245794, 0.049322034074728835, + 0.04696459911988671, 0.043962790687246206, 0.04054110976944769, + 0.036912091403024316, 0.03325418433123545, 0.029702330065998872, + 0.026348364680459868, 0.02324712520287775, 0.020424728514969015, + 0.01788673828888998, 0.015625091735571195, 0.013623446212733203, + }; + Eigen::Map u_delta_x(u.delta_x.data(), u.delta_x.size()); + REQUIRE(u_delta_x.isApprox(Eigen::Map( + u_delta_x_ref_vec.data(), u_delta_x_ref_vec.size()))); + + REQUIRE(u(0.3) == Approx(0.8209004724107448)); + int i; + double x_tilde; + std::tie(i, x_tilde) = u.split(0.3); + REQUIRE(i == 10); + REQUIRE(x_tilde == Approx(0.9304866288710429)); +} + TEST_CASE("basis.u(x)", "[basis]") { /* # Julia implementation diff --git a/test/sampling.cxx b/test/sampling.cxx index f85ec7f..38045a9 100644 --- a/test/sampling.cxx +++ b/test/sampling.cxx @@ -3,6 +3,7 @@ #include +#include // for Approx #include #include @@ -12,6 +13,7 @@ using namespace sparseir; using namespace std; +using Catch::Approx; using ComplexF64 = std::complex; @@ -39,11 +41,11 @@ TEST_CASE("TauSampling test", "[sampling]") { }; Eigen::MatrixXd mat_ref = Eigen::Map(mat_ref_vec.data(), 1, 19); Eigen::MatrixXd mat = eval_matrix(&tau_sampling, basis, x); - std::cout << "mat: \n" << mat << std::endl; - std::cout << "mat_ref: \n" << mat_ref << std::endl; + + REQUIRE(basis->u[0](0.3) == Approx(0.8209004724107448)); + REQUIRE(basis->u(x).transpose().isApprox(mat)); - // REQUIRE(basis->u(x).transpose().isApprox(mat_ref)); - // Rest of the test case... + REQUIRE(basis->u(x).transpose().isApprox(mat_ref)); } TEST_CASE("complex_test") { @@ -157,7 +159,7 @@ TEST_CASE("complex_test") { } auto tau_sampling = TauSampling(basis); - + /* { int dim = 0; Eigen::Tensor gl = movedim(originalgl, 0, dim); @@ -174,6 +176,7 @@ TEST_CASE("complex_test") { std::cout << "gl_from_tau(0, 0), gl(0, 0): " << gl_from_tau(0, 0) << ", " << gl(0, 0) << std::endl; std::cout << "gl_from_tau(0, 1), gl(0, 1): " << gl_from_tau(0, 1) << ", " << gl(0, 1) << std::endl; } + */ /* { int dim = 1; diff --git a/test/utils.cxx b/test/utils.cxx index f8a0b96..107d051 100644 --- a/test/utils.cxx +++ b/test/utils.cxx @@ -43,10 +43,34 @@ TEST_CASE("issorted", "[utils]") REQUIRE(sparseir::issorted(ev1) == true); } - using ComplexF64 = std::complex; +TEST_CASE("diff", "[utils]") +{ + SECTION("std::vector") + { + std::vector xs = {1.0, 3.0, 9.0, 1.0, 2.0}; + REQUIRE(sparseir::diff(xs) == std::vector({2.0, 6.0, -8.0, 1.0})); + } + SECTION("Eigen::VectorXd") + { + Eigen::VectorXd ev1(5); + ev1 << 1.0, 3.0, 9.0, 1.0, 2.0; + Eigen::VectorXd expected(4); + expected << 2.0, 6.0, -8.0, 1.0; + REQUIRE(sparseir::diff(ev1) == expected); + } + SECTION("Eigen::VectorX") + { + Eigen::VectorX ev1(5); + ev1 << 1.0, 3.0, 9.0, 1.0, 2.0; + Eigen::VectorX expected(4); + expected << 2.0, 6.0, -8.0, 1.0; + REQUIRE(sparseir::diff(ev1) == expected); + } +} +using ComplexF64 = std::complex; TEST_CASE("movedim", "[utils]"){ const Eigen::Index s_size = 12;