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/Parity.hpp"
       9             : #include "pairinteraction/enums/TransformationType.hpp"
      10             : #include "pairinteraction/ket/KetAtom.hpp"
      11             : #include "pairinteraction/ket/KetPair.hpp"
      12             : #include "pairinteraction/utils/TaskControl.hpp"
      13             : #include "pairinteraction/utils/eigen_assertion.hpp"
      14             : #include "pairinteraction/utils/eigen_compat.hpp"
      15             : 
      16             : #include <cassert>
      17             : #include <set>
      18             : 
      19             : namespace pairinteraction {
      20             : 
      21             : template <typename Scalar>
      22             : class BasisAtom;
      23             : 
      24             : template <typename Derived>
      25        1824 : void Basis<Derived>::perform_sorter_checks(const std::vector<TransformationType> &labels) const {
      26             :     // Check if the labels are valid sorting labels
      27        3658 :     for (const auto &label : labels) {
      28        1844 :         if (label == TransformationType::CANONICAL_ORDER) {
      29           0 :             throw std::invalid_argument("Sorting by canonical order is not supported.");
      30             :         }
      31        1844 :         if (!utils::is_sorting(label)) {
      32           0 :             throw std::invalid_argument("One of the labels is not a valid sorting label.");
      33             :         }
      34             :     }
      35        1818 : }
      36             : 
      37             : template <typename Derived>
      38         624 : void Basis<Derived>::perform_blocks_checks(
      39             :     const std::set<TransformationType> &unique_labels) const {
      40             :     // Check if the states are sorted by the requested labels
      41         624 :     std::set<TransformationType> unique_labels_present;
      42        1233 :     for (const auto &label : get_transformation().transformation_type) {
      43         653 :         if (!utils::is_sorting(label) || unique_labels_present.size() >= unique_labels.size()) {
      44          51 :             break;
      45             :         }
      46         607 :         unique_labels_present.insert(label);
      47             :     }
      48         625 :     if (unique_labels != unique_labels_present) {
      49           0 :         throw std::invalid_argument("The states are not sorted by the requested labels.");
      50             :     }
      51             : 
      52             :     // Throw a meaningful error if getting the blocks by energy is requested as this might be a
      53             :     // common mistake
      54         622 :     if (unique_labels.contains(TransformationType::SORT_BY_ENERGY)) {
      55           0 :         throw std::invalid_argument("States do not store the energy and thus no energy blocks can "
      56             :                                     "be obtained. Use an energy operator instead.");
      57             :     }
      58         617 : }
      59             : 
      60             : template <typename Derived>
      61        2054 : Basis<Derived>::Basis(ketvec_t &&kets)
      62        4108 :     : kets(std::move(kets)), coefficients{{static_cast<Eigen::Index>(this->kets.size()),
      63        2054 :                                            static_cast<Eigen::Index>(this->kets.size())},
      64        4108 :                                           {TransformationType::CANONICAL_ORDER}} {
      65        2054 :     if (this->kets.empty()) {
      66           0 :         throw std::invalid_argument("The basis must contain at least one element.");
      67             :     }
      68        2054 :     state_index_to_quantum_number_f.reserve(this->kets.size());
      69        2054 :     state_index_to_quantum_number_m.reserve(this->kets.size());
      70        2054 :     state_index_to_parity.reserve(this->kets.size());
      71     2611496 :     for (const auto &ket : this->kets) {
      72     2609442 :         real_t f = std::numeric_limits<real_t>::max();
      73     2609442 :         real_t m = std::numeric_limits<real_t>::max();
      74     2609442 :         Parity p = Parity::UNKNOWN;
      75             :         // TODO: this is a workaround, and should be fixed, once we restructure the quantum number
      76             :         // handling of the Basis class
      77             :         if constexpr (requires { ket->get_quantum_number(std::string{}); }) {
      78       42793 :             f = ket->get_quantum_number("f");
      79       42793 :             m = ket->get_quantum_number("m");
      80       42793 :             p = static_cast<Parity>(static_cast<int>(ket->get_quantum_number("parity")));
      81             :         } else {
      82     2566649 :             m = ket->get_quantum_number_m();
      83             :         }
      84     2609442 :         state_index_to_quantum_number_f.push_back(f);
      85     2609442 :         state_index_to_quantum_number_m.push_back(m);
      86     2609442 :         state_index_to_parity.push_back(p);
      87     2609442 :         if (f == std::numeric_limits<real_t>::max()) {
      88     2566649 :             _has_quantum_number_f = false;
      89             :         }
      90     2609442 :         if (m == std::numeric_limits<real_t>::max()) {
      91           0 :             _has_quantum_number_m = false;
      92             :         }
      93     2609442 :         if (p == Parity::UNKNOWN) {
      94     2566649 :             _has_parity = false;
      95             :         }
      96             :     }
      97        2054 :     coefficients.matrix.setIdentity();
      98        2054 : }
      99             : 
     100             : template <typename Derived>
     101        1217 : bool Basis<Derived>::has_quantum_number_f() const {
     102        1217 :     return _has_quantum_number_f;
     103             : }
     104             : 
     105             : template <typename Derived>
     106     5262703 : bool Basis<Derived>::has_quantum_number_m() const {
     107     5262703 :     return _has_quantum_number_m;
     108             : }
     109             : 
     110             : template <typename Derived>
     111        1224 : bool Basis<Derived>::has_parity() const {
     112        1224 :     return _has_parity;
     113             : }
     114             : 
     115             : template <typename Derived>
     116        5109 : const Derived &Basis<Derived>::derived() const {
     117        5109 :     return static_cast<const Derived &>(*this);
     118             : }
     119             : 
     120             : template <typename Derived>
     121        1189 : const typename Basis<Derived>::ketvec_t &Basis<Derived>::get_kets() const {
     122        1189 :     return kets;
     123             : }
     124             : 
     125             : template <typename Derived>
     126             : const Eigen::SparseMatrix<typename Basis<Derived>::scalar_t, Eigen::RowMajor> &
     127       14421 : Basis<Derived>::get_coefficients() const {
     128       14421 :     return coefficients.matrix;
     129             : }
     130             : 
     131             : template <typename Derived>
     132        1060 : std::shared_ptr<const Derived> Basis<Derived>::copy_with_coefficients(
     133             :     const Eigen::SparseMatrix<scalar_t, Eigen::RowMajor> &values) const {
     134        1060 :     if (values.rows() != coefficients.matrix.rows()) {
     135           0 :         throw std::invalid_argument("Incompatible number of rows.");
     136             :     }
     137        1060 :     if (values.cols() != coefficients.matrix.cols()) {
     138           0 :         throw std::invalid_argument("Incompatible number of columns.");
     139             :     }
     140             : 
     141             :     // Create a copy of the current object and update the coefficients of the copy
     142        1060 :     auto result = std::make_shared<Derived>(derived());
     143             : 
     144        1060 :     result->coefficients.matrix = values;
     145        1060 :     result->coefficients.transformation_type = {TransformationType::ARBITRARY};
     146             : 
     147        1060 :     std::fill(result->state_index_to_quantum_number_f.begin(),
     148        1060 :               result->state_index_to_quantum_number_f.end(), std::numeric_limits<real_t>::max());
     149        1060 :     std::fill(result->state_index_to_quantum_number_m.begin(),
     150        1060 :               result->state_index_to_quantum_number_m.end(), std::numeric_limits<real_t>::max());
     151        1060 :     std::fill(result->state_index_to_parity.begin(), result->state_index_to_parity.end(),
     152        1060 :               Parity::UNKNOWN);
     153        1060 :     result->_has_quantum_number_f = false;
     154        1060 :     result->_has_quantum_number_m = false;
     155        1060 :     result->_has_parity = false;
     156             : 
     157        2120 :     return result;
     158        1060 : }
     159             : 
     160             : template <typename Derived>
     161           0 : typename Basis<Derived>::real_t Basis<Derived>::get_quantum_number_f(size_t state_index) const {
     162           0 :     real_t quantum_number_f = state_index_to_quantum_number_f.at(state_index);
     163           0 :     if (quantum_number_f == std::numeric_limits<real_t>::max()) {
     164           0 :         throw std::invalid_argument("The state does not have a well-defined quantum number f.");
     165             :     }
     166           0 :     return quantum_number_f;
     167             : }
     168             : 
     169             : template <typename Derived>
     170     5261526 : typename Basis<Derived>::real_t Basis<Derived>::get_quantum_number_m(size_t state_index) const {
     171     5261526 :     real_t quantum_number_m = state_index_to_quantum_number_m.at(state_index);
     172     5261526 :     if (quantum_number_m == std::numeric_limits<real_t>::max()) {
     173           0 :         throw std::invalid_argument("The state does not have a well-defined quantum number m.");
     174             :     }
     175     5261526 :     return quantum_number_m;
     176             : }
     177             : 
     178             : template <typename Derived>
     179       32795 : Parity Basis<Derived>::get_parity(size_t state_index) const {
     180       32795 :     Parity parity = state_index_to_parity.at(state_index);
     181       32795 :     if (parity == Parity::UNKNOWN) {
     182           0 :         throw std::invalid_argument("The state does not have a well-defined parity.");
     183             :     }
     184       32795 :     return parity;
     185             : }
     186             : 
     187             : template <typename Derived>
     188        2028 : std::shared_ptr<const Derived> Basis<Derived>::get_state(size_t state_index) const {
     189             :     // Create a copy of the current object
     190        2028 :     auto restricted = std::make_shared<Derived>(derived());
     191             : 
     192             :     // Restrict the copy to the single requested state
     193        2028 :     restricted->coefficients.matrix = restricted->coefficients.matrix.col(state_index);
     194             : 
     195        2028 :     restricted->state_index_to_quantum_number_f = {state_index_to_quantum_number_f[state_index]};
     196        2028 :     restricted->state_index_to_quantum_number_m = {state_index_to_quantum_number_m[state_index]};
     197        2028 :     restricted->state_index_to_parity = {state_index_to_parity[state_index]};
     198             : 
     199        4056 :     restricted->_has_quantum_number_f =
     200        2028 :         restricted->state_index_to_quantum_number_f[0] != std::numeric_limits<real_t>::max();
     201        4056 :     restricted->_has_quantum_number_m =
     202        2028 :         restricted->state_index_to_quantum_number_m[0] != std::numeric_limits<real_t>::max();
     203        2028 :     restricted->_has_parity = restricted->state_index_to_parity[0] != Parity::UNKNOWN;
     204             : 
     205        4056 :     return restricted;
     206        2028 : }
     207             : 
     208             : template <typename Derived>
     209             : std::shared_ptr<const typename Basis<Derived>::ket_t>
     210      172252 : Basis<Derived>::get_ket(size_t ket_index) const {
     211      172252 :     return kets[ket_index];
     212             : }
     213             : 
     214             : template <typename Derived>
     215           4 : typename Basis<Derived>::Iterator Basis<Derived>::begin() const {
     216           4 :     return kets.begin();
     217             : }
     218             : 
     219             : template <typename Derived>
     220           4 : typename Basis<Derived>::Iterator Basis<Derived>::end() const {
     221           4 :     return kets.end();
     222             : }
     223             : 
     224             : template <typename Derived>
     225           8 : Basis<Derived>::Iterator::Iterator(typename ketvec_t::const_iterator it) : it{std::move(it)} {}
     226             : 
     227             : template <typename Derived>
     228         240 : bool Basis<Derived>::Iterator::operator!=(const Iterator &other) const {
     229         240 :     return other.it != it;
     230             : }
     231             : 
     232             : template <typename Derived>
     233         236 : std::shared_ptr<const typename Basis<Derived>::ket_t> Basis<Derived>::Iterator::operator*() const {
     234         236 :     return *it;
     235             : }
     236             : 
     237             : template <typename Derived>
     238         236 : typename Basis<Derived>::Iterator &Basis<Derived>::Iterator::operator++() {
     239         236 :     ++it;
     240         236 :     return *this;
     241             : }
     242             : 
     243             : template <typename Derived>
     244        6602 : size_t Basis<Derived>::get_number_of_states() const {
     245        6602 :     return coefficients.matrix.cols();
     246             : }
     247             : 
     248             : template <typename Derived>
     249       48376 : size_t Basis<Derived>::get_number_of_kets() const {
     250       48376 :     return coefficients.matrix.rows();
     251             : }
     252             : 
     253             : template <typename Derived>
     254             : const Transformation<typename Basis<Derived>::scalar_t> &
     255         729 : Basis<Derived>::get_transformation() const {
     256         729 :     return coefficients;
     257             : }
     258             : 
     259             : template <typename Derived>
     260           1 : Sorting Basis<Derived>::get_sorter(const std::vector<TransformationType> &labels) const {
     261           1 :     perform_sorter_checks(labels);
     262             : 
     263             :     // Throw a meaningful error if sorting by energy is requested as this might be a common mistake
     264           1 :     if (std::find(labels.begin(), labels.end(), TransformationType::SORT_BY_ENERGY) !=
     265           2 :         labels.end()) {
     266           0 :         throw std::invalid_argument("States do not store the energy and thus can not be sorted by "
     267             :                                     "the energy. Use an energy operator instead.");
     268             :     }
     269             : 
     270             :     // Initialize transformation
     271           1 :     Sorting transformation;
     272           1 :     transformation.matrix.resize(coefficients.matrix.cols());
     273           1 :     transformation.matrix.setIdentity();
     274             : 
     275             :     // Get the sorter
     276           1 :     get_sorter_without_checks(labels, transformation);
     277             : 
     278             :     // Check if all labels have been used for sorting
     279           1 :     if (labels != transformation.transformation_type) {
     280           0 :         throw std::invalid_argument("The states could not be sorted by all the requested labels.");
     281             :     }
     282             : 
     283           1 :     return transformation;
     284           0 : }
     285             : 
     286             : template <typename Derived>
     287             : std::vector<IndicesOfBlock>
     288           1 : Basis<Derived>::get_indices_of_blocks(const std::vector<TransformationType> &labels) const {
     289           1 :     perform_sorter_checks(labels);
     290             : 
     291           1 :     std::set<TransformationType> unique_labels(labels.begin(), labels.end());
     292           1 :     perform_blocks_checks(unique_labels);
     293             : 
     294             :     // Get the blocks
     295           1 :     IndicesOfBlocksCreator blocks_creator({0, static_cast<size_t>(coefficients.matrix.cols())});
     296           1 :     get_indices_of_blocks_without_checks(unique_labels, blocks_creator);
     297             : 
     298           2 :     return blocks_creator.create();
     299           1 : }
     300             : 
     301             : template <typename Derived>
     302         573 : void Basis<Derived>::get_sorter_without_checks(const std::vector<TransformationType> &labels,
     303             :                                                Sorting &transformation) const {
     304         573 :     constexpr real_t numerical_precision = 100 * std::numeric_limits<real_t>::epsilon();
     305             : 
     306         573 :     int *perm_begin = transformation.matrix.indices().data();
     307         578 :     int *perm_end = perm_begin + coefficients.matrix.cols();
     308         570 :     const int *perm_back = perm_end - 1;
     309             : 
     310             :     // Sort the vector based on the requested labels
     311         570 :     set_task_status("Sorting basis states...");
     312      609110 :     std::stable_sort(perm_begin, perm_end, [&](int a, int b) {
     313      268683 :         for (const auto &label : labels) {
     314      173543 :             switch (label) {
     315       27028 :             case TransformationType::SORT_BY_PARITY:
     316       27028 :                 if (state_index_to_parity[a] != state_index_to_parity[b]) {
     317       55156 :                     return state_index_to_parity[a] < state_index_to_parity[b];
     318             :                 }
     319       17851 :                 break;
     320      148151 :             case TransformationType::SORT_BY_QUANTUM_NUMBER_M:
     321      296239 :                 if (std::abs(state_index_to_quantum_number_m[a] -
     322      296671 :                              state_index_to_quantum_number_m[b]) > numerical_precision) {
     323       46049 :                     return state_index_to_quantum_number_m[a] < state_index_to_quantum_number_m[b];
     324             :                 }
     325      103015 :                 break;
     326           0 :             case TransformationType::SORT_BY_QUANTUM_NUMBER_F:
     327           0 :                 if (std::abs(state_index_to_quantum_number_f[a] -
     328           0 :                              state_index_to_quantum_number_f[b]) > numerical_precision) {
     329           0 :                     return state_index_to_quantum_number_f[a] < state_index_to_quantum_number_f[b];
     330             :                 }
     331           0 :                 break;
     332           0 :             default:
     333           0 :                 std::abort(); // Can't happen because of previous checks
     334             :             }
     335             :         }
     336       94867 :         return false; // Elements are equal
     337             :     });
     338             : 
     339             :     // Check for invalid values and add transformation types
     340        1186 :     for (const auto &label : labels) {
     341         611 :         switch (label) {
     342          41 :         case TransformationType::SORT_BY_PARITY:
     343          41 :             if (state_index_to_parity[*perm_back] == Parity::UNKNOWN) {
     344           0 :                 throw std::invalid_argument(
     345             :                     "States cannot be labeled and thus not sorted by the parity.");
     346             :             }
     347          41 :             transformation.transformation_type.push_back(TransformationType::SORT_BY_PARITY);
     348          41 :             break;
     349         571 :         case TransformationType::SORT_BY_QUANTUM_NUMBER_M:
     350         571 :             if (state_index_to_quantum_number_m[*perm_back] == std::numeric_limits<real_t>::max()) {
     351           0 :                 throw std::invalid_argument(
     352             :                     "States cannot be labeled and thus not sorted by the quantum number m.");
     353             :             }
     354         575 :             transformation.transformation_type.push_back(
     355         575 :                 TransformationType::SORT_BY_QUANTUM_NUMBER_M);
     356         569 :             break;
     357           0 :         case TransformationType::SORT_BY_QUANTUM_NUMBER_F:
     358           0 :             if (state_index_to_quantum_number_f[*perm_back] == std::numeric_limits<real_t>::max()) {
     359           0 :                 throw std::invalid_argument(
     360             :                     "States cannot be labeled and thus not sorted by the quantum number f.");
     361             :             }
     362           0 :             transformation.transformation_type.push_back(
     363           0 :                 TransformationType::SORT_BY_QUANTUM_NUMBER_F);
     364           0 :             break;
     365           2 :         default:
     366           2 :             std::abort(); // Can't happen because of previous checks
     367             :         }
     368             :     }
     369         574 : }
     370             : 
     371             : template <typename Derived>
     372         572 : void Basis<Derived>::get_indices_of_blocks_without_checks(
     373             :     const std::set<TransformationType> &unique_labels,
     374             :     IndicesOfBlocksCreator &blocks_creator) const {
     375         572 :     constexpr real_t numerical_precision = 100 * std::numeric_limits<real_t>::epsilon();
     376             : 
     377         572 :     auto last_quantum_number_f = state_index_to_quantum_number_f[0];
     378         563 :     auto last_quantum_number_m = state_index_to_quantum_number_m[0];
     379         566 :     auto last_parity = state_index_to_parity[0];
     380             : 
     381         568 :     set_task_status("Identifying basis blocks...");
     382       32530 :     for (int i = 0; i < coefficients.matrix.cols(); ++i) {
     383       68998 :         for (auto label : unique_labels) {
     384       37792 :             if (label == TransformationType::SORT_BY_QUANTUM_NUMBER_F &&
     385           0 :                 std::abs(state_index_to_quantum_number_f[i] - last_quantum_number_f) >
     386             :                     numerical_precision) {
     387           0 :                 blocks_creator.add(i);
     388           0 :                 break;
     389             :             }
     390       69464 :             if (label == TransformationType::SORT_BY_QUANTUM_NUMBER_M &&
     391       31723 :                 std::abs(state_index_to_quantum_number_m[i] - last_quantum_number_m) >
     392             :                     numerical_precision) {
     393         569 :                 blocks_creator.add(i);
     394         565 :                 break;
     395             :             }
     396       43493 :             if (label == TransformationType::SORT_BY_PARITY &&
     397        6326 :                 state_index_to_parity[i] != last_parity) {
     398         144 :                 blocks_creator.add(i);
     399         143 :                 break;
     400             :             }
     401             :         }
     402       31064 :         last_quantum_number_f = state_index_to_quantum_number_f[i];
     403       31857 :         last_quantum_number_m = state_index_to_quantum_number_m[i];
     404       31945 :         last_parity = state_index_to_parity[i];
     405             :     }
     406         543 : }
     407             : 
     408             : template <typename Derived>
     409         168 : std::shared_ptr<const Derived> Basis<Derived>::canonicalized() const {
     410         168 :     auto result = std::make_shared<Derived>(derived());
     411             : 
     412         168 :     size_t n = kets.size();
     413             : 
     414         168 :     result->coefficients.matrix.resize(n, n);
     415         168 :     result->coefficients.matrix.setIdentity();
     416         168 :     result->coefficients.transformation_type = {TransformationType::CANONICAL_ORDER};
     417             : 
     418         168 :     result->state_index_to_quantum_number_f.resize(n);
     419         168 :     result->state_index_to_quantum_number_m.resize(n);
     420         168 :     result->state_index_to_parity.resize(n);
     421         168 :     result->_has_quantum_number_f = true;
     422         168 :     result->_has_quantum_number_m = true;
     423         168 :     result->_has_parity = true;
     424             : 
     425       15430 :     for (size_t i = 0; i < n; ++i) {
     426       15262 :         real_t f = std::numeric_limits<real_t>::max();
     427       15262 :         real_t m = std::numeric_limits<real_t>::max();
     428       15262 :         Parity p = Parity::UNKNOWN;
     429             :         // TODO: this is a workaround, and should be fixed, once we restructure the quantum number
     430             :         // handling of the Basis class
     431             :         if constexpr (requires { kets[i]->get_quantum_number(std::string{}); }) {
     432       14558 :             f = kets[i]->get_quantum_number("f");
     433       14558 :             m = kets[i]->get_quantum_number("m");
     434       14558 :             p = static_cast<Parity>(static_cast<int>(kets[i]->get_quantum_number("parity")));
     435             :         } else {
     436         704 :             m = kets[i]->get_quantum_number_m();
     437             :         }
     438       15262 :         result->state_index_to_quantum_number_f[i] = f;
     439       15262 :         result->state_index_to_quantum_number_m[i] = m;
     440       15262 :         result->state_index_to_parity[i] = p;
     441       15262 :         if (f == std::numeric_limits<real_t>::max()) {
     442         704 :             result->_has_quantum_number_f = false;
     443             :         }
     444       15262 :         if (m == std::numeric_limits<real_t>::max()) {
     445           0 :             result->_has_quantum_number_m = false;
     446             :         }
     447       15262 :         if (p == Parity::UNKNOWN) {
     448         704 :             result->_has_parity = false;
     449             :         }
     450             :     }
     451             : 
     452         336 :     return result;
     453         168 : }
     454             : 
     455             : template <typename Derived>
     456         108 : bool Basis<Derived>::is_canonical() const {
     457         108 :     return get_transformation().transformation_type ==
     458         216 :         std::vector<TransformationType>{TransformationType::CANONICAL_ORDER};
     459             : }
     460             : 
     461             : template <typename Derived>
     462        1210 : std::shared_ptr<const Derived> Basis<Derived>::transformed(const Sorting &transformation) const {
     463             :     // Create a copy of the current object
     464        1210 :     auto transformed = std::make_shared<Derived>(derived());
     465             : 
     466        1209 :     if (coefficients.matrix.cols() == 0) {
     467           0 :         return transformed;
     468             :     }
     469             : 
     470             :     // Apply the transformation
     471        1192 :     set_task_status("Applying basis sorting...");
     472        1212 :     transformed->coefficients.matrix = coefficients.matrix * transformation.matrix;
     473        1210 :     transformed->coefficients.transformation_type = transformation.transformation_type;
     474             : 
     475        1199 :     transformed->state_index_to_quantum_number_f.resize(transformation.matrix.size());
     476        1205 :     transformed->state_index_to_quantum_number_m.resize(transformation.matrix.size());
     477        1207 :     transformed->state_index_to_parity.resize(transformation.matrix.size());
     478             : 
     479        1206 :     set_task_status("Relabeling sorted basis states...");
     480       63272 :     for (int i = 0; i < transformation.matrix.size(); ++i) {
     481       61985 :         transformed->state_index_to_quantum_number_f[i] =
     482       61986 :             state_index_to_quantum_number_f[transformation.matrix.indices()[i]];
     483       61953 :         transformed->state_index_to_quantum_number_m[i] =
     484       62077 :             state_index_to_quantum_number_m[transformation.matrix.indices()[i]];
     485       62072 :         transformed->state_index_to_parity[i] =
     486       62059 :             state_index_to_parity[transformation.matrix.indices()[i]];
     487             :     }
     488             : 
     489        1174 :     return transformed;
     490        1208 : }
     491             : 
     492             : template <typename Derived>
     493             : std::shared_ptr<const Derived>
     494         649 : Basis<Derived>::transformed(const Transformation<scalar_t> &transformation) const {
     495             :     // TODO why is "numerical_precision = 100 * std::sqrt(coefficients.matrix.rows()) *
     496             :     // std::numeric_limits<real_t>::epsilon()" too small for figuring out whether m is conserved?
     497         649 :     real_t numerical_precision = 0.001;
     498             : 
     499             :     // Create a copy of the current object
     500         649 :     auto transformed = std::make_shared<Derived>(derived());
     501             : 
     502         645 :     if (coefficients.matrix.cols() == 0) {
     503           0 :         return transformed;
     504             :     }
     505             : 
     506             :     // Apply the transformation
     507             :     // If a quantum number turns out to be conserved by the transformation, it will be
     508             :     // rounded to the nearest half integer to avoid loss of numerical_precision.
     509         643 :     set_task_status("Applying basis transformation...");
     510         649 :     transformed->coefficients.matrix = coefficients.matrix * transformation.matrix;
     511         649 :     transformed->coefficients.transformation_type = transformation.transformation_type;
     512             : 
     513         649 :     Eigen::SparseMatrix<real_t> probs = transformation.matrix.cwiseAbs2().transpose();
     514             : 
     515         649 :     set_task_status("Updating transformed quantum numbers...");
     516             :     {
     517        1297 :         auto map = Eigen::Map<const Eigen::VectorX<real_t>>(state_index_to_quantum_number_f.data(),
     518         648 :                                                             state_index_to_quantum_number_f.size());
     519         646 :         Eigen::VectorX<real_t> val = probs * map;
     520         648 :         Eigen::VectorX<real_t> sq = probs * map.cwiseAbs2();
     521         646 :         Eigen::VectorX<real_t> diff = (val.cwiseAbs2() - sq).cwiseAbs();
     522         642 :         transformed->state_index_to_quantum_number_f.resize(probs.rows());
     523             : 
     524       63363 :         for (size_t i = 0; i < transformed->state_index_to_quantum_number_f.size(); ++i) {
     525       62763 :             if (diff[i] < numerical_precision) {
     526        4898 :                 transformed->state_index_to_quantum_number_f[i] = std::round(val[i] * 2) / 2;
     527             :             } else {
     528       57874 :                 transformed->state_index_to_quantum_number_f[i] =
     529       57829 :                     std::numeric_limits<real_t>::max();
     530       57883 :                 transformed->_has_quantum_number_f = false;
     531             :             }
     532             :         }
     533         645 :     }
     534             : 
     535             :     {
     536        1294 :         auto map = Eigen::Map<const Eigen::VectorX<real_t>>(state_index_to_quantum_number_m.data(),
     537         646 :                                                             state_index_to_quantum_number_m.size());
     538         645 :         Eigen::VectorX<real_t> val = probs * map;
     539         648 :         Eigen::VectorX<real_t> sq = probs * map.cwiseAbs2();
     540         649 :         Eigen::VectorX<real_t> diff = (val.cwiseAbs2() - sq).cwiseAbs();
     541         647 :         transformed->state_index_to_quantum_number_m.resize(probs.rows());
     542             : 
     543       63483 :         for (size_t i = 0; i < transformed->state_index_to_quantum_number_m.size(); ++i) {
     544       62895 :             if (diff[i] < numerical_precision) {
     545       60213 :                 transformed->state_index_to_quantum_number_m[i] = std::round(val[i] * 2) / 2;
     546             :             } else {
     547        2646 :                 transformed->state_index_to_quantum_number_m[i] =
     548        2640 :                     std::numeric_limits<real_t>::max();
     549        2629 :                 transformed->_has_quantum_number_m = false;
     550             :             }
     551             :         }
     552         649 :     }
     553             : 
     554             :     {
     555             :         using utype = std::underlying_type<Parity>::type;
     556         647 :         Eigen::VectorX<real_t> map(state_index_to_parity.size());
     557      101344 :         for (size_t i = 0; i < state_index_to_parity.size(); ++i) {
     558      100749 :             map[i] = static_cast<utype>(state_index_to_parity[i]);
     559             :         }
     560         612 :         Eigen::VectorX<real_t> val = probs * map;
     561         649 :         Eigen::VectorX<real_t> sq = probs * map.cwiseAbs2();
     562         647 :         Eigen::VectorX<real_t> diff = (val.cwiseAbs2() - sq).cwiseAbs();
     563         648 :         transformed->state_index_to_parity.resize(probs.rows());
     564             : 
     565       63397 :         for (size_t i = 0; i < transformed->state_index_to_parity.size(); ++i) {
     566       62789 :             if (diff[i] < numerical_precision) {
     567       56652 :                 transformed->state_index_to_parity[i] = static_cast<Parity>(std::lround(val[i]));
     568             :             } else {
     569        6027 :                 transformed->state_index_to_parity[i] = Parity::UNKNOWN;
     570        6078 :                 transformed->_has_parity = false;
     571             :             }
     572             :         }
     573         649 :     }
     574             : 
     575         648 :     return transformed;
     576         645 : }
     577             : 
     578             : // Explicit instantiations
     579             : template class Basis<BasisAtom<double>>;
     580             : template class Basis<BasisAtom<std::complex<double>>>;
     581             : template class Basis<BasisPair<double>>;
     582             : template class Basis<BasisPair<std::complex<double>>>;
     583             : } // namespace pairinteraction