Skip to content

Commit

Permalink
Merge pull request #95 from SpM-lab/terasaki/implement-evaluate
Browse files Browse the repository at this point in the history 
Implement evaluate
  • Loading branch information
@terasakisatoshi
terasakisatoshi authored Feb 1, 2025
2 parents 4585c97 + c144077 commit f48e72d
Show file tree
Hide file tree
Showing 6 changed files with 216 additions and 75 deletions.
3 changes: 0 additions & 3 deletions include/sparseir/basis.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,8 +96,6 @@ 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
Expand Down Expand Up @@ -219,7 +217,6 @@ 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
202 changes: 170 additions & 32 deletions include/sparseir/sampling.hpp
View file
Original file line number Diff line number Diff line change
@@ -1,22 +1,186 @@
#pragma once

#include <Eigen/Dense>
#include <unsupported/Eigen/CXX11/Tensor>

#include <Eigen/SVD>
#include <memory>
#include <stdexcept>
#include <vector>
#include <functional>

namespace sparseir {

template <int N>
Eigen::array<int, N> getperm(int src, int dst) {
Eigen::array<int, N> perm;
if (src == dst) {
for (int i = 0; i < N; ++i) {
perm[i] = i;
}
return perm;
}

int pos = 0;
for (int i = 0; i < N; ++i) {
if (i == dst) {
perm[i] = src;
} else {
// src の位置をスキップ
if (pos == src)
++pos;
perm[i] = pos;
++pos;
}
}
return perm;
}

// movedim: テンソル arr の次元 src を次元 dst に移動する（他の次元の順序はそのまま）
template<typename T, int N>
Eigen::Tensor<T, N> movedim(const Eigen::Tensor<T, N>& arr, int src, int dst) {
if (src == dst) {
return arr;
}
auto perm = getperm<N>(src, dst);
return arr.shuffle(perm);
}

/*
A possible C++11/Eigen port of the Julia code that applies a matrix operation
along a chosen dimension of an N-dimensional array.
- matop_along_dim(buffer, mat, arr, dim, op) moves the requested dimension
to the first or last as needed, then calls matop.
- matop(buffer, mat, arr, op, dim) reshapes the tensor into a 2D view and
calls the operation function.
For simplicity, we assume:
1) T is either float/double or a complex type.
2) mat is an Eigen::MatrixXd (real) or possibly Eigen::MatrixXcd if needed.
3) The user-provided op has the signature:
void op(Eigen::Ref<Eigen::MatrixX<T>> out,
const Eigen::Ref<const Eigen::MatrixXd>& mat,
const Eigen::Ref<const Eigen::MatrixX<T>>& in);
…or a variant suitable to your problem.
*/

template <typename T, int N>
Eigen::Tensor<T, N>& matop(
Eigen::Tensor<T, N>& buffer,
const Eigen::MatrixXd& mat,
const Eigen::Tensor<T, N>& arr,
const std::function<void(
Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>>& outMap,
const Eigen::MatrixXd& mat,
const Eigen::Map<const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>>& inMap
)>& op,
int dim)
{
if (dim != 0 && dim != N - 1) {
throw std::domain_error("Dimension must be 0 or N-1 for matop.");
}
if (dim == 0) {
Eigen::Index rowDim = arr.dimension(0);
Eigen::Index colDim = arr.size() / rowDim;

using MatrixType = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>;

Eigen::Map<const MatrixType> inMap(arr.data(), rowDim, colDim);
Eigen::Map<MatrixType> outMap(buffer.data(), mat.rows(), colDim);

outMap = mat * inMap;
for (int i = 0; i < outMap.size(); ++i) {
buffer.data()[i] = outMap.data()[i];
}
} else {
Eigen::Index rowDim = arr.size() / arr.dimension(N - 1);
Eigen::Index colDim = arr.dimension(N - 1);
using MatrixType = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>;

Eigen::Map<const MatrixType> inMap(arr.data(), rowDim, colDim);
Eigen::Map<MatrixType> outMap(buffer.data(), rowDim, mat.rows());

outMap = inMap * mat.transpose();

for (int i = 0; i < outMap.size(); ++i) {
buffer.data()[i] = outMap.data()[i];
}
}
return buffer;
}

template <typename T, int N>
Eigen::Tensor<T, N>& matop_along_dim(
Eigen::Tensor<T, N>& buffer,
const Eigen::MatrixXd& mat,
const Eigen::Tensor<T, N>& arr,
int dim,
const std::function<void(
Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>>& out,
const Eigen::MatrixXd& mat,
const Eigen::Map<const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>>& in
)>& op)
{
if (dim < 0 || dim >= N) {
throw std::domain_error("Dimension must be in [0, N).");
}

if (dim == 0) {
return matop<T, N>(buffer, mat, arr, op, 0);
} else if (dim != N - 1) {
auto perm = getperm<N>(dim, 0);
Eigen::Tensor<T, N> arr_perm = arr.shuffle(perm).eval();
Eigen::Tensor<T, N> buffer_perm = buffer.shuffle(perm).eval();

matop<T, N>(buffer_perm, mat, arr_perm, op, 0);

auto inv_perm = getperm<N>(0, dim);
buffer = buffer_perm.shuffle(inv_perm).eval();
} else {
return matop<T, N>(buffer, mat, arr, op, N - 1);
}

return buffer;
}

template <typename S>
class AbstractSampling {
public:
virtual ~AbstractSampling() = default;

// Evaluate the basis coefficients at sampling points
virtual Eigen::VectorXd evaluate(
const Eigen::VectorXd& al,
const Eigen::VectorXd* points = nullptr) const = 0;
template<typename T, int N>
Eigen::Tensor<T, N> evaluate(const Eigen::Tensor<T, N>& al, int dim = 1) const {
if (dim < 0 || dim >= N) {
throw std::runtime_error(
"evaluate: dimension must be in [0..N). Got dim=" + std::to_string(dim));
}

if (get_matrix().cols() != al.dimension(dim)) {
throw std::runtime_error(
"Mismatch: matrix.cols()=" + std::to_string(get_matrix().cols())
+ ", but al.dimension(" + std::to_string(dim) + ")="
+ std::to_string(al.dimension(dim)));
}

Eigen::Tensor<T, N> buffer(al.dimensions());

std::function<void(
Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>>&,
const Eigen::MatrixXd&,
const Eigen::Map<const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>>&
)> op = [](
Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>>& out,
const Eigen::MatrixXd& mat,
const Eigen::Map<const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>>& in
) {
out.noalias() = mat * in;
};

matop_along_dim<T, N>(buffer, get_matrix(), al, dim, op);
return buffer;
}

// Fit values at sampling points to basis coefficients
virtual Eigen::VectorXd fit(
Expand All @@ -25,6 +189,7 @@ class AbstractSampling {

// Get the sampling points
virtual const Eigen::VectorXd& sampling_points() const = 0;
virtual Eigen::MatrixXd get_matrix() const = 0;
};
// Helper function declarations
// Forward declarations
Expand All @@ -40,7 +205,6 @@ inline Eigen::MatrixXd eval_matrix(const TauSampling<S>* tau_sampling,

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

Expand All @@ -63,7 +227,6 @@ class TauSampling : public AbstractSampling<S> {

// 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()) {
Expand All @@ -84,27 +247,9 @@ class TauSampling : public AbstractSampling<S> {
}
}

Eigen::VectorXd evaluate(
const Eigen::VectorXd& al,
const Eigen::VectorXd* points = nullptr) const override {
if (points) {
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::MatrixXd get_matrix() const override {
return matrix_;
}

Eigen::VectorXd fit(
Expand All @@ -118,13 +263,6 @@ class TauSampling : public AbstractSampling<S> {
);
return local_svd.solve(ax);
}

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

Expand Down
1 change: 0 additions & 1 deletion include/sparseir/sve.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,6 @@ class SVEResult {
inline std::tuple<double, std::string, std::string>
choose_accuracy(double epsilon, const std::string &Twork)
{
std::cout << "Twork: " << Twork << std::endl;
if (Twork != "Float64" && Twork != "Float64x2") {
throw std::invalid_argument("Twork must be either 'Float64' or 'Float64x2'");
}
Expand Down
35 changes: 0 additions & 35 deletions include/sparseir/utils.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,39 +43,4 @@ bool issorted(const Eigen::MatrixBase<Derived>& vec) {
return true;
}

template <int N>
Eigen::array<int, N> getperm(int src, int dst) {
Eigen::array<int, N> perm;
if (src == dst) {
for (int i = 0; i < N; ++i) {
perm[i] = i;
}
return perm;
}

int pos = 0;
for (int i = 0; i < N; ++i) {
if (i == dst) {
perm[i] = src;
} else {
// src の位置をスキップ
if (pos == src)
++pos;
perm[i] = pos;
++pos;
}
}
return perm;
}

// movedim: テンソル arr の次元 src を次元 dst に移動する（他の次元の順序はそのまま）
template<typename T, int N>
Eigen::Tensor<T, N> movedim(const Eigen::Tensor<T, N>& arr, int src, int dst) {
if (src == dst) {
return arr;
}
auto perm = getperm<N>(src, dst);
return arr.shuffle(perm);
}

}
Loading

0 comments on commit f48e72d

Please sign in to comment.