.. _program_listing_file_SeQuant_domain_mbpt_biorthogonalization.cpp:

Program Listing for File biorthogonalization.cpp
================================================

|exhale_lsh| :ref:`Return to documentation for file <file_SeQuant_domain_mbpt_biorthogonalization.cpp>` (``SeQuant/domain/mbpt/biorthogonalization.cpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   #include <SeQuant/domain/mbpt/biorthogonalization.hpp>
   
   #include <SeQuant/core/container.hpp>
   #include <SeQuant/core/expr.hpp>
   #include <SeQuant/core/index.hpp>
   #include <SeQuant/core/math.hpp>
   #include <SeQuant/core/utility/expr.hpp>
   #include <SeQuant/core/utility/indices.hpp>
   #include <SeQuant/core/utility/macros.hpp>
   #include <SeQuant/core/utility/permutation.hpp>
   
   #include <Eigen/Eigenvalues>
   
   #include <libperm/Permutation.hpp>
   #include <libperm/Rank.hpp>
   #include <libperm/Utils.hpp>
   
   #include <algorithm>
   
   namespace sequant {
   
   template <typename T>
   struct compare_first_less {
     bool operator()(const T& lhs, const T& rhs) const {
       return lhs.first < rhs.first;
     }
   };
   
   using IndexPair = std::pair<Index, Index>;
   using ParticlePairings = container::svector<IndexPair>;
   
   ResultExpr biorthogonal_transform_copy(const ResultExpr& expr,
                                          double threshold) {
     container::svector<ResultExpr> wrapper = {expr.clone()};
   
     biorthogonal_transform(wrapper, threshold);
   
     return wrapper.front();
   }
   
   container::svector<ResultExpr> biorthogonal_transform_copy(
       const container::svector<ResultExpr>& exprs, double threshold) {
     container::svector<ResultExpr> copy;
     copy.reserve(exprs.size());
   
     std::transform(exprs.begin(), exprs.end(), std::back_inserter(copy),
                    [](const ResultExpr& expr) { return expr.clone(); });
   
     biorthogonal_transform(copy, threshold);
   
     return copy;
   }
   
   void biorthogonal_transform(ResultExpr& expr, double threshold) {
     // TODO: avoid copy
     expr = biorthogonal_transform_copy(expr, threshold);
   }
   
   Eigen::MatrixXd permutational_overlap_matrix(std::size_t n_particles) {
     const auto n = boost::numeric_cast<Eigen::Index>(factorial(n_particles));
   
     // The matrix only contains integer entries but all operations we want to do
     // with the matrix will (in general) require non-integer scalars which in
     // Eigen only works if you start from a non-integer matrix.
     Eigen::MatrixXd M(n, n);
     M.setZero();
   
     // TODO: Can we fill the entire matrix only by knowing the entries of one
     // row/column? For n_particles < 4, every consecutive col/row is only rotated
     // by one compared to the one before
     for (std::size_t row = 0; row < n; ++row) {
       perm::Permutation ref = perm::unrank(row, n_particles);
       ref->invert();
   
       // The identity permutation always has as many disjoint cycles as the number
       // of elements it acts on
       M(row, row) = std::pow(-2, n_particles);
   
       for (std::size_t col = row + 1; col < n; ++col) {
         // Get permutation that transforms the permutation of rank1 into the one
         // of current rank i
         perm::Permutation current = perm::unrank(col, n_particles);
         current->postMultiply(ref);
   
         auto cycles = current->toDisjointCycles(n_particles);
         std::size_t n_cycles = std::distance(cycles.begin(), cycles.end());
   
         auto entry = std::pow(-2, n_cycles);
   
         M(row, col) = entry;
         M(col, row) = entry;
       }
     }
   
     if (n_particles % 2 != 0) {
       M *= -1;
     }
   
     SEQUANT_ASSERT(M.isApprox(M.transpose()));
   
     return M;
   }
   
   Eigen::MatrixXd compute_biorth_coeffs(std::size_t n_particles,
                                         double threshold) {
     auto perm_ovlp_mat = permutational_overlap_matrix(n_particles);
     SEQUANT_ASSERT(perm_ovlp_mat.rows() == perm_ovlp_mat.cols());
     SEQUANT_ASSERT(perm_ovlp_mat.isApprox(perm_ovlp_mat.transpose()));
   
     // Find Pseudo Inverse
     auto decomp =
         Eigen::CompleteOrthogonalDecomposition<decltype(perm_ovlp_mat)>();
     decomp.setThreshold(threshold);
     decomp.compute(perm_ovlp_mat);
   
     Eigen::MatrixXd pinv = decomp.pseudoInverse();
     // The pseudo inverse of a symmetric matrix should also be symmetric
     SEQUANT_ASSERT(pinv.isApprox(pinv.transpose()));
   
     // We need to normalize to the amount of non-zero eigenvalues via
     // normalization = #eigenvalues / #non-zero eigenvalues
     // Since perm_ovlp_mat is symmetric, it is diagonalizable and for every
     // diagonalizable matrix, its rank equals the amount of non-zero eigenvalues.
     double normalization =
         static_cast<double>(perm_ovlp_mat.rows()) / decomp.rank();
   
     pinv *= normalization;
   
     return pinv;
   }
   
   void sort_pairings(ParticlePairings& pairing) {
     std::stable_sort(pairing.begin(), pairing.end(),
                      compare_first_less<IndexPair>{});
   }
   
   std::size_t rank_transformation_perms(const ParticlePairings& reference,
                                         const ParticlePairings& current) {
     SEQUANT_ASSERT(reference.size() == current.size());
     SEQUANT_ASSERT(std::is_sorted(reference.begin(), reference.end(),
                                   compare_first_less<IndexPair>{}));
     SEQUANT_ASSERT(std::is_sorted(current.begin(), current.end(),
                                   compare_first_less<IndexPair>{}));
   
     perm::Permutation perm = perm::computeTransformationPermutation(
         reference, current, [](const IndexPair& lhs, const IndexPair& rhs) {
           return lhs.second < rhs.second;
         });
   
     return perm::rank(perm, reference.size());
   }
   
   ExprPtr create_expr_for(const ParticlePairings& ref_pairing,
                           const perm::Permutation& perm,
                           const container::svector<ParticlePairings>& pairings,
                           const container::svector<ExprPtr>& base_exprs) {
     // Note: perm only applies to the p->second for every pair p in ref_pairing
   
     // assert that all pairings are sorted w.r.t. first
     SEQUANT_ASSERT(std::all_of(
         pairings.begin(), pairings.end(), [](const ParticlePairings& pairing) {
           return std::is_sorted(pairing.begin(), pairing.end(),
                                 compare_first_less<IndexPair>{});
         }));
     SEQUANT_ASSERT(std::is_sorted(ref_pairing.begin(), ref_pairing.end(),
                                   compare_first_less<IndexPair>{}));
   
     container::set<std::pair<IndexSpace, IndexSpace>> ref_space_pairing;
     ref_space_pairing.reserve(ref_pairing.size());
     for (std::size_t i = 0; i < ref_pairing.size(); ++i) {
       ref_space_pairing.emplace(ref_pairing[i].first.space(),
                                 ref_pairing[perm->image(i)].second.space());
     }
   
     // Look for a ParticlePairings object that pairs indices belonging to index
     // spaces compatible with ref_space_pairing
     auto it = std::find_if(
         pairings.begin(), pairings.end(), [&](const ParticlePairings& p) {
           SEQUANT_ASSERT(p.size() == ref_pairing.size());
   
           for (const IndexPair& pair : p) {
             if (ref_space_pairing.find(
                     std::make_pair(pair.first.space(), pair.second.space())) ==
                 ref_space_pairing.end()) {
               return false;
             }
           }
   
           return true;
         });
   
     if (it == pairings.end()) {
       throw std::runtime_error(
           "Missing explicit expression for a required index pairing in "
           "biorthogonalization");
     }
   
     auto idx = std::distance(pairings.begin(), it);
     const ParticlePairings& base = *it;
   
     SEQUANT_ASSERT(base.size() == ref_pairing.size());
   
     container::map<Index, Index> replacements;
     for (std::size_t i = 0; i < base.size(); ++i) {
       std::size_t ref_idx = perm->image(i);
   
       // Remember that all index pairings are sorted w.r.t. first and hence we are
       // only looking for permutations in second
       SEQUANT_ASSERT(base[i].first == ref_pairing[i].first);
       const bool differs_in_second =
           base[i].second != ref_pairing[ref_idx].second;
   
       if (!differs_in_second) {
         // This particle pairing is identical
         continue;
       }
   
       SEQUANT_ASSERT(differs_in_second);
   
       // Note: we may only permute indices belonging to the same space
       // (otherwise, we would produce non-sensical expressions)
       if (base[i].second.space() == ref_pairing[ref_idx].second.space()) {
         // base and ref_pairing differ in the second index of the current
         // pairing and their index space matches -> can just permute them
         replacements.emplace(base[i].second, ref_pairing[ref_idx].second);
       } else {
         // Index spaces of the differing index (second) in the pairings are
         // different as well. Since the tensors are assumed to be
         // particle-symmetric, we can instead permute the first indices in the
         // pairings, which are of the same space (that's guaranteed by the way we
         // chose base).
         SEQUANT_ASSERT(base[i].first.space() ==
                        ref_pairing[ref_idx].first.space());
         replacements.emplace(base[i].first, ref_pairing[ref_idx].first);
       }
     }
   
     ExprPtr expr = base_exprs.at(idx)->clone();
   
     if (!replacements.empty()) {
   #ifdef SEQUANT_ASSERT_ENABLED
       for (const auto& [first, second] : replacements) {
         SEQUANT_ASSERT(first.space() == second.space());
       }
   #endif
       expr = transform_expr(expr, replacements);
     }
   
     return expr;
   }
   
   void biorthogonal_transform(container::svector<ResultExpr>& result_exprs,
                               double threshold) {
     if (result_exprs.empty()) {
       return;
     }
   
     // We expect all ResultExpr objects to be equal except for the permutation of
     // indices
     // Also, we are assuming that all given ResultExpr objects are
     // particle-symmetric
     SEQUANT_ASSERT(std::all_of(
         result_exprs.begin(), result_exprs.end(), [&](const ResultExpr& expr) {
           return expr.has_label() == result_exprs.front().has_label() &&
                  (!expr.has_label() ||
                   expr.label() == result_exprs.front().label());
         }));
     SEQUANT_ASSERT(std::all_of(
         result_exprs.begin(), result_exprs.end(), [&](const ResultExpr& expr) {
           return expr.symmetry() == result_exprs.front().symmetry();
         }));
     SEQUANT_ASSERT(std::all_of(
         result_exprs.begin(), result_exprs.end(), [&](const ResultExpr& expr) {
           return expr.braket_symmetry() == result_exprs.front().braket_symmetry();
         }));
     SEQUANT_ASSERT(std::all_of(
         result_exprs.begin(), result_exprs.end(), [&](const ResultExpr& expr) {
           return expr.column_symmetry() == result_exprs.front().column_symmetry();
         }));
     SEQUANT_ASSERT(std::all_of(
         result_exprs.begin(), result_exprs.end(), [&](const ResultExpr& expr) {
           return expr.bra().size() == result_exprs.front().bra().size() &&
                  std::is_permutation(expr.bra().begin(), expr.bra().end(),
                                      result_exprs.front().bra().begin());
         }));
     SEQUANT_ASSERT(std::all_of(
         result_exprs.begin(), result_exprs.end(), [&](const ResultExpr& expr) {
           return expr.ket().size() == result_exprs.front().ket().size() &&
                  std::is_permutation(expr.ket().begin(), expr.ket().end(),
                                      result_exprs.front().ket().begin());
         }));
     SEQUANT_ASSERT(std::all_of(
         result_exprs.begin(), result_exprs.end(), [&](const ResultExpr& expr) {
           return expr.aux().size() == result_exprs.front().aux().size() &&
                  std::is_permutation(expr.aux().begin(), expr.aux().end(),
                                      result_exprs.front().aux().begin());
         }));
     SEQUANT_ASSERT(std::all_of(
         result_exprs.begin(), result_exprs.end(), [](const ResultExpr& res) {
           return res.column_symmetry() == ColumnSymmetry::Symm;
         }));
   
     // Furthermore, we expect that there is no symmetrization operator present in
     // the expressions as that would imply transforming also the symmetrization
     // operator, which is incorrect. This is because the idea during
     // biorthogonalization is that we project onto e.g.
     // \tilde{E}^{IJ}_{AB} = c_1 E^{IJ}_{AB} + c_2 E^{JI}_{AB}
     // instead of E^{IJ}_{AB} directly. In either case though, the result looks
     // like R^{IJ}_{AB} and the index pairing of the result is what determines
     // the required symmetrization. Hence, the symmetrization operator must not
     // be changed when transforming from one representation into the other.
     assert(std::all_of(
         result_exprs.begin(), result_exprs.end(), [](const ResultExpr& res) {
           bool found = false;
           res.expression()->visit(
               [&](const ExprPtr& expr) {
                 if (expr->is<Tensor>() && (expr->as<Tensor>().label() == L"S" ||
                                            expr->as<Tensor>().label() == L"A")) {
                   found = true;
                 };
               },
               true);
           return !found;
         }));
   
     auto externals = result_exprs |
                      ranges::views::transform([](const ResultExpr& expr) {
                        return expr.index_particle_grouping<IndexPair>();
                      }) |
                      ranges::to<container::svector<ParticlePairings>>();
     ranges::for_each(externals, sort_pairings);
   
     auto ranks = externals | ranges::views::transform([&](const auto& p) {
                    return rank_transformation_perms(externals.front(), p);
                  }) |
                  ranges::to<container::svector<std::size_t>>();
   
     const std::size_t n_particles = externals.front().size();
   
     Eigen::MatrixXd coefficients = compute_biorth_coeffs(n_particles, threshold);
   
     auto num_perms = factorial(n_particles);
     SEQUANT_ASSERT(num_perms == coefficients.rows());
     SEQUANT_ASSERT(num_perms == coefficients.cols());
   
     auto original_exprs = result_exprs |
                           ranges::views::transform([](const ResultExpr& res) {
                             return res.expression();
                           }) |
                           ranges::to<container::svector<ExprPtr>>();
   
     for (std::size_t i = 0; i < result_exprs.size(); ++i) {
       result_exprs.at(i).expression() = ex<Constant>(0);
       perm::Permutation reference = perm::unrank(ranks.at(i), n_particles);
       reference->invert();
   
       for (std::size_t rank = 0; rank < num_perms; ++rank) {
         perm::Permutation perm = perm::unrank(rank, n_particles);
         perm->postMultiply(reference);
   
         result_exprs.at(i).expression() +=
             ex<Constant>(
                 to_rational(coefficients(ranks.at(i), rank), threshold)) *
             create_expr_for(externals.at(i), perm, externals, original_exprs);
       }
   
       simplify(result_exprs.at(i).expression());
     }
   }
   
   ExprPtr biorthogonal_transform(
       const sequant::ExprPtr& expr,
       const container::svector<container::svector<sequant::Index>>&
           ext_index_groups,
       const double threshold) {
     ResultExpr res(
         bra(ext_index_groups | ranges::views::transform([](const auto& pair) {
               return get_ket_idx(pair);
             }) |
             ranges::to<container::svector<Index>>()),
         ket(ext_index_groups | ranges::views::transform([](const auto& pair) {
               return get_bra_idx(pair);
             }) |
             ranges::to<container::svector<Index>>()),
         aux(IndexList{}), Symmetry::Nonsymm, BraKetSymmetry::Nonsymm,
         ColumnSymmetry::Symm, {}, expr);
   
     biorthogonal_transform(res, threshold);
   
     return res.expression();
   }
   
   }  // namespace sequant