diff --git a/CMakeLists.txt b/CMakeLists.txt index b04b376fc9..f119f10dbc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -306,16 +306,28 @@ set(SeQuant_symb_src SeQuant/core/context.hpp SeQuant/core/expressions/abstract_tensor.cpp SeQuant/core/expressions/abstract_tensor.hpp + SeQuant/core/expressions/constant.cpp + SeQuant/core/expressions/constant.hpp SeQuant/core/expressions/expr.cpp SeQuant/core/expressions/expr.hpp SeQuant/core/expressions/expr_algorithms.cpp SeQuant/core/expressions/expr_algorithms.hpp SeQuant/core/expressions/expr_operators.hpp + SeQuant/core/expressions/expr_ptr.cpp + SeQuant/core/expressions/expr_ptr.hpp SeQuant/core/expressions/expr_range.hpp + SeQuant/core/expressions/power.cpp + SeQuant/core/expressions/power.hpp + SeQuant/core/expressions/product.cpp + SeQuant/core/expressions/product.hpp SeQuant/core/expressions/result_expr.cpp SeQuant/core/expressions/result_expr.hpp + SeQuant/core/expressions/sum.cpp + SeQuant/core/expressions/sum.hpp SeQuant/core/expressions/tensor.cpp SeQuant/core/expressions/tensor.hpp + SeQuant/core/expressions/variable.cpp + SeQuant/core/expressions/variable.hpp SeQuant/core/hash.cpp SeQuant/core/hash.hpp SeQuant/core/hugenholtz.hpp diff --git a/SeQuant/core/expressions/constant.cpp b/SeQuant/core/expressions/constant.cpp new file mode 100644 index 0000000000..04bb759c54 --- /dev/null +++ b/SeQuant/core/expressions/constant.cpp @@ -0,0 +1,77 @@ +#include +#include +#include +#include +#include + +namespace sequant { + +std::wstring Constant::to_latex() const { + return L"{" + io::latex::to_string(value()) + L"}"; +} + +Expr::type_id_type Constant::type_id() const { return get_type_id(); } + +bool Constant::is_scalar() const { return true; } + +ExprPtr Constant::clone() const { return ex(this->value()); } + +void Constant::adjoint() { + value_ = conj(value_); + reset_hash_value(); +} + +Constant &Constant::operator*=(const Expr &that) { + if (that.is()) { + value_ *= that.as().value(); + } else { + throw Exception("Constant::operator*=(that): not valid for that"); + } + + reset_hash_value(); + + return *this; +} + +Constant &Constant::operator+=(const Expr &that) { + if (that.is()) { + value_ += that.as().value(); + } else { + throw Exception("Constant::operator+=(that): not valid for that"); + } + + reset_hash_value(); + + return *this; +} + +Constant &Constant::operator-=(const Expr &that) { + if (that.is()) { + value_ -= that.as().value(); + } else { + throw Exception("Constant::operator-=(that): not valid for that"); + } + + reset_hash_value(); + + return *this; +} + +bool Constant::is_zero(scalar_type v) { return v.is_zero(); } + +bool Constant::is_zero() const { return is_zero(this->value()); } + +Expr::hash_type Constant::memoizing_hash() const { + if (!hash_value_) { + hash_value_ = hash::value(value_); + } else { + SEQUANT_ASSERT(*hash_value_ == hash::value(value_)); + } + return *hash_value_; +} + +bool Constant::static_equal(const Expr &that) const { + return value() == static_cast(that).value(); +} + +} // namespace sequant diff --git a/SeQuant/core/expressions/constant.hpp b/SeQuant/core/expressions/constant.hpp index 9ae67a9cb2..9bade3b542 100644 --- a/SeQuant/core/expressions/constant.hpp +++ b/SeQuant/core/expressions/constant.hpp @@ -3,8 +3,6 @@ #include #include -#include -#include #include #include @@ -14,6 +12,8 @@ namespace sequant { +class ExprPtr; + // implementation details of Constant; prefer sequant::detail over an unnamed // namespace in a header (see CppCoreGuidelines SF.21) namespace detail { @@ -67,68 +67,37 @@ class Constant : public Expr { throw Exception("Constant::value: cannot convert value to type T"); } - std::wstring to_latex() const override { - return L"{" + io::latex::to_string(value()) + L"}"; - } + std::wstring to_latex() const override; - type_id_type type_id() const override { return get_type_id(); } + type_id_type type_id() const override; - bool is_scalar() const override { return true; } + bool is_scalar() const override; - ExprPtr clone() const override { return ex(this->value()); } + ExprPtr clone() const override; /// @brief adjoint of a Constant is its complex conjugate virtual void adjoint() override; - virtual Expr &operator*=(const Expr &that) override { - if (that.is()) { - value_ *= that.as().value(); - } else { - throw Exception("Constant::operator*=(that): not valid for that"); - } - return *this; - } + Constant &operator*=(const Expr &that); - virtual Expr &operator+=(const Expr &that) override { - if (that.is()) { - value_ += that.as().value(); - } else { - throw Exception("Constant::operator+=(that): not valid for that"); - } - return *this; - } + Constant &operator+=(const Expr &that); - virtual Expr &operator-=(const Expr &that) override { - if (that.is()) { - value_ -= that.as().value(); - } else { - throw Exception("Constant::operator-=(that): not valid for that"); - } - return *this; - } + Constant &operator-=(const Expr &that); /// @param[in] v a scalar /// @return true if this is zero - static bool is_zero(scalar_type v) { return v.is_zero(); } + static bool is_zero(scalar_type v); /// @return `Constant::is_zero(this->value())` - bool is_zero() const final { return is_zero(this->value()); } + bool is_zero() const final; private: scalar_type value_; - hash_type memoizing_hash() const override { - if (!hash_value_) { - hash_value_ = hash::value(value_); - } else { - SEQUANT_ASSERT(*hash_value_ == hash::value(value_)); - } - return *hash_value_; - } + hash_type memoizing_hash() const override; + + bool static_equal(const Expr &that) const override; - bool static_equal(const Expr &that) const override { - return value() == static_cast(that).value(); - } }; // class Constant } // namespace sequant diff --git a/SeQuant/core/expressions/expr.cpp b/SeQuant/core/expressions/expr.cpp index a8ce23862e..cb9655e12b 100644 --- a/SeQuant/core/expressions/expr.cpp +++ b/SeQuant/core/expressions/expr.cpp @@ -2,32 +2,11 @@ // Created by Eduard Valeyev on 2019-02-06. // -#include -#include #include -#include +#include #include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include +#include +#include namespace sequant { @@ -77,493 +56,8 @@ ExprPtr &Expr::back() { return at(size() - 1); } const ExprPtr &Expr::back() const { return at(size() - 1); } -ExprPtr ExprPtr::clone() const & { - if (!*this) return {}; - return ExprPtr(as_shared_ptr()->clone()); -} - -ExprPtr ExprPtr::clone() && noexcept { return std::move(*this); } - -ExprPtr::base_type &ExprPtr::as_shared_ptr() & { - return static_cast(*this); -} -const ExprPtr::base_type &ExprPtr::as_shared_ptr() const & { - return static_cast(*this); -} -ExprPtr::base_type &&ExprPtr::as_shared_ptr() && { - return static_cast(*this); -} - -Expr &ExprPtr::operator*() & { - SEQUANT_ASSERT(this->operator bool()); - return *(this->get()); -} - -const Expr &ExprPtr::operator*() const & { - SEQUANT_ASSERT(this->operator bool()); - return *(this->get()); -} - -Expr &&ExprPtr::operator*() && { - SEQUANT_ASSERT(this->operator bool()); - return std::move(*(this->get())); -} - -ExprPtr &ExprPtr::operator+=(const ExprPtr &other) { - if (!other) return *this; - - if (!*this) { - *this = other.clone(); - } else if (as_shared_ptr()->is()) { - as_shared_ptr()->operator+=(*other); - } else if (as_shared_ptr()->is() && other->is()) { - *this = ex(this->as().value() + - other->as().value()); - } else { - *this = ex(ExprPtrList{*this, other}); - } - return *this; -} - -ExprPtr &ExprPtr::operator-=(const ExprPtr &other) { - if (!other) return *this; - - if (!*this) { - *this = ex(-1) * other.clone(); - } else if (as_shared_ptr()->is()) { - as_shared_ptr()->operator-=(*other); - } else if (as_shared_ptr()->is() && other->is()) { - *this = ex(this->as().value() - - other->as().value()); - } else { - *this = ex(ExprPtrList{*this, ex(-1, ExprPtrList{other})}); - } - return *this; -} - -ExprPtr &ExprPtr::operator*=(const ExprPtr &other) { - if (!other) return *this; - - if (!*this) { - *this = other.clone(); - } else if (as_shared_ptr()->is()) { - as_shared_ptr()->operator*=(*other); - } else if (as_shared_ptr()->is() && other->is()) { - *this = ex(this->as().value() * - other->as().value()); - } else { - *this = ex(ExprPtrList{*this, other}); - } - return *this; -} - -std::size_t ExprPtr::size() const { return this->get()->size(); } - -std::wstring ExprPtr::to_latex() const { return as_shared_ptr()->to_latex(); } - -Exception Expr::not_implemented(const char *fn) const { - std::ostringstream oss; - oss << "Expr::" << fn - << " not implemented in this derived class (type_name=" << type_name() - << ")"; - return Exception(oss.str()); -} - -std::wstring Expr::to_latex() const { throw not_implemented("to_latex"); } - -ExprPtr Expr::clone() const { throw not_implemented("clone"); } - -void Expr::adjoint() { throw not_implemented("adjoint"); } - -Expr &Expr::operator*=(const Expr &) { throw not_implemented("operator*="); } - -Expr &Expr::operator^=(const Expr &) { throw not_implemented("operator^="); } - -Expr &Expr::operator+=(const Expr &) { throw not_implemented("operator+="); } - -Expr &Expr::operator-=(const Expr &) { throw not_implemented("operator-="); } - -ExprPtr adjoint(const ExprPtr &expr) { - auto result = expr->clone(); - result->adjoint(); - return result; -} - -void Constant::adjoint() { - value_ = conj(value_); - reset_hash_value(); -} - -std::wstring_view Variable::label() const { return label_; } - -void Variable::set_label(std::wstring label) { - label_ = std::move(label); - reset_hash_value(); -} - -void Variable::conjugate() { conjugated_ = !conjugated_; } - -bool Variable::conjugated() const { return conjugated_; } - -std::wstring Variable::to_latex() const { - std::wstring result = L"{" + io::latex::utf_to_string(label_) + L"}"; - if (conjugated_) result = L"{" + result + L"^*" + L"}"; - return result; -} - -ExprPtr Variable::clone() const { return ex(*this); } - -void Variable::adjoint() { conjugate(); } - -bool Product::is_commutative() const { - bool result = true; - const auto nfactors = size(); - for (size_t f = 0; f != nfactors; ++f) { - for (size_t s = 1; result && s != nfactors; ++s) { - result &= factors_[f]->commutes_with(*factors_[s]); - } - } - return result; -} - -ExprPtr Product::canonicalize_impl(CanonicalizeOptions opts) { - // recursively canonicalize non-tensor subfactors (tensors will be - // canonicalized as part of the TN built of all tensor factors of this) ... - ranges::for_each(factors_, [this, opts](auto &factor) { - if (factor.template is()) { - return; - } - auto bp = factor->canonicalize(opts); - if (bp) { - SEQUANT_ASSERT(bp->template is()); - this->scalar_ *= std::static_pointer_cast(bp)->value(); - } - }); - - if (Logger::instance().canonicalize) { - std::wcout << "Product canonicalization(" << to_wstring(opts.method) - << ") input: " << to_latex() << std::endl; - } - - // pull out all scalar factors to the front - auto is_scalar = [](const auto &factor) { return factor->is_scalar(); }; - auto scalars = - factors_ | ranges::views::filter(is_scalar) | ranges::to_vector; - // scalars commute, so we can reorder them freely - ranges::sort(scalars, [](const auto &first, const auto &second) { - return *first < *second; - }); - - factors_ = factors_ | ranges::views::filter([&is_scalar](const auto &factor) { - return !is_scalar(factor); - }) | - ranges::to; - - // if there are no factors, insert scalars back and return - if (factors_.empty()) { - factors_.insert(factors_.begin(), scalars.begin(), scalars.end()); - return {}; - } - - auto contains_nontensors = ranges::any_of(factors_, [](const auto &factor) { - return std::dynamic_pointer_cast(factor) == nullptr; - }); - if (!contains_nontensors) { // tensor network canonization is a special case - // that's done in - // TensorNetwork - auto make_canonical_tn = [this, &opts](auto *tn_null_ptr) { - using TN = std::decay_t>; - ExprPtr canon_factor; - TN tn(this->factors_); - if constexpr (TN::version() == 3) { - canon_factor = tn.canonicalize( - TensorCanonicalizer::cardinal_tensor_labels(), opts); - } else { - using NamedIndexSet = tensor_network::NamedIndexSet; - std::shared_ptr named_indices = - !opts.named_indices - ? nullptr - : std::make_shared(opts.named_indices->begin(), - opts.named_indices->end()); - canon_factor = tn.canonicalize( - TensorCanonicalizer::cardinal_tensor_labels(), - opts.method == CanonicalizationMethod::Rapid, named_indices.get()); - } - return std::pair{std::move(tn), canon_factor}; - }; - using TN = TensorNetwork; - auto [tn, canon_factor] = make_canonical_tn(static_cast(nullptr)); - - const auto &tensors = tn.tensors(); - using std::size; - SEQUANT_ASSERT(size(tensors) == size(factors_)); - using std::begin; - using std::end; - std::transform(begin(tensors), end(tensors), begin(factors_), - [](const auto &tptr) { - auto exprptr = std::dynamic_pointer_cast(tptr); - SEQUANT_ASSERT(exprptr); - return exprptr; - }); - if (canon_factor) scalar_ *= canon_factor->template as().value(); - this->reset_hash_value(); - } else { // if contains non-tensors, do commutation-checking resort - - // comparer that respects cardinal tensor labels - auto &cardinal_tensor_labels = - TensorCanonicalizer::cardinal_tensor_labels(); - auto local_compare = [&cardinal_tensor_labels](const ExprPtr &first, - const ExprPtr &second) { - if (first->is() && second->is()) { - const auto first_label = first->as().label(); - const auto second_label = second->as().label(); - if (first_label == second_label) return *first < *second; - const auto first_is_cardinal_it = ranges::find_if( - cardinal_tensor_labels, - [&first_label](const std::wstring &l) { return l == first_label; }); - const auto first_is_cardinal = - first_is_cardinal_it != ranges::end(cardinal_tensor_labels); - const auto second_is_cardinal_it = ranges::find_if( - cardinal_tensor_labels, [&second_label](const std::wstring &l) { - return l == second_label; - }); - const auto second_is_cardinal = - second_is_cardinal_it != ranges::end(cardinal_tensor_labels); - if (first_is_cardinal && second_is_cardinal) - return first_is_cardinal_it < second_is_cardinal_it; - else if (first_is_cardinal && !second_is_cardinal) - return true; - else if (!first_is_cardinal && second_is_cardinal) - return false; - else { - SEQUANT_ASSERT(!first_is_cardinal && !second_is_cardinal); - return *first < *second; - } - } else - return *first < *second; - }; - - // ... then resort, respecting commutativity - using std::begin; - using std::end; - if (static_commutativity()) { - if (is_commutative()) { - std::stable_sort(begin(factors_), end(factors_), local_compare); - } - } else { - // must do bubble sort if not commuting to avoid swapping elements across - // a noncommuting element - bubble_sort( - begin(factors_), end(factors_), - [&local_compare](const ExprPtr &first, const ExprPtr &second) { - bool result = (first->commutes_with(*second)) - ? local_compare(first, second) - : false; - return result; - }); - } - } - // reinsert scalar factors at the front - factors_.insert(factors_.begin(), scalars.begin(), scalars.end()); - - // TODO evaluate product of Tensors (turn this into Products of Products) - - if (Logger::instance().canonicalize) - std::wcout << "Product canonicalization(" << to_wstring(opts.method) - << ") result: " << to_latex() << std::endl; - - return {}; // side effects are absorbed into the scalar_ -} - -void Product::adjoint() { - SEQUANT_ASSERT(static_commutativity() == false); // assert no slicing - auto adj_scalar = conj(scalar()); - using namespace ranges; - auto adj_factors = - factors() | views::reverse | - views::transform([](auto &expr) { return ::sequant::adjoint(expr); }); - using std::swap; - *this = - Product(adj_scalar, ranges::begin(adj_factors), ranges::end(adj_factors)); -} - -ExprPtr Product::canonicalize(CanonicalizeOptions opt) { - return this->canonicalize_impl(opt); -} - -ExprPtr Product::rapid_canonicalize(CanonicalizeOptions opt) { - SEQUANT_ASSERT(opt.method == CanonicalizationMethod::Rapid); - return this->canonicalize_impl(opt); -} - -void CProduct::adjoint() { - auto adj_scalar = conj(scalar()); - using namespace ranges; - // no need to reverse for commutative product - auto adj_factors = factors() | views::transform([](auto &&expr) { - return ::sequant::adjoint(expr); - }); - *this = CProduct(adj_scalar, ranges::begin(adj_factors), - ranges::end(adj_factors)); -} - -void NCProduct::adjoint() { - auto adj_scalar = conj(scalar()); - using namespace ranges; - // no need to reverse for commutative product - auto adj_factors = - factors() | views::reverse | - views::transform([](auto &&expr) { return ::sequant::adjoint(expr); }); - *this = NCProduct(adj_scalar, ranges::begin(adj_factors), - ranges::end(adj_factors)); -} - -void Sum::adjoint() { - using namespace ranges; - auto adj_summands = summands() | views::transform([](auto &&expr) { - return ::sequant::adjoint(expr); - }); - *this = Sum(ranges::begin(adj_summands), ranges::end(adj_summands)); -} - -ExprPtr Sum::canonicalize_impl(bool multipass, CanonicalizeOptions opts) { - if (Logger::instance().canonicalize) - std::wcout << "Sum::canonicalize_impl: input = " - << to_latex_align(shared_from_this()) << std::endl; - - const auto npasses = multipass ? 2 : 1; - for (auto pass = 0; pass != npasses; ++pass) { - const auto rapid = (pass % 2 == 0); - - // canonicalizing TNs in a sum requires treating named indices as - // meaningful/distinct - auto opts_copy = opts; - opts_copy.ignore_named_index_labels = - CanonicalizeOptions::IgnoreNamedIndexLabel::No; - if (rapid) { - opts_copy.method = CanonicalizationMethod::Lexicographic; - } else - opts_copy.method = opts.method | CanonicalizationMethod::Topological; - - // recursively canonicalize summands ... - // using for_each and direct access to summands - sequant::for_each(summands_, [pass, &opts_copy, &rapid](ExprPtr &summand) { - ExprPtr bp; - if (rapid) { - bp = summand->rapid_canonicalize(opts_copy); - } else { - bp = summand->canonicalize(opts_copy); - } - if (bp) { - SEQUANT_ASSERT(bp->template is()); - summand = ex(std::static_pointer_cast(bp)->value(), - ExprPtrList{summand}); - } - }); - if (Logger::instance().canonicalize) - std::wcout << "Sum::canonicalize_impl (pass=" << pass - << "): after canonicalizing summands = " - << to_latex_align(shared_from_this()) << std::endl; - - HashingAccumulator acc; - for (auto &summand : summands_) { - acc.append(summand); - } - - // last pass? sort by hash then by Expr::operator< - // N.B. no point in differentiating between canonicalization methods here - // since need to sort in both cases - auto new_sum = - (pass == npasses - 1) ? acc.make_canonicalized_sum() : acc.make_sum(); - this->swap(*new_sum); - - if (Logger::instance().canonicalize) - std::wcout << "Sum::canonicalize_impl (pass=" << pass - << "): after reducing summands = " - << to_latex_align(shared_from_this()) << std::endl; - } - - return {}; // side effects are absorbed into summands -} - -HashingAccumulator &HashingAccumulator::append(ExprPtr summand, bool flatten) { - // flatten, if needed - if (flatten && summand.is()) { - for (auto &subsummand : summand.as().summands()) { - this->append(subsummand, flatten); - } - return *this; - } - - // process summand as a whole - auto it = summands_.find(summand); - if (it == summands_.end()) { - summands_.emplace(summand); - } else { // found existing term with the same hash - auto existing_summand = *it; - if (summand.template is()) { - if (existing_summand.is()) { - // both are products - add them - existing_summand.as().add_identical( - summand.template as()); - } else { - // convert existing term to product and add - auto product_copy = std::make_shared(summand->clone()); - product_copy->add_identical(existing_summand); - summands_.erase(it); - summands_.emplace(std::move(product_copy)); - } - } else { - if (existing_summand.is()) { - existing_summand.as().add_identical(summand); - } else { - // neither is a product - create new product - auto product_form = std::make_shared(); - product_form->append(2, summand.template as()); - summands_.erase(it); - summands_.emplace(std::move(product_form)); - } - } - } - - return *this; -} - -SumPtr HashingAccumulator::make_sum_impl(bool canonicalize) { - Sum::summands_type summands; - summands.reserve(summands_.size()); - for (auto summand : summands_) { - if (!summand->is_zero()) { - summands.push_back(summand); - } - } - - if (canonicalize) { - ranges::sort(summands, [](const auto &e1, const auto &e2) { - if (e1->hash_value() == e2->hash_value()) { - return e1 < e2; - } else { - return e1->hash_value() < e2->hash_value(); - } - }); - } - - return std::make_shared(std::move(summands), Sum::move_only_tag{}); -} - -SumPtr HashingAccumulator::make_sum() { return make_sum_impl(false); } - -SumPtr HashingAccumulator::make_canonicalized_sum() { - return make_sum_impl(true); -} - -ExprPtr HashingAccumulator::make_expr(bool canonicalize) { - if (summands_.size() == 0) { - return ex(0); - } else if (summands_.size() == 1) - return *(summands_.begin()); - else - return make_sum_impl(canonicalize); +std::wstring Expr::to_latex() const { + throw Exception("to_latex not implemented for " + type_name()); } bool proportional_to::operator()(const ExprPtr &expr1, diff --git a/SeQuant/core/expressions/expr.hpp b/SeQuant/core/expressions/expr.hpp index 7bee7e8bc3..e1f7ba1870 100644 --- a/SeQuant/core/expressions/expr.hpp +++ b/SeQuant/core/expressions/expr.hpp @@ -75,9 +75,7 @@ class Expr : public std::enable_shared_from_this { virtual std::wstring to_latex() const; /// @return a clone of this object, i.e. an object that is equal to @c this - /// @note - must be overridden in the derived class. - /// - the default implementation throws an exception - virtual ExprPtr clone() const; + virtual ExprPtr clone() const = 0; /// like Expr::shared_from_this, but returns ExprPtr /// @return a shared_ptr to this object wrapped into ExprPtr, if this object @@ -228,9 +226,7 @@ class Expr : public std::enable_shared_from_this { } /// @brief changes this to its adjoint - /// @note base implementation throws, must be reimplemented in the derived - /// class - virtual void adjoint(); + virtual void adjoint() = 0; /// Computes and returns the hash value. If default @p hasher is used then the /// value will be memoized, otherwise @p hasher will be used to compute the @@ -247,17 +243,9 @@ class Expr : public std::enable_shared_from_this { } /// Computes and returns the derived type identifier - /// @note this function must be overridden in the derived class /// @sa Expr::get_type_id /// @return the hash value for this Expr - virtual type_id_type type_id() const -#if __GNUG__ - { - abort(); - } -#else - = 0; -#endif + virtual type_id_type type_id() const = 0; friend inline bool operator==(const Expr &a, const Expr &b); @@ -330,38 +318,6 @@ class Expr : public std::enable_shared_from_this { return boost::core::demangle(typeid(*this).name()); } - /** @name in-place arithmetic operators - * Virtual in-place arithmetic operators to be overridden in expressions for - * which these make sense. - */ - ///@{ - - /// @brief in-place multiply @c *this by @c that - /// @return reference to @c *this - /// @throw Exception if not implemented for this class, or cannot be - /// implemented for the particular @c that - virtual Expr &operator*=(const Expr &that); - - /// @brief in-place non-commutatively-multiply @c *this by @c that - /// @return reference to @c *this - /// @throw Exception if not implemented for this class, or cannot be - /// implemented for the particular @c that - virtual Expr &operator^=(const Expr &that); - - /// @brief in-place add @c that to @c *this - /// @return reference to @c *this - /// @throw Exception if not implemented for this class, or cannot be - /// implemented for the particular @c that - virtual Expr &operator+=(const Expr &that); - - /// @brief in-place subtract @c that from @c *this - /// @return reference to @c *this - /// @throw Exception if not implemented for this class, or cannot be - /// implemented for the particular @c that - virtual Expr &operator-=(const Expr &that); - - ///@} - ExprIterator begin(); ExprIterator end(); ConstExprIterator begin() const; @@ -441,14 +397,7 @@ class Expr : public std::enable_shared_from_this { /// @note @c that is guaranteed to be of same type as @c *this, hence can be /// statically cast /// @return true if @c that is equivalent to *this - virtual bool static_equal([[maybe_unused]] const Expr &that) const -#if __GNUG__ - { - abort(); - } -#else - = 0; -#endif + virtual bool static_equal(const Expr &that) const = 0; /// @param that an Expr object /// @note @c that is guaranteed to be of same type as @c *this, hence can be @@ -483,12 +432,6 @@ class Expr : public std::enable_shared_from_this { static type_id_type type_id = get_next_type_id(); return type_id; } - - private: - /// @input[in] fn the name of function that is missing in this class - /// @return an Exception object containing a message describing that @p - /// fn is missing from this type - Exception not_implemented(const char *fn) const; }; // class Expr static_assert(std::ranges::sized_range); diff --git a/SeQuant/core/expressions/expr_operators.hpp b/SeQuant/core/expressions/expr_operators.hpp index e1e0967e3c..2f83df592a 100644 --- a/SeQuant/core/expressions/expr_operators.hpp +++ b/SeQuant/core/expressions/expr_operators.hpp @@ -14,96 +14,9 @@ #include #include -#include namespace sequant { -inline bool operator==(const ExprPtr &left, const ExprPtr &right) { - return *left == *right; -} - -inline ExprPtr operator*(const ExprPtr &left, const ExprPtr &right) { - if (left.is() && right.is()) { - auto c_ = left->clone(); - auto &c = c_.as(); - c *= right.as(); - return c_; - } - - auto left_is_product = left->is(); - auto right_is_product = right->is(); - if (!left_is_product && !right_is_product) { - return ex(ExprPtrList{left, right}); - } else if (left_is_product) { - auto result = std::static_pointer_cast(left->clone()); - result->append(1, right); - return result; - } else { // right_is_product - auto result = std::static_pointer_cast(right->clone()); - result->prepend(1, left); - return result; - } - - SEQUANT_UNREACHABLE; -} - -/// Unlike @code operator*(const ExprPtr&, const ExprPtr&) @endcode this -/// produces a non-commutative product (i.e. NCProduct) -inline ExprPtr operator^(const ExprPtr &left, const ExprPtr &right) { - auto left_is_product = left->is(); - auto right_is_product = right->is(); - if (!left_is_product && !right_is_product) { - return ex(ExprPtrList{left, right}); - } else if (left_is_product) { - auto result = std::make_shared(left->clone().as()); - result->append(1, right); - return result; - } else { // right_is_product - auto result = std::make_shared(right->clone().as()); - result->prepend(1, left); - return result; - } - - SEQUANT_UNREACHABLE; -} - -inline ExprPtr operator+(const ExprPtr &left, const ExprPtr &right) { - auto left_is_sum = left->is(); - auto right_is_sum = right->is(); - if (!left_is_sum && !right_is_sum) { - return ex(ExprPtrList{left, right}); - } else if (left_is_sum) { - auto result = std::static_pointer_cast(left->clone()); - result->append(right); - return result; - } else { // right_is_sum - auto result = std::static_pointer_cast(right->clone()); - result->prepend(left); - return result; - } - - SEQUANT_UNREACHABLE; -} - -inline ExprPtr operator-(const ExprPtr &left, const ExprPtr &right) { - auto left_is_sum = left->is(); - if (!left_is_sum) { - return ex(ExprPtrList{ - left, - (right->is() ? ex(-right->as().value()) - : ex(-1, ExprPtrList{right}))}); - } else if (left_is_sum) { - auto result = std::static_pointer_cast(left->clone()); - if (right->is()) - result->append(ex(-right->as().value())); - else - result->append(ex(-1, ExprPtrList{right})); - return result; - } - - SEQUANT_UNREACHABLE; -} - template requires(std::constructible_from) ExprPtr operator+(const ExprPtr &lhs, T &&rhs) { diff --git a/SeQuant/core/expressions/expr_ptr.cpp b/SeQuant/core/expressions/expr_ptr.cpp new file mode 100644 index 0000000000..b1e21a0bb5 --- /dev/null +++ b/SeQuant/core/expressions/expr_ptr.cpp @@ -0,0 +1,190 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace sequant { + +ExprPtr ExprPtr::clone() const & { + if (!*this) return {}; + return ExprPtr(as_shared_ptr()->clone()); +} + +ExprPtr ExprPtr::clone() && noexcept { return std::move(*this); } + +ExprPtr::base_type &ExprPtr::as_shared_ptr() & { + return static_cast(*this); +} +const ExprPtr::base_type &ExprPtr::as_shared_ptr() const & { + return static_cast(*this); +} +ExprPtr::base_type &&ExprPtr::as_shared_ptr() && { + return static_cast(*this); +} + +Expr &ExprPtr::operator*() & { + SEQUANT_ASSERT(this->operator bool()); + return *(this->get()); +} + +const Expr &ExprPtr::operator*() const & { + SEQUANT_ASSERT(this->operator bool()); + return *(this->get()); +} + +Expr &&ExprPtr::operator*() && { + SEQUANT_ASSERT(this->operator bool()); + return std::move(*(this->get())); +} + +ExprPtr &ExprPtr::operator+=(const ExprPtr &other) { + if (!other) return *this; + + if (!*this) { + *this = other.clone(); + } else if (as_shared_ptr()->is()) { + as() += *other; + } else if (as_shared_ptr()->is() && other->is()) { + *this = ex(this->as().value() + + other->as().value()); + } else { + *this = ex(ExprPtrList{*this, other}); + } + return *this; +} + +ExprPtr &ExprPtr::operator-=(const ExprPtr &other) { + if (!other) return *this; + + if (!*this) { + *this = ex(-1) * other.clone(); + } else if (as_shared_ptr()->is()) { + as() -= *other; + } else if (as_shared_ptr()->is() && other->is()) { + *this = ex(this->as().value() - + other->as().value()); + } else { + *this = ex(ExprPtrList{*this, ex(-1, ExprPtrList{other})}); + } + return *this; +} + +ExprPtr &ExprPtr::operator*=(const ExprPtr &other) { + if (!other) return *this; + + if (!*this) { + *this = other.clone(); + } else if (as_shared_ptr()->is()) { + as() *= *other; + } else if (as_shared_ptr()->is() && other->is()) { + *this = ex(this->as().value() * + other->as().value()); + } else { + *this = ex(ExprPtrList{*this, other}); + } + return *this; +} + +std::size_t ExprPtr::size() const { return this->get()->size(); } + +std::wstring ExprPtr::to_latex() const { return as_shared_ptr()->to_latex(); } + +ExprPtr adjoint(const ExprPtr &expr) { + auto result = expr->clone(); + result->adjoint(); + return result; +} + +bool operator==(const ExprPtr &left, const ExprPtr &right) { + return *left == *right; +} + +ExprPtr operator*(const ExprPtr &left, const ExprPtr &right) { + if (left.is() && right.is()) { + auto c_ = left->clone(); + auto &c = c_.as(); + c *= right.as(); + return c_; + } + + auto left_is_product = left->is(); + auto right_is_product = right->is(); + if (!left_is_product && !right_is_product) { + return ex(ExprPtrList{left, right}); + } else if (left_is_product) { + auto result = std::static_pointer_cast(left->clone()); + result->append(1, right); + return result; + } else { // right_is_product + auto result = std::static_pointer_cast(right->clone()); + result->prepend(1, left); + return result; + } + + SEQUANT_UNREACHABLE; +} + +/// Unlike @code operator*(const ExprPtr&, const ExprPtr&) @endcode this +/// produces a non-commutative product (i.e. NCProduct) +ExprPtr operator^(const ExprPtr &left, const ExprPtr &right) { + auto left_is_product = left->is(); + auto right_is_product = right->is(); + if (!left_is_product && !right_is_product) { + return ex(ExprPtrList{left, right}); + } else if (left_is_product) { + auto result = std::make_shared(left->clone().as()); + result->append(1, right); + return result; + } else { // right_is_product + auto result = std::make_shared(right->clone().as()); + result->prepend(1, left); + return result; + } + + SEQUANT_UNREACHABLE; +} + +ExprPtr operator+(const ExprPtr &left, const ExprPtr &right) { + auto left_is_sum = left->is(); + auto right_is_sum = right->is(); + if (!left_is_sum && !right_is_sum) { + return ex(ExprPtrList{left, right}); + } else if (left_is_sum) { + auto result = std::static_pointer_cast(left->clone()); + result->append(right); + return result; + } else { // right_is_sum + auto result = std::static_pointer_cast(right->clone()); + result->prepend(left); + return result; + } + + SEQUANT_UNREACHABLE; +} + +ExprPtr operator-(const ExprPtr &left, const ExprPtr &right) { + auto left_is_sum = left->is(); + if (!left_is_sum) { + return ex(ExprPtrList{ + left, + (right->is() ? ex(-right->as().value()) + : ex(-1, ExprPtrList{right}))}); + } else if (left_is_sum) { + auto result = std::static_pointer_cast(left->clone()); + if (right->is()) + result->append(ex(-right->as().value())); + else + result->append(ex(-1, ExprPtrList{right})); + return result; + } + + SEQUANT_UNREACHABLE; +} + +} // namespace sequant diff --git a/SeQuant/core/expressions/expr_ptr.hpp b/SeQuant/core/expressions/expr_ptr.hpp index 5eb532d0d0..39b9ebbe40 100644 --- a/SeQuant/core/expressions/expr_ptr.hpp +++ b/SeQuant/core/expressions/expr_ptr.hpp @@ -157,6 +157,15 @@ using ExprPtrVector = container::svector; /// @return the adjoint of @p expr ExprPtr adjoint(const ExprPtr &expr); +ExprPtr operator*(const ExprPtr &left, const ExprPtr &right); + +/// Unlike @code operator*(const ExprPtr&, const ExprPtr&) @endcode this +/// produces a non-commutative product (i.e. NCProduct) +ExprPtr operator^(const ExprPtr &left, const ExprPtr &right); + +ExprPtr operator+(const ExprPtr &left, const ExprPtr &right); +ExprPtr operator-(const ExprPtr &left, const ExprPtr &right); + } // namespace sequant #endif // SEQUANT_EXPRESSIONS_EXPR_PTR_HPP diff --git a/SeQuant/core/expressions/power.cpp b/SeQuant/core/expressions/power.cpp new file mode 100644 index 0000000000..f3920dac24 --- /dev/null +++ b/SeQuant/core/expressions/power.cpp @@ -0,0 +1,190 @@ +#include +#include +#include +#include +#include +#include + +namespace sequant { + +Power::Power(ExprPtr base, exponent_type exponent) + : base_{}, exponent_{std::move(exponent)} { + SEQUANT_ASSERT(base); + SEQUANT_ASSERT(base->is() || base->is()); + // clone on construction so that external + // mutations of the input cannot invalidate our memoized hash + base_ = base->clone(); + // 0^n is defined only for n >= 0 (0^0 = 1 by convention) + SEQUANT_ASSERT(!base_->is() || !base_->as().is_zero() || + exponent_ >= 0); +} + +const ExprPtr& Power::base() const { return base_; } + +const Power::exponent_type& Power::exponent() const { return exponent_; } + +bool Power::conjugated() const { return conjugated_; } + +void Power::conjugate() { + conjugated_ = !conjugated_; + reset_hash_value(); +} + +bool Power::is_zero() const { + return exponent_ > 0 && base_->is() && + base_->as().is_zero(); +} + +void Power::flatten(ExprPtr& expr) { + if (!expr || !expr->is()) return; + const auto& pw = expr->as(); + + // b^1 = b and conjugate if needed + if (pw.exponent_ == 1) { + auto lifted = pw.base_->clone(); + if (pw.conjugated_) lifted->adjoint(); + expr = std::move(lifted); + return; + } + // b^0 = 1 for any base (the ctor rejects 0^(negative) + if (pw.exponent_ == 0) { + expr = ex(Constant::scalar_type{1}); + return; + } + if (!pw.base_->is()) return; + + using scalar_type = Constant::scalar_type; + const auto& base_val = pw.base_->as().value(); + + // 1^k = 1 for any rational k. + if (base_val == scalar_type{1}) { + expr = ex(scalar_type{1}); + return; + } + + // Both remaining fold cases share one shape — `rational base raised to + // an integer exponent` — so we normalize to that shape and run a single + // exp-by-squaring loop. Anything else is left untouched. + // + // Case A: integer exponent (any Constant base, real or complex). + // `base` is just `base_val`. + // Case B: half-integer exponent on a non-negative real rational base + // `p/q` with both `p` and `q` perfect squares. Then + // (p/q)^(m/2) = (sqrt(p)/sqrt(q))^m, + // so we replace `base` with `sqrt(p)/sqrt(q)` (still a rational) and + // keep `exp_int = m`. + + // initialize the base + scalar_type base{0}; + auto exp_nr = numerator(pw.exponent_); // numerator of exponent + + if (denominator(pw.exponent_) == 1) { + base = base_val; + } else if (denominator(pw.exponent_) == 2 && base_val.imag() == 0 && + base_val.real() >= 0) { + intmax_t p = numerator(base_val.real()); + intmax_t q = denominator(base_val.real()); // > 0 by Boost's convention, + // sign is with the numerator + + // check for perfect squares + intmax_t p_rem{0}, q_rem{0}; + intmax_t p_root = boost::multiprecision::sqrt(p, p_rem); + intmax_t q_root = boost::multiprecision::sqrt(q, q_rem); + // fold if p and q are perfect squares, else return + if (p_rem != 0 || q_rem != 0) return; + base = scalar_type{rational(p_root) / rational(q_root)}; + } else { + return; + } + + // Standard exp-by-squaring; for negative exponents we power the + // magnitude and invert at the end. + const bool negate = exp_nr < 0; + if (negate) exp_nr = -exp_nr; + scalar_type value{1}; + scalar_type b = base; + while (exp_nr > 0) { + if (exp_nr % 2 != 0) value *= b; + exp_nr /= 2; + if (exp_nr > 0) b *= b; + } + if (negate) value = scalar_type{1} / value; + + if (pw.conjugated_) value = conj(value); + expr = ex(std::move(value)); +} + +Expr::type_id_type Power::type_id() const { return get_type_id(); } + +bool Power::is_scalar() const { return true; } + +ExprPtr Power::clone() const { + auto cloned = ex(base_, exponent_); + if (conjugated_) cloned->as().conjugate(); + return cloned; +} + +void Power::adjoint() { conjugate(); } + +Power& Power::operator*=(const Expr& that) { + // b^e1 *= b^e2 -> b^(e1+e2) + if (that.is()) { + const auto& other = that.as(); + if (conjugated_ == other.conjugated_ && *base_ == *other.base_) { + exponent_ += other.exponent_; + reset_hash_value(); + return *this; + } + } + // (b^e)* *= b* -> (b^(e+1))* + else if (base_->is() && that.is()) { + // check effective conjugation of Variable in this and that, if valid + // operation iff they match + const auto& base_var = base_->as(); + const auto& that_var = that.as(); + if (base_var.label() == that_var.label() && + (base_var.conjugated() ^ conjugated_) == that_var.conjugated()) { + exponent_ += rational{1}; + reset_hash_value(); + return *this; + } + } + // C^e *= C -> C^(e+1) + else if (!conjugated_ && *base_ == that) { + exponent_ += rational{1}; + reset_hash_value(); + return *this; + } + throw Exception("Power::operator*=(that): not valid for that"); +} + +Expr::hash_type Power::memoizing_hash() const { + auto compute_hash = [this]() { + if (exponent_ == 1 && !conjugated_) return hash::value(*base_); + auto val = hash::value(*base_); + hash::combine(val, hash::value(exponent_)); + hash::combine(val, conjugated_); + return val; + }; + + if (!hash_value_) { + hash_value_ = compute_hash(); + } else { + SEQUANT_ASSERT(*hash_value_ == compute_hash()); + } + return *hash_value_; +} + +bool Power::static_equal(const Expr& that) const { + const auto& other = static_cast(that); + return exponent_ == other.exponent_ && conjugated_ == other.conjugated_ && + *base_ == *other.base_; +} + +bool Power::static_less_than(const Expr& that) const { + const auto& other = static_cast(that); + if (*base_ != *other.base_) return *base_ < *other.base_; + if (exponent_ != other.exponent_) return exponent_ < other.exponent_; + return conjugated_ < other.conjugated_; +} +} // namespace sequant diff --git a/SeQuant/core/expressions/power.hpp b/SeQuant/core/expressions/power.hpp index f5731dabf3..29fb69e505 100644 --- a/SeQuant/core/expressions/power.hpp +++ b/SeQuant/core/expressions/power.hpp @@ -5,10 +5,7 @@ #include #include #include -#include -#include #include -#include namespace sequant { @@ -27,17 +24,7 @@ class Power : public Expr { /// @param[in] base the base expression; must be a Constant or Variable. /// @param[in] exponent rational exponent - Power(ExprPtr base, exponent_type exponent) - : base_{}, exponent_{std::move(exponent)} { - SEQUANT_ASSERT(base); - SEQUANT_ASSERT(base->is() || base->is()); - // clone on construction so that external - // mutations of the input cannot invalidate our memoized hash - base_ = base->clone(); - // 0^n is defined only for n >= 0 (0^0 = 1 by convention) - SEQUANT_ASSERT(!base_->is() || !base_->as().is_zero() || - exponent_ >= 0); - } + Power(ExprPtr base, exponent_type exponent); /// @overload constructs a `Variable` base from @p label template @@ -55,29 +42,23 @@ class Power : public Expr { : Power(ex(std::forward(value)), std::move(exponent)) {} /// @return the base expression - const ExprPtr& base() const { return base_; } + const ExprPtr& base() const; /// @return the rational exponent - const exponent_type& exponent() const { return exponent_; } + const exponent_type& exponent() const; /// @return whether this Power has been complex-conjugated via adjoint() /// @note Conjugation is tracked as a flag because, in general, /// `conj(base^exponent) != conj(base)^exponent` - bool conjugated() const { return conjugated_; } + bool conjugated() const; /// @brief toggles the conjugation flag - void conjugate() { - conjugated_ = !conjugated_; - reset_hash_value(); - } + void conjugate(); /// @return true if the base is zero and the exponent is positive /// @note Construction rejects all undefined 0^n cases; 0^0 is legal and /// treated as 1. - bool is_zero() const override { - return exponent_ > 0 && base_->is() && - base_->as().is_zero(); - } + bool is_zero() const override; /// @brief Attempts to flatten a Power, mutating @p expr in place. Folds /// when @p expr holds a Power and any of: @@ -93,97 +74,16 @@ class Power : public Expr { /// @note Only square-root exponents are folded (that is the only /// case needed in practice right now). Extending to general n-th roots only /// requires replacing the integer-square-root step with an integer n-th-root. - static void flatten(ExprPtr& expr) { - if (!expr || !expr->is()) return; - const auto& pw = expr->as(); + static void flatten(ExprPtr& expr); - // b^1 = b and conjugate if needed - if (pw.exponent_ == 1) { - auto lifted = pw.base_->clone(); - if (pw.conjugated_) lifted->adjoint(); - expr = std::move(lifted); - return; - } - // b^0 = 1 for any base (the ctor rejects 0^(negative) - if (pw.exponent_ == 0) { - expr = ex(Constant::scalar_type{1}); - return; - } - if (!pw.base_->is()) return; + type_id_type type_id() const override; - using scalar_type = Constant::scalar_type; - const auto& base_val = pw.base_->as().value(); + bool is_scalar() const override; - // 1^k = 1 for any rational k. - if (base_val == scalar_type{1}) { - expr = ex(scalar_type{1}); - return; - } - - // Both remaining fold cases share one shape — `rational base raised to - // an integer exponent` — so we normalize to that shape and run a single - // exp-by-squaring loop. Anything else is left untouched. - // - // Case A: integer exponent (any Constant base, real or complex). - // `base` is just `base_val`. - // Case B: half-integer exponent on a non-negative real rational base - // `p/q` with both `p` and `q` perfect squares. Then - // (p/q)^(m/2) = (sqrt(p)/sqrt(q))^m, - // so we replace `base` with `sqrt(p)/sqrt(q)` (still a rational) and - // keep `exp_int = m`. - - // initialize the base - scalar_type base{0}; - auto exp_nr = numerator(pw.exponent_); // numerator of exponent - - if (denominator(pw.exponent_) == 1) { - base = base_val; - } else if (denominator(pw.exponent_) == 2 && base_val.imag() == 0 && - base_val.real() >= 0) { - intmax_t p = numerator(base_val.real()); - intmax_t q = denominator(base_val.real()); // > 0 by Boost's convention, - // sign is with the numerator - - // check for perfect squares - intmax_t p_rem{0}, q_rem{0}; - intmax_t p_root = boost::multiprecision::sqrt(p, p_rem); - intmax_t q_root = boost::multiprecision::sqrt(q, q_rem); - // fold if p and q are perfect squares, else return - if (p_rem != 0 || q_rem != 0) return; - base = scalar_type{rational(p_root) / rational(q_root)}; - } else { - return; - } - - // Standard exp-by-squaring; for negative exponents we power the - // magnitude and invert at the end. - const bool negate = exp_nr < 0; - if (negate) exp_nr = -exp_nr; - scalar_type value{1}; - scalar_type b = base; - while (exp_nr > 0) { - if (exp_nr % 2 != 0) value *= b; - exp_nr /= 2; - if (exp_nr > 0) b *= b; - } - if (negate) value = scalar_type{1} / value; - - if (pw.conjugated_) value = conj(value); - expr = ex(std::move(value)); - } - - type_id_type type_id() const override { return get_type_id(); } - - bool is_scalar() const override { return true; } - - ExprPtr clone() const override { - auto cloned = ex(base_, exponent_); - if (conjugated_) cloned->as().conjugate(); - return cloned; - } + ExprPtr clone() const override; /// @brief adjoint of Power: flips the conjugation flag. - void adjoint() override { conjugate(); } + void adjoint() override; /// @brief Combines exponents when effective bases match: /// - `b^e1 *= b^e2` → `b^(e1+e2)` when this and @p that share the same @@ -193,37 +93,7 @@ class Power : public Expr { /// and the effective conjugation parities align. For a Constant base /// only the fully unconjugated case combines. /// @throw Exception if @p that is not combinable. - Expr& operator*=(const Expr& that) override { - // b^e1 *= b^e2 -> b^(e1+e2) - if (that.is()) { - const auto& other = that.as(); - if (conjugated_ == other.conjugated_ && *base_ == *other.base_) { - exponent_ += other.exponent_; - reset_hash_value(); - return *this; - } - } - // (b^e)* *= b* -> (b^(e+1))* - else if (base_->is() && that.is()) { - // check effective conjugation of Variable in this and that, if valid - // operation iff they match - const auto& base_var = base_->as(); - const auto& that_var = that.as(); - if (base_var.label() == that_var.label() && - (base_var.conjugated() ^ conjugated_) == that_var.conjugated()) { - exponent_ += rational{1}; - reset_hash_value(); - return *this; - } - } - // C^e *= C -> C^(e+1) - else if (!conjugated_ && *base_ == that) { - exponent_ += rational{1}; - reset_hash_value(); - return *this; - } - throw Exception("Power::operator*=(that): not valid for that"); - } + Power& operator*=(const Expr& that); private: ExprPtr base_; @@ -232,35 +102,11 @@ class Power : public Expr { /// @return hash of this Power /// @note when exponent is 1 and not conjugated the hash matches the base's - hash_type memoizing_hash() const override { - auto compute_hash = [this]() { - if (exponent_ == 1 && !conjugated_) return hash::value(*base_); - auto val = hash::value(*base_); - hash::combine(val, hash::value(exponent_)); - hash::combine(val, conjugated_); - return val; - }; - - if (!hash_value_) { - hash_value_ = compute_hash(); - } else { - SEQUANT_ASSERT(*hash_value_ == compute_hash()); - } - return *hash_value_; - } + hash_type memoizing_hash() const override; - bool static_equal(const Expr& that) const override { - const auto& other = static_cast(that); - return exponent_ == other.exponent_ && conjugated_ == other.conjugated_ && - *base_ == *other.base_; - } + bool static_equal(const Expr& that) const override; - bool static_less_than(const Expr& that) const override { - const auto& other = static_cast(that); - if (*base_ != *other.base_) return *base_ < *other.base_; - if (exponent_ != other.exponent_) return exponent_ < other.exponent_; - return conjugated_ < other.conjugated_; - } + bool static_less_than(const Expr& that) const override; }; } // namespace sequant diff --git a/SeQuant/core/expressions/product.cpp b/SeQuant/core/expressions/product.cpp new file mode 100644 index 0000000000..48ac8a2684 --- /dev/null +++ b/SeQuant/core/expressions/product.cpp @@ -0,0 +1,400 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +namespace sequant { + +Product::Product(ExprPtrList factors, Flatten flatten_tag) { + using std::begin; + using std::end; + for (auto it = begin(factors); it != end(factors); ++it) + append(1, *it, flatten_tag); +} + +Product &Product::append(ExprPtr factor, Flatten flatten_tag) { + return this->append(1, factor, flatten_tag); +} + +const Product::scalar_type &Product::scalar() const { return scalar_; } + +bool Product::is_zero() const { return Constant::is_zero(this->scalar()); } + +const Product::factors_type &Product::factors() const { return factors_; } +Product::factors_type &Product::factors() { return factors_; } + +const ExprPtr &Product::factor(size_t i) const { return factors_.at(i); } + +bool Product::empty() const { return factors_.empty(); } + +bool Product::is_commutative() const { + bool result = true; + const auto nfactors = size(); + for (size_t f = 0; f != nfactors; ++f) { + for (size_t s = 1; result && s != nfactors; ++s) { + result &= factors_[f]->commutes_with(*factors_[s]); + } + } + return result; +} + +ExprPtr Product::canonicalize_impl(CanonicalizeOptions opts) { + // recursively canonicalize non-tensor subfactors (tensors will be + // canonicalized as part of the TN built of all tensor factors of this) ... + ranges::for_each(factors_, [this, opts](auto &factor) { + if (factor.template is()) { + return; + } + auto bp = factor->canonicalize(opts); + if (bp) { + SEQUANT_ASSERT(bp->template is()); + this->scalar_ *= std::static_pointer_cast(bp)->value(); + } + }); + + if (Logger::instance().canonicalize) { + std::wcout << "Product canonicalization(" << to_wstring(opts.method) + << ") input: " << to_latex() << std::endl; + } + + // pull out all scalar factors to the front + auto is_scalar = [](const auto &factor) { return factor->is_scalar(); }; + auto scalars = + factors_ | ranges::views::filter(is_scalar) | ranges::to_vector; + // scalars commute, so we can reorder them freely + ranges::sort(scalars, [](const auto &first, const auto &second) { + return *first < *second; + }); + + factors_ = factors_ | ranges::views::filter([&is_scalar](const auto &factor) { + return !is_scalar(factor); + }) | + ranges::to; + + // if there are no factors, insert scalars back and return + if (factors_.empty()) { + factors_.insert(factors_.begin(), scalars.begin(), scalars.end()); + return {}; + } + + auto contains_nontensors = ranges::any_of(factors_, [](const auto &factor) { + return std::dynamic_pointer_cast(factor) == nullptr; + }); + if (!contains_nontensors) { // tensor network canonization is a special case + // that's done in + // TensorNetwork + auto make_canonical_tn = [this, &opts](auto *tn_null_ptr) { + using TN = std::decay_t>; + ExprPtr canon_factor; + TN tn(this->factors_); + if constexpr (TN::version() == 3) { + canon_factor = tn.canonicalize( + TensorCanonicalizer::cardinal_tensor_labels(), opts); + } else { + using NamedIndexSet = tensor_network::NamedIndexSet; + std::shared_ptr named_indices = + !opts.named_indices + ? nullptr + : std::make_shared(opts.named_indices->begin(), + opts.named_indices->end()); + canon_factor = tn.canonicalize( + TensorCanonicalizer::cardinal_tensor_labels(), + opts.method == CanonicalizationMethod::Rapid, named_indices.get()); + } + return std::pair{std::move(tn), canon_factor}; + }; + using TN = TensorNetwork; + auto [tn, canon_factor] = make_canonical_tn(static_cast(nullptr)); + + const auto &tensors = tn.tensors(); + using std::size; + SEQUANT_ASSERT(size(tensors) == size(factors_)); + using std::begin; + using std::end; + std::transform(begin(tensors), end(tensors), begin(factors_), + [](const auto &tptr) { + auto exprptr = std::dynamic_pointer_cast(tptr); + SEQUANT_ASSERT(exprptr); + return exprptr; + }); + if (canon_factor) scalar_ *= canon_factor->template as().value(); + this->reset_hash_value(); + } else { // if contains non-tensors, do commutation-checking resort + + // comparer that respects cardinal tensor labels + auto &cardinal_tensor_labels = + TensorCanonicalizer::cardinal_tensor_labels(); + auto local_compare = [&cardinal_tensor_labels](const ExprPtr &first, + const ExprPtr &second) { + if (first->is() && second->is()) { + const auto first_label = first->as().label(); + const auto second_label = second->as().label(); + if (first_label == second_label) return *first < *second; + const auto first_is_cardinal_it = ranges::find_if( + cardinal_tensor_labels, + [&first_label](const std::wstring &l) { return l == first_label; }); + const auto first_is_cardinal = + first_is_cardinal_it != ranges::end(cardinal_tensor_labels); + const auto second_is_cardinal_it = ranges::find_if( + cardinal_tensor_labels, [&second_label](const std::wstring &l) { + return l == second_label; + }); + const auto second_is_cardinal = + second_is_cardinal_it != ranges::end(cardinal_tensor_labels); + if (first_is_cardinal && second_is_cardinal) + return first_is_cardinal_it < second_is_cardinal_it; + else if (first_is_cardinal && !second_is_cardinal) + return true; + else if (!first_is_cardinal && second_is_cardinal) + return false; + else { + SEQUANT_ASSERT(!first_is_cardinal && !second_is_cardinal); + return *first < *second; + } + } else + return *first < *second; + }; + + // ... then resort, respecting commutativity + using std::begin; + using std::end; + if (static_commutativity()) { + if (is_commutative()) { + std::stable_sort(begin(factors_), end(factors_), local_compare); + } + } else { + // must do bubble sort if not commuting to avoid swapping elements across + // a noncommuting element + bubble_sort( + begin(factors_), end(factors_), + [&local_compare](const ExprPtr &first, const ExprPtr &second) { + bool result = (first->commutes_with(*second)) + ? local_compare(first, second) + : false; + return result; + }); + } + } + // reinsert scalar factors at the front + factors_.insert(factors_.begin(), scalars.begin(), scalars.end()); + + // TODO evaluate product of Tensors (turn this into Products of Products) + + if (Logger::instance().canonicalize) + std::wcout << "Product canonicalization(" << to_wstring(opts.method) + << ") result: " << to_latex() << std::endl; + + return {}; // side effects are absorbed into the scalar_ +} + +void Product::adjoint() { + SEQUANT_ASSERT(static_commutativity() == false); // assert no slicing + auto adj_scalar = conj(scalar()); + using namespace ranges; + auto adj_factors = + factors() | views::reverse | + views::transform([](auto &expr) { return ::sequant::adjoint(expr); }); + using std::swap; + *this = + Product(adj_scalar, ranges::begin(adj_factors), ranges::end(adj_factors)); +} + +ExprPtr Product::canonicalize(CanonicalizeOptions opt) { + return this->canonicalize_impl(opt); +} + +ExprPtr Product::rapid_canonicalize(CanonicalizeOptions opt) { + SEQUANT_ASSERT(opt.method == CanonicalizationMethod::Rapid); + return this->canonicalize_impl(opt); +} + +bool Product::static_commutativity() const { return false; } + +std::wstring Product::to_latex() const { return to_latex(false); } + +std::wstring Product::to_latex(bool negate) const { + std::wstring result; + result = L"{"; + if (!scalar().is_zero()) { + const auto scal = negate ? -scalar() : scalar(); + if (!scal.is_identity()) { + // replace -1 prefactor by - + if (!(negate ? scalar() : -scalar()).is_identity()) { + result += io::latex::to_string(scal); + } else { + result += L"{-}"; + } + } + for (const auto &i : factors()) { + if (i->is()) + result += L"\\bigl(" + i->to_latex() + L"\\bigr)"; + else + result += i->to_latex(); + } + } + result += L"}"; + return result; +} + +Product::type_id_type Product::type_id() const { + return get_type_id(); +}; + +/// @return an identical clone of this Product (a deep copy allocated on the +/// heap) +/// @note this does not flatten the product +ExprPtr Product::clone() const { return ex(this->deep_copy()); } + +Product Product::deep_copy() const { + auto cloned_factors = + factors() | ranges::views::transform([](const ExprPtr &ptr) { + return ptr ? ptr->clone() : nullptr; + }); + Product result(this->scalar(), ExprPtrList{}); + ranges::for_each(cloned_factors, [&](const auto &cloned_factor) { + result.append(1, std::move(cloned_factor), Flatten::No); + }); + return result; +} + +Product &Product::operator*=(const Expr &that) { + if (!that.is()) { + this->append(1, const_cast(that).shared_from_this()); + } else { + scalar_ *= that.as().value(); + } + return *this; +} + +void Product::add_identical(const Product &other) { + SEQUANT_ASSERT(ranges::equal(this->factors(), other.factors())); + scalar_ += other.scalar_; +} + +void Product::add_identical(const std::shared_ptr &other) { + SEQUANT_ASSERT(ranges::equal(this->factors(), other->factors())); + scalar_ += other->scalar_; +} + +void Product::add_identical(const ExprPtr &other) { + if (other.is()) return this->add_identical(other.as()); + + // only makes sense if this has a single factor + SEQUANT_ASSERT(this->factors_.size() == 1 && this->factors_[0] == other); + scalar_ += 1; +} + +ExprIterator Product::begin_subexpr() { + if (!factors_.empty()) { + reset_hash_value(); + } + + return ExprIterator{factors_.data()}; +} + +ExprIterator Product::end_subexpr() { + return ExprIterator{factors_.data() + factors_.size()}; +} + +ConstExprIterator Product::begin_subexpr() const { + return ConstExprIterator{factors_.data()}; +} + +ConstExprIterator Product::end_subexpr() const { + return ConstExprIterator{factors_.data() + factors_.size()}; +} + +/// @return the hash of this object, by hashing only the factors, +/// not the scalar to make possible rapid finding of Products that only +/// differ by a factor +/// @note this ensures that hash of a Product involving a single factor is +/// identical to the hash of the factor itself. +Expr::hash_type Product::memoizing_hash() const { + auto compute_hash = [this]() { + if (factors().size() == 1) + return factors_[0]->hash_value(); + else { + auto deref_factors = + factors() | + ranges::views::transform( + [](const ExprPtr &ptr) -> const Expr & { return *ptr; }); + auto value = + hash::range(ranges::begin(deref_factors), ranges::end(deref_factors)); + return value; + } + }; + + if (!hash_value_) { + hash_value_ = compute_hash(); + } else { + SEQUANT_ASSERT(*hash_value_ == compute_hash()); + } + + return *hash_value_; +} + +bool Product::static_equal(const Expr &that) const { + const auto &that_cast = static_cast(that); + if (scalar() == that_cast.scalar() && + factors().size() == that_cast.factors().size()) { + if (this->empty()) return true; + // compare hash values first + if (this->hash_value() == + that.hash_value()) // hash values agree -> do full comparison + return std::equal(begin_subexpr(), end_subexpr(), that.begin_subexpr(), + expr_ptr_comparer); + else + return false; + } else + return false; +} + +CProduct::CProduct(const Product &other) : Product(other) {} +CProduct::CProduct(Product &&other) : Product(std::move(other)) {} + +bool CProduct::is_commutative() const { return true; } + +void CProduct::adjoint() { + auto adj_scalar = conj(scalar()); + using namespace ranges; + // no need to reverse for commutative product + auto adj_factors = factors() | views::transform([](auto &&expr) { + return ::sequant::adjoint(expr); + }); + *this = CProduct(adj_scalar, ranges::begin(adj_factors), + ranges::end(adj_factors)); +} + +bool CProduct::static_commutativity() const { return true; } + +NCProduct::NCProduct(const Product &other) : Product(other) {} +NCProduct::NCProduct(Product &&other) : Product(std::move(other)) {} + +bool NCProduct::is_commutative() const { return false; } + +void NCProduct::adjoint() { + auto adj_scalar = conj(scalar()); + using namespace ranges; + // no need to reverse for commutative product + auto adj_factors = + factors() | std::views::reverse | std::views::transform([](auto &&expr) { + return ::sequant::adjoint(expr); + }); + *this = NCProduct(adj_scalar, ranges::begin(adj_factors), + ranges::end(adj_factors)); +} + +bool NCProduct::static_commutativity() const { return true; } + +} // namespace sequant diff --git a/SeQuant/core/expressions/product.hpp b/SeQuant/core/expressions/product.hpp index 8d9e1a629a..1c12ed2728 100644 --- a/SeQuant/core/expressions/product.hpp +++ b/SeQuant/core/expressions/product.hpp @@ -4,17 +4,12 @@ #include #include #include -#include #include #include -#include #include #include -#include -#include #include -#include #include #include @@ -35,6 +30,7 @@ class Product : public Expr { enum class Flatten { Once, Recursively, Yes = Recursively, No }; using scalar_type = Constant::scalar_type; + using factors_type = container::svector; Product() = default; virtual ~Product() = default; @@ -46,12 +42,7 @@ class Product : public Expr { /// construct a Product out of zero or more factors (multiplied by 1) /// @param factors the factors /// @param flatten_tag if Flatten::Yes, flatten the factors - Product(ExprPtrList factors, Flatten flatten_tag = Flatten::Yes) { - using std::begin; - using std::end; - for (auto it = begin(factors); it != end(factors); ++it) - append(1, *it, flatten_tag); - } + Product(ExprPtrList factors, Flatten flatten_tag = Flatten::Yes); /// construct a Product out of zero or more factors (multiplied by 1) /// @param rng a range of factors; if rng is a Product, it will be flattened @@ -198,9 +189,7 @@ class Product : public Expr { /// @param factor a factor by which to multiply the product /// @param flatten_tag specifies whether (and how) to flatten the argument(s) /// @return @c *this - Product &append(ExprPtr factor, Flatten flatten_tag = Flatten::Yes) { - return this->append(1, factor, flatten_tag); - } + Product &append(ExprPtr factor, Flatten flatten_tag = Flatten::Yes); /// (post-)multiplies the product by @c factor /// @param factor a factor by which to multiply the product @@ -268,13 +257,13 @@ class Product : public Expr { flatten_tag); } - const auto &scalar() const { return scalar_; } + const scalar_type &scalar() const; /// @return `Constant::is_zero(this->scalar())` - bool is_zero() const override { return Constant::is_zero(this->scalar()); } + bool is_zero() const override; - const auto &factors() const { return factors_; } - auto &factors() { return factors_; } + const factors_type &factors() const; + factors_type &factors(); /// @brief View view of factors that are scalars (anything for which /// Expr::is_scalar() returns true). @@ -294,10 +283,10 @@ class Product : public Expr { /// Factor accessor /// @param i factor index /// @return ith factor - const ExprPtr &factor(size_t i) const { return factors_.at(i); } + const ExprPtr &factor(size_t i) const; /// @return true if the number of factors is zero - bool empty() const { return factors_.empty(); } + bool empty() const; /// @brief checks commutativity recursively /// @return true if definitely commutative, false definitely not commutative @@ -309,175 +298,73 @@ class Product : public Expr { /// factors, with complex-conjugated scalar virtual void adjoint() override; - private: - /// @return true if commutativity is decidable statically - /// @sa CProduct::static_commutativity() and NCProduct::static_commutativity() - virtual bool static_commutativity() const { return false; } - - public: - std::wstring to_latex() const override { return to_latex(false); } + std::wstring to_latex() const override; /// just like Expr::to_latex() , but can negate before conversion /// @param[in] negate if true, scalar will be before conversion - std::wstring to_latex(bool negate) const { - std::wstring result; - result = L"{"; - if (!scalar().is_zero()) { - const auto scal = negate ? -scalar() : scalar(); - if (!scal.is_identity()) { - // replace -1 prefactor by - - if (!(negate ? scalar() : -scalar()).is_identity()) { - result += io::latex::to_string(scal); - } else { - result += L"{-}"; - } - } - for (const auto &i : factors()) { - if (i->is()) - result += L"\\bigl(" + i->to_latex() + L"\\bigr)"; - else - result += i->to_latex(); - } - } - result += L"}"; - return result; - } + std::wstring to_latex(bool negate) const; - type_id_type type_id() const override { return get_type_id(); }; + type_id_type type_id() const override; /// @return an identical clone of this Product (a deep copy allocated on the /// heap) /// @note this does not flatten the product - ExprPtr clone() const override { return ex(this->deep_copy()); } - - Product deep_copy() const { - auto cloned_factors = - factors() | ranges::views::transform([](const ExprPtr &ptr) { - return ptr ? ptr->clone() : nullptr; - }); - Product result(this->scalar(), ExprPtrList{}); - ranges::for_each(cloned_factors, [&](const auto &cloned_factor) { - result.append(1, std::move(cloned_factor), Flatten::No); - }); - return result; - } + ExprPtr clone() const override; - virtual Expr &operator*=(const Expr &that) override { - if (!that.is()) { - this->append(1, const_cast(that).shared_from_this()); - } else { - scalar_ *= that.as().value(); - } - return *this; - } + Product deep_copy() const; - void add_identical(const Product &other) { - SEQUANT_ASSERT(ranges::equal(this->factors(), other.factors())); - scalar_ += other.scalar_; - } + Product &operator*=(const Expr &that); - void add_identical(const std::shared_ptr &other) { - SEQUANT_ASSERT(ranges::equal(this->factors(), other->factors())); - scalar_ += other->scalar_; - } + void add_identical(const Product &other); - void add_identical(const ExprPtr &other) { - if (other.is()) return this->add_identical(other.as()); + void add_identical(const std::shared_ptr &other); - // only makes sense if this has a single factor - SEQUANT_ASSERT(this->factors_.size() == 1 && this->factors_[0] == other); - scalar_ += 1; - } + void add_identical(const ExprPtr &other); - ExprIterator begin_subexpr() override { - if (!factors_.empty()) { - reset_hash_value(); - } + ExprIterator begin_subexpr() override; - return ExprIterator{factors_.data()}; - } + ExprIterator end_subexpr() override; - ExprIterator end_subexpr() override { - return ExprIterator{factors_.data() + factors_.size()}; - } + ConstExprIterator begin_subexpr() const override; - ConstExprIterator begin_subexpr() const override { - return ConstExprIterator{factors_.data()}; - } + ConstExprIterator end_subexpr() const override; - ConstExprIterator end_subexpr() const override { - return ConstExprIterator{factors_.data() + factors_.size()}; - } + virtual ExprPtr canonicalize( + CanonicalizeOptions opt = + CanonicalizeOptions::default_options()) override; + + virtual ExprPtr rapid_canonicalize( + CanonicalizeOptions opts = + CanonicalizeOptions::default_options().copy_and_set( + CanonicalizationMethod::Rapid)) override; private: scalar_type scalar_ = {1, 0}; - container::svector factors_{}; + factors_type factors_{}; + + /// @return true if commutativity is decidable statically + /// @sa CProduct::static_commutativity() and NCProduct::static_commutativity() + virtual bool static_commutativity() const; /// @return the hash of this object, by hashing only the factors, /// not the scalar to make possible rapid finding of Products that only /// differ by a factor /// @note this ensures that hash of a Product involving a single factor is /// identical to the hash of the factor itself. - hash_type memoizing_hash() const override { - auto compute_hash = [this]() { - if (factors().size() == 1) - return factors_[0]->hash_value(); - else { - auto deref_factors = - factors() | - ranges::views::transform( - [](const ExprPtr &ptr) -> const Expr & { return *ptr; }); - auto value = hash::range(ranges::begin(deref_factors), - ranges::end(deref_factors)); - return value; - } - }; - - if (!hash_value_) { - hash_value_ = compute_hash(); - } else { - SEQUANT_ASSERT(*hash_value_ == compute_hash()); - } - - return *hash_value_; - } + hash_type memoizing_hash() const override; ExprPtr canonicalize_impl(CanonicalizeOptions); - public: - virtual ExprPtr canonicalize( - CanonicalizeOptions opt = - CanonicalizeOptions::default_options()) override; - virtual ExprPtr rapid_canonicalize( - CanonicalizeOptions opts = - CanonicalizeOptions::default_options().copy_and_set( - CanonicalizationMethod::Rapid)) override; - - private: - bool static_equal(const Expr &that) const override { - const auto &that_cast = static_cast(that); - if (scalar() == that_cast.scalar() && - factors().size() == that_cast.factors().size()) { - if (this->empty()) return true; - // compare hash values first - if (this->hash_value() == - that.hash_value()) // hash values agree -> do full comparison - return std::equal(begin_subexpr(), end_subexpr(), that.begin_subexpr(), - expr_ptr_comparer); - else - return false; - } else - return false; - } + bool static_equal(const Expr &that) const override; }; // class Product class CProduct : public Product { public: using Product::Product; - CProduct(const Product &other) : Product(other) {} - CProduct(Product &&other) : Product(other) {} + CProduct(const Product &other); + CProduct(Product &&other); - bool is_commutative() const override { return true; } + bool is_commutative() const override; /// @brief adjoint of a CProduct is a product of adjoints of its factors, with /// complex-conjugated scalar @@ -485,23 +372,23 @@ class CProduct : public Product { virtual void adjoint() override; private: - bool static_commutativity() const override { return true; } + bool static_commutativity() const override; }; // class CProduct class NCProduct : public Product { public: using Product::Product; - NCProduct(const Product &other) : Product(other) {} - NCProduct(Product &&other) : Product(other) {} + NCProduct(const Product &other); + NCProduct(Product &&other); - bool is_commutative() const override { return false; } + bool is_commutative() const override; /// @brief adjoint of a NCProduct is a reserved product of adjoints of its /// factors, with complex-conjugated scalar virtual void adjoint() override; private: - bool static_commutativity() const override { return true; } + bool static_commutativity() const override; }; // class NCProduct } // namespace sequant diff --git a/SeQuant/core/expressions/sum.cpp b/SeQuant/core/expressions/sum.cpp new file mode 100644 index 0000000000..880c47d2dc --- /dev/null +++ b/SeQuant/core/expressions/sum.cpp @@ -0,0 +1,386 @@ +#include +#include +#include +#include +#include + +namespace sequant { + +Sum::Sum(ExprPtrList summands) { + // use append to flatten out Sum summands + for (auto &&summand : summands) { + append(std::forward(summand)); + } +} + +Sum::Sum(summands_type &&summands, move_only_tag) + : summands_(std::move(summands)) { + std::size_t pos = 0; + for (auto it = summands_.begin(); it != summands_.end(); ++it) { + auto &summand = *it; + bool do_erase = false; + if (summand->is_zero()) { + do_erase = true; + } else if (summand->is()) { + auto summand_constant = summand.as_shared_ptr(); + if (constant_summand_idx_) { // add up to the existing constant ... + SEQUANT_ASSERT(summands_.at(*constant_summand_idx_)->is()); + summands_[*constant_summand_idx_].as() += *summand_constant; + do_erase = true; + } else { // or memorize the position of the constant + constant_summand_idx_ = pos; + } + } + + // erase if needed + if (do_erase) { + summands_.erase(it); + it = summands_.begin(); + std::advance(it, pos); + } else + ++pos; + } +} + +Sum &Sum::append(ExprPtr summand) { + SEQUANT_ASSERT(summand); + if (!summand->is()) { + if (!summand->is_zero()) { // exclude zeros + if (summand->is()) { // add up constants + // immediately, if possible + auto summand_constant = summand.as_shared_ptr(); + if (constant_summand_idx_) { + SEQUANT_ASSERT(summands_.at(*constant_summand_idx_)->is()); + summands_[*constant_summand_idx_].as() += *summand; + } else { + summands_.push_back(summand->clone()); + constant_summand_idx_ = summands_.size() - 1; + } + } else { + summands_.push_back(summand->clone()); + } + reset_hash_value(); + } + } else { // this recursively flattens Sum summands + for (auto &subsummand : *summand) this->append(subsummand); + } + return *this; +} + +Sum &Sum::prepend(ExprPtr summand) { + SEQUANT_ASSERT(summand); + if (!summand->is()) { + if (!summand->is_zero()) { + // exclude zeros + if (summand->is()) { + auto summand_constant = summand.as_shared_ptr(); + if (constant_summand_idx_) { // add up to the existing constant ... + SEQUANT_ASSERT(summands_.at(*constant_summand_idx_)->is()); + summands_[*constant_summand_idx_].as() += *summand_constant; + } else { // or include the nonzero constant and update + // constant_summand_idx_ + summands_.insert(summands_.begin(), summand->clone()); + constant_summand_idx_ = 0; + } + } else { + summands_.insert(summands_.begin(), summand->clone()); + if (constant_summand_idx_) // if have a constant, update its position + ++*constant_summand_idx_; + } + reset_hash_value(); + } + } else { // this recursively flattens Sum summands + for (auto &subsummand : *summand) this->prepend(subsummand); + } + return *this; +} + +const Sum::summands_type &Sum::summands() const { return summands_; } + +const ExprPtr &Sum::summand(size_t i) const { return summands_.at(i); } + +ExprPtr Sum::take_n(size_t count) const { + const auto e = (count >= summands_.size() ? summands_.end() + : (summands_.begin() + count)); + return ex(summands_.begin(), e); +} + +ExprPtr Sum::take_n(size_t offset, size_t count) const { + const auto offset_plus_count = offset + count; + const auto b = (offset >= summands_.size() ? summands_.end() + : (summands_.begin() + offset)); + const auto e = (offset_plus_count >= summands_.size() + ? summands_.end() + : (summands_.begin() + offset_plus_count)); + return ex(b, e); +} + +bool Sum::empty() const { return summands_.empty(); } + +std::size_t Sum::size() const { return summands_.size(); } + +std::wstring Sum::to_latex() const { + std::wstring result; + result = L"{ \\bigl("; + std::size_t counter = 0; + for (const auto &i : summands()) { + const auto i_is_product = i->is(); + if (!i_is_product) { + result += (counter == 0) ? i->to_latex() : (L" + " + i->to_latex()); + } else { // i_is_product + const auto i_prod = i->as(); + const auto scalar = i_prod.scalar(); + if (scalar.real() < 0 || (scalar.real() == 0 && scalar.imag() < 0)) { + result += L" - " + i_prod.to_latex(true); + } else { + result += (counter == 0) ? i->to_latex() : (L" + " + i->to_latex()); + } + } + ++counter; + } + result += L"\\bigr) }"; + return result; +} + +Expr::type_id_type Sum::type_id() const { return Expr::get_type_id(); }; + +ExprPtr Sum::clone() const { + auto cloned_summands = + summands() | + ranges::views::transform([](const ExprPtr &ptr) { return ptr->clone(); }); + return ex(ranges::begin(cloned_summands), ranges::end(cloned_summands)); +} + +void Sum::adjoint() { + using namespace ranges; + auto adj_summands = summands() | views::transform([](auto &&expr) { + return ::sequant::adjoint(expr); + }); + *this = Sum(ranges::begin(adj_summands), ranges::end(adj_summands)); +} + +ExprPtr Sum::canonicalize_impl(bool multipass, CanonicalizeOptions opts) { + if (Logger::instance().canonicalize) + std::wcout << "Sum::canonicalize_impl: input = " + << to_latex_align(shared_from_this()) << std::endl; + + const auto npasses = multipass ? 2 : 1; + for (auto pass = 0; pass != npasses; ++pass) { + const auto rapid = (pass % 2 == 0); + + // canonicalizing TNs in a sum requires treating named indices as + // meaningful/distinct + auto opts_copy = opts; + opts_copy.ignore_named_index_labels = + CanonicalizeOptions::IgnoreNamedIndexLabel::No; + if (rapid) { + opts_copy.method = CanonicalizationMethod::Lexicographic; + } else + opts_copy.method = opts.method | CanonicalizationMethod::Topological; + + // recursively canonicalize summands ... + // using for_each and direct access to summands + sequant::for_each(summands_, [&opts_copy, &rapid](ExprPtr &summand) { + ExprPtr bp; + if (rapid) { + bp = summand->rapid_canonicalize(opts_copy); + } else { + bp = summand->canonicalize(opts_copy); + } + if (bp) { + SEQUANT_ASSERT(bp->template is()); + summand = ex(std::static_pointer_cast(bp)->value(), + ExprPtrList{summand}); + } + }); + if (Logger::instance().canonicalize) + std::wcout << "Sum::canonicalize_impl (pass=" << pass + << "): after canonicalizing summands = " + << to_latex_align(shared_from_this()) << std::endl; + + HashingAccumulator acc; + for (auto &summand : summands_) { + acc.append(summand); + } + + // last pass? sort by hash then by Expr::operator< + // N.B. no point in differentiating between canonicalization methods here + // since need to sort in both cases + auto new_sum = + (pass == npasses - 1) ? acc.make_canonicalized_sum() : acc.make_sum(); + using std::swap; + swap(*this, *new_sum); + + if (Logger::instance().canonicalize) + std::wcout << "Sum::canonicalize_impl (pass=" << pass + << "): after reducing summands = " + << to_latex_align(shared_from_this()) << std::endl; + } + + return {}; // side effects are absorbed into summands +} + +Sum &Sum::operator+=(const Expr &that) { + this->append(const_cast(that).shared_from_this()); + return *this; +} + +Sum &Sum::operator-=(const Expr &that) { + if (that.is()) + this->append(ex(-that.as().value())); + else + this->append(ex( + -1, ExprPtrList{const_cast(that).shared_from_this()})); + return *this; +} + +ExprIterator Sum::begin_subexpr() { + if (!summands_.empty()) { + reset_hash_value(); + } + + return ExprIterator{summands_.data()}; +} + +ExprIterator Sum::end_subexpr() { + return ExprIterator{summands_.data() + summands_.size()}; +} + +ConstExprIterator Sum::begin_subexpr() const { + return ConstExprIterator{summands_.data()}; +} + +ConstExprIterator Sum::end_subexpr() const { + return ConstExprIterator{summands_.data() + summands_.size()}; +} + +Expr::hash_type Sum::memoizing_hash() const { + auto compute_hash = [this]() { + if (summands_.size() == 1) + return summands_[0]->hash_value(); + else { + auto deref_summands = + summands() | + ranges::views::transform( + [](const ExprPtr &ptr) -> const Expr & { return *ptr; }); + auto value = hash::range(ranges::begin(deref_summands), + ranges::end(deref_summands)); + return value; + } + }; + + if (!hash_value_) { + hash_value_ = compute_hash(); + } else { + SEQUANT_ASSERT(*hash_value_ == compute_hash()); + } + + return *hash_value_; +} + +ExprPtr Sum::canonicalize(CanonicalizeOptions opt) { + return canonicalize_impl(true, opt); +} +ExprPtr Sum::rapid_canonicalize(CanonicalizeOptions opts) { + SEQUANT_ASSERT(opts.method == CanonicalizationMethod::Rapid); + return canonicalize_impl(false, opts); +} + +bool Sum::static_equal(const Expr &that) const { + const auto &that_cast = static_cast(that); + if (summands().size() == that_cast.summands().size()) { + if (this->empty()) return true; + // compare hash values first + if (this->hash_value() == + that.hash_value()) // hash values agree -> do full comparison + return std::equal(begin_subexpr(), end_subexpr(), that.begin_subexpr(), + expr_ptr_comparer); + else + return false; + } else + return false; +} + +HashingAccumulator &HashingAccumulator::append(ExprPtr summand, bool flatten) { + // flatten, if needed + if (flatten && summand.is()) { + for (auto &subsummand : summand.as().summands()) { + this->append(subsummand, flatten); + } + return *this; + } + + // process summand as a whole + auto it = summands_.find(summand); + if (it == summands_.end()) { + summands_.emplace(summand); + } else { // found existing term with the same hash + auto existing_summand = *it; + if (summand.template is()) { + if (existing_summand.is()) { + // both are products - add them + existing_summand.as().add_identical( + summand.template as()); + } else { + // convert existing term to product and add + auto product_copy = std::make_shared(summand->clone()); + product_copy->add_identical(existing_summand); + summands_.erase(it); + summands_.emplace(std::move(product_copy)); + } + } else { + if (existing_summand.is()) { + existing_summand.as().add_identical(summand); + } else { + // neither is a product - create new product + auto product_form = std::make_shared(); + product_form->append(2, summand.template as()); + summands_.erase(it); + summands_.emplace(std::move(product_form)); + } + } + } + + return *this; +} + +SumPtr HashingAccumulator::make_sum_impl(bool canonicalize) { + Sum::summands_type summands; + summands.reserve(summands_.size()); + for (auto summand : summands_) { + if (!summand->is_zero()) { + summands.push_back(summand); + } + } + + if (canonicalize) { + ranges::sort(summands, [](const auto &e1, const auto &e2) { + if (e1->hash_value() == e2->hash_value()) { + return e1 < e2; + } else { + return e1->hash_value() < e2->hash_value(); + } + }); + } + + return std::make_shared(std::move(summands), Sum::move_only_tag{}); +} + +SumPtr HashingAccumulator::make_sum() { return make_sum_impl(false); } + +SumPtr HashingAccumulator::make_canonicalized_sum() { + return make_sum_impl(true); +} + +ExprPtr HashingAccumulator::make_expr(bool canonicalize) { + if (summands_.size() == 0) { + return ex(0); + } else if (summands_.size() == 1) + return *(summands_.begin()); + else + return make_sum_impl(canonicalize); +} + +bool HashingAccumulator::empty() const { return summands_.empty(); } + +} // namespace sequant diff --git a/SeQuant/core/expressions/sum.hpp b/SeQuant/core/expressions/sum.hpp index 4cfcda2c8b..a0f3186c45 100644 --- a/SeQuant/core/expressions/sum.hpp +++ b/SeQuant/core/expressions/sum.hpp @@ -35,20 +35,9 @@ class Sum : public Expr { Sum &operator=(const Sum &) = default; Sum &operator=(Sum &&) = default; - void swap(Sum &other) { - Sum tmp = std::move(other); - other = std::move(*this); - *this = std::move(tmp); - } - /// construct a Sum out of zero or more summands /// @param summands an initializer list of summands - Sum(ExprPtrList summands) { - // use append to flatten out Sum summands - for (auto &&summand : summands) { - append(std::forward(summand)); - } - } + Sum(ExprPtrList summands); /// construct a Sum out of a range of summands /// @param begin the begin iterator @@ -93,120 +82,30 @@ class Sum : public Expr { /// construct a Sum by moving in the summands, no flattening is performed, /// but zeros will be omitted and constants added up /// @param summands the summands to move in - explicit Sum(summands_type &&summands, move_only_tag) - : summands_(std::move(summands)) { - std::size_t pos = 0; - for (auto it = summands_.begin(); it != summands_.end(); ++it) { - auto &summand = *it; - bool do_erase = false; - if (summand->is_zero()) { - do_erase = true; - } else if (summand->is()) { - auto summand_constant = summand.as_shared_ptr(); - if (constant_summand_idx_) { // add up to the existing constant ... - SEQUANT_ASSERT(summands_.at(*constant_summand_idx_)->is()); - *summands_[*constant_summand_idx_] += *summand_constant; - do_erase = true; - } else { // or memorize the position of the constant - constant_summand_idx_ = pos; - } - } - - // erase if needed - if (do_erase) { - summands_.erase(it); - it = summands_.begin(); - std::advance(it, pos); - } else - ++pos; - } - } + explicit Sum(summands_type &&summands, move_only_tag); /// append a summand to the sum /// @param summand the summand - Sum &append(ExprPtr summand) { - SEQUANT_ASSERT(summand); - if (!summand->is()) { - if (!summand->is_zero()) { // exclude zeros - if (summand->is()) { // add up constants - // immediately, if possible - auto summand_constant = summand.as_shared_ptr(); - if (constant_summand_idx_) { - SEQUANT_ASSERT( - summands_.at(*constant_summand_idx_)->is()); - *(summands_[*constant_summand_idx_]) += *summand; - } else { - summands_.push_back(summand->clone()); - constant_summand_idx_ = summands_.size() - 1; - } - } else { - summands_.push_back(summand->clone()); - } - reset_hash_value(); - } - } else { // this recursively flattens Sum summands - for (auto &subsummand : *summand) this->append(subsummand); - } - return *this; - } + Sum &append(ExprPtr summand); /// prepend a summand to the sum /// @param summand the summand - Sum &prepend(ExprPtr summand) { - SEQUANT_ASSERT(summand); - if (!summand->is()) { - if (!summand->is_zero()) { - // exclude zeros - if (summand->is()) { - auto summand_constant = summand.as_shared_ptr(); - if (constant_summand_idx_) { // add up to the existing constant ... - SEQUANT_ASSERT( - summands_.at(*constant_summand_idx_)->is()); - *summands_[*constant_summand_idx_] += *summand_constant; - } else { // or include the nonzero constant and update - // constant_summand_idx_ - summands_.insert(summands_.begin(), summand->clone()); - constant_summand_idx_ = 0; - } - } else { - summands_.insert(summands_.begin(), summand->clone()); - if (constant_summand_idx_) // if have a constant, update its position - ++*constant_summand_idx_; - } - reset_hash_value(); - } - } else { // this recursively flattens Sum summands - for (auto &subsummand : *summand) this->prepend(subsummand); - } - return *this; - } + Sum &prepend(ExprPtr summand); /// Summands accessor - const auto &summands() const { return summands_; } + const summands_type &summands() const; /// Summand accessor /// @param i summand index /// @return ith summand - const ExprPtr &summand(size_t i) const { return summands_.at(i); } + const ExprPtr &summand(size_t i) const; /// Takes the first @c count elements of the sum - ExprPtr take_n(size_t count) const { - const auto e = (count >= summands_.size() ? summands_.end() - : (summands_.begin() + count)); - return ex(summands_.begin(), e); - } + ExprPtr take_n(size_t count) const; /// Takes the first @c count elements of the sum starting with element @c /// offset - ExprPtr take_n(size_t offset, size_t count) const { - const auto offset_plus_count = offset + count; - const auto b = (offset >= summands_.size() ? summands_.end() - : (summands_.begin() + offset)); - const auto e = (offset_plus_count >= summands_.size() - ? summands_.end() - : (summands_.begin() + offset_plus_count)); - return ex(b, e); - } + ExprPtr take_n(size_t offset, size_t count) const; /// @tparam Filter a boolean predicate type, such `Filter(const ExprPtr&)` /// evaluates to true @@ -218,82 +117,31 @@ class Sum : public Expr { } /// @return true if the number of factors is zero - bool empty() const { return summands_.empty(); } + bool empty() const; /// @return the number of summands in a Sum - std::size_t size() const { return summands_.size(); } - - std::wstring to_latex() const override { - std::wstring result; - result = L"{ \\bigl("; - std::size_t counter = 0; - for (const auto &i : summands()) { - const auto i_is_product = i->is(); - if (!i_is_product) { - result += (counter == 0) ? i->to_latex() : (L" + " + i->to_latex()); - } else { // i_is_product - const auto i_prod = i->as(); - const auto scalar = i_prod.scalar(); - if (scalar.real() < 0 || (scalar.real() == 0 && scalar.imag() < 0)) { - result += L" - " + i_prod.to_latex(true); - } else { - result += (counter == 0) ? i->to_latex() : (L" + " + i->to_latex()); - } - } - ++counter; - } - result += L"\\bigr) }"; - return result; - } + std::size_t size() const; - Expr::type_id_type type_id() const override { - return Expr::get_type_id(); - }; + std::wstring to_latex() const override; - ExprPtr clone() const override { - auto cloned_summands = - summands() | ranges::views::transform( - [](const ExprPtr &ptr) { return ptr->clone(); }); - return ex(ranges::begin(cloned_summands), - ranges::end(cloned_summands)); - } + Expr::type_id_type type_id() const override; + + ExprPtr clone() const override; /// @brief adjoint of a Sum is a sum of adjoints of its factors virtual void adjoint() override; - virtual Expr &operator+=(const Expr &that) override { - this->append(const_cast(that).shared_from_this()); - return *this; - } + Sum &operator+=(const Expr &that); - virtual Expr &operator-=(const Expr &that) override { - if (that.is()) - this->append(ex(-that.as().value())); - else - this->append(ex( - -1, ExprPtrList{const_cast(that).shared_from_this()})); - return *this; - } + Sum &operator-=(const Expr &that); - ExprIterator begin_subexpr() override { - if (!summands_.empty()) { - reset_hash_value(); - } + ExprIterator begin_subexpr() override; - return ExprIterator{summands_.data()}; - } + ExprIterator end_subexpr() override; - ExprIterator end_subexpr() override { - return ExprIterator{summands_.data() + summands_.size()}; - } - - ConstExprIterator begin_subexpr() const override { - return ConstExprIterator{summands_.data()}; - } + ConstExprIterator begin_subexpr() const override; - ConstExprIterator end_subexpr() const override { - return ConstExprIterator{summands_.data() + summands_.size()}; - } + ConstExprIterator end_subexpr() const override; private: summands_type summands_{}; @@ -304,61 +152,21 @@ class Sum : public Expr { /// @return the hash of this object /// @note this ensures that hash of a Sum of a single summand is /// identical to the hash of the summand itself. - hash_type memoizing_hash() const override { - auto compute_hash = [this]() { - if (summands_.size() == 1) - return summands_[0]->hash_value(); - else { - auto deref_summands = - summands() | - ranges::views::transform( - [](const ExprPtr &ptr) -> const Expr & { return *ptr; }); - auto value = hash::range(ranges::begin(deref_summands), - ranges::end(deref_summands)); - return value; - } - }; - - if (!hash_value_) { - hash_value_ = compute_hash(); - } else { - SEQUANT_ASSERT(*hash_value_ == compute_hash()); - } - - return *hash_value_; - } + hash_type memoizing_hash() const override; /// @param multipass if true, will do a multipass canonicalization, with extra /// cleanup pass after the deep canonization pass ExprPtr canonicalize_impl(bool multipass, CanonicalizeOptions opt); - virtual ExprPtr canonicalize( - CanonicalizeOptions opt = - CanonicalizeOptions::default_options()) override { - return canonicalize_impl(true, opt); - } - virtual ExprPtr rapid_canonicalize( + ExprPtr canonicalize(CanonicalizeOptions opt = + CanonicalizeOptions::default_options()) override; + + ExprPtr rapid_canonicalize( CanonicalizeOptions opts = CanonicalizeOptions::default_options().copy_and_set( - CanonicalizationMethod::Rapid)) override { - SEQUANT_ASSERT(opts.method == CanonicalizationMethod::Rapid); - return canonicalize_impl(false, opts); - } + CanonicalizationMethod::Rapid)) override; - bool static_equal(const Expr &that) const override { - const auto &that_cast = static_cast(that); - if (summands().size() == that_cast.summands().size()) { - if (this->empty()) return true; - // compare hash values first - if (this->hash_value() == - that.hash_value()) // hash values agree -> do full comparison - return std::equal(begin_subexpr(), end_subexpr(), that.begin_subexpr(), - expr_ptr_comparer); - else - return false; - } else - return false; - } + bool static_equal(const Expr &that) const override; }; // class Sum /// @brief utility for eagerly accumulating summands in a hash table @@ -380,7 +188,7 @@ class HashingAccumulator { /// zero summands), or the lone summand itself ExprPtr make_expr(bool canonicalize = true); - bool empty() const { return summands_.empty(); } + bool empty() const; private: /// @brief Common implementation for make_sum and make_canonicalized_sum diff --git a/SeQuant/core/expressions/variable.cpp b/SeQuant/core/expressions/variable.cpp new file mode 100644 index 0000000000..58b15d00e9 --- /dev/null +++ b/SeQuant/core/expressions/variable.cpp @@ -0,0 +1,67 @@ +#include +#include +#include +#include +#include + +#include +#include + +namespace sequant { + +Variable::Variable(std::wstring label) + : label_(std::move(label)), conjugated_(false) {} + +Variable::Variable(const std::string &label) + : label_(sequant::toUtf16(label)), conjugated_(false) {} + +Expr::type_id_type Variable::type_id() const { return get_type_id(); } + +bool Variable::is_scalar() const { return true; } + +Expr::hash_type Variable::memoizing_hash() const { + auto compute_hash = [this]() { + auto val = hash::value(label_); + hash::combine(val, conjugated_); + return val; + }; + + if (!hash_value_) { + hash_value_ = compute_hash(); + } else { + SEQUANT_ASSERT(*hash_value_ == compute_hash()); + } + + return *hash_value_; +} + +bool Variable::static_equal(const Expr &that) const { + return label_ == static_cast(that).label_ && + conjugated_ == static_cast(that).conjugated_; +} + +std::wstring_view Variable::label() const { return label_; } + +void Variable::set_label(std::wstring label) { + label_ = std::move(label); + reset_hash_value(); +} + +void Variable::conjugate() { + conjugated_ = !conjugated_; + reset_hash_value(); +} + +bool Variable::conjugated() const { return conjugated_; } + +std::wstring Variable::to_latex() const { + std::wstring result = L"{" + io::latex::utf_to_string(label_) + L"}"; + if (conjugated_) result = L"{" + result + L"^*" + L"}"; + return result; +} + +ExprPtr Variable::clone() const { return ex(*this); } + +void Variable::adjoint() { conjugate(); } + +} // namespace sequant diff --git a/SeQuant/core/expressions/variable.hpp b/SeQuant/core/expressions/variable.hpp index 650ce64a39..4943845ac4 100644 --- a/SeQuant/core/expressions/variable.hpp +++ b/SeQuant/core/expressions/variable.hpp @@ -2,17 +2,16 @@ #define SEQUANT_EXPRESSIONS_VARIABLE_HPP #include -#include #include -#include #include -#include #include #include namespace sequant { +class ExprPtr; + /// This is represented as a "run-time" complex rational number class Variable : public Expr, public MutatableLabeled { public: @@ -29,10 +28,9 @@ class Variable : public Expr, public MutatableLabeled { std::constructible_from) explicit Variable(U &&label) : label_(std::forward(label)) {} - Variable(std::wstring label) : label_(std::move(label)), conjugated_(false) {} + Variable(std::wstring label); - Variable(const std::string &label) - : label_(sequant::toUtf16(label)), conjugated_(false) {} + Variable(const std::string &label); /// @return variable label /// @warning conjugation does not change it @@ -48,9 +46,9 @@ class Variable : public Expr, public MutatableLabeled { std::wstring to_latex() const override; - type_id_type type_id() const override { return get_type_id(); } + type_id_type type_id() const override; - bool is_scalar() const override { return true; } + bool is_scalar() const override; ExprPtr clone() const override; @@ -61,26 +59,10 @@ class Variable : public Expr, public MutatableLabeled { std::wstring label_; bool conjugated_ = false; - hash_type memoizing_hash() const override { - auto compute_hash = [this]() { - auto val = hash::value(label_); - hash::combine(val, conjugated_); - return val; - }; - - if (!hash_value_) { - hash_value_ = compute_hash(); - } else { - SEQUANT_ASSERT(*hash_value_ == compute_hash()); - } - - return *hash_value_; - } - - bool static_equal(const Expr &that) const override { - return label_ == static_cast(that).label_ && - conjugated_ == static_cast(that).conjugated_; - } + hash_type memoizing_hash() const override; + + bool static_equal(const Expr &that) const override; + }; // class Variable } // namespace sequant diff --git a/SeQuant/core/op.hpp b/SeQuant/core/op.hpp index 5225429122..4cc47e993d 100644 --- a/SeQuant/core/op.hpp +++ b/SeQuant/core/op.hpp @@ -955,6 +955,13 @@ class NormalOperator : public Operator, } }; +template <> +const container::svector & +NormalOperator::labels(); +template <> +const container::svector & +NormalOperator::labels(); + static_assert( is_tensor>, "The NormalOperator class does not fulfill the " @@ -1061,6 +1068,8 @@ class NormalOperatorSequence : public container::svector>, return Expr::get_type_id(); }; + ExprPtr clone() const override { return ex(*this); } + friend bool operator==(const NormalOperatorSequence &nopseq1, const NormalOperatorSequence &nopseq2) { return nopseq1.vacuum() == nopseq2.vacuum() && diff --git a/SeQuant/core/tensor_network/v1.cpp b/SeQuant/core/tensor_network/v1.cpp index cee3abce17..f351250420 100644 --- a/SeQuant/core/tensor_network/v1.cpp +++ b/SeQuant/core/tensor_network/v1.cpp @@ -458,7 +458,7 @@ ExprPtr TensorNetworkV1::canonicalize( nondefault_canonizer_ptr ? nondefault_canonizer_ptr.get() : &default_tensor_canonizer; auto bp = tensor_canonizer->apply(*tensor); - if (bp) *canon_byproduct *= *bp; + if (bp) canon_byproduct.as() *= *bp; } } edges_.clear(); diff --git a/SeQuant/core/wick.impl.hpp b/SeQuant/core/wick.impl.hpp index 89042c6f1f..f57579ee13 100644 --- a/SeQuant/core/wick.impl.hpp +++ b/SeQuant/core/wick.impl.hpp @@ -831,7 +831,7 @@ ExprPtr WickTheorem::compute(const bool count_only, nopseq->push_back(factor->template as>()); } else { SEQUANT_ASSERT(factor->is_cnumber()); - *prefactor *= *factor; + prefactor.as() *= *factor; } } init_input(nopseq); diff --git a/tests/unit/test_expr.cpp b/tests/unit/test_expr.cpp index b9112c7ee6..e9a1e0924d 100644 --- a/tests/unit/test_expr.cpp +++ b/tests/unit/test_expr.cpp @@ -38,6 +38,7 @@ struct Dummy : public sequant::Expr { std::wstring to_latex() const override { return L"{\\text{Dummy}}"; } type_id_type type_id() const override { return get_type_id(); }; sequant::ExprPtr clone() const override { return sequant::ex(); } + void adjoint() override {} bool static_equal(const sequant::Expr &) const override { return true; } }; @@ -69,6 +70,8 @@ struct VecExpr : public std::vector, public sequant::Expr { type_id_type type_id() const override { return get_type_id>(); }; + void adjoint() override {} + sequant::ConstExprIterator begin_subexpr() const override { if constexpr (sequant::Expr::is_shared_ptr_of_expr::value) { return sequant::ConstExprIterator{base_type::data()}; @@ -413,10 +416,6 @@ TEST_CASE("expr", "[elements]") { } SECTION("adjoint") { - { // not implemented by default - const auto e = std::make_shared(); - REQUIRE_THROWS_AS(e->adjoint(), Exception); - } { // implemented in Adjointable const auto e = std::make_shared(); REQUIRE_NOTHROW(e->adjoint());