From 0f1e02a917ccfe08c0f27d27cb364cb3ccaae703 Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Fri, 22 Nov 2024 14:09:35 +0900 Subject: [PATCH 01/15] [WIP] Porting Python to C++ --- include/sparseir/kernel.hpp | 666 ++++++++++++++++++++++ include/sparseir/sparseir-header-only.hpp | 1 + test/CMakeLists.txt | 1 + test/kernel.cxx | 142 +++++ 4 files changed, 810 insertions(+) create mode 100644 include/sparseir/kernel.hpp create mode 100644 test/kernel.cxx diff --git a/include/sparseir/kernel.hpp b/include/sparseir/kernel.hpp new file mode 100644 index 0000000..42126c4 --- /dev/null +++ b/include/sparseir/kernel.hpp @@ -0,0 +1,666 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +namespace sparseir +{ + // Forward declaration of ReducedKernel + class ReducedKernel; + + /** + * @brief Abstract base class for an integral kernel K(x, y). + * + * AbstractKernel represents a real binary function K(x, y) used in a Fredholm + * integral equation of the first kind: + * + * u(x) = ∫ K(x, y) v(y) dy + * + * where x ∈ [xmin, xmax] and y ∈ [ymin, ymax]. + * + * In general, the kernel is applied to a scaled spectral function ρ'(y) as: + * + * ∫ K(x, y) ρ'(y) dy, + * + * where ρ'(y) = w(y) ρ(y). + */ + class AbstractKernel + { + public: + /** + * @brief Evaluate kernel at point (x, y). + * + * For given x, y, return the value of K(x, y). The parameters x_plus and + * x_minus, if given, shall contain the values of x - xmin and xmax - x, + * respectively. This is useful if either difference is to be formed and + * cancellation is expected. + * + * @param x The x value. + * @param y The y value. + * @param x_plus Optional. x - xmin. + * @param x_minus Optional. xmax - x. + * @return The value of K(x, y). + */ + virtual double operator()(double x, double y, double x_plus = std::numeric_limits::quiet_NaN(), + double x_minus = std::numeric_limits::quiet_NaN()) const = 0; + + /** + * @brief Return symmetrized kernel K(x, y) + sign * K(x, -y). + * + * Should be overridden by derived classes if they support symmetrization. + * + * @param sign The sign (+1 or -1). + * @return A shared pointer to the symmetrized kernel. + */ + virtual std::shared_ptr get_symmetrized(int sign) const + { + throw std::runtime_error("get_symmetrized not implemented in base class"); + } + + /** + * @brief Return tuple (xmin, xmax) delimiting the range of allowed x values. + * + * @return A pair containing xmin and xmax. + */ + virtual std::pair xrange() const + { + return std::make_pair(-1.0, 1.0); + } + + /** + * @brief Return tuple (ymin, ymax) delimiting the range of allowed y values. + * + * @return A pair containing ymin and ymax. + */ + virtual std::pair yrange() const + { + return std::make_pair(-1.0, 1.0); + } + + /** + * @brief Check if the kernel is centrosymmetric. + * + * Returns true if and only if K(x, y) == K(-x, -y) for all values of x and y. + * This allows the kernel to be block-diagonalized, speeding up the singular + * value expansion by a factor of 4. Defaults to false. + * + * @return True if the kernel is centrosymmetric, false otherwise. + */ + virtual bool is_centrosymmetric() const + { + return false; + } + + /** + * @brief Power with which the y coordinate scales. + * + * @return The power with which y scales. + */ + virtual int ypower() const + { + return 0; + } + + /** + * @brief Convergence radius of the Matsubara basis asymptotic model. + * + * For improved relative numerical accuracy, the IR basis functions on the + * Matsubara axis can be evaluated from an asymptotic expression for + * abs(n) > conv_radius. If conv_radius is infinity, then the asymptotics are + * unused (the default). + * + * @return The convergence radius. + */ + virtual double conv_radius() const + { + return std::numeric_limits::infinity(); + } + + /** + * @brief Return the weight function for given statistics. + * + * @param statistics 'F' for fermions or 'B' for bosons. + * @return A function representing the weight function w(y). + */ + virtual std::function weight_func(char statistics) const + { + if (statistics != 'F' && statistics != 'B') + { + throw std::invalid_argument("statistics must be 'F' for fermions or 'B' for bosons"); + } + return [](double /*x*/) + { return 1.0; }; + } + + virtual ~AbstractKernel() {} + }; + + /** + * @brief Fermionic/bosonic analytical continuation kernel. + * + * In dimensionless variables x = 2τ/β - 1, y = βω/Λ, the integral kernel is a function on [-1, 1] × [-1, 1]: + * + * K(x, y) = exp(-Λ y (x + 1) / 2) / (1 + exp(-Λ y)) + * + * LogisticKernel is a fermionic analytic continuation kernel. + * Nevertheless, one can model the τ dependence of a bosonic correlation function as follows: + * + * ∫ [exp(-Λ y (x + 1) / 2) / (1 - exp(-Λ y))] ρ(y) dy = ∫ K(x, y) ρ'(y) dy, + * + * with ρ'(y) = w(y) ρ(y), where the weight function is given by w(y) = 1 / tanh(Λ y / 2). + */ + class LogisticKernel : public AbstractKernel + { + public: + double lambda_; ///< The kernel cutoff Λ. + + /** + * @brief Constructor for LogisticKernel. + * + * @param lambda The kernel cutoff Λ. + */ + explicit LogisticKernel(double lambda) : lambda_(lambda) + { + if (lambda_ < 0) + { + throw std::domain_error("Kernel cutoff Λ must be non-negative"); + } + } + + /** + * @brief Evaluate the kernel at point (x, y). + * + * @param x The x value. + * @param y The y value. + * @param x_plus Optional. x - xmin. + * @param x_minus Optional. xmax - x. + * @return The value of K(x, y). + */ + double operator()(double x, double y, double x_plus = std::numeric_limits::quiet_NaN(), + double x_minus = std::numeric_limits::quiet_NaN()) const override + { + // Check that x and y are within the valid ranges + std::pair x_range = xrange(); + double xmin = x_range.first; + double xmax = x_range.second; + if (x < xmin || x > xmax) + { + throw std::out_of_range("x value not in range [-1, 1]"); + } + std::pair y_range = yrange(); + double ymin = y_range.first; + double ymax = y_range.second; + if (y < ymin || y > ymax) + { + throw std::out_of_range("y value not in range [-1, 1]"); + } + + std::tuple uv_values = compute_uv(x, y, x_plus, x_minus); + double u_plus = std::get<0>(uv_values); + double u_minus = std::get<1>(uv_values); + double v = std::get<2>(uv_values); + + return compute(u_plus, u_minus, v); + } + + /** + * @brief Check if the kernel is centrosymmetric. + * + * @return True, since LogisticKernel is centrosymmetric. + */ + bool is_centrosymmetric() const override + { + return true; + } + + /** + * @brief Convergence radius of the Matsubara basis asymptotic model. + * + * For LogisticKernel, conv_radius = 40 * Λ. + * + * @return The convergence radius. + */ + double conv_radius() const override + { + return 40.0 * lambda_; + } + + /** + * @brief Return the weight function for given statistics. + * + * @param statistics 'F' for fermions or 'B' for bosons. + * @return A function representing the weight function w(y). + */ + std::function weight_func(char statistics) const override + { + if (statistics == 'F') + { + // Fermionic weight function: w(y) == 1 + return [](double /*y*/) + { return 1.0; }; + } + else if (statistics == 'B') + { + // Bosonic weight function: w(y) == 1 / tanh(Λ*y/2) + double lambda = lambda_; + return [lambda](double y) + { + return 1.0 / std::tanh(0.5 * lambda * y); + }; + } + else + { + throw std::invalid_argument("statistics must be 'F' for fermions or 'B' for bosons"); + } + } + + private: + /** + * @brief Compute the variables u_plus, u_minus, v. + * + * @param x The x value. + * @param y The y value. + * @param x_plus x - xmin. + * @param x_minus xmax - x. + * @return A tuple containing u_plus, u_minus, and v. + */ + std::tuple compute_uv(double x, double y, double x_plus, double x_minus) const + { + // Compute u_plus, u_minus, v + if (std::isnan(x_plus)) + { + x_plus = 1.0 + x; + } + if (std::isnan(x_minus)) + { + x_minus = 1.0 - x; + } + double u_plus = 0.5 * x_plus; + double u_minus = 0.5 * x_minus; + double v = lambda_ * y; + return std::make_tuple(u_plus, u_minus, v); + } + + /** + * @brief Compute the kernel value using u_plus, u_minus, and v. + * + * @param u_plus Computed u_plus. + * @param u_minus Computed u_minus. + * @param v Computed v. + * @return The value of K(x, y). + */ + double compute(double u_plus, double u_minus, double v) const + { + double abs_v = std::abs(v); + + double numerator; + double denominator; + + if (v >= 0) + { + numerator = std::exp(-u_plus * abs_v); + } + else + { + numerator = std::exp(-u_minus * abs_v); + } + + denominator = 1.0 + std::exp(-abs_v); + + return numerator / denominator; + } + }; + + /** + * @brief Regularized bosonic analytical continuation kernel. + * + * In dimensionless variables x = 2τ/β - 1, y = βω/Λ, the integral kernel is a function on [-1, 1] × [-1, 1]: + * + * K(x, y) = y * exp(-Λ y (x + 1) / 2) / (exp(-Λ y) - 1) + * + * Care has to be taken in evaluating this expression around y = 0. + */ + class RegularizedBoseKernel : public AbstractKernel + { + public: + double lambda_; ///< The kernel cutoff Λ. + + /** + * @brief Constructor for RegularizedBoseKernel. + * + * @param lambda The kernel cutoff Λ. + */ + explicit RegularizedBoseKernel(double lambda) : lambda_(lambda) + { + if (lambda_ < 0) + { + throw std::domain_error("Kernel cutoff Λ must be non-negative"); + } + } + + /** + * @brief Evaluate the kernel at point (x, y). + * + * @param x The x value. + * @param y The y value. + * @param x_plus Optional. x - xmin. + * @param x_minus Optional. xmax - x. + * @return The value of K(x, y). + */ + double operator()(double x, double y, double x_plus = std::numeric_limits::quiet_NaN(), + double x_minus = std::numeric_limits::quiet_NaN()) const override + { + // Check that x and y are within the valid ranges + std::pair xrange_values = xrange(); + double xmin = xrange_values.first; + double xmax = xrange_values.second; + if (x < xmin || x > xmax) + { + throw std::out_of_range("x value not in range [-1, 1]"); + } + std::pair yrange_values = yrange(); + double ymin = yrange_values.first; + double ymax = yrange_values.second; + if (y < ymin || y > ymax) + { + throw std::out_of_range("y value not in range [-1, 1]"); + } + + std::tuple uv_values = compute_uv(x, y, x_plus, x_minus); + double u_plus = std::get<0>(uv_values); + double u_minus = std::get<1>(uv_values); + double v = std::get<2>(uv_values); + + return compute(u_plus, u_minus, v); + } + + /** + * @brief Check if the kernel is centrosymmetric. + * + * @return True, since RegularizedBoseKernel is centrosymmetric. + */ + bool is_centrosymmetric() const override + { + return true; + } + + /** + * @brief Returns the power with which y scales. + * + * For RegularizedBoseKernel, ypower = 1. + * + * @return The power with which y scales. + */ + int ypower() const override + { + return 1; + } + + /** + * @brief Convergence radius of the Matsubara basis asymptotic model. + * + * For RegularizedBoseKernel, conv_radius = 40 * Λ. + * + * @return The convergence radius. + */ + double conv_radius() const override + { + return 40.0 * lambda_; + } + + /** + * @brief Return the weight function for given statistics. + * + * @param statistics 'F' for fermions or 'B' for bosons. + * @return A function representing the weight function w(y). + */ + std::function weight_func(char statistics) const override + { + if (statistics == 'F') + { + throw std::runtime_error("Kernel is designed for bosonic functions"); + } + else if (statistics == 'B') + { + // Bosonic weight function: w(y) == 1 / y + return [](double y) + { return 1.0 / y; }; + } + else + { + throw std::invalid_argument("statistics must be 'F' for fermions or 'B' for bosons"); + } + } + + private: + /** + * @brief Compute the variables u_plus, u_minus, v. + * + * @param x The x value. + * @param y The y value. + * @param x_plus x - xmin. + * @param x_minus xmax - x. + * @return A tuple containing u_plus, u_minus, and v. + */ + std::tuple compute_uv(double x, double y, double x_plus, double x_minus) const + { + // Compute u_plus, u_minus, v + if (std::isnan(x_plus)) + { + x_plus = 1.0 + x; + } + if (std::isnan(x_minus)) + { + x_minus = 1.0 - x; + } + double u_plus = 0.5 * x_plus; + double u_minus = 0.5 * x_minus; + double v = lambda_ * y; + return std::make_tuple(u_plus, u_minus, v); + } + + /** + * @brief Compute the kernel value using u_plus, u_minus, and v. + * + * @param u_plus Computed u_plus. + * @param u_minus Computed u_minus. + * @param v Computed v. + * @return The value of K(x, y). + */ + double compute(double u_plus, double u_minus, double v) const + { + double absv = std::abs(v); + + double numerator; + double denominator; + + if (v >= 0) + { + numerator = std::exp(-u_plus * absv); + } + else + { + numerator = std::exp(-u_minus * absv); + } + + // Handle small values of absv to avoid division by zero + double value; + + if (absv > 1e-200) + { + denominator = std::expm1(-absv); // exp(-absv) - 1 + value = -1.0 / lambda_ * numerator * (absv / denominator); + } + else + { + // Limit as absv -> 0 + value = -1.0 / lambda_ * numerator * (-1.0); + } + + return value; + } + }; + + /** + * @brief Restriction of centrosymmetric kernel to positive interval. + * + * For a kernel K on [-1, 1] × [-1, 1] that is centrosymmetric, i.e., + * K(x, y) = K(-x, -y), it is straightforward to show that the left/right + * singular vectors can be chosen as either odd or even functions. + * + * Consequently, they are singular functions of a reduced kernel K_red on + * [0, 1] × [0, 1] that is given as either: + * + * K_red(x, y) = K(x, y) ± K(x, -y) + * + * This kernel is what this class represents. The full singular functions can be + * reconstructed by (anti-)symmetrically continuing them to the negative axis. + */ + class ReducedKernel : public AbstractKernel + { + public: + std::shared_ptr inner_kernel_; ///< The inner kernel K. + int sign_; ///< The sign (+1 or -1). + + /** + * @brief Constructor for ReducedKernel. + * + * @param inner_kernel The inner kernel K. + * @param sign The sign (+1 or -1). Must satisfy abs(sign) == 1. + */ + ReducedKernel(std::shared_ptr inner_kernel, int sign) + : inner_kernel_(std::move(inner_kernel)), sign_(sign) + { + if (!inner_kernel_->is_centrosymmetric()) + { + throw std::invalid_argument("Inner kernel must be centrosymmetric"); + } + if (sign != 1 && sign != -1) + { + throw std::invalid_argument("sign must be -1 or 1"); + } + } + + /** + * @brief Evaluate the reduced kernel at point (x, y). + * + * @param x The x value. + * @param y The y value. + * @param x_plus Optional. x - xmin. + * @param x_minus Optional. xmax - x. + * @return The value of K_red(x, y). + */ + double operator()(double x, double y, double x_plus = std::numeric_limits::quiet_NaN(), + double x_minus = std::numeric_limits::quiet_NaN()) const override + { + return call_reduced(x, y, x_plus, x_minus); + } + + /** + * @brief Return tuple (xmin, xmax) delimiting the range of allowed x values. + * + * For ReducedKernel, xrange is modified to [0, xmax_inner]. + * + * @return A pair containing xmin and xmax. + */ + std::pair xrange() const override + { + auto range = inner_kernel_->xrange(); + return std::make_pair(0.0, range.second); + } + + /** + * @brief Return tuple (ymin, ymax) delimiting the range of allowed y values. + * + * For ReducedKernel, yrange is modified to [0, ymax_inner]. + * + * @return A pair containing ymin and ymax. + */ + std::pair yrange() const override + { + auto range = inner_kernel_->yrange(); + return std::make_pair(0.0, range.second); + } + + /** + * @brief Check if the kernel is centrosymmetric. + * + * @return False, since ReducedKernel cannot be symmetrized further. + */ + bool is_centrosymmetric() const override + { + return false; + } + + /** + * @brief Attempting to symmetrize a ReducedKernel will result in an error. + * + * @param sign The sign (+1 or -1). + * @throws std::runtime_error Cannot symmetrize twice. + */ + std::shared_ptr get_symmetrized(int /*sign*/) const override + { + throw std::runtime_error("Cannot symmetrize twice"); + } + + /** + * @brief Returns the power with which y scales. + * + * @return The ypower of the inner kernel. + */ + int ypower() const override + { + return inner_kernel_->ypower(); + } + + /** + * @brief Convergence radius of the Matsubara basis asymptotic model. + * + * @return The convergence radius of the inner kernel. + */ + double conv_radius() const override + { + return inner_kernel_->conv_radius(); + } + + private: + /** + * @brief Evaluate the reduced kernel. + * + * @param x The x value. + * @param y The y value. + * @param x_plus x - xmin. + * @param x_minus xmax - x. + * @return The value of K_red(x, y). + */ + double call_reduced(double x, double y, double x_plus, double x_minus) const + { + // The reduced kernel is defined only over the interval [0, 1], which + // means we must add one to get the x_plus for the inner kernels. + if (std::isnan(x_plus)) + { + x_plus = 1.0 + x; + } + // x_minus remains the same + + // Evaluate inner kernel at (x, y) and (x, -y) + double K_plus = inner_kernel_->operator()(x, y, x_plus, x_minus); + double K_minus = inner_kernel_->operator()(x, -y, x_plus, x_minus); + + if (sign_ == 1) + { + return K_plus + K_minus; + } + else + { + return K_plus - K_minus; + } + } + }; + +} // namespace sparseir \ No newline at end of file diff --git a/include/sparseir/sparseir-header-only.hpp b/include/sparseir/sparseir-header-only.hpp index 8077a34..86d5095 100644 --- a/include/sparseir/sparseir-header-only.hpp +++ b/include/sparseir/sparseir-header-only.hpp @@ -7,3 +7,4 @@ #include "./gauss.hpp" #include "./freq.hpp" #include "./poly.hpp" +#include "./kernel.hpp" \ No newline at end of file diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 1ac9408..1282b4d 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -19,6 +19,7 @@ add_executable(libsparseirtests gauss.cxx freq.cxx poly.cxx + kernel.cxx ) target_link_libraries(libsparseirtests PRIVATE Catch2::Catch2WithMain) diff --git a/test/kernel.cxx b/test/kernel.cxx new file mode 100644 index 0000000..19b81a9 --- /dev/null +++ b/test/kernel.cxx @@ -0,0 +1,142 @@ +#include +#include +#include +#include +#include +#include + +#include + +#include + +using xprec::DDouble; + +TEST_CASE("kernel.hpp") +{ + // Use a vector of shared pointers to AbstractKernel objects + std::vector> kernels = { + std::make_shared(9.0), + std::make_shared(8.0), + std::make_shared(120000.0), + std::make_shared(127500.0), + // Obtain symmetrized kernels + //sparseir::get_symmetrized(std::make_shared(40000.0), -1), + //sparseir::get_symmetrized(std::make_shared(35000.0), -1) + }; + + for (const auto& K_ptr : kernels) + { + const auto& K = *K_ptr; // Dereference the shared_ptr to get the kernel + + using T = float; + using T_x = double; + + // Convert Rule to type T + sparseir::Rule _rule = sparseir::legendre(10); + sparseir::Rule rule = sparseir::convert(_rule); + + // Obtain SVE hints for the kernel + // TODO: Implement sve_hints + // auto hints = sparseir::sve_hints(K, std::numeric_limits::epsilon()); + + // Generate piecewise Gaussian quadrature rules for x and y + // TODO: Implement piecewise and segments_x and segments_y + // auto gauss_x = sparseir::piecewise(rule, sparseir::segments_x(hints)); + // auto gauss_y = sparseir::piecewise(rule, sparseir::segments_y(hints)); + + T epsilon = std::numeric_limits::epsilon(); + T tiny = std::numeric_limits::min() / epsilon; + + // Compute the matrix from Gaussian quadrature + // TODO: Implement matrix_from_gauss + // Eigen::MatrixX result = sparseir::matrix_from_gauss(K, gauss_x, gauss_y); + + // Convert gauss_x and gauss_y to higher precision T_x + // TODO: Implement convert + // auto gauss_x_Tx = sparseir::convert(gauss_x); + // auto gauss_y_Tx = sparseir::convert(gauss_y); + + // Compute the matrix in higher precision + // TODO: Implement matrix_from_gauss + // Eigen::MatrixX result_x = sparseir::matrix_from_gauss(K, gauss_x_Tx, gauss_y_Tx); + + // T_x magn = result_x.cwiseAbs().maxCoeff(); + + // Check that the difference is within tolerance + // REQUIRE((result.template cast() - result_x).cwiseAbs().maxCoeff() <= 2 * magn * epsilon); + + // Eigen::Array reldiff = (result.cwiseAbs().array() < tiny) + // .select(T(1.0), result.array() / result_x.template cast().array()); + + // REQUIRE((reldiff - T(1.0)).cwiseAbs().maxCoeff() <= 100 * epsilon); + } + + // Second test set: "singularity with Λ = Λ" + std::vector lambdas = {10.0, 42.0, 10000.0}; + std::mt19937 rng; + std::uniform_real_distribution dist(-1.0, 1.0); + + for (double Lambda : lambdas) + { + std::vector xs(10); + for (auto& x : xs) + { + x = dist(rng); + } + + auto K_ptr = std::make_shared(Lambda); + const auto& K = *K_ptr; + const double tol = 1e-6; + + for (double x : xs) + { + REQUIRE(std::abs(K(x, 0.0) - 1.0 / Lambda) <= tol); + } + } + + // Third test set: "unit tests" + { + auto K_ptr = std::make_shared(42.0); + const auto& K = *K_ptr; + auto K_symm_ptr = sparseir::get_symmetrized(K_ptr, 1); + const auto& K_symm = *K_symm_ptr; + + REQUIRE(!sparseir::iscentrosymmetric(K_symm)); + + REQUIRE_THROWS_AS(sparseir::get_symmetrized(K_symm_ptr, -1), std::runtime_error); + + double expected_value = 1.0 / std::tanh(0.5 * 42.0 * 1e-8); + double computed_value = sparseir::weight_func(K, sparseir::Bosonic())(1e-8); + REQUIRE(computed_value == Approx(expected_value)); + + REQUIRE(sparseir::weight_func(K, sparseir::Fermionic())(482) == Approx(1.0)); + + REQUIRE(sparseir::weight_func(K_symm, sparseir::Bosonic())(1e-3) == Approx(1.0)); + + REQUIRE(sparseir::weight_func(K_symm, sparseir::Fermionic())(482) == Approx(1.0)); + + auto K99_ptr = std::make_shared(99.0); + const auto& K99 = *K99_ptr; + auto hints = sparseir::sve_hints(K99, 1e-6); + + REQUIRE(sparseir::nsvals(hints) == 56); + + REQUIRE(sparseir::ngauss(hints) == 10); + + REQUIRE(sparseir::ypower(K99) == 1); + + REQUIRE(sparseir::ypower(*sparseir::get_symmetrized(K99_ptr, -1)) == 1); + + REQUIRE(sparseir::ypower(*sparseir::get_symmetrized(K99_ptr, 1)) == 1); + + REQUIRE(sparseir::conv_radius(K99) == Approx(40.0 * 99.0)); + + REQUIRE(sparseir::conv_radius(*sparseir::get_symmetrized(K99_ptr, -1)) == Approx(40.0 * 99.0)); + + REQUIRE(sparseir::conv_radius(*sparseir::get_symmetrized(K99_ptr, 1)) == Approx(40.0 * 99.0)); + + REQUIRE(sparseir::weight_func(K99, sparseir::Bosonic())(482) == Approx(1.0 / 482.0)); + + REQUIRE_THROWS_AS(sparseir::weight_func(K99, sparseir::Fermionic())(1.0), std::runtime_error); + } +} \ No newline at end of file From 45388e779b7b097eaede62254bd303537b058e18 Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Fri, 22 Nov 2024 14:24:46 +0900 Subject: [PATCH 02/15] Add conversion function --- include/sparseir/gauss.hpp | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/include/sparseir/gauss.hpp b/include/sparseir/gauss.hpp index 4fe30ea..432a83e 100644 --- a/include/sparseir/gauss.hpp +++ b/include/sparseir/gauss.hpp @@ -192,4 +192,17 @@ Matrix legendre_collocation(const Rule& rule, int n = -1 return res; } +template +Rule convert(const Rule &rule) +{ + std::vector x(rule.x.begin(), rule.x.end()); + std::vector w(rule.w.begin(), rule.w.end()); + TargetType a = static_cast(rule.a); + TargetType b = static_cast(rule.b); + std::vector x_forward(rule.x_forward.begin(), rule.x_forward.end()); + std::vector x_backward(rule.x_backward.begin(), rule.x_backward.end()); + + return Rule(x, w, a, b, x_forward, x_backward); +} + } // namespace sparseir \ No newline at end of file From 235c7dc0c04d7906df64ca675adb0d9a095849e3 Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Fri, 22 Nov 2024 14:40:40 +0900 Subject: [PATCH 03/15] Update test for kernel --- test/kernel.cxx | 122 +++++++++++++++--------------------------------- 1 file changed, 37 insertions(+), 85 deletions(-) diff --git a/test/kernel.cxx b/test/kernel.cxx index 19b81a9..789a719 100644 --- a/test/kernel.cxx +++ b/test/kernel.cxx @@ -1,142 +1,94 @@ +#include #include -#include #include #include #include +#include #include -#include - +// Include your sparseir headers #include using xprec::DDouble; -TEST_CASE("kernel.hpp") +TEST_CASE("Kernel Accuracy Test") { - // Use a vector of shared pointers to AbstractKernel objects + using T = float; + using T_x = double; + + // List of kernels to test std::vector> kernels = { std::make_shared(9.0), std::make_shared(8.0), std::make_shared(120000.0), std::make_shared(127500.0), - // Obtain symmetrized kernels - //sparseir::get_symmetrized(std::make_shared(40000.0), -1), - //sparseir::get_symmetrized(std::make_shared(35000.0), -1) + // Symmetrized kernels + std::make_shared(40000.0)->get_symmetrized(-1), + std::make_shared(35000.0)->get_symmetrized(-1), }; - for (const auto& K_ptr : kernels) + for (const auto &K_ptr : kernels) { - const auto& K = *K_ptr; // Dereference the shared_ptr to get the kernel - - using T = float; - using T_x = double; + const auto &K = *K_ptr; // Convert Rule to type T sparseir::Rule _rule = sparseir::legendre(10); sparseir::Rule rule = sparseir::convert(_rule); // Obtain SVE hints for the kernel - // TODO: Implement sve_hints - // auto hints = sparseir::sve_hints(K, std::numeric_limits::epsilon()); + auto hints = K.sve_hints(2.2e-16); // Generate piecewise Gaussian quadrature rules for x and y - // TODO: Implement piecewise and segments_x and segments_y - // auto gauss_x = sparseir::piecewise(rule, sparseir::segments_x(hints)); - // auto gauss_y = sparseir::piecewise(rule, sparseir::segments_y(hints)); + auto gauss_x = sparseir::piecewise(rule, hints.template segments_x()); + auto gauss_y = sparseir::piecewise(rule, hints.template segments_y()); T epsilon = std::numeric_limits::epsilon(); T tiny = std::numeric_limits::min() / epsilon; // Compute the matrix from Gaussian quadrature - // TODO: Implement matrix_from_gauss - // Eigen::MatrixX result = sparseir::matrix_from_gauss(K, gauss_x, gauss_y); + Eigen::Matrix result = sparseir::matrix_from_gauss(K, gauss_x, gauss_y); // Convert gauss_x and gauss_y to higher precision T_x - // TODO: Implement convert - // auto gauss_x_Tx = sparseir::convert(gauss_x); - // auto gauss_y_Tx = sparseir::convert(gauss_y); + auto gauss_x_Tx = sparseir::convert(gauss_x); + auto gauss_y_Tx = sparseir::convert(gauss_y); // Compute the matrix in higher precision - // TODO: Implement matrix_from_gauss - // Eigen::MatrixX result_x = sparseir::matrix_from_gauss(K, gauss_x_Tx, gauss_y_Tx); + Eigen::Matrix result_x = sparseir::matrix_from_gauss(K, gauss_x_Tx, gauss_y_Tx); - // T_x magn = result_x.cwiseAbs().maxCoeff(); + T_x magn = result_x.cwiseAbs().maxCoeff(); // Check that the difference is within tolerance - // REQUIRE((result.template cast() - result_x).cwiseAbs().maxCoeff() <= 2 * magn * epsilon); + REQUIRE((result.template cast() - result_x).cwiseAbs().maxCoeff() <= 2 * magn * epsilon); - // Eigen::Array reldiff = (result.cwiseAbs().array() < tiny) - // .select(T(1.0), result.array() / result_x.template cast().array()); + Eigen::Array reldiff = (result.cwiseAbs().array() < tiny) + .select(T(1.0), result.array() / result_x.template cast().array()); - // REQUIRE((reldiff - T(1.0)).cwiseAbs().maxCoeff() <= 100 * epsilon); + REQUIRE((reldiff - T(1.0)).cwiseAbs().maxCoeff() <= 100 * epsilon); } +} - // Second test set: "singularity with Λ = Λ" +TEST_CASE("Kernel Singularity Test") +{ std::vector lambdas = {10.0, 42.0, 10000.0}; std::mt19937 rng; std::uniform_real_distribution dist(-1.0, 1.0); - for (double Lambda : lambdas) + for (double lambda : lambdas) { - std::vector xs(10); - for (auto& x : xs) + // Generate random x values + std::vector x_values(1000); + for (double &x : x_values) { x = dist(rng); } - auto K_ptr = std::make_shared(Lambda); - const auto& K = *K_ptr; - const double tol = 1e-6; + sparseir::RegularizedBoseKernel K(lambda); - for (double x : xs) + for (double x : x_values) { - REQUIRE(std::abs(K(x, 0.0) - 1.0 / Lambda) <= tol); + double expected = 1.0 / lambda; + double computed = K(x, 0.0); + REQUIRE(isApprox(computed, expected)); } } - - // Third test set: "unit tests" - { - auto K_ptr = std::make_shared(42.0); - const auto& K = *K_ptr; - auto K_symm_ptr = sparseir::get_symmetrized(K_ptr, 1); - const auto& K_symm = *K_symm_ptr; - - REQUIRE(!sparseir::iscentrosymmetric(K_symm)); - - REQUIRE_THROWS_AS(sparseir::get_symmetrized(K_symm_ptr, -1), std::runtime_error); - - double expected_value = 1.0 / std::tanh(0.5 * 42.0 * 1e-8); - double computed_value = sparseir::weight_func(K, sparseir::Bosonic())(1e-8); - REQUIRE(computed_value == Approx(expected_value)); - - REQUIRE(sparseir::weight_func(K, sparseir::Fermionic())(482) == Approx(1.0)); - - REQUIRE(sparseir::weight_func(K_symm, sparseir::Bosonic())(1e-3) == Approx(1.0)); - - REQUIRE(sparseir::weight_func(K_symm, sparseir::Fermionic())(482) == Approx(1.0)); - - auto K99_ptr = std::make_shared(99.0); - const auto& K99 = *K99_ptr; - auto hints = sparseir::sve_hints(K99, 1e-6); - - REQUIRE(sparseir::nsvals(hints) == 56); - - REQUIRE(sparseir::ngauss(hints) == 10); - - REQUIRE(sparseir::ypower(K99) == 1); - - REQUIRE(sparseir::ypower(*sparseir::get_symmetrized(K99_ptr, -1)) == 1); - - REQUIRE(sparseir::ypower(*sparseir::get_symmetrized(K99_ptr, 1)) == 1); - - REQUIRE(sparseir::conv_radius(K99) == Approx(40.0 * 99.0)); - - REQUIRE(sparseir::conv_radius(*sparseir::get_symmetrized(K99_ptr, -1)) == Approx(40.0 * 99.0)); - - REQUIRE(sparseir::conv_radius(*sparseir::get_symmetrized(K99_ptr, 1)) == Approx(40.0 * 99.0)); - - REQUIRE(sparseir::weight_func(K99, sparseir::Bosonic())(482) == Approx(1.0 / 482.0)); - - REQUIRE_THROWS_AS(sparseir::weight_func(K99, sparseir::Fermionic())(1.0), std::runtime_error); - } } \ No newline at end of file From 48c1358a0b9d6d34f74dc61b80bf324c15235611 Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Fri, 22 Nov 2024 15:01:29 +0900 Subject: [PATCH 04/15] Fix constructor --- include/sparseir/gauss.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/sparseir/gauss.hpp b/include/sparseir/gauss.hpp index 432a83e..a6576d4 100644 --- a/include/sparseir/gauss.hpp +++ b/include/sparseir/gauss.hpp @@ -202,7 +202,7 @@ Rule convert(const Rule &rule) std::vector x_forward(rule.x_forward.begin(), rule.x_forward.end()); std::vector x_backward(rule.x_backward.begin(), rule.x_backward.end()); - return Rule(x, w, a, b, x_forward, x_backward); + return Rule(x, w, x_forward, x_backward, a, b); } } // namespace sparseir \ No newline at end of file From 01b1f3d2600ab13a19e6fc8d4d5217e08dc66cc3 Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Fri, 22 Nov 2024 15:01:37 +0900 Subject: [PATCH 05/15] Fix poly --- test/poly.cxx | 53 +++++++++++++++++++-------------------------------- 1 file changed, 20 insertions(+), 33 deletions(-) diff --git a/test/poly.cxx b/test/poly.cxx index bf876ea..6a987f6 100644 --- a/test/poly.cxx +++ b/test/poly.cxx @@ -21,16 +21,6 @@ #include #include - -// Helper function for approximate equality -bool isApprox(const Eigen::MatrixXd& a, const Eigen::MatrixXd& b, double tol = 1e-12) { - return ((a - b).array().abs() < tol).all(); -} - -bool isApprox(const double a, const double b, double tol = 1e-12) { - return (std::abs(a - b) < tol); -} - TEST_CASE("StableRNG") { // Initialize data directly with the given values Eigen::MatrixXd data(3, 3); @@ -42,9 +32,6 @@ TEST_CASE("StableRNG") { knots << 0.507134318967235, 0.5766150365607372, 0.7126662232433161, 0.7357313003784003; - // Check that data matches expected values - REQUIRE(isApprox(data, data)); - // Check that knots are sorted REQUIRE(std::is_sorted(knots.data(), knots.data() + knots.size())); @@ -56,7 +43,7 @@ TEST_CASE("StableRNG") { 0.13124858727993338, 0.2193663343416914, 0.7756615110113394; // Check that ddata matches expected values - REQUIRE(isApprox(ddata, ddata)); + REQUIRE(ddata.isApprox(ddata)); } TEST_CASE("sparseir::PiecewiseLegendrePoly(data::Matrix, knots::Vector, l::Int)") { @@ -76,10 +63,10 @@ TEST_CASE("sparseir::PiecewiseLegendrePoly(data::Matrix, knots::Vector, l::Int)" sparseir::PiecewiseLegendrePoly pwlp(data, knots, l); // Test that the object is initialized correctly - REQUIRE(isApprox(pwlp.data, data)); - REQUIRE(isApprox(pwlp.xmin, knots[0])); - REQUIRE(isApprox(pwlp.xmax, knots[knots.size() - 1])); - REQUIRE(isApprox(pwlp.knots, knots)); + REQUIRE(pwlp.data.isApprox(data)); + REQUIRE(pwlp.xmin == knots[0]); + REQUIRE(pwlp.xmax == knots[knots.size() - 1]); + REQUIRE(pwlp.knots.isApprox(knots)); REQUIRE(pwlp.polyorder == data.rows()); REQUIRE(pwlp.symm == 0); } @@ -110,19 +97,19 @@ TEST_CASE("PiecewiseLegendrePoly(data, p::PiecewiseLegendrePoly; symm=symm(p))") sparseir::PiecewiseLegendrePoly ddata_pwlp(ddata, pwlp, randsymm); // Test that ddata_pwlp is initialized correctly - REQUIRE(isApprox(ddata_pwlp.data, ddata)); + REQUIRE(ddata_pwlp.data.isApprox(ddata)); REQUIRE(ddata_pwlp.symm == randsymm); // Check that other fields match between pwlp and ddata_pwlp REQUIRE(pwlp.polyorder == ddata_pwlp.polyorder); REQUIRE(pwlp.xmin == ddata_pwlp.xmin); REQUIRE(pwlp.xmax == ddata_pwlp.xmax); - REQUIRE(isApprox(pwlp.knots, ddata_pwlp.knots)); - REQUIRE(isApprox(pwlp.delta_x, ddata_pwlp.delta_x)); + REQUIRE(pwlp.knots.isApprox(ddata_pwlp.knots)); + REQUIRE(pwlp.delta_x == ddata_pwlp.delta_x); REQUIRE(pwlp.l == ddata_pwlp.l); - REQUIRE(isApprox(pwlp.xm, ddata_pwlp.xm)); - REQUIRE(isApprox(pwlp.inv_xs, ddata_pwlp.inv_xs)); - REQUIRE(isApprox(pwlp.norms, ddata_pwlp.norms)); + REQUIRE(pwlp.xm.isApprox(ddata_pwlp.xm)); + REQUIRE(pwlp.inv_xs.isApprox(ddata_pwlp.inv_xs)); + REQUIRE(pwlp.norms.isApprox(ddata_pwlp.norms)); } TEST_CASE("sparseir::PiecewiseLegendrePolyVector") { @@ -213,10 +200,10 @@ TEST_CASE("sparseir::PiecewiseLegendrePolyVector") { // Test properties REQUIRE(polys.xmin() == pwlp1.xmin); REQUIRE(polys.xmax() == pwlp1.xmax); - REQUIRE(isApprox(polys.get_knots(), pwlp1.knots)); - REQUIRE(isApprox(polys.get_delta_x(), pwlp1.delta_x)); + REQUIRE(polys.get_knots().isApprox(pwlp1.knots)); + REQUIRE(polys.get_delta_x() == pwlp1.delta_x); REQUIRE(polys.get_polyorder() == pwlp1.polyorder); - REQUIRE(isApprox(polys.get_norms(), pwlp1.norms)); + REQUIRE(polys.get_norms().isApprox(pwlp1.norms)); // Test symm std::vector expected_symm = {pwlp1.symm, pwlp2.symm, pwlp3.symm}; @@ -274,19 +261,19 @@ TEST_CASE("Deriv") { sparseir::PiecewiseLegendrePoly deriv_pwlp = pwlp.deriv(); // Test that derivative data matches - REQUIRE(isApprox(deriv_pwlp.data, ddata)); + REQUIRE(deriv_pwlp.data.isApprox(ddata)); REQUIRE(deriv_pwlp.symm == 0); // Check that other fields match REQUIRE(pwlp.polyorder == deriv_pwlp.polyorder); REQUIRE(pwlp.xmin == deriv_pwlp.xmin); REQUIRE(pwlp.xmax == deriv_pwlp.xmax); - REQUIRE(isApprox(pwlp.knots, deriv_pwlp.knots)); - REQUIRE(isApprox(pwlp.delta_x, deriv_pwlp.delta_x)); + REQUIRE(pwlp.knots.isApprox(deriv_pwlp.knots)); + REQUIRE(pwlp.delta_x == deriv_pwlp.delta_x); REQUIRE(pwlp.l == deriv_pwlp.l); - REQUIRE(isApprox(pwlp.xm, deriv_pwlp.xm)); - REQUIRE(isApprox(pwlp.inv_xs, deriv_pwlp.inv_xs)); - REQUIRE(isApprox(pwlp.norms, deriv_pwlp.norms)); + REQUIRE(pwlp.xm.isApprox(deriv_pwlp.xm)); + REQUIRE(pwlp.inv_xs.isApprox(deriv_pwlp.inv_xs)); + REQUIRE(pwlp.norms.isApprox(deriv_pwlp.norms)); } TEST_CASE("Overlap") { From 87472b2de40263a39b8edde62dc153c5c83c2375 Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Fri, 22 Nov 2024 15:01:49 +0900 Subject: [PATCH 06/15] Fix conversion rule --- test/kernel.cxx | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/kernel.cxx b/test/kernel.cxx index 789a719..e38b833 100644 --- a/test/kernel.cxx +++ b/test/kernel.cxx @@ -32,8 +32,9 @@ TEST_CASE("Kernel Accuracy Test") const auto &K = *K_ptr; // Convert Rule to type T - sparseir::Rule _rule = sparseir::legendre(10); - sparseir::Rule rule = sparseir::convert(_rule); + sparseir::Rule ddouble_rule = sparseir::legendre(10); + sparseir::Rule double_rule = sparseir::convert(ddouble_rule); + sparseir::Rule rule = sparseir::convert(double_rule); // Obtain SVE hints for the kernel auto hints = K.sve_hints(2.2e-16); From 08d5336783d0ef1731c42773d0ba19ab4ebe2cee Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Fri, 22 Nov 2024 15:02:54 +0900 Subject: [PATCH 07/15] Use isApprox in Eigen --- test/kernel.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/kernel.cxx b/test/kernel.cxx index e38b833..a311a76 100644 --- a/test/kernel.cxx +++ b/test/kernel.cxx @@ -89,7 +89,7 @@ TEST_CASE("Kernel Singularity Test") { double expected = 1.0 / lambda; double computed = K(x, 0.0); - REQUIRE(isApprox(computed, expected)); + REQUIRE(Eigen::internal::isApprox(computed, expected)); } } } \ No newline at end of file From 48604935c059733afd17db2f724bbecf1de6046a Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Fri, 22 Nov 2024 15:38:53 +0900 Subject: [PATCH 08/15] [WIP] --- include/sparseir/kernel.hpp | 177 ++++++++++++++++++++++++++++++++++++ 1 file changed, 177 insertions(+) diff --git a/include/sparseir/kernel.hpp b/include/sparseir/kernel.hpp index 42126c4..421bbaf 100644 --- a/include/sparseir/kernel.hpp +++ b/include/sparseir/kernel.hpp @@ -8,10 +8,19 @@ #include #include +#pragma once + +#include +#include + namespace sparseir { // Forward declaration of ReducedKernel class ReducedKernel; + class SVEHints; + class SVEHintsLogistic; + class SVEHintsRegularizedBose; + /** * @brief Abstract base class for an integral kernel K(x, y). @@ -49,6 +58,8 @@ namespace sparseir virtual double operator()(double x, double y, double x_plus = std::numeric_limits::quiet_NaN(), double x_minus = std::numeric_limits::quiet_NaN()) const = 0; + virtual std::shared_ptr sve_hints(double epsilon) const = 0; + /** * @brief Return symmetrized kernel K(x, y) + sign * K(x, -y). * @@ -208,6 +219,14 @@ namespace sparseir return compute(u_plus, u_minus, v); } + /* + // Inside class LogisticKernel definition + std::shared_ptr sve_hints(double epsilon) const override + { + return std::make_shared(*this, epsilon); + } + */ + /** * @brief Check if the kernel is centrosymmetric. * @@ -379,6 +398,14 @@ namespace sparseir return compute(u_plus, u_minus, v); } + /* + // Inside class RegularizedBoseKernel definition + std::shared_ptr sve_hints(double epsilon) const override + { + return std::make_shared(*this, epsilon); + } + */ + /** * @brief Check if the kernel is centrosymmetric. * @@ -663,4 +690,154 @@ namespace sparseir } }; +} // namespace sparseir + +namespace sparseir +{ + + class SVEHints + { + public: + virtual ~SVEHints() = default; + + // Functions to compute segments for x and y + virtual std::vector segments_x() const = 0; + virtual std::vector segments_x_ddouble() const = 0; + + virtual std::vector segments_y() const = 0; + virtual std::vector segments_y_ddouble() const = 0; + + // Additional methods if needed + virtual int nsvals() const = 0; + virtual int ngauss() const = 0; + }; + + class SVEHintsLogistic : public SVEHints + { + public: + SVEHintsLogistic(const LogisticKernel &kernel, double epsilon) + : kernel_(kernel), epsilon_(epsilon) {} + + std::vector segments_x() const override; + std::vector segments_y() const override; + + std::vector segments_x_ddouble() const override; + std::vector segments_y_ddouble() const override; + + int nsvals() const override; + int ngauss() const override; + + private: + const LogisticKernel &kernel_; + double epsilon_; + }; + + class SVEHintsRegularizedBose : public SVEHints + { + public: + SVEHintsRegularizedBose(const RegularizedBoseKernel &kernel, double epsilon) + : kernel_(kernel), epsilon_(epsilon) {} + + template + std::vector segments_x() const; + + template + std::vector segments_y() const; + + int nsvals() const override; + int ngauss() const override; + + std::vector segments_x_ddouble() const override; + std::vector segments_y_ddouble() const override; + + private: + const RegularizedBoseKernel &kernel_; + double epsilon_; + }; + + template + std::vector SVEHintsLogistic::segments_x() const + { + T h = static_cast(1.0 / kernel_.lambda_); + std::vector segments; + + if (std::isinf(h) || std::isnan(h)) + { + segments.push_back(static_cast(0)); + } + else + { + segments = {static_cast(-1), -h, static_cast(0), h, static_cast(1)}; + } + + return segments; + } + + template + std::vector SVEHintsLogistic::segments_y() const + { + return {static_cast(-1), static_cast(0), static_cast(1)}; + } + + int SVEHintsLogistic::nsvals() const + { + // Implement as needed + return 0; // Placeholder + } + + int SVEHintsLogistic::ngauss() const + { + // Implement as needed + return 0; // Placeholder + } + + std::vector SVEHintsLogistic::segments_x() const + { + int nzeros = std::max(static_cast(std::round(15 * std::log10(kernel_.lambda_))), 1); + + std::vector temp(nzeros); + double coeff = 0.143; + for (int i = 0; i < nzeros; ++i) + { + temp[i] = coeff * static_cast(i); + } + + std::vector diffs(nzeros); + for (int i = 0; i < nzeros; ++i) + { + diffs[i] = 1.0 / std::cosh(temp[i]); + } + + std::vector zeros(nzeros); + zeros[0] = diffs[0]; + for (int i = 1; i < nzeros; ++i) + { + zeros[i] = zeros[i - 1] + diffs[i]; + } + + // Normalize zeros + double last_zero = zeros.back(); + for (double &z : zeros) + { + z /= last_zero; + } + + // Create segments: [-reverse(zeros); 0; zeros] + std::vector segments; + segments.reserve(2 * nzeros + 1); + for (auto it = zeros.rbegin(); it != zeros.rend(); ++it) + { + segments.push_back(-(*it)); + } + segments.push_back(0.0); + segments.insert(segments.end(), zeros.begin(), zeros.end()); + + return segments; + } + + std::vector SVEHintsLogistic::segments_y() const + { + return {-1.0, 0.0, 1.0}; + } + } // namespace sparseir \ No newline at end of file From e213c0931acd75aef6886e1a29197290d2f99c42 Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Fri, 22 Nov 2024 15:41:51 +0900 Subject: [PATCH 09/15] Fix --- include/sparseir/kernel.hpp | 27 ++------------------------- 1 file changed, 2 insertions(+), 25 deletions(-) diff --git a/include/sparseir/kernel.hpp b/include/sparseir/kernel.hpp index 421bbaf..6dd6a14 100644 --- a/include/sparseir/kernel.hpp +++ b/include/sparseir/kernel.hpp @@ -755,28 +755,9 @@ namespace sparseir double epsilon_; }; - template - std::vector SVEHintsLogistic::segments_x() const - { - T h = static_cast(1.0 / kernel_.lambda_); - std::vector segments; - - if (std::isinf(h) || std::isnan(h)) - { - segments.push_back(static_cast(0)); - } - else - { - segments = {static_cast(-1), -h, static_cast(0), h, static_cast(1)}; - } - - return segments; - } - - template - std::vector SVEHintsLogistic::segments_y() const + std::vector SVEHintsLogistic::segments_y() const { - return {static_cast(-1), static_cast(0), static_cast(1)}; + return {-1.0, 0.0, 1.0}; } int SVEHintsLogistic::nsvals() const @@ -835,9 +816,5 @@ namespace sparseir return segments; } - std::vector SVEHintsLogistic::segments_y() const - { - return {-1.0, 0.0, 1.0}; - } } // namespace sparseir \ No newline at end of file From 2a4fb1d39818420d5a862b4ae2a5a519d9b32653 Mon Sep 17 00:00:00 2001 From: SatoshiTerasaki Date: Thu, 28 Nov 2024 12:18:29 +0900 Subject: [PATCH 10/15] Fix C++ implementation by hand --- include/sparseir/kernel.hpp | 189 +++++++++++++++++++++--------------- 1 file changed, 109 insertions(+), 80 deletions(-) diff --git a/include/sparseir/kernel.hpp b/include/sparseir/kernel.hpp index 6dd6a14..ed4a9f6 100644 --- a/include/sparseir/kernel.hpp +++ b/include/sparseir/kernel.hpp @@ -17,7 +17,7 @@ namespace sparseir { // Forward declaration of ReducedKernel class ReducedKernel; - class SVEHints; + class AbstractSVEHints; class SVEHintsLogistic; class SVEHintsRegularizedBose; @@ -58,7 +58,7 @@ namespace sparseir virtual double operator()(double x, double y, double x_plus = std::numeric_limits::quiet_NaN(), double x_minus = std::numeric_limits::quiet_NaN()) const = 0; - virtual std::shared_ptr sve_hints(double epsilon) const = 0; + // virtual auto sve_hints(double epsilon) const = 0; /** * @brief Return symmetrized kernel K(x, y) + sign * K(x, -y). @@ -694,37 +694,128 @@ namespace sparseir namespace sparseir { + /* + AbstractSVEHints - class SVEHints + Discretization hints for singular value expansion of a given kernel. + */ + class AbstractSVEHints { public: - virtual ~SVEHints() = default; + virtual ~AbstractSVEHints() = default; // Functions to compute segments for x and y virtual std::vector segments_x() const = 0; - virtual std::vector segments_x_ddouble() const = 0; - virtual std::vector segments_y() const = 0; - virtual std::vector segments_y_ddouble() const = 0; // Additional methods if needed virtual int nsvals() const = 0; virtual int ngauss() const = 0; }; - class SVEHintsLogistic : public SVEHints + class SVEHintsLogistic final : public AbstractSVEHints { public: SVEHintsLogistic(const LogisticKernel &kernel, double epsilon) : kernel_(kernel), epsilon_(epsilon) {} - std::vector segments_x() const override; - std::vector segments_y() const override; + std::vector segments_x() const override{ + int nzeros = std::max(static_cast(std::round(15 * std::log10(kernel_.lambda_))), 1); + + // Create a range of values + std::vector temp(nzeros); + for (int i = 0; i < nzeros; ++i) { + temp[i] = (double)0.143 * i; + } - std::vector segments_x_ddouble() const override; - std::vector segments_y_ddouble() const override; + // Calculate diffs using the inverse hyperbolic cosine + std::vector diffs(nzeros); + for (int i = 0; i < nzeros; ++i) { + diffs[i] = 1.0 / std::cosh(temp[i]); + } - int nsvals() const override; + // Calculate cumulative sum of diffs + std::vector zeros(nzeros); + zeros[0] = diffs[0]; + for (int i = 1; i < nzeros; ++i) { + zeros[i] = zeros[i - 1] + diffs[i]; + } + + // Normalize zeros + double last_zero = zeros.back(); + for (int i = 0; i < nzeros; ++i) { + zeros[i] /= last_zero; + } + + // Create the final segments vector + std::vector segments; + segments.reserve(2 * nzeros + 1); + + // Add reversed zeros, zero, and zeros to segments + segments.insert(segments.end(), zeros.rbegin(), zeros.rend()); + segments.push_back(0.0); + segments.insert(segments.end(), zeros.begin(), zeros.end()); + + return segments; + }; + std::vector segments_y() const override { + // Calculate the number of zeros + int nzeros = std::max(static_cast(std::round(20 * std::log10(kernel_.lambda_))), 2); + + // Initial differences + std::vector diffs = { + 0.01523, 0.03314, 0.04848, 0.05987, 0.06703, 0.07028, 0.07030, + 0.06791, 0.06391, 0.05896, 0.05358, 0.04814, 0.04288, 0.03795, + 0.03342, 0.02932, 0.02565, 0.02239, 0.01951, 0.01699 + }; + + // Truncate diffs if necessary + if (nzeros < diffs.size()) { + diffs.resize(nzeros); + } + + // Calculate trailing differences + for (int i = 20; i < nzeros; ++i) { + double x = (double)0.141 * i; + diffs.push_back(0.25 * std::exp(-x)); + } + + // Calculate cumulative sum of diffs + std::vector zeros(nzeros); + zeros[0] = diffs[0]; + for (int i = 1; i < nzeros; ++i) { + zeros[i] = zeros[i - 1] + diffs[i]; + } + + // Normalize zeros + double last_zero = zeros.back(); + for (int i = 0; i < nzeros; ++i) { + zeros[i] /= last_zero; + } + + // Adjust zeros + for (int i = 0; i < nzeros; ++i) { + zeros[i] -= 1.0; + } + + // Create the final segments vector + std::vector segments; + segments.reserve(2 * nzeros + 3); + + // Add -1, zeros, 0, reversed zeros, and 1 to segments + segments.push_back(-1.0); + segments.insert(segments.end(), zeros.begin(), zeros.end()); + segments.push_back(0.0); + segments.insert(segments.end(), zeros.rbegin(), zeros.rend()); + segments.push_back(1.0); + + return segments; + } + + int nsvals() const override { + double log10_Lambda = std::max(1.0, std::log10(kernel_.lambda_)); + return static_cast(std::round((25 + log10_Lambda) * log10_Lambda)); + } int ngauss() const override; private: @@ -732,7 +823,7 @@ namespace sparseir double epsilon_; }; - class SVEHintsRegularizedBose : public SVEHints + class SVEHintsRegularizedBose : public AbstractSVEHints { public: SVEHintsRegularizedBose(const RegularizedBoseKernel &kernel, double epsilon) @@ -744,77 +835,15 @@ namespace sparseir template std::vector segments_y() const; - int nsvals() const override; + int nsvals() const override{ + double log10_Lambda = std::max(1.0, std::log10(kernel_.lambda_)); + return static_cast(std::round(28 * log10_Lambda)); + } int ngauss() const override; - std::vector segments_x_ddouble() const override; - std::vector segments_y_ddouble() const override; - private: const RegularizedBoseKernel &kernel_; double epsilon_; }; - std::vector SVEHintsLogistic::segments_y() const - { - return {-1.0, 0.0, 1.0}; - } - - int SVEHintsLogistic::nsvals() const - { - // Implement as needed - return 0; // Placeholder - } - - int SVEHintsLogistic::ngauss() const - { - // Implement as needed - return 0; // Placeholder - } - - std::vector SVEHintsLogistic::segments_x() const - { - int nzeros = std::max(static_cast(std::round(15 * std::log10(kernel_.lambda_))), 1); - - std::vector temp(nzeros); - double coeff = 0.143; - for (int i = 0; i < nzeros; ++i) - { - temp[i] = coeff * static_cast(i); - } - - std::vector diffs(nzeros); - for (int i = 0; i < nzeros; ++i) - { - diffs[i] = 1.0 / std::cosh(temp[i]); - } - - std::vector zeros(nzeros); - zeros[0] = diffs[0]; - for (int i = 1; i < nzeros; ++i) - { - zeros[i] = zeros[i - 1] + diffs[i]; - } - - // Normalize zeros - double last_zero = zeros.back(); - for (double &z : zeros) - { - z /= last_zero; - } - - // Create segments: [-reverse(zeros); 0; zeros] - std::vector segments; - segments.reserve(2 * nzeros + 1); - for (auto it = zeros.rbegin(); it != zeros.rend(); ++it) - { - segments.push_back(-(*it)); - } - segments.push_back(0.0); - segments.insert(segments.end(), zeros.begin(), zeros.end()); - - return segments; - } - - } // namespace sparseir \ No newline at end of file From a30f7cc22328456ef039a128690cabb6b431a6cc Mon Sep 17 00:00:00 2001 From: SatoshiTerasaki Date: Thu, 28 Nov 2024 13:36:23 +0900 Subject: [PATCH 11/15] Fix segments_x --- include/sparseir/gauss.hpp | 5 +-- include/sparseir/kernel.hpp | 67 ++++++++++++++++++++++++++----------- test/kernel.cxx | 28 +++++++++------- 3 files changed, 66 insertions(+), 34 deletions(-) diff --git a/include/sparseir/gauss.hpp b/include/sparseir/gauss.hpp index a6576d4..6a7eb6c 100644 --- a/include/sparseir/gauss.hpp +++ b/include/sparseir/gauss.hpp @@ -75,13 +75,14 @@ class Rule { return Rule(x, new_w, x_forward, x_backward, a, b); } - Rule piecewise(const std::vector& edges) const { + template + Rule piecewise(const std::vector& edges) const { if (!is_sorted(edges.begin(), edges.end())) { throw std::invalid_argument("segments ends must be ordered ascendingly"); } std::vector> rules; for (size_t i = 0; i < edges.size() - 1; ++i) { - rules.push_back(reseat(edges[i], edges[i + 1])); + rules.push_back(reseat(T(edges[i]), T(edges[i + 1]))); } return join(rules); } diff --git a/include/sparseir/kernel.hpp b/include/sparseir/kernel.hpp index ed4a9f6..ae3f3bc 100644 --- a/include/sparseir/kernel.hpp +++ b/include/sparseir/kernel.hpp @@ -7,10 +7,6 @@ #include #include #include - -#pragma once - -#include #include namespace sparseir @@ -58,7 +54,12 @@ namespace sparseir virtual double operator()(double x, double y, double x_plus = std::numeric_limits::quiet_NaN(), double x_minus = std::numeric_limits::quiet_NaN()) const = 0; - // virtual auto sve_hints(double epsilon) const = 0; + /* + sve_hints(double epsilon) const + { + return std::make_shared(*this, epsilon); + } + */ /** * @brief Return symmetrized kernel K(x, y) + sign * K(x, -y). @@ -219,9 +220,9 @@ namespace sparseir return compute(u_plus, u_minus, v); } - /* // Inside class LogisticKernel definition - std::shared_ptr sve_hints(double epsilon) const override + /* + std::shared_ptr sve_hints(double epsilon) const { return std::make_shared(*this, epsilon); } @@ -398,9 +399,9 @@ namespace sparseir return compute(u_plus, u_minus, v); } - /* // Inside class RegularizedBoseKernel definition - std::shared_ptr sve_hints(double epsilon) const override + /* + std::shared_ptr sve_hints(double epsilon) const { return std::make_shared(*this, epsilon); } @@ -748,14 +749,11 @@ namespace sparseir } // Create the final segments vector - std::vector segments; - segments.reserve(2 * nzeros + 1); - - // Add reversed zeros, zero, and zeros to segments - segments.insert(segments.end(), zeros.rbegin(), zeros.rend()); - segments.push_back(0.0); - segments.insert(segments.end(), zeros.begin(), zeros.end()); - + std::vector segments(2 * nzeros + 1, 0); + for (int i = 0; i < nzeros; ++i) { + segments[i] = -zeros[nzeros - i - 1]; + segments[nzeros + i + 1] = zeros[i]; + } return segments; }; std::vector segments_y() const override { @@ -816,7 +814,9 @@ namespace sparseir double log10_Lambda = std::max(1.0, std::log10(kernel_.lambda_)); return static_cast(std::round((25 + log10_Lambda) * log10_Lambda)); } - int ngauss() const override; + int ngauss() const override{ + return epsilon_ >= std::sqrt(std::numeric_limits::epsilon()) ? 10 : 16; + }; private: const LogisticKernel &kernel_; @@ -846,4 +846,33 @@ namespace sparseir double epsilon_; }; -} // namespace sparseir \ No newline at end of file +} // namespace sparseir + +namespace sparseir{ + + // Function to provide SVE hints + inline SVEHintsLogistic sve_hints(const LogisticKernel& kernel, double epsilon) { + return SVEHintsLogistic(kernel, epsilon); + } + +/* +function ngauss end +ngauss(hints::SVEHintsLogistic) = hints.ε ≥ sqrt(eps()) ? 10 : 16 +ngauss(hints::SVEHintsRegularizedBose) = hints.ε ≥ sqrt(eps()) ? 10 : 16 +ngauss(hints::SVEHintsReduced) = ngauss(hints.inner_hints) +*/ + + /* + SVEHintsRegularizedBose sve_hints(const RegularizedBoseKernel& kernel, double epsilon) { + return SVEHintsRegularizedBose(kernel, epsilon); + } + */ + + /* + std::shared_ptr sve_hints(const AbstractReducedKernel& kernel, double epsilon) { + // Assume kernel.inner is a method that returns the inner kernel + auto innerHints = sve_hints(kernel.inner(), epsilon); + return std::make_shared(innerHints); + } + */ +} \ No newline at end of file diff --git a/test/kernel.cxx b/test/kernel.cxx index a311a76..50da4d1 100644 --- a/test/kernel.cxx +++ b/test/kernel.cxx @@ -4,6 +4,7 @@ #include #include #include +#include #include // Include your sparseir headers @@ -17,31 +18,31 @@ TEST_CASE("Kernel Accuracy Test") using T_x = double; // List of kernels to test - std::vector> kernels = { - std::make_shared(9.0), - std::make_shared(8.0), - std::make_shared(120000.0), - std::make_shared(127500.0), + std::vector kernels = { + sparseir::LogisticKernel(9.0), + //sparseir::RegularizedBoseKernel(8.0), + sparseir::LogisticKernel(120000.0), + //sparseir::RegularizedBoseKernel(127500.0), // Symmetrized kernels - std::make_shared(40000.0)->get_symmetrized(-1), - std::make_shared(35000.0)->get_symmetrized(-1), + //sparseir::LogisticKernel(40000.0)->get_symmetrized(-1), + //std::make_shared(35000.0)->get_symmetrized(-1), }; - for (const auto &K_ptr : kernels) + for (const auto &K : kernels) { - const auto &K = *K_ptr; - // Convert Rule to type T sparseir::Rule ddouble_rule = sparseir::legendre(10); sparseir::Rule double_rule = sparseir::convert(ddouble_rule); sparseir::Rule rule = sparseir::convert(double_rule); // Obtain SVE hints for the kernel - auto hints = K.sve_hints(2.2e-16); + auto hints = sparseir::sve_hints(K, 2.2e-16); // Generate piecewise Gaussian quadrature rules for x and y - auto gauss_x = sparseir::piecewise(rule, hints.template segments_x()); - auto gauss_y = sparseir::piecewise(rule, hints.template segments_y()); + auto gauss_x = rule.piecewise(hints.segments_x()); + /* + auto gauss_y = rule.piecewise(hints.segments_y()); + T epsilon = std::numeric_limits::epsilon(); T tiny = std::numeric_limits::min() / epsilon; @@ -65,6 +66,7 @@ TEST_CASE("Kernel Accuracy Test") .select(T(1.0), result.array() / result_x.template cast().array()); REQUIRE((reldiff - T(1.0)).cwiseAbs().maxCoeff() <= 100 * epsilon); + */ } } From f851e3beeeaa318566348cef06884ca2ca69affd Mon Sep 17 00:00:00 2001 From: SatoshiTerasaki Date: Thu, 28 Nov 2024 13:52:25 +0900 Subject: [PATCH 12/15] Fix gauss implementation --- include/sparseir/kernel.hpp | 22 ++++++++++++---------- test/kernel.cxx | 3 +-- 2 files changed, 13 insertions(+), 12 deletions(-) diff --git a/include/sparseir/kernel.hpp b/include/sparseir/kernel.hpp index ae3f3bc..191baed 100644 --- a/include/sparseir/kernel.hpp +++ b/include/sparseir/kernel.hpp @@ -790,6 +790,10 @@ namespace sparseir for (int i = 0; i < nzeros; ++i) { zeros[i] /= last_zero; } + zeros.pop_back(); + + // updated nzeros + nzeros = zeros.size(); // Adjust zeros for (int i = 0; i < nzeros; ++i) { @@ -797,16 +801,14 @@ namespace sparseir } // Create the final segments vector - std::vector segments; - segments.reserve(2 * nzeros + 3); - - // Add -1, zeros, 0, reversed zeros, and 1 to segments - segments.push_back(-1.0); - segments.insert(segments.end(), zeros.begin(), zeros.end()); - segments.push_back(0.0); - segments.insert(segments.end(), zeros.rbegin(), zeros.rend()); - segments.push_back(1.0); - + std::vector segments(2 * nzeros + 3, 0); + for (int i = 0; i < nzeros; ++i) { + segments[1+i] = zeros[i]; + segments[1+nzeros + 1 + i] = -zeros[nzeros-i-1]; + } + segments[0] = -1.0; + segments[1 + nzeros] = 0.0; + segments[2*nzeros + 2] = 1.0; return segments; } diff --git a/test/kernel.cxx b/test/kernel.cxx index 50da4d1..de871bf 100644 --- a/test/kernel.cxx +++ b/test/kernel.cxx @@ -40,9 +40,8 @@ TEST_CASE("Kernel Accuracy Test") // Generate piecewise Gaussian quadrature rules for x and y auto gauss_x = rule.piecewise(hints.segments_x()); - /* auto gauss_y = rule.piecewise(hints.segments_y()); - + /* T epsilon = std::numeric_limits::epsilon(); T tiny = std::numeric_limits::min() / epsilon; From fb9d0b1e4a2a93485199beff51ee00bab8e36730 Mon Sep 17 00:00:00 2001 From: SatoshiTerasaki Date: Thu, 28 Nov 2024 14:04:58 +0900 Subject: [PATCH 13/15] Done!!! for LogisticKernel --- include/sparseir/kernel.hpp | 28 ++++++++++++++++++++++++++++ test/kernel.cxx | 3 --- 2 files changed, 28 insertions(+), 3 deletions(-) diff --git a/include/sparseir/kernel.hpp b/include/sparseir/kernel.hpp index 191baed..05935c9 100644 --- a/include/sparseir/kernel.hpp +++ b/include/sparseir/kernel.hpp @@ -852,6 +852,34 @@ namespace sparseir namespace sparseir{ + // Function to compute matrix from Gauss rules + template + Eigen::MatrixX matrix_from_gauss( + const AbstractKernel &kernel, const Rule& gauss_x, const Rule& gauss_y) { + + size_t n = gauss_x.x.size(); + size_t m = gauss_y.x.size(); + Eigen::Matrix res(n, m); + + // Parallelize using threads + std::vector threads; + for (size_t i = 0; i < n; ++i) { + threads.emplace_back([&, i]() { + for (size_t j = 0; j < m; ++j) { + res(i, j) = kernel(gauss_x.x[i], gauss_y.x[j], + gauss_x.x_forward[i], gauss_x.x_backward[i]); + } + }); + } + + // Join threads + for (auto& thread : threads) { + thread.join(); + } + + return res; + } + // Function to provide SVE hints inline SVEHintsLogistic sve_hints(const LogisticKernel& kernel, double epsilon) { return SVEHintsLogistic(kernel, epsilon); diff --git a/test/kernel.cxx b/test/kernel.cxx index de871bf..a2451f6 100644 --- a/test/kernel.cxx +++ b/test/kernel.cxx @@ -41,7 +41,6 @@ TEST_CASE("Kernel Accuracy Test") // Generate piecewise Gaussian quadrature rules for x and y auto gauss_x = rule.piecewise(hints.segments_x()); auto gauss_y = rule.piecewise(hints.segments_y()); - /* T epsilon = std::numeric_limits::epsilon(); T tiny = std::numeric_limits::min() / epsilon; @@ -55,7 +54,6 @@ TEST_CASE("Kernel Accuracy Test") // Compute the matrix in higher precision Eigen::Matrix result_x = sparseir::matrix_from_gauss(K, gauss_x_Tx, gauss_y_Tx); - T_x magn = result_x.cwiseAbs().maxCoeff(); // Check that the difference is within tolerance @@ -65,7 +63,6 @@ TEST_CASE("Kernel Accuracy Test") .select(T(1.0), result.array() / result_x.template cast().array()); REQUIRE((reldiff - T(1.0)).cwiseAbs().maxCoeff() <= 100 * epsilon); - */ } } From f365e58de88f9daccc8a613a7a46edb398eee7ab Mon Sep 17 00:00:00 2001 From: SatoshiTerasaki Date: Thu, 28 Nov 2024 14:16:46 +0900 Subject: [PATCH 14/15] Implement ngauss for RegularizedBoseKernel --- include/sparseir/kernel.hpp | 16 ++++++++------ test/kernel.cxx | 43 +++++++++++++++++++++---------------- 2 files changed, 35 insertions(+), 24 deletions(-) diff --git a/include/sparseir/kernel.hpp b/include/sparseir/kernel.hpp index 05935c9..9ab20ef 100644 --- a/include/sparseir/kernel.hpp +++ b/include/sparseir/kernel.hpp @@ -841,7 +841,9 @@ namespace sparseir double log10_Lambda = std::max(1.0, std::log10(kernel_.lambda_)); return static_cast(std::round(28 * log10_Lambda)); } - int ngauss() const override; + int ngauss() const override { + return epsilon_ >= std::sqrt(std::numeric_limits::epsilon()) ? 10 : 16; + }; private: const RegularizedBoseKernel &kernel_; @@ -885,6 +887,12 @@ namespace sparseir{ return SVEHintsLogistic(kernel, epsilon); } + /* + inline SVEHintsRegularizedBose sve_hints(const RegularizedBoseKernel& kernel, double epsilon) { + return SVEHintsRegularizedBose(kernel, epsilon); + } + */ + /* function ngauss end ngauss(hints::SVEHintsLogistic) = hints.ε ≥ sqrt(eps()) ? 10 : 16 @@ -892,11 +900,7 @@ ngauss(hints::SVEHintsRegularizedBose) = hints.ε ≥ sqrt(eps()) ? 10 : 16 ngauss(hints::SVEHintsReduced) = ngauss(hints.inner_hints) */ - /* - SVEHintsRegularizedBose sve_hints(const RegularizedBoseKernel& kernel, double epsilon) { - return SVEHintsRegularizedBose(kernel, epsilon); - } - */ + /* std::shared_ptr sve_hints(const AbstractReducedKernel& kernel, double epsilon) { diff --git a/test/kernel.cxx b/test/kernel.cxx index a2451f6..8dd2775 100644 --- a/test/kernel.cxx +++ b/test/kernel.cxx @@ -12,24 +12,10 @@ using xprec::DDouble; -TEST_CASE("Kernel Accuracy Test") -{ - using T = float; - using T_x = double; - - // List of kernels to test - std::vector kernels = { - sparseir::LogisticKernel(9.0), - //sparseir::RegularizedBoseKernel(8.0), - sparseir::LogisticKernel(120000.0), - //sparseir::RegularizedBoseKernel(127500.0), - // Symmetrized kernels - //sparseir::LogisticKernel(40000.0)->get_symmetrized(-1), - //std::make_shared(35000.0)->get_symmetrized(-1), - }; - - for (const auto &K : kernels) - { +template +bool kernel_accuracy_test(Kernel &K){ + using T = float; + using T_x = double; // Convert Rule to type T sparseir::Rule ddouble_rule = sparseir::legendre(10); sparseir::Rule double_rule = sparseir::convert(ddouble_rule); @@ -63,6 +49,27 @@ TEST_CASE("Kernel Accuracy Test") .select(T(1.0), result.array() / result_x.template cast().array()); REQUIRE((reldiff - T(1.0)).cwiseAbs().maxCoeff() <= 100 * epsilon); + return true; +} + + +TEST_CASE("Kernel Accuracy Test") +{ + { + // List of kernels to test + std::vector kernels = { + sparseir::LogisticKernel(9.0), + //sparseir::RegularizedBoseKernel(8.0), + sparseir::LogisticKernel(120000.0), + //sparseir::RegularizedBoseKernel(127500.0), + // Symmetrized kernels + //sparseir::LogisticKernel(40000.0)->get_symmetrized(-1), + //std::make_shared(35000.0)->get_symmetrized(-1), + }; + for (const auto &K : kernels) + { + REQUIRE(kernel_accuracy_test(K)); + } } } From 53cc9d3e88b4c41febb90e57a8134e0bad0381cf Mon Sep 17 00:00:00 2001 From: SatoshiTerasaki Date: Thu, 28 Nov 2024 15:27:09 +0900 Subject: [PATCH 15/15] Pass tests for RegularizedBoseKernel --- include/sparseir/gauss.hpp | 2 +- include/sparseir/kernel.hpp | 114 ++++++++++++++++++++++++++++-------- test/kernel.cxx | 91 +++++++++++++++++++--------- 3 files changed, 153 insertions(+), 54 deletions(-) diff --git a/include/sparseir/gauss.hpp b/include/sparseir/gauss.hpp index 6a7eb6c..ec026cf 100644 --- a/include/sparseir/gauss.hpp +++ b/include/sparseir/gauss.hpp @@ -77,7 +77,7 @@ class Rule { template Rule piecewise(const std::vector& edges) const { - if (!is_sorted(edges.begin(), edges.end())) { + if (!std::is_sorted(edges.begin(), edges.end())) { throw std::invalid_argument("segments ends must be ordered ascendingly"); } std::vector> rules; diff --git a/include/sparseir/kernel.hpp b/include/sparseir/kernel.hpp index 9ab20ef..b834cb6 100644 --- a/include/sparseir/kernel.hpp +++ b/include/sparseir/kernel.hpp @@ -706,8 +706,10 @@ namespace sparseir virtual ~AbstractSVEHints() = default; // Functions to compute segments for x and y - virtual std::vector segments_x() const = 0; - virtual std::vector segments_y() const = 0; + //template + //virtual std::vector segments_x() const = 0; + //template + //virtual std::vector segments_y() const = 0; // Additional methods if needed virtual int nsvals() const = 0; @@ -720,48 +722,51 @@ namespace sparseir SVEHintsLogistic(const LogisticKernel &kernel, double epsilon) : kernel_(kernel), epsilon_(epsilon) {} - std::vector segments_x() const override{ + template + std::vector segments_x() const { int nzeros = std::max(static_cast(std::round(15 * std::log10(kernel_.lambda_))), 1); // Create a range of values - std::vector temp(nzeros); + std::vector temp(nzeros); for (int i = 0; i < nzeros; ++i) { - temp[i] = (double)0.143 * i; + temp[i] = (T)0.143 * i; } // Calculate diffs using the inverse hyperbolic cosine - std::vector diffs(nzeros); + std::vector diffs(nzeros); for (int i = 0; i < nzeros; ++i) { diffs[i] = 1.0 / std::cosh(temp[i]); } // Calculate cumulative sum of diffs - std::vector zeros(nzeros); + std::vector zeros(nzeros); zeros[0] = diffs[0]; for (int i = 1; i < nzeros; ++i) { zeros[i] = zeros[i - 1] + diffs[i]; } // Normalize zeros - double last_zero = zeros.back(); + T last_zero = zeros.back(); for (int i = 0; i < nzeros; ++i) { zeros[i] /= last_zero; } // Create the final segments vector - std::vector segments(2 * nzeros + 1, 0); + std::vector segments(2 * nzeros + 1, 0); for (int i = 0; i < nzeros; ++i) { segments[i] = -zeros[nzeros - i - 1]; segments[nzeros + i + 1] = zeros[i]; } return segments; }; - std::vector segments_y() const override { + + template + std::vector segments_y() const { // Calculate the number of zeros int nzeros = std::max(static_cast(std::round(20 * std::log10(kernel_.lambda_))), 2); // Initial differences - std::vector diffs = { + std::vector diffs = { 0.01523, 0.03314, 0.04848, 0.05987, 0.06703, 0.07028, 0.07030, 0.06791, 0.06391, 0.05896, 0.05358, 0.04814, 0.04288, 0.03795, 0.03342, 0.02932, 0.02565, 0.02239, 0.01951, 0.01699 @@ -774,19 +779,19 @@ namespace sparseir // Calculate trailing differences for (int i = 20; i < nzeros; ++i) { - double x = (double)0.141 * i; + T x = (T)0.141 * i; diffs.push_back(0.25 * std::exp(-x)); } // Calculate cumulative sum of diffs - std::vector zeros(nzeros); + std::vector zeros(nzeros); zeros[0] = diffs[0]; for (int i = 1; i < nzeros; ++i) { zeros[i] = zeros[i - 1] + diffs[i]; } // Normalize zeros - double last_zero = zeros.back(); + T last_zero = zeros.back(); for (int i = 0; i < nzeros; ++i) { zeros[i] /= last_zero; } @@ -801,14 +806,14 @@ namespace sparseir } // Create the final segments vector - std::vector segments(2 * nzeros + 3, 0); + std::vector segments(2 * nzeros + 3, 0); for (int i = 0; i < nzeros; ++i) { segments[1+i] = zeros[i]; segments[1+nzeros + 1 + i] = -zeros[nzeros-i-1]; } - segments[0] = -1.0; - segments[1 + nzeros] = 0.0; - segments[2*nzeros + 2] = 1.0; + segments[0] = -T(1.0); + segments[1 + nzeros] = T( 0.0); + segments[2*nzeros + 2] = T(1.0); return segments; } @@ -831,11 +836,74 @@ namespace sparseir SVEHintsRegularizedBose(const RegularizedBoseKernel &kernel, double epsilon) : kernel_(kernel), epsilon_(epsilon) {} - template - std::vector segments_x() const; + template + std::vector segments_x() const { + int nzeros = std::max(static_cast(std::round(15 * std::log10(kernel_.lambda_))), 15); + std::vector temp(nzeros); + std::vector diffs(nzeros); + std::vector zeros(nzeros); + + for (int i = 0; i < nzeros; ++i) { + temp[i] = T(0.18) * i; + diffs[i] = 1.0 / std::cosh(temp[i]); + } + + std::partial_sum(diffs.begin(), diffs.end(), zeros.begin()); + T last_zero = zeros.back(); + std::transform(zeros.begin(), zeros.end(), zeros.begin(), [last_zero](T z) { return z / last_zero; }); - template - std::vector segments_y() const; + std::vector result(2 * nzeros + 1, T(0)); + + for (int i = 0; i < nzeros; ++i) { + result[i] = -zeros[nzeros - i - 1]; + result[nzeros + i + 1] = zeros[i]; + } + + return result; + } + + template + std::vector segments_y() const { + int nzeros = std::max(static_cast(std::round(20 * std::log10(kernel_.lambda_))), 20); + std::vector diffs(nzeros); + + for (int j = 0; j < nzeros; ++j) { + diffs[j] = T(0.12) / std::exp(0.0337 * j * std::log(j + 1)); + } + + // Calculate cumulative sum of diffs + std::vector zeros(nzeros); + zeros[0] = diffs[0]; + for (int i = 1; i < nzeros; ++i) { + zeros[i] = zeros[i - 1] + diffs[i]; + } + + // Normalize zeros + T last_zero = zeros.back(); + for (int i = 0; i < nzeros; ++i) { + zeros[i] /= last_zero; + } + zeros.pop_back(); + + // updated nzeros + nzeros = zeros.size(); + + // Adjust zeros + for (int i = 0; i < nzeros; ++i) { + zeros[i] -= 1.0; + } + + std::vector result(2 * nzeros + 3, T(0)); + + for (int i = 0; i < nzeros; ++i) { + result[1+i] = zeros[i]; + result[1+nzeros + 1 + i] = -zeros[nzeros-i-1]; + } + result[0] = T(-1); + result[1 + nzeros] = T(0); + result[2*nzeros + 2] = T(1); + return result; + } int nsvals() const override{ double log10_Lambda = std::max(1.0, std::log10(kernel_.lambda_)); @@ -887,11 +955,9 @@ namespace sparseir{ return SVEHintsLogistic(kernel, epsilon); } - /* inline SVEHintsRegularizedBose sve_hints(const RegularizedBoseKernel& kernel, double epsilon) { return SVEHintsRegularizedBose(kernel, epsilon); } - */ /* function ngauss end diff --git a/test/kernel.cxx b/test/kernel.cxx index 8dd2775..d1790db 100644 --- a/test/kernel.cxx +++ b/test/kernel.cxx @@ -13,45 +13,52 @@ using xprec::DDouble; template -bool kernel_accuracy_test(Kernel &K){ - using T = float; - using T_x = double; - // Convert Rule to type T - sparseir::Rule ddouble_rule = sparseir::legendre(10); - sparseir::Rule double_rule = sparseir::convert(ddouble_rule); - sparseir::Rule rule = sparseir::convert(double_rule); +bool kernel_accuracy_test(Kernel &K) { + using T = float; + using T_x = double; - // Obtain SVE hints for the kernel - auto hints = sparseir::sve_hints(K, 2.2e-16); + // Convert Rule to type T + auto ddouble_rule = sparseir::legendre(10); + auto double_rule = sparseir::convert(ddouble_rule); + auto rule = sparseir::convert(double_rule); - // Generate piecewise Gaussian quadrature rules for x and y - auto gauss_x = rule.piecewise(hints.segments_x()); - auto gauss_y = rule.piecewise(hints.segments_y()); + // Obtain SVE hints for the kernel + auto hints = sparseir::sve_hints(K, 2.2e-16); - T epsilon = std::numeric_limits::epsilon(); - T tiny = std::numeric_limits::min() / epsilon; + // Generate piecewise Gaussian quadrature rules for x and y - // Compute the matrix from Gaussian quadrature - Eigen::Matrix result = sparseir::matrix_from_gauss(K, gauss_x, gauss_y); + auto sx = hints.segments_x(); + auto sy = hints.segments_y(); - // Convert gauss_x and gauss_y to higher precision T_x - auto gauss_x_Tx = sparseir::convert(gauss_x); - auto gauss_y_Tx = sparseir::convert(gauss_y); + REQUIRE(std::is_sorted(sx.begin(), sx.end())); + REQUIRE(std::is_sorted(sy.begin(), sy.end())); - // Compute the matrix in higher precision - Eigen::Matrix result_x = sparseir::matrix_from_gauss(K, gauss_x_Tx, gauss_y_Tx); - T_x magn = result_x.cwiseAbs().maxCoeff(); + auto gauss_x = rule.piecewise(sx); + auto gauss_y = rule.piecewise(sy); - // Check that the difference is within tolerance - REQUIRE((result.template cast() - result_x).cwiseAbs().maxCoeff() <= 2 * magn * epsilon); + T epsilon = std::numeric_limits::epsilon(); + T tiny = std::numeric_limits::min() / epsilon; - Eigen::Array reldiff = (result.cwiseAbs().array() < tiny) - .select(T(1.0), result.array() / result_x.template cast().array()); + // Compute the matrix from Gaussian quadrature + auto result = sparseir::matrix_from_gauss(K, gauss_x, gauss_y); - REQUIRE((reldiff - T(1.0)).cwiseAbs().maxCoeff() <= 100 * epsilon); - return true; -} + // Convert gauss_x and gauss_y to higher precision T_x + auto gauss_x_Tx = sparseir::convert(gauss_x); + auto gauss_y_Tx = sparseir::convert(gauss_y); + + // Compute the matrix in higher precision + auto result_x = sparseir::matrix_from_gauss(K, gauss_x_Tx, gauss_y_Tx); + T_x magn = result_x.cwiseAbs().maxCoeff(); + + // Check that the difference is within tolerance + REQUIRE((result.template cast() - result_x).cwiseAbs().maxCoeff() <= 2 * magn * epsilon); + auto reldiff = (result.cwiseAbs().array() < tiny) + .select(T(1.0), result.array() / result_x.template cast().array()); + + REQUIRE((reldiff - T(1.0)).cwiseAbs().maxCoeff() <= 100 * epsilon); + return true; +} TEST_CASE("Kernel Accuracy Test") { @@ -71,6 +78,32 @@ TEST_CASE("Kernel Accuracy Test") REQUIRE(kernel_accuracy_test(K)); } } + + { + // List of kernels to test + std::vector kernels = { + sparseir::RegularizedBoseKernel(8.0), + sparseir::RegularizedBoseKernel(127500.0), + }; + for (const auto &K : kernels) + { + REQUIRE(kernel_accuracy_test(K)); + } + } + /* + { + // List of kernels to test + std::vector kernels = { + // Symmetrized kernels + sparseir::LogisticKernel(40000.0).get_symmetrized(-1), + sparseir::RegularizedBoseKernel(35000.0).get_symmetrized(-1), + }; + for (const auto &K : kernels) + { + REQUIRE(kernel_accuracy_test(K)); + } + } + */ } TEST_CASE("Kernel Singularity Test")