Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Shinaoka/debug sve #75

Merged
merged 11 commits into from
Dec 25, 2024
Merged

Shinaoka/debug sve #75

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
016c219
WIP
shinaoka Dec 24, 2024
fe2950c
WIP
shinaoka Dec 25, 2024
24fe860
WIP
shinaoka Dec 25, 2024
e61a823
WIP: Kernel depends on T
shinaoka Dec 25, 2024
fdbb3d2
Bug fixed
shinaoka Dec 25, 2024
8b99d20
LogisticKernelOdd works
shinaoka Dec 25, 2024
0715da9
Add more tests to sve
shinaoka Dec 25, 2024
e016108
Pass basis.cxx
shinaoka Dec 25, 2024
b850a9b
Fix test
shinaoka Dec 25, 2024
a8d09a8
Fix test
shinaoka Dec 25, 2024
1edb456
fix CI
shinaoka Dec 25, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
71 changes: 37 additions & 34 deletions include/sparseir/basis.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,8 @@

namespace sparseir {

Eigen::VectorXd default_sampling_points(const PiecewiseLegendrePolyVector &u, int L);

// Abstract class with S = Fermionic or Bosonic
template <typename S>
class AbstractBasis {
Expand Down Expand Up @@ -160,11 +162,11 @@ class AbstractBasis {

namespace sparseir {

template <typename S, typename K=LogisticKernel>
template <typename S, typename K=LogisticKernel<DDouble>>
class FiniteTempBasis : public AbstractBasis<S> {
public:
K kernel;
std::shared_ptr<SVEResult<K>> sve_result;
std::shared_ptr<K> kernel;
std::shared_ptr<SVEResult> sve_result;
double accuracy;
double beta;
PiecewiseLegendrePolyVector u;
Expand All @@ -174,7 +176,7 @@ class FiniteTempBasis : public AbstractBasis<S> {
PiecewiseLegendreFTVector<S> uhat_full;

FiniteTempBasis(double beta, double omega_max, double epsilon,
K kernel, SVEResult<K> sve_result, int max_size = -1)
std::shared_ptr<K> kernel, SVEResult sve_result, int max_size = -1)
{
if (sve_result.s.size() == 0) {
throw std::runtime_error("SVE result sve_result.s is empty");
Expand All @@ -189,9 +191,9 @@ class FiniteTempBasis : public AbstractBasis<S> {
}
this->beta = beta;
this->kernel = kernel;
this->sve_result = std::make_shared<SVEResult<K>>(sve_result);
this->sve_result = std::make_shared<SVEResult>(sve_result);

double wmax = this->kernel.lambda_ / beta;
double wmax = this->kernel->lambda_ / beta;

auto part_result = sve_result.part(epsilon, max_size);
PiecewiseLegendrePolyVector u_ = std::get<0>(part_result);
Expand All @@ -205,15 +207,15 @@ class FiniteTempBasis : public AbstractBasis<S> {
this->accuracy = sve_result.s(s_.size() - 1) / sve_result_s0;
}

this->s = (std::sqrt(beta / 2 * wmax) * std::pow(wmax, -(this->kernel.ypower()))) * s_;
this->s = (std::sqrt(beta / 2 * wmax) * std::pow(wmax, -(this->kernel->ypower()))) * s_;

Eigen::Tensor<double, 3> udata3d = sve_result.u.get_data();
PiecewiseLegendrePolyVector uhat_base_full =
PiecewiseLegendrePolyVector(sqrt(beta) * udata3d, sve_result.u);
S statistics = S();

this->uhat_full = PiecewiseLegendreFTVector<S>(
uhat_base_full, statistics, kernel.conv_radius());
uhat_base_full, statistics, kernel->conv_radius());

std::vector<PiecewiseLegendreFT<S>> uhat_polyvec;
for (int i = 0; i < this->s.size(); ++i) {
Expand All @@ -223,20 +225,19 @@ class FiniteTempBasis : public AbstractBasis<S> {
}

// Delegating constructor 1
FiniteTempBasis(double beta, double omega_max,
double epsilon, int max_size=-1)
: FiniteTempBasis(beta, omega_max, epsilon,
LogisticKernel(beta * omega_max),
compute_sve<LogisticKernel>(LogisticKernel(beta * omega_max), epsilon),
FiniteTempBasis(double beta, double omega_max, double epsilon, int max_size=-1)
: FiniteTempBasis(beta, omega_max, epsilon, std::make_shared<K>(beta * omega_max),
compute_sve<typename K::ScalarT>(std::make_shared<K>(beta * omega_max), epsilon),
max_size){}

// Delegating constructor 2
FiniteTempBasis(double beta, double omega_max,
double epsilon, K kernel)
double epsilon, std::shared_ptr<K> kernel)
: FiniteTempBasis(beta, omega_max, epsilon,
kernel,
compute_sve<K>(kernel, epsilon),
compute_sve<typename K::ScalarT>(kernel, epsilon),
-1){}

// Overload operator[] for indexing (get a subset of the basis)
FiniteTempBasis<S, K> operator[](const std::pair<int, int> &range) const
{
Expand All @@ -252,22 +253,22 @@ class FiniteTempBasis : public AbstractBasis<S> {
double get_accuracy() const { return accuracy; }

// Getter for ωmax
double get_wmax() const { return kernel.lambda_ / beta; }
double get_wmax() const { return kernel->lambda_ / beta; }

// Getter for SVEResult
const SVEResult<K> &getSVEResult() const { return sve_result; }
std::shared_ptr<const SVEResult> getSVEResult() const { return sve_result; }

// Getter for kernel
const K &getKernel() const { return kernel; }
std::shared_ptr<const K> getKernel() const { return kernel; }

// Getter for Λ
double Lambda() const { return kernel.lambda_; }
double Lambda() const { return kernel->lambda_; }

// Default τ sampling points
Eigen::VectorXd defaultTauSamplingPoints() const
{
Eigen::VectorXd x =
default_sampling_points(sve_result.u, static_cast<int>(s.size()));
default_sampling_points(sve_result->u, static_cast<int>(s.size()));
return (beta / 2.0) * (x.array() + 1.0);
}

Expand All @@ -283,14 +284,14 @@ class FiniteTempBasis : public AbstractBasis<S> {
Eigen::VectorXd defaultOmegaSamplingPoints() const
{
Eigen::VectorXd y =
default_sampling_points(sve_result.v, static_cast<int>(s.size()));
default_sampling_points(sve_result->v, static_cast<int>(s.size()));
return get_wmax() * y.array();
}

// Rescale function
FiniteTempBasis<S, K> rescale(double new_beta) const
{
double new_omega_max = kernel.lambda_ / new_beta;
double new_omega_max = kernel->lambda_ / new_beta;
return FiniteTempBasis<S, K>(
new_beta,
new_omega_max,
Expand Down Expand Up @@ -396,34 +397,36 @@ inline Eigen::VectorXd default_tau_sampling_points(std::shared_ptr<FiniteTempBas
return (basis.beta / 2.0) * (x.array() + 1.0);
}

inline std::pair<FiniteTempBasis<Fermionic, LogisticKernel>,
FiniteTempBasis<Bosonic, LogisticKernel>>
template<typename T>
std::pair<FiniteTempBasis<Fermionic, LogisticKernel<T>>,
FiniteTempBasis<Bosonic, LogisticKernel<T>>>
finite_temp_bases(
double beta, double omega_max,
double epsilon,
SVEResult<LogisticKernel> sve_result
SVEResult sve_result
)
{
LogisticKernel kernel(beta * omega_max);
auto basis_f = FiniteTempBasis<Fermionic, LogisticKernel>(
auto kernel = std::make_shared<LogisticKernel<T>>(beta * omega_max);
auto basis_f = FiniteTempBasis<Fermionic, LogisticKernel<T>>(
beta, omega_max, epsilon, kernel, sve_result);
auto basis_b = FiniteTempBasis<Bosonic, LogisticKernel>(
auto basis_b = FiniteTempBasis<Bosonic, LogisticKernel<T>>(
beta, omega_max, epsilon, kernel, sve_result);
return std::make_pair(basis_f, basis_b);
}

inline std::pair<FiniteTempBasis<Fermionic, LogisticKernel>,
FiniteTempBasis<Bosonic, LogisticKernel>>
template<typename T>
std::pair<FiniteTempBasis<Fermionic, LogisticKernel<T>>,
FiniteTempBasis<Bosonic, LogisticKernel<T>>>
finite_temp_bases(
double beta, double omega_max,
double epsilon = std::numeric_limits<double>::quiet_NaN()
)
{
LogisticKernel kernel(beta * omega_max);
SVEResult<LogisticKernel> sve_result = compute_sve<LogisticKernel>(kernel, epsilon);
auto basis_f = FiniteTempBasis<Fermionic, LogisticKernel>(
auto kernel = std::make_shared<LogisticKernel<T>>(beta * omega_max);
SVEResult sve_result = compute_sve<T>(kernel, epsilon);
auto basis_f = FiniteTempBasis<Fermionic, LogisticKernel<T>>(
beta, omega_max, epsilon, kernel, sve_result);
auto basis_b = FiniteTempBasis<Bosonic, LogisticKernel>(
auto basis_b = FiniteTempBasis<Bosonic, LogisticKernel<T>>(
beta, omega_max, epsilon, kernel, sve_result);
return std::make_pair(basis_f, basis_b);
}
Expand Down
8 changes: 4 additions & 4 deletions include/sparseir/basis_set.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ class FiniteTempBasisSet {

// Constructors
FiniteTempBasisSet(Scalar beta, Scalar omega_max, Scalar epsilon = std::numeric_limits<Scalar>::quiet_NaN(),
const SVEResult<Scalar>& sve_result = SVEResult<Scalar>())
const SVEResult& sve_result = SVEResult())
: beta_(beta), omega_max_(omega_max), epsilon_(epsilon) {
initialize(sve_result);
}
Expand All @@ -38,10 +38,10 @@ class FiniteTempBasisSet {
const std::vector<int>& wn_f() const { return smpl_wn_f_->sampling_frequencies(); }
const std::vector<int>& wn_b() const { return smpl_wn_b_->sampling_frequencies(); }

const SVEResult<Scalar>& sve_result() const { return sve_result_; }
const SVEResult& sve_result() const { return sve_result_; }

private:
void initialize(const SVEResult<Scalar>& sve_result_input) {
void initialize(const SVEResult& sve_result_input) {
if (std::isnan(epsilon_)) {
epsilon_ = std::numeric_limits<Scalar>::epsilon();
}
Expand Down Expand Up @@ -75,7 +75,7 @@ class FiniteTempBasisSet {
MatsubaraSamplingPtr smpl_wn_f_;
MatsubaraSamplingPtr smpl_wn_b_;

SVEResult<Scalar> sve_result_;
SVEResult sve_result_;
};

} // namespace sparseir
44 changes: 27 additions & 17 deletions include/sparseir/gauss.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,22 +77,29 @@ class Rule {
T b = 1)
: x(x), w(w), a(a), b(b)
{
this->x_forward = x_forward.size() == 0
? Eigen::VectorX<T>::Zero(x.size())
: x_forward;
this->x_backward = x_backward.size() == 0
? Eigen::VectorX<T>::Zero(x.size())
: x_backward;
if (x_forward.size() == 0) {
std::transform(x.data(), x.data() + x.size(),
this->x_forward.data(),
[a](T xi) { return xi - a; });
}
if (x_backward.size() == 0) {
std::transform(x.data(), x.data() + x.size(),
this->x_backward.data(),
[b](T xi) { return b - xi; });
}
//this->x_forward = x_forward.size() == 0
//? Eigen::VectorX<T>::Zero(x.size())
//: x_forward;
//this->x_backward = x_backward.size() == 0
//? Eigen::VectorX<T>::Zero(x.size())
//: x_backward;
this->x_forward = Eigen::VectorX<T>::Zero(x.size());
this->x_backward = Eigen::VectorX<T>::Zero(x.size());
std::transform(this->x.data(), this->x.data() + this->x.size(),
this->x_forward.data(), [a](T xi) { return xi - a; });
std::transform(this->x.data(), this->x.data() + this->x.size(),
this->x_backward.data(), [b](T xi) { return b - xi; });
//
//if (x_forward.size() == 0) {
//std::transform(x.data(), x.data() + x.size(),
//this->x_forward.data(),
//[a](T xi) { return xi - a; });
//}
//if (x_backward.size() == 0) {
//std::transform(x.data(), x.data() + x.size(),
//this->x_backward.data(),
//[b](T xi) { return b - xi; });
//}
}

template <typename U>
Expand Down Expand Up @@ -144,7 +151,8 @@ class Rule {
}
std::vector<Rule<T>> rules;
for (size_t i = 0; i < edges.size() - 1; ++i) {
rules.push_back(reseat(T(edges[i]), T(edges[i + 1])));
auto rule_ = reseat(T(edges[i]), T(edges[i + 1]));
rules.push_back(rule_);
}
return join(rules);
}
Expand All @@ -166,6 +174,7 @@ class Rule {
T prev_b = a;
Eigen::VectorX<T> x, w, x_forward, x_backward;

int counter = 0;
for (const auto &curr : gauss_list) {
if (curr.a != prev_b) {
throw std::invalid_argument("Gauss rules must be ascending");
Expand All @@ -191,6 +200,7 @@ class Rule {
w.tail(curr.w.size()) = curr.w;
x_forward.tail(curr_x_forward.size()) = curr_x_forward;
x_backward.tail(curr_x_backward.size()) = curr_x_backward;
counter ++;
}

return Rule<T>(x, w, x_forward, x_backward, a, b);
Expand Down
Loading
Loading