Merge pull request #73 from SpM-lab/terasaki/default-samplint-points

Write tests for default_sampling_points
SpM-lab · Dec 24, 2024 · b1d2839 · b1d2839
2 parents 9a43fb9 + 7cf2429
commit b1d2839
Show file tree

Hide file tree

Showing 2 changed files with 45 additions and 34 deletions.
diff --git a/include/sparseir/basis.hpp b/include/sparseir/basis.hpp
@@ -224,7 +224,7 @@ class FiniteTempBasis : public AbstractBasis<S> {
 
     // Delegating constructor 1
     FiniteTempBasis(double beta, double omega_max,
-                    double epsilon, int max_size)
+                    double epsilon, int max_size=-1)
         : FiniteTempBasis(beta, omega_max, epsilon,
                         LogisticKernel(beta * omega_max),
                         compute_sve<LogisticKernel>(LogisticKernel(beta * omega_max), epsilon),
@@ -305,39 +305,6 @@ class FiniteTempBasis : public AbstractBasis<S> {
     // Placeholder statistics function
     std::shared_ptr<S> statistics() const { return std::make_shared<S>(); }
 
-    // Default sampling points function
-    Eigen::VectorXd
-    default_sampling_points(const PiecewiseLegendrePolyVector &u, int L) const
-    {
-        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;
-        }
-    }
-
     // Default Matsubara sampling points function
     Eigen::VectorXd defaultMatsubaraSamplingPoints(
         const PiecewiseLegendreFTVector<S> &u_hat_full, int L,
@@ -391,6 +358,37 @@ 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;
+        }
+    }
+
 inline std::pair<FiniteTempBasis<Fermionic, LogisticKernel>,
                  FiniteTempBasis<Bosonic, LogisticKernel>>
     finite_temp_bases(

diff --git a/test/basis.cxx b/test/basis.cxx
@@ -99,4 +99,17 @@ TEST_CASE("FiniteTempBasis consistency tests", "[basis]") {
         REQUIRE(rescaled_basis.sve_result->s.size() == basis.sve_result->s.size());
         REQUIRE(rescaled_basis.get_wmax() == 6.0);
     }
+
+    SECTION("default_sampling_points") {
+        auto basis = sparseir::FiniteTempBasis<sparseir::Fermionic, sparseir::LogisticKernel>(3., 4., 1e-6);
+        auto sve = basis.sve_result;
+        auto s = sve->s;
+        //std::cout << "Singular values: " << s.transpose() << std::endl;
+        int L = 10;
+        Eigen::VectorXd pts_L = default_sampling_points(sve->u, L);
+        REQUIRE(pts_L.size() == L);
+        Eigen::VectorXd pts_100 = default_sampling_points(sve->u, 100);
+        REQUIRE(pts_100.size() == 24);
+    }
+
 }