Line data    Source code

       1             : // SPDX-FileCopyrightText: 2025 PairInteraction Developers
       2             : // SPDX-License-Identifier: LGPL-3.0-or-later
       3             : 
       4             : #include "pairinteraction/system/GreenTensor.hpp"
       5             : 
       6             : #include "pairinteraction/utils/eigen_assertion.hpp"
       7             : #include "pairinteraction/utils/eigen_compat.hpp"
       8             : #include "pairinteraction/utils/spherical.hpp"
       9             : #include "pairinteraction/utils/traits.hpp"
      10             : 
      11             : #include <Eigen/Dense>
      12             : #include <Eigen/Sparse>
      13             : #include <complex>
      14             : #include <map>
      15             : #include <spdlog/spdlog.h>
      16             : #include <unsupported/Eigen/Splines>
      17             : 
      18             : namespace pairinteraction {
      19             : template <typename Scalar>
      20        1632 : GreenTensor<Scalar>::ConstantEntry::ConstantEntry(int row, int col, Scalar val)
      21        1632 :     : row_(row), col_(col), val_(val) {}
      22             : 
      23             : template <typename Scalar>
      24        1668 : Scalar GreenTensor<Scalar>::ConstantEntry::val() const {
      25        1668 :     return val_;
      26             : }
      27             : 
      28             : template <typename Scalar>
      29        3297 : int GreenTensor<Scalar>::ConstantEntry::row() const noexcept {
      30        3297 :     return row_;
      31             : }
      32             : 
      33             : template <typename Scalar>
      34        3297 : int GreenTensor<Scalar>::ConstantEntry::col() const noexcept {
      35        3297 :     return col_;
      36             : }
      37             : 
      38             : template <typename Scalar>
      39          18 : GreenTensor<Scalar>::OmegaDependentEntry::OmegaDependentEntry(int row, int col,
      40             :                                                               Eigen::Spline<real_t, 1> real_spline,
      41             :                                                               Eigen::Spline<real_t, 1> imag_spline)
      42          18 :     : row_(row), col_(col), real_spline(std::move(real_spline)),
      43          18 :       imag_spline(std::move(imag_spline)) {}
      44             : 
      45             : template <typename Scalar>
      46          54 : Scalar GreenTensor<Scalar>::OmegaDependentEntry::val(double omega) const {
      47             :     if constexpr (traits::NumTraits<Scalar>::is_complex_v) {
      48          27 :         return {real_spline(omega)(0), imag_spline(omega)(0)};
      49             :     } else {
      50          27 :         return real_spline(omega)(0);
      51             :     }
      52             : }
      53             : 
      54             : template <typename Scalar>
      55          54 : int GreenTensor<Scalar>::OmegaDependentEntry::row() const noexcept {
      56          54 :     return row_;
      57             : }
      58             : 
      59             : template <typename Scalar>
      60          54 : int GreenTensor<Scalar>::OmegaDependentEntry::col() const noexcept {
      61          54 :     return col_;
      62             : }
      63             : 
      64             : template <typename Scalar>
      65         538 : void GreenTensor<Scalar>::create_entries_from_cartesian(
      66             :     int kappa1, int kappa2, const Eigen::MatrixX<Scalar> &tensor_in_cartesian_coordinates) {
      67             : 
      68         538 :     const real_t scale = tensor_in_cartesian_coordinates.norm();
      69         538 :     const real_t numerical_precision = 100 * scale * std::numeric_limits<real_t>::epsilon();
      70             : 
      71        1076 :     Eigen::SparseMatrix<complex_t> tensor =
      72         538 :         (spherical::get_transformator<complex_t>(kappa1) * tensor_in_cartesian_coordinates *
      73         538 :          spherical::get_transformator<complex_t>(kappa2).adjoint())
      74         538 :             .sparseView(1, numerical_precision);
      75             : 
      76         538 :     std::vector<Entry> entries;
      77        2212 :     for (int k = 0; k < tensor.outerSize(); ++k) {
      78        3306 :         for (typename Eigen::SparseMatrix<complex_t>::InnerIterator it(tensor, k); it; ++it) {
      79             :             if constexpr (traits::NumTraits<Scalar>::is_complex_v) {
      80         807 :                 entries.emplace_back(ConstantEntry(it.row(), it.col(), it.value()));
      81             :             } else {
      82         825 :                 entries.emplace_back(ConstantEntry(it.row(), it.col(), it.value().real()));
      83         825 :                 assert(abs(it.value().imag()) < numerical_precision);
      84             :             }
      85             :         }
      86             :     }
      87         538 :     entries_map[{kappa1, kappa2}] = std::move(entries);
      88         538 : }
      89             : 
      90             : template <typename Scalar>
      91           6 : void GreenTensor<Scalar>::create_entries_from_cartesian(
      92             :     int kappa1, int kappa2,
      93             :     const std::vector<Eigen::MatrixX<Scalar>> &tensors_in_cartesian_coordinates,
      94             :     const std::vector<double> &omegas) {
      95             : 
      96           6 :     if (tensors_in_cartesian_coordinates.size() != omegas.size()) {
      97           0 :         throw std::invalid_argument("The number of tensors and omegas must match.");
      98             :     }
      99             : 
     100           6 :     if (tensors_in_cartesian_coordinates.size() < 4) {
     101           0 :         throw std::invalid_argument(
     102             :             "At least 4 tensors are required for the applied cubic spline interpolation.");
     103             :     }
     104             : 
     105           6 :     auto num_knots = static_cast<int>(omegas.size());
     106           6 :     Eigen::Map<const Eigen::RowVectorXd> knots(omegas.data(), num_knots);
     107             : 
     108           6 :     constexpr int spline_degree = 3; // cubic spline interpolation
     109             : 
     110             :     // Temporary storage wih key = (row, col) and value = vector of one double per omega
     111           6 :     std::map<std::pair<int, int>, std::pair<Eigen::RowVectorXd, Eigen::RowVectorXd>> temp_map;
     112          30 :     for (int idx = 0; idx < num_knots; ++idx) {
     113             : 
     114          24 :         const real_t scale = tensors_in_cartesian_coordinates[idx].norm();
     115          24 :         const real_t numerical_precision = 100 * scale * std::numeric_limits<real_t>::epsilon();
     116             : 
     117          48 :         Eigen::SparseMatrix<complex_t> tensor =
     118          24 :             (spherical::get_transformator<complex_t>(kappa1) *
     119          24 :              tensors_in_cartesian_coordinates[idx] *
     120          24 :              spherical::get_transformator<complex_t>(kappa2).adjoint())
     121          24 :                 .sparseView(1, numerical_precision);
     122             : 
     123          96 :         for (int k = 0; k < tensor.outerSize(); ++k) {
     124         144 :             for (typename Eigen::SparseMatrix<complex_t>::InnerIterator it(tensor, k); it; ++it) {
     125          72 :                 std::pair<int, int> key{it.row(), it.col()};
     126          72 :                 auto &[vec_real, vec_imag] =
     127             :                     temp_map
     128         144 :                         .try_emplace(key, Eigen::RowVectorXd::Zero(num_knots),
     129             :                                      Eigen::RowVectorXd::Zero(num_knots))
     130          72 :                         .first->second;
     131          72 :                 vec_real(idx) = it.value().real();
     132             :                 if constexpr (traits::NumTraits<Scalar>::is_complex_v) {
     133          36 :                     vec_imag(idx) = it.value().imag();
     134             :                 } else {
     135          36 :                     assert(abs(it.value().imag()) < numerical_precision);
     136             :                 }
     137             :             }
     138             :         }
     139             :     }
     140             : 
     141             :     // Set the green tensor entries with spline interpolation
     142           6 :     std::vector<Entry> entries;
     143           6 :     entries.reserve(temp_map.size());
     144          24 :     for (const auto &[key, value] : temp_map) {
     145          18 :         const auto &[vec_real, vec_imag] = value;
     146          18 :         const auto &[row, col] = key;
     147             : 
     148          18 :         Eigen::Spline<real_t, 1> real_spline =
     149             :             Eigen::SplineFitting<Eigen::Spline<real_t, 1>>::Interpolate(vec_real, spline_degree,
     150             :                                                                         knots);
     151             : 
     152          18 :         Eigen::Spline<real_t, 1> imag_spline;
     153             :         if constexpr (traits::NumTraits<Scalar>::is_complex_v) {
     154           9 :             imag_spline = Eigen::SplineFitting<Eigen::Spline<real_t, 1>>::Interpolate(
     155             :                 vec_imag, spline_degree, knots);
     156             :         }
     157             : 
     158          18 :         entries.emplace_back(
     159          36 :             OmegaDependentEntry(row, col, std::move(real_spline), std::move(imag_spline)));
     160             :     }
     161           6 :     entries_map[{kappa1, kappa2}] = std::move(entries);
     162           6 : }
     163             : 
     164             : template <typename Scalar>
     165             : const std::vector<typename GreenTensor<Scalar>::Entry> &
     166        5241 : GreenTensor<Scalar>::get_spherical_entries(int kappa1, int kappa2) const {
     167        5241 :     if (auto it = entries_map.find({kappa1, kappa2}); it != entries_map.end()) {
     168        1602 :         return it->second;
     169             :     }
     170        3639 :     static const std::vector<Entry> empty_entries;
     171        3639 :     return empty_entries;
     172             : }
     173             : 
     174             : // Explicit instantiations
     175             : template class GreenTensor<double>;
     176             : template class GreenTensor<std::complex<double>>;
     177             : } // namespace pairinteraction