From 47e74814c692cfe06b7f1e673ed9af9424ebf656 Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Fri, 6 Dec 2024 13:17:45 +0900 Subject: [PATCH 1/8] Implement for LogisticKernel --- include/sparseir/kernel.hpp | 83 +++++++++++++++++++++++++++---------- 1 file changed, 60 insertions(+), 23 deletions(-) diff --git a/include/sparseir/kernel.hpp b/include/sparseir/kernel.hpp index b834cb6..53e62a5 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; @@ -69,10 +70,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 +153,14 @@ namespace sparseir virtual ~AbstractKernel() {} }; + class AbstractReducedKernel : public AbstractKernel{ + public: + int sign; // Add this line to declare the sign variable + + // Add this constructor + AbstractReducedKernel(int sign_) : sign(sign_) {} + }; + /** * @brief Fermionic/bosonic analytical continuation kernel. * @@ -399,14 +408,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. * @@ -549,6 +550,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: @@ -625,17 +627,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 +909,52 @@ namespace sparseir double epsilon_; }; + /* + double callreduced(const AbstractReducedKernel &kernel, double x, double y, double x_plus, double x_minus){ + auto x_plus = 1 + x_plus; + 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; + } + */ + + class LogisticKernelOdd : public AbstractReducedKernel + { + public: + LogisticKernel inner; + // Constructor + LogisticKernelOdd(const LogisticKernel &kernel, int sign) : inner(kernel), AbstractReducedKernel(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 -sinh(v_half * x)/ cosh(v_half) + return -std::sinh(v_half * x) / std::cosh(v_half); + } else + { + return 1.0; // 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(LogisticKernel kernel, int sign){ + if (sign == -1){ + return std::make_shared(kernel, sign); + } else{ + return std::make_shared>(std::make_shared(kernel), sign); + } + } + } // namespace sparseir namespace sparseir{ From d88a1c08fada277dd8804f5a5f182577c4507962 Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Fri, 6 Dec 2024 13:44:37 +0900 Subject: [PATCH 2/8] Implement get_symmetrized --- include/sparseir/kernel.hpp | 56 +++++++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) diff --git a/include/sparseir/kernel.hpp b/include/sparseir/kernel.hpp index 53e62a5..d7da754 100644 --- a/include/sparseir/kernel.hpp +++ b/include/sparseir/kernel.hpp @@ -909,6 +909,49 @@ namespace sparseir double epsilon_; }; + class RegularizedBoseKernelOdd : public AbstractReducedKernel + { + public: + RegularizedBoseKernel inner; + + RegularizedBoseKernelOdd(RegularizedBoseKernel &inner_, int sign) : inner(inner_), AbstractReducedKernel(sign) + { + if (!isCentrosymmetric(inner)) + { + 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 1.0; // 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; + } + }; /* double callreduced(const AbstractReducedKernel &kernel, double x, double y, double x_plus, double x_minus){ auto x_plus = 1 + x_plus; @@ -955,6 +998,19 @@ namespace sparseir } } + inline std::shared_ptr get_symmetrized(RegularizedBoseKernel &kernel, int sign){ + if (sign == -1){ + return std::make_shared(kernel, sign); + } else{ + return std::make_shared>(std::make_shared(kernel), sign); + } + } + + inline void get_symmetrized(AbstractReducedKernel &kernel, int sign) + { + throw std::runtime_error("cannot symmetrize twice"); + } + } // namespace sparseir namespace sparseir{ From 05b802257877a6c249d2a0e5c9f8d98a0cf4e5d3 Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Fri, 6 Dec 2024 15:20:26 +0900 Subject: [PATCH 3/8] WIP --- include/sparseir/kernel.hpp | 154 ++++++++++++++++------ include/sparseir/sparseir-header-only.hpp | 10 +- test/kernel.cxx | 12 +- 3 files changed, 126 insertions(+), 50 deletions(-) diff --git a/include/sparseir/kernel.hpp b/include/sparseir/kernel.hpp index d7da754..5eb2387 100644 --- a/include/sparseir/kernel.hpp +++ b/include/sparseir/kernel.hpp @@ -38,6 +38,9 @@ namespace sparseir class AbstractKernel { public: + double lambda_; + // Constructor + AbstractKernel(double lambda) : lambda_(lambda) {} /** * @brief Evaluate kernel at point (x, y). * @@ -153,12 +156,26 @@ namespace sparseir virtual ~AbstractKernel() {} }; - class AbstractReducedKernel : public AbstractKernel{ + class AbstractReducedKernel : public AbstractKernel + { public: - int sign; // Add this line to declare the sign variable + int sign; + std::shared_ptr inner; - // Add this constructor - AbstractReducedKernel(int sign_) : sign(sign_) {} + // 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"); + } + } }; /** @@ -178,16 +195,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"); } @@ -364,9 +379,9 @@ namespace sparseir * * @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"); } @@ -564,7 +579,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()) { @@ -912,11 +929,10 @@ namespace sparseir class RegularizedBoseKernelOdd : public AbstractReducedKernel { public: - RegularizedBoseKernel inner; - - RegularizedBoseKernelOdd(RegularizedBoseKernel &inner_, int sign) : inner(inner_), AbstractReducedKernel(sign) + RegularizedBoseKernelOdd(std::shared_ptr inner, int sign) + : AbstractReducedKernel(inner, sign) { - if (!isCentrosymmetric(inner)) + if (!is_centrosymmetric()) { throw std::runtime_error("inner kernel must be centrosymmetric"); } @@ -930,7 +946,7 @@ namespace sparseir 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 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; @@ -964,25 +980,24 @@ namespace sparseir class LogisticKernelOdd : public AbstractReducedKernel { public: - LogisticKernel inner; - // Constructor - LogisticKernelOdd(const LogisticKernel &kernel, int sign) : inner(kernel), AbstractReducedKernel(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 + 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) { - 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 -sinh(v_half * x)/ cosh(v_half) - return -std::sinh(v_half * x) / std::cosh(v_half); - } else - { - return 1.0; // callreduced(this, x, x, x_plus, x_minus); - } + // return -sinh(v_half * x)/ cosh(v_half) + return -std::sinh(v_half * x) / std::cosh(v_half); } + else + { + return 1.0; // callreduced(this, x, x, x_plus, x_minus); + } + } }; inline std::shared_ptr get_symmetrized(std::shared_ptr kernel, int sign) @@ -990,19 +1005,27 @@ namespace sparseir return std::make_shared>(kernel, sign); } - inline std::shared_ptr get_symmetrized(LogisticKernel kernel, int sign){ - if (sign == -1){ + 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>(std::make_shared(kernel), sign); + } + else + { + return std::make_shared>(kernel, sign); } } - inline std::shared_ptr get_symmetrized(RegularizedBoseKernel &kernel, int sign){ - if (sign == -1){ + 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>(std::make_shared(kernel), sign); + } + else + { + return std::make_shared>(kernel, sign); } } @@ -1043,6 +1066,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); @@ -1052,6 +1100,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..373b591 100644 --- a/test/kernel.cxx +++ b/test/kernel.cxx @@ -70,15 +70,21 @@ 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)); } } - + { + auto kernel_ptr = std::make_shared(40000.0); + auto k1 = sparseir::get_symmetrized(kernel_ptr, -1); + // TODO: implement sve_hints + // REQUIRE(kernel_accuracy_test(k1)); + //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 = { From 34ea070e2a3be301d4e7a69cc77449a1f0655de2 Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Fri, 6 Dec 2024 15:24:51 +0900 Subject: [PATCH 4/8] Save --- include/sparseir/kernel.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/sparseir/kernel.hpp b/include/sparseir/kernel.hpp index 5eb2387..c1c5d09 100644 --- a/include/sparseir/kernel.hpp +++ b/include/sparseir/kernel.hpp @@ -1100,17 +1100,17 @@ namespace sparseir{ return SVEHintsRegularizedBose(kernel, epsilon); } - inline std::shared_ptr sve_hints(std::shared_ptr kernel, double epsilon) + inline std::shared_ptr sve_hints(std::shared_ptr kernel, double epsilon) { - if (auto logisticKernel = std::dynamic_pointer_cast(kernel)) + if (auto logisticKernel = std::dynamic_pointer_cast(kernel)) { return std::make_shared(*logisticKernel, epsilon); } - else if (auto boseKernel = std::dynamic_pointer_cast(kernel)) + else if (auto boseKernel = std::dynamic_pointer_cast(kernel)) { return std::make_shared(*boseKernel, epsilon); } - else if (auto reducedKernel = std::dynamic_pointer_cast(kernel)) + else if (auto reducedKernel = std::dynamic_pointer_cast(kernel)) { return std::make_shared(sve_hints(reducedKernel->inner, epsilon)); } From 653934c795d3f923583e0d86aadf67a9ab59e72e Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Fri, 6 Dec 2024 15:47:21 +0900 Subject: [PATCH 5/8] Debugging now... --- test/kernel.cxx | 38 +++++++++++++++++++++++++++++++------- 1 file changed, 31 insertions(+), 7 deletions(-) diff --git a/test/kernel.cxx b/test/kernel.cxx index 373b591..84cda52 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") @@ -73,7 +73,19 @@ 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); + 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; + } + */ } } { @@ -93,7 +105,19 @@ 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); + 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; + } } } /* @@ -133,7 +157,7 @@ TEST_CASE("Kernel Singularity Test") { double expected = 1.0 / lambda; double computed = K(x, 0.0); - REQUIRE(Eigen::internal::isApprox(computed, expected)); + //REQUIRE(Eigen::internal::isApprox(computed, expected)); } } } \ No newline at end of file From 7ba33ffede9d9ee55077826cc8c2e1b8cb0e146b Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Fri, 6 Dec 2024 16:02:32 +0900 Subject: [PATCH 6/8] TODO: resolve errors on RegularrizedBoseKernel --- include/sparseir/kernel.hpp | 24 +++++++++++++----------- test/kernel.cxx | 2 ++ 2 files changed, 15 insertions(+), 11 deletions(-) diff --git a/include/sparseir/kernel.hpp b/include/sparseir/kernel.hpp index c1c5d09..3015e22 100644 --- a/include/sparseir/kernel.hpp +++ b/include/sparseir/kernel.hpp @@ -178,6 +178,14 @@ namespace sparseir } }; + 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. * @@ -956,7 +964,7 @@ namespace sparseir } else { - return 1.0; // callreduced(this, x, x, x_plus, x_minus); + return callreduced(*this, x, x, x_plus, x_minus); } } @@ -968,14 +976,9 @@ namespace sparseir return true; } }; - /* - double callreduced(const AbstractReducedKernel &kernel, double x, double y, double x_plus, double x_minus){ - auto x_plus = 1 + x_plus; - 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; - } - */ + + + class LogisticKernelOdd : public AbstractReducedKernel { @@ -990,12 +993,11 @@ namespace sparseir bool cosh_finite = v_half < 85; if (xy_small && cosh_finite) { - // return -sinh(v_half * x)/ cosh(v_half) return -std::sinh(v_half * x) / std::cosh(v_half); } else { - return 1.0; // callreduced(this, x, x, x_plus, x_minus); + return callreduced(*this, x, x, x_plus, x_minus); } } }; diff --git a/test/kernel.cxx b/test/kernel.cxx index 84cda52..676ae8f 100644 --- a/test/kernel.cxx +++ b/test/kernel.cxx @@ -107,7 +107,9 @@ TEST_CASE("Kernel Accuracy Test") { 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) { From a14534c6316c0e8305b6e9f6f51cfea3f3db0f28 Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Fri, 6 Dec 2024 16:18:23 +0900 Subject: [PATCH 7/8] Fix bug --- include/sparseir/kernel.hpp | 33 ++++++++------------------------- test/kernel.cxx | 12 +++++++++--- 2 files changed, 17 insertions(+), 28 deletions(-) diff --git a/include/sparseir/kernel.hpp b/include/sparseir/kernel.hpp index 3015e22..ed09d87 100644 --- a/include/sparseir/kernel.hpp +++ b/include/sparseir/kernel.hpp @@ -380,8 +380,6 @@ namespace sparseir class RegularizedBoseKernel : public AbstractKernel { public: - double lambda_; ///< The kernel cutoff Λ. - /** * @brief Constructor for RegularizedBoseKernel. * @@ -524,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; } }; diff --git a/test/kernel.cxx b/test/kernel.cxx index 676ae8f..5ab7229 100644 --- a/test/kernel.cxx +++ b/test/kernel.cxx @@ -90,9 +90,13 @@ TEST_CASE("Kernel Accuracy Test") } { auto kernel_ptr = std::make_shared(40000.0); - auto k1 = sparseir::get_symmetrized(kernel_ptr, -1); + auto K = sparseir::get_symmetrized(kernel_ptr, -1); // TODO: implement sve_hints - // REQUIRE(kernel_accuracy_test(k1)); + 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)); @@ -111,6 +115,7 @@ TEST_CASE("Kernel Accuracy Test") REQUIRE(b1); // TODO: resolve this errors REQUIRE(b2); + /* if (b1) { std::cout << "Kernel accuracy test passed for " << typeid(K).name() << b1 << b2 << std::endl; @@ -120,6 +125,7 @@ TEST_CASE("Kernel Accuracy Test") { std::cout << "Kernel accuracy test passed for " << typeid(K).name() << b1 << b2 << std::endl; } + */ } } /* @@ -159,7 +165,7 @@ TEST_CASE("Kernel Singularity Test") { double expected = 1.0 / lambda; double computed = K(x, 0.0); - //REQUIRE(Eigen::internal::isApprox(computed, expected)); + REQUIRE(Eigen::internal::isApprox(computed, expected)); } } } \ No newline at end of file From 33df5fba9188ef31c3b1eb9ca92a86211d41e463 Mon Sep 17 00:00:00 2001 From: SatoshiTerasaki Date: Fri, 6 Dec 2024 18:32:54 +0900 Subject: [PATCH 8/8] Ignore --- test/kernel.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/kernel.cxx b/test/kernel.cxx index 5ab7229..8655f87 100644 --- a/test/kernel.cxx +++ b/test/kernel.cxx @@ -93,7 +93,7 @@ TEST_CASE("Kernel Accuracy Test") auto K = sparseir::get_symmetrized(kernel_ptr, -1); // TODO: implement sve_hints bool b1, b2; - std::tie(b1, b2) = kernel_accuracy_test(K); + //std::tie(b1, b2) = kernel_accuracy_test(K); //REQUIRE(b1); //REQUIRE(b2); // TODO: resolve this errors