Skip to content

Commit

Permalink
Merge pull request #53 from SpM-lab/terasaki/porting-basis
Browse files Browse the repository at this point in the history 
Update implementation in poly and basis
  • Loading branch information
@terasakisatoshi
terasakisatoshi authored Dec 16, 2024
2 parents f6a8854 + 4275fa3 commit f88e6a1
Show file tree
Hide file tree
Showing 6 changed files with 515 additions and 274 deletions.
191 changes: 140 additions & 51 deletions include/sparseir/basis.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,8 @@

namespace sparseir {

template <typename T>
// Abstract class with S = Fermionic or Bosonic
template <typename S>
class AbstractBasis {
public:
virtual ~AbstractBasis() { }
Expand All @@ -28,7 +29,7 @@ class AbstractBasis {
* @param tau Imaginary time variable.
* @return Value of the l-th basis function at time tau.
*/
virtual T u(int l, T tau) const = 0;
// virtual S u(int l, double tau) const = 0;

/**
* @brief Basis functions on the reduced Matsubara frequency axis.
Expand All @@ -46,14 +47,14 @@ class AbstractBasis {
* @param wn Reduced Matsubara frequency (integer multiplier).
* @return Value of the l-th basis function at frequency wn.
*/
virtual T uhat(int l, int wn) const = 0;
// virtual S uhat(int l, int wn) const = 0;

/**
* @brief Quantum statistic ("F" for fermionic, "B" for bosonic).
*
* @return Character representing the quantum statistics.
*/
virtual char statistics() const = 0;
// virtual S statistics() const = 0;

/**
* @brief Access basis functions/singular values for given index/indices.
Expand All @@ -64,21 +65,21 @@ class AbstractBasis {
* @param index Index or range of indices.
* @return Pointer to the truncated basis (implementation-defined).
*/
virtual AbstractBasis<T> *operator[](int index) const = 0;
// virtual AbstractBasis<S> *operator[](int index) const = 0;

/**
* @brief Shape of the basis function set.
*
* @return Pair representing the shape (rows, columns).
*/
virtual std::pair<int, int> shape() const = 0;
// virtual std::pair<int, int> shape() const = 0;

/**
* @brief Number of basis functions / singular values.
*
* @return Size of the basis function set.
*/
virtual int size() const = 0;
// virtual int size() const = 0;

/**
* @brief Significances of the basis functions.
Expand All @@ -89,7 +90,7 @@ class AbstractBasis {
*
* @return Vector of significance values.
*/
virtual const std::vector<T> &significance() const = 0;
virtual const Eigen::VectorXd significance() const = 0;

/**
* @brief Accuracy of the basis.
Expand All @@ -99,41 +100,41 @@ class AbstractBasis {
*
* @return Accuracy value.
*/
virtual T accuracy() const
virtual double accuracy() const
{
const auto &sig = significance();
return !sig.empty() ? sig.back() : static_cast<T>(0);
const Eigen::VectorXd sig = significance();
return sig.size() > 0 ? sig(sig.size() - 1) : static_cast<double>(0);
}

/**
* @brief Basis cutoff parameter, Λ == β * wmax, or NaN if not present.
*
* @return Cutoff parameter Λ.
*/
virtual T lambda() const = 0;
// virtual double lambda() const = 0;

/**
* @brief Inverse temperature.
*
* @return Value of β.
*/
virtual T beta() const = 0;
// virtual double beta() const = 0;

/**
* @brief Real frequency cutoff or NaN if not present.
*
* @return Maximum real frequency wmax.
*/
virtual T wmax() const = 0;
// virtual double wmax() const = 0;

/**
* @brief Default sampling points on the imaginary time axis.
*
* @param npoints Minimum number of sampling points to return.
* @return Vector of sampling points on the τ-axis.
*/
virtual std::vector<T>
default_tau_sampling_points(int npoints = 0) const = 0;
// virtual Eigen::VectorXd
// default_tau_sampling_points(int npoints = 0) const = 0;

/**
* @brief Default sampling points on the imaginary frequency axis.
Expand All @@ -142,36 +143,117 @@ class AbstractBasis {
* @param positive_only If true, only return non-negative frequencies.
* @return Vector of sampling points on the Matsubara axis.
*/
virtual std::vector<int>
default_matsubara_sampling_points(int npoints = 0,
bool positive_only = false) const = 0;
// virtual Eigen::VectorXd
// default_matsubara_sampling_points(int npoints = 0,
// bool positive_only = false) const = 0;

/**
* @brief Returns true if the sampling is expected to be well-conditioned.
*
* @return True if well-conditioned.
*/
virtual bool is_well_conditioned() const { return true; }
// virtual bool is_well_conditioned() const { return true; }
};

} // namespace sparseir

namespace sparseir {

template <typename S, typename K>
template <typename S, typename K=LogisticKernel>
class FiniteTempBasis : public AbstractBasis<S> {
public:
K kernel;
SVEResult<K> sve_result;
double accuracy;
double beta; // β
double beta;
PiecewiseLegendrePolyVector u;
PiecewiseLegendrePolyVector v;
std::vector<double> s;
Eigen::VectorXd 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)
{
if (beta <= 0.0) {
throw std::domain_error(
"Inverse temperature beta must be positive");
}
if (omega_max < 0.0) {
throw std::domain_error(
"Frequency cutoff omega_max must be non-negative");
}

this->beta = beta;
this->kernel = kernel;
this->sve_result = sve_result;

// 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);

this->accuracy = (sve_result.s.size() > s_.size())
? sve_result.s[s_.size()] / sve_result.s[0]
: sve_result.s[s.size() - 1] / sve_result.s[0];

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;
int u_symm = u_[0].symm;
int v_symm = v_[0].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_;

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());

std::vector<PiecewiseLegendreFT<S>> uhat_polyvec;
for (int i = 0; i < 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>())
Expand Down Expand Up @@ -222,40 +304,23 @@ class FiniteTempBasis : public AbstractBasis<S> {
kernel.convRadius());
uhat = uhat_full.slice(0, L);
}

// Show function (for printing)
void show() const
{
std::cout << s.size() << "-element FiniteTempBasis<" << typeid(S).name()
<< "> with "
<< "β = " << beta << ", ωmax = " << getOmegaMax()
<< " and singular values:\n";
for (size_t i = 0; i < s.size() - 1; ++i)
std::cout << " " << s[i] << "\n";
std::cout << " " << s.back() << "\n";
}

*/
// Overload operator[] for indexing (get a subset of the basis)
FiniteTempBasis<S, K> operator[](const std::pair<int, int> &range) const
{
int new_size = range.second - range.first + 1;
return FiniteTempBasis<S, K>(statistics(), beta, getOmegaMax(), 0.0,
return FiniteTempBasis<S, K>(statistics(), beta, get_wmax(), 0.0,
new_size, kernel, sve_result);
}

// Calculate significance
Eigen::VectorXd significance() const
{
Eigen::VectorXd s_vec =
Eigen::Map<const Eigen::VectorXd>(s.data(), s.size());
return s_vec / s_vec[0];
}
const Eigen::VectorXd significance() const override { return s / s[0]; }

// Getter for accuracy
double getAccuracy() const { return accuracy; }
double get_accuracy() const { return accuracy; }

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

// Getter for SVEResult
const SVEResult<K> &getSVEResult() const { return sve_result; }
Expand All @@ -264,7 +329,7 @@ class FiniteTempBasis : public AbstractBasis<S> {
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
Expand All @@ -287,16 +352,15 @@ class FiniteTempBasis : public AbstractBasis<S> {
{
Eigen::VectorXd y =
default_sampling_points(sve_result.v, static_cast<int>(s.size()));
return getOmegaMax() * y.array();
return get_wmax() * y.array();
}

// Rescale function
FiniteTempBasis<S, K> rescale(double new_beta) const
{
double new_omega_max = kernel.Lambda() / new_beta;
return FiniteTempBasis<S, K>(statistics(), new_beta, new_omega_max, 0.0,
static_cast<int>(s.size()), kernel,
sve_result);
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()));
}

private:
Expand Down Expand Up @@ -389,4 +453,29 @@ class FiniteTempBasis : public AbstractBasis<S> {
}
};

/*
inline std::pair<FiniteTempBasis<Fermionic, LogisticKernel>, FiniteTempBasis<Bosonic,
LogisticKernel>> finite_temp_bases( double beta, double omega_max, double
epsilon = std::numeric_limits<double>::quiet_NaN()
)
{
return std::make_pair<FiniteTempBasis<Fermionic, LogisticKernel>>(beta,
omega_max, epsilon), FiniteTempBasis<Bosonic, LogisticKernel>(beta, omega_max,
epsilon);
}
*/

inline std::pair<FiniteTempBasis<Fermionic, LogisticKernel>,
FiniteTempBasis<Bosonic, LogisticKernel>>
finite_temp_bases(
double beta, double omega_max,
double epsilon = std::numeric_limits<double>::quiet_NaN(),
SVEResult<LogisticKernel> sve_result = SVEResult<LogisticKernel>())
{
LogisticKernel kernel(beta * omega_max);
return std::make_pair(FiniteTempBasis<Fermionic, LogisticKernel>(
beta, omega_max, epsilon, kernel, sve_result),
FiniteTempBasis<Bosonic, LogisticKernel>(
beta, omega_max, epsilon, kernel, sve_result));
}
} // namespace sparseir
6 changes: 5 additions & 1 deletion include/sparseir/kernel.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,8 @@ class AbstractKernel {
public:
double lambda_;
// Constructor
AbstractKernel(double lambda) : lambda_(lambda) { }
AbstractKernel(){}
AbstractKernel(double lambda) : lambda_(lambda) {}

/**
* @brief Evaluate kernel at point (x, y).
Expand Down Expand Up @@ -204,6 +205,9 @@ class LogisticKernel : public AbstractKernel {
public:
double lambda_; ///< The kernel cutoff Λ.

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

/**
* @brief Constructor for LogisticKernel.
*
Expand Down
Loading

0 comments on commit f88e6a1

Please sign in to comment.