From 84e7cf27cdae94de300e9572b8c661d00e45ed34 Mon Sep 17 00:00:00 2001 From: SatoshiTerasaki Date: Tue, 24 Dec 2024 20:09:38 +0900 Subject: [PATCH] Update --- include/sparseir/augment.hpp | 7 ++++ include/sparseir/basis.hpp | 7 ++++ include/sparseir/sampling.hpp | 77 ++++++++++++++++------------------- test/sampling.cxx | 2 +- test/sve.cxx | 7 ++++ 5 files changed, 56 insertions(+), 44 deletions(-) diff --git a/include/sparseir/augment.hpp b/include/sparseir/augment.hpp index b5f807a..5625137 100644 --- a/include/sparseir/augment.hpp +++ b/include/sparseir/augment.hpp @@ -231,4 +231,11 @@ class AugmentedBasis : public AbstractBasis { }; +template +inline Eigen::VectorXd default_tau_sampling_points(AugmentedBasis 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 index 908ed56..6f6a04e 100644 --- a/include/sparseir/basis.hpp +++ b/include/sparseir/basis.hpp @@ -389,6 +389,13 @@ inline Eigen::VectorXd default_sampling_points(const PiecewiseLegendrePolyVector } } +template +inline Eigen::VectorXd default_tau_sampling_points(std::shared_ptr> 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); +} + inline std::pair, FiniteTempBasis> finite_temp_bases( diff --git a/include/sparseir/sampling.hpp b/include/sparseir/sampling.hpp index 92ba910..a493c52 100644 --- a/include/sparseir/sampling.hpp +++ b/include/sparseir/sampling.hpp @@ -7,19 +7,19 @@ namespace sparseir { -template +template class AbstractSampling { public: virtual ~AbstractSampling() = default; // Evaluate the basis coefficients at sampling points - virtual Eigen::Matrix evaluate( - const Eigen::Matrix& al, + virtual Eigen::Matrix evaluate( + const Eigen::Matrix& al, const Eigen::VectorXd* points = nullptr) const = 0; // Fit values at sampling points to basis coefficients - virtual Eigen::Matrix fit( - const Eigen::Matrix& ax, + virtual Eigen::Matrix fit( + const Eigen::Matrix& ax, const Eigen::VectorXd* points = nullptr) const = 0; // Condition number of the transformation matrix @@ -29,34 +29,34 @@ class AbstractSampling { virtual const Eigen::VectorXd& sampling_points() const = 0; // Get the transformation matrix - virtual const Eigen::Matrix& matrix() const = 0; + virtual const Eigen::Matrix& matrix() const = 0; // Get the associated basis - virtual const std::shared_ptr>& basis() const = 0; + //virtual const std::shared_ptr>& basis() const = 0; protected: // Create a new sampling object for different sampling points - virtual std::shared_ptr> for_sampling_points( + virtual std::shared_ptr for_sampling_points( const Eigen::VectorXd& x) const = 0; }; -template -class TauSampling : public AbstractSampling { +template +class TauSampling : public AbstractSampling { public: TauSampling( - const std::shared_ptr>& basis, + const std::shared_ptr>& basis, const Eigen::VectorXd* sampling_points = nullptr) : basis_(basis) { if (sampling_points) { sampling_points_ = *sampling_points; } else { - sampling_points_ = basis_->default_tau_sampling_points(); + sampling_points_ = default_tau_sampling_points(basis_); } construct_matrix(); } - Eigen::Matrix evaluate( - const Eigen::Matrix& al, + Eigen::Matrix evaluate( + const Eigen::Matrix& al, const Eigen::VectorXd* points = nullptr) const override { if (points) { auto sampling = for_sampling_points(*points); @@ -65,8 +65,9 @@ class TauSampling : public AbstractSampling { return matrix_ * al; } - Eigen::Matrix fit( - const Eigen::Matrix& ax, + /* + Eigen::Matrix fit( + const Eigen::Matrix& ax, const Eigen::VectorXd* points = nullptr) const override { if (points) { auto sampling = for_sampling_points(*points); @@ -74,6 +75,7 @@ class TauSampling : public AbstractSampling { } return solve(ax); } + */ double cond() const override { return cond_; @@ -83,18 +85,19 @@ class TauSampling : public AbstractSampling { return sampling_points_; } - const Eigen::Matrix& matrix() const override { + const Eigen::Matrix& matrix() const override { return matrix_; } - - const std::shared_ptr>& basis() const override { + /* + const std::shared_ptr>& basis() const override { return basis_; } + */ protected: - std::shared_ptr> for_sampling_points( + std::shared_ptr> for_sampling_points( const Eigen::VectorXd& x) const override { - return std::make_shared>(basis_, &x); + return std::make_shared(basis_, &x); } private: @@ -107,24 +110,24 @@ class TauSampling : public AbstractSampling { } } - Eigen::Matrix solve( - const Eigen::Matrix& ax) const { + Eigen::Matrix solve( + const Eigen::Matrix& ax) const { return solver_.solve(ax); } double compute_condition_number( - const Eigen::Matrix& mat) const { - Eigen::JacobiSVD> svd( + const Eigen::Matrix& mat) const { + Eigen::JacobiSVD> svd( mat, Eigen::ComputeThinU | Eigen::ComputeThinV); double cond = svd.singularValues()(0) / svd.singularValues()(svd.singularValues().size() - 1); return cond; } - std::shared_ptr> basis_; + std::shared_ptr> basis_; Eigen::VectorXd sampling_points_; - Eigen::Matrix matrix_; - mutable Eigen::ColPivHouseholderQR> solver_; + Eigen::Matrix matrix_; + mutable Eigen::ColPivHouseholderQR> solver_; double cond_; }; @@ -144,25 +147,19 @@ class MatsubaraSampling : public AbstractSampling> { construct_matrix(); } + /* Eigen::Matrix, Eigen::Dynamic, 1> evaluate( const Eigen::Matrix& al, const Eigen::VectorXi* points = nullptr) const override { if (points) { + // TODO: Implement for_sampling_points auto sampling = for_sampling_points(*points); return sampling->matrix() * al; } return matrix_ * al; } + */ - Eigen::Matrix fit( - const Eigen::Matrix, Eigen::Dynamic, 1>& ax, - const Eigen::VectorXi* points = nullptr) const override { - if (points) { - auto sampling = for_sampling_points(*points); - return sampling->solve(ax); - } - return solve(ax); - } double cond() const override { return cond_; @@ -180,12 +177,6 @@ class MatsubaraSampling : public AbstractSampling> { return basis_; } -protected: - std::shared_ptr>> for_sampling_points( - const Eigen::VectorXi& n) const override { - return std::make_shared>(basis_, positive_only_, &n); - } - private: void construct_matrix() { matrix_ = basis_->uhat(sampling_points_); diff --git a/test/sampling.cxx b/test/sampling.cxx index 5827487..5cfd5b6 100644 --- a/test/sampling.cxx +++ b/test/sampling.cxx @@ -17,7 +17,7 @@ TEST_CASE("sampling", "[Sampling]") { auto basis = std::make_shared>( beta, omega_max, std::numeric_limits::quiet_NaN(), sparseir::LogisticKernel(beta * omega_max), sve_result); - //sparseir::TauSampling tau_sampling(basis); + // sparseir::TauSampling TauSampling(basis); // Here we can check the type or properties of tau_sampling if needed REQUIRE(true); // Placeholder } diff --git a/test/sve.cxx b/test/sve.cxx index 1072d3b..76d3c3e 100644 --- a/test/sve.cxx +++ b/test/sve.cxx @@ -113,6 +113,13 @@ TEST_CASE("CentrosymmSVE", "[CentrosymmSVE]") 1e-6); } +TEST_CASE("compute_sve", "[compute_sve]"){ + sparseir::LogisticKernel lk(12.0); + auto sve = sparseir::compute_sve(lk); + auto s = sve.s; + //std::cout << "S values: \n" << s << std::endl; +} + TEST_CASE("sve.cpp", "[compute_sve]") {