Support non-covariant indices in evaluation expressions.

SeQuant/core/algorithm.hpp

@@ -138,34 +138,6 @@ bool next_permutation_parity(int& parity, BidirIt first, BidirIt last,
   }
 }
 
-///
-/// Given a size of a container @c n, returns a range of bitsets. The bitsets
-/// represent the sieve that can be used to construct the subsequences of the
-/// size-n container.
-///
-inline auto subset_indices(size_t n) {
-  using ranges::views::iota;
-  using ranges::views::transform;
-  return iota(size_t{1}, size_t(1 << n)) |
-         transform([n](auto i) { return boost::dynamic_bitset<>(n, i); });
-}
-
-///
-/// Given an iterable @c it and a bitset @c bs, select the elements in the
-/// iterable that correspond to an 'on' bit.
-///
-template <typename Iterable>
-auto subsequence(Iterable const& it, boost::dynamic_bitset<> const& bs) {
-  using ranges::views::filter;
-  using ranges::views::iota;
-  using ranges::views::transform;
-  using ranges::views::zip;
-
-  auto bits = iota(size_t{0}) | transform([&bs](auto i) { return bs.at(i); });
-  return zip(it, bits) | filter([](auto&& kv) { return std::get<1>(kv); }) |
-         transform([](auto&& kv) { return std::get<0>(kv); });
-}
-
 ///
 /// All elements in the vector belong to the integral range [-1,N)
 /// where N is the length of the [Expr] (ie. the iterable of expressions)
@@ -191,6 +163,120 @@ auto inits(Rng const& rng) {
          });
 }
 
+namespace bits {
+
+///
+/// @brief Generates all ordered bipartitions of the set bits in the mask.
+///
+/// For a given mask representing a set S, this function returns a view of pairs
+/// (A, B) such that A ∪ B = S, A ∩ B = {}, and A, B are subsets of S. The
+/// iteration order results in B descending and A ascending.
+///
+/// @tparam T The unsigned integer type of the mask.
+/// @param mask The bitmask representing the set.
+/// @return A view of pairs of disjoint submasks.
+///
+template <std::unsigned_integral T>
+constexpr auto bipartitions_ordered(T mask) {
+  // number of possible bipartitions
+  size_t const nparts = (1 << std::popcount(mask));
+  return std::views::iota(size_t{0}, nparts) |
+         std::views::transform([mask, p = mask](auto) mutable {
+           auto p_ = p;
+           p = (p - 1) & mask;
+           return std::pair<T, T>{p_ ^ mask, p_};
+         });
+}
+
+///
+/// @brief Generates unordered bipartitions of the mask.
+///
+/// This function returns a subset of the ordered bipartitions such that for
+/// every pair {A, B}, the pair {B, A} is excluded (unless A == B, which is
+/// handled naturally). Specifically, it takes the first half of the ordered
+/// bipartitions.
+///
+/// @param mask The bitmask representing the set.
+/// @return A view of unordered bipartitions.
+///
+constexpr auto bipartitions_unordered(std::unsigned_integral auto mask) {
+  auto bps = bipartitions_ordered(mask);
+  size_t const nparts = std::ranges::distance(bps) / 2;
+  return bps | std::views::take(nparts);
+}
+
+///
+/// @brief Alias for bipartitions_unordered.
+///
+/// @param mask The bitmask representing the set.
+/// @return A view of unordered bipartitions.
+///
+constexpr auto bipartitions(std::unsigned_integral auto mask) {
+  return bipartitions_unordered(mask);
+}
+
+///
+/// @brief Generates all subsets of the mask in ascending numerical order.
+///
+/// @param mask The bitmask representing the set.
+/// @return A view of submasks sorted ascending.
+///
+constexpr auto subsets_ascending(std::unsigned_integral auto mask) {
+  return bipartitions_ordered(mask) | std::views::elements<0>;
+}
+
+///
+/// @brief Generates all subsets of the mask in descending numerical order.
+///
+/// @param mask The bitmask representing the set.
+/// @return A view of submasks sorted descending.
+///
+constexpr auto subsets_descending(std::unsigned_integral auto mask) {
+  return bipartitions_ordered(mask) | std::views::elements<1>;
+}
+
+///
+/// @brief Generates all subsets of the mask in ascending numerical order.
+///
+/// @param mask The bitmask representing the set.
+/// @return A view of submasks.
+///
+constexpr auto subsets(std::unsigned_integral auto mask) {
+  return subsets_ascending(mask);
+}
+
+///
+/// @brief Generates indices of set bits in the mask.
+/// @param mask The bitmask representing the set.
+/// @return A view of indices of set bits.
+///
+constexpr auto on_bits_index(std::unsigned_integral auto mask) {
+  size_t const nbits = std::popcount(mask);
+  return std::views::iota(size_t{0}, nbits) |
+         std::views::transform([m = mask](auto) mutable {
+           size_t i = std::countr_zero(m);
+           m &= (m - 1);
+           return i;
+         });
+}
+
+///
+/// @brief Creates a view adaptor that selects elements from the input range
+/// based on indices.
+/// This function returns a closure object that, when applied to a range of
+/// indices, maps those indices to elements in the source range `rng`.
+/// @param rng The source range from which elements are selected.
+/// @return A view adaptor that transforms indices into elements from `rng`.
+///
+constexpr auto sieve(std::ranges::viewable_range auto&& rng) {
+  auto elems = std::views::all(std::forward<decltype(rng)>(rng));
+  return std::views::transform([elems](auto i) {
+    return *std::ranges::next(std::ranges::begin(elems), i);
+  });
+}
+
+}  // namespace bits
+
 }  // namespace sequant
 
 #endif  // SEQUANT_ALGORITHM_HPP

