-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #62 from SpM-lab/terasaki/porting-augment
Porting augment
- Loading branch information
Showing
6 changed files
with
300 additions
and
2 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,166 @@ | ||
#pragma once | ||
|
||
#include <Eigen/Dense> | ||
#include <cmath> | ||
#include <complex> | ||
#include <limits> | ||
#include <memory> | ||
#include <stdexcept> | ||
#include <vector> | ||
|
||
namespace sparseir { | ||
|
||
template <typename T> | ||
class AbstractAugmentation { | ||
public: | ||
virtual T operator()(T tau) const = 0; // Ensure this is virtual | ||
virtual T deriv(T tau, int n = 1) const = 0; // Add this line | ||
virtual std::complex<T> hat(int n) const = 0; // Add this line | ||
virtual ~AbstractAugmentation() = default; | ||
virtual std::unique_ptr<AbstractAugmentation<T>> clone() const = 0; | ||
}; | ||
|
||
template <typename S> | ||
class AugmentedBasis : public AbstractBasis<S> { | ||
public: | ||
AugmentedBasis(std::shared_ptr<AbstractBasis<S>> basis, | ||
const std::vector<std::unique_ptr<AbstractAugmentation<S>>>& augmentations) | ||
: _basis(basis), _augmentations(augmentations), _naug(augmentations.size()) { | ||
//Error Handling: Check for null basis pointer | ||
if (!_basis) { | ||
throw std::invalid_argument("Basis cannot be null"); | ||
} | ||
//Check for valid augmentations | ||
for (const auto& aug : _augmentations) { | ||
if (!aug) { | ||
throw std::invalid_argument("Augmentation cannot be null"); | ||
} | ||
} | ||
} | ||
|
||
size_t size() const { return _basis->size() + _naug; } | ||
|
||
Eigen::VectorXd u(const Eigen::VectorXd& tau) const { | ||
Eigen::VectorXd result(size()); | ||
for (size_t i = 0; i < _naug; ++i) { | ||
result(i) = (*_augmentations[i])(tau(i)); | ||
} | ||
for (size_t i = _naug; i < size(); ++i) { | ||
result(i) = _basis->u(tau(i - _naug))(i - _naug); | ||
} | ||
return result; | ||
} | ||
|
||
Eigen::VectorXcf uhat(const Eigen::VectorXcf& wn) const { | ||
Eigen::VectorXcf result(size()); | ||
for (size_t i = 0; i < _naug; ++i) { | ||
result(i) = (*_augmentations[i]).hat(wn(i)); | ||
} | ||
for (size_t i = _naug; i < size(); ++i) { | ||
result(i) = _basis->uhat(wn(i - _naug))(i - _naug); | ||
} | ||
return result; | ||
} | ||
|
||
Eigen::VectorXd v(const Eigen::VectorXd& w) const { | ||
return _basis->v(w); | ||
} | ||
|
||
Eigen::VectorXd s() const { return _basis->s(); } | ||
|
||
double beta() const { return _basis->beta(); } | ||
|
||
double wmax() const { return _basis->wmax(); } | ||
|
||
std::shared_ptr<Statistics> statistics() const { | ||
return _basis->statistics(); // Assuming _basis also returns a shared_ptr | ||
} | ||
|
||
private: | ||
std::shared_ptr<AbstractBasis<S>> _basis; | ||
std::vector<std::unique_ptr<AbstractAugmentation<S>>> _augmentations; | ||
size_t _naug; | ||
}; | ||
|
||
|
||
template <typename T> | ||
class TauConst : public AbstractAugmentation<T> { | ||
public: | ||
TauConst(T beta) : beta_(beta), norm_(1.0 / std::sqrt(beta)) { | ||
if (beta_ <= 0) { | ||
throw std::invalid_argument("beta must be positive"); | ||
} | ||
} | ||
|
||
T operator()(T tau) const override { | ||
check_tau_range(tau, beta_); | ||
return norm_; | ||
} | ||
T deriv(T tau, int n = 1) const { | ||
if (n == 0) return (*this)(tau); | ||
return 0.0; | ||
} | ||
std::complex<T> hat(int n) const { | ||
return norm_ * std::sqrt(beta_); | ||
} | ||
std::unique_ptr<AbstractAugmentation<T>> clone() const override { | ||
return std::make_unique<TauConst<T>>(*this); | ||
} | ||
|
||
private: | ||
T beta_; | ||
T norm_; | ||
}; | ||
|
||
template <typename T> | ||
class TauLinear : public AbstractAugmentation<T> { | ||
public: | ||
TauLinear(T beta) : beta_(beta), norm_(std::sqrt(3.0 / beta)) { | ||
if (beta_ <= 0) { | ||
throw std::invalid_argument("beta must be positive"); | ||
} | ||
} | ||
|
||
T operator()(T tau) const override { | ||
check_tau_range(tau, beta_); | ||
return norm_ * (2.0 * tau / beta_ - 1.0); | ||
} | ||
T deriv(T tau, int n = 1) const override { | ||
if (n == 1) return norm_ * 2.0 / beta_; | ||
return 0.0; | ||
} | ||
std::complex<T> hat(int n) const override { | ||
if (n == 0) return 0.0; | ||
return norm_ * 2.0 * beta_ / (n * M_PI * std::complex<T>(0, 1)); | ||
} | ||
std::unique_ptr<AbstractAugmentation<T>> clone() const override { | ||
return std::make_unique<TauLinear<T>>(*this); | ||
} | ||
|
||
private: | ||
T beta_; | ||
T norm_; | ||
}; | ||
|
||
template <typename T> | ||
class MatsubaraConst : public AbstractAugmentation<T> { | ||
public: | ||
MatsubaraConst(T beta) : beta_(beta) { | ||
if (beta_ <= 0) { | ||
throw std::invalid_argument("beta must be positive"); | ||
} | ||
} | ||
|
||
T operator()(T tau) const override { return std::numeric_limits<T>::quiet_NaN(); } | ||
T deriv(T tau, int n = 1) const override { return std::numeric_limits<T>::quiet_NaN(); } | ||
std::complex<T> hat(int n) const override { return 1.0; } | ||
std::unique_ptr<AbstractAugmentation<T>> clone() const override { | ||
return std::make_unique<MatsubaraConst<T>>(*this); | ||
} | ||
|
||
private: | ||
T beta_; | ||
}; | ||
|
||
|
||
} // namespace sparseir |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -11,4 +11,4 @@ | |
#include "./kernel.hpp" | ||
#include "./sve.hpp" | ||
#include "./basis.hpp" | ||
|
||
#include "./augment.hpp" |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,129 @@ | ||
#include <Eigen/Dense> | ||
#include <algorithm> | ||
#include <cmath> | ||
#include <complex> | ||
#include <limits> | ||
#include <memory> | ||
#include <stdexcept> | ||
#include <vector> | ||
|
||
#include <catch2/catch_test_macros.hpp> | ||
|
||
#include <sparseir/sparseir-header-only.hpp> | ||
#include <xprec/ddouble-header-only.hpp> | ||
|
||
using namespace sparseir; | ||
using namespace std; | ||
|
||
TEST_CASE("Augmented bosonic basis") { | ||
using T = double; | ||
T beta = 1000.0; | ||
T wmax = 2.0; | ||
|
||
// Create bosonic basis | ||
auto kernel = LogisticKernel(beta * wmax); | ||
auto sve_result = compute_sve(kernel, 1e-6); | ||
auto basis = make_shared<FiniteTempBasis<Bosonic>>(beta, wmax, 1e-6, kernel, sve_result); | ||
// Create augmented basis with TauConst and TauLinear | ||
/* | ||
vector<unique_ptr<AbstractAugmentation<T>>> augmentations; | ||
//augmentations.push_back(make_unique<TauConst<T>>(beta)); | ||
//augmentations.push_back(make_unique<TauLinear<T>>(beta)); | ||
AugmentedBasis<T> basis_aug(basis, augmentations); | ||
REQUIRE(basis_aug.size() == basis->size() + 2); | ||
// Define G(τ) = c - exp(-τ * pole) / (1 - exp(-β * pole)) | ||
T pole = 1.0; | ||
T c = 1e-2; | ||
// Create tau sampling points | ||
auto tau_sampling = make_tau_sampling(basis_aug); | ||
auto tau = tau_sampling.tau; | ||
Eigen::VectorX<T> gtau(tau.size()); | ||
for (size_t i = 0; i < tau.size(); ++i) { | ||
gtau(i) = c - exp(-tau(i) * pole) / (1 - exp(-beta * pole)); | ||
} | ||
T magn = gtau.maxCoeff(); | ||
// This illustrates that "naive" fitting is a problem if the fitting matrix | ||
// is not well-conditioned. | ||
Eigen::MatrixXd tau_matrix = tau_sampling.matrix; | ||
Eigen::VectorX<T> gl_fit_bad = tau_matrix.completeOrthogonalDecomposition().solve(gtau); | ||
Eigen::VectorX<T> gtau_reconst_bad = tau_matrix * gl_fit_bad; | ||
REQUIRE(!gtau_reconst_bad.isApprox(gtau, 1e-13 * magn)); | ||
REQUIRE(gtau_reconst_bad.isApprox(gtau, 5e-16 * tau_matrix.norm() * magn)); | ||
REQUIRE(tau_matrix.norm() > 1e7); | ||
REQUIRE(tau_matrix.rows() == basis_aug.size()); | ||
REQUIRE(tau_matrix.cols() == tau.size()); | ||
// Now do the fit properly | ||
Eigen::VectorX<T> gl_fit = tau_matrix.colPivHouseholderQr().solve(gtau); | ||
Eigen::VectorX<T> gtau_reconst = tau_matrix * gl_fit; | ||
REQUIRE(gtau_reconst.isApprox(gtau, 1e-14 * magn)); | ||
*/ | ||
} | ||
|
||
|
||
TEST_CASE("Vertex basis with stat = $stat", "[augment]") { | ||
/* | ||
for (const shared_ptr<Statistics>& stat : {make_shared<Fermionic>(), make_shared<Bosonic>()}) { | ||
using T = double; | ||
T beta = 1000.0; | ||
T wmax = 2.0; | ||
auto basis = make_shared<FiniteTempBasis<T>>(stat, beta, wmax, 1e-6); | ||
vector<unique_ptr<AbstractAugmentation<T>>> augmentations; | ||
augmentations.push_back(make_unique<MatsubaraConst<T>>(beta)); | ||
AugmentedBasis<T> basis_aug(basis, augmentations); | ||
REQUIRE(!basis_aug.uhat.empty()); | ||
// G(iν) = c + 1 / (iν - pole) | ||
T pole = 1.0; | ||
T c = 1.0; | ||
auto matsu_sampling = make_matsubara_sampling(basis_aug); | ||
Eigen::VectorXcf gi_n(matsu_sampling.wn.size()); | ||
for (size_t i = 0; i < matsu_sampling.wn.size(); ++i) { | ||
complex<T> iwn(0, matsu_sampling.wn(i)); | ||
gi_n(i) = c + 1.0 / (iwn - pole); | ||
} | ||
Eigen::VectorXcf gl = fit(matsu_sampling, gi_n); | ||
Eigen::VectorXcf gi_n_reconst = evaluate(matsu_sampling, gl); | ||
REQUIRE(gi_n_reconst.isApprox(gi_n, gi_n.maxCoeff() * 1e-7)); | ||
} | ||
} | ||
TEST_CASE("unit tests", "[augment]") { | ||
using T = double; | ||
T beta = 1000.0; | ||
T wmax = 2.0; | ||
auto basis = make_shared<FiniteTempBasis<T>>(Bosonic(), beta, wmax, 1e-6); | ||
vector<unique_ptr<AbstractAugmentation<T>>> augmentations; | ||
augmentations.push_back(make_unique<TauConst<T>>(beta)); | ||
augmentations.push_back(make_unique<TauLinear<T>>(beta)); | ||
AugmentedBasis<T> basis_aug(basis, augmentations); | ||
SECTION("getindex") { | ||
REQUIRE(basis_aug.u.size() == basis_aug.size()); | ||
REQUIRE(basis_aug.u[0]->operator()(0.0) == basis_aug[0](0.0)); | ||
REQUIRE(basis_aug.u[1]->operator()(0.0) == basis_aug[1](0.0)); | ||
} | ||
size_t len_basis = basis->size(); | ||
size_t len_aug = len_basis + 2; | ||
REQUIRE(basis_aug.size() == len_aug); | ||
// REQUIRE(basis_aug.accuracy == basis->accuracy); | ||
REQUIRE(basis_aug.Lambda() == beta * wmax); | ||
REQUIRE(basis_aug.wmax() == wmax); | ||
REQUIRE(basis_aug.u.size() == len_aug); | ||
REQUIRE(basis_aug.u(0.8).size() == len_aug); | ||
//REQUIRE(basis_aug.uhat(MatsubaraFreq<T>(4.0)).size() == len_aug); | ||
REQUIRE(basis_aug.u.minCoeff() == 0.0); | ||
REQUIRE(basis_aug.u.maxCoeff() == beta); | ||
//Further tests omitted for brevity, adapt as needed from Julia code. | ||
*/ | ||
} |