Merge pull request #258 from ValeevGroup/evaleev/feature/canonicalize…

…-slots-phase

`TN::canonicalize_slots` produces phase associated with reordering antisymmetric bra/ket slots bundles

CMakeLists.txt

@@ -262,6 +262,7 @@ set(SeQuant_src
         SeQuant/core/utility/string.hpp
         SeQuant/core/utility/string.cpp
         SeQuant/core/utility/tuple.hpp
+        SeQuant/core/utility/swap.hpp
         SeQuant/core/wick.hpp
         SeQuant/core/wick.impl.hpp
         SeQuant/core/wolfram.hpp
@@ -415,7 +416,7 @@ if (SEQUANT_IWYU)
 endif()
 
 if (SEQUANT_BENCHMARKS)
-	add_subdirectory(benchmarks)
+    add_subdirectory(benchmarks)
 endif()
 
 ### unit tests

SeQuant/core/algorithm.hpp

@@ -22,8 +22,9 @@ using suitable_call_operator =
     decltype(std::declval<Callable>()(std::declval<Args>()...));
 
 /// @brief bubble sort that uses swap exclusively
-template <typename ForwardIter, typename Sentinel, typename Compare>
-void bubble_sort(ForwardIter begin, Sentinel end, Compare comp) {
+template <typename ForwardIter, typename Sentinel,
+          typename Compare = std::less<>>
+void bubble_sort(ForwardIter begin, Sentinel end, Compare comp = {}) {
   if (begin == end) return;  // no-op for empty range
   bool swapped;
   do {

SeQuant/core/eval_expr.cpp

@@ -265,9 +265,9 @@ size_t hash_imed(EvalExpr const& left, EvalExpr const& right,
   return h;
 }
 
-std::pair<container::svector<Index>,  // bra
-          container::svector<Index>   // ket
-          >
+[[maybe_unused]] std::pair<container::svector<Index>,  // bra
+                           container::svector<Index>   // ket
+                           >
 target_braket(Tensor const& t1, Tensor const& t2) noexcept {
   using ranges::contains;
   using ranges::views::concat;

SeQuant/core/index.hpp

@@ -6,14 +6,11 @@
 #define SEQUANT_INDEX_H
 
 #include <SeQuant/core/container.hpp>
-#include <SeQuant/core/hash.hpp>
-#include <iostream>
-// #include <SeQuant/core/space.hpp>
 #include <SeQuant/core/context.hpp>
+#include <SeQuant/core/hash.hpp>
 #include <SeQuant/core/tag.hpp>
 #include <SeQuant/core/utility/string.hpp>
-// Only needed due to a (likely) compiler bug in Apple Clang
-// #include <SeQuant/core/attr.hpp>
+#include <SeQuant/core/utility/swap.hpp>
 
 #include <algorithm>
 #include <atomic>
@@ -23,6 +20,7 @@
 #include <cwchar>
 #include <functional>
 #include <initializer_list>
+#include <iostream>
 #include <iterator>
 #include <map>
 #include <mutex>
@@ -868,28 +866,10 @@ void Index::canonicalize_proto_indices() {
     std::stable_sort(begin(proto_indices_), end(proto_indices_));
 }
 
-class IndexSwapper {
- public:
-  IndexSwapper() : even_num_of_swaps_(true) {}
-  static IndexSwapper &thread_instance() {
-    static thread_local IndexSwapper instance_{};
-    return instance_;
-  }
-
-  bool even_num_of_swaps() const { return even_num_of_swaps_; }
-  void reset() { even_num_of_swaps_ = true; }
-
- private:
-  std::atomic<bool> even_num_of_swaps_;
-  void toggle() { even_num_of_swaps_ = !even_num_of_swaps_; }
-
-  friend inline void swap(Index &, Index &);
-};
-
 /// swap operator helps tracking # of swaps
 inline void swap(Index &first, Index &second) {
   std::swap(first, second);
-  IndexSwapper::thread_instance().toggle();
+  detail::count_swap<Index>();
 }
 
 /// Generates temporary indices

SeQuant/core/tensor_canonicalizer.hpp

@@ -184,14 +184,13 @@ class DefaultTensorCanonicalizer : public TensorCanonicalizer {
         auto _bra = bra_range(t);
         auto _ket = ket_range(t);
         //      std::wcout << "canonicalizing " << to_latex(t);
-        IndexSwapper::thread_instance().reset();
+        reset_ts_swap_counter<Index>();
         // std::{stable_}sort does not necessarily use swap! so must implement
         // sort outselves .. thankfully ranks will be low so can stick with
         // bubble
         bubble_sort(begin(_bra), end(_bra), idxcmp);
         bubble_sort(begin(_ket), end(_ket), idxcmp);
-        if (is_antisymm)
-          even = IndexSwapper::thread_instance().even_num_of_swaps();
+        if (is_antisymm) even = ts_swap_counter_is_even<Index>();
         //      std::wcout << " is " << (even ? "even" : "odd") << " and
         //      produces " << to_latex(t) << std::endl;
       } break;

SeQuant/core/tensor_network.cpp

@@ -14,6 +14,7 @@
 #include <SeQuant/core/latex.hpp>
 #include <SeQuant/core/logger.hpp>
 #include <SeQuant/core/tag.hpp>
+#include <SeQuant/core/tensor.hpp>
 #include <SeQuant/core/tensor_network.hpp>
 #include <SeQuant/core/tensor_network/vertex_painter.hpp>
 #include <SeQuant/core/tensor_network_v2.hpp>
@@ -212,7 +213,7 @@ ExprPtr TensorNetwork::canonicalize(
       container::multimap<size_t, std::pair<size_t, size_t>> color2idx;
       // collect colors and anonymous indices sorted by colors
       size_t idx_ord = 0;
-      for (auto &&ttpair : edges_) {
+      for ([[maybe_unused]] auto &&ttpair : edges_) {
         if (is_anonymous_index_ord(idx_ord)) {
           auto color = vcolors[idx_ord];
           if (colors.find(color) == colors.end()) colors.insert(color);
@@ -511,7 +512,7 @@ TensorNetwork::make_bliss_graph(const named_indices_t *named_indices_ptr,
     assert(!bundle.empty());
     ++nv;  // each symmetric protoindex bundle is a vertex
     std::wstring spbundle_label = L"<";
-    std::size_t pi_count = 0;
+    [[maybe_unused]] std::size_t pi_count = 0;
     const auto end = bundle.end();
     auto it = bundle.begin();
     spbundle_label += it->full_label();
@@ -621,7 +622,7 @@ TensorNetwork::make_bliss_graph(const named_indices_t *named_indices_ptr,
   index_cnt = 0;
   ranges::for_each(edges_, [&](const Edge &edge) {
     assert(edge.size() > 0);
-    const auto edge_connected = edge.size() == 2;
+    [[maybe_unused]] const auto edge_connected = edge.size() == 2;
     for (int t = 0; t != edge.size(); ++t) {
       const auto &terminal = edge[t];
       const auto tensor_ord = terminal.tensor_ord;
@@ -872,7 +873,10 @@ TensorNetwork::SlotCanonicalizationMetadata TensorNetwork::canonicalize_slots(
   // distinct_named_indices = false
   auto [graph, vlabels, vcolors, vtypes] =
       make_bliss_graph(&named_indices, /* distinct_named_indices = */ false);
-  //    graph->write_dot(std::wcout, vlabels);
+  if (Logger::instance().canonicalize_dot) {
+    std::wcout << "Input graph for canonicalization:\n";
+    graph->write_dot(std::wcout, vlabels);
+  }
 
   // canonize the graph
   bliss::Stats stats;
@@ -890,6 +894,7 @@ TensorNetwork::SlotCanonicalizationMetadata TensorNetwork::canonicalize_slots(
       return pvector;
     };
 
+    std::wcout << "Canonicalized graph:\n";
     auto cvlabels = permute(vlabels, cl);
     metadata.graph->write_dot(std::wcout, cvlabels);
   }
@@ -961,6 +966,115 @@ TensorNetwork::SlotCanonicalizationMetadata TensorNetwork::canonicalize_slots(
 
   }  // named indices resort to canonical order
 
+  // compute the phase associated with *slot* canonicalization ...
+  // - choose an index-independent order of slots: loop over bra then ket then
+  // aux slots of each tensor ... record each Index in the order of its
+  // appearance when iterating over slots. This is the input ordinal of this
+  // Index. NB We can't just use Index::full_label() to look up Index in edges_
+  // since this would produce label-dependent ordinals
+  // - canonicalization reorders slots according to the topology-based order of
+  // indices.
+  // - for each antisymmetric bra/ket bundle determine its parity according to
+  // the input ordinals and canonical ordinals, the product is the phase due to
+  // slot canonicalization
+  {
+    metadata.phase = 1;
+
+    // iterate over tensor bra/ket/aux index slots in the canonical input order,
+    // for each antisymmetric bra/ket compute phase
+    container::map<Index, std::size_t, FullLabelCompare>
+        idx_inord;  // Index -> input ordinal; helps with computing the ordinals
+                    // during the traversal
+    for (auto &_t : tensors_) {
+      assert(std::dynamic_pointer_cast<Tensor>(_t));
+      auto t = std::static_pointer_cast<Tensor>(_t);
+
+      // returns an iterator to {Index,inord} pair
+      auto index_inord_it = [&](const Index &idx) {
+        auto it = idx_inord.find(idx.full_label());
+        if (it == idx_inord.end()) {
+          const auto inord = idx_inord.size();
+          bool inserted;
+          std::tie(it, inserted) = idx_inord.emplace(idx, inord);
+          assert(inserted);
+        }
+        return it;
+      };
+
+      // computes parity of a bra/ket bundle due to the reordering of its slots
+      // involved in the TN canonicalization
+      auto input_to_canonical_parity = [&](const auto &idx_rng) {
+        using ranges::size;
+        const auto sz = size(idx_rng);
+        if (sz < 2) {  // no phase for 1-index bundles, but still process the
+                       // indices to ensure ordinals are correct
+          using ranges::begin;
+          index_inord_it(*begin(idx_rng));
+          return 1;
+        }
+
+        auto parity = [&](auto &rng) {
+          reset_ts_swap_counter<std::size_t>();
+          using ranges::begin;
+          using ranges::end;
+          bubble_sort(begin(rng), end(rng));
+          return ts_swap_counter_is_even<std::size_t>() ? +1 : -1;
+        };
+
+        container::vector<SwapCountable<std::size_t>> ordinals_input;
+        container::vector<SwapCountable<std::size_t>> ordinals_canonical;
+        ordinals_input.reserve(sz);
+        ordinals_canonical.reserve(sz);
+        for (const auto &idx : idx_rng) {
+          // input ordinal
+          std::size_t inord;
+          std::tie(std::ignore, inord) = *index_inord_it(idx);
+          ordinals_input.emplace_back(inord);
+
+          // use canonical ordinal of the index vertex as the canonical index
+          // ordinal to find it need the ordinal of the corresponding vertex in
+          // the input graph
+          auto input_vertex_it = this->edges_.find(idx.full_label());
+          assert(input_vertex_it != this->edges_.end());
+          const auto inord_vertex = input_vertex_it - this->edges_.begin();
+          const auto canord_vertex = cl[inord_vertex];
+
+          // to which edge did this get mapped?
+          //          std::wcout << "vlabels[ord_original=" << inord_vertex
+          //                     << "]=" << vlabels[inord_vertex]
+          //                     << " -> ord_canonical=" << canord_vertex <<
+          //                     std::endl;
+          ordinals_canonical.emplace_back(canord_vertex);
+        }
+
+        // parity = parity(original range) * parity(canonical range)
+        const auto parity_original = parity(ordinals_input);
+        const auto parity_canonical = parity(ordinals_canonical);
+        return parity_original * parity_canonical;
+      };
+
+      // canonical order of slots: bra, then ket, then aux. We only care about
+      // indices in tensor slots here, so no need to worry about protoindex
+      // slots
+      if (t->symmetry() == Symmetry::antisymm) {
+        // bra first, then ket
+        metadata.phase *= input_to_canonical_parity(t->bra());
+        metadata.phase *= input_to_canonical_parity(t->ket());
+      } else {  // although we don't need to worry about phases, we still need
+                // to process all indices so that the input ordinals are correct
+        for (auto &&idx : t->bra()) {
+          index_inord_it(idx);
+        }
+        for (auto &&idx : t->ket()) {
+          index_inord_it(idx);
+        }
+      }
+      for (auto &&idx : t->aux()) {
+        index_inord_it(idx);
+      }
+    }
+  }
+
   return metadata;
 }
 

SeQuant/core/tensor_network.hpp

@@ -216,8 +216,7 @@ class TensorNetwork {
         } else {
           throw std::logic_error(
               "TensorNetwork::TensorNetwork: non-tensors in the given "
-              "expression "
-              "range");
+              "expression range");
         }
       }
       return;
@@ -284,6 +283,11 @@ class TensorNetwork {
     /// canonicalized colored graph, use graph->cmp to compare against another
     /// to detect equivalence
     std::shared_ptr<bliss::Graph> graph;
+
+    /// if tensor network contains tensors with antisymmetric bra/ket this
+    /// reports the phase change due to permutation of slots relative to their
+    /// input order
+    std::int8_t phase = +1;  // +1 or -1
   };
 
   /// Like canonicalize(), but only use graph-based canonicalization to
@@ -336,6 +340,15 @@ class TensorNetwork {
     bool operator()(const std::wstring_view &first, const Edge &second) const {
       return first < second.idx().full_label();
     }
+    bool operator()(const Index &first, const Index &second) const {
+      return first.full_label() < second.full_label();
+    }
+    bool operator()(const Index &first, const std::wstring_view &second) const {
+      return first.full_label() < second;
+    }
+    bool operator()(const std::wstring_view &first, const Index &second) const {
+      return first < second.full_label();
+    }
   };
   // Index -> Edge, sorted by full label
   using edges_t = container::set<Edge, FullLabelCompare>;

SeQuant/core/tensor_network_v2.hpp

@@ -264,6 +264,11 @@ class TensorNetworkV2 {
     /// canonicalized colored graph, use graph->cmp to compare against another
     /// to detect equivalence
     std::shared_ptr<bliss::Graph> graph;
+
+    /// if tensor network contains tensors with antisymmetric bra/ket this
+    /// reports the phase change due to permutation of slots relative to their
+    /// input order
+    std::int8_t phase = +1;  // +1 or -1
   };
 
   /// Like canonicalize(), but only use graph-based canonicalization to