SeQuant/core/meta.hpp

@@ -9,6 +9,7 @@
 #include <memory>
 #include <range/v3/range/access.hpp>
 #include <range/v3/range/traits.hpp>
+#include <ranges>
 #include <type_traits>
 
 namespace sequant {
@@ -506,6 +507,23 @@ struct std_array_size<std::array<T, N>>
 template <typename T>
 constexpr inline std::size_t std_array_size_v = std_array_size<T>::value;
 
+///
+/// True if @tparam Rng is a range whose value type is convertible
+/// to @tparam V .
+///
+template <typename Rng, typename V>
+concept range_of = std::ranges::range<Rng> &&
+                   std::is_convertible_v<std::ranges::range_value_t<Rng>, V>;
+
+///
+/// True if @tparam Rng is a range whose value type is convertible
+/// to @tparam V and the size (std::ranges::distance) is known in constant time.
+///
+/// @see https://en.cppreference.com/w/cpp/ranges/sized_range.html
+///
+template <typename Rng, typename V>
+concept sized_range_of = std::ranges::sized_range<Rng> && range_of<Rng, V>;
+
 }  // namespace meta
 }  // namespace sequant
 

SeQuant/core/utility/indices.hpp

@@ -566,6 +566,46 @@ auto get_ket_idx(T&& container)
   }
 }
 
+auto subexpr_index_counts(meta::sized_range_of<Tensor> auto const& tnsrs) {
+  size_t const N = ranges::distance(tnsrs);
+  container::vector<container::map<Index, size_t, Index::FullLabelCompare>>
+      result((1 << N));
+  for (auto i = 1; i < result.size(); ++i) {
+    auto bs = boost::dynamic_bitset<>(N, i);
+    for (auto b = 0; b < N; ++b)
+      if (bs[b])
+        for (auto&& ix : ranges::at(tnsrs, b).indices())
+          if (auto [it, yn] = result[i].try_emplace(ix, 1); !yn)  //
+            ++(it->second);
+  }
+  return result;
+}
+
+template <meta::sized_range_of<Tensor> Rng, meta::range_of<Index> Ixs>
+auto subexpr_target_indices(Rng const& tnsrs, Ixs const& tixs) {
+  using IndexSet = container::set<Index, Index::FullLabelCompare>;
+  size_t const N = ranges::distance(tnsrs);
+  container::vector<IndexSet> result((1 << N));
+
+  if (result.empty()) return result;
+
+  for (auto i = 0; i < N; ++i) {
+    auto&& ixs = ranges::at(tnsrs, i).indices();
+    result[(1 << i)] = IndexSet(ranges::begin(ixs), ranges::end(ixs));
+  }
+
+  *result.rbegin() = IndexSet(ranges::begin(tixs), ranges::end(tixs));
+
+  auto counts = subexpr_index_counts(tnsrs);
+
+  for (auto i = 0; i < result.size(); ++i)
+    for (auto&& [k, v] : counts[i])
+      if (v == 1 || (v > 0 && counts.at(counts.size() - i - 1).contains(k)))
+        result[i].emplace(k);
+
+  return result;
+}
+
 }  // namespace sequant
 
 #endif  // SEQUANT_CORE_UTILITY_INDICES_HPP