Skip to content

Commit

Permalink
Merge pull request #105 from SpM-lab/terasaki/fix-poly
Browse files Browse the repository at this point in the history 
Fix constructor of FiniteTempBasis
  • Loading branch information
@terasakisatoshi
terasakisatoshi authored Feb 5, 2025
2 parents bffa894 + ae59f08 commit e16ad7e
Show file tree
Hide file tree
Showing 7 changed files with 217 additions and 62 deletions.
35 changes: 23 additions & 12 deletions include/sparseir/basis.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -95,22 +95,33 @@ class FiniteTempBasis : public AbstractBasis<S> {
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<int> 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<int> 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<Eigen::VectorXi>(u_symm_vec.data(), u_symm_vec.size());
Expand Down Expand Up @@ -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
Expand Down
18 changes: 9 additions & 9 deletions include/sparseir/poly.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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();

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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<double(double)> f) const
Expand Down Expand Up @@ -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));
}
}

Expand Down
64 changes: 30 additions & 34 deletions include/sparseir/sve.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -325,16 +325,8 @@ class SamplingSVE : public AbstractSVE<T> {
}
}

// Manually compute differences for dsegs_x and dsegs_y
Eigen::VectorX<T> 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<T> 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) {
Expand All @@ -356,6 +348,17 @@ class SamplingSVE : public AbstractSVE<T> {

std::vector<PiecewiseLegendrePoly> polyvec_u;
std::vector<PiecewiseLegendrePoly> polyvec_v;
std::vector<double> segs_x_double(segs_x.size());
std::vector<double> segs_y_double(segs_y.size());

for (int i = 0; i < segs_x.size(); ++i) {
segs_x_double[i] = static_cast<double>(segs_x[i]);
}
for (int i = 0; i < segs_y.size(); ++i) {
segs_y_double[i] = static_cast<double>(segs_y[i]);
}
Eigen::VectorXd knots_x = Eigen::Map<Eigen::VectorXd>(segs_x_double.data(), segs_x_double.size());
Eigen::VectorXd knots_y = Eigen::Map<Eigen::VectorXd>(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),
Expand All @@ -365,16 +368,15 @@ class SamplingSVE : public AbstractSVE<T> {
slice_double(j, k) = static_cast<double>(u_data(j, k, i));
}
}
std::vector<double> segs_x_double;
segs_x_double.reserve(segs_x.size());
for (const auto &x : segs_x) {
segs_x_double.push_back(static_cast<double>(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
Expand All @@ -387,22 +389,16 @@ class SamplingSVE : public AbstractSVE<T> {
}
}

std::vector<double> segs_y_double;
segs_y_double.reserve(segs_y.size());
for (const auto &y : segs_y) {
segs_y_double.push_back(static_cast<double>(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);
Expand Down
20 changes: 19 additions & 1 deletion include/sparseir/utils.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,4 +61,22 @@ bool tensorIsApprox(const Eigen::Tensor<T, N> &a, const Eigen::Tensor<T, N> &b,
return true;
}

}
template <typename T>
inline std::vector<T> diff(const std::vector<T> &xs) {
std::vector<T> diff(xs.size() - 1);
for (size_t i = 0; i < xs.size() - 1; ++i) {
diff[i] = xs[i+1] - xs[i];
}
return diff;
}

template <typename T>
inline Eigen::VectorX<T> diff(const Eigen::VectorX<T> &xs) {
Eigen::VectorX<T> 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
103 changes: 103 additions & 0 deletions test/basis.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,109 @@ using namespace std;

using ComplexF64 = std::complex<double>;

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<FiniteTempBasis<Bosonic>>(beta, Lambda, 1e-15,
kernel, sve_result);

std::vector<double> 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<Eigen::VectorXd> s(basis->s.data(), basis->s.size());
REQUIRE(s.size() == s_ref_vec.size());
Eigen::VectorXd s_ref =
Eigen::Map<Eigen::VectorXd>(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<double> 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<Eigen::VectorXd> u_knots(u.knots.data(), u.knots.size());
REQUIRE(u_knots.isApprox(Eigen::Map<Eigen::VectorXd>(
u_knot_ref_vec.data(), u_knot_ref_vec.size())));

std::vector<double> 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<Eigen::VectorXd> u_delta_x(u.delta_x.data(), u.delta_x.size());
REQUIRE(u_delta_x.isApprox(Eigen::Map<Eigen::VectorXd>(
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
Expand Down
13 changes: 8 additions & 5 deletions test/sampling.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@

#include <catch2/catch_test_macros.hpp>

#include <catch2/catch_approx.hpp> // for Approx
#include <sparseir/sparseir-header-only.hpp>
#include <xprec/ddouble-header-only.hpp>

Expand All @@ -12,6 +13,7 @@

using namespace sparseir;
using namespace std;
using Catch::Approx;

using ComplexF64 = std::complex<double>;

Expand Down Expand Up @@ -39,11 +41,11 @@ TEST_CASE("TauSampling test", "[sampling]") {
};
Eigen::MatrixXd mat_ref = Eigen::Map<Eigen::MatrixXd>(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") {
Expand Down Expand Up @@ -157,7 +159,7 @@ TEST_CASE("complex_test") {
}

auto tau_sampling = TauSampling<Bosonic>(basis);

/*
{
int dim = 0;
Eigen::Tensor<ComplexF64, 2> gl = movedim(originalgl, 0, dim);
Expand All @@ -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;
Expand Down
Loading

0 comments on commit e16ad7e

Please sign in to comment.