diff --git a/CMakeLists.txt b/CMakeLists.txt index b04b376fc9..1083c18082 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -306,16 +306,30 @@ 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_container.cpp + SeQuant/core/expressions/expr_container.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/export/export.hpp b/SeQuant/core/export/export.hpp index 4b3627218a..86f627fecd 100644 --- a/SeQuant/core/export/export.hpp +++ b/SeQuant/core/export/export.hpp @@ -539,10 +539,10 @@ void track_usage(const EvalNode &node, PreprocessResult &result) { handle_variable(expr.as()); } else if (expr.is()) { const Power &power = expr.as(); - if (power.base().is()) { - handle_tensor(power.base().as()); - } else if (power.base().is()) { - handle_variable(power.base().as()); + if (power.base()->is()) { + handle_tensor(power.base()->as()); + } else if (power.base()->is()) { + handle_variable(power.base()->as()); } } } @@ -762,10 +762,10 @@ class PreprocessVisitor { handle_variable(node.left()->as_variable()); } else if (node.left()->is_power()) { const Power &power = node.left()->as_power(); - if (power.base().is()) { - handle_tensor(power.base().as()); - } else if (power.base().is()) { - handle_variable(power.base().as()); + if (power.base()->is()) { + handle_tensor(power.base()->as()); + } else if (power.base()->is()) { + handle_variable(power.base()->as()); } } @@ -775,10 +775,10 @@ class PreprocessVisitor { handle_variable(node.right()->as_variable()); } else if (node.right()->is_power()) { const Power &power = node.right()->as_power(); - if (power.base().is()) { - handle_tensor(power.base().as()); - } else if (power.base().is()) { - handle_variable(power.base().as()); + if (power.base()->is()) { + handle_tensor(power.base()->as()); + } else if (power.base()->is()) { + handle_variable(power.base()->as()); } } } diff --git a/SeQuant/core/export/itf.hpp b/SeQuant/core/export/itf.hpp index 3c1d0803fb..a7db26cb07 100644 --- a/SeQuant/core/export/itf.hpp +++ b/SeQuant/core/export/itf.hpp @@ -215,7 +215,7 @@ class ItfGenerator : public Generator { } std::string represent(const Power &power, const Context &ctx) const override { - const ExprPtr &base = power.base(); + const ExprContainer &base = power.base(); // ITF can only express powers of Constants if (!base->is()) { throw Exception( diff --git a/SeQuant/core/export/julia_tensor_operations.hpp b/SeQuant/core/export/julia_tensor_operations.hpp index 70342fa2dc..8140ee2c66 100644 --- a/SeQuant/core/export/julia_tensor_operations.hpp +++ b/SeQuant/core/export/julia_tensor_operations.hpp @@ -130,7 +130,7 @@ class JuliaTensorOperationsGenerator : public Generator { } std::string represent(const Power &power, const Context &ctx) const override { - const ExprPtr &base = power.base(); + const ExprContainer &base = power.base(); std::string base_str = to_julia_expr(*base, ctx); if (base->is() && base->as().conjugated()) { base_str = wrap_conj(std::move(base_str)); diff --git a/SeQuant/core/export/python_einsum.hpp b/SeQuant/core/export/python_einsum.hpp index ca1a6e8501..49cf3f84f3 100644 --- a/SeQuant/core/export/python_einsum.hpp +++ b/SeQuant/core/export/python_einsum.hpp @@ -198,7 +198,7 @@ class PythonEinsumGeneratorBase : public Generator { } std::string represent(const Power &power, const Context &ctx) const override { - const ExprPtr &base = power.base(); + const ExprContainer &base = power.base(); std::string base_str = stringify_scalar(*base, ctx); if (base->is() && base->as().conjugated()) { base_str = wrap_conj(std::move(base_str)); diff --git a/SeQuant/core/export/text_generator.hpp b/SeQuant/core/export/text_generator.hpp index 0a478fc5dd..c99d9f3d64 100644 --- a/SeQuant/core/export/text_generator.hpp +++ b/SeQuant/core/export/text_generator.hpp @@ -121,7 +121,7 @@ class TextGenerator : public Generator { } std::string represent(const Power &power, const Context &ctx) const override { - const ExprPtr &base = power.base(); + const ExprContainer &base = power.base(); std::string base_str = stringify(*base, ctx); if (base->is() && base->as().conjugated()) { base_str = wrap_conj(std::move(base_str)); diff --git a/SeQuant/core/export/utils.cpp b/SeQuant/core/export/utils.cpp index 182cc1977b..fd0841e6cc 100644 --- a/SeQuant/core/export/utils.cpp +++ b/SeQuant/core/export/utils.cpp @@ -4,8 +4,7 @@ #include -#include -#include +#include #include #include @@ -30,7 +29,7 @@ std::string format_power_exponent(const Power::exponent_type &exponent, return ss.str(); } -std::string format_power_base(const ExprPtr &base, std::string base_str) { +std::string format_power_base(const ExprContainer &base, std::string base_str) { if (base->is()) { const auto &v = base->as().value(); if (v.imag() == 0 && diff --git a/SeQuant/core/export/utils.hpp b/SeQuant/core/export/utils.hpp index 14798f42b5..103f217ac0 100644 --- a/SeQuant/core/export/utils.hpp +++ b/SeQuant/core/export/utils.hpp @@ -5,9 +5,7 @@ #ifndef SEQUANT_CORE_EXPORT_UTILS_HPP #define SEQUANT_CORE_EXPORT_UTILS_HPP -#include -#include -#include +#include #include @@ -27,7 +25,7 @@ std::string format_power_exponent(const Power::exponent_type &exponent, /// @param base_str @p base already rendered to a string by the caller /// @return @p base_str, wrapped in parens iff @p base is a Constant whose /// value is a non-integer or negative real -std::string format_power_base(const ExprPtr &base, std::string base_str); +std::string format_power_base(const ExprContainer &base, std::string base_str); } // namespace sequant::detail diff --git a/SeQuant/core/expr.hpp b/SeQuant/core/expr.hpp index 9c8e4ed763..31a09b779a 100644 --- a/SeQuant/core/expr.hpp +++ b/SeQuant/core/expr.hpp @@ -9,6 +9,7 @@ #include #include #include +#include #include #include #include diff --git a/SeQuant/core/expr_fwd.hpp b/SeQuant/core/expr_fwd.hpp index 1dbb36493c..05575986fe 100644 --- a/SeQuant/core/expr_fwd.hpp +++ b/SeQuant/core/expr_fwd.hpp @@ -12,6 +12,7 @@ namespace sequant { class Expr; class ResultExpr; class ExprPtr; +class ExprContainer; class Labeled; class Constant; diff --git a/SeQuant/core/expressions/constant.cpp b/SeQuant/core/expressions/constant.cpp new file mode 100644 index 0000000000..7bf3854572 --- /dev/null +++ b/SeQuant/core/expressions/constant.cpp @@ -0,0 +1,117 @@ +#include +#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; } + +void Constant::adjoint() { + value_ = conj(value_); + reset_hash_value(); +} + +Constant &Constant::operator*=(const Constant &that) { + value_ *= that.value(); + + reset_hash_value(); + + return *this; +} + +Constant &Constant::operator*=(const Expr &that) { + if (!that.is()) { + throw Exception("Constant::operator*=(that): not valid for that"); + } + + return *this *= that.as(); +} + +Constant &Constant::operator+=(const Constant &that) { + value_ += that.value(); + + reset_hash_value(); + + return *this; +} + +Constant &Constant::operator+=(const Expr &that) { + if (!that.is()) { + throw Exception("Constant::operator+=(that): not valid for that"); + } + + return *this += that.as(); +} + +Constant &Constant::operator-=(const Constant &that) { + value_ -= that.value(); + + reset_hash_value(); + + return *this; +} + +Constant &Constant::operator-=(const Expr &that) { + if (!that.is()) { + throw Exception("Constant::operator-=(that): not valid for that"); + } + + return *this -= that.as(); +} + +bool Constant::is_zero(scalar_type v) { return v.is_zero(); } + +bool Constant::is_zero() const { return is_zero(this->value()); } + +std::unique_ptr Constant::unique_copy() const { + return std::make_unique(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(); +} + +Constant operator*(const Constant &lhs, const Constant &rhs) { + Constant result(lhs); + + result *= rhs; + + return result; +} + +Constant operator+(const Constant &lhs, const Constant &rhs) { + Constant result(lhs); + + result += rhs; + + return result; +} + +Constant operator-(const Constant &lhs, const Constant &rhs) { + Constant result(lhs); + + result -= rhs; + + return result; +} + +} // namespace sequant diff --git a/SeQuant/core/expressions/constant.hpp b/SeQuant/core/expressions/constant.hpp index 9ae67a9cb2..9d433f32d2 100644 --- a/SeQuant/core/expressions/constant.hpp +++ b/SeQuant/core/expressions/constant.hpp @@ -3,17 +3,18 @@ #include #include -#include -#include #include #include #include +#include #include 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,70 +68,47 @@ 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"}"; - } - - type_id_type type_id() const override { return get_type_id(); } + std::wstring to_latex() const override; - bool is_scalar() const override { return true; } + type_id_type type_id() const override; - ExprPtr clone() const override { return ex(this->value()); } + bool is_scalar() 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 Constant &that); + 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 Constant &that); + 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 Constant &that); + 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; + + protected: + std::unique_ptr unique_copy() const override; 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 +Constant operator*(const Constant &lhs, const Constant &rhs); +Constant operator+(const Constant &lhs, const Constant &rhs); +Constant operator-(const Constant &lhs, const Constant &rhs); + } // namespace sequant #endif // SEQUANT_EXPRESSIONS_CONSTANT_HPP diff --git a/SeQuant/core/expressions/expr.cpp b/SeQuant/core/expressions/expr.cpp index a8ce23862e..59fdd2c2d3 100644 --- a/SeQuant/core/expressions/expr.cpp +++ b/SeQuant/core/expressions/expr.cpp @@ -2,32 +2,13 @@ // 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 + +#include namespace sequant { @@ -77,506 +58,38 @@ 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()); +std::wstring Expr::to_latex() const { + throw Exception("to_latex not implemented for " + type_name()); } -const Expr &ExprPtr::operator*() const & { - SEQUANT_ASSERT(this->operator bool()); - return *(this->get()); -} +ExprPtr Expr::clone() const { return unique_copy(); } -Expr &&ExprPtr::operator*() && { - SEQUANT_ASSERT(this->operator bool()); - return std::move(*(this->get())); +std::shared_ptr Expr::shared_from_this() { + return std::enable_shared_from_this::shared_from_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::shared_ptr Expr::shared_from_this() const { + return std::enable_shared_from_this::shared_from_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; +std::weak_ptr Expr::weak_from_this() { + return std::enable_shared_from_this::weak_from_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::weak_ptr Expr::weak_from_this() const { + return std::enable_shared_from_this::weak_from_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); -} - -bool proportional_to::operator()(const ExprPtr &expr1, - const ExprPtr &expr2) const { - if (expr1->type_id() != - expr2->type_id()) { // if expr1 is a Product with single factor == expr2, - // or vice versa +bool proportional_to::operator()(const Expr &expr1, const Expr &expr2) const { + if (expr1.type_id() != expr2.type_id()) { + // if expr1 is a Product with single factor == expr2, + // or vice versa if (expr1.is()) { return expr1.as().factors().size() == 1 && - expr1.as().factors().front() == expr2; + *expr1.as().factors().front() == expr2; } else if (expr2.is()) { return expr2.as().factors().size() == 1 && - expr2.as().factors().front() == expr1; + *expr2.as().factors().front() == expr1; } else return false; } @@ -587,10 +100,15 @@ bool proportional_to::operator()(const ExprPtr &expr1, return true; } if (expr1.is()) { - return expr1->hash_value() == expr2->hash_value() && + return expr1.hash_value() == expr2.hash_value() && expr1.as().factors() == expr2.as().factors(); } return expr1 == expr2; } +bool proportional_to::operator()(const ExprPtr &expr1, + const ExprPtr &expr2) const { + return (*this)(*expr1, *expr2); +} + } // namespace sequant diff --git a/SeQuant/core/expressions/expr.hpp b/SeQuant/core/expressions/expr.hpp index 7bee7e8bc3..ebbc21072a 100644 --- a/SeQuant/core/expressions/expr.hpp +++ b/SeQuant/core/expressions/expr.hpp @@ -59,6 +59,8 @@ static const wchar_t adjoint_label = L'\u207A'; /// @endcode class Expr : public std::enable_shared_from_this { public: + friend class ExprContainer; + using hash_type = std::size_t; using type_id_type = int; // to speed up comparisons @@ -75,31 +77,48 @@ 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; + ExprPtr clone() const; + + [[deprecated("Expr objects may no longer be managed by shared_ptr")]] std:: + shared_ptr + shared_from_this(); + [[deprecated("Expr objects may no longer be managed by shared_ptr")]] std:: + shared_ptr + shared_from_this() const; + [[deprecated("Expr objects may no longer be managed by shared_ptr")]] std:: + weak_ptr + weak_from_this(); + [[deprecated("Expr objects may no longer be managed by shared_ptr")]] std:: + weak_ptr + weak_from_this() const; /// like Expr::shared_from_this, but returns ExprPtr /// @return a shared_ptr to this object wrapped into ExprPtr, if this object /// is already managed by a shared_ptr, else returns a shared_ptr to a clone /// of this object wrapped into ExprPtr - ExprPtr exprptr_from_this() { + [[deprecated("Expr objects may no longer be managed by shared_ptr")]] ExprPtr + exprptr_from_this() { + SEQUANT_PRAGMA_IGNORE_DEPRECATED_BEGIN if (weak_from_this().use_count() == 0) return this->clone(); else return static_cast(this->shared_from_this()); + SEQUANT_PRAGMA_IGNORE_DEPRECATED_END } /// like Expr::shared_from_this, but returns ExprPtr /// @return a shared_ptr to this object wrapped into ExprPtr, if this object /// is already managed by a shared_ptr, else returns a shared_ptr to a clone /// of this object wrapped into ExprPtr - ExprPtr exprptr_from_this() const { + [[deprecated("Expr objects may no longer be managed by shared_ptr")]] ExprPtr + exprptr_from_this() const { + SEQUANT_PRAGMA_IGNORE_DEPRECATED_BEGIN if (weak_from_this().use_count() == 0) return this->clone(); else return static_cast( std::const_pointer_cast(this->shared_from_this())); + SEQUANT_PRAGMA_IGNORE_DEPRECATED_END } /// Canonicalizes @c this and returns the byproduct of canonicalization (e.g. @@ -228,9 +247,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 @@ -241,23 +258,22 @@ class Expr : public std::enable_shared_from_this { /// Expr::memoizing_hash /// @return the hash value for this Expr hash_type hash_value( - std::function &)> hasher = {}) - const { - return hasher ? hasher(shared_from_this()) : memoizing_hash(); + std::function hasher = {}) const { + return hasher ? hasher(*this) : memoizing_hash(); + } + + [[deprecated( + "Use a hashing function that takes a const Expr & instead of " + "shared_ptr")]] hash_type + hash_value(std::function &)> + hasher) const { + return hasher ? hasher(this->clone()) : memoizing_hash(); } /// 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 +346,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; @@ -395,9 +379,10 @@ class Expr : public std::enable_shared_from_this { typename E, typename Visitor, typename = std::enable_if_t, Expr>>> static bool visit_impl(E &&expr, Visitor &&visitor, const bool atoms_only) { - if (expr.weak_from_this().use_count() == 0) - throw Exception( - "Expr::visit: cannot visit expressions not managed by shared_ptr"); + constexpr bool visitor_uses_exprptr = + std::is_invocable_r_v, + ExprPtr &>; + for (auto &subexpr_ptr : expr.expr()) { const auto subexpr_is_an_atom = subexpr_ptr->is_atom(); const auto need_to_visit_subexpr = !atoms_only || subexpr_is_an_atom; @@ -406,18 +391,32 @@ class Expr : public std::enable_shared_from_this { visited = visit_impl(*subexpr_ptr, std::forward(visitor), atoms_only); // call on the subexpression itself, if not yet done so - if (need_to_visit_subexpr && !visited) visitor(subexpr_ptr); + if (need_to_visit_subexpr && !visited) { + if constexpr (visitor_uses_exprptr) { + visitor(subexpr_ptr); + } else { + visitor(*subexpr_ptr); + } + } } + // N.B. can only visit itself if visitor is nonmutating! bool this_visited = false; if constexpr (std::is_invocable_r_v, const ExprPtr &>) { if (!atoms_only || expr.is_atom()) { - const ExprPtr this_exprptr = expr.exprptr_from_this(); - visitor(this_exprptr); + visitor(expr.clone()); + this_visited = true; + } + } else if constexpr (std::is_invocable_r_v, + const Expr &>) { + if (!atoms_only || expr.is_atom()) { + visitor(std::as_const(expr)); this_visited = true; } } + return this_visited; } @@ -441,14 +440,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 @@ -469,6 +461,8 @@ class Expr : public std::enable_shared_from_this { return true; } + virtual std::unique_ptr unique_copy() const = 0; + private: /// @return returns next type id in the grand class list static type_id_type get_next_type_id() { @@ -483,12 +477,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); @@ -514,6 +502,7 @@ struct proportional_to { /// @param[in] expr1 /// @param[in] expr2 /// @return true if @p expr1 is proportional to @p expr2 + bool operator()(const Expr &expr1, const Expr &expr2) const; bool operator()(const ExprPtr &expr1, const ExprPtr &expr2) const; }; diff --git a/SeQuant/core/expressions/expr_algorithms.cpp b/SeQuant/core/expressions/expr_algorithms.cpp index 982f2f09dd..df21424cbc 100644 --- a/SeQuant/core/expressions/expr_algorithms.cpp +++ b/SeQuant/core/expressions/expr_algorithms.cpp @@ -19,10 +19,10 @@ namespace sequant { -std::wstring to_latex_align(const ExprPtr& exprptr, size_t max_lines_per_align, +std::wstring to_latex_align(const Expr& expr, size_t max_lines_per_align, size_t max_terms_per_line) { - std::wstring result = io::latex::to_string(exprptr); - if (exprptr->is()) { + std::wstring result = io::latex::to_string(expr); + if (expr.is()) { result.erase(0, 7); // remove leading "{ \bigl" result.replace(result.size() - 8, 8, L")"); // replace trailing "\bigr) }" with ")" @@ -78,6 +78,10 @@ std::wstring to_latex_align(const ExprPtr& exprptr, size_t max_lines_per_align, result += L"\n\\end{align}"; return result; } +std::wstring to_latex_align(const ExprPtr& exprptr, size_t max_lines_per_align, + size_t max_terms_per_line) { + return to_latex_align(*exprptr, max_lines_per_align, max_terms_per_line); +} std::size_t size(const Expr& expr) { return ranges::size(expr); } diff --git a/SeQuant/core/expressions/expr_algorithms.hpp b/SeQuant/core/expressions/expr_algorithms.hpp index 1bd0454e54..c21a9deade 100644 --- a/SeQuant/core/expressions/expr_algorithms.hpp +++ b/SeQuant/core/expressions/expr_algorithms.hpp @@ -22,6 +22,8 @@ namespace sequant { /// @param max_lines_per_align the maximum number of lines in the align before /// starting new align block (if zero, will produce single align block) /// @param max_terms_per_line the maximum number of terms per line +std::wstring to_latex_align(const Expr& expr, size_t max_lines_per_align = 0, + size_t max_terms_per_line = 1); std::wstring to_latex_align(const ExprPtr& exprptr, size_t max_lines_per_align = 0, size_t max_terms_per_line = 1); diff --git a/SeQuant/core/expressions/expr_container.cpp b/SeQuant/core/expressions/expr_container.cpp new file mode 100644 index 0000000000..7aaefc3e8c --- /dev/null +++ b/SeQuant/core/expressions/expr_container.cpp @@ -0,0 +1,239 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace sequant { + +ExprPtr to_expr_ptr(ExprContainer &&container) { + return std::move(container).take_expr(); +} + +ExprContainer::ExprContainer(const ExprContainer &container) + : ExprContainer(container->unique_copy()) {} + +ExprContainer::ExprContainer(std::unique_ptr expr) + : expr_(std::move(expr)) { + SEQUANT_ASSERT(expr_ != nullptr); +} + +ExprContainer::ExprContainer(const Expr &expr) + : ExprContainer(expr.unique_copy()) {} + +ExprContainer::ExprContainer(Expr &&expr) + : ExprContainer(std::move(expr).unique_copy()) {} + +ExprContainer::ExprContainer(const ExprPtr &ptr) + : ExprContainer(ptr->unique_copy()) {} + +ExprContainer &ExprContainer::operator=(const ExprContainer &container) { + // Copy-and-swap + ExprContainer copy(container); + + swap(*this, copy); + + return *this; +} + +ExprContainer &ExprContainer::operator=(const Expr &expr) { + // Copy-and-swap + ExprContainer copy(expr); + + swap(*this, copy); + + return *this; +} + +ExprContainer &ExprContainer::operator=(Expr &&expr) { + // Copy-and-swap + ExprContainer copy(std::move(expr)); + + swap(*this, copy); + + return *this; +} + +ExprContainer ExprContainer::copy() const { return expr_->unique_copy(); } + +std::unique_ptr ExprContainer::take_expr() && { return std::move(expr_); } + +ExprIterator ExprContainer::begin() { + SEQUANT_ASSERT(expr_); + return expr_->begin(); +} + +ExprIterator ExprContainer::end() { + SEQUANT_ASSERT(expr_); + return expr_->end(); +} + +ConstExprIterator ExprContainer::begin() const { + SEQUANT_ASSERT(expr_); + return std::as_const(*expr_).begin(); +} + +ConstExprIterator ExprContainer::end() const { + SEQUANT_ASSERT(expr_); + return std::as_const(*expr_).end(); +} + +ConstExprIterator ExprContainer::cbegin() const { return begin(); } + +ConstExprIterator ExprContainer::cend() const { return end(); } + +ExprContainer::operator const Expr &() const { return *expr_; } + +ExprContainer::operator Expr &() & { return *expr_; } +ExprContainer::operator Expr &&() && { return std::move(*expr_); } + +const Expr &ExprContainer::operator*() const & { return *expr_; } + +Expr &ExprContainer::operator*() & { return *expr_; } + +Expr &&ExprContainer::operator*() && { return std::move(*expr_); } + +const Expr *ExprContainer::operator->() const { return expr_.get(); } + +Expr *ExprContainer::operator->() { return expr_.get(); } + +ExprContainer &ExprContainer::operator+=(const Expr &expr) { + if (expr_->is()) { + expr_->as() += expr; + } else if (expr_->is() && expr.is()) { + *this = expr_->as() + expr.as(); + } else { + *this = Sum(ExprPtrList{to_expr_ptr(std::move(*this)), expr.clone()}); + } + + return *this; +} + +ExprContainer &ExprContainer::operator-=(const Expr &expr) { + if (expr_->is()) { + expr_->as() -= expr; + } else if (expr_->is() && expr.is()) { + *this = expr_->as() - expr.as(); + } else { + *this = Sum(ExprPtrList{to_expr_ptr(std::move(*this)), + ex(-1, ExprPtrList{expr.clone()})}); + } + + return *this; +} + +ExprContainer &ExprContainer::operator*=(const Expr &expr) { + if (expr_->is()) { + expr_->as() *= expr; + } else if (expr_->is() && expr.is()) { + *this = expr_->as() * expr.as(); + } else { + *this = Product(ExprPtrList{to_expr_ptr(std::move(*this)), expr.clone()}); + } + return *this; +} + +ExprContainer &ExprContainer::operator^=(const Expr &other) { + auto this_is_product = expr_->is(); + auto other_is_product = other.is(); + if (!this_is_product && !other_is_product) { + *this = + NCProduct(ExprPtrList{to_expr_ptr(std::move(*this)), other.clone()}); + } else if (this_is_product) { + *this = NCProduct(std::move(expr_->as())); + expr_->as().append(1, other.clone()); + } else { // other_is_product + NCProduct result(other.clone().as()); + result.prepend(1, to_expr_ptr(std::move(*this))); + *this = std::move(result); + } + + return *this; +} + +void swap(ExprContainer &lhs, ExprContainer &rhs) { + std::swap(lhs.expr_, rhs.expr_); +} + +ExprContainer operator+(const Expr &lhs, const Expr &rhs) { + ExprContainer cont(lhs); + cont += rhs; + + return cont; +} + +ExprContainer operator+(const ExprContainer &lhs, const Expr &rhs) { + return static_cast(lhs) + rhs; +} + +ExprContainer operator+(const Expr &lhs, const ExprContainer &rhs) { + return lhs + static_cast(rhs); +} + +ExprContainer operator-(const Expr &lhs, const Expr &rhs) { + ExprContainer cont(lhs); + cont -= rhs; + + return cont; +} + +ExprContainer operator-(const ExprContainer &lhs, const Expr &rhs) { + return static_cast(lhs) - rhs; +} + +ExprContainer operator-(const Expr &lhs, const ExprContainer &rhs) { + return lhs - static_cast(rhs); +} + +ExprContainer operator*(const Expr &lhs, const Expr &rhs) { + ExprContainer cont(lhs); + cont *= rhs; + + return cont; +} + +ExprContainer operator*(const ExprContainer &lhs, const Expr &rhs) { + return static_cast(lhs) * rhs; +} + +ExprContainer operator*(const Expr &lhs, const ExprContainer &rhs) { + return lhs * static_cast(rhs); +} + +ExprContainer operator^(const Expr &lhs, const Expr &rhs) { + ExprContainer cont(lhs); + cont ^= rhs; + + return cont; +} + +ExprContainer operator^(const ExprContainer &lhs, const Expr &rhs) { + return static_cast(lhs) ^ rhs; +} + +ExprContainer operator^(const Expr &lhs, const ExprContainer &rhs) { + return lhs ^ static_cast(rhs); +} + +bool operator==(const ExprContainer &lhs, const ExprPtr &rhs) { + return *lhs == *rhs; +} + +bool operator==(const ExprPtr &lhs, const ExprContainer &rhs) { + return *lhs == *rhs; +} + +ExprContainer adjoint(const ExprContainer &cont) { + ExprContainer copy = cont.copy(); + copy->adjoint(); + return copy; +} + +} // namespace sequant diff --git a/SeQuant/core/expressions/expr_container.hpp b/SeQuant/core/expressions/expr_container.hpp new file mode 100644 index 0000000000..23d97c8424 --- /dev/null +++ b/SeQuant/core/expressions/expr_container.hpp @@ -0,0 +1,92 @@ +#ifndef SEQUANT_EXPRESSIONS_EXPR_CONTAINER_HPP +#define SEQUANT_EXPRESSIONS_EXPR_CONTAINER_HPP + +#include + +#include +#include + +namespace sequant { + +class Expr; +class ExprPtr; + +class ExprContainer { + public: + explicit ExprContainer(const ExprContainer &container); + ExprContainer(ExprContainer &&container) = default; + + explicit ExprContainer(const Expr &expr); + ExprContainer(Expr &&expr); + + explicit ExprContainer(const ExprPtr &ptr); + + ExprContainer &operator=(const ExprContainer &container); + ExprContainer &operator=(ExprContainer &&container) = default; + + ExprContainer &operator=(const Expr &expr); + ExprContainer &operator=(Expr &&expr); + + ~ExprContainer() = default; + + ExprContainer copy() const; + + std::unique_ptr take_expr() &&; + + ExprIterator begin(); + ExprIterator end(); + ConstExprIterator begin() const; + ConstExprIterator end() const; + ConstExprIterator cbegin() const; + ConstExprIterator cend() const; + + operator const Expr &() const; + operator Expr &() &; + operator Expr &&() &&; + + const Expr &operator*() const &; + Expr &operator*() &; + Expr &&operator*() &&; + + const Expr *operator->() const; + Expr *operator->(); + + ExprContainer &operator+=(const Expr &expr); + ExprContainer &operator-=(const Expr &expr); + ExprContainer &operator*=(const Expr &expr); + ExprContainer &operator^=(const Expr &expr); + + friend void swap(ExprContainer &, ExprContainer &); + + private: + std::unique_ptr expr_; + + ExprContainer(std::unique_ptr expr); +}; + +using ExprContainerList = std::initializer_list; + +ExprContainer operator+(const Expr &lhs, const Expr &rhs); +ExprContainer operator+(const ExprContainer &lhs, const Expr &rhs); +ExprContainer operator+(const Expr &lhs, const ExprContainer &rhs); + +ExprContainer operator-(const Expr &lhs, const Expr &rhs); +ExprContainer operator-(const ExprContainer &lhs, const Expr &rhs); +ExprContainer operator-(const Expr &lhs, const ExprContainer &rhs); + +ExprContainer operator*(const Expr &lhs, const Expr &rhs); +ExprContainer operator*(const ExprContainer &lhs, const ExprContainer &rhs); +ExprContainer operator*(const Expr &lhs, const Expr &rhs); + +ExprContainer operator^(const Expr &lhs, const Expr &rhs); +ExprContainer operator^(const ExprContainer &lhs, const Expr &rhs); +ExprContainer operator^(const Expr &lhs, const ExprContainer &rhs); + +bool operator==(const ExprContainer &lhs, const ExprPtr &rhs); +bool operator==(const ExprPtr &lhs, const ExprContainer &rhs); + +ExprContainer adjoint(const ExprContainer &expr); + +} // namespace sequant + +#endif // SEQUANT_EXPRESSIONS_EXPR_CONTAINER_HPP diff --git a/SeQuant/core/expressions/expr_operators.hpp b/SeQuant/core/expressions/expr_operators.hpp index e1e0967e3c..4db7740876 100644 --- a/SeQuant/core/expressions/expr_operators.hpp +++ b/SeQuant/core/expressions/expr_operators.hpp @@ -14,175 +14,89 @@ #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) { +ExprPtr operator+(const std::same_as auto &lhs, T &&rhs) { return lhs + ex(std::forward(rhs)); } template requires(std::constructible_from) -ExprPtr operator+(T &&lhs, const ExprPtr &rhs) { +ExprPtr operator+(T &&lhs, const std::same_as auto &rhs) { return ex(std::forward(lhs)) + rhs; } template requires(std::constructible_from) -ExprPtr operator-(const ExprPtr &lhs, T &&rhs) { +ExprPtr operator-(const std::same_as auto &lhs, T &&rhs) { return lhs - ex(std::forward(rhs)); } template requires(std::constructible_from) -ExprPtr operator-(T &&lhs, const ExprPtr &rhs) { +ExprPtr operator-(T &&lhs, const std::same_as auto &rhs) { return ex(std::forward(lhs)) - rhs; } template requires(std::constructible_from) -ExprPtr operator*(const ExprPtr &lhs, T &&rhs) { +ExprPtr operator*(const std::same_as auto &lhs, T &&rhs) { return lhs * ex(std::forward(rhs)); } template requires(std::constructible_from) -ExprPtr operator*(T &&lhs, const ExprPtr &rhs) { +ExprPtr operator*(T &&lhs, const std::same_as auto &rhs) { return ex(std::forward(lhs)) * rhs; } template requires(std::is_arithmetic_v) -ExprPtr operator/(const ExprPtr &lhs, T &&rhs) { +ExprPtr operator/(const std::same_as auto &lhs, T &&rhs) { return lhs * ex(rational(1, std::forward(rhs))); } -inline ExprPtr operator/(const ExprPtr &lhs, const Constant &rhs) { +inline ExprPtr operator/(const std::same_as auto &lhs, + const Constant &rhs) { return lhs * ex(1.0 / rhs.value()); } template requires(std::constructible_from) -ExprPtr operator+(T &&lhs, const ExprPtr &rhs) { +ExprPtr operator+(T &&lhs, const std::same_as auto &rhs) { return ex(std::forward(lhs)) + rhs; } template requires(std::constructible_from) -ExprPtr operator+(const ExprPtr &lhs, T &&rhs) { +ExprPtr operator+(const std::same_as auto &lhs, T &&rhs) { return lhs + ex(std::forward(rhs)); } template requires(std::constructible_from) -ExprPtr operator-(T &&lhs, const ExprPtr &rhs) { +ExprPtr operator-(T &&lhs, const std::same_as auto &rhs) { return ex(std::forward(lhs)) - rhs; } template requires(std::constructible_from) -ExprPtr operator-(const ExprPtr &lhs, T &&rhs) { +ExprPtr operator-(const std::same_as auto &lhs, T &&rhs) { return lhs - ex(std::forward(rhs)); } template requires(std::constructible_from) -ExprPtr operator*(T &&lhs, const ExprPtr &rhs) { +ExprPtr operator*(T &&lhs, const std::same_as auto &rhs) { return ex(std::forward(lhs)) * rhs; } template requires(std::constructible_from) -ExprPtr operator*(const ExprPtr &lhs, T &&rhs) { +ExprPtr operator*(const std::same_as auto &lhs, T &&rhs) { return lhs * ex(std::forward(rhs)); } diff --git a/SeQuant/core/expressions/expr_ptr.cpp b/SeQuant/core/expressions/expr_ptr.cpp new file mode 100644 index 0000000000..9276b92db7 --- /dev/null +++ b/SeQuant/core/expressions/expr_ptr.cpp @@ -0,0 +1,196 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace sequant { + +ExprPtr::ExprPtr(const ExprContainer &container) : ExprPtr(container.copy()) {} + +ExprPtr::ExprPtr(ExprContainer &&container) + : ExprPtr(std::move(container).take_expr()) {} + +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..b453c16f6a 100644 --- a/SeQuant/core/expressions/expr_ptr.hpp +++ b/SeQuant/core/expressions/expr_ptr.hpp @@ -10,6 +10,8 @@ namespace sequant { +class ExprContainer; + /// @brief ExprPtr is a multiple-owner smart pointer to Expr /// It can be used mostly interchangeably with `std::shared_ptr`, but @@ -23,6 +25,8 @@ class ExprPtr : public std::shared_ptr { ExprPtr() = default; ExprPtr(const ExprPtr &) = default; ExprPtr(ExprPtr &&) = default; + explicit ExprPtr(const ExprContainer &container); + ExprPtr(ExprContainer &&container); template , Expr> || std::is_base_of_v>>> @@ -157,6 +161,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..ea8eeb259d --- /dev/null +++ b/SeQuant/core/expressions/power.cpp @@ -0,0 +1,212 @@ +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace sequant { + +Power::Power(const ExprPtr& base, exponent_type exponent) + : Power(ExprContainer(base), std::move(exponent)) {} + +Power::Power(ExprContainer base, exponent_type exponent) + : base_{std::move(base)}, exponent_{std::move(exponent)} { + SEQUANT_ASSERT(base_->is() || base_->is()); + // 0^n is defined only for n >= 0 (0^0 = 1 by convention) + SEQUANT_ASSERT(!base_->is() || !base_->as().is_zero() || + exponent_ >= 0); +} + +const ExprContainer& 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(); +} + +template +void flatt_impl(E& expr) { + auto create_constant = [](const auto& val) { + if constexpr (std::same_as, ExprContainer>) { + return Constant(val); + } else { + return ex(val); + } + }; + + const auto& pw = expr->template as(); + + // b^1 = b and conjugate if needed + if (pw.exponent() == 1) { + auto lifted = pw.base().copy(); + 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 = create_constant(Constant::scalar_type{1}); + return; + } + if (!pw.base()->template is()) return; + + using scalar_type = Constant::scalar_type; + const auto& base_val = pw.base()->template as().value(); + + // 1^k = 1 for any rational k. + if (base_val == scalar_type{1}) { + expr = create_constant(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()); + + 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 = create_constant(std::move(value)); +} + +void Power::flatten(ExprPtr& expr) { + if (!expr || !expr->is()) return; + + flatt_impl(expr); +} + +void Power::flatten(ExprContainer& expr) { + if (!expr->is()) return; + + flatt_impl(expr); +} + +Expr::type_id_type Power::type_id() const { return get_type_id(); } + +bool Power::is_scalar() const { return true; } + +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"); +} + +std::unique_ptr Power::unique_copy() const { + auto copy = std::make_unique(base_.copy(), exponent_); + if (conjugated_) copy->as().conjugate(); + return copy; +} + +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..86b9840156 100644 --- a/SeQuant/core/expressions/power.hpp +++ b/SeQuant/core/expressions/power.hpp @@ -3,12 +3,14 @@ #include #include +#include #include +#include #include -#include -#include #include -#include + +#include +#include namespace sequant { @@ -27,57 +29,40 @@ 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(const ExprPtr& base, exponent_type exponent); + Power(ExprContainer base, exponent_type exponent); /// @overload constructs a `Variable` base from @p label template - requires std::constructible_from && - (!std::convertible_to) + requires(std::constructible_from && !expr_holder) Power(L&& label, exponent_type exponent) : Power(ex(std::forward(label)), std::move(exponent)) {} /// @overload constructs a `Constant` base from scalar @p value template - requires(!std::constructible_from && - !std::convertible_to && + requires(!std::constructible_from && !expr_holder && std::constructible_from) Power(V&& value, exponent_type exponent) : Power(ex(std::forward(value)), std::move(exponent)) {} /// @return the base expression - const ExprPtr& base() const { return base_; } + const ExprContainer& 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 +78,15 @@ 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(); - - // 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 + static void flatten(ExprPtr& expr); + static void flatten(ExprContainer& expr); - // 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; - } + type_id_type type_id() const override; - // 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; - } + bool is_scalar() 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,74 +96,23 @@ 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); + + protected: + std::unique_ptr unique_copy() const override; private: - ExprPtr base_; + ExprContainer base_; exponent_type exponent_; bool conjugated_ = false; /// @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..e904eba3c0 --- /dev/null +++ b/SeQuant/core/expressions/product.cpp @@ -0,0 +1,399 @@ +#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::unique_ptr Product::unique_copy() const { + return std::make_unique(deep_copy()); +} + +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(); +}; + +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, that.clone()); + } 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..d40c957a61 100644 --- a/SeQuant/core/expressions/product.hpp +++ b/SeQuant/core/expressions/product.hpp @@ -4,18 +4,14 @@ #include #include #include -#include #include #include -#include #include #include -#include -#include #include -#include +#include #include #include @@ -35,6 +31,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 +43,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 @@ -69,7 +61,7 @@ class Product : public Expr { if constexpr (rng_is_expr || rng_is_exprptr) { ExprPtr rng_as_exprptr; if constexpr (rng_is_expr) { - rng_as_exprptr = rng.exprptr_from_this(); + rng_as_exprptr = rng.clone(); } else { rng_as_exprptr = rng; } @@ -188,19 +180,17 @@ class Product : public Expr { typename = std::enable_if_t>> Product &append(T scalar, Factor &&factor, Flatten flatten_tag = Flatten::Yes) { - return this->append(scalar, - std::static_pointer_cast( - std::forward(factor).shared_from_this()), - flatten_tag); + return this->append( + scalar, + std::static_pointer_cast(std::forward(factor).clone()), + flatten_tag); } /// (post-)multiplies the product by@c factor /// @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 @@ -209,9 +199,9 @@ class Product : public Expr { /// @warning if @p factor is a Product, it is flattened recursively template >> Product &append(Factor &&factor, Flatten flatten_tag = Flatten::Yes) { - return this->append(std::static_pointer_cast( - std::forward(factor).shared_from_this()), - flatten_tag); + return this->append( + std::static_pointer_cast(std::forward(factor).clone()), + flatten_tag); } /// (pre-)multiplies the product by @c scalar times @c factor @@ -262,19 +252,19 @@ class Product : public Expr { typename = std::enable_if_t>> Product &prepend(T scalar, Factor &&factor, Flatten flatten_tag = Flatten::Yes) { - return this->prepend(scalar, - std::static_pointer_cast( - std::forward(factor).shared_from_this()), - flatten_tag); + return this->prepend( + scalar, + std::static_pointer_cast(std::forward(factor).clone()), + 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 +284,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 +299,71 @@ 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(); }; - - /// @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; - } + type_id_type type_id() 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; + + protected: + std::unique_ptr unique_copy() const 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 +371,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..51fbaa7c81 --- /dev/null +++ b/SeQuant/core/expressions/sum.cpp @@ -0,0 +1,388 @@ +#include +#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(); }; + +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(*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(*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(*this) + << std::endl; + } + + return {}; // side effects are absorbed into summands +} + +Sum &Sum::operator+=(const Expr &that) { + this->append(that.clone()); + return *this; +} + +Sum &Sum::operator-=(const Expr &that) { + if (that.is()) + this->append(ex(-that.as().value())); + else + this->append(ex(-1, ExprPtrList{that.clone()})); + 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()}; +} + +std::unique_ptr Sum::unique_copy() const { + auto cloned_summands = + summands() | + ranges::views::transform([](const ExprPtr &ptr) { return ptr->clone(); }); + return std::make_unique(ranges::begin(cloned_summands), + ranges::end(cloned_summands)); +} + +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..b70fe70686 100644 --- a/SeQuant/core/expressions/sum.hpp +++ b/SeQuant/core/expressions/sum.hpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include @@ -16,7 +17,10 @@ #include #include +#include +#include #include +#include #include namespace sequant { @@ -35,20 +39,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 @@ -63,19 +56,18 @@ class Sum : public Expr { /// construct a Sum out of a range of summands /// @param rng a range - template - requires(meta::is_range_v> && - !meta::is_same_v, ExprPtrList>) + template + requires(!std::same_as, ExprPtrList>) explicit Sum(Range &&rng) { // N.B. use append to flatten out Sum summands - constexpr auto rng_is_expr = - meta::is_base_of_v>; + constexpr auto rng_is_expr = is_an_expr_v>; constexpr auto rng_is_exprptr = - meta::is_same_v>; + std::same_as>; + if constexpr (rng_is_expr || rng_is_exprptr) { ExprPtr rng_as_exprptr; if constexpr (rng_is_expr) { - rng_as_exprptr = rng.exprptr_from_this(); + rng_as_exprptr = rng.clone(); } else { rng_as_exprptr = rng; } @@ -93,207 +85,67 @@ 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 - /// @param f an object of Filter type - /// Selects elements {`e`} for which `f(e)` is true - template + /// @param f Boolean predicate + /// @returns A sum containing only the summands for which f was true. + template Filter> ExprPtr filter(Filter &&f) const { - return ex(summands_ | ranges::views::filter(f)); + return ex(summands_ | + ranges::views::transform([](const auto &e) { return *e; }) | + ranges::views::filter(f)); } /// @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; /// @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; - ConstExprIterator begin_subexpr() const override { - return ConstExprIterator{summands_.data()}; - } + ConstExprIterator end_subexpr() const override; - ConstExprIterator end_subexpr() const override { - return ConstExprIterator{summands_.data() + summands_.size()}; - } + protected: + std::unique_ptr unique_copy() const override; private: summands_type summands_{}; @@ -304,61 +156,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 +192,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/tensor.cpp b/SeQuant/core/expressions/tensor.cpp index 0495e2dd22..ea546c8c43 100644 --- a/SeQuant/core/expressions/tensor.cpp +++ b/SeQuant/core/expressions/tensor.cpp @@ -45,4 +45,8 @@ ExprPtr Tensor::canonicalize(CanonicalizeOptions) { return canonicalizer_ptr ? canonicalizer_ptr->apply(*this) : ExprPtr{}; } +std::unique_ptr Tensor::unique_copy() const { + return std::make_unique(*this); +} + } // namespace sequant diff --git a/SeQuant/core/expressions/tensor.hpp b/SeQuant/core/expressions/tensor.hpp index 07b3b7d5c2..9023bc0af8 100644 --- a/SeQuant/core/expressions/tensor.hpp +++ b/SeQuant/core/expressions/tensor.hpp @@ -754,8 +754,6 @@ class Tensor : public Expr, public AbstractTensor, public MutatableLabeled { type_id_type type_id() const override { return get_type_id(); }; - ExprPtr clone() const override { return ex(*this); } - void reset_tags() const { ranges::for_each(slots(), [](const auto &idx) { idx.reset_tag(); }); } @@ -775,6 +773,9 @@ class Tensor : public Expr, public AbstractTensor, public MutatableLabeled { return false; // TODO do we compare typeid? labels? probably the latter } + protected: + std::unique_ptr unique_copy() const override; + private: std::wstring label_{}; sequant::bra bra_{}; diff --git a/SeQuant/core/expressions/traits.hpp b/SeQuant/core/expressions/traits.hpp index efeb818d82..76c793bf81 100644 --- a/SeQuant/core/expressions/traits.hpp +++ b/SeQuant/core/expressions/traits.hpp @@ -4,6 +4,9 @@ #include #include +#include +#include + namespace sequant { template @@ -36,6 +39,10 @@ constexpr bool is_a_power_v = meta::is_base_of_v; template constexpr bool is_power_v = meta::is_same_v; +template +concept expr_holder = std::same_as, ExprPtr> || + std::same_as, ExprContainer>; + } // namespace sequant #endif // SEQUANT_EXPRESSIONS_TRAITS_HPP diff --git a/SeQuant/core/expressions/variable.cpp b/SeQuant/core/expressions/variable.cpp new file mode 100644 index 0000000000..15d614f478 --- /dev/null +++ b/SeQuant/core/expressions/variable.cpp @@ -0,0 +1,70 @@ +#include +#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; +} + +std::unique_ptr Variable::unique_copy() const { + return std::make_unique(*this); +} + +void Variable::adjoint() { conjugate(); } + +} // namespace sequant diff --git a/SeQuant/core/expressions/variable.hpp b/SeQuant/core/expressions/variable.hpp index 650ce64a39..de13956ae3 100644 --- a/SeQuant/core/expressions/variable.hpp +++ b/SeQuant/core/expressions/variable.hpp @@ -2,17 +2,17 @@ #define SEQUANT_EXPRESSIONS_VARIABLE_HPP #include -#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 +29,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,39 +47,24 @@ class Variable : public Expr, public MutatableLabeled { std::wstring to_latex() const override; - type_id_type type_id() const override { return get_type_id(); } - - bool is_scalar() const override { return true; } + type_id_type type_id() const override; - ExprPtr clone() const override; + bool is_scalar() const override; /// @brief adjoint of a Variable is its complex conjugate virtual void adjoint() override; + protected: + std::unique_ptr unique_copy() const override; + private: 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..b304f39853 100644 --- a/SeQuant/core/op.hpp +++ b/SeQuant/core/op.hpp @@ -389,7 +389,10 @@ class Operator : public container::svector>, public Expr { type_id_type type_id() const override { return get_type_id(); }; - ExprPtr clone() const override { return std::make_shared(*this); } + protected: + std::unique_ptr unique_copy() const override { + return std::make_unique(*this); + } private: base_type make_ops(Action action, IndexList indices) { @@ -711,10 +714,6 @@ class NormalOperator : public Operator, return Expr::get_type_id(); }; - ExprPtr clone() const override { - return std::make_shared(*this); - } - virtual void adjoint() override { // same as base adjoint(), but updates extra state Operator::adjoint(); @@ -739,6 +738,11 @@ class NormalOperator : public Operator, return mutated; } + protected: + std::unique_ptr unique_copy() const override { + return std::make_unique(*this); + } + private: Vacuum vacuum_; std::size_t ncreators_ = 0; @@ -955,6 +959,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 " @@ -1068,6 +1079,11 @@ class NormalOperatorSequence : public container::svector>, static_cast(nopseq2); } + protected: + std::unique_ptr unique_copy() const override { + return std::make_unique(*this); + } + private: Vacuum vacuum_ = Vacuum::Physical; /// ensures that all operators use same vacuum, and sets 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/SeQuant/domain/mbpt/op.cpp b/SeQuant/domain/mbpt/op.cpp index 50f7b295d4..171981c375 100644 --- a/SeQuant/domain/mbpt/op.cpp +++ b/SeQuant/domain/mbpt/op.cpp @@ -580,14 +580,14 @@ ExprPtr OpMaker::operator()( if (!dep && csv) { if (opclass == OpClass::Ex) { if constexpr (assert_enabled()) { - for (auto&& s : cre_spaces_) { + for ([[maybe_unused]] const auto& s : cre_spaces_) { SEQUANT_ASSERT(isr->contains_unoccupied(s)); } } dep = UseDepIdx::Bra; } else if (opclass == OpClass::Deex) { if constexpr (assert_enabled()) { - for (auto&& s : ann_spaces_) { + for ([[maybe_unused]] const auto& s : ann_spaces_) { SEQUANT_ASSERT(isr->contains_unoccupied(s)); } } diff --git a/SeQuant/domain/mbpt/op.hpp b/SeQuant/domain/mbpt/op.hpp index 658e1bacba..45e909f841 100644 --- a/SeQuant/domain/mbpt/op.hpp +++ b/SeQuant/domain/mbpt/op.hpp @@ -43,7 +43,7 @@ #include #include #include -#include +#include #include #include #include @@ -950,6 +950,9 @@ class Operator : public Operator { /// @brief returns the perturbation order of this operator [[nodiscard]] size_t order() const { return order_; } + protected: + std::unique_ptr unique_copy() const override; + private: std::function qn_action_; @@ -963,8 +966,6 @@ class Operator : public Operator { Expr::type_id_type type_id() const override; - ExprPtr clone() const override; - std::wstring to_latex() const override; Expr::hash_type memoizing_hash() const override; diff --git a/SeQuant/domain/mbpt/op.ipp b/SeQuant/domain/mbpt/op.ipp index 8f28062b51..2030799552 100644 --- a/SeQuant/domain/mbpt/op.ipp +++ b/SeQuant/domain/mbpt/op.ipp @@ -12,6 +12,8 @@ #include #include +#include + namespace sequant { namespace mbpt { @@ -192,8 +194,8 @@ Expr::type_id_type Operator::type_id() const { }; template -ExprPtr Operator::clone() const { - return ex(*this); +std::unique_ptr Operator::unique_copy() const { + return std::make_unique(*this); } // Expresses general operators in human interpretable form. for example: diff --git a/tests/unit/test_expr.cpp b/tests/unit/test_expr.cpp index b9112c7ee6..9fbf6424e5 100644 --- a/tests/unit/test_expr.cpp +++ b/tests/unit/test_expr.cpp @@ -37,7 +37,10 @@ struct Dummy : public sequant::Expr { virtual ~Dummy() = default; 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(); } + std::unique_ptr unique_copy() const override { + return std::make_unique(); + } + void adjoint() override {} bool static_equal(const sequant::Expr &) const override { return true; } }; @@ -69,6 +72,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()}; @@ -107,8 +112,8 @@ struct VecExpr : public std::vector, public sequant::Expr { static_cast(static_cast(that)); } - sequant::ExprPtr clone() const override { - return sequant::ex(this->begin(), this->end()); + std::unique_ptr unique_copy() const override { + return std::make_unique(this->begin(), this->end()); } }; @@ -120,8 +125,8 @@ struct Adjointable : public sequant::Expr { return L"{\\text{Adjointable}{" + std::to_wstring(v) + L"}}"; } type_id_type type_id() const override { return get_type_id(); }; - sequant::ExprPtr clone() const override { - return sequant::ex(v); + std::unique_ptr unique_copy() const override { + return std::make_unique(v); } bool static_equal(const sequant::Expr &that) const override { return v == that.as().v; @@ -282,16 +287,16 @@ TEST_CASE("expr", "[elements]") { const auto c2 = ex(rational{1, 2}); const auto vx = ex(L"x"); - { // constructors + SECTION("constructors") { REQUIRE_NOTHROW(Power(c2, rational{1, 2})); REQUIRE_NOTHROW(Power(vx, rational{3, 1})); // convenience ctors: REQUIRE(Power(L"x", 2) == Power(vx, rational{2})); REQUIRE(Power(L"x", rational{1, 2}) == Power(vx, rational{1, 2})); - REQUIRE(Power(2, 3) == Power(ex(2), rational{3})); + REQUIRE(Power(2, 3) == Power(Constant(2), rational{3})); REQUIRE(Power(rational{2, 3}, 2) == - Power(ex(rational{2, 3}), rational{2})); + Power(Constant(rational{2, 3}), rational{2})); if constexpr (sequant::assert_behavior() == sequant::AssertBehavior::Throw) { // base must be a Constant or Variable; Power-of-Power is not allowed @@ -299,25 +304,25 @@ TEST_CASE("expr", "[elements]") { REQUIRE_THROWS(Power(inner, rational{2, 3})); // 0^n is defined only for n >= 0 - REQUIRE_THROWS(Power(ex(0), rational{-1})); + REQUIRE_THROWS(Power(Constant(0), rational{-1})); } } - { // accessors + SECTION("accessors") { Power p(c2, rational{1, 2}); - REQUIRE(p.base() == ex(rational{1, 2})); + REQUIRE(p.base() == Constant(rational{1, 2})); REQUIRE(p.exponent() == rational{1, 2}); // is_zero: base == 0, exponent > 0 - Power pz(ex(0), rational{2}); + Power pz(Constant(0), rational{2}); REQUIRE(pz.is_zero()); // 0^0 is not zero by our convention - Power pz2(ex(0), rational{0}); + Power pz2(Constant(0), rational{0}); REQUIRE(!pz2.is_zero()); REQUIRE(!p.is_zero()); } - { // comparison + SECTION("comparison") { Power p1(c2, rational{1, 2}); Power p2(c2, rational{1, 2}); REQUIRE(p1 == p2); @@ -336,7 +341,7 @@ TEST_CASE("expr", "[elements]") { REQUIRE(!(plt_a < plt_c)); } - { // operator*= + SECTION("operator*=") { // b^e1 *= b^e2 -> b^(e1+e2) Power pa(vx, rational{1, 2}); Power pb(vx, rational{1, 3}); @@ -366,14 +371,14 @@ TEST_CASE("expr", "[elements]") { REQUIRE(pc_conj2.exponent() == rational{3, 2}); // 2^{1/2} * 2^{1/2} = 2 - Power pe(ex(2), rational{1, 2}); - Power pf(ex(2), rational{1, 2}); + Power pe(Constant(2), rational{1, 2}); + Power pf(Constant(2), rational{1, 2}); pe *= pf; REQUIRE(pe.exponent() == rational{1}); REQUIRE(to_latex(pe) == Constant(2).to_latex()); } - { // Power should NOT be absorbed into Product::scalar_ + SECTION("Don't absorb into Product::scalar") { auto p = ex(vx, rational{1, 2}); auto prod = ex(Product{}); prod->as().append(1, p, Product::Flatten::Yes); @@ -413,21 +418,17 @@ 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()); REQUIRE_NOTHROW(adjoint(e)); // check free-function adjoint } - { // Constant + SECTION("Constant") { const auto e = std::make_shared(Constant::scalar_type{1, 2}); REQUIRE_NOTHROW(e->adjoint()); REQUIRE(e->value() == Constant::scalar_type{1, -2}); } - { // Variable + SECTION("Variabkle") { const auto e = std::make_shared(L"q"); REQUIRE(e->conjugated() == false); REQUIRE_NOTHROW(e->adjoint()); @@ -436,7 +437,7 @@ TEST_CASE("expr", "[elements]") { REQUIRE_NOTHROW(e->adjoint()); REQUIRE(e->conjugated() == false); } - { // Product + SECTION("Product") { // Product const auto e = std::make_shared(); e->append(Constant::scalar_type{2, -1}, ex()); e->append(1, ex(-2)); @@ -445,7 +446,7 @@ TEST_CASE("expr", "[elements]") { REQUIRE(e->factors()[0]->as().v == 2); REQUIRE(e->factors()[1]->as().v == -1); } - { // CProduct + SECTION("CProduct") { const auto e = std::make_shared(); e->append(Constant::scalar_type{2, -1}, ex()); e->append(1, ex(-2)); @@ -454,7 +455,7 @@ TEST_CASE("expr", "[elements]") { REQUIRE(e->factors()[0]->as().v == -1); REQUIRE(e->factors()[1]->as().v == 2); } - { // NCProduct + SECTION("NCProduct") { const auto e = std::make_shared(); e->append(Constant::scalar_type{2, -1}, ex()); e->append(1, ex(-2)); @@ -463,7 +464,7 @@ TEST_CASE("expr", "[elements]") { REQUIRE(e->factors()[0]->as().v == 2); REQUIRE(e->factors()[1]->as().v == -1); } - { // Sum + SECTION("Sum") { const auto e = std::make_shared(); e->append(ex()); e->append(ex(-2)); @@ -471,7 +472,8 @@ TEST_CASE("expr", "[elements]") { REQUIRE(e->summands()[0]->as().v == -1); REQUIRE(e->summands()[1]->as().v == 2); } - { // Power: adjoint flips the conjugation flag; base/exponent unchanged + SECTION("Power") { + // adjoint flips the conjugation flag; base/exponent unchanged Power pv(ex(L"z"), rational{1, 2}); REQUIRE(!pv.conjugated()); pv.adjoint(); @@ -910,10 +912,14 @@ TEST_CASE("expr", "[elements]") { REQUIRE(hash_value(ex(1)) == hash_value(ex(1))); - auto hasher = [](const std::shared_ptr &) -> unsigned int { + auto hasher1 = [](const std::shared_ptr &) -> unsigned int { return 0; }; - REQUIRE_NOTHROW(ex(1)->hash_value(hasher) == 0); + auto hasher2 = [](const Expr &) -> unsigned int { return 2; }; + SEQUANT_PRAGMA_IGNORE_DEPRECATED_BEGIN + REQUIRE_NOTHROW(ex(1)->hash_value(hasher1) == 0); + SEQUANT_PRAGMA_IGNORE_DEPRECATED_END + REQUIRE_NOTHROW(ex(1)->hash_value(hasher2) == 2); { // Power const auto c2 = ex(rational{1, 2}); @@ -1118,6 +1124,157 @@ TEST_CASE("expr", "[elements]") { } } + SECTION("ExprContainer") { + SECTION("Constructors") { + SECTION("from conrete") { + ExprContainer cont1(Constant(1)); + REQUIRE(cont1->is()); + REQUIRE(cont1->as() == Constant(1)); + + ExprContainer cont2(Variable("bla")); + REQUIRE(cont2->is()); + REQUIRE(cont2->as() == Variable("bla")); + + Product prod(ExprPtrList{ex("bla"), ex(2)}); + ExprContainer cont3(prod); + REQUIRE(cont3->is()); + REQUIRE(cont3->as() == prod); + + Sum sum(ExprPtrList{ex("bla"), ex(2)}); + ExprContainer cont4(sum); + REQUIRE(cont4->is()); + REQUIRE(cont4->as() == sum); + } + SECTION("from base") { + ExprPtr expr = ex(42); + + ExprContainer cont(*expr); + REQUIRE(cont->is()); + REQUIRE(cont->as().value() == 42); + } + SECTION("conversion via 'assignment'") { + ExprContainer cont1 = Constant(1); + REQUIRE(cont1->is()); + REQUIRE(cont1->as() == Constant(1)); + + ExprContainer cont2 = Variable("bla"); + REQUIRE(cont2->is()); + REQUIRE(cont2->as() == Variable("bla")); + + // Note: copy-ctor is explicit so in order for this "assignment" to + // work, we need to assign rvalues + Product prod(ExprPtrList{ex("bla"), ex(2)}); + ExprContainer cont3 = Product(prod); + REQUIRE(cont3->is()); + REQUIRE(cont3->as() == prod); + + Sum sum(ExprPtrList{ex("bla"), ex(2)}); + ExprContainer cont4 = Sum(sum); + REQUIRE(cont4->is()); + REQUIRE(cont4->as() == sum); + } + } + SECTION("Assignment") { + ExprContainer cont = Constant(1); + REQUIRE(cont->is()); + REQUIRE(cont->as() == Constant(1)); + + cont = Variable("bla"); + REQUIRE(cont->is()); + REQUIRE(cont->as() == Variable("bla")); + + Product prod(ExprPtrList{ex("bla"), ex(2)}); + cont = Product(prod); + REQUIRE(cont->is()); + REQUIRE(cont->as() == prod); + + Sum sum(ExprPtrList{ex("bla"), ex(2)}); + cont = Sum(sum); + REQUIRE(cont->is()); + REQUIRE(cont->as() == sum); + } + SECTION("value semantics") { + ExprContainer cont = Constant(1); + ExprContainer copy(cont); + copy = Variable("test"); + + REQUIRE(copy->is()); + REQUIRE(cont->is()); + } + SECTION("conversion to Expr &") { + bool passed1 = false; + + auto func1 = [&passed1](Expr &) { passed1 = true; }; + + ExprContainer expr = Constant(5); + func1(expr); + + REQUIRE(passed1); + + bool passed2 = false; + auto func2 = [&passed2](const Expr &) { passed2 = true; }; + + func2(std::as_const(expr)); + + REQUIRE(passed2); + + bool passed3 = false; + auto func3 = [&passed3](Expr &&) { passed3 = true; }; + + func3(std::move(expr)); + + REQUIRE(passed3); + } + SECTION("freestanding Expr arithmetic") { + // This allows to use arithmetic directly on Expr & instances (instead of + // requiring ExprPtr or ExprContainer wrappers) + ExprContainer res = Constant(1) + Variable("One"); + REQUIRE_THAT(res, EquivalentTo("1 + One")); + + res = Variable("A") * Tensor("T", bra({"a1"}), ket()) - Constant(42); + REQUIRE_THAT(res, EquivalentTo("A * T{a1} - 42")); + } + SECTION("In-place ExprContainer arithmetic") { + ExprContainer res = Constant(1); + res += Variable("A"); + res -= Variable("B"); + res *= Constant(3); + + REQUIRE_THAT(res, EquivalentTo("(1 + A - B) * 3")); + } + SECTION("Conversion to ExprPtr") { + ExprContainer cont = Constant(3); + + ExprPtr ptr = std::move(cont); + REQUIRE(ptr->is()); + REQUIRE(ptr->as().value() == 3); + + cont = Variable("A"); + // Copy-conversion-ctor is explicit + ptr = ExprPtr(cont); + REQUIRE(ptr->is()); + REQUIRE(ptr->as().label() == L"A"); + + // Ensure that ptr actually points to a copy + ptr->as().set_label(L"B"); + REQUIRE(cont->as().label() == L"A"); + REQUIRE(ptr->as().label() == L"B"); + } + SECTION("Conversion from ExprPtr") { + // This conversion is always explicit as it always has to perform a copy. + // Even a moved-from ExpPtr might point to an object that is co-owned by + // another ExprPtr and thus "resource stealing" is not possible. + ExprPtr ptr = ex(2); + ExprContainer cont(ptr); + REQUIRE(cont->is()); + REQUIRE(cont->as().value() == 2); + + ptr->as() = Constant(3); + REQUIRE(ptr->as().value() == 3); + REQUIRE(cont->as().value() == 2); + } + } + SECTION("ResultExpr") { SECTION("accessors") { SECTION("as_variable") {