Skip to content

Commit

Permalink
Merge pull request #91 from SpM-lab/terasaki/update-sampling-
Browse files Browse the repository at this point in the history 
Update sampling
  • Loading branch information
@terasakisatoshi
terasakisatoshi authored Jan 30, 2025
2 parents 4f7aebd + 132b795 commit 645823c
Show file tree
Hide file tree
Showing 5 changed files with 215 additions and 250 deletions.
33 changes: 33 additions & 0 deletions include/sparseir/basis.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,39 @@ class FiniteTempBasis : public AbstractBasis<S> {
auto uk = u_.polyvec[0].knots;
std::cout << "uk.size(): " << uk.size() << std::endl;
std::cout << "uk: \n" << uk << std::endl;
/*
Port the following Julia code to C++:
# The polynomials are scaled to the new variables by transforming
the # knots according to: tau = β/2 * (x + 1), w = ωmax * y. Scaling # the
data is not necessary as the normalization is inferred. ωmax = Λ(kernel) / β
u_knots = (β / 2) .* (knots(u_) .+ 1)
v_knots = ωmax .* knots(v_)
u = PiecewiseLegendrePolyVector(u_, u_knots; Δx=(β / 2) .* Δx(u_),
symm=symm(u_)) v = PiecewiseLegendrePolyVector(v_, v_knots; Δx=ωmax .*
Δx(v_), symm=symm(v_))
*/

Eigen::VectorXd u_knots = (beta / 2) * (uk.array() + 1);
Eigen::VectorXd v_knots = wmax * uk;

Eigen::VectorXd deltax4u = (beta / 2) * u_.polyvec[0].get_delta_x();
Eigen::VectorXd deltax4v = wmax * v_.polyvec[0].get_delta_x();
std::vector<int> u_symm_vec;
for (int i = 0; i < u_.size(); ++i) {
u_symm_vec.push_back(u_.polyvec[i].get_symm());
}
std::vector<int> v_symm_vec;
for (int i = 0; i < v_.size(); ++i) {
v_symm_vec.push_back(v_.polyvec[i].get_symm());
}

Eigen::VectorXi u_symm = Eigen::Map<Eigen::VectorXi>(u_symm_vec.data(), u_symm_vec.size());
Eigen::VectorXi v_symm = Eigen::Map<Eigen::VectorXi>(v_symm_vec.data(), v_symm_vec.size());

PiecewiseLegendrePolyVector u_new = PiecewiseLegendrePolyVector(u_, u_knots, deltax4u, u_symm);
PiecewiseLegendrePolyVector v_new = PiecewiseLegendrePolyVector(v_, v_knots, deltax4v, v_symm);
this->u = u_new;
this->v = v_new;
this->s = (std::sqrt(beta / 2 * wmax) * std::pow(wmax, -(kernel.ypower()))) * s_;

Eigen::Tensor<double, 3> udata3d = sve_result.u.get_data();
Expand Down Expand Up @@ -187,6 +219,7 @@ class FiniteTempBasis : public AbstractBasis<S> {
const Eigen::VectorXd default_tau_sampling_points() const override {
int sz = this->sve_result->s.size();
auto x = default_sampling_points(this->sve_result->u, sz);
std::cout << "x=" << x << std::endl;
return (this->beta / 2.0) * (x.array() + 1.0);
}

Expand Down
242 changes: 85 additions & 157 deletions include/sparseir/sampling.hpp
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
#pragma once

#include <Eigen/Dense>
#include <Eigen/SVD>
#include <memory>
#include <stdexcept>
#include <vector>
Expand All @@ -13,72 +14,115 @@ class AbstractSampling {
virtual ~AbstractSampling() = default;

// Evaluate the basis coefficients at sampling points
virtual Eigen::Matrix<double, Eigen::Dynamic, 1> evaluate(
const Eigen::Matrix<double, Eigen::Dynamic, 1>& al,
virtual Eigen::VectorXd evaluate(
const Eigen::VectorXd& al,
const Eigen::VectorXd* points = nullptr) const = 0;

// Fit values at sampling points to basis coefficients
virtual Eigen::Matrix<double, Eigen::Dynamic, 1> fit(
const Eigen::Matrix<double, Eigen::Dynamic, 1>& ax,
virtual Eigen::VectorXd fit(
const Eigen::VectorXd& ax,
const Eigen::VectorXd* points = nullptr) const = 0;

// Condition number of the transformation matrix
virtual double cond() const = 0;

// Get the sampling points
virtual const Eigen::VectorXd& sampling_points() const = 0;
};
// Helper function declarations
// Forward declarations
template <typename S>
class TauSampling;

// Get the transformation matrix
virtual const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& matrix() const = 0;
template <typename S>
inline Eigen::MatrixXd eval_matrix(const TauSampling<S>* tau_sampling,
const std::shared_ptr<FiniteTempBasis<S>>& basis,
const Eigen::VectorXd& x){
// Initialize matrix with correct dimensions
Eigen::MatrixXd matrix(x.size(), basis->size());

// Get the associated basis
//virtual const std::shared_ptr<AbstractBasis<S>>& basis() const = 0;
// Evaluate basis functions at sampling points
auto u_eval = basis->u(x);
std::cout << "u_eval = " << u_eval << std::endl;
// Transpose and scale by singular values
matrix = u_eval.transpose();

protected:
// Create a new sampling object for different sampling points
virtual std::shared_ptr<AbstractSampling> for_sampling_points(
const Eigen::VectorXd& x) const = 0;
};
return matrix;
}

template <typename S>
class TauSampling : public AbstractSampling<S> {
public:
TauSampling(
const std::shared_ptr<AbstractBasis<S>>& basis,
const Eigen::VectorXd* sampling_points = nullptr)
: basis_(basis) {
if (sampling_points) {
sampling_points_ = *sampling_points;
} else {
sampling_points_ = default_tau_sampling_points(basis_);
const std::shared_ptr<FiniteTempBasis<S>>& basis) : basis_(basis) {
// Get default sampling points from basis
sampling_points_ = basis_->default_tau_sampling_points();

// Ensure matrix dimensions are correct
if (sampling_points_.size() == 0) {
throw std::runtime_error("No sampling points generated");
}

// Initialize evaluation matrix with correct dimensions
matrix_ = eval_matrix(this, basis_, sampling_points_);
std::cout << "matrix_ = " << matrix_ << std::endl;
// Check matrix dimensions
if (matrix_.rows() != sampling_points_.size() ||
matrix_.cols() != basis_->size()) {
throw std::runtime_error(
"Matrix dimensions mismatch: got " +
std::to_string(matrix_.rows()) + "x" +
std::to_string(matrix_.cols()) +
", expected " + std::to_string(sampling_points_.size()) +
"x" + std::to_string(basis_->size()));
}
construct_matrix();

// Initialize SVD
matrix_svd_ = Eigen::JacobiSVD<Eigen::MatrixXd>(
matrix_,
Eigen::ComputeFullU | Eigen::ComputeFullV
);
}

Eigen::Matrix<double, Eigen::Dynamic, 1> evaluate(
const Eigen::Matrix<double, Eigen::Dynamic, 1>& al,
Eigen::VectorXd evaluate(
const Eigen::VectorXd& al,
const Eigen::VectorXd* points = nullptr) const override {
if (points) {
auto sampling = for_sampling_points(*points);
return sampling->matrix() * al;
auto eval_mat = eval_matrix(this, basis_, *points);
if (eval_mat.cols() != al.size()) {
throw std::runtime_error(
"Input vector size mismatch: got " +
std::to_string(al.size()) +
", expected " + std::to_string(eval_mat.cols()));
}
return eval_mat * al;
}

if (matrix_.cols() != al.size()) {
throw std::runtime_error(
"Input vector size mismatch: got " +
std::to_string(al.size()) +
", expected " + std::to_string(matrix_.cols()));
}
return matrix_ * al;
}

/*
Eigen::Matrix<double, Eigen::Dynamic, 1> fit(
const Eigen::Matrix<double, Eigen::Dynamic, 1>& ax,
Eigen::VectorXd fit(
const Eigen::VectorXd& ax,
const Eigen::VectorXd* points = nullptr) const override {
if (points) {
auto sampling = for_sampling_points(*points);
return sampling->solve(ax);
auto eval_mat = eval_matrix(this, basis_, *points);
Eigen::JacobiSVD<Eigen::MatrixXd> local_svd(
eval_mat,
Eigen::ComputeFullU | Eigen::ComputeFullV
);
return local_svd.solve(ax);
}
return solve(ax);
}
*/

double cond() const override {
return cond_;
if (ax.size() != matrix_.rows()) {
throw std::runtime_error(
"Input vector size mismatch: got " +
std::to_string(ax.size()) +
", expected " + std::to_string(matrix_.rows()));
}
return matrix_svd_.solve(ax);
}

const Eigen::VectorXd& sampling_points() const override {
Expand All @@ -89,128 +133,12 @@ class TauSampling : public AbstractSampling<S> {
return sampling_points_;
}

const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& matrix() const override {
return matrix_;
}
/*
const std::shared_ptr<AbstractBasis<double>>& basis() const override {
return basis_;
}
*/

protected:
std::shared_ptr<AbstractSampling<S>> for_sampling_points(
const Eigen::VectorXd& x) const override {
return std::make_shared<TauSampling>(basis_, &x);
}

private:
void construct_matrix() {
matrix_ = basis_->u(sampling_points_).transpose();
cond_ = compute_condition_number(matrix_);
solver_.compute(matrix_);
if (solver_.info() != Eigen::Success) {
throw std::runtime_error("Matrix decomposition failed.");
}
}

Eigen::Matrix<double, Eigen::Dynamic, 1> solve(
const Eigen::Matrix<double, Eigen::Dynamic, 1>& ax) const {
return solver_.solve(ax);
}

double compute_condition_number(
const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& mat) const {
Eigen::JacobiSVD<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>> svd(
mat, Eigen::ComputeThinU | Eigen::ComputeThinV);
double cond = svd.singularValues()(0) /
svd.singularValues()(svd.singularValues().size() - 1);
return cond;
}

std::shared_ptr<AbstractBasis<S>> basis_;
std::shared_ptr<FiniteTempBasis<S>> basis_;
Eigen::VectorXd sampling_points_;
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> matrix_;
mutable Eigen::ColPivHouseholderQR<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>> solver_;
double cond_;
Eigen::MatrixXd matrix_;
Eigen::JacobiSVD<Eigen::MatrixXd> matrix_svd_;
};

template <typename T>
class MatsubaraSampling : public AbstractSampling<std::complex<T>> {
public:
MatsubaraSampling(
const std::shared_ptr<AbstractBasis<T>>& basis,
bool positive_only = false,
const Eigen::VectorXi* sampling_points = nullptr)
: basis_(basis), positive_only_(positive_only) {
if (sampling_points) {
sampling_points_ = *sampling_points;
} else {
sampling_points_ = basis_->default_matsubara_sampling_points(positive_only);
}
construct_matrix();
}

/*
Eigen::Matrix<std::complex<T>, Eigen::Dynamic, 1> evaluate(
const Eigen::Matrix<T, Eigen::Dynamic, 1>& al,
const Eigen::VectorXi* points = nullptr) const override {
if (points) {
// TODO: Implement for_sampling_points
auto sampling = for_sampling_points(*points);
return sampling->matrix() * al;
}
return matrix_ * al;
}
*/


double cond() const override {
return cond_;
}

const Eigen::VectorXi& sampling_points() const override {
return sampling_points_;
}

const Eigen::Matrix<std::complex<T>, Eigen::Dynamic, Eigen::Dynamic>& matrix() const override {
return matrix_;
}

const std::shared_ptr<AbstractBasis<T>>& basis() const override {
return basis_;
}

private:
void construct_matrix() {
matrix_ = basis_->uhat(sampling_points_);
cond_ = compute_condition_number(matrix_);
solver_.compute(matrix_);
if (solver_.info() != Eigen::Success) {
throw std::runtime_error("Matrix decomposition failed.");
}
}

Eigen::Matrix<T, Eigen::Dynamic, 1> solve(
const Eigen::Matrix<std::complex<T>, Eigen::Dynamic, 1>& ax) const {
return solver_.solve(ax.real());
}

double compute_condition_number(
const Eigen::Matrix<std::complex<T>, Eigen::Dynamic, Eigen::Dynamic>& mat) const {
Eigen::JacobiSVD<Eigen::Matrix<std::complex<T>, Eigen::Dynamic, Eigen::Dynamic>> svd(
mat, Eigen::ComputeThinU | Eigen::ComputeThinV);
double cond = svd.singularValues()(0).real() /
svd.singularValues()(svd.singularValues().size() - 1).real();
return cond;
}

std::shared_ptr<AbstractBasis<T>> basis_;
bool positive_only_;
Eigen::VectorXi sampling_points_;
Eigen::Matrix<std::complex<T>, Eigen::Dynamic, Eigen::Dynamic> matrix_;
mutable Eigen::ColPivHouseholderQR<Eigen::Matrix<std::complex<T>, Eigen::Dynamic, Eigen::Dynamic>> solver_;
double cond_;
};

} // namespace sparseir
2 changes: 1 addition & 1 deletion test/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ add_executable(libsparseirtests
sve.cxx
basis.cxx
augment.cxx
#sampling.cxx
sampling.cxx
#basis_set.cxx
)

Expand Down
14 changes: 9 additions & 5 deletions test/augment.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -169,11 +169,15 @@ TEST_CASE("Augmented bosonic basis") {


TEST_CASE("Vertex basis with stat = $stat", "[augment]") {

for (const shared_ptr<Statistics>& stat : {static_pointer_cast<Statistics>(make_shared<Fermionic>()), static_pointer_cast<Statistics>(make_shared<Bosonic>())}) {
using T = double;
T beta = 1000.0;
T wmax = 2.0;
// Fix: Use proper template parameters for Statistics
std::vector<std::shared_ptr<Statistics>> stats = {
std::make_shared<Fermionic>(),
std::make_shared<Bosonic>()
};

for (const auto& stat : stats) {
double beta = 1000.0;
double wmax = 2.0;
/*
auto basis = make_shared<FiniteTempBasis<T>>(stat, beta, wmax, 1e-6);
vector<shared_ptr<AbstractAugmentation<T>>> augmentations;
Expand Down
Loading

0 comments on commit 645823c

Please sign in to comment.