diff --git a/CMakeLists.txt b/CMakeLists.txt index b04b376fc9..2918b52087 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -367,6 +367,8 @@ set(SeQuant_symb_src SeQuant/core/tensor_network/v2.hpp SeQuant/core/tensor_network/v3.cpp SeQuant/core/tensor_network/v3.hpp + SeQuant/core/tree_index.cpp + SeQuant/core/tree_index.hpp SeQuant/core/utility/aggregate.hpp SeQuant/core/utility/atomic.hpp SeQuant/core/utility/context.hpp diff --git a/SeQuant/core/attr.hpp b/SeQuant/core/attr.hpp index ad39ae8efb..ef8cc8be33 100644 --- a/SeQuant/core/attr.hpp +++ b/SeQuant/core/attr.hpp @@ -11,6 +11,7 @@ #include #include #include +#include namespace sequant { @@ -130,10 +131,10 @@ enum class BraKetPos { /// /// @note This does not include slot bundles, like braket, etc. enum class SlotType { - Bra, - Ket, - Aux, - Proto, + Bra = 0b1, + Ket = 0b10, + Aux = 0b100, + Proto = 0b1000, }; template @@ -156,6 +157,33 @@ std::basic_ostream& operator<<( return stream; } +struct SlotTypes { + std::underlying_type_t active = 0; + + constexpr SlotTypes(SlotType type) + : active(static_cast(type)) {} + constexpr SlotTypes(decltype(active) val) : active(val) {} + + constexpr bool operator==(const SlotTypes&) const = default; + constexpr auto operator<=>(const SlotTypes&) const = default; + + constexpr bool operator&(SlotType type) const { + return active & static_cast(type); + } + + constexpr SlotTypes operator|(SlotType type) const { + return {active | static_cast(type)}; + } +}; + +constexpr SlotTypes operator|(SlotType lhs, SlotType rhs) { + using IntType = std::underlying_type_t; + return {static_cast(lhs) | static_cast(rhs)}; +} + +static constexpr const SlotTypes AnySlotType = + SlotType::Bra | SlotType::Ket | SlotType::Aux | SlotType::Proto; + enum class Statistics { FermiDirac, BoseEinstein, diff --git a/SeQuant/core/export/export.hpp b/SeQuant/core/export/export.hpp index 4b3627218a..2dda1066a3 100644 --- a/SeQuant/core/export/export.hpp +++ b/SeQuant/core/export/export.hpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -798,7 +799,8 @@ void preprocess_and_maybe_log(ExportNode &tree, PreprocessResult &result, std::cout << "Tree before preprocessing:\n" << tree.tikz( [](const ExportNode &node) { - return "$" + toUtf8(to_latex(node->expr())) + "$"; + return "$" + toUtf8(io::latex::to_string(node->expr())) + + "$"; }, [](const ExportNode) -> std::string { return ""; }) << "\n"; @@ -811,7 +813,8 @@ void preprocess_and_maybe_log(ExportNode &tree, PreprocessResult &result, std::cout << "Tree after pre-processing:\n" << tree.tikz( [](const ExportNode &node) { - return "$" + toUtf8(to_latex(node->expr())) + "$"; + return "$" + toUtf8(io::latex::to_string(node->expr())) + + "$"; }, [](const ExportNode) -> std::string { return ""; }) << "\n"; diff --git a/SeQuant/core/expressions/expr.cpp b/SeQuant/core/expressions/expr.cpp index a8ce23862e..bbc30d66cc 100644 --- a/SeQuant/core/expressions/expr.cpp +++ b/SeQuant/core/expressions/expr.cpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include @@ -65,6 +66,14 @@ const ExprPtr &Expr::operator[](std::size_t idx) const { return begin()[idx]; } +ExprPtr &ExprPtr::operator[](const TreeIndex &idx) { + return idx.select_from(*this); +} + +const ExprPtr &ExprPtr::operator[](const TreeIndex &idx) const { + return idx.select_from(*this); +} + ExprPtr &Expr::at(std::size_t idx) { return (*this)[idx]; } const ExprPtr &Expr::at(std::size_t idx) const { return (*this)[idx]; } @@ -183,6 +192,12 @@ Expr &Expr::operator+=(const Expr &) { throw not_implemented("operator+="); } Expr &Expr::operator-=(const Expr &) { throw not_implemented("operator-="); } +Expr &Expr::operator[](const TreeIndex &idx) { return idx.select_from(*this); } + +const Expr &Expr::operator[](const TreeIndex &idx) const { + return idx.select_from(*this); +} + ExprPtr adjoint(const ExprPtr &expr) { auto result = expr->clone(); result->adjoint(); diff --git a/SeQuant/core/expressions/expr.hpp b/SeQuant/core/expressions/expr.hpp index 7bee7e8bc3..372ad50087 100644 --- a/SeQuant/core/expressions/expr.hpp +++ b/SeQuant/core/expressions/expr.hpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include @@ -381,6 +382,10 @@ class Expr : public std::enable_shared_from_this { ExprPtr &operator[](std::size_t idx); const ExprPtr &operator[](std::size_t idx) const; + /// @return The subexpression identified by the given index + Expr &operator[](const TreeIndex &idx); + const Expr &operator[](const TreeIndex &idx) const; + ExprPtr &at(std::size_t idx); const ExprPtr &at(std::size_t idx) const; diff --git a/SeQuant/core/expressions/expr_ptr.hpp b/SeQuant/core/expressions/expr_ptr.hpp index 5eb532d0d0..fa140ed86c 100644 --- a/SeQuant/core/expressions/expr_ptr.hpp +++ b/SeQuant/core/expressions/expr_ptr.hpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include @@ -111,6 +112,10 @@ class ExprPtr : public std::shared_ptr { /// @return reference to @c *this ExprPtr &operator*=(const ExprPtr &); + /// @return The subexpression identified by the given index + ExprPtr &operator[](const TreeIndex &idx); + const ExprPtr &operator[](const TreeIndex &idx) const; + /// @tparam T an Expr type /// @return true if this object is of type @c T template diff --git a/SeQuant/core/io/latex/latex.cpp b/SeQuant/core/io/latex/latex.cpp index 74e37f60a4..00bf83508d 100644 --- a/SeQuant/core/io/latex/latex.cpp +++ b/SeQuant/core/io/latex/latex.cpp @@ -5,12 +5,15 @@ #include #include #include +#include #include #include #include #include #include +#include +#include #include namespace sequant::io::latex { @@ -64,6 +67,19 @@ std::wstring to_string(const Power& power) { return result; } +std::wstring to_string(const ResultExpr& expr) { + std::wstringstream stream; + if (expr.produces_tensor()) { + stream << to_string(expr.result_as_tensor()); + } else { + stream << to_string(expr.result_as_variable()); + } + + stream << " = " << to_string(expr.expression()); + + return stream.str(); +} + namespace detail { template diff --git a/SeQuant/core/io/latex/latex.hpp b/SeQuant/core/io/latex/latex.hpp index 1a0dbfa237..bf18fcfffa 100644 --- a/SeQuant/core/io/latex/latex.hpp +++ b/SeQuant/core/io/latex/latex.hpp @@ -19,7 +19,8 @@ namespace sequant { class Power; -} +class ResultExpr; +} // namespace sequant namespace sequant::io::latex { @@ -105,6 +106,8 @@ std::wstring to_string(const rational& num); std::wstring to_string(const Power& power); +std::wstring to_string(const ResultExpr& expr); + namespace detail { template diff --git a/SeQuant/core/meta.hpp b/SeQuant/core/meta.hpp index 6d57715fc0..404e76e3f7 100644 --- a/SeQuant/core/meta.hpp +++ b/SeQuant/core/meta.hpp @@ -9,6 +9,7 @@ #include #include #include +#include #include #include @@ -570,6 +571,27 @@ static_assert( std::same_as, const float &>); static_assert(std::same_as, float &>); +#if defined(__cpp_lib_forward_like) && __cpp_lib_forward_like >= 202207L +using std::forward_like; +#else +template +// Taken from cppreference https://en.cppreference.com/cpp/utility/forward_like +constexpr auto &&forward_like(U &&x) noexcept { + constexpr bool is_adding_const = std::is_const_v>; + if constexpr (std::is_lvalue_reference_v) { + if constexpr (is_adding_const) + return std::as_const(x); + else + return static_cast(x); + } else { + if constexpr (is_adding_const) + return std::move(std::as_const(x)); + else + return std::move(x); + } +} +#endif + /// /// True if @p T is a range of rank @p Rank whose value type is convertible to /// @p V. diff --git a/SeQuant/core/tree_index.cpp b/SeQuant/core/tree_index.cpp new file mode 100644 index 0000000000..81ed9e9932 --- /dev/null +++ b/SeQuant/core/tree_index.cpp @@ -0,0 +1,71 @@ +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace sequant { + +TreeIndex::TreeIndex(container::svector positions) + : positions_(std::move(positions)) {} + +TreeIndex::TreeIndex(std::initializer_list positions) + : TreeIndex(container::svector(std::move(positions))) {} + +template +ExprType &select_impl(ExprType &expr, + const container::svector &positions) { + constexpr bool is_expr_ptr = + std::same_as, ExprPtr>; + + ExprType *selected = &expr; + + for (std::size_t current : positions) { + if (current >= selected->size()) { + throw Exception("Position " + std::to_string(current) + + " is out of bounds for an expression of dimension " + + std::to_string(selected->size())); + } + + using std::ranges::begin; + + auto it = begin(*selected) + current; + + if constexpr (is_expr_ptr) { + selected = &(*it); + } else { + selected = &(*(*it)); + } + } + + SEQUANT_ASSERT(selected); + + return *selected; +} + +ExprPtr &TreeIndex::select_from(ExprPtr &expr) const { + return select_impl(expr, positions_); +} + +const ExprPtr &TreeIndex::select_from(const ExprPtr &expr) const { + return select_impl(expr, positions_); +} + +Expr &TreeIndex::select_from(Expr &expr) const { + return select_impl(expr, positions_); +} + +const Expr &TreeIndex::select_from(const Expr &expr) const { + return select_impl(expr, positions_); +} + +std::size_t TreeIndex::depth() const { return positions_.size(); } + +} // namespace sequant diff --git a/SeQuant/core/tree_index.hpp b/SeQuant/core/tree_index.hpp new file mode 100644 index 0000000000..ddb9c521be --- /dev/null +++ b/SeQuant/core/tree_index.hpp @@ -0,0 +1,75 @@ +#ifndef SEQUANT_TREEINDEX_H +#define SEQUANT_TREEINDEX_H + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace sequant { + +/// Represents an index (position) into a tree structure +/// @note Using this to select a subexpression from an expression tree takes +/// O(N) time where N = depth() +class TreeIndex { + public: + TreeIndex() = default; + + TreeIndex(container::svector positions); + + TreeIndex(std::initializer_list positions); + + template + requires(std::integral>) + TreeIndex(Positions &&positions) { + if constexpr (std::ranges::sized_range) { + positions_.reserve(std::ranges::size(positions)); + } + + for (auto val : positions) { + SEQUANT_ASSERT(val >= 0); + positions_.push_back(val); + } + } + + ExprPtr &select_from(ExprPtr &expr) const; + const ExprPtr &select_from(const ExprPtr &expr) const; + Expr &select_from(Expr &expr) const; + const Expr &select_from(const Expr &expr) const; + + std::size_t depth() const; + + bool operator==(const TreeIndex &) const = default; + std::strong_ordering operator<=>(const TreeIndex &) const = default; + + template + friend std::basic_ostream &operator<<( + std::basic_ostream &stream, const TreeIndex &idx) { + stream << "{"; + + for (std::size_t i = 0; i < idx.positions_.size(); ++i) { + stream << " " << idx.positions_[i] << " "; + + if (i + 1 < idx.positions_.size()) { + stream << "->"; + } + } + + stream << "}"; + + return stream; + } + + private: + container::svector positions_; +}; + +} // namespace sequant + +#endif // SEQUANT_TREEINDEX_H diff --git a/SeQuant/core/utility/expr.hpp b/SeQuant/core/utility/expr.hpp index 95e6061b33..7ba3795e67 100644 --- a/SeQuant/core/utility/expr.hpp +++ b/SeQuant/core/utility/expr.hpp @@ -71,7 +71,7 @@ bool is_valid(const ResultExpr &expr, std::string *msg = nullptr); /// /// @param expression The expression to modify /// @param label The label of the tensor that shall be removed -/// @returns The removed tensor, if any occurrance has been found +/// @returns The removed tensor, if any occurrence has been found std::optional pop_tensor(ExprPtr &expression, std::wstring_view label); /// Replaces a given target expression by a given replacement diff --git a/SeQuant/core/utility/indices.hpp b/SeQuant/core/utility/indices.hpp index f25b1ecad2..27f6eae317 100644 --- a/SeQuant/core/utility/indices.hpp +++ b/SeQuant/core/utility/indices.hpp @@ -364,19 +364,54 @@ Map get_used_indices_with_counts(const ExprPtr& expr) { /// @returns A set of all unique indices used in the provided expression /// @note includes pure protoindices -template > +/// @tparam types The set of SlotTypes of indices to consider. Any index that +/// never occurs in these kind of slots at least once is filtered out +template , SlotTypes types = AnySlotType> Set get_used_indices(const Expr& expr) { - return get_used_indices_with_counts(expr) | - std::views::transform( - [](const auto& idx_count) { return idx_count.first; }) | - ranges::to; + if constexpr (types == AnySlotType) { + return get_used_indices_with_counts(expr) | + std::views::transform( + [](const auto& idx_count) { return idx_count.first; }) | + ranges::to; + } else { + return get_used_indices_with_counts(expr) | + std::views::filter([](const auto& pair) -> bool { + if constexpr (types & SlotType::Bra) { + if (pair.second.bra > 0) { + return true; + } + } + if constexpr (types & SlotType::Ket) { + if (pair.second.ket > 0) { + return true; + } + } + if constexpr (types & SlotType::Aux) { + if (pair.second.aux > 0) { + return true; + } + } + if constexpr (types & SlotType::Proto) { + if (pair.second.proto > 0) { + return true; + } + } + + return false; + }) | + std::views::transform( + [](const auto& idx_count) { return idx_count.first; }) | + ranges::to; + } } /// @returns A set of all unique indices used in the provided expression /// @note includes pure protoindices -template > +/// @tparam types The set of SlotTypes of indices to consider. Any index that +/// never occurs in these kind of slots at least once is filtered out +template , SlotTypes types = AnySlotType> Set get_used_indices(const ExprPtr& expr) { - return get_used_indices(*expr); + return get_used_indices(*expr); } template , typename Rng> @@ -719,8 +754,11 @@ decltype(auto) get_ket_idx(Group&& group) { /// (based on slot type) or assert that they already are in canonical order, /// if group is const. decltype(auto) as_index_group_view(SlottedIndexGroup auto&& group) { - static_assert(static_cast(SlotType::Bra) == 0); - static_assert(static_cast(SlotType::Ket) == 1); + static_assert(SlotType::Bra < SlotType::Ket); + static_assert(SlotType::Bra < SlotType::Aux); + static_assert(SlotType::Bra < SlotType::Proto); + static_assert(SlotType::Ket < SlotType::Aux); + static_assert(SlotType::Ket < SlotType::Proto); // We have to ensure a unique order of indices if we're getting rid of the // SlotType tag diff --git a/SeQuant/domain/mbpt/spin.cpp b/SeQuant/domain/mbpt/spin.cpp index a182e40374..9dec9a9aec 100644 --- a/SeQuant/domain/mbpt/spin.cpp +++ b/SeQuant/domain/mbpt/spin.cpp @@ -44,6 +44,8 @@ #include #include #include +#include +#include #include #include #include @@ -619,51 +621,63 @@ ExprPtr symmetrize_expr(const ProductPtr& product) { return factor->is() && factor->as().label() == reserved::antisymm_label(); }); - if (it == ranges::end(factors)) return product; - const auto& A_tensor = (*it)->as(); + + if (it == ranges::end(factors)) { + return product; + } + + const Tensor& A_tensor = (*it)->as(); SEQUANT_ASSERT(A_tensor.label() == reserved::antisymm_label()); auto A_is_nconserving = A_tensor.bra_rank() == A_tensor.ket_rank(); - if (A_is_nconserving && A_tensor.bra_rank() == 1) - return remove_tensor(product, reserved::antisymm_label()); + SEQUANT_ASSERT(A_tensor.aux_rank() == 0); - SEQUANT_ASSERT(A_tensor.rank() > 1); + if (A_is_nconserving && A_tensor.rank() == 1) { + return remove_tensor(product, reserved::antisymm_label()); + } - auto S = Tensor{}; + std::optional S; if (A_is_nconserving) { + SEQUANT_ASSERT(A_tensor.rank() > 1); S = Tensor(reserved::symm_label(), A_tensor.bra(), A_tensor.ket(), A_tensor.aux(), Symmetry::Nonsymm); - } else { // A is N-nonconserving + } else { auto n = std::min(A_tensor.bra_rank(), A_tensor.ket_rank()); - container::svector bra_list(A_tensor.bra().begin(), - A_tensor.bra().begin() + n); - container::svector ket_list(A_tensor.ket().begin(), - A_tensor.ket().begin() + n); - S = Tensor(reserved::symm_label(), bra(std::move(bra_list)), - ket(std::move(ket_list)), A_tensor.aux(), Symmetry::Nonsymm); + + if (n > 0) { + container::svector bra_list(A_tensor.bra().begin(), + A_tensor.bra().begin() + n); + container::svector ket_list(A_tensor.ket().begin(), + A_tensor.ket().begin() + n); + S = Tensor(reserved::symm_label(), bra(std::move(bra_list)), + ket(std::move(ket_list)), A_tensor.aux(), Symmetry::Nonsymm); + } } - const auto nf = rational{1, factorial(S.ket_rank())}; + + const auto nf = S ? rational{1, factorial(S->ket_rank())} : 1; // Generate replacement maps from a list of Index type (could be a bra or a // ket) // Uses a permuted list of int to generate permutations // TODO factor out for reuse auto maps_from_list = [](const container::svector& list) { - container::svector int_list(list.size()); + container::svector int_list(list.size()); std::iota(int_list.begin(), int_list.end(), 0); container::svector> result; do { container::map map; auto list_ptr = list.begin(); - for (auto&& i : int_list) { + for (std::size_t i : int_list) { map.emplace(*list_ptr, list[i]); list_ptr++; } result.push_back(map); - } while (std::next_permutation(int_list.begin(), int_list.end())); + } while (std::ranges::next_permutation(int_list).found); + SEQUANT_ASSERT(result.size() == boost::numeric_cast(factorial(list.size()))); + return result; }; @@ -671,7 +685,9 @@ ExprPtr symmetrize_expr(const ProductPtr& product) { // TODO factor out for reuse auto get_phase = [](const container::map& map) { container::svector idx_list; - for (const auto& [key, val] : map) idx_list.push_back(val); + auto indices = map | std::ranges::views::values; + idx_list.insert(idx_list.end(), indices.begin(), indices.end()); + reset_ts_swap_counter(); bubble_sort(std::begin(idx_list), std::end(idx_list)); return ts_swap_counter_is_even() ? 1 : -1; @@ -683,15 +699,21 @@ ExprPtr symmetrize_expr(const ProductPtr& product) { maps = maps_from_list(A_tensor.bra()); } else { SEQUANT_ASSERT(A_tensor.bra_rank() != A_tensor.ket_rank()); + maps = A_tensor.bra_rank() > A_tensor.ket_rank() ? maps_from_list(A_tensor.bra()) : maps_from_list(A_tensor.ket()); } + SEQUANT_ASSERT(!maps.empty()); + for (auto&& map : maps) { Product new_product{}; - new_product.scale(product->scalar() * nf); - new_product.append(get_phase(map), ex(S)); + new_product.scale(product->scalar() * nf * get_phase(map)); + if (S) { + new_product.append(ex(S.value())); + } + auto temp_product = remove_tensor(product, reserved::antisymm_label()); for (auto&& term : *temp_product) { if (term->is()) { @@ -705,8 +727,12 @@ ExprPtr symmetrize_expr(const ProductPtr& product) { term->type_name()); } } + result->append(ex(new_product)); - } // map + } + + detail::reset_idx_tags(result); + return result; } @@ -938,24 +964,19 @@ ExprPtr closed_shell_spintrace_impl(const ExprPtr& expression, // non-symmetric. // full_expansion: it fully expands the antisymmetrizer directly (can be used // for v2 eqs, however it is not an optimized way). - auto partially_or_fully_expand = [&full_expansion](const ExprPtr& expr) { - auto temp = expr; - if (has_tensor(temp, reserved::antisymm_label())) { - if (full_expansion) { - temp = expand_A_op(temp); - } else { - temp = symmetrize_expr(temp); - } + ExprPtr expr = expression; + expand(expr); + + if (has_tensor(expr, reserved::antisymm_label())) { + if (full_expansion) { + expr = expand_A_op(expr); + } else { + expr = symmetrize_expr(expr); } - temp = expand_antisymm(temp); - rapid_simplify(temp); - return temp; - }; - ExprPtr expr = partially_or_fully_expand(expression); + } + + expr = expand_antisymm(expr); - // Index tags are cleaned prior to calling the fast canonicalizer - detail::reset_idx_tags(expr); // This call is REQUIRED - expand(expr); // This call is REQUIRED simplify(expr); // full simplify to combine terms before count_cycles // Lambda for spin-tracing a product term @@ -1376,7 +1397,8 @@ std::vector open_shell_spintrace_impl( // Grand index list contains both internal and external indices container::set grand_idxlist = - get_used_indices(expr); + get_used_indices(expr); container::set ext_idxlist; for (const auto& idxgrp : ext_index_groups) { @@ -1646,7 +1668,7 @@ std::vector open_shell_CC_spintrace(const ExprPtr& expr) { return expr_vec; } -template +template ExprPtr spintrace_impl(const ExprPtr& expression, IdxGroups&& ext_index_groups, bool spinfree_index_spaces) { // Escape immediately if expression is a constant @@ -1654,7 +1676,7 @@ ExprPtr spintrace_impl(const ExprPtr& expression, IdxGroups&& ext_index_groups, return expression; } - if constexpr (assert_enabled()) { + if constexpr (check_ext_indices && assert_enabled()) { // Verify that the number of external indices matches the number of indices // in ext_index_groups, UNLESS user overrode external definitions in default // context @@ -1727,7 +1749,8 @@ ExprPtr spintrace_impl(const ExprPtr& expression, IdxGroups&& ext_index_groups, ExprPtr expr = product->clone(); // List of all indices in the expression container::set grand_idxlist = - get_used_indices(expr); + get_used_indices(expr); // List of external indices, i.e. indices that are not summed over Einstein // style (indices that are not repeated in an expression) @@ -1829,6 +1852,8 @@ ExprPtr spintrace_impl(const ExprPtr& expression, IdxGroups&& ext_index_groups, // Expand antisymmetrizer operator (A) if present in the expression ExprPtr expr = expression; + expand(expr); + if (has_tensor(expr, reserved::antisymm_label())) expr = expand_A_op(expr); if (expr->is()) expr = ex(1) * expr; @@ -1849,12 +1874,16 @@ ExprPtr spintrace_impl(const ExprPtr& expression, IdxGroups&& ext_index_groups, result_sum->append(term); result = result_sum; } - return result; } else { throw Exception("Invalid Expr type in spintrace: " + expr->type_name()); } + SEQUANT_ASSERT(result); + detail::reset_idx_tags(result); + + simplify(result); + return result; } @@ -1862,32 +1891,40 @@ ExprPtr spintrace(const ExprPtr& expression, const container::svector>& ext_index_groups, bool spinfree_index_spaces) { - return spintrace_impl(expression, as_view_of_index_groups(ext_index_groups), - spinfree_index_spaces); + return spintrace_impl(expression, + as_view_of_index_groups(ext_index_groups), + spinfree_index_spaces); } ExprPtr spintrace(const ExprPtr& expression, EmptyInitializerList, bool spinfree_index_spaces) { - return spintrace_impl(expression, - container::svector>{}, - spinfree_index_spaces); + return spintrace_impl(expression, + container::svector>{}, + spinfree_index_spaces); } ExprPtr spintrace( const ExprPtr& expression, const container::svector>& ext_index_groups, bool spinfree_index_spaces) { - return spintrace_impl(expression, ext_index_groups, spinfree_index_spaces); + return spintrace_impl(expression, ext_index_groups, + spinfree_index_spaces); +} + +ExprPtr resultexpr_spintrace_delegate( + const ExprPtr& expression, + const container::svector>& + ext_index_groups, + bool spinfree_index_spaces) { + return spintrace_impl(expression, + as_view_of_index_groups(ext_index_groups), + spinfree_index_spaces); } container::svector spintrace(const ResultExpr& expr, bool spinfree_index_spaces) { - using TraceFunction = ExprPtr (*)( - const ExprPtr&, - const container::svector>&, bool); - return detail::wrap_trace>( - expr, static_cast(&spintrace), spinfree_index_spaces); + expr, &resultexpr_spintrace_delegate, spinfree_index_spaces); } } // namespace sequant::mbpt diff --git a/tests/unit/catch2_sequant.hpp b/tests/unit/catch2_sequant.hpp index 3dc0185ef1..1346c25b5c 100644 --- a/tests/unit/catch2_sequant.hpp +++ b/tests/unit/catch2_sequant.hpp @@ -16,6 +16,7 @@ #include #include +#include #include #include #include @@ -182,7 +183,9 @@ ExprVar to_expression(T &&expression) { using std::begin; using std::end; - if constexpr (std::is_convertible_v) { + using BaseT = std::remove_cvref_t; + + if constexpr (std::is_convertible_v) { std::wstring string = sequant::toUtf16(std::forward(expression)); if (std::find(begin(string), end(string), L'=') != end(string)) { @@ -194,7 +197,7 @@ ExprVar to_expression(T &&expression) { std::string(std::forward(expression)), {.def_perm_symm = sequant::Symmetry::Nonsymm}); } - } else if constexpr (std::is_convertible_v) { + } else if constexpr (std::is_convertible_v) { if (std::find(begin(expression), end(expression), L'=') != end(expression)) { return sequant::deserialize( @@ -205,13 +208,13 @@ ExprVar to_expression(T &&expression) { std::wstring(std::forward(expression)), {.def_perm_symm = sequant::Symmetry::Nonsymm}); } - } else if constexpr (std::is_convertible_v) { + } else if constexpr (std::same_as) { return expression; - } else if constexpr (std::is_convertible_v) { + } else if constexpr (std::same_as) { // Clone in order to not have to worry about later modification return expression.clone(); } else { - static_assert(std::is_convertible_v, + static_assert(std::same_as, "Invalid type for expression"); // Clone in order to not have to worry about later modification diff --git a/tests/unit/test_expr.cpp b/tests/unit/test_expr.cpp index b9112c7ee6..4a9e807e82 100644 --- a/tests/unit/test_expr.cpp +++ b/tests/unit/test_expr.cpp @@ -13,6 +13,7 @@ #include #include #include +#include #include #include @@ -25,6 +26,7 @@ #include #include #include +#include #include #include #include @@ -1152,4 +1154,44 @@ TEST_CASE("expr", "[elements]") { REQUIRE_THAT(pairings, ::Catch::Matchers::UnorderedRangeEquals(expected)); } } + + SECTION("TreeIndex") { + std::vector> tests = { + {"Var", {}, "Var"}, {"A * B", {}, "A * B"}, + {"A + B", {0}, "A"}, {"A + B", {1}, "B"}, + {"A + B", {}, "A + B"}, {"A + B", {0}, "A"}, + {"A + B", {1}, "B"}, {"A + B * C", {0}, "A"}, + {"A + B * C", {1}, "B * C"}, {"A + B * C", {1, 0}, "B"}, + {"A + B * C", {1, 1}, "C"}, {"(A + B) * C", {0}, "A + B"}, + {"(A + B) * C", {0, 1}, "B"}, {"(A + B) * C", {1}, "C"}, + }; + + for (const auto &[expr_string, idx, expected_str] : tests) { + CAPTURE(expr_string); + CAPTURE(idx); + CAPTURE(expected_str); + + ExprPtr expr = deserialize(expr_string); + ExprPtr expected = deserialize(expected_str); + + REQUIRE_THAT(idx.select_from(expr), EquivalentTo(expected)); + REQUIRE_THAT(idx.select_from(std::as_const(expr)), + EquivalentTo(expected)); + REQUIRE_THAT(idx.select_from(*expr), EquivalentTo(*expected)); + REQUIRE_THAT(idx.select_from(std::as_const(*expr)), + EquivalentTo(*expected)); + + REQUIRE_THAT(expr[idx], EquivalentTo(*expected)); + REQUIRE_THAT(std::as_const(expr)[idx], EquivalentTo(*expected)); + REQUIRE_THAT((*expr)[idx], EquivalentTo(*expected)); + REQUIRE_THAT(std::as_const(*expr)[idx], EquivalentTo(*expected)); + } + + SECTION("Out-of-bounds exception") { + const ExprPtr expr = deserialize("A * B"); + + REQUIRE_THROWS_AS(TreeIndex({2}).select_from(expr), Exception); + REQUIRE_THROWS_AS(TreeIndex({0, 1}).select_from(expr), Exception); + } + } } diff --git a/tests/unit/test_index.cpp b/tests/unit/test_index.cpp index 998f3c8c48..2cff127a76 100644 --- a/tests/unit/test_index.cpp +++ b/tests/unit/test_index.cpp @@ -64,8 +64,7 @@ TEST_CASE("index", "[elements][index]") { Index i2_alt(isr->retrieve("i1"), 2); REQUIRE(i2 == i2_alt); - // examples with proto indices - { + SECTION("with proto-indices") { REQUIRE_NOTHROW(Index(isr->retrieve(L"i"), 3, {i1, i2})); Index i3(isr->retrieve(L"i"), 3, {i1, i2}); REQUIRE(i3.label() == L"i_3"); @@ -142,8 +141,7 @@ TEST_CASE("index", "[elements][index]") { } } - // can use bytestrings also - { + SECTION("from bytestring") { Index i1("i_1"); REQUIRE(i1.label() == L"i_1"); REQUIRE(i1.space() == isr->retrieve("i")); @@ -165,8 +163,7 @@ TEST_CASE("index", "[elements][index]") { get_default_context().index_space_registry()->retrieve("α_1")); } - // with default Context label defines Index - { + SECTION("default Context label defines Index") { const auto ctx_resetter = set_scoped_default_context(Context{}); REQUIRE_NOTHROW(Index(L"i")); Index i(L"i"); diff --git a/tests/unit/test_spin.cpp b/tests/unit/test_spin.cpp index 24091c6a53..ce102c7b7a 100644 --- a/tests/unit/test_spin.cpp +++ b/tests/unit/test_spin.cpp @@ -270,6 +270,16 @@ TEST_CASE("spin", "[spin]") { REQUIRE_THAT(result, EquivalentTo("-1/2 g{p1,p2;p4,p3} + g{p1,p2;p3,p4}")); } + { + // Ensure auxiliary indices are ignored + const ExprPtr expr = deserialize("1/4 g{p1,p2;p3,p4;p5} d{;;p5}", + {.def_perm_symm = Symmetry::Antisymm}); + auto result = spintrace(expr, IdxGroupList{{"p1", "p3"}, {"p2", "p4"}}); + REQUIRE_THAT( + result, + EquivalentTo( + "-1/2 g{p1,p2;p4,p3;p5} d{;;p5} + g{p1,p2;p3,p4;p5} d{;;p5}")); + } { // Note the provided external index pairings which is different from the // way this tensor is written down Also, the prefactor of -1 is important diff --git a/tests/unit/test_utilities.cpp b/tests/unit/test_utilities.cpp index a8b7130bc3..980102780c 100644 --- a/tests/unit/test_utilities.cpp +++ b/tests/unit/test_utilities.cpp @@ -397,6 +397,45 @@ TEST_CASE("utilities", "[utilities]") { REQUIRE_THAT(get_used_indices(expr), Catch::Matchers::UnorderedRangeEquals(indices)); } + + SECTION("restricted") { + const ExprPtr expr = deserialize("A{a1;a2;a3} B{i1;i3} C{i3}"); + + std::vector expected = {"a1", Index("i1", {"i2"}), "i3"}; + std::vector actual = + get_used_indices, SlotType::Bra>(expr); + REQUIRE_THAT(actual, Catch::Matchers::UnorderedRangeEquals(expected)); + + expected = {"a2", "i3"}; + actual = get_used_indices, SlotType::Ket>(expr); + REQUIRE_THAT(actual, Catch::Matchers::UnorderedRangeEquals(expected)); + + expected = {"a3"}; + actual = get_used_indices, SlotType::Aux>(expr); + REQUIRE_THAT(actual, Catch::Matchers::UnorderedRangeEquals(expected)); + + expected = {"i2"}; + actual = get_used_indices, SlotType::Proto>(expr); + REQUIRE_THAT(actual, Catch::Matchers::UnorderedRangeEquals(expected)); + + expected = {"i2", "a3"}; + actual = + get_used_indices, SlotType::Proto | SlotType::Aux>( + expr); + REQUIRE_THAT(actual, Catch::Matchers::UnorderedRangeEquals(expected)); + + expected = {"i2", "a3", "a2", "i3"}; + actual = + get_used_indices, SlotType::Proto | SlotType::Aux | + SlotType::Ket>(expr); + REQUIRE_THAT(actual, Catch::Matchers::UnorderedRangeEquals(expected)); + + expected = {"i2", "a3", "a2", "i3", "a1", Index("i1", {"i2"})}; + actual = get_used_indices, + SlotType::Proto | SlotType::Aux | + SlotType::Ket | SlotType::Bra>(expr); + REQUIRE_THAT(actual, Catch::Matchers::UnorderedRangeEquals(expected)); + } } SECTION("replace") { diff --git a/utilities/external-interface/CMakeLists.txt b/utilities/external-interface/CMakeLists.txt index cb1071235d..fe8d063a5f 100644 --- a/utilities/external-interface/CMakeLists.txt +++ b/utilities/external-interface/CMakeLists.txt @@ -8,6 +8,22 @@ add_executable(external_interface external_interface.cpp processing.cpp utils.cpp + + canonicalize_step.cpp + density_fitting_step.cpp + execution_context.cpp + executor.cpp + export_step.cpp + optimization_step.cpp + output_step.cpp + processing_step.cpp + processing_step_factory.cpp + projection_step.cpp + read_input_step.cpp + simplify_step.cpp + spintracing_step.cpp + to_export_tree_step.cpp + validate_step.cpp ) set_target_properties(external_interface PROPERTIES CXX_SCAN_FOR_MODULES OFF) diff --git a/utilities/external-interface/canonicalize_step.cpp b/utilities/external-interface/canonicalize_step.cpp new file mode 100644 index 0000000000..d8b5311780 --- /dev/null +++ b/utilities/external-interface/canonicalize_step.cpp @@ -0,0 +1,54 @@ +#include "canonicalize_step.hpp" +#include "processing_data.hpp" +#include "processing_step_factory.hpp" + +#include + +#include + +#include +#include + +namespace sequant::util::extint { + +SEQUANT_EXTINT_REGISTER_STEP_TYPE(CanonicalizeStep, "canonicalize"); + +std::string CanonicalizeStep::kind() const { return "canonicalize"; } + +bool CanonicalizeStep::accepts_options() const { return false; } + +bool CanonicalizeStep::requires_options() const { return false; } + +void CanonicalizeStep::set_options(const nlohmann::json &) { + throw Exception(kind() + " doesn't take any options"); +} + +std::size_t CanonicalizeStep::process(std::string_view id_prefix, + std::size_t id_start, + ExecutionContext &ctx, + const ExpressionData &data) { + std::vector outputs; + + for (const ResultExpr &expr : data.expressions) { + ResultExpr clone = expr.clone(); + canonicalize(clone); + + outputs.emplace_back(std::move(clone)); + } + + if (outputs != data.expressions) { + ExpressionData data_obj; + data_obj.expressions.insert(data_obj.expressions.end(), + std::make_move_iterator(outputs.begin()), + std::make_move_iterator(outputs.end())); + ctx.set_data(id_prefix, id_start, std::move(data_obj)); + + return 1; + } + + return 0; +} + +bool CanonicalizeStep::alias_unchanged_inputs() const { return true; } + +} // namespace sequant::util::extint diff --git a/utilities/external-interface/canonicalize_step.hpp b/utilities/external-interface/canonicalize_step.hpp new file mode 100644 index 0000000000..8495ae51f5 --- /dev/null +++ b/utilities/external-interface/canonicalize_step.hpp @@ -0,0 +1,35 @@ +#ifndef SEQUANT_EXTERNAL_INTERFACE_CANONICALIZESTEP_HPP +#define SEQUANT_EXTERNAL_INTERFACE_CANONICALIZESTEP_HPP + +#include "execution_context.hpp" +#include "processing_data.hpp" +#include "processing_step.hpp" + +#include + +#include + +#include +#include + +namespace sequant::util::extint { + +class CanonicalizeStep : public OneByOneProcessingStep { + public: + std::string kind() const override; + + bool accepts_options() const override; + bool requires_options() const override; + void set_options(const nlohmann::json &options) override; + + protected: + std::size_t process(std::string_view id_prefix, std::size_t id_start, + ExecutionContext &ctx, + const ExpressionData &data) override; + + bool alias_unchanged_inputs() const override; +}; + +} // namespace sequant::util::extint + +#endif // SEQUANT_EXTERNAL_INTERFACE_CANONICALIZESTEP_HPP diff --git a/utilities/external-interface/density_fitting_step.cpp b/utilities/external-interface/density_fitting_step.cpp new file mode 100644 index 0000000000..704dca0f92 --- /dev/null +++ b/utilities/external-interface/density_fitting_step.cpp @@ -0,0 +1,93 @@ +#include "density_fitting_step.hpp" +#include "processing_data.hpp" +#include "processing_step_factory.hpp" + +#include +#include +#include +#include +#include + +#include + +namespace sequant::util::extint { + +SEQUANT_EXTINT_REGISTER_STEP_TYPE(DensityFittingStep, "density_fitting"); + +std::string DensityFittingStep::kind() const { return "density_fitting"; } + +bool DensityFittingStep::accepts_options() const { return true; } + +bool DensityFittingStep::requires_options() const { return true; } + +void DensityFittingStep::set_options(const nlohmann::json &options) { + if (!options.is_object()) { + throw Exception(kind() + " expects a JSON object for its options!"); + } + + for (const auto &[key, value] : options.items()) { + if (key == "auxiliary_space") { + if (!value.is_string()) { + throw Exception("Option '" + key + "' for " + kind() + + " requires string argument"); + } + + aux_space_ = get_default_context().index_space_registry()->retrieve( + value.get()); + } else if (key == "integral_label") { + if (!value.is_string()) { + throw Exception("Option '" + key + "' for " + kind() + + " requires string argument"); + } + + two_elec_int_label_ = value.get(); + } else if (key == "df_tensor_label") { + if (!value.is_string()) { + throw Exception("Option '" + key + "' for " + kind() + + " requires string argument"); + } + + df_label_ = value.get(); + } else { + throw Exception("Unknown option key for " + kind() + ": '" + key + "'"); + } + } + + if (aux_space_ == IndexSpace::null) { + throw Exception(kind() + " requires the auxiliary_space option to be set"); + } +} + +std::size_t DensityFittingStep::process(std::string_view id_prefix, + std::size_t id_start, + ExecutionContext &ctx, + const ExpressionData &data) { + bool had_effect = false; + std::vector modified; + for (const ResultExpr &expr : data.expressions) { + ExprPtr result = + mbpt::density_fit(expr.expression(), aux_space_, + toUtf16(two_elec_int_label_), toUtf16(df_label_)); + + had_effect |= result != expr.expression(); + if (expr.produces_tensor()) { + modified.emplace_back( + ResultExpr(expr.result_as_tensor(), std::move(result))); + } else { + modified.emplace_back( + ResultExpr(expr.result_as_variable(), std::move(result))); + } + } + + if (had_effect) { + ctx.set_data(id_prefix, id_start, + ExpressionData{.expressions = std::move(modified)}); + return 1; + } + + return 0; +} + +bool DensityFittingStep::alias_unchanged_inputs() const { return true; } + +} // namespace sequant::util::extint diff --git a/utilities/external-interface/density_fitting_step.hpp b/utilities/external-interface/density_fitting_step.hpp new file mode 100644 index 0000000000..bac88b0250 --- /dev/null +++ b/utilities/external-interface/density_fitting_step.hpp @@ -0,0 +1,41 @@ +#ifndef SEQUANT_EXTERNAL_INTERFACE_DENSITYFITTINGSTEP_HPP +#define SEQUANT_EXTERNAL_INTERFACE_DENSITYFITTINGSTEP_HPP + +#include "execution_context.hpp" +#include "processing_data.hpp" +#include "processing_step.hpp" + +#include + +#include + +#include +#include +#include + +namespace sequant::util::extint { + +class DensityFittingStep : public OneByOneProcessingStep { + public: + std::string kind() const override; + + bool accepts_options() const override; + bool requires_options() const override; + void set_options(const nlohmann::json &options) override; + + protected: + std::size_t process(std::string_view id_prefix, std::size_t id_start, + ExecutionContext &ctx, + const ExpressionData &data) override; + + bool alias_unchanged_inputs() const override; + + private: + IndexSpace aux_space_ = IndexSpace::null; + std::string two_elec_int_label_ = "g"; + std::string df_label_ = "DF"; +}; + +} // namespace sequant::util::extint + +#endif // SEQUANT_EXTERNAL_INTERFACE_DENSITYFITTINGSTEP_HPP diff --git a/utilities/external-interface/examples/nevpt2_new.json b/utilities/external-interface/examples/nevpt2_new.json new file mode 100644 index 0000000000..dd89ab8151 --- /dev/null +++ b/utilities/external-interface/examples/nevpt2_new.json @@ -0,0 +1,115 @@ +{ + "driver_format_version": 2, + "index_spaces": [ + { + "label": "a", + "size": 1000, + "real_valued": true, + "meta": { + "name": "External", + "tag": "e" + } + }, + { + "label": "u", + "size": 5, + "real_valued": true, + "meta": { + "name": "Active", + "tag": "a" + } + }, + { + "label": "i", + "size": 80, + "real_valued": true, + "meta": { + "name": "Closed", + "tag": "c" + } + }, + { + "label": "F", + "size": 100, + "real_valued": true, + "meta": { + "name": "Auxiliary", + "tag": "F" + } + } + ], + "steps": [ + { + "id": "input", + "kind": "read_input", + "options": { + "file_path": [ + "nevpt2/nevpt2_en0.inp", + "nevpt2/nevpt2_en.inp", + + "nevpt2/nevpt2_res1_i1.inp", + "nevpt2/nevpt2_res1_s0.inp", + "nevpt2/nevpt2_res2_s1_singles.inp", + "nevpt2/nevpt2_res1_s1.inp", + + "nevpt2/nevpt2_res2_p0.inp", + "nevpt2/nevpt2_res2_p2.inp", + "nevpt2/nevpt2_res2_i2.inp", + "nevpt2/nevpt2_res2_p1.inp", + "nevpt2/nevpt2_res2_s1.inp", + "nevpt2/nevpt2_res2_s2.inp" + ], + "default_symmetry": "antisymmetric" + }, + "outputs": { + "ecc0": "0", + "ecc": "1", + "en": "0-1", + "res1": "2-5", + "res2": "6-11", + "res": "2-11" + } + }, + { + "kind": "validate", + "inputs": "input" + }, + { + "id": "DF", + "kind": "density_fitting", + "inputs": "input", + "options": { + "auxiliary_space": "F" + } + }, + { + "id": "traced", + "kind": "spintracing", + "inputs": "DF", + "options": { + "algorithm": "closed_shell" + } + }, + { + "id": "biorth", + "kind": "project", + "inputs": "traced.res", + "options": { + "method": "biorthogonal" + } + }, + { + "id": "opt", + "kind": "optimize", + "inputs": [ + "traced.en", + "biorth" + ] + }, + { + "id": "treeify", + "kind": "to_export_tree", + "inputs": "opt" + } + ] +} diff --git a/utilities/external-interface/execution_context.cpp b/utilities/external-interface/execution_context.cpp new file mode 100644 index 0000000000..849cd53b57 --- /dev/null +++ b/utilities/external-interface/execution_context.cpp @@ -0,0 +1,352 @@ +#include "execution_context.hpp" +#include "processing_data.hpp" + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace sequant::util::extint { + +void ExecutionContext::set_data(std::string_view prefix, std::size_t counter, + ProcessingData data) { + set_data(prefix, std::ranges::single_view{counter}, + std::ranges::single_view{std::move(data)}); +} + +void ExecutionContext::add_data_alias(std::string_view id, std::string alias) { + add_data_alias(std::ranges::single_view(std::move(id)), std::move(alias)); +} + +bool ExecutionContext::has_data(std::string_view id) const { + return id_to_data_indices_.find(id) != id_to_data_indices_.end(); +} + +std::size_t ExecutionContext::dataset_size(std::string_view id) const { + if (!is_valid_id(id, false)) { + throw Exception("Invalid ID '" + std::string(id) + "'"); + } + + auto it = id_to_data_indices_.find(id); + + return it == id_to_data_indices_.end() ? 0 : it->second.size(); +} + +template +std::vector< + ExecutionContext::Data>> +get_data_impl(DataVec &&data, ID2DataIdxMap &&id2dataidx, + DataIdx2IDMap &&dataidx2id, std::string_view id) { + if (!ExecutionContext::is_valid_id(id, true)) { + throw Exception("Invalid id '" + std::string(id) + "'"); + } + + using RefT = meta::mimic_constness_t; + + std::vector> selected_data; + + for (const std::string_view current : ExecutionContext::expand_id(id)) { + auto it = id2dataidx.find(current); + + if (it == id2dataidx.end()) { + throw Exception("No data available for ID '" + std::string(current) + + "'"); + } + + for (std::size_t idx : it->second) { + ExecutionContext::Data ret_data{.data = data.at(idx)}; + + for (std::string_view alias : dataidx2id.at(idx)) { +#if defined(__cpp_lib_associative_heterogeneous_insertion) && \ + __cpp_lib_associative_heterogeneous_insertion >= 202311L + if (id2dataidx.at(alias).size() > 1) { +#else + // No heterogenous lookup overload for at() until C++26 + if (id2dataidx.at(std::string(alias)).size() > 1) { +#endif + ret_data.associated_group_ids.emplace_back(std::move(alias)); + } else { + ret_data.associated_ids.emplace_back(std::move(alias)); + } + } + + // Ensure that we order the IDs such that the generic index-IDs come after + // potentially manually assigned (that will have more meaningful names) + auto not_ends_with_num = [](std::string_view id) -> bool { + SEQUANT_ASSERT(!id.empty()); + if (id.empty()) { + return false; + } + + auto it = id.rfind('.'); + if (it == std::string_view::npos) { + // The generic number-IDs always contain a period + return true; + } + + std::string_view suffix = id.substr(it + 1); + return suffix.find_first_not_of("0123456789") != std::string_view::npos; + }; + + std::ranges::stable_partition(ret_data.associated_ids, not_ends_with_num); + std::ranges::stable_partition(ret_data.associated_group_ids, + not_ends_with_num); + + selected_data.emplace_back(std::move(ret_data)); + } + } + + return selected_data; +} + +std::vector> +ExecutionContext::get_data(std::string_view id) const { + return get_data_impl(data_, id_to_data_indices_, data_idx_to_ids_, id); +} + +std::vector> ExecutionContext::get_data( + std::string_view id) { + return get_data_impl(data_, id_to_data_indices_, data_idx_to_ids_, id); +} + +bool ExecutionContext::is_valid_id(std::string_view id, bool allow_selectors) { + auto validate_non_selector = [](std::string_view part) -> bool { + for (char c : part) { + if (!std::isalnum(c) & c != '_' && c != '.') { + return false; + } + } + + return true; + }; + + auto validate_range_component = [](std::string_view comp) -> bool { + for (char c : comp) { + if (!std::isdigit(c)) { + return false; + } + } + + return true; + }; + + auto validate_selector = [&validate_non_selector, &validate_range_component]( + std::string_view selector) -> bool { + for (auto &&comp : selector | std::ranges::views::split(',')) { + std::string_view part(comp.begin(), comp.end()); + + if (validate_non_selector(part)) { + continue; + } + + auto dash_pos = part.find('-'); + if (dash_pos == std::string_view::npos || + part.find('-', dash_pos + 1) != std::string_view::npos) { + // Either no dash or more than one dash in single selector component + return false; + } + + std::string_view prefix = part.substr(0, dash_pos); + std::string_view suffix = part.substr(dash_pos + 1); + + if (prefix.empty() || suffix.empty()) { + return false; + } + + if (!validate_range_component(prefix) || + !validate_range_component(suffix)) { + return false; + } + } + + return true; + }; + + if (id.empty()) { + return false; + } + + std::string_view::size_type bracket_begin = 0; + std::string_view::size_type prev_pos = 0; + do { + bracket_begin = id.find('[', prev_pos); + + if (!validate_non_selector(id.substr(prev_pos, bracket_begin - prev_pos))) { + return false; + } + + if (bracket_begin == std::string_view::npos) { + continue; + } + + if (!allow_selectors) { + return false; + } + + std::string_view::size_type bracket_end = id.find(']', bracket_begin); + + if (bracket_end == std::string_view::npos) { + return false; + } + if (bracket_begin + 1 == bracket_end) { + return false; + } + + if (!validate_selector( + id.substr(bracket_begin + 1, bracket_end - bracket_begin - 1))) { + return false; + } + + prev_pos = bracket_end + 1; + } while (bracket_begin != std::string_view::npos && prev_pos < id.size()); + + return true; +} + +std::vector expand_selector(std::string_view selector) { + std::vector processed; + + for (auto &&comp : selector | std::ranges::views::split(',')) { + std::string_view current(comp.begin(), comp.end()); + // trim whitespace + while (!current.empty() && current.front() == ' ') { + current.remove_prefix(1); + } + while (!current.empty() && current.back() == ' ') { + current.remove_suffix(1); + } + + if (current.empty()) { + throw Exception("Empty selector component in '[" + std::string(selector) + + "']"); + } + + if (auto dash_pos = current.find('-'); dash_pos != std::string::npos) { + // Numeric ranges such as "1-3" + std::size_t from = string_to(current.substr(0, dash_pos)); + std::size_t to = string_to(current.substr(dash_pos + 1)); + + if (from > to) { + std::swap(from, to); + } + + for (std::size_t val : std::ranges::views::iota(from, to + 1)) { + processed.emplace_back(std::to_string(val)); + } + } else { + processed.emplace_back(std::move(current)); + } + } + + SEQUANT_ASSERT(std::none_of(processed.begin(), processed.end(), + [](const auto &p) { return p.empty(); })); + + return processed; +} + +std::vector> create_id_partitions( + std::string_view id) { + std::vector> partitions; + + std::size_t begin = -1; + std::size_t prev_pos = 0; + do { + begin = id.find('[', begin + 1); + + std::vector part = {std::string(id.substr(prev_pos, begin))}; + SEQUANT_ASSERT(!part.back().empty()); + partitions.emplace_back(std::move(part)); + + if (begin == std::string_view::npos) { + continue; + } + + auto end = id.find(']', begin); + + if (end == std::string_view::npos) { + throw Exception("Unbalanced brackets in selector '" + std::string(id) + + "'"); + } + if (begin + 1 == end) { + throw Exception("Empty selector in '" + std::string(id) + "'"); + } + + std::string_view selector = id.substr(begin + 1, end - begin - 1); + partitions.emplace_back(expand_selector(selector)); + + prev_pos = end + 1; + } while (begin != std::string_view::npos && prev_pos < id.size()); + + SEQUANT_ASSERT(std::none_of(partitions.begin(), partitions.end(), + [](const auto &p) { return p.empty(); })); + + return partitions; +} + +std::vector ExecutionContext::expand_id(std::string_view id) { + if (!is_valid_id(id, true)) { + throw Exception("Invalid id '" + std::string(id) + "'"); + } + + std::vector> partitions = create_id_partitions(id); + std::vector indices(partitions.size(), 0); + + auto has_more = [&partitions, &indices]() { + for (std::size_t i = 0; i < indices.size(); ++i) { + SEQUANT_ASSERT(partitions[i].size() > 0); + if (indices[i] >= partitions[i].size()) { + return false; + } + } + + return true; + }; + + auto increment = [&partitions, &indices]() mutable { + for (std::size_t i = 0; i < indices.size(); ++i) { + if (indices[i] < partitions[i].size() - 1) { + ++indices[i]; + return; + } + indices[i] = 0; + } + + // indicate end has been reached + std::fill(indices.begin(), indices.end(), + std::numeric_limits::max()); + }; + + SEQUANT_ASSERT(has_more()); + + std::vector expanded; + + while (has_more()) { + std::stringstream stream; + for (std::size_t i = 0; i < indices.size(); ++i) { + stream << partitions[i].at(indices[i]); + + if (i + 1 < indices.size()) { + stream << "."; + } + } + + expanded.emplace_back(stream.str()); + + increment(); + } + + SEQUANT_ASSERT(!expanded.empty()); + + return expanded; +} + +} // namespace sequant::util::extint diff --git a/utilities/external-interface/execution_context.hpp b/utilities/external-interface/execution_context.hpp new file mode 100644 index 0000000000..46781050a8 --- /dev/null +++ b/utilities/external-interface/execution_context.hpp @@ -0,0 +1,139 @@ +#ifndef SEQUANT_EXTERNAL_INTERFACE_EXECUTIONCONTEXT_HPP +#define SEQUANT_EXTERNAL_INTERFACE_EXECUTIONCONTEXT_HPP + +#include "processing_data.hpp" + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace sequant::util::extint { + +class ExecutionContext { + public: + template + struct Data { + /// The actual data object (reference) + std::reference_wrapper data; + /// List of IDs that this data object is assigned to + std::vector associated_ids = {}; + /// List of IDs that are associated with a group of data objects which the + /// present one is a part of + std::vector associated_group_ids = {}; + + operator DataT &() { return data.get(); } + operator std::add_const_t &() const { return data.get(); } + }; + + ExecutionContext() = default; + + void set_data(std::string_view prefix, std::size_t counter, + ProcessingData data); + + template + requires(std::integral> && + std::same_as, ProcessingData>) + void set_data(std::string_view prefix, CounterValues &&counters, + Data &&data) { + if (std::ranges::empty(counters)) { + throw Exception("Attempted to add data without specifying any counter"); + } + if (std::ranges::empty(data)) { + throw Exception("Attempted to register empty dataset"); + } + + std::vector ids; + ids.reserve(std::ranges::size(counters)); + + for (auto &¤t : counters) { + ids.emplace_back(std::string(prefix) + "." + std::to_string(current)); + + if (id_to_data_indices_.find(ids.back()) != id_to_data_indices_.end()) { + throw Exception("Duplicate data ID '" + ids.back() + "'"); + } + if (!is_valid_id(ids.back(), false)) { + throw Exception("Invalid id '" + ids.back() + "'"); + } + } + + std::vector data_indices(std::ranges::size(data)); + std::iota(data_indices.begin(), data_indices.end(), data_.size()); + + for (std::string ¤t : ids) { + set_id_idx_assoc(std::move(current), data_indices); + } + + for (auto &¤t : data) { + data_.emplace_back(std::move(current)); + } + } + + void add_data_alias(std::string_view id, std::string alias); + + template + requires(std::is_convertible_v, + std::string_view>) + void add_data_alias(IDs &&ids, std::string alias) { + if (id_to_data_indices_.find(alias) != id_to_data_indices_.end()) { + throw Exception("Alias '" + alias + "' already exists as a data ID"); + } + + for (const auto ¤t_id : ids) { + for (const std::string &expanded : expand_id(current_id)) { + auto it = id_to_data_indices_.find(expanded); + + if (it == id_to_data_indices_.end()) { + throw Exception("Attempted to alias non-existent ID '" + expanded + + "'"); + } + + SEQUANT_ASSERT(!it->second.empty()); + set_id_idx_assoc(alias, it->second); + } + } + } + + bool has_data(std::string_view id) const; + + std::size_t dataset_size(std::string_view id) const; + + std::vector> get_data(std::string_view id) const; + std::vector> get_data(std::string_view id); + + static bool is_valid_id(std::string_view id, bool allow_selectors); + + static std::vector expand_id(std::string_view id); + + private: + std::deque data_; + std::map, std::less<>> + id_to_data_indices_; + std::map> data_idx_to_ids_; + + template + void set_id_idx_assoc(std::string id, Indices &&indices) { + if constexpr (!std::ranges::range) { + set_id_idx_assoc(std::move(id), + std::ranges::single_view(std::move(indices))); + } else { + for (std::size_t idx : indices) { + id_to_data_indices_[id].emplace_back(idx); + data_idx_to_ids_[idx].emplace_back(id); + } + } + } +}; + +} // namespace sequant::util::extint + +#endif // SEQUANT_EXTERNAL_INTERFACE_EXECUTIONCONTEXT_HPP diff --git a/utilities/external-interface/executor.cpp b/utilities/external-interface/executor.cpp new file mode 100644 index 0000000000..7d8eeff757 --- /dev/null +++ b/utilities/external-interface/executor.cpp @@ -0,0 +1,112 @@ +#include "executor.hpp" +#include "processing_step.hpp" +#include "processing_step_factory.hpp" + +#include + +#include +#include +#include +#include + +namespace sequant::util::extint { + +void Executor::execute(const nlohmann::json &steps) { + if (!steps.is_array()) { + throw Exception("Steps object must be an array"); + } + + std::size_t step_id_counter = 0; + + for (const nlohmann::json &step : steps) { + const std::string_view kind = step.at("kind").get(); + + std::vector inputs; + + if (step.contains("inputs")) { + const nlohmann::json &inps = step.at("inputs"); + if (inps.is_string()) { + inputs.emplace_back(inps.get()); + } else if (inps.is_array()) { + for (const auto ¤t : inps) { + if (!current.is_string()) { + throw Exception("Entries in inputs array must be strings"); + } + + inputs.emplace_back(current.get()); + } + } else { + throw Exception("inputs field must be either a string or an array"); + } + } + + std::unique_ptr proc_step = + ProcessingStepFactory::instance().instantiate(kind); + + if (step.contains("options")) { + if (!proc_step->accepts_options()) { + throw Exception("Processing step '" + std::string(kind) + + "' does not take options but some where given"); + } + + proc_step->set_options(step.at("options")); + } else if (proc_step->requires_options()) { + throw Exception("Processing step '" + std::string(kind) + + "' requires options but none where given"); + } + + std::string step_id = step.contains("id") + ? step.at("id").get() + : "step" + std::to_string(step_id_counter) + "." + + proc_step->kind(); + + if (num_outputs_.find(step_id) != num_outputs_.end()) { + throw Exception("Duplicate step ID '" + step_id + "'"); + } + + std::cout << "Executing '" << step_id << "' (" << kind << ")... "; + std::cout.flush(); + + std::chrono::steady_clock::time_point start = + std::chrono::steady_clock::now(); + + const std::size_t produced_outputs = + proc_step->run(step_id, context_, inputs); + + std::chrono::steady_clock::duration delta = + std::chrono::steady_clock::now() - start; + if (delta > std::chrono::minutes(1)) { + std::cout << std::chrono::duration_cast(delta) + << std::endl; + } else if (delta > std::chrono::seconds(1)) { + std::cout << std::chrono::duration_cast(delta) + << std::endl; + } else { + std::cout << std::chrono::duration_cast(delta) + << std::endl; + } + + num_outputs_.emplace(step_id, produced_outputs); + + if (step.contains("outputs")) { + const nlohmann::json &outputs = step.at("outputs"); + if (!outputs.is_object()) { + throw Exception("outputs field must be an object"); + } + + for (const auto &[key, val] : outputs.items()) { + context_.add_data_alias(step_id + "[" + val.get() + "]", + step_id + "." + key); + } + } + + ++step_id_counter; + } +} + +void Executor::reset() { + num_outputs_.clear(); + context_ = {}; +} + +} // namespace sequant::util::extint diff --git a/utilities/external-interface/executor.hpp b/utilities/external-interface/executor.hpp new file mode 100644 index 0000000000..0f5405649d --- /dev/null +++ b/utilities/external-interface/executor.hpp @@ -0,0 +1,30 @@ +#ifndef SEQUANT_EXTERNAL_INTERFACE_EXECUTOR_HPP +#define SEQUANT_EXTERNAL_INTERFACE_EXECUTOR_HPP + +#include "execution_context.hpp" + +#include + +#include +#include +#include +#include + +namespace sequant::util::extint { + +class Executor { + public: + Executor() = default; + + void execute(const nlohmann::json &steps); + + void reset(); + + private: + ExecutionContext context_; + std::map> num_outputs_; +}; + +} // namespace sequant::util::extint + +#endif // SEQUANT_EXTERNAL_INTERFACE_EXECUTOR_HPP diff --git a/utilities/external-interface/export_step.cpp b/utilities/external-interface/export_step.cpp new file mode 100644 index 0000000000..65aa3918c4 --- /dev/null +++ b/utilities/external-interface/export_step.cpp @@ -0,0 +1,223 @@ +#include "export_step.hpp" +#include "processing_data.hpp" +#include "processing_step_factory.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace sequant::util::extint { + +SEQUANT_EXTINT_REGISTER_STEP_TYPE(ExportStep, "export"); + +std::string ExportStep::kind() const { return "export"; } + +bool ExportStep::accepts_options() const { return true; } + +bool ExportStep::requires_options() const { return true; } + +void ExportStep::set_options(const nlohmann::json &options) { + if (!options.is_object()) { + throw Exception(kind() + " expects a JSON object for its options!"); + } + + for (const auto &[key, value] : options.items()) { + if (key == "language") { + if (!value.is_string()) { + throw Exception("Value for " + kind() + " option '" + key + + "' must be a string"); + } + + language_ = value.get(); + } else if (key == "optimize") { + if (!value.is_boolean()) { + throw Exception("Value for " + kind() + " option '" + key + + "' must be a boolean"); + } + + optimize_ = value.get(); + } else if (key == "output") { + if (!value.is_string()) { + throw Exception("Value for " + kind() + " option '" + key + + "' must be a string"); + } + + filepath_ = value.get(); + } else { + throw Exception("Unknown option key for " + kind() + ": '" + key + "'"); + } + + if (language_.empty()) { + throw Exception("The 'language' option for " + kind() + " is mandatory"); + } + + if (filepath_.empty()) { + throw Exception("The 'output' option for " + kind() + " is mandatory"); + } + } +} + +struct ExportTreeDataCompare { + bool operator()(const ExportTreeData &lhs, const ExportTreeData &rhs) const { + if (lhs.entries.size() != rhs.entries.size()) { + return lhs.entries.size() < rhs.entries.size(); + } + + for (std::size_t i = 0; i < lhs.entries.size(); ++i) { + if (!equal(lhs.entries.at(i), rhs.entries.at(i))) { + return (*this)(lhs.entries.at(i), rhs.entries.at(i)); + } + } + + return false; + } + + bool equal(const ExportTreeData::Entry &lhs, + const ExportTreeData::Entry &rhs) const { + if (lhs.symm_contribution_target != rhs.symm_contribution_target) { + return false; + } + + if (lhs.tree->is_tensor() != rhs.tree->is_tensor()) { + return false; + } + + return lhs.tree->is_tensor() + ? lhs.tree->as_tensor() == rhs.tree->as_tensor() + : lhs.tree->as_variable() == rhs.tree->as_variable(); + } + + bool operator()(const ExportTreeData::Entry &lhs, + const ExportTreeData::Entry &rhs) const { + if (lhs.symm_contribution_target != rhs.symm_contribution_target) { + return lhs.symm_contribution_target < rhs.symm_contribution_target; + } + + if (lhs.tree->is_tensor() != rhs.tree->is_tensor()) { + return rhs.tree->is_tensor(); + } + + return lhs.tree->is_tensor() + ? lhs.tree->as_tensor() < rhs.tree->as_tensor() + : lhs.tree->as_variable() < rhs.tree->as_variable(); + } +}; + +std::size_t ExportStep::run(std::string_view, ExecutionContext &exctx, + const std::vector &inputs) { + std::vector< + std::reference_wrapper>> + data; + for (std::string_view current_input : inputs) { + for (const ExecutionContext::Data ¤t : + exctx.get_data(current_input)) { + data.push_back(std::cref(current)); + } + } + + auto to_export_data = [&](std::size_t idx) { + return convert_data(data.at(idx).get().data.get()); + }; + + std::vector access_order(data.size()); + std::iota(access_order.begin(), access_order.end(), 0); + + // Sort inputs in partitions of contributions to the same result + std::ranges::sort(access_order, ExportTreeDataCompare{}, to_export_data); + + std::vector> groups; + for (std::size_t idx : access_order) { + const ExportTreeData ¤t_data = to_export_data(idx); + + if (std::ranges::any_of(current_data.entries, [&](const auto &entry) { + return !ExportTreeDataCompare{}.equal(entry, + current_data.entries.front()); + })) { + throw Exception( + "Expected all entries in an ExportTreeData object to contribute to " + "the same result"); + } + + // TODO: determine group name from metadata + std::string name = "We'll see"; + + auto it = + std::ranges::find_if(groups, [&name](const ExpressionGroup<> &group) { + return group.is_named() && group.name() == name; + }); + if (it == groups.end()) { + groups.emplace_back(std::move(name)); + it = groups.end() - 1; + } + + ExpressionGroup<> &group = *it; + + for (const ExportTreeData::Entry ¤t : current_data.entries) { + group.add(current.tree); + } + + // TODO: if current_data was last in its partition (based on above sorting) + // we need to check if symmetrization is required and if so, generate the + // code snippet that does it and add it to the group + } + + auto perform_export = [&](auto &&generator, const auto &genctx) + requires( + std::derived_from, + ExportContext> && + std::derived_from, + Generator>>) + { + // TODO: setup load/store/create/import strategies + + // TODO: handle index batching + + export_groups(groups, generator, genctx); + + return generator.get_generated_code(); + }; + + std::string generated; + + if (language_ == "itf") { + ItfGenerator generator; + ItfContext genctx; + + if (optimize_) { + generated = perform_export( + GenerationOptimizer(std::move(generator)), + genctx); + } else { + generated = perform_export(generator, genctx); + } + } else { + throw Exception("Unsupported export language '" + language_ + "'"); + } + + SEQUANT_ASSERT(!generated.empty()); + + std::ofstream stream(filepath_); + stream << generated; + + return 0; +} + +} // namespace sequant::util::extint diff --git a/utilities/external-interface/export_step.hpp b/utilities/external-interface/export_step.hpp new file mode 100644 index 0000000000..efff06825f --- /dev/null +++ b/utilities/external-interface/export_step.hpp @@ -0,0 +1,35 @@ +#ifndef SEQUANT_EXTERNAL_INTERFACE_EXPORTSTEP_HPP +#define SEQUANT_EXTERNAL_INTERFACE_EXPORTSTEP_HPP + +#include "execution_context.hpp" +#include "processing_data.hpp" +#include "processing_step.hpp" + +#include + +#include +#include +#include + +namespace sequant::util::extint { + +class ExportStep : public ProcessingStep { + public: + std::string kind() const override; + + bool accepts_options() const override; + bool requires_options() const override; + void set_options(const nlohmann::json &options) override; + + std::size_t run(std::string_view step_id, ExecutionContext &ctx, + const std::vector &inputs = {}) override; + + protected: + std::string language_; + bool optimize_ = true; + std::filesystem::path filepath_; +}; + +} // namespace sequant::util::extint + +#endif // SEQUANT_EXTERNAL_INTERFACE_EXPORTSTEP_HPP diff --git a/utilities/external-interface/external_interface.cpp b/utilities/external-interface/external_interface.cpp index c7a8a622d6..c026d59769 100644 --- a/utilities/external-interface/external_interface.cpp +++ b/utilities/external-interface/external_interface.cpp @@ -1,3 +1,4 @@ +#include "executor.hpp" #include "format_support.hpp" #include "processing.hpp" #include "utils.hpp" @@ -36,22 +37,25 @@ #include #include #include +#include #include #include #include #include #include +#include using nlohmann::json; -using namespace sequant; template <> -struct std::hash { - std::size_t operator()(const Tensor &tensor) const { +struct std::hash { + std::size_t operator()(const sequant::Tensor &tensor) const { return tensor.hash_value(); } }; +namespace sequant::util::extint { + class ItfExportContext : public ItfContext { public: ItfExportContext(const IndexSpaceMeta &meta) : m_meta(&meta) {} @@ -604,7 +608,8 @@ void generateCode(const json &details, const IndexSpaceMeta &spaceMeta) { } } -void registerIndexSpaces(const json &spaces, IndexSpaceMeta &meta) { +void registerIndexSpaces(const json &spaces, IndexSpaceMeta &meta, + std::size_t version) { IndexSpaceRegistry ®istry = *get_default_context().mutable_index_space_registry(); @@ -621,8 +626,13 @@ void registerIndexSpaces(const json &spaces, IndexSpaceMeta &meta) { } IndexSpaceMeta::Entry entry; - entry.name = current.at("name").get(); - entry.tag = current.at("tag").get(); + if (version == 1) { + entry.name = current.at("name").get(); + entry.tag = current.at("tag").get(); + } else { + entry.name = current.at("meta").at("name").get(); + entry.tag = current.at("meta").at("tag").get(); + } std::wstring label = toUtf16(current.at("label").get()); Field field = @@ -649,12 +659,34 @@ void process(const json &driver, IndexSpaceMeta &spaceMeta) { throw Exception("Missing index_spaces definition"); } - registerIndexSpaces(driver.at("index_spaces"), spaceMeta); + std::size_t version = 1; + if (driver.contains("driver_format_version")) { + version = driver.at("driver_format_version").get(); + } + + if (version == 0) { + throw Exception("driver_format_version has a minimum value of 1"); + } + + registerIndexSpaces(driver.at("index_spaces"), spaceMeta, version); + + if (version == 1) { + if (driver.contains("code_generation")) { + const json &details = driver.at("code_generation"); - if (driver.contains("code_generation")) { - const json &details = driver.at("code_generation"); + generateCode(details, spaceMeta); + } + } else if (version == 2) { + if (!driver.contains("steps")) { + throw Exception("Missing steps specification"); + } - generateCode(details, spaceMeta); + Executor executor; + executor.execute(driver.at("steps")); + } else { + throw Exception( + "Requested driver_format_version too recent for this implementation: " + + std::to_string(version)); } } @@ -663,7 +695,10 @@ void generalSetup() { mbpt::cardinal_tensor_labels()); } +} // namespace sequant::util::extint + int main(int argc, char **argv) { + using namespace sequant; set_locale(); Context ctx({.index_space_registry = IndexSpaceRegistry(), .vacuum = Vacuum::SingleProduct}); @@ -674,7 +709,7 @@ int main(int argc, char **argv) { // to use the new names. ctx.set(CanonicalizeOptions{.method = CanonicalizationMethod::Complete}); set_default_context(ctx); - generalSetup(); + util::extint::generalSetup(); CLI::App app( "Interface for reading in equations generated outside of SeQuant"); @@ -717,7 +752,7 @@ int main(int argc, char **argv) { json::parse(in, /*callback*/ nullptr, /*allow_exceptions*/ true, /*skip_comments*/ true); - process(driver_info, spaceMeta); + util::extint::process(driver_info, spaceMeta); } catch (const std::exception &e) { spdlog::error("Unexpected error: {}", e.what()); return 1; diff --git a/utilities/external-interface/optimization_step.cpp b/utilities/external-interface/optimization_step.cpp new file mode 100644 index 0000000000..9d9e13088f --- /dev/null +++ b/utilities/external-interface/optimization_step.cpp @@ -0,0 +1,117 @@ +#include "optimization_step.hpp" +#include "processing_data.hpp" +#include "processing_step_factory.hpp" +#include "utils.hpp" + +#include +#include + +#include + +#include +#include + +namespace sequant::util::extint { + +SEQUANT_EXTINT_REGISTER_STEP_TYPE(OptimizationStep, "optimize"); + +std::string OptimizationStep::kind() const { return "optimize"; } + +bool OptimizationStep::accepts_options() const { return true; } + +bool OptimizationStep::requires_options() const { return false; } + +void OptimizationStep::set_options(const nlohmann::json &options) { + if (!options.is_object()) { + throw Exception(kind() + " expects a JSON object for its options!"); + } + + for (const auto &[key, value] : options.items()) { + if (key == "objective") { + if (!value.is_string()) { + throw Exception("Value for " + kind() + " option '" + key + + "' must be a string"); + } + + if (value == "DenseFLOPs") { + options_.objective_function = ObjectiveFunction::DenseFLOPs; + } else if (value == "DenseSize") { + options_.objective_function = ObjectiveFunction::DenseSize; + } else if (value == "DensePeakSize") { + options_.objective_function = ObjectiveFunction::DensePeakSize; + } else if (value == "DensePeakSizeBatched") { + options_.objective_function = ObjectiveFunction::DensePeakSizeBatched; + } else { + throw Exception("Invalid value for " + kind() + " option '" + key + + "': '" + value.get() + "'"); + } + } else if (key == "reorder_sums") { + if (!value.is_boolean()) { + throw Exception("Value for " + kind() + " option '" + key + + "' must be a boolean"); + } + + options_.reorder = + value.get() ? ReorderSum::Reorder : ReorderSum::NoReorder; + } else if (key == "cse") { + if (!value.is_string()) { + throw Exception("Value for " + kind() + " option '" + key + + "' must be a string"); + } + + if (value == "none") { + options_.CSE.subnet = false; + } else if (value == "subnet") { + options_.CSE.subnet = true; + } else { + throw Exception("Invalid value for " + kind() + " option '" + key + + "': '" + value.get() + "'"); + } + } else if (key == "intermediate_size_penalty") { + if (!value.is_number()) { + throw Exception("Value for " + kind() + " option '" + key + + "' must be a number"); + } + + options_.footprint_weight = value.get(); + } else if (key == "prune_outer_products") { + if (!value.is_boolean()) { + throw Exception("Value for " + kind() + " option '" + key + + "' must be a boolean"); + } + + options_.prune_outer_products = value.get(); + } else { + throw Exception("Unknown option key for " + kind() + ": '" + key + "'"); + } + } +} + +std::size_t OptimizationStep::process(std::string_view id_prefix, + std::size_t id_start, + ExecutionContext &ctx, + const ExpressionData &data) { + ExpressionData result; + result.expressions.reserve(data.expressions.size()); + for (const ResultExpr &input : data.expressions) { + result.expressions.emplace_back(input.clone()); + + std::optional symmetrizer = + pop_symmetrizer(result.expressions.back()); + + optimize(result.expressions.back(), options_); + + if (symmetrizer.has_value()) { + result.expressions.back().expression() = ex( + ExprPtrList{std::move(symmetrizer.value()), + std::move(result.expressions.back().expression())}, + Product::Flatten::No); + } + } + + ctx.set_data(id_prefix, id_start, std::move(result)); + + return 1; +} + +} // namespace sequant::util::extint diff --git a/utilities/external-interface/optimization_step.hpp b/utilities/external-interface/optimization_step.hpp new file mode 100644 index 0000000000..e3476fdb48 --- /dev/null +++ b/utilities/external-interface/optimization_step.hpp @@ -0,0 +1,35 @@ +#ifndef SEQUANT_EXTERNAL_INTERFACE_OPTIMIZATIONSTEP_HPP +#define SEQUANT_EXTERNAL_INTERFACE_OPTIMIZATIONSTEP_HPP + +#include "execution_context.hpp" +#include "processing_data.hpp" +#include "processing_step.hpp" + +#include + +#include + +#include +#include + +namespace sequant::util::extint { + +class OptimizationStep : public OneByOneProcessingStep { + public: + std::string kind() const override; + + bool accepts_options() const override; + bool requires_options() const override; + void set_options(const nlohmann::json &options) override; + + protected: + OptimizeOptions options_; + + std::size_t process(std::string_view id_prefix, std::size_t id_start, + ExecutionContext &ctx, + const ExpressionData &data) override; +}; + +} // namespace sequant::util::extint + +#endif // SEQUANT_EXTERNAL_INTERFACE_OPTIMIZATIONSTEP_HPP diff --git a/utilities/external-interface/output_step.cpp b/utilities/external-interface/output_step.cpp new file mode 100644 index 0000000000..013c948894 --- /dev/null +++ b/utilities/external-interface/output_step.cpp @@ -0,0 +1,74 @@ +#include "output_step.hpp" +#include "processing_step_factory.hpp" + +#include +#include +#include +#include + +#include + +#include +#include + +namespace sequant::util::extint { + +SEQUANT_EXTINT_REGISTER_STEP_TYPE(OutputStep, "output"); + +std::string OutputStep::kind() const { return "Output"; } + +bool OutputStep::accepts_options() const { return true; } + +bool OutputStep::requires_options() const { return false; } + +void OutputStep::set_options(const nlohmann::json &options) { + if (!options.is_object()) { + throw Exception(kind() + " expects a JSON object for its options!"); + } + + for (const auto &[key, value] : options.items()) { + if (key == "format") { + if (!value.is_string()) { + throw Exception("Option 'format' for " + kind() + + " is expected to be a string"); + } + + if (value == "latex") { + latex_ = true; + } else if (value == "serialize") { + latex_ = false; + } else { + throw Exception("Invalid output format for " + kind() + ": '" + + value.get() + "'"); + } + } else if (key == "annotate_symmetry") { + if (!value.is_boolean()) { + throw Exception("Option 'annotate_symmetry' for " + kind() + + " is expected to be a " + "boolean"); + } + + annot_symm_ = value.get(); + } else { + throw Exception("Unknown option key for " + kind() + ": '" + key + "'"); + } + } +} + +std::size_t OutputStep::process(std::string_view, std::size_t, + ExecutionContext &, + const ExpressionData &data) { + for (const auto &expr : data.expressions) { + if (latex_) { + std::wcout << io::latex::to_string(expr) << std::endl; + } else { + std::wcout << io::serialization::to_string(expr, + {.annot_symm = annot_symm_}) + << std::endl; + } + } + + return 0; +} + +} // namespace sequant::util::extint diff --git a/utilities/external-interface/output_step.hpp b/utilities/external-interface/output_step.hpp new file mode 100644 index 0000000000..b9a5280791 --- /dev/null +++ b/utilities/external-interface/output_step.hpp @@ -0,0 +1,35 @@ +#ifndef SEQUANT_EXTERNAL_INTERFACE_OUTPUTSTEP_HPP +#define SEQUANT_EXTERNAL_INTERFACE_OUTPUTSTEP_HPP + +#include "execution_context.hpp" +#include "processing_data.hpp" +#include "processing_step.hpp" + +#include + +#include +#include + +namespace sequant::util::extint { + +class OutputStep : public OneByOneProcessingStep { + public: + std::string kind() const override; + + bool accepts_options() const override; + bool requires_options() const override; + void set_options(const nlohmann::json &options) override; + + protected: + std::size_t process(std::string_view id_prefix, std::size_t id_start, + ExecutionContext &ctx, + const ExpressionData &data) override; + + private: + bool latex_ = false; + bool annot_symm_ = true; +}; + +} // namespace sequant::util::extint + +#endif // SEQUANT_EXTERNAL_INTERFACE_OUTPUTSTEP_HPP diff --git a/utilities/external-interface/processing_data.hpp b/utilities/external-interface/processing_data.hpp new file mode 100644 index 0000000000..cc9150a36b --- /dev/null +++ b/utilities/external-interface/processing_data.hpp @@ -0,0 +1,59 @@ +#ifndef SEQUANT_EXTERNAL_INTERFACE_PROCESSINGDATA_HPP +#define SEQUANT_EXTERNAL_INTERFACE_PROCESSINGDATA_HPP + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace sequant::util::extint { + +struct ExpressionData { + static constexpr std::string_view name{"ExpressionData"}; + + std::vector expressions; +}; + +struct ExportTreeData { + static constexpr std::string_view name{"ExportTreeData"}; + + struct Entry { + ExportNode<> tree; + /// The result the current tree contributes to _after symmetrizing_ over the + /// external indices. This implies that if this is set, symmetrization is + /// required. + std::optional symm_contribution_target; + }; + + std::vector entries; +}; + +using ProcessingData = std::variant; + +template + requires(std::is_assignable_v, Target>) +decltype(auto) convert_data(Input &&input) { + using RetType = decltype(meta::forward_like(std::declval())); + + return std::visit( + [](auto &&inp) -> RetType { + using Current = std::remove_cvref_t; + if constexpr (std::is_same_v) { + return std::forward(inp); + } + + throw Exception("Can't convert " + std::string(Current::name) + " to " + + std::string(Target::name)); + }, + std::forward(input)); +} + +} // namespace sequant::util::extint + +#endif // SEQUANT_EXTERNAL_INTERFACE_PROCESSINGDATA_HPP diff --git a/utilities/external-interface/processing_step.cpp b/utilities/external-interface/processing_step.cpp new file mode 100644 index 0000000000..502cf561a7 --- /dev/null +++ b/utilities/external-interface/processing_step.cpp @@ -0,0 +1,7 @@ +#include "processing_step.hpp" + +namespace sequant::util::extint { + +ProcessingStep::~ProcessingStep() = default; + +} // namespace sequant::util::extint diff --git a/utilities/external-interface/processing_step.hpp b/utilities/external-interface/processing_step.hpp new file mode 100644 index 0000000000..5ed759f7a0 --- /dev/null +++ b/utilities/external-interface/processing_step.hpp @@ -0,0 +1,167 @@ +#ifndef SEQUANT_EXTERNAL_INTERFACE_PROCESSINGSTEP_HPP +#define SEQUANT_EXTERNAL_INTERFACE_PROCESSINGSTEP_HPP + +#include "execution_context.hpp" + +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +namespace sequant::util::extint { + +class ProcessingStep { + public: + ProcessingStep() = default; + virtual ~ProcessingStep(); + + virtual std::string kind() const = 0; + + virtual bool accepts_options() const = 0; + virtual bool requires_options() const = 0; + virtual void set_options(const nlohmann::json &options) = 0; + + virtual std::size_t run(std::string_view step_id, ExecutionContext &ctx, + const std::vector &inputs = {}) = 0; +}; + +namespace detail { + +template +class OBOProcessFunctionMember { + protected: + virtual std::size_t process(std::string_view id_prefix, std::size_t id_start, + ExecutionContext &ctx, const DataType &data) = 0; +}; + +template <> +class OBOProcessFunctionMember { + protected: + virtual std::size_t process(std::string_view id_prefix, std::size_t id_start, + ExecutionContext &ctx) = 0; +}; + +} // namespace detail + +template +class OneByOneProcessingStep + : public ProcessingStep, + public detail::OBOProcessFunctionMember { + public: + std::size_t run(std::string_view step_id, ExecutionContext &ctx, + const std::vector &inputs = {}) override { + std::size_t total_outputs = 0; + + if constexpr (std::same_as) { + if (!inputs.empty()) { + throw Exception(kind() + " doesn't take any inputs"); + } + + try { + total_outputs = this->process(step_id, total_outputs, ctx); + } catch (const std::exception &e) { + throw Exception("Error while processing " + kind() + ": " + e.what()); + } + } else { + std::map> group_assocs; + + for (std::string_view current_input : inputs) { + for (const ExecutionContext::Data ¤t : + ctx.get_data(current_input)) { + const DataType &data = convert_data(current.data.get()); + + SEQUANT_ASSERT(!current.associated_ids.empty()); + + std::size_t new_outputs; + try { + new_outputs = this->process(step_id, total_outputs, ctx, data); + } catch (const std::exception &e) { + throw Exception("Error in " + kind() + " on input " + + std::string(current.associated_ids.front()) + ": " + + e.what()); + } + + if (new_outputs == 0 && alias_unchanged_inputs()) { + new_outputs = 1; + ctx.add_data_alias( + current.associated_ids.front(), + std::string(step_id) + "." + std::to_string(total_outputs)); + } + + if (new_outputs > 0) { + std::string created_outputs = + std::string(step_id) + "[" + std::to_string(total_outputs) + + "-" + std::to_string(total_outputs + new_outputs - 1) + "]"; + + for (std::string_view assoc_id : current.associated_ids) { + if (auto pos = assoc_id.find('.'); + pos != std::string_view::npos) { + // Remove step ID of previous step + assoc_id = assoc_id.substr(pos); + } + + if (assoc_id.empty() || + assoc_id.substr(1).find_first_not_of("0123456789") == + std::string_view::npos) { + // Skip "trivial" (autogenerated) IDs + continue; + } + + ctx.add_data_alias(created_outputs, + std::string(step_id) + std::string(assoc_id)); + } + + for (std::string_view group_id : current.associated_group_ids) { + group_assocs[group_id].emplace_back(created_outputs); + } + } + + total_outputs += new_outputs; + } + } + + for (auto &[group_id, assoc_outputs] : group_assocs) { + const std::size_t expected_size = ctx.dataset_size(group_id); + SEQUANT_ASSERT(expected_size >= assoc_outputs.size()); + + if (expected_size != assoc_outputs.size()) { + continue; + } + + std::string out_group_id; + if (auto pos = group_id.find('.'); pos != std::string_view::npos) { + // Remove previous step ID + out_group_id = + std::string(step_id) + std::string(group_id.substr(pos)); + } else { + out_group_id = std::string(step_id) + "." + std::string(group_id); + } + + ctx.add_data_alias(assoc_outputs, std::move(out_group_id)); + } + } + + ctx.add_data_alias( + std::ranges::views::iota(std::size_t(0), total_outputs) | + std::ranges::views::transform([&step_id](std::size_t num) { + return std::string(step_id) + "." + std::to_string(num); + }), + std::string(step_id)); + + return total_outputs; + } + + protected: + virtual bool alias_unchanged_inputs() const { return false; } +}; + +} // namespace sequant::util::extint + +#endif // SEQUANT_EXTERNAL_INTERFACE_PROCESSINGSTEP_HPP diff --git a/utilities/external-interface/processing_step_factory.cpp b/utilities/external-interface/processing_step_factory.cpp new file mode 100644 index 0000000000..ddba35ad16 --- /dev/null +++ b/utilities/external-interface/processing_step_factory.cpp @@ -0,0 +1,43 @@ +#include "processing_step_factory.hpp" +#include "processing_step.hpp" + +#include + +#include +#include +#include +#include + +namespace sequant::util::extint { + +ProcessingStepFactory &ProcessingStepFactory::instance() { + static ProcessingStepFactory factory; + + return factory; +} + +std::unique_ptr ProcessingStepFactory::instantiate( + std::string_view type) const { + auto it = instatiators_.find(type); + + if (it == instatiators_.end()) { + throw Exception( + "Attempted to instantiated unknown processing step of kind '" + + std::string(type) + "'"); + } + + return it->second(); +} + +bool ProcessingStepFactory::register_class(std::string type, + InstatiateFunc func) { + auto [it, inserted] = instatiators_.emplace(std::move(type), std::move(func)); + + if (!inserted) { + throw Exception("Duplicate processing step type name '" + it->first + "'"); + } + + return inserted; +} + +} // namespace sequant::util::extint diff --git a/utilities/external-interface/processing_step_factory.hpp b/utilities/external-interface/processing_step_factory.hpp new file mode 100644 index 0000000000..8bb2a969df --- /dev/null +++ b/utilities/external-interface/processing_step_factory.hpp @@ -0,0 +1,44 @@ +#ifndef SEQUANT_EXTERNAL_INTERFACE_PROCESSING_STEP_FACTORY_HPP +#define SEQUANT_EXTERNAL_INTERFACE_PROCESSING_STEP_FACTORY_HPP + +#include "processing_step.hpp" + +#include + +#include +#include +#include +#include +#include + +namespace sequant::util::extint { + +class ProcessingStepFactory { + public: + using InstatiateFunc = std::function()>; + + static ProcessingStepFactory &instance(); + + std::unique_ptr instantiate(std::string_view type) const; + + bool register_class(std::string type, InstatiateFunc func); + + private: + std::map> instatiators_; + + ProcessingStepFactory() = default; +}; + +#define SEQUANT_EXTINT_REGISTER_STEP_TYPE(class_name, type_string) \ + namespace { \ + static const bool SEQUANT_CONCAT(SEQUANT_CONCAT(registered_step_type_, \ + class_name), \ + __LINE__) = \ + ProcessingStepFactory::instance().register_class(type_string, []() { \ + return std::make_unique(); \ + }); \ + } + +} // namespace sequant::util::extint + +#endif // SEQUANT_EXTERNAL_INTERFACE_PROCESSING_STEP_FACTORY_HPP diff --git a/utilities/external-interface/projection_step.cpp b/utilities/external-interface/projection_step.cpp new file mode 100644 index 0000000000..20bd7e893e --- /dev/null +++ b/utilities/external-interface/projection_step.cpp @@ -0,0 +1,85 @@ +#include "projection_step.hpp" +#include "processing_data.hpp" +#include "processing_step_factory.hpp" +#include "utils.hpp" + +#include +#include +#include + +#include + +#include +#include + +namespace sequant::util::extint { + +SEQUANT_EXTINT_REGISTER_STEP_TYPE(ProjectionStep, "project"); + +std::string ProjectionStep::kind() const { return "project"; } + +bool ProjectionStep::accepts_options() const { return true; } + +bool ProjectionStep::requires_options() const { return true; } + +void ProjectionStep::set_options(const nlohmann::json &options) { + if (!options.is_object()) { + throw Exception(kind() + " expects a JSON object for its options!"); + } + + for (const auto &[key, value] : options.items()) { + if (key == "method") { + if (!value.is_string()) { + throw Exception("Value for " + kind() + " option '" + key + + "' must be a string"); + } + if (value == "biorthogonal") { + // Since we don't support any other method for now, we don't have to + // store this option + } else { + throw Exception("Invalid value for " + kind() + " option '" + key + + "': '" + value.get() + "'"); + } + } else { + throw Exception("Unknown option key for " + kind() + ": '" + key + "'"); + } + } +} + +std::size_t ProjectionStep::process(std::string_view id_prefix, + std::size_t id_start, ExecutionContext &ctx, + const ExpressionData &data) { + container::svector transformed; + + transformed.insert(transformed.end(), data.expressions.begin(), + data.expressions.end()); + + container::svector> symmetrizers; + symmetrizers.reserve(transformed.size()); + for (ResultExpr ¤t : transformed) { + symmetrizers.push_back(pop_symmetrizer(current)); + } + + mbpt::biorthogonal_transform(transformed); + + for (std::size_t i = 0; i < transformed.size(); ++i) { + if (symmetrizers.at(i).has_value()) { + transformed.at(i).expression() = + ex(ExprPtrList{std::move(symmetrizers.at(i).value()), + std::move(transformed.at(i).expression())}, + Product::Flatten::No); + } + + simplify(transformed.at(i)); + } + + ExpressionData data_obj; + data_obj.expressions.insert(data_obj.expressions.end(), + std::make_move_iterator(transformed.begin()), + std::make_move_iterator(transformed.end())); + ctx.set_data(id_prefix, id_start, std::move(data_obj)); + + return 1; +} + +} // namespace sequant::util::extint diff --git a/utilities/external-interface/projection_step.hpp b/utilities/external-interface/projection_step.hpp new file mode 100644 index 0000000000..300ddba97e --- /dev/null +++ b/utilities/external-interface/projection_step.hpp @@ -0,0 +1,31 @@ +#ifndef SEQUANT_EXTERNAL_INTERFACE_PROJECTIONSTEP_HPP +#define SEQUANT_EXTERNAL_INTERFACE_PROJECTIONSTEP_HPP + +#include "execution_context.hpp" +#include "processing_data.hpp" +#include "processing_step.hpp" + +#include + +#include +#include + +namespace sequant::util::extint { + +class ProjectionStep : public OneByOneProcessingStep { + public: + std::string kind() const override; + + bool accepts_options() const override; + bool requires_options() const override; + void set_options(const nlohmann::json &options) override; + + protected: + std::size_t process(std::string_view id_prefix, std::size_t id_start, + ExecutionContext &ctx, + const ExpressionData &data) override; +}; + +} // namespace sequant::util::extint + +#endif // SEQUANT_EXTERNAL_INTERFACE_PROJECTIONSTEP_HPP diff --git a/utilities/external-interface/read_input_step.cpp b/utilities/external-interface/read_input_step.cpp new file mode 100644 index 0000000000..33a46a63b1 --- /dev/null +++ b/utilities/external-interface/read_input_step.cpp @@ -0,0 +1,98 @@ +#include "read_input_step.hpp" +#include "processing_step_factory.hpp" + +#include +#include +#include + +#include + +#include +#include +#include +#include + +namespace sequant::util::extint { + +SEQUANT_EXTINT_REGISTER_STEP_TYPE(ReadInputStep, "read_input"); + +std::string ReadInputStep::kind() const { return "read_input"; } + +bool ReadInputStep::accepts_options() const { return true; } + +bool ReadInputStep::requires_options() const { return true; } + +void ReadInputStep::set_options(const nlohmann::json &options) { + if (!options.is_object()) { + throw Exception(kind() + " expects a JSON object for its options!"); + } + + for (const auto &[key, value] : options.items()) { + if (key == "file_path") { + if (value.is_array()) { + for (const nlohmann::json ¤t : value) { + if (!current.is_string()) { + throw Exception("Array entry for " + kind() + " option '" + key + + "' must be strings"); + } + + input_paths_.emplace_back(current.get()); + } + } else if (value.is_string()) { + input_paths_.emplace_back(value.get()); + } else { + throw Exception("Invalid data type for " + kind() + " option '" + key + + "'"); + } + } else if (key == "default_symmetry") { + if (!value.is_string()) { + throw Exception("Value for " + kind() + " option '" + key + + "' must be a string"); + } + if (value == "none") { + options_.def_perm_symm = Symmetry::Nonsymm; + options_.def_col_symm = ColumnSymmetry::Nonsymm; + options_.def_braket_symm = BraKetSymmetry::Nonsymm; + } else if (value == "antisymmetric") { + options_.def_perm_symm = Symmetry::Antisymm; + options_.def_col_symm = ColumnSymmetry::Symm; + options_.def_braket_symm = BraKetSymmetry::Nonsymm; + } else if (value == "symmetric") { + options_.def_perm_symm = Symmetry::Symm; + options_.def_col_symm = ColumnSymmetry::Symm; + options_.def_braket_symm = BraKetSymmetry::Nonsymm; + } else { + throw Exception("Invalid value for " + kind() + " option '" + key + + "': '" + value.get() + "'"); + } + } else { + throw Exception("Unknown option key for " + kind() + ": '" + key + "'"); + } + } +} + +std::size_t ReadInputStep::process(std::string_view id_prefix, + std::size_t id_start, + ExecutionContext &ctx) { + std::size_t counter = 0; + + for (const std::filesystem::path ¤t : input_paths_) { + if (!std::filesystem::exists(current)) { + throw Exception("Input file '" + current.string() + "' does not exist"); + } + + // Read input file + std::ifstream in(current); + const std::string contents(std::istreambuf_iterator(in), {}); + + ResultExpr expr = + io::serialization::from_string(contents, options_); + + ctx.set_data(id_prefix, id_start + counter++, + ExpressionData{.expressions = {std::move(expr)}}); + } + + return counter; +} + +} // namespace sequant::util::extint diff --git a/utilities/external-interface/read_input_step.hpp b/utilities/external-interface/read_input_step.hpp new file mode 100644 index 0000000000..833df5de87 --- /dev/null +++ b/utilities/external-interface/read_input_step.hpp @@ -0,0 +1,37 @@ +#ifndef SEQUANT_EXTERNAL_INTERFACE_READINPUTSTEP_HPP +#define SEQUANT_EXTERNAL_INTERFACE_READINPUTSTEP_HPP + +#include "execution_context.hpp" +#include "processing_step.hpp" + +#include + +#include + +#include +#include +#include +#include + +namespace sequant::util::extint { + +class ReadInputStep : public OneByOneProcessingStep { + public: + std::string kind() const override; + + bool accepts_options() const override; + bool requires_options() const override; + void set_options(const nlohmann::json &options) override; + + protected: + std::size_t process(std::string_view id_prefix, std::size_t id_start, + ExecutionContext &ctx) override; + + private: + std::vector input_paths_; + io::serialization::DeserializationOptions options_; +}; + +} // namespace sequant::util::extint + +#endif // SEQUANT_EXTERNAL_INTERFACE_READINPUTSTEP_HPP diff --git a/utilities/external-interface/simplify_step.cpp b/utilities/external-interface/simplify_step.cpp new file mode 100644 index 0000000000..425efce8c1 --- /dev/null +++ b/utilities/external-interface/simplify_step.cpp @@ -0,0 +1,53 @@ +#include "simplify_step.hpp" +#include "processing_data.hpp" +#include "processing_step_factory.hpp" + +#include + +#include + +#include +#include + +namespace sequant::util::extint { + +SEQUANT_EXTINT_REGISTER_STEP_TYPE(SimplifyStep, "simplify"); + +std::string SimplifyStep::kind() const { return "simplify"; } + +bool SimplifyStep::accepts_options() const { return false; } + +bool SimplifyStep::requires_options() const { return false; } + +void SimplifyStep::set_options(const nlohmann::json &) { + throw Exception(kind() + " doesn't take any options"); +} + +std::size_t SimplifyStep::process(std::string_view id_prefix, + std::size_t id_start, ExecutionContext &ctx, + const ExpressionData &data) { + std::vector outputs; + + for (const ResultExpr &expr : data.expressions) { + ResultExpr clone = expr.clone(); + simplify(clone); + + outputs.emplace_back(std::move(clone)); + } + + if (outputs != data.expressions) { + ExpressionData data_obj; + data_obj.expressions.insert(data_obj.expressions.end(), + std::make_move_iterator(outputs.begin()), + std::make_move_iterator(outputs.end())); + ctx.set_data(id_prefix, id_start, std::move(data_obj)); + + return 1; + } + + return 0; +} + +bool SimplifyStep::alias_unchanged_inputs() const { return true; } + +} // namespace sequant::util::extint diff --git a/utilities/external-interface/simplify_step.hpp b/utilities/external-interface/simplify_step.hpp new file mode 100644 index 0000000000..b10a13483c --- /dev/null +++ b/utilities/external-interface/simplify_step.hpp @@ -0,0 +1,35 @@ +#ifndef SEQUANT_EXTERNAL_INTERFACE_SIMPLIFYSTEP_HPP +#define SEQUANT_EXTERNAL_INTERFACE_SIMPLIFYSTEP_HPP + +#include "execution_context.hpp" +#include "processing_data.hpp" +#include "processing_step.hpp" + +#include + +#include + +#include +#include + +namespace sequant::util::extint { + +class SimplifyStep : public OneByOneProcessingStep { + public: + std::string kind() const override; + + bool accepts_options() const override; + bool requires_options() const override; + void set_options(const nlohmann::json &options) override; + + protected: + std::size_t process(std::string_view id_prefix, std::size_t id_start, + ExecutionContext &ctx, + const ExpressionData &data) override; + + bool alias_unchanged_inputs() const override; +}; + +} // namespace sequant::util::extint + +#endif // SEQUANT_EXTERNAL_INTERFACE_SIMPLIFYSTEP_HPP diff --git a/utilities/external-interface/spintracing_step.cpp b/utilities/external-interface/spintracing_step.cpp new file mode 100644 index 0000000000..7f29e80c01 --- /dev/null +++ b/utilities/external-interface/spintracing_step.cpp @@ -0,0 +1,77 @@ +#include "spintracing_step.hpp" +#include "processing_data.hpp" +#include "processing_step_factory.hpp" + +#include +#include +#include + +#include + +#include +#include + +namespace sequant::util::extint { + +SEQUANT_EXTINT_REGISTER_STEP_TYPE(SpintracingStep, "spintracing"); + +std::string SpintracingStep::kind() const { return "spintracing"; } + +bool SpintracingStep::accepts_options() const { return true; } + +bool SpintracingStep::requires_options() const { return false; } + +void SpintracingStep::set_options(const nlohmann::json &options) { + if (!options.is_object()) { + throw Exception(kind() + " expects a JSON object for its options!"); + } + + for (const auto &[key, value] : options.items()) { + if (key == "algorithm") { + if (!value.is_string()) { + throw Exception("Option '" + key + "' for " + kind() + + " requires string argument"); + } + + if (value == "rigorous") { + use_closed_shell_algo_ = false; + } else if (value == "closed_shell") { + use_closed_shell_algo_ = true; + } else { + throw Exception("Unknown value '" + value.get() + + "' for option " + key + " of " + kind()); + } + } else { + throw Exception("Unknown option key for " + kind() + ": '" + key + "'"); + } + } +} + +std::size_t SpintracingStep::process(std::string_view id_prefix, + std::size_t id_start, + ExecutionContext &ctx, + const ExpressionData &data) { + std::size_t num_outputs = 0; + + for (const ResultExpr &expr : data.expressions) { + container::svector result; + if (use_closed_shell_algo_) { + result = mbpt::closed_shell_spintrace(expr); + } else { + result = mbpt::spintrace(expr); + } + + ExpressionData output; + output.expressions.insert(output.expressions.end(), + std::make_move_iterator(result.begin()), + std::make_move_iterator(result.end())); + + ctx.set_data(id_prefix, id_start + num_outputs, std::move(output)); + + num_outputs += 1; + } + + return num_outputs; +} + +} // namespace sequant::util::extint diff --git a/utilities/external-interface/spintracing_step.hpp b/utilities/external-interface/spintracing_step.hpp new file mode 100644 index 0000000000..9335123b6e --- /dev/null +++ b/utilities/external-interface/spintracing_step.hpp @@ -0,0 +1,36 @@ +#ifndef SEQUANT_EXTERNAL_INTERFACE_SPINTRACINGSTEP_HPP +#define SEQUANT_EXTERNAL_INTERFACE_SPINTRACINGSTEP_HPP + +#include "execution_context.hpp" +#include "processing_data.hpp" +#include "processing_step.hpp" + +#include + +#include + +#include +#include + +namespace sequant::util::extint { + +class SpintracingStep : public OneByOneProcessingStep { + public: + std::string kind() const override; + + bool accepts_options() const override; + bool requires_options() const override; + void set_options(const nlohmann::json &options) override; + + protected: + std::size_t process(std::string_view id_prefix, std::size_t id_start, + ExecutionContext &ctx, + const ExpressionData &data) override; + + private: + bool use_closed_shell_algo_ = false; +}; + +} // namespace sequant::util::extint + +#endif // SEQUANT_EXTERNAL_INTERFACE_SPINTRACINGSTEP_HPP diff --git a/utilities/external-interface/to_export_tree_step.cpp b/utilities/external-interface/to_export_tree_step.cpp new file mode 100644 index 0000000000..2d8115d10d --- /dev/null +++ b/utilities/external-interface/to_export_tree_step.cpp @@ -0,0 +1,74 @@ +#include "to_export_tree_step.hpp" +#include "processing_data.hpp" +#include "processing_step_factory.hpp" +#include "utils.hpp" + +#include +#include +#include +#include + +#include + +#include +#include + +namespace sequant::util::extint { + +SEQUANT_EXTINT_REGISTER_STEP_TYPE(ToExportTreeStep, "to_export_tree"); + +std::string ToExportTreeStep::kind() const { return "to_export_tree"; } + +bool ToExportTreeStep::accepts_options() const { return false; } + +bool ToExportTreeStep::requires_options() const { return false; } + +void ToExportTreeStep::set_options(const nlohmann::json &) { + throw Exception(kind() + " doesn't take any options"); +} + +std::size_t ToExportTreeStep::process(std::string_view id_prefix, + std::size_t id_start, + ExecutionContext &ctx, + const ExpressionData &data) { + ExportTreeData output; + + for (const ResultExpr ¤t : data.expressions) { + std::optional symmetrized_result; + + ExportNode<> tree = [&]() { + if (!needsSymmetrization(current.expression())) { + return to_export_tree(current); + } + + ResultExpr copy = current.clone(); + [[maybe_unused]] std::optional symmetrizer = + pop_symmetrizer(copy); + + SEQUANT_ASSERT(copy.produces_tensor()); + SEQUANT_ASSERT(copy.has_label()); + + SEQUANT_ASSERT(symmetrizer.has_value()); + SEQUANT_ASSERT(symmetrizer->is()); + SEQUANT_ASSERT(symmetrizer->as().bra() == copy.ket()); + SEQUANT_ASSERT(symmetrizer->as().ket() == copy.bra()); + + symmetrized_result = copy.result_as_tensor(); + + copy.result_as_tensor().set_label( + std::wstring(symmetrized_result->label()) + L"u"); + + return to_export_tree(copy); + }(); + + output.entries.push_back( + {.tree = std::move(tree), + .symm_contribution_target = std::move(symmetrized_result)}); + } + + ctx.set_data(id_prefix, id_start, std::move(output)); + + return 1; +} + +} // namespace sequant::util::extint diff --git a/utilities/external-interface/to_export_tree_step.hpp b/utilities/external-interface/to_export_tree_step.hpp new file mode 100644 index 0000000000..323aadf082 --- /dev/null +++ b/utilities/external-interface/to_export_tree_step.hpp @@ -0,0 +1,31 @@ +#ifndef SEQUANT_EXTERNAL_INTERFACE_TO_EXPORT_TREE_STEP_HPP +#define SEQUANT_EXTERNAL_INTERFACE_TO_EXPORT_TREE_STEP_HPP + +#include "execution_context.hpp" +#include "processing_data.hpp" +#include "processing_step.hpp" + +#include + +#include +#include + +namespace sequant::util::extint { + +class ToExportTreeStep : public OneByOneProcessingStep { + public: + std::string kind() const override; + + bool accepts_options() const override; + bool requires_options() const override; + void set_options(const nlohmann::json &options) override; + + protected: + std::size_t process(std::string_view id_prefix, std::size_t id_start, + ExecutionContext &ctx, + const ExpressionData &data) override; +}; + +} // namespace sequant::util::extint + +#endif // SEQUANT_EXTERNAL_INTERFACE_TO_EXPORT_TREE_STEP_HPP diff --git a/utilities/external-interface/utils.cpp b/utilities/external-interface/utils.cpp index b6fb576336..7aa88e7646 100644 --- a/utilities/external-interface/utils.cpp +++ b/utilities/external-interface/utils.cpp @@ -10,7 +10,6 @@ #include #include -#include #include #include diff --git a/utilities/external-interface/validate_step.cpp b/utilities/external-interface/validate_step.cpp new file mode 100644 index 0000000000..35fab92d31 --- /dev/null +++ b/utilities/external-interface/validate_step.cpp @@ -0,0 +1,41 @@ +#include "validate_step.hpp" +#include "processing_step_factory.hpp" + +#include +#include +#include + +#include +#include + +namespace sequant::util::extint { + +SEQUANT_EXTINT_REGISTER_STEP_TYPE(ValidateStep, "validate"); + +std::string ValidateStep::kind() const { return "validate"; } + +bool ValidateStep::accepts_options() const { return false; } + +bool ValidateStep::requires_options() const { return false; } + +void ValidateStep::set_options(const nlohmann::json &) { + throw Exception("validate doesn't take any options"); +} + +std::size_t ValidateStep::process(std::string_view, std::size_t, + ExecutionContext &, + const ExpressionData &data) { + std::size_t expr_counter = 1; + for (const ResultExpr &expr : data.expressions) { + std::string msg; + if (!is_valid(expr, &msg)) { + throw Exception("Expression #" + std::to_string(expr_counter) + + " is invalid: " + msg); + } + ++expr_counter; + } + + return 0; +} + +} // namespace sequant::util::extint diff --git a/utilities/external-interface/validate_step.hpp b/utilities/external-interface/validate_step.hpp new file mode 100644 index 0000000000..65ec929737 --- /dev/null +++ b/utilities/external-interface/validate_step.hpp @@ -0,0 +1,30 @@ +#ifndef SEQUANT_EXTERNAL_INTERFACE_VALIDATESTEP_HPP +#define SEQUANT_EXTERNAL_INTERFACE_VALIDATESTEP_HPP + +#include "execution_context.hpp" +#include "processing_data.hpp" +#include "processing_step.hpp" + +#include + +#include +#include + +namespace sequant::util::extint { + +class ValidateStep : public OneByOneProcessingStep { + public: + std::string kind() const override; + + bool accepts_options() const override; + bool requires_options() const override; + void set_options(const nlohmann::json &options) override; + + std::size_t process(std::string_view id_prefix, std::size_t id_start, + ExecutionContext &ctx, + const ExpressionData &data) override; +}; + +} // namespace sequant::util::extint + +#endif // SEQUANT_EXTERNAL_INTERFACE_VALIDATESTEP_HPP