Skip to content

Commit

Permalink
Merge pull request #75 from SpM-lab/shinaoka/debug_sve
Browse files Browse the repository at this point in the history 
Shinaoka/debug sve
  • Loading branch information
@shinaoka
shinaoka authored Dec 25, 2024
2 parents a089c33 + 1edb456 commit d83f1f5
Show file tree
Hide file tree
Showing 10 changed files with 717 additions and 550 deletions.
71 changes: 37 additions & 34 deletions include/sparseir/basis.hpp
View file
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
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
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

0 comments on commit d83f1f5

Please sign in to comment.