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/eval/eval_expr.cpp

@@ -30,6 +30,11 @@
 
 namespace sequant {
 
+using EvalExprNode = FullBinaryNode<EvalExpr>;
+
+using IndexSet = container::set<Index, Index::FullLabelCompare>;
+using IndexCount = container::map<Index, size_t, Index::FullLabelCompare>;
+
 namespace {
 
 size_t hash_terminal_tensor(Tensor const&) noexcept;
@@ -300,7 +305,39 @@ struct ExprWithHash {
   size_t hash;
 };
 
-using EvalExprNode = FullBinaryNode<EvalExpr>;
+void index_count(IndexCount& counts, ExprPtr const& expr) {
+  if (!expr) return;
+  if (expr->is<Tensor>()) {
+    for (auto&& i : expr->as<Tensor>().const_indices()) {
+      auto&& [it, yn] = counts.emplace(i, 1);
+      if (!yn) it->second++;
+    }
+  } else if (expr->is<Sum>() && !expr->empty())
+    index_count(counts, expr->front());
+  else if (expr->is<Product>())
+    for (auto&& fac : *expr) index_count(counts, fac);
+}
+
+IndexCount index_count(meta::range_of<ExprPtr> auto const& facs) {
+  IndexCount result;
+  for (auto&& f : facs) index_count(result, f);
+  return result;
+}
+
+IndexSet target_indices(meta::range_of<ExprPtr> auto const& facs,
+                        IndexSet const& tidxs) {
+  IndexSet result;
+  for (auto&& [k, v] : index_count(facs))
+    if (v == 1 || tidxs.contains(k)) result.emplace(k);
+  return result;
+}
+
+auto imed_target_indices(Product const& prod, IndexSet const& tidxs) {
+  return inits(prod.factors()) |
+         ranges::views::transform(
+             [&tidxs](auto&& facs) { return target_indices(facs, tidxs); }) |
+         ranges::to_vector;
+}
 
 ///
 /// \brief Collect tensors appearing as a factor at the leaf node of a product

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,52 @@ 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 @p T is a range of rank @p Rank whose value type is convertible to
+/// @p V.
+///
+/// A rank of 0 means @p T is directly convertible to @p V.
+/// A rank of 1 means @p T is a range of values convertible to @p V.
+/// A rank of N means @p T is a range of ranges ... (N times) of values
+/// convertible to @p V.
+///
+template <typename, typename, size_t>
+constexpr bool range_of_rank{};
+
+template <typename T, typename V>
+constexpr bool range_of_rank<T, V, 0> = std::is_convertible_v<T, V>;
+
+template <std::ranges::range Rng, typename V, size_t Rank>
+constexpr bool range_of_rank<Rng, V, Rank> =
+    range_of_rank<std::ranges::range_value_t<Rng>, V, (Rank - 1)>;
+
+///
+/// True if @tparam Rng is a range whose value type is convertible
+/// to @tparam V .
+///
+/// A rank of 0 means @tparam Rng is directly convertible to @tparam V.
+/// A rank of 1 means @tparam Rng is a range of values convertible to @tparam V.
+/// A rank of N means @tparam Rng is a range of ranges ... (N times) of values
+/// convertible to @tparam V.
+///
+template <typename Rng, typename V, size_t Rank = 1>
+concept range_of = range_of_rank<Rng, V, Rank>;
+
+///
+/// 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.
+///
+/// A rank of 0 means @tparam Rng is directly convertible to @tparam V.
+/// A rank of 1 means @tparam Rng is a range of values convertible to @tparam V.
+/// A rank of N means @tparam Rng is a range of ranges ... (N times) of values
+/// convertible to @tparam V.
+///
+/// @see https://en.cppreference.com/w/cpp/ranges/sized_range.html
+///
+template <typename Rng, typename V, size_t Rank = 1>
+concept sized_range_of =
+    std::ranges::sized_range<Rng> && range_of<Rng, V, Rank>;
+
 }  // namespace meta
 }  // namespace sequant