Line data Source code
1 : // SPDX-FileCopyrightText: 2024 PairInteraction Developers
2 : // SPDX-License-Identifier: LGPL-3.0-or-later
3 :
4 : #include "pairinteraction/basis/Basis.hpp"
5 :
6 : #include "pairinteraction/basis/BasisAtom.hpp"
7 : #include "pairinteraction/basis/BasisPair.hpp"
8 : #include "pairinteraction/enums/SorterType.hpp"
9 : #include "pairinteraction/ket/KetAtom.hpp"
10 : #include "pairinteraction/ket/KetPair.hpp"
11 : #include "pairinteraction/utils/TaskControl.hpp"
12 : #include "pairinteraction/utils/eigen_assertion.hpp"
13 : #include "pairinteraction/utils/eigen_compat.hpp"
14 :
15 : #include <algorithm>
16 : #include <cassert>
17 : #include <cmath>
18 : #include <cstdint>
19 : #include <map>
20 : #include <set>
21 : #include <string>
22 : #include <vector>
23 :
24 : namespace pairinteraction {
25 :
26 : template <typename Scalar>
27 : class BasisAtom;
28 :
29 : template <typename Derived>
30 666 : void Basis<Derived>::perform_blocks_checks(const std::set<SorterType> &unique_labels) const {
31 : // Currently, since states in a basis dont store the energy, they cannot be sorted by energy.
32 666 : if (unique_labels.contains(SorterType::ENERGY)) {
33 1 : throw std::invalid_argument("Blocks cannot be obtained by the energy. Note that sorting "
34 : "a system by the energy is supported nevertheless.");
35 : }
36 :
37 : // Check if the states are labeled by the requested labels
38 634 : std::vector<const std::vector<real_t> *> quantum_numbers_of_labels;
39 634 : quantum_numbers_of_labels.reserve(unique_labels.size());
40 1294 : for (const auto &label : unique_labels) {
41 621 : quantum_numbers_of_labels.push_back(&get_quantum_numbers_of_states(label));
42 : }
43 :
44 : // Check if the states are sorted by the requested labels, i.e., if states that share the same
45 : // labels are contiguous. Otherwise, states belonging to the same block would be scattered over
46 : // several blocks and couplings between them would be lost.
47 : // The quantum numbers are (half-)integers, thus doubling them yields an exact integer
48 : using label_t = std::vector<int64_t>;
49 37642 : auto label_of = [&](Eigen::Index i) {
50 36993 : label_t label;
51 36993 : label.reserve(quantum_numbers_of_labels.size());
52 75236 : for (const auto *quantum_numbers : quantum_numbers_of_labels) {
53 38250 : label.push_back(std::llround(2 * (*quantum_numbers)[i]));
54 : }
55 36880 : return label;
56 0 : };
57 :
58 : // Collect the labels of the contiguous blocks of states
59 649 : std::vector<label_t> labels_of_blocks;
60 37541 : for (Eigen::Index i = 0; i < coefficients.cols(); ++i) {
61 37040 : label_t label = label_of(i);
62 36884 : if (labels_of_blocks.empty() || labels_of_blocks.back() != label) {
63 1378 : labels_of_blocks.push_back(label);
64 : }
65 : }
66 :
67 : // If the states are sorted, no two blocks share the same labels
68 666 : std::sort(labels_of_blocks.begin(), labels_of_blocks.end());
69 647 : if (std::adjacent_find(labels_of_blocks.begin(), labels_of_blocks.end()) !=
70 1286 : labels_of_blocks.end()) {
71 2 : throw std::invalid_argument("The states are not sorted by the requested labels.");
72 : }
73 645 : }
74 :
75 : template <typename Derived>
76 2218 : Basis<Derived>::Basis(ketvec_t &&kets)
77 4436 : : kets(std::move(kets)), coefficients(static_cast<Eigen::Index>(this->kets.size()),
78 6654 : static_cast<Eigen::Index>(this->kets.size())) {
79 2218 : if (this->kets.empty()) {
80 0 : throw std::invalid_argument("The basis must contain at least one element.");
81 : }
82 2218 : make_canonical();
83 2218 : }
84 :
85 : template <typename Derived>
86 2386 : void Basis<Derived>::make_canonical() {
87 2386 : auto n = static_cast<Eigen::Index>(kets.size());
88 2386 : if (coefficients.rows() != n || coefficients.cols() != n) {
89 0 : throw std::invalid_argument("The coefficient matrix must be of shape n x n, where n is "
90 : "the number of kets.");
91 : }
92 :
93 : // Make the i-th state equal to the i-th ket
94 2386 : coefficients.setIdentity();
95 :
96 : // Assign each state the quantum numbers of the ket it is equal to
97 9544 : for (const auto &[label, name] : Derived::sorter_type_to_quantum_number_name) {
98 7158 : std::vector<real_t> quantum_numbers;
99 7158 : quantum_numbers.reserve(kets.size());
100 8020923 : for (const auto &ket : kets) {
101 8013765 : quantum_numbers.push_back(ket->has_quantum_number(name)
102 13232767 : ? static_cast<real_t>(ket->get_quantum_number(name))
103 5219002 : : std::numeric_limits<real_t>::max());
104 : }
105 7158 : quantum_numbers_of_states[name] = std::move(quantum_numbers);
106 : }
107 :
108 : // Assign each state the energy of the ket it is equal to
109 2386 : energy_of_states.clear();
110 2386 : energy_of_states.reserve(kets.size());
111 2673641 : for (const auto &ket : kets) {
112 2671255 : energy_of_states.push_back(static_cast<real_t>(ket->get_energy()));
113 : }
114 2386 : }
115 :
116 : template <typename Derived>
117 6235 : bool Basis<Derived>::has_quantum_number(const std::string &name) const {
118 6235 : auto it = quantum_numbers_of_states.find(name);
119 6197 : if (it == quantum_numbers_of_states.end()) {
120 0 : return false;
121 : }
122 457486 : return std::none_of(it->second.begin(), it->second.end(), [](real_t quantum_number) {
123 451308 : return quantum_number == std::numeric_limits<real_t>::max();
124 6233 : });
125 : }
126 :
127 : template <typename Derived>
128 1602 : bool Basis<Derived>::has_energy() const {
129 1602 : return std::none_of(energy_of_states.begin(), energy_of_states.end(),
130 192025 : [](real_t energy) { return energy == std::numeric_limits<real_t>::max(); });
131 : }
132 :
133 : template <typename Derived>
134 505 : bool Basis<Derived>::is_sorted_by_energy() const {
135 505 : return has_energy() && std::is_sorted(energy_of_states.begin(), energy_of_states.end());
136 : }
137 :
138 : template <typename Derived>
139 6084 : const Derived &Basis<Derived>::derived() const {
140 6084 : return static_cast<const Derived &>(*this);
141 : }
142 :
143 : template <typename Derived>
144 : const std::vector<typename Basis<Derived>::real_t> &
145 1885 : Basis<Derived>::get_quantum_numbers_of_states(SorterType label) const {
146 1885 : auto it = Derived::sorter_type_to_quantum_number_name.find(label);
147 1846 : if (it == Derived::sorter_type_to_quantum_number_name.end()) {
148 0 : throw std::invalid_argument(
149 : "The states of this basis cannot be labeled by the requested sorter type.");
150 : }
151 1852 : const std::string &name = it->second;
152 1850 : if (!has_quantum_number(name)) {
153 2 : throw std::invalid_argument(
154 : "States cannot be labeled and thus not sorted by the quantum number " + name + ".");
155 : }
156 3800 : return quantum_numbers_of_states.at(name);
157 : }
158 :
159 : template <typename Derived>
160 2357 : const typename Basis<Derived>::ketvec_t &Basis<Derived>::get_kets() const {
161 2357 : return kets;
162 : }
163 :
164 : template <typename Derived>
165 : const Eigen::SparseMatrix<typename Basis<Derived>::scalar_t, Eigen::RowMajor> &
166 15495 : Basis<Derived>::get_coefficients() const {
167 15495 : return coefficients;
168 : }
169 :
170 : template <typename Derived>
171 1129 : std::shared_ptr<const Derived> Basis<Derived>::copy_with_coefficients(
172 : const Eigen::SparseMatrix<scalar_t, Eigen::RowMajor> &values) const {
173 1129 : if (values.rows() != coefficients.rows()) {
174 0 : throw std::invalid_argument("Incompatible number of rows.");
175 : }
176 1129 : if (values.cols() != coefficients.cols()) {
177 0 : throw std::invalid_argument("Incompatible number of columns.");
178 : }
179 :
180 : // Create a copy of the current object and update the coefficients of the copy
181 1129 : auto result = std::make_shared<Derived>(derived());
182 :
183 1129 : result->coefficients = values;
184 :
185 4516 : for (auto &[name, quantum_numbers] : result->quantum_numbers_of_states) {
186 3387 : std::fill(quantum_numbers.begin(), quantum_numbers.end(),
187 6774 : std::numeric_limits<real_t>::max());
188 : }
189 1129 : std::fill(result->energy_of_states.begin(), result->energy_of_states.end(),
190 1129 : std::numeric_limits<real_t>::max());
191 :
192 2258 : return result;
193 1129 : }
194 :
195 : template <typename Derived>
196 5389722 : typename Basis<Derived>::real_t Basis<Derived>::get_quantum_number(const std::string &name,
197 : size_t state_index) const {
198 5389722 : auto it = quantum_numbers_of_states.find(name);
199 5389722 : if (it == quantum_numbers_of_states.end()) {
200 0 : throw std::invalid_argument("The states are not labeled by the quantum number " + name +
201 : ".");
202 : }
203 5389722 : real_t quantum_number = it->second.at(state_index);
204 5389722 : if (quantum_number == std::numeric_limits<real_t>::max()) {
205 1 : throw std::invalid_argument("The state does not have a well-defined quantum number " +
206 : name + ".");
207 : }
208 5389721 : return quantum_number;
209 : }
210 :
211 : template <typename Derived>
212 112446 : typename Basis<Derived>::real_t Basis<Derived>::get_energy(size_t state_index) const {
213 112446 : real_t energy = energy_of_states.at(state_index);
214 112446 : if (energy == std::numeric_limits<real_t>::max()) {
215 2 : throw std::invalid_argument("The state does not have a well-defined energy.");
216 : }
217 112444 : return energy;
218 : }
219 :
220 : template <typename Derived>
221 2624 : std::shared_ptr<const Derived> Basis<Derived>::get_state(size_t state_index) const {
222 : // Create a copy of the current object
223 2624 : auto restricted = std::make_shared<Derived>(derived());
224 :
225 : // Restrict the copy to the single requested state
226 2624 : restricted->coefficients = restricted->coefficients.col(state_index);
227 :
228 10496 : for (auto &[name, quantum_numbers] : restricted->quantum_numbers_of_states) {
229 7872 : quantum_numbers = {quantum_numbers_of_states.at(name)[state_index]};
230 : }
231 2624 : restricted->energy_of_states = {energy_of_states.at(state_index)};
232 :
233 5248 : return restricted;
234 2624 : }
235 :
236 : template <typename Derived>
237 : std::shared_ptr<const typename Basis<Derived>::ket_t>
238 186464 : Basis<Derived>::get_ket(size_t ket_index) const {
239 186464 : return kets[ket_index];
240 : }
241 :
242 : template <typename Derived>
243 5 : typename Basis<Derived>::Iterator Basis<Derived>::begin() const {
244 5 : return kets.begin();
245 : }
246 :
247 : template <typename Derived>
248 5 : typename Basis<Derived>::Iterator Basis<Derived>::end() const {
249 5 : return kets.end();
250 : }
251 :
252 : template <typename Derived>
253 10 : Basis<Derived>::Iterator::Iterator(typename ketvec_t::const_iterator it) : it{std::move(it)} {}
254 :
255 : template <typename Derived>
256 497 : bool Basis<Derived>::Iterator::operator!=(const Iterator &other) const {
257 497 : return other.it != it;
258 : }
259 :
260 : template <typename Derived>
261 492 : std::shared_ptr<const typename Basis<Derived>::ket_t> Basis<Derived>::Iterator::operator*() const {
262 492 : return *it;
263 : }
264 :
265 : template <typename Derived>
266 492 : typename Basis<Derived>::Iterator &Basis<Derived>::Iterator::operator++() {
267 492 : ++it;
268 492 : return *this;
269 : }
270 :
271 : template <typename Derived>
272 7626 : size_t Basis<Derived>::get_number_of_states() const {
273 7626 : return coefficients.cols();
274 : }
275 :
276 : template <typename Derived>
277 49383 : size_t Basis<Derived>::get_number_of_kets() const {
278 49383 : return coefficients.rows();
279 : }
280 :
281 : template <typename Derived>
282 : Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic>
283 6 : Basis<Derived>::get_sorter(const std::vector<SorterType> &labels) const {
284 : // Note that sorting by the energy is rejected by get_sorter_without_checks
285 :
286 : // Initialize the sorter
287 6 : Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic> sorter(coefficients.cols());
288 6 : sorter.setIdentity();
289 :
290 6 : get_sorter_without_checks(labels, sorter);
291 :
292 6 : return sorter;
293 0 : }
294 :
295 : template <typename Derived>
296 : std::vector<IndicesOfBlock>
297 7 : Basis<Derived>::get_indices_of_blocks(const std::vector<SorterType> &labels) const {
298 7 : std::set<SorterType> unique_labels(labels.begin(), labels.end());
299 7 : perform_blocks_checks(unique_labels);
300 :
301 : // Get the blocks
302 3 : IndicesOfBlocksCreator blocks_creator({0, static_cast<size_t>(coefficients.cols())});
303 3 : get_indices_of_blocks_without_checks(unique_labels, blocks_creator);
304 :
305 6 : return blocks_creator.create();
306 7 : }
307 :
308 : template <typename Derived>
309 567 : void Basis<Derived>::get_sorter_without_checks(
310 : const std::vector<SorterType> &labels,
311 : Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic> &sorter) const {
312 567 : constexpr real_t numerical_precision = 100 * std::numeric_limits<real_t>::epsilon();
313 :
314 : // Currently, since states in a basis dont store the energy, they cannot be sorted by energy.
315 : // Checking this upfront also guarantees that only labels that refer to a quantum number are
316 : // encountered below.
317 567 : if (std::find(labels.begin(), labels.end(), SorterType::ENERGY) != labels.end()) {
318 0 : throw std::invalid_argument(
319 : "States in a basis do not store the energy and thus can not be sorted by it. "
320 : "Note that sorting a system by the energy is supported nevertheless.");
321 : }
322 :
323 576 : if (coefficients.cols() == 0) {
324 0 : return;
325 : }
326 :
327 574 : std::vector<const std::vector<real_t> *> quantum_numbers_of_labels;
328 574 : quantum_numbers_of_labels.reserve(labels.size());
329 1200 : for (const auto &label : labels) {
330 621 : quantum_numbers_of_labels.push_back(&get_quantum_numbers_of_states(label));
331 : }
332 :
333 570 : int *perm_begin = sorter.indices().data();
334 564 : int *perm_end = perm_begin + coefficients.cols();
335 :
336 : // Sort the vector based on the requested labels
337 569 : set_task_status("Sorting basis states...");
338 146582 : std::stable_sort(perm_begin, perm_end, [&](int a, int b) {
339 264350 : for (const auto *quantum_numbers : quantum_numbers_of_labels) {
340 169636 : if (std::abs((*quantum_numbers)[a] - (*quantum_numbers)[b]) > numerical_precision) {
341 54768 : return (*quantum_numbers)[a] < (*quantum_numbers)[b];
342 : }
343 : }
344 93528 : return false; // Elements are equal
345 : });
346 605 : }
347 :
348 : template <typename Derived>
349 568 : void Basis<Derived>::get_indices_of_blocks_without_checks(
350 : const std::set<SorterType> &unique_labels, IndicesOfBlocksCreator &blocks_creator) const {
351 568 : constexpr real_t numerical_precision = 100 * std::numeric_limits<real_t>::epsilon();
352 :
353 568 : if (coefficients.cols() == 0) {
354 0 : return;
355 : }
356 :
357 560 : std::vector<const std::vector<real_t> *> quantum_numbers_of_labels;
358 560 : quantum_numbers_of_labels.reserve(unique_labels.size());
359 1208 : for (const auto &label : unique_labels) {
360 625 : quantum_numbers_of_labels.push_back(&get_quantum_numbers_of_states(label));
361 : }
362 :
363 : // A new block starts whenever a state differs from its predecessor in one of the labels
364 587 : set_task_status("Identifying basis blocks...");
365 31657 : for (int i = 1; i < coefficients.cols(); ++i) {
366 66987 : for (const auto *quantum_numbers : quantum_numbers_of_labels) {
367 36467 : if (std::abs((*quantum_numbers)[i] - (*quantum_numbers)[i - 1]) > numerical_precision) {
368 736 : blocks_creator.add(i);
369 735 : break;
370 : }
371 : }
372 : }
373 548 : }
374 :
375 : template <typename Derived>
376 168 : std::shared_ptr<const Derived> Basis<Derived>::canonicalized() const {
377 168 : auto result = std::make_shared<Derived>(derived());
378 168 : size_t n = kets.size();
379 168 : result->coefficients.resize(n, n);
380 168 : result->make_canonical();
381 336 : return result;
382 168 : }
383 :
384 : template <typename Derived>
385 112 : bool Basis<Derived>::is_canonical() const {
386 112 : constexpr real_t numerical_precision = 100 * std::numeric_limits<real_t>::epsilon();
387 :
388 : // The basis is canonical if and only if its coefficient matrix is the identity matrix, i.e.,
389 : // if the i-th state is the i-th ket
390 112 : if (coefficients.rows() != coefficients.cols()) {
391 0 : return false;
392 : }
393 :
394 112 : Eigen::Index num_ones = 0;
395 2064 : for (Eigen::Index row = 0; row < coefficients.outerSize(); ++row) {
396 1954 : for (typename Eigen::SparseMatrix<scalar_t, Eigen::RowMajor>::InnerIterator it(coefficients,
397 : row);
398 3914 : it; ++it) {
399 3915 : if (it.row() == it.col() &&
400 2923 : std::abs(it.value() - static_cast<scalar_t>(1)) <= numerical_precision) {
401 1952 : ++num_ones;
402 10 : } else if (std::abs(it.value()) > numerical_precision) {
403 2 : return false;
404 : }
405 : }
406 : }
407 110 : return num_ones == coefficients.rows();
408 : }
409 :
410 : template <typename Derived>
411 1256 : std::shared_ptr<const Derived> Basis<Derived>::transformed(
412 : const Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic> &sorter) const {
413 : // Create a copy of the current object
414 1256 : auto transformed = std::make_shared<Derived>(derived());
415 :
416 1260 : if (coefficients.cols() == 0) {
417 0 : return transformed;
418 : }
419 :
420 : // Apply the sorting
421 1248 : set_task_status("Applying basis sorting...");
422 1264 : transformed->coefficients = coefficients * sorter;
423 :
424 1258 : set_task_status("Relabeling sorted basis states...");
425 4828 : for (auto &[name, transformed_quantum_numbers] : transformed->quantum_numbers_of_states) {
426 3784 : const auto &quantum_numbers = quantum_numbers_of_states.at(name);
427 3761 : transformed_quantum_numbers.resize(sorter.size());
428 193312 : for (int i = 0; i < sorter.size(); ++i) {
429 189183 : transformed_quantum_numbers[i] = quantum_numbers[sorter.indices()[i]];
430 : }
431 : }
432 1260 : transformed->energy_of_states.resize(sorter.size());
433 65095 : for (int i = 0; i < sorter.size(); ++i) {
434 63635 : transformed->energy_of_states[i] = energy_of_states[sorter.indices()[i]];
435 : }
436 :
437 1235 : return transformed;
438 1261 : }
439 :
440 : template <typename Derived>
441 : std::shared_ptr<const Derived>
442 916 : Basis<Derived>::transformed(const Eigen::SparseMatrix<scalar_t, Eigen::RowMajor> &transformation,
443 : const std::vector<real_t> &energy_of_transformed_states) const {
444 1790 : if (!energy_of_transformed_states.empty() &&
445 884 : energy_of_transformed_states.size() != static_cast<size_t>(transformation.cols())) {
446 1 : throw std::invalid_argument(
447 : "The number of energies must match the number of transformed states.");
448 : }
449 :
450 913 : const real_t numerical_precision =
451 897 : 100 * std::sqrt(coefficients.rows()) * std::numeric_limits<real_t>::epsilon();
452 :
453 : // Create a copy of the current object
454 913 : auto transformed = std::make_shared<Derived>(derived());
455 :
456 924 : if (coefficients.cols() == 0) {
457 0 : return transformed;
458 : }
459 :
460 : // Apply the transformation
461 922 : set_task_status("Applying basis transformation...");
462 926 : transformed->coefficients = coefficients * transformation;
463 :
464 925 : Eigen::SparseMatrix<real_t> probs = transformation.cwiseAbs2().transpose();
465 : // Drop explicitly stored zeros
466 920 : probs.prune(
467 902887 : [](Eigen::Index /*row*/, Eigen::Index /*col*/, const real_t &value) { return value != 0; });
468 925 : Eigen::VectorX<real_t> norm = probs * Eigen::VectorX<real_t>::Ones(probs.cols());
469 :
470 920 : set_task_status("Updating transformed quantum numbers...");
471 3656 : for (auto &[name, transformed_quantum_numbers] : transformed->quantum_numbers_of_states) {
472 2753 : const auto &quantum_numbers = quantum_numbers_of_states.at(name);
473 5508 : auto map = Eigen::Map<const Eigen::VectorX<real_t>>(
474 2750 : quantum_numbers.data(), static_cast<Eigen::Index>(quantum_numbers.size()));
475 2734 : Eigen::VectorX<real_t> val = probs * map;
476 2776 : Eigen::VectorX<real_t> sq = probs * map.cwiseAbs2();
477 2771 : transformed_quantum_numbers.resize(probs.rows());
478 :
479 284180 : for (size_t i = 0; i < transformed_quantum_numbers.size(); ++i) {
480 281463 : const real_t mean = val[i] / norm[i];
481 281240 : const real_t mean_sq = sq[i] / norm[i];
482 281287 : const real_t variance = std::abs(mean_sq - mean * mean);
483 281506 : if (variance < numerical_precision * std::max(static_cast<real_t>(1), mean_sq)) {
484 157249 : transformed_quantum_numbers[i] = std::round(mean * 2) / 2;
485 : } else {
486 124478 : transformed_quantum_numbers[i] = std::numeric_limits<real_t>::max();
487 : }
488 : }
489 : }
490 :
491 : // The transformed states only have well-defined energies if they are explicitly specified,
492 : // e.g., as eigenenergies
493 923 : if (energy_of_transformed_states.empty()) {
494 48 : transformed->energy_of_states.assign(static_cast<size_t>(transformation.cols()),
495 48 : std::numeric_limits<real_t>::max());
496 : } else {
497 897 : transformed->energy_of_states = energy_of_transformed_states;
498 : }
499 :
500 909 : return transformed;
501 918 : }
502 :
503 : // Explicit instantiations
504 : template class Basis<BasisAtom<double>>;
505 : template class Basis<BasisAtom<std::complex<double>>>;
506 : template class Basis<BasisPair<double>>;
507 : template class Basis<BasisPair<std::complex<double>>>;
508 : } // namespace pairinteraction
|
