diff --git a/include/sparseir/kernel.hpp b/include/sparseir/kernel.hpp index b834cb6..ed09d87 100644 --- a/include/sparseir/kernel.hpp +++ b/include/sparseir/kernel.hpp @@ -12,6 +12,7 @@ namespace sparseir { // Forward declaration of ReducedKernel + template class ReducedKernel; class AbstractSVEHints; class SVEHintsLogistic; @@ -37,6 +38,9 @@ namespace sparseir class AbstractKernel { public: + double lambda_; + // Constructor + AbstractKernel(double lambda) : lambda_(lambda) {} /** * @brief Evaluate kernel at point (x, y). * @@ -69,10 +73,10 @@ namespace sparseir * @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"); - } + //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. @@ -152,6 +156,36 @@ namespace sparseir virtual ~AbstractKernel() {} }; + class AbstractReducedKernel : public AbstractKernel + { + public: + int sign; + std::shared_ptr inner; + + // Constructor + AbstractReducedKernel(std::shared_ptr inner_kernel, int sign) + : AbstractKernel(inner_kernel->lambda_), inner(std::move(inner_kernel)), sign(sign) + { + // Validate inputs + if (!inner->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"); + } + } + }; + + inline double callreduced(const AbstractReducedKernel &kernel, double x, double y, double x_plus, double x_minus) + { + x_plus += 1; + auto K_plus = (*kernel.inner)(x, +y, x_plus, x_minus); + auto K_minus = (*kernel.inner)(x, -y, x_minus, x_plus); + return K_plus + kernel.sign * K_minus; + } + /** * @brief Fermionic/bosonic analytical continuation kernel. * @@ -169,16 +203,14 @@ namespace sparseir 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) + LogisticKernel(double lambda) : AbstractKernel(lambda) { - if (lambda_ < 0) + if (lambda < 0) { throw std::domain_error("Kernel cutoff Λ must be non-negative"); } @@ -348,16 +380,14 @@ namespace sparseir 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) + explicit RegularizedBoseKernel(double lambda) : AbstractKernel(lambda) { - if (lambda_ < 0) + if (lambda < 0) { throw std::domain_error("Kernel cutoff Λ must be non-negative"); } @@ -399,14 +429,6 @@ namespace sparseir 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. * @@ -500,37 +522,22 @@ namespace sparseir * @param v Computed v. * @return The value of K(x, y). */ - double compute(double u_plus, double u_minus, double v) const - { + double compute(double u_plus, double u_minus, double v) const{ double absv = std::abs(v); + double enum_val = std::exp(-absv * (v >= 0 ? u_plus : u_minus)); - 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) + // Handle the tricky expression v / (exp(v) - 1) + double denom; + if (absv >= 1e-200) { - denominator = std::expm1(-absv); // exp(-absv) - 1 - value = -1.0 / lambda_ * numerator * (absv / denominator); + denom = absv / std::expm1(-absv); } else { - // Limit as absv -> 0 - value = -1.0 / lambda_ * numerator * (-1.0); + denom = -1; // Assuming T is a floating-point type } - return value; + return -1 / static_cast(this->lambda_) * enum_val * denom; } }; @@ -549,6 +556,7 @@ namespace sparseir * This kernel is what this class represents. The full singular functions can be * reconstructed by (anti-)symmetrically continuing them to the negative axis. */ + template class ReducedKernel : public AbstractKernel { public: @@ -562,7 +570,9 @@ namespace sparseir * @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) + : AbstractKernel(inner_kernel->lambda_), // Initialize base class + inner_kernel_(std::move(inner_kernel)), + sign_(sign) { if (!inner_kernel_->is_centrosymmetric()) { @@ -625,17 +635,6 @@ namespace sparseir 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. * @@ -918,6 +917,108 @@ namespace sparseir double epsilon_; }; + class RegularizedBoseKernelOdd : public AbstractReducedKernel + { + public: + RegularizedBoseKernelOdd(std::shared_ptr inner, int sign) + : AbstractReducedKernel(inner, sign) + { + if (!is_centrosymmetric()) + { + throw std::runtime_error("inner kernel must be centrosymmetric"); + } + if (std::abs(sign) != 1) + { + throw std::domain_error("sign must be -1 or 1"); + } + } + + // Implement the pure virtual function from the parent class + double operator()(double x, double y, double x_plus = std::numeric_limits::quiet_NaN(), + double x_minus = std::numeric_limits::quiet_NaN()) const override + { + double v_half = inner->lambda_ * 0.5 * y; + double xv_half = x * v_half; + bool xy_small = xv_half < 1; + bool sinh_range = 1e-200 < v_half && v_half < 85; + if (xy_small && sinh_range) + { + return y * std::sinh(xv_half) / std::sinh(v_half); + } + else + { + return callreduced(*this, x, x, x_plus, x_minus); + } + } + + // You'll need to implement the isCentrosymmetric function + // Here's a placeholder + bool isCentrosymmetric(RegularizedBoseKernel &kernel) + { + // Implement this function + return true; + } + }; + + + + + class LogisticKernelOdd : public AbstractReducedKernel + { + public: + LogisticKernelOdd(std::shared_ptr inner, int sign) + : AbstractReducedKernel(inner, sign){} + // Implement the pure virtual function from the parent class + double operator()(double x, double y, double x_plus = std::numeric_limits::quiet_NaN(), + double x_minus = std::numeric_limits::quiet_NaN()) const override{ + double v_half = inner->lambda_ * 0.5 * y; + bool xy_small = x * v_half < 1; + bool cosh_finite = v_half < 85; + if (xy_small && cosh_finite) + { + return -std::sinh(v_half * x) / std::cosh(v_half); + } + else + { + return callreduced(*this, x, x, x_plus, x_minus); + } + } + }; + + inline std::shared_ptr get_symmetrized(std::shared_ptr kernel, int sign) + { + return std::make_shared>(kernel, sign); + } + + inline std::shared_ptr get_symmetrized(std::shared_ptr kernel, int sign) + { + if (sign == -1) + { + return std::make_shared(kernel, sign); + } + else + { + return std::make_shared>(kernel, sign); + } + } + + inline std::shared_ptr get_symmetrized(std::shared_ptr kernel, int sign) + { + if (sign == -1) + { + return std::make_shared(kernel, sign); + } + else + { + return std::make_shared>(kernel, sign); + } + } + + inline void get_symmetrized(AbstractReducedKernel &kernel, int sign) + { + throw std::runtime_error("cannot symmetrize twice"); + } + } // namespace sparseir namespace sparseir{ @@ -950,6 +1051,31 @@ namespace sparseir{ return res; } + class SVEHintsReduced : public AbstractSVEHints + { + public: + SVEHintsReduced(std::shared_ptr inner_hints) + : inner(inner_hints) {} + + // Implement required methods + int nsvals() const override + { + // Implement this function + // For example, you can delegate the call to the inner object + return inner->nsvals(); + } + + int ngauss() const override + { + // Implement this function + // For example, you can delegate the call to the inner object + return inner->ngauss(); + } + + private: + std::shared_ptr inner; + }; + // Function to provide SVE hints inline SVEHintsLogistic sve_hints(const LogisticKernel& kernel, double epsilon) { return SVEHintsLogistic(kernel, epsilon); @@ -959,6 +1085,30 @@ namespace sparseir{ return SVEHintsRegularizedBose(kernel, epsilon); } + inline std::shared_ptr sve_hints(std::shared_ptr kernel, double epsilon) + { + if (auto logisticKernel = std::dynamic_pointer_cast(kernel)) + { + return std::make_shared(*logisticKernel, epsilon); + } + else if (auto boseKernel = std::dynamic_pointer_cast(kernel)) + { + return std::make_shared(*boseKernel, epsilon); + } + else if (auto reducedKernel = std::dynamic_pointer_cast(kernel)) + { + return std::make_shared(sve_hints(reducedKernel->inner, epsilon)); + } + else + { + throw std::invalid_argument("Unsupported kernel type for SVE hints"); + } + } + + inline SVEHintsReduced sve_hints(const AbstractReducedKernel& kernel, double epsilon) { + return SVEHintsReduced(sve_hints(kernel.inner, epsilon)); + } + /* function ngauss end ngauss(hints::SVEHintsLogistic) = hints.ε ≥ sqrt(eps()) ? 10 : 16 diff --git a/include/sparseir/sparseir-header-only.hpp b/include/sparseir/sparseir-header-only.hpp index 31e5193..9b21cfc 100644 --- a/include/sparseir/sparseir-header-only.hpp +++ b/include/sparseir/sparseir-header-only.hpp @@ -1,12 +1,10 @@ #pragma once +#include "./_linalg.hpp" #include "./_specfuncs.hpp" - #include "./_root.hpp" -#include "./_linalg.hpp" - -#include "./gauss.hpp" #include "./freq.hpp" +#include "./svd.hpp" +#include "./gauss.hpp" #include "./poly.hpp" -#include "./kernel.hpp" -#include "./svd.hpp" \ No newline at end of file +#include "./kernel.hpp" \ No newline at end of file diff --git a/test/kernel.cxx b/test/kernel.cxx index d1790db..8655f87 100644 --- a/test/kernel.cxx +++ b/test/kernel.cxx @@ -13,7 +13,7 @@ using xprec::DDouble; template -bool kernel_accuracy_test(Kernel &K) { +std::tuple kernel_accuracy_test(Kernel &K) { using T = float; using T_x = double; @@ -51,13 +51,13 @@ bool kernel_accuracy_test(Kernel &K) { 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); + bool b1 = (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; + bool b2 = (reldiff - T(1.0)).cwiseAbs().maxCoeff() <= 100 * epsilon; + return std::make_tuple(b1, b2); } TEST_CASE("Kernel Accuracy Test") @@ -70,15 +70,37 @@ TEST_CASE("Kernel Accuracy Test") 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)); + bool b1, b2; + std::tie(b1, b2) = kernel_accuracy_test(K); + REQUIRE(b1); + REQUIRE(b2); + /* + if (b1){ + std::cout << "Kernel accuracy test passed for " << typeid(K).name() << b1 << b2 << std::endl; + } + + if (b2){ + std::cout << "Kernel accuracy test passed for " << typeid(K).name() << b1 << b2 << std::endl; + } + */ } } - + { + auto kernel_ptr = std::make_shared(40000.0); + auto K = sparseir::get_symmetrized(kernel_ptr, -1); + // TODO: implement sve_hints + bool b1, b2; + //std::tie(b1, b2) = kernel_accuracy_test(K); + //REQUIRE(b1); + //REQUIRE(b2); + // TODO: resolve this errors + //auto k2 = sparseir::get_symmetrized(sparseir::RegularizedBoseKernel(40000.0), -1); + // TODO: implement sve_hints + // REQUIRE(kernel_accuracy_test(k2)); + } { // List of kernels to test std::vector kernels = { @@ -87,7 +109,23 @@ TEST_CASE("Kernel Accuracy Test") }; for (const auto &K : kernels) { - REQUIRE(kernel_accuracy_test(K)); + bool b1, b2; + std::tie(b1, b2) = kernel_accuracy_test(K); + // TODO: resolve this errors + REQUIRE(b1); + // TODO: resolve this errors + REQUIRE(b2); + /* + if (b1) + { + std::cout << "Kernel accuracy test passed for " << typeid(K).name() << b1 << b2 << std::endl; + } + + if (b2) + { + std::cout << "Kernel accuracy test passed for " << typeid(K).name() << b1 << b2 << std::endl; + } + */ } } /*