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