.. _program_listing_file_SeQuant_domain_mbpt_antisymmetrizer.hpp:

Program Listing for File antisymmetrizer.hpp
============================================

|exhale_lsh| :ref:`Return to documentation for file <file_SeQuant_domain_mbpt_antisymmetrizer.hpp>` (``SeQuant/domain/mbpt/antisymmetrizer.hpp``)

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

.. code-block:: cpp

   //
   // created by Conner Masteran 06/1/2021
   //
   
   #ifndef SEQUANT_DOMAIN_MBPT_ANTISYMMETRIZER_HPP
   #define SEQUANT_DOMAIN_MBPT_ANTISYMMETRIZER_HPP
   
   #include <SeQuant/core/expr.hpp>
   #include <SeQuant/core/math.hpp>
   #include <SeQuant/core/op.hpp>
   #include <SeQuant/core/tensor.hpp>
   #include <algorithm>
   #include <iostream>
   #include <vector>
   
   namespace sequant {
   
   class antisymm_element {
     using IndexGroup = std::pair<size_t, size_t>;
   
    private:
     std::vector<IndexGroup>
         index_group;  // where each tensor begins and ends. assumes particle
                       // conserving ops. needed to keep track of canonical
                       // ordering
     std::vector<std::pair<int, std::vector<Index>>>
         unique_bras_list;  // list of unique bra orderings with associated integer
                            // for the sign
     std::vector<std::pair<int, std::vector<Index>>>
         unique_kets_list;  // list of unique ket orderings with associated integer
                            // for the sign
     ExprPtr current_product;  // used to keep track of the original expression
                               // recieved by the constructor
   
     // generates all possible permutations while observing the canonical ordering
     // of each tensor/operator function kept general to work with other data
     // types, algorithm does not require sequant::Index objects
     // @param a list of T objects/indices in the canonical ordering
     // @return a list of all permutations each with an associated sign
     template <typename T>
     std::vector<std::pair<int, std::vector<T>>> gen_antisymm_unique(
         std::vector<T> ordered_indices) {
       std::vector<std::pair<int, std::vector<T>>> result;
       // next_permutation needs an ordered list. works for integers quite well.
       // Easiest to map T to an integer corresponding to its position in a vector
       // and then go back.
       std::vector<int> ordered_numbers;
       for (int i = 0; i < ordered_indices.size(); i++) {
         ordered_numbers.push_back(i);
       }
   
       // return {next_exists, nswaps for this permutation} N.B.
       // nswaps_relative_to_input != nswaps_relative_to_оriginal
       auto swapcounting_tracking_next_permutation =
           [](auto first, auto last) -> std::pair<bool, int> {
         int nswaps = 0;
         using BidirIt = decltype(first);
         if (first == last) return std::make_pair(false, 0);
         BidirIt i = last;
         if (first == --i) return std::make_pair(false, 0);
         while (true) {
           BidirIt i1, i2;
           i1 = i;
           if (*--i < *i1) {
             i2 = last;
             while (!(*i < *--i2))
               ;
             std::iter_swap(i, i2);
             ++nswaps;
             std::reverse(i1, last);
             nswaps += (std::distance(i1, last)) /
                       2;  // logic from
                           // https://en.cppreference.com/w/cpp/algorithm/reverse
             return std::make_pair(true, nswaps);
           }
           if (i == first) {
             std::reverse(first, last);
             nswaps += (std::distance(first, last)) / 2;
             return std::make_pair(false, nswaps);
           }
         }
       };
   
       std::vector<int> numbers = ordered_numbers;
       result.push_back({1, ordered_indices});
   
       int total_swaps = 0;  // even # swaps produces positive and odd # of swaps
                             // produce negative
       [[maybe_unused]] int counter = 0;
   
       bool do_next_perm = true;
   
       while (do_next_perm) {
         auto do_swaps =
             swapcounting_tracking_next_permutation(begin(numbers), end(numbers));
         do_next_perm = do_swaps.first;
         total_swaps += do_swaps.second;
         if (!do_next_perm) {
           break;
         }
         auto is_canonical_sign =
             [this, &total_swaps](const auto& indices) -> std::pair<bool, bool> {
           for (auto& group :
                this->index_group) {  // There is only one sorted possibility in a
                                      // set (tensor) considering that no index
                                      // label should be the same.
             if (!is_sorted(indices.begin() + group.first,
                            indices.begin() + group.second))
               return {
                   false,
                   total_swaps % 2 ==
                       0};  // total swaps divisible by 2 means true = positive.
           }
           return {true, total_swaps % 2 == 0};
         };
         // sieve out non-canonical terms
         auto [is_canonical, sign] = is_canonical_sign(numbers);
         if (is_canonical) {
           std::vector<T> return_vec(ordered_indices.size());
           for (int i = 0; i < numbers.size(); i++) {
             return_vec[i] = ordered_indices[numbers[i]];
           }
           if (sign) {
             result.push_back({1, return_vec});
           } else {
             result.push_back({-1, return_vec});
           }
         }
         counter += 1;
       }
       return result;
     }
   
    public:
     // takes a sequant::ExprPtr and generates all antisymmetric unique
     // permutations of that expression. requires that ex_ is a product expression
     // at this point
     // @param ex_ as product
     // populates a result ExprPtr that the user can grab. result is in general a
     // Sum.
     antisymm_element(ExprPtr ex_) {
       current_product = ex_;
       assert(ex_->is<Product>());
       Constant::scalar_type starting_constant = 1;
       unsigned long begining_index = 0;
       for (auto it = begin(*ex_); it != end(*ex_); it++) {
         if (it->get()->is<Tensor>()) {
           auto factor = it->get()->as<Tensor>();
           index_group.push_back(
               {begining_index, begining_index + factor.bra_rank()});
           begining_index += factor.bra_rank();
           assert(factor.bra_rank() == factor.ket_rank());
           for (int i = 0; i < factor.rank(); i++) {
             sorted_bra_indices.push_back(factor.bra()[i]);
             sorted_ket_indices.push_back(factor.ket()[i]);
           }
         } else if (it->get()->is<Constant>()) {
           starting_constant *= it->get()->as<Constant>().value();
         }
   
         else if (it->get()->is<FNOperator>()) {
           auto factor = it->get()->as<FNOperator>();
           index_group.push_back(
               {begining_index, begining_index + factor.nannihilators()});
           begining_index += factor.nannihilators();
           assert(factor.ncreators() == factor.nannihilators());
           for (int i = 0; i < factor.rank(); i++) {
             sorted_bra_indices.push_back(factor.annihilators()[i].index());
             sorted_ket_indices.push_back(factor.creators()[i].index());
           }
         } else {
           throw " unknown type of product, factor is not tensor, constant, or NormalOperator with FermiDirac statistics";
         }
       }
   
       unique_bras_list = gen_antisymm_unique(sorted_bra_indices);
       unique_kets_list = gen_antisymm_unique(sorted_ket_indices);
   
       auto new_sum = ex<Constant>(0);
   
       auto summand_exists =
           [&new_sum](ExprPtr ex) {  // check whether a summand has already been
                                     // generated to screen out same terms.
             bool value = false;
             for (auto summand = new_sum->begin_subexpr();
                  summand != new_sum->end_subexpr(); summand++) {
               value = ex.get()->as<Product>() ==
                       summand->get()
                           ->as<Product>();  // ensure that this equality is
                                             // mathematical and not hash based.
               if (value == true) {
                 return value;
               }
             }
             return value;
           };
   
       for (int i = 0; i < unique_bras_list.size(); i++) {
         for (int j = 0; j < unique_kets_list.size(); j++) {  // product level
   
           auto new_product =
               ex<Constant>(unique_kets_list[i].first * unique_bras_list[j].first *
                            starting_constant);
   
           int index_label_pos = 0;
           for (auto it = begin(*ex_); it != end(*ex_); it++) {  // factor level
             if (it->get()->is<Constant>()) {
             }  // constant already captured in the first loop.
             else if (it->get()->is<Tensor>()) {
               auto old_tensor = it->get()->as<Tensor>();
               auto label = old_tensor.label();
               std::vector<Index> new_bras;
               std::vector<Index> new_kets;
               for (auto k = 0; k < old_tensor.rank(); k++) {  // index level
                 new_bras.push_back(unique_bras_list[i].second[index_label_pos]);
                 new_kets.push_back(unique_kets_list[j].second[index_label_pos]);
                 index_label_pos++;
               }
               auto new_tensor = ex<Tensor>(label, bra(std::move(new_bras)),
                                            ket(std::move(new_kets)));
               new_product = new_tensor * new_product;
               new_product->canonicalize();
             }
   
             else if (it->get()->is<FNOperator>()) {
               auto old_Nop = it->get()->as<FNOperator>();
               std::vector<Index> new_anni;
               std::vector<Index> new_crea;
               for (auto k = 0; k < old_Nop.rank(); k++) {
                 new_anni.push_back(unique_bras_list[i].second[index_label_pos]);
                 new_crea.push_back(unique_kets_list[j].second[index_label_pos]);
                 index_label_pos++;
               }
               auto new_Nop = ex<FNOperator>(cre(new_crea), ann(new_anni));
               new_product = new_product * new_Nop;
               // std::wcout << "product:  " << to_latex(new_product) << std::endl;
               new_product->canonicalize();
             }
   
             else {
               throw " unknown type of product, factor is not tensor, constant, or NormalOperator with FermiDirac statistics";
             }
           }
           if (!summand_exists(
                   new_product)) {  // since products are canonicalized, repeats
                                    // can be found.
             new_sum = new_sum + new_product;
           }
         }
       }
   
       result = new_sum;
     }
   
     std::vector<Index> sorted_bra_indices;  // The original order of the upper
                                             // indices on a given term
     std::vector<Index> sorted_ket_indices;  // the original order of the lower
                                             // indices on a given term
     ExprPtr result;
   };
   
   class antisymmetrize {
    public:
     ExprPtr result = ex<Constant>(0);
   
     antisymmetrize(ExprPtr s) {
       if (s->is<Sum>()) {
         for (auto&& product :
              s->as<Sum>().summands()) {  // for each element in the sum
           if (product->is<Product>()) {
             antisymm_element ele(
                 product);  // calculate the sum of all the valid permutations for
                            // each term. each object here should only live until
                            // this loop ends.
             result = result + ele.result;  // append that to the final list;
           } else {
             result = result + product;
           }
         }
       } else if (s->is<Product>()) {
         antisymm_element answer(ex<Product>(s->as<Product>()));
         result = answer.result;
       } else {
         result = s;
       }
     }
   };
   
   namespace antisymm {
   // function which counts number of closed loops from the starting order of upper
   // and lower indices and the contracted indices. since the ordering of the new
   // contracted indices is arbitrary, the algorithm searched for the upper index
   // which would connect to the lower index checks if contracted lower index
   // closes the loop, if not, continue searching until the corresponding upper
   // index is not present, or the loop closes. keep track of which indices are
   // used so that loops are not double counted
   //@param1 initial order of upper indices before contraction, @param 2 initial
   // order of lower indices before contraction, @param 3 set of contracted upper
   // indices, @param 4 set of lower contracted indices.
   // @out the number of loops.
   // TODO Test this function extensively and add more asserts
   int num_closed_loops(std::vector<Index> init_upper,
                        std::vector<Index> init_lower,
                        std::vector<Index> new_upper,
                        std::vector<Index> new_lower) {
     int result = 0;
     auto equal_indices = [](Index i1, Index i2) {
       return i1.label() == i2.label();
     };
     auto where_same_ele = [&](Index i1, std::vector<Index> ref_list) {
       int hit_counter = 0;
       int where;
       bool in_list = false;
       for (int i = 0; i < ref_list.size(); i++) {
         if (equal_indices(i1, ref_list[i])) {
           hit_counter += 1;
           where = i;
           in_list = true;
         }
       }
       assert(hit_counter < 2);
       std::pair<int, bool> result(where, in_list);
       return result;
     };
   
     std::vector<Index> in_loop;  // lower indices already in a loop.
     for (int i = 0; i < new_upper.size(); i++) {
       assert(new_upper.size() == new_lower.size());
       if (!where_same_ele(new_lower[i], in_loop).second) {
         auto starting_point = where_same_ele(new_upper[i], init_upper).first;
         auto chain_end = init_lower[starting_point];
   
         bool closed = equal_indices(chain_end, new_lower[i]);
         auto lower_index = new_lower[i];
         while (!closed) {
           auto where_exist = where_same_ele(
               init_upper[where_same_ele(lower_index, init_lower).first],
               new_upper);  // if the  initial upper index is the same particle as
                            // the new lower index in question, find it in the new
                            // upper list.
           if (!where_exist
                    .second) {  // if the upper particle index originally connected
                                // to the lower index in question is not part of a
                                // contraction, there is no loop with these.
             break;
           } else {
             lower_index = new_lower[where_exist.first];  // the lower index below
                                                          // the found upper index
             in_loop.push_back(
                 lower_index);  // this lower index is part of one loop so it
                                // cannot be part of another.
             closed = equal_indices(chain_end,
                                    lower_index);  // is the lower index the same
                                                   // as the end of the chain?
           }
         }
         if (closed) {
           result += 1;
         }
         in_loop.push_back(
             new_lower[i]);  // put the starting lower index in the list
       }
     }
     assert(in_loop.size() <= init_lower.size());
     assert(result <= init_lower.size());
     return result;
   }
   
   // for the mnemonic rules of spin summing to work, each individual
   // tensor/FNOperator needs to maximally resemble the original tensor, so indices
   // may need to be swapped (and the sign changed).
   //  this does not matter for spin-orbital expressions, but the rules become
   //  incredibly simple if things stay maximally the same.
   //  @param1 original upper indices
   //  @param2 original lower indices
   //  @param3 starting expression
   ExprPtr max_similarity(const std::vector<Index>& original_upper,
                          const std::vector<Index>& original_lower,
                          ExprPtr expression) {
     // index pairing is originally understood as a position in the original
     // vectors, but for this case, a map may do better.
     std::map<std::wstring, std::wstring> original_map;
     for (int i = 0; i < original_upper.size(); i++) {
       original_map.emplace(original_upper[i].to_latex(),
                            original_lower[i].to_latex());
     }
     for (auto&& product : expression->as<Sum>().summands()) {
       for (auto&& factor : product->as<Product>().factors()) {
         int og_pairs = 0;
         int new_pairs = 0;
         if (factor->is<Tensor>()) {
           std::vector<Index> current_upper;
           std::vector<Index> current_lower;
           if (factor->as<Tensor>().rank() == 2) {
             for (int i = 0; i < 2; i++) {
               assert(
                   original_map.find(factor->as<Tensor>().ket()[i].to_latex()) !=
                   original_map.end());
               if (original_map.find(factor->as<Tensor>().ket()[i].to_latex())
                       ->second == factor->as<Tensor>().bra()[i].to_latex()) {
                 og_pairs += 1;
               }
               current_upper.push_back(factor->as<Tensor>().ket()[i]);
               current_lower.push_back(factor->as<Tensor>().bra()[i]);
             }
             std::iter_swap(current_lower.begin(), current_lower.begin() + 1);
             for (int i = 0; i < 2; i++) {
               assert(original_map.find(current_upper[i].to_latex()) !=
                      original_map.end());
               if (original_map.find(current_upper[i].to_latex())->second ==
                   current_lower[i].to_latex()) {
                 new_pairs += 1;
               }
             }
           }
           if (new_pairs > og_pairs) {
             factor = ex<Constant>(-1) * ex<Tensor>(factor->as<Tensor>().label(),
                                                    bra(std::move(current_lower)),
                                                    ket(std::move(current_upper)));
           }
         } else if (factor->is<FNOperator>()) {
           std::vector<Index> current_upper;
           std::vector<Index> current_lower;
           if (factor->as<FNOperator>().rank() == 2) {
             for (int i = 0; i < 2; i++) {
               assert(original_map.find(factor->as<FNOperator>()
                                            .creators()[i]
                                            .index()
                                            .to_latex()) != original_map.end());
               if (original_map
                       .find(factor->as<FNOperator>()
                                 .creators()[i]
                                 .index()
                                 .to_latex())
                       ->second ==
                   factor->as<FNOperator>().annihilators()[i].index().to_latex()) {
                 og_pairs += 1;
               }
               current_upper.push_back(
                   factor->as<FNOperator>().creators()[i].index());
               current_lower.push_back(
                   factor->as<FNOperator>().annihilators()[i].index());
             }
             std::iter_swap(current_lower.begin(), current_lower.begin() + 1);
             for (int i = 0; i < 2; i++) {
               assert(original_map.find(current_upper[i].to_latex()) !=
                      original_map.end());
               if (original_map.find(current_upper[i].to_latex())->second ==
                   current_lower[i].to_latex()) {
                 new_pairs += 1;
               }
             }
           }
           if (new_pairs > og_pairs) {
             factor =
                 ex<Constant>(-1) * ex<FNOperator>(cre(std::move(current_upper)),
                                                   ann(std::move(current_lower)));
           }
         }
       }
     }
     simplify(expression);
     return expression;
   }
   
   // not a general spin-summing procedure, implementation for a known singlet
   // state for the prefactor rules to apply.
   // @param1 an antisymm_element object strictly so the original ordering of the
   // indices is known
   // @param2 bool singlet_state? the looping rules and contraction prefactors are
   // a direct result of the singlet state approximation to densities.
   // TODO: use a generalized spin summation for non-singlet states
   ExprPtr spin_sum(std::vector<Index> original_upper,
                    std::vector<Index> original_lower, ExprPtr expression,
                    bool singlet_state = true) {
     if (singlet_state) {
       // std::wcout << "before spin sum! :" << std::endl <<
       // to_latex_align(expression) << std::endl;
       auto init_upper = original_upper;
       auto init_lower = original_lower;
       max_similarity(original_upper, original_lower, expression);
       // may need to add separate loop if the result is a single product or
       // Operator/Tensor
       auto return_val = ex<Constant>(0);
       for (auto&& product : expression->as<Sum>().summands()) {
         auto prefactor = ex<Constant>(
             rational{1, 8});  // each term in the 3 body decomp, should have
                               // (1 /2^n) prefactor where n is rank of product
                               // so product rank always 3 for 3 body decomp
         std::vector<Index> new_upper;
         std::vector<Index> new_lower;
         for (auto&& factor : product->as<Product>().factors()) {
           if (factor->is<Tensor>() && factor->as<Tensor>().label() == L"γ") {
             // prefactor = ex<Constant>(-0.5) *
             // ex<Constant>(factor->as<Tensor>().rank()) * prefactor;
             for (int i = 0; i < factor->as<Tensor>().rank(); i++) {
               new_upper.push_back(factor->as<Tensor>().ket()[i]);
               new_lower.push_back(factor->as<Tensor>().bra()[i]);
             }
             factor = ex<Tensor>(L"Γ", factor->as<Tensor>().bra(),
                                 factor->as<Tensor>().ket());
           } else if (factor->is<FNOperator>()) {
             // prefactor = ex<Constant>(-0.5) *
             // ex<Constant>(factor->as<Tensor>().rank()) * prefactor;
             for (int i = 0; i < factor->as<FNOperator>().rank(); i++) {
               new_upper.push_back(factor->as<FNOperator>().creators()[i].index());
               new_lower.push_back(
                   factor->as<FNOperator>().annihilators()[i].index());
             }
           }
         }
         // std::wcout << to_latex_align(product) << std::endl;
         int nloops =
             num_closed_loops(init_upper, init_lower, new_upper, new_lower);
         // std::wcout  << "nloops: " << nloops << std::endl;
         if (nloops == 0) {
         } else {
           prefactor = ex<Constant>(pow2(nloops)) * prefactor;
         }
         return_val = (product * prefactor) + return_val;
       }
       return_val->canonicalize();
       return return_val;
     } else {
       throw " non-singlet states not yet supported";
     }
   }
   }  // namespace antisymm
   }  // namespace sequant
   
   #endif  // SEQUANT_DOMAIN_MBPT_ANTISYMMETRIZER_HPP