Skip to content

Commit

Permalink
Merge pull request #60 from SpM-lab/terasaki/fix-s-size-zero
Browse files Browse the repository at this point in the history 
Terasaki/fix s size zero
  • Loading branch information
@terasakisatoshi
terasakisatoshi authored Dec 17, 2024
2 parents 9fa8f27 + 375ef5d commit d6092c8
Show file tree
Hide file tree
Showing 5 changed files with 106 additions and 187 deletions.
146 changes: 35 additions & 111 deletions include/sparseir/basis.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
#include <iostream>
#include <stdexcept>
#include <utility>
#include <memory> // for std::shared_ptr
#include <vector>

namespace sparseir {
Expand Down Expand Up @@ -163,7 +164,7 @@ template <typename S, typename K=LogisticKernel>
class FiniteTempBasis : public AbstractBasis<S> {
public:
K kernel;
SVEResult<K> sve_result;
std::shared_ptr<SVEResult<K>> sve_result;
double accuracy;
double beta;
PiecewiseLegendrePolyVector u;
Expand All @@ -172,30 +173,12 @@ class FiniteTempBasis : public AbstractBasis<S> {
PiecewiseLegendreFTVector<S> uhat;
PiecewiseLegendreFTVector<S> uhat_full;

// Constructor
FiniteTempBasis(double beta, double omega_max,
double epsilon = std::numeric_limits<double>::quiet_NaN(),
int max_size = -1)
{
LogisticKernel kernel = LogisticKernel(beta * omega_max);
SVEResult<LogisticKernel> sve_result =
compute_sve<LogisticKernel>(kernel, epsilon);
FiniteTempBasis<S, LogisticKernel>(beta, omega_max, epsilon, kernel,
sve_result, max_size);
}

FiniteTempBasis(double beta, double omega_max, double epsilon,
K kernel)
{
SVEResult<K> sve_result = compute_sve<K>(kernel, epsilon);
int max_size = -1;
FiniteTempBasis<S, K>(beta, omega_max, epsilon, kernel, sve_result,
max_size);
}

FiniteTempBasis(double beta, double omega_max, double epsilon,
K kernel, SVEResult<K> sve_result, int max_size = -1)
K kernel, SVEResult<K> sve_result, int max_size = -1)
{
if (sve_result.s.size() == 0) {
throw std::runtime_error("SVE result sve_result.s is empty");
}
if (beta <= 0.0) {
throw std::domain_error(
"Inverse temperature beta must be positive");
Expand All @@ -204,21 +187,16 @@ class FiniteTempBasis : public AbstractBasis<S> {
throw std::domain_error(
"Frequency cutoff omega_max must be non-negative");
}

this->beta = beta;
this->kernel = kernel;
this->sve_result = sve_result;
this->sve_result = std::make_shared<SVEResult<K>>(sve_result);

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

// std::tuple<
// PiecewiseLegendrePolyVector,
// Eigen::VectorXd,
// PiecewiseLegendrePolyVector
// >
auto part_result = sve_result.part(epsilon, max_size);
PiecewiseLegendrePolyVector u_ = std::get<0>(part_result);
Eigen::VectorXd s_ = std::get<1>(part_result);
PiecewiseLegendrePolyVector v_ = std::get<2>(part_result);

double sve_result_s0 = sve_result.s(0);

if (sve_result.s.size() > s_.size()) {
Expand All @@ -227,25 +205,7 @@ class FiniteTempBasis : public AbstractBasis<S> {
this->accuracy = sve_result.s(s_.size() - 1) / sve_result_s0;
}

double wmax = kernel.lambda_ / beta;
Eigen::VectorXd u_knots = u_[0].knots;
Eigen::VectorXd v_knots = v_[0].knots;
Eigen::VectorXd u_delta_x = u_[0].delta_x;
Eigen::VectorXd v_delta_x = v_[0].delta_x;
Eigen::VectorXi u_symm = Eigen::VectorXi::Zero(u_.size());
Eigen::VectorXi v_symm = Eigen::VectorXi::Zero(v_.size());
for (int i = 0; i < u_.size(); ++i) {
u_symm[i] = u_[i].symm;
v_symm[i] = v_[i].symm;
}

u_knots = (beta / 2) * (u_knots.array() + 1);
v_knots = wmax * v_knots;
this->u = PiecewiseLegendrePolyVector(u_, u_knots, u_delta_x, u_symm);
this->v = PiecewiseLegendrePolyVector(v_, v_knots, v_delta_x, v_symm);

this->s =
std::sqrt(beta / 2 * wmax) * std::pow(wmax, -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 =
Expand All @@ -256,63 +216,27 @@ class FiniteTempBasis : public AbstractBasis<S> {
uhat_base_full, statistics, kernel.conv_radius());

std::vector<PiecewiseLegendreFT<S>> uhat_polyvec;
for (int i = 0; i < s.size(); ++i) {
for (int i = 0; i < this->s.size(); ++i) {
uhat_polyvec.push_back(this->uhat_full[i]);
}
}
this->uhat = PiecewiseLegendreFTVector<S>(uhat_polyvec);
}
/*
FiniteTempBasis(S statistics, double beta, double omega_max,
double epsilon = 0.0, int max_size = -1, K kernel = K(),
SVEResult<K> sve_result = SVEResult<K>())
: kernel(kernel), sve_result(sve_result), beta(beta)
{
if (beta <= 0.0)
throw std::domain_error("Inverse temperature β must be positive");
if (omega_max < 0.0)
throw std::domain_error(
"Frequency cutoff ωmax must be non-negative");
// Partition the SVE result
auto part_result = sve_result.part(epsilon, max_size);
auto u_ = std::get<0>(part_result);
auto s_ = std::get<1>(part_result);
auto v_ = std::get<2>(part_result);

int L = static_cast<int>(s_.size());
// Calculate accuracy
if (sve_result.s.size() > s_.size()) {
accuracy = sve_result.s[s_.size()] / sve_result.s[0];
} else {
accuracy = sve_result.s.back() / sve_result.s[0];
}
// Delegating constructor 1
FiniteTempBasis(double beta, double omega_max,
double epsilon, int max_size)
: FiniteTempBasis(beta, omega_max, epsilon,
LogisticKernel(beta * omega_max),
compute_sve<LogisticKernel>(LogisticKernel(beta * omega_max), epsilon),
max_size){}

// Scaling variables
omega_max = kernel.Lambda() / beta;
Eigen::VectorXd u_knots = (beta / 2.0) * (u_.knots.array() + 1.0);
Eigen::VectorXd v_knots = omega_max * v_.knots;
u = PiecewiseLegendrePolyVector(u_, u_knots, (beta / 2.0) * u_.delta_x,
u_.symmetry);
v = PiecewiseLegendrePolyVector(v_, v_knots, omega_max * v_.delta_x,
v_.symmetry);
// Scale singular values
double scale_factor = std::sqrt(beta / 2.0 * omega_max) *
std::pow(omega_max, -kernel.getYPower());
s.resize(L);
for (int i = 0; i < L; ++i)
s[i] = scale_factor * s_[i];
// Fourier transforms scaling
PiecewiseLegendrePolyVector u_base_full(
std::sqrt(beta) * sve_result.u.data, sve_result.u);
uhat_full = PiecewiseLegendreFTVector<S>(u_base_full, statistics,
kernel.convRadius());
uhat = uhat_full.slice(0, L);
}
*/
// Delegating constructor 2
FiniteTempBasis(double beta, double omega_max,
double epsilon, K kernel)
: FiniteTempBasis(beta, omega_max, epsilon,
kernel,
compute_sve<K>(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 Down Expand Up @@ -367,8 +291,14 @@ class FiniteTempBasis : public AbstractBasis<S> {
FiniteTempBasis<S, K> rescale(double new_beta) const
{
double new_omega_max = kernel.lambda_ / new_beta;
return FiniteTempBasis<S, K>(new_beta, new_omega_max, 0.0, kernel,
sve_result, static_cast<int>(s.size()));
return FiniteTempBasis<S, K>(
new_beta,
new_omega_max,
std::numeric_limits<double>::quiet_NaN(),
kernel,
*sve_result,
static_cast<int>(s.size())
);
}

private:
Expand Down Expand Up @@ -470,13 +400,10 @@ inline std::pair<FiniteTempBasis<Fermionic, LogisticKernel>,
)
{
LogisticKernel kernel(beta * omega_max);
std::cout << "sve_result.s.size() = " << sve_result.s.size() << std::endl;
auto basis_f = FiniteTempBasis<Fermionic, LogisticKernel>(
beta, omega_max, epsilon, kernel, sve_result);
std::cout << "basis_f.s.size() = " << basis_f.s.size() << std::endl;
auto basis_b = FiniteTempBasis<Bosonic, LogisticKernel>(
beta, omega_max, epsilon, kernel, sve_result);
std::cout << "basis_b.s.size() = " << basis_b.s.size() << std::endl;
return std::make_pair(basis_f, basis_b);
}

Expand All @@ -488,14 +415,11 @@ inline std::pair<FiniteTempBasis<Fermionic, LogisticKernel>,
)
{
LogisticKernel kernel(beta * omega_max);
SVEResult<LogisticKernel> sve_result = compute_sve(kernel, epsilon);
std::cout << "sve_result.s.size() = " << sve_result.s.size() << std::endl;
SVEResult<LogisticKernel> sve_result = compute_sve<LogisticKernel>(kernel, epsilon);
auto basis_f = FiniteTempBasis<Fermionic, LogisticKernel>(
beta, omega_max, epsilon, kernel, sve_result);
std::cout << "basis_f.s.size() = " << basis_f.s.size() << std::endl;
auto basis_b = FiniteTempBasis<Bosonic, LogisticKernel>(
beta, omega_max, epsilon, kernel, sve_result);
std::cout << "basis_b.s.size() = " << basis_b.s.size() << std::endl;
return std::make_pair(basis_f, basis_b);
}
} // namespace sparseir
2 changes: 0 additions & 2 deletions include/sparseir/kernel.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -203,8 +203,6 @@ inline double callreduced(const AbstractReducedKernel &kernel, double x,
*/
class LogisticKernel : public AbstractKernel {
public:
double lambda_; ///< The kernel cutoff Λ.

// Default constructor
LogisticKernel() : AbstractKernel() {}

Expand Down
27 changes: 13 additions & 14 deletions include/sparseir/sve.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -216,8 +216,8 @@ class SamplingSVE : public AbstractSVE<K, T> {
const Eigen::MatrixX<T> &u = u_list[0];
const Eigen::VectorX<T> &s_ = s_list[0];
const Eigen::MatrixX<T> &v = v_list[0];
Eigen::VectorXd s = s_.template cast<double>();

Eigen::VectorXd s = s_.template cast<double>();
Eigen::VectorX<T> gauss_x_w =
Eigen::VectorX<T>::Map(gauss_x.w.data(), gauss_x.w.size());
Eigen::VectorX<T> gauss_y_w =
Expand Down Expand Up @@ -257,7 +257,6 @@ class SamplingSVE : public AbstractSVE<K, T> {
}
}
}

// Manually compute differences for dsegs_x and dsegs_y
Eigen::VectorX<T> dsegs_x(segs_x.size() - 1);
for (int i = 0; i < segs_x.size() - 1; ++i) {
Expand Down Expand Up @@ -298,31 +297,30 @@ class SamplingSVE : public AbstractSVE<K, T> {
for (int i = 0; i < u_data.dimension(2); ++i) {
Eigen::MatrixXd slice_double(u_data.dimension(0),
u_data.dimension(1));
for (int j = 0; j < u_data.dimension(1); ++j) {
for (int k = 0; k < u_data.dimension(2); ++k) {
for (int j = 0; j < u_data.dimension(0); ++j) {
for (int k = 0; k < u_data.dimension(1); ++k) {
slice_double(j, k) = static_cast<double>(u_data(j, k, i));
}
}

std::vector<double> segs_x_double;
segs_x_double.reserve(segs_x.size());
for (const auto &x : segs_x) {
segs_x_double.push_back(static_cast<double>(x));
}

polyvec_u.push_back(PiecewiseLegendrePoly(
slice_double,
Eigen::VectorXd::Map(segs_x_double.data(),
segs_x_double.size()),
i));
}


// Repeat similar changes for v_data
for (int i = 0; i < v_data.dimension(2); ++i) {
Eigen::MatrixXd slice_double(v_data.dimension(0),
v_data.dimension(1));
for (int j = 0; j < v_data.dimension(1); ++j) {
for (int k = 0; k < v_data.dimension(2); ++k) {
for (int j = 0; j < v_data.dimension(0); ++j) {
for (int k = 0; k < v_data.dimension(1); ++k) {
slice_double(j, k) = static_cast<double>(v_data(j, k, i));
}
}
Expand All @@ -340,6 +338,9 @@ class SamplingSVE : public AbstractSVE<K, T> {
i));
}




PiecewiseLegendrePolyVector ulx(polyvec_u);
PiecewiseLegendrePolyVector vly(polyvec_v);
canonicalize(ulx, vly);
Expand Down Expand Up @@ -559,8 +560,6 @@ auto pre_postprocess(K &kernel, double safe_epsilon, int n_gauss,
T cutoff_actual = std::isnan(cutoff)
? 2 * T(std::numeric_limits<double>::epsilon())
: T(cutoff);
std::cout << "Cutoff: " << cutoff_actual << std::endl;

std::vector<Eigen::MatrixX<T>> u_list_truncated;
std::vector<Eigen::VectorX<T>> s_list_truncated;
std::vector<Eigen::MatrixX<T>> v_list_truncated;
Expand All @@ -583,12 +582,12 @@ SVEResult<K> compute_sve(K kernel, double epsilon = std::numeric_limits<double>:
double safe_epsilon;
std::string Twork_actual;
std::string svd_strategy_actual;
std::cout << "Twork: " << Twork << std::endl;
std::cout << "svd_strat: " << svd_strat << std::endl;
//std::cout << "Twork: " << Twork << std::endl;
//std::cout << "svd_strat: " << svd_strat << std::endl;
std::tie(safe_epsilon, Twork_actual, svd_strategy_actual) =
choose_accuracy(epsilon, Twork, svd_strat);
std::cout << "Twork_actual: " << Twork_actual << std::endl;
std::cout << "svd_strategy_actual: " << svd_strategy_actual << std::endl;
//std::cout << "Twork_actual: " << Twork_actual << std::endl;
//std::cout << "svd_strategy_actual: " << svd_strategy_actual << std::endl;

if (Twork_actual == "Float64") {
return pre_postprocess<K, double>(kernel, safe_epsilon, n_gauss, cutoff,
Expand Down
20 changes: 8 additions & 12 deletions test/basis.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,19 +77,15 @@ TEST_CASE("FiniteTempBasis consistency tests", "[basis]") {
sparseir::LogisticKernel(beta * omega_max));
sparseir::SVEResult<sparseir::LogisticKernel> sve =
sparseir::compute_sve(sparseir::LogisticKernel(beta * omega_max));

REQUIRE(sve.s.size() > 0);
REQUIRE(basis.s.size() > 0);
double scale = std::sqrt(beta / 2.0 * omega_max);
/*
std::cout << "scale = " << scale << std::endl;
// Ensure the correct function or member is used for singular values
Eigen::VectorXd scaled_s_eigen = sve.s * scale;
std::cout << "scaled_s_eigen " << scaled_s_eigen << std::endl;
std::cout << "basis.s.size()" << basis.s.size() << std::endl;
REQUIRE(basis.s.size() == sve.s.size());
REQUIRE(basis.s.isApprox(scaled_s_eigen));
// Access accuracy as a member variable if it's not a function
REQUIRE(basis.accuracy == sve.s(sve.s.size() - 1) / sve.s(0));
*/
REQUIRE(std::abs(basis.accuracy - sve.s(sve.s.size() - 1) / sve.s(0)) < 1e-10);
}

SECTION("Rescaling test") {
Expand All @@ -98,9 +94,9 @@ TEST_CASE("FiniteTempBasis consistency tests", "[basis]") {
double epsilon = 1e-6;

// Specify both template parameters
//sparseir::FiniteTempBasis<sparseir::Fermionic, sparseir::LogisticKernel> basis(beta, omega_max, epsilon, sparseir::LogisticKernel(beta * omega_max));
//sparseir::FiniteTempBasis<sparseir::Fermionic, sparseir::LogisticKernel> rescaled_basis = basis.rescale(2.0);

//REQUIRE(rescaled_basis.get_wmax() == 6.0);
sparseir::FiniteTempBasis<sparseir::Fermionic, sparseir::LogisticKernel> basis(beta, omega_max, epsilon, sparseir::LogisticKernel(beta * omega_max));
sparseir::FiniteTempBasis<sparseir::Fermionic, sparseir::LogisticKernel> rescaled_basis = basis.rescale(2.0);
REQUIRE(rescaled_basis.sve_result->s.size() == basis.sve_result->s.size());
REQUIRE(rescaled_basis.get_wmax() == 6.0);
}
}
Loading

0 comments on commit d6092c8

Please sign in to comment.