Implement evaluate

include/sparseir/basis.hpp

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

include/sparseir/sampling.hpp

@@ -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(
@@ -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
@@ -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();
 
@@ -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()) {
@@ -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(
@@ -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);
     }
 

include/sparseir/sve.hpp

@@ -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'");
     }

include/sparseir/utils.hpp

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