diff --git a/include/sparseir/gauss.hpp b/include/sparseir/gauss.hpp index 4fe30ea..ec026cf 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 { - if (!is_sorted(edges.begin(), edges.end())) { + template + Rule piecewise(const std::vector& edges) const { + if (!std::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); } @@ -192,4 +193,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, x_forward, x_backward, a, b); +} + } // namespace sparseir \ No newline at end of file diff --git a/include/sparseir/kernel.hpp b/include/sparseir/kernel.hpp new file mode 100644 index 0000000..b834cb6 --- /dev/null +++ b/include/sparseir/kernel.hpp @@ -0,0 +1,978 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace sparseir +{ + // Forward declaration of ReducedKernel + class ReducedKernel; + class AbstractSVEHints; + class SVEHintsLogistic; + class SVEHintsRegularizedBose; + + + /** + * @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; + + /* + sve_hints(double epsilon) const + { + return std::make_shared(*this, epsilon); + } + */ + + /** + * @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); + } + + // Inside class LogisticKernel definition + /* + std::shared_ptr sve_hints(double epsilon) const + { + return std::make_shared(*this, epsilon); + } + */ + + /** + * @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); + } + + // Inside class RegularizedBoseKernel definition + /* + std::shared_ptr sve_hints(double epsilon) const + { + return std::make_shared(*this, epsilon); + } + */ + + /** + * @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 + +namespace sparseir +{ + /* + AbstractSVEHints + + Discretization hints for singular value expansion of a given kernel. + */ + class AbstractSVEHints + { + public: + virtual ~AbstractSVEHints() = default; + + // Functions to compute segments for x and y + //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; + virtual int ngauss() const = 0; + }; + + class SVEHintsLogistic final : public AbstractSVEHints + { + public: + SVEHintsLogistic(const LogisticKernel &kernel, double epsilon) + : kernel_(kernel), epsilon_(epsilon) {} + + 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); + for (int i = 0; i < nzeros; ++i) { + temp[i] = (T)0.143 * i; + } + + // 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]); + } + + // 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; + } + + // Create the final segments vector + 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; + }; + + 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 = { + 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) { + T x = (T)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 + 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; + } + + // Create the final segments vector + 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] = -T(1.0); + segments[1 + nzeros] = T( 0.0); + segments[2*nzeros + 2] = T(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{ + return epsilon_ >= std::sqrt(std::numeric_limits::epsilon()) ? 10 : 16; + }; + + private: + const LogisticKernel &kernel_; + double epsilon_; + }; + + class SVEHintsRegularizedBose : public AbstractSVEHints + { + public: + SVEHintsRegularizedBose(const RegularizedBoseKernel &kernel, double epsilon) + : kernel_(kernel), epsilon_(epsilon) {} + + 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; }); + + 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_)); + return static_cast(std::round(28 * log10_Lambda)); + } + int ngauss() const override { + return epsilon_ >= std::sqrt(std::numeric_limits::epsilon()) ? 10 : 16; + }; + + private: + const RegularizedBoseKernel &kernel_; + double epsilon_; + }; + +} // 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); + } + + inline SVEHintsRegularizedBose sve_hints(const RegularizedBoseKernel& kernel, double epsilon) { + return SVEHintsRegularizedBose(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) +*/ + + + + /* + 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/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..d1790db --- /dev/null +++ b/test/kernel.cxx @@ -0,0 +1,133 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +// Include your sparseir headers +#include + +using xprec::DDouble; + +template +bool kernel_accuracy_test(Kernel &K) { + using T = float; + using T_x = double; + + // Convert Rule to type T + auto ddouble_rule = sparseir::legendre(10); + auto double_rule = sparseir::convert(ddouble_rule); + auto rule = sparseir::convert(double_rule); + + // Obtain SVE hints for the kernel + auto hints = sparseir::sve_hints(K, 2.2e-16); + + // Generate piecewise Gaussian quadrature rules for x and y + + auto sx = hints.segments_x(); + auto sy = hints.segments_y(); + + REQUIRE(std::is_sorted(sx.begin(), sx.end())); + REQUIRE(std::is_sorted(sy.begin(), sy.end())); + + auto gauss_x = rule.piecewise(sx); + auto gauss_y = rule.piecewise(sy); + + T epsilon = std::numeric_limits::epsilon(); + T tiny = std::numeric_limits::min() / epsilon; + + // Compute the matrix from Gaussian quadrature + auto result = sparseir::matrix_from_gauss(K, gauss_x, gauss_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); + + // 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") +{ + { + // 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)); + } + } + + { + // 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") +{ + 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) + { + // Generate random x values + std::vector x_values(1000); + for (double &x : x_values) + { + x = dist(rng); + } + + sparseir::RegularizedBoseKernel K(lambda); + + for (double x : x_values) + { + double expected = 1.0 / lambda; + double computed = K(x, 0.0); + REQUIRE(Eigen::internal::isApprox(computed, expected)); + } + } +} \ No newline at end of file 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") {