Skip to content

Commit

Permalink
Merge pull request #74 from SpM-lab/terasaki/default-tau-samplint-points
Browse files Browse the repository at this point in the history 
Implement default-tau-samplint-points
  • Loading branch information
@terasakisatoshi
terasakisatoshi authored Dec 24, 2024
2 parents b1d2839 + 84e7cf2 commit a089c33
Show file tree
Hide file tree
Showing 5 changed files with 56 additions and 44 deletions.
7 changes: 7 additions & 0 deletions include/sparseir/augment.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -231,4 +231,11 @@ class AugmentedBasis : public AbstractBasis<S> {
};


template <typename S, typename B, typename F, typename FHAT>
inline Eigen::VectorXd default_tau_sampling_points(AugmentedBasis<S, B, F, FHAT> 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
7 changes: 7 additions & 0 deletions include/sparseir/basis.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -389,6 +389,13 @@ inline Eigen::VectorXd default_sampling_points(const PiecewiseLegendrePolyVector
}
}

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);
}

inline std::pair<FiniteTempBasis<Fermionic, LogisticKernel>,
FiniteTempBasis<Bosonic, LogisticKernel>>
finite_temp_bases(
Expand Down
77 changes: 34 additions & 43 deletions include/sparseir/sampling.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,19 +7,19 @@

namespace sparseir {

template <typename T>
template <typename S>
class AbstractSampling {
public:
virtual ~AbstractSampling() = default;

// Evaluate the basis coefficients at sampling points
virtual Eigen::Matrix<T, Eigen::Dynamic, 1> evaluate(
const Eigen::Matrix<T, Eigen::Dynamic, 1>& al,
virtual Eigen::Matrix<double, Eigen::Dynamic, 1> evaluate(
const Eigen::Matrix<double, Eigen::Dynamic, 1>& al,
const Eigen::VectorXd* points = nullptr) const = 0;

// Fit values at sampling points to basis coefficients
virtual Eigen::Matrix<T, Eigen::Dynamic, 1> fit(
const Eigen::Matrix<T, Eigen::Dynamic, 1>& ax,
virtual Eigen::Matrix<double, Eigen::Dynamic, 1> fit(
const Eigen::Matrix<double, Eigen::Dynamic, 1>& ax,
const Eigen::VectorXd* points = nullptr) const = 0;

// Condition number of the transformation matrix
Expand All @@ -29,34 +29,34 @@ class AbstractSampling {
virtual const Eigen::VectorXd& sampling_points() const = 0;

// Get the transformation matrix
virtual const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& matrix() const = 0;
virtual const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& matrix() const = 0;

// Get the associated basis
virtual const std::shared_ptr<AbstractBasis<T>>& basis() const = 0;
//virtual const std::shared_ptr<AbstractBasis<S>>& basis() const = 0;

protected:
// Create a new sampling object for different sampling points
virtual std::shared_ptr<AbstractSampling<T>> for_sampling_points(
virtual std::shared_ptr<AbstractSampling> for_sampling_points(
const Eigen::VectorXd& x) const = 0;
};

template <typename T>
class TauSampling : public AbstractSampling<T> {
template <typename S>
class TauSampling : public AbstractSampling<S> {
public:
TauSampling(
const std::shared_ptr<AbstractBasis<T>>& basis,
const std::shared_ptr<AbstractBasis<S>>& 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<T, Eigen::Dynamic, 1> evaluate(
const Eigen::Matrix<T, Eigen::Dynamic, 1>& al,
Eigen::Matrix<double, Eigen::Dynamic, 1> evaluate(
const Eigen::Matrix<double, Eigen::Dynamic, 1>& al,
const Eigen::VectorXd* points = nullptr) const override {
if (points) {
auto sampling = for_sampling_points(*points);
Expand All @@ -65,15 +65,17 @@ class TauSampling : public AbstractSampling<T> {
return matrix_ * al;
}

Eigen::Matrix<T, Eigen::Dynamic, 1> fit(
const Eigen::Matrix<T, Eigen::Dynamic, 1>& ax,
/*
Eigen::Matrix<double, Eigen::Dynamic, 1> fit(
const Eigen::Matrix<double, Eigen::Dynamic, 1>& ax,
const Eigen::VectorXd* 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_;
Expand All @@ -83,18 +85,19 @@ class TauSampling : public AbstractSampling<T> {
return sampling_points_;
}

const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& matrix() const override {
const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& matrix() const override {
return matrix_;
}

const std::shared_ptr<AbstractBasis<T>>& basis() const override {
/*
const std::shared_ptr<AbstractBasis<double>>& basis() const override {
return basis_;
}
*/

protected:
std::shared_ptr<AbstractSampling<T>> for_sampling_points(
std::shared_ptr<AbstractSampling<S>> for_sampling_points(
const Eigen::VectorXd& x) const override {
return std::make_shared<TauSampling<T>>(basis_, &x);
return std::make_shared<TauSampling>(basis_, &x);
}

private:
Expand All @@ -107,24 +110,24 @@ class TauSampling : public AbstractSampling<T> {
}
}

Eigen::Matrix<T, Eigen::Dynamic, 1> solve(
const Eigen::Matrix<T, Eigen::Dynamic, 1>& ax) const {
Eigen::Matrix<double, Eigen::Dynamic, 1> solve(
const Eigen::Matrix<double, Eigen::Dynamic, 1>& ax) const {
return solver_.solve(ax);
}

double compute_condition_number(
const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& mat) const {
Eigen::JacobiSVD<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>> svd(
const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& mat) const {
Eigen::JacobiSVD<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>> svd(
mat, Eigen::ComputeThinU | Eigen::ComputeThinV);
double cond = svd.singularValues()(0) /
svd.singularValues()(svd.singularValues().size() - 1);
return cond;
}

std::shared_ptr<AbstractBasis<T>> basis_;
std::shared_ptr<AbstractBasis<S>> basis_;
Eigen::VectorXd sampling_points_;
Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> matrix_;
mutable Eigen::ColPivHouseholderQR<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>> solver_;
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> matrix_;
mutable Eigen::ColPivHouseholderQR<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>> solver_;
double cond_;
};

Expand All @@ -144,25 +147,19 @@ class MatsubaraSampling : public AbstractSampling<std::complex<T>> {
construct_matrix();
}

/*
Eigen::Matrix<std::complex<T>, Eigen::Dynamic, 1> evaluate(
const Eigen::Matrix<T, Eigen::Dynamic, 1>& 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<T, Eigen::Dynamic, 1> fit(
const Eigen::Matrix<std::complex<T>, 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_;
Expand All @@ -180,12 +177,6 @@ class MatsubaraSampling : public AbstractSampling<std::complex<T>> {
return basis_;
}

protected:
std::shared_ptr<AbstractSampling<std::complex<T>>> for_sampling_points(
const Eigen::VectorXi& n) const override {
return std::make_shared<MatsubaraSampling<T>>(basis_, positive_only_, &n);
}

private:
void construct_matrix() {
matrix_ = basis_->uhat(sampling_points_);
Expand Down
2 changes: 1 addition & 1 deletion test/sampling.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ TEST_CASE("sampling", "[Sampling]") {
auto basis = std::make_shared<sparseir::FiniteTempBasis<sparseir::Fermionic, sparseir::LogisticKernel>>(
beta, omega_max, std::numeric_limits<double>::quiet_NaN(), sparseir::LogisticKernel(beta * omega_max), sve_result);

//sparseir::TauSampling<double> tau_sampling(basis);
// sparseir::TauSampling<sparseir::Fermionic> TauSampling(basis);
// Here we can check the type or properties of tau_sampling if needed
REQUIRE(true); // Placeholder
}
Expand Down
7 changes: 7 additions & 0 deletions test/sve.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<sparseir::LogisticKernel>(lk);
auto s = sve.s;
//std::cout << "S values: \n" << s << std::endl;
}

TEST_CASE("sve.cpp", "[compute_sve]")
{

Expand Down

0 comments on commit a089c33

Please sign in to comment.