Merge pull request #90 from SpM-lab/terasaki/port-sampling-points

Define default_tau_sampling_points as a method
SpM-lab · Jan 30, 2025 · 4f7aebd · 4f7aebd
2 parents a3ff815 + 3937bca
commit 4f7aebd
Show file tree

Hide file tree

Showing 2 changed files with 40 additions and 47 deletions.
diff --git a/include/sparseir/augment.hpp b/include/sparseir/augment.hpp
@@ -300,6 +300,12 @@ class AugmentedBasis : public AbstractBasis<S> {
     const AugmentedMatsubaraFunction<S>& uhat() const { return *uhat_; }
     size_t nAug() const { return augmentations_.size(); }
 
+    const Eigen::VectorXd default_tau_sampling_points() const override {
+        int sz = this->basis_->sve_result->s.size() + this->augmentations_.size();
+        auto x = default_sampling_points(this->basis_->sve_result->u, sz);
+        return (this->basis_->get_beta() / 2.0) * (x.array() + 1.0);
+    }
+
     // Factory method
     static std::shared_ptr<AugmentedBasis<S>> create(
         std::shared_ptr<FiniteTempBasis<S>> basis,
@@ -308,11 +314,4 @@ class AugmentedBasis : public AbstractBasis<S> {
     }
 };
 
-template <typename S>
-inline Eigen::VectorXd default_tau_sampling_points(AugmentedBasis<S> basis){
-    int sz = basis.basis.sve_result.s.size() + basis.augmentations.size();
-    auto x = default_samplint_points(basis.basis.sve_result.u, sz);
-    return (basis.beta / 2.0) * (x.array() + 1.0);
-}
-
 } // namespace sparseir
diff --git a/include/sparseir/basis.hpp b/include/sparseir/basis.hpp
@@ -42,6 +42,8 @@ class AbstractBasis {
     virtual double get_accuracy() const = 0;
     virtual double get_wmax() const = 0;
     virtual size_t size() const = 0;
+    virtual const Eigen::VectorXd default_tau_sampling_points() const = 0;
+
 };
 
 } // namespace sparseir
@@ -157,14 +159,6 @@ class FiniteTempBasis : public AbstractBasis<S> {
     // Getter for Λ
     //double Lambda() const { return lambda; }
 
-    // Default τ sampling points
-    Eigen::VectorXd defaultTauSamplingPoints() const
-    {
-        Eigen::VectorXd x =
-            default_sampling_points(sve_result->u, static_cast<int>(s.size()));
-        return (this->get_beta() / 2.0) * (x.array() + 1.0);
-    }
-
     // Default Matsubara sampling points
     Eigen::VectorXd
     defaultMatsubaraSamplingPoints(bool positive_only = false) const
@@ -190,6 +184,12 @@ class FiniteTempBasis : public AbstractBasis<S> {
             lambda, *sve_result, static_cast<int>(s.size()));
     }
 
+    const Eigen::VectorXd default_tau_sampling_points() const override {
+        int sz = this->sve_result->s.size();
+        auto x = default_sampling_points(this->sve_result->u, sz);
+        return (this->beta / 2.0) * (x.array() + 1.0);
+    }
+
 private:
     // Placeholder statistics function
     std::shared_ptr<S> statistics() const { return std::make_shared<S>(); }
@@ -249,40 +249,33 @@ class FiniteTempBasis : public AbstractBasis<S> {
 
     // Default sampling points function
 inline Eigen::VectorXd default_sampling_points(const PiecewiseLegendrePolyVector &u, int L) {
-        if (u.xmin() != -1.0 || u.xmax() != 1.0)
-            throw std::runtime_error("Expecting unscaled functions here.");
-
-        if (L < u.size()) {
-            // TODO: Resolve this errors.
-            return u.polyvec[L].roots();
-        } else {
-            // Approximate roots by extrema
-            // TODO: resolve this error
-            PiecewiseLegendrePoly poly = u.polyvec.back();
-            Eigen::VectorXd maxima = poly.deriv().roots();
-
-            double left = (maxima[0] + poly.xmin) / 2.0;
-            double right = (maxima[maxima.size() - 1] + poly.xmax) / 2.0;
-
-            Eigen::VectorXd x0(maxima.size() + 2);
-            x0[0] = left;
-            x0.tail(maxima.size()) = maxima;
-            x0[x0.size() - 1] = right;
-
-            if (x0.size() != L) {
-                std::cerr << "Warning: Expected " << L
-                          << " sampling points, got " << x0.size() << ".\n";
-            }
-
-            return x0;
+    if (u.xmin() != -1.0 || u.xmax() != 1.0)
+        throw std::runtime_error("Expecting unscaled functions here.");
+
+    if (L < u.size()) {
+        // TODO: Resolve this errors.
+        return u.polyvec[L].roots();
+    } else {
+        // Approximate roots by extrema
+        // TODO: resolve this error
+        PiecewiseLegendrePoly poly = u.polyvec.back();
+        Eigen::VectorXd maxima = poly.deriv().roots();
+
+        double left = (maxima[0] + poly.xmin) / 2.0;
+        double right = (maxima[maxima.size() - 1] + poly.xmax) / 2.0;
+
+        Eigen::VectorXd x0(maxima.size() + 2);
+        x0[0] = left;
+        x0.segment(1, maxima.size()) = maxima;
+        x0[x0.size() - 1] = right;
+
+        if (x0.size() != L) {
+            std::cerr << "Warning: Expected " << L
+                        << " sampling points, got " << x0.size() << ".\n";
         }
-    }
 
-template <typename S>
-inline Eigen::VectorXd default_tau_sampling_points(std::shared_ptr<FiniteTempBasis<S>> basis){
-    int sz = basis.sve_result.s.size();
-    auto x = default_samplint_points(basis.sve_result.u, sz);
-    return (basis.beta / 2.0) * (x.array() + 1.0);
+        return x0;
+    }
 }
 
 inline std::pair<FiniteTempBasis<Fermionic>,
@@ -316,4 +309,5 @@ inline std::pair<FiniteTempBasis<Fermionic>,
         beta, omega_max, epsilon, kernel, sve_result);
     return std::make_pair(basis_f, basis_b);
 }
+
 } // namespace sparseir