Skip to content

Commit

Permalink
Merge pull request #70 from SpM-lab/terasaki/matsubara-freq-constructor
Browse files Browse the repository at this point in the history 
Improve Matsubara freq constructor
  • Loading branch information
@terasakisatoshi
terasakisatoshi authored Dec 24, 2024
2 parents e50791f + 79201c2 commit 8f86cdd
Show file tree
Hide file tree
Showing 3 changed files with 102 additions and 43 deletions.
4 changes: 2 additions & 2 deletions include/sparseir/augment.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -82,8 +82,8 @@ class TauLinear : public AbstractAugmentation {
}

std::complex<double> operator()(MatsubaraFreq<Fermionic> n) const override {
std::invalid_argument("TauConst is not a Fermionic basis.");
return std::numeric_limits<double>::quiet_NaN();
throw std::invalid_argument("TauConst is not a Fermionic basis.");
return std::numeric_limits<std::complex<double>>::quiet_NaN();
}

std::function<double(double)> deriv(int order = 1) const override {
Expand Down
31 changes: 23 additions & 8 deletions include/sparseir/freq.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ inline std::unique_ptr<Statistics> create_statistics(int zeta)
throw std::domain_error("Unknown statistics type");
}

// Matsubara frequency class template
// MatsubaraFreq class template
template <typename S>
class MatsubaraFreq {
static_assert(std::is_base_of<Statistics, S>::value,
Expand All @@ -50,21 +50,36 @@ class MatsubaraFreq {
public:
int n;

inline MatsubaraFreq(int n) : n(n)
{
// コンストラクタ
inline MatsubaraFreq(int n) : n(n) {
static_assert(std::is_same<S, Fermionic>::value || std::is_same<S, Bosonic>::value,
"S must be Fermionic or Bosonic");
S stat;
if (!stat.allowed(n))
if (!stat.allowed(n)) {
throw std::domain_error("Frequency is not allowed for this type");
}
instance_ = std::make_shared<S>(stat); // 適切な型のインスタンスを保持
}

inline double value(double beta) const { return n * M_PI / beta; }
inline std::complex<double> value_im(double beta) const
{
// 値の計算
inline double value(double beta) const {
return n * M_PI / beta;
}

// 複素数の値
inline std::complex<double> value_im(double beta) const {
return std::complex<double>(0, value(beta));
}

inline S statistics() const { return S(); }
// 統計型の取得
inline S statistics() const { return *std::static_pointer_cast<S>(instance_); }

// n の取得
inline int get_n() const { return n; }

private:
// インスタンスを保持する共有ポインタ
std::shared_ptr<Statistics> instance_;
};

// Typedefs for convenience
Expand Down
110 changes: 77 additions & 33 deletions test/augment.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,48 +7,92 @@
#include <stdexcept>
#include <vector>

#include <catch2/catch_test_macros.hpp>
#include <catch2/catch_approx.hpp> // for Approx

#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("AbstractAugmentation") {
SECTION("TauConst") {
double beta = 1000.0;
auto tc = TauConst(beta);
double tau = 0.5;
double y = tc(tau);
auto dtc = tc.deriv();
double x = 2.0;
REQUIRE(tc.beta == beta);
REQUIRE(dtc(x) == 0.0);
TEST_CASE("SparseIR Basis Functions", "[SparseIR]")
{
using Catch::Approx;
using namespace sparseir;

SECTION("TauConst")
{
REQUIRE_THROWS_AS(TauConst(-34), std::domain_error);

TauConst tc(123);
REQUIRE(tc.beta == 123.0);

REQUIRE_THROWS_AS(tc(-3), std::domain_error);
REQUIRE_THROWS_AS(tc(321), std::domain_error);

REQUIRE(tc(100) == 1 / std::sqrt(123));
//REQUIRE(tc(MatsubaraFreq(0)) == std::sqrt(123));
//REQUIRE(tc(MatsubaraFreq(92)) == 0.0);
//REQUIRE_THROWS_AS(tc(MatsubaraFreq(93)), std::runtime_error);

//REQUIRE(sparseir::deriv(tc)(4.2) == 0.0);
//REQUIRE(sparseir::deriv(tc, 0) == tc);
}
SECTION("TauLinear") {
double beta = 1000.0;
double tau_0 = 0.5;
double tau_1 = 1.0;
auto tl = TauLinear(beta);
double tau = 0.75;
double y = tl(tau);
auto dtl = tl.deriv();
double x = 2.0;
REQUIRE(true);
//REQUIRE(dtl(x) == -beta * (tau - tau_0) / pow(tau_1 - tau_0, 2));

SECTION("TauLinear")
{
REQUIRE_THROWS_AS(TauLinear(-34), std::domain_error);

TauLinear tl(123);
REQUIRE(tl.beta == Approx(123.0));

REQUIRE_THROWS_AS(tl(-3), std::domain_error);
REQUIRE_THROWS_AS(tl(321), std::domain_error);
REQUIRE(tl.norm == Approx(std::sqrt(3.0 / 123.0)));
double tau = 100;
REQUIRE(tl(tau) == std::sqrt(3.0 / 123.0) * (2.0 / 123. * tau - 1.));

MatsubaraFreq<Bosonic> freq0(0);
REQUIRE(tl(freq0) == 0.0);
MatsubaraFreq<Bosonic> freq92(92);
// Calculate the expected complex value
std::complex<double> expected_value = std::sqrt(3. / 123.) * 2. / std::complex<double>(0, 1) * 123. / (92. * M_PI);
// Get the actual value from the function
std::complex<double> actual_value = tl(freq92);

REQUIRE(actual_value.real() == Approx(expected_value.real()));
REQUIRE(actual_value.imag() == Approx(expected_value.imag()));

MatsubaraFreq<Fermionic> freq93(93);
REQUIRE_THROWS_AS(tl(freq93), std::invalid_argument);

//REQUIRE(sparseir::deriv(tl, 0) == tl);
//REQUIRE(sparseir::deriv(tl)(4.2) ==
// Approx(std::sqrt(3 / 123) * 2 / 123));
//REQUIRE(sparseir::deriv(tl, 2)(4.2) == Approx(0.0));
}
SECTION("MatsubaraConst") {
double beta = 1000.0;
double w_0 = 0.5;
double w_1 = 1.0;
auto mc = MatsubaraConst(beta);
double w = 0.75;
double y = mc(w);
auto dmc = mc.deriv();
double x = 2.0;
REQUIRE(true);
//REQUIRE(dmc(x) == -beta * (w - w_0) / pow(w_1 - w_0, 2));

SECTION("MatsubaraConst")
{
REQUIRE_THROWS_AS(MatsubaraConst(-34), std::domain_error);

MatsubaraConst mc(123);
REQUIRE(mc.beta == Approx(123.0));

REQUIRE_THROWS_AS(mc(-3), std::domain_error);
REQUIRE_THROWS_AS(mc(321), std::domain_error);

REQUIRE(std::isnan(mc(100)));
MatsubaraFreq<Bosonic> freq0(0);
REQUIRE(mc(freq0) == 1.0);
MatsubaraFreq<Bosonic> freq92(0);
REQUIRE(mc(freq92) == 1.0);
MatsubaraFreq<Fermionic> freq93(93);
REQUIRE(mc(freq93) == 1.0);

//REQUIRE(sparseir::deriv(mc) == mc);
//REQUIRE(sparseir::deriv(mc, 0) == mc);
}
}

Expand Down

0 comments on commit 8f86cdd

Please sign in to comment.