From 7210f61eec28241fb2b57c01d86a4d1883f14ffb Mon Sep 17 00:00:00 2001 From: Alexander Condello Date: Wed, 26 Aug 2026 15:33:22 -0700 Subject: [PATCH 1/2] Add interval class --- .../include/dwave-optimization/interval.hpp | 195 ++++++++++++++++++ dwave/optimization/src/interval.cpp | 49 +++++ meson.build | 1 + .../notes/interval-ad3c6f2bd14d3747.yaml | 3 + tests/cpp/meson.build | 1 + tests/cpp/test_interval.cpp | 146 +++++++++++++ 6 files changed, 395 insertions(+) create mode 100644 dwave/optimization/include/dwave-optimization/interval.hpp create mode 100644 dwave/optimization/src/interval.cpp create mode 100644 releasenotes/notes/interval-ad3c6f2bd14d3747.yaml create mode 100644 tests/cpp/test_interval.cpp diff --git a/dwave/optimization/include/dwave-optimization/interval.hpp b/dwave/optimization/include/dwave-optimization/interval.hpp new file mode 100644 index 000000000..bee59eb32 --- /dev/null +++ b/dwave/optimization/include/dwave-optimization/interval.hpp @@ -0,0 +1,195 @@ +// Copyright 2026 D-Wave +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#pragma once + +#include +#include +#include +#include +#include + +#include "dwave-optimization/typing.hpp" + +namespace dwave::optimization { + +/// An interval encodes a range of possible values. +/// +/// Note that this class does not (yet) implement outward rounding. +template +struct interval { + /// Construct an empty interval. + constexpr interval() = default; + + /// Construct an interval of values between inf and sup (inclusive). + /// When ``sup < inf`` the interval is treated as empty. + constexpr interval(T inf, T sup) noexcept : infimum(inf), supremum(sup) {} + + /// Copy constructor. + interval(const interval&) = default; + + /// Create an ``interval`` from another interval. + template + requires(std::same_as) + interval(const interval& other) noexcept : interval(other.infimum, other.supremum) {} + // dev note: we could expand this. E.g., we could support all promotions + // allowed by NumPy promotion. + + /// Move constructor. + interval(interval&&) = default; + + /// Copy assignment operator. + interval& operator=(const interval&) = default; + + /// Move assignment operator. + interval& operator=(interval&&) = default; + + /// Destructor. + ~interval() = default; + + /// An interval evalutes to `true` if it is not empty. + explicit constexpr operator bool() const noexcept { return infimum <= supremum; } + + /// Two intervals are treated as equal if they are the same type and have the same endpoints + /// or if they are both null. + constexpr bool operator==(const interval& rhs) const { + if (not static_cast(*this) and not static_cast(rhs)) return true; // both null + return infimum == rhs.infimum and supremum == rhs.supremum; + } + // dev note: we could support other type combinations in the future + + /// Comparison operators <, <=, >=, > are used for strict subset, subset, superset, and strict + /// superset respectively. + constexpr std::partial_ordering operator<=>(const interval& rhs) const { + // If we're equal then we're equivalent + if (*this == rhs) return std::partial_ordering::equivalent; + + // If lhs != rhs then at most one can be empty + if (not static_cast(*this)) return std::partial_ordering::less; // empty < not-empty + if (not static_cast(rhs)) return std::partial_ordering::greater; // not-empty > empty + + // Ok, neither are empty + + // If lhs <= rhs and lhs != rhs then lhs < rhs + if (rhs.infimum <= infimum and supremum <= rhs.supremum) { + return std::partial_ordering::less; + } + + // If lhs >= rhs and lhs != rhs then lhs > rhs + if (infimum <= rhs.infimum and rhs.supremum <= supremum) { + return std::partial_ordering::greater; + } + + // Otherwise we're not comparable + return std::partial_ordering::unordered; + } + // dev note: we could support other type combinations in the future + + /// Negate and swap the values in the interval. + /// For boolean intervals, negation is treated as logical not. + constexpr interval operator-() const { + // For bool, we overload this to be negation + if constexpr (std::same_as) return interval(not supremum, not infimum); + + // -INT_MIN is undefined. Under the assumption that if the user is using + // INT_MIN/INT_MAX they probably are trying to say "unbounded" we do a + // weird thing and just define -INT_MIN := INT_MAX and -INT_MAX := INT_MIN + // even though that's wrong and leads to some slightly weird outcomes + if constexpr (std::integral) { + using limits = std::numeric_limits; + if (infimum == limits::lowest() and supremum == limits::max()) { + return *this; + } else if (infimum == limits::lowest()) { + return interval(-supremum, limits::max()); + } else if (supremum == limits::max()) { + return interval(limits::lowest(), -infimum); + } + } + + return interval(-supremum, -infimum); + } + + /// Intersection with ``rhs``. + constexpr interval& operator&=(const interval& rhs) { + // If lhs is an empty interval, then the intersection is just lhs + if (not static_cast(*this)) return *this; + + // If rhs is an empty interval, then the intersection is just rhs + if (not static_cast(rhs)) return *this = rhs; + + if (infimum < rhs.infimum) infimum = rhs.infimum; + if (rhs.supremum < supremum) supremum = rhs.supremum; + + return *this; + } + + /// Union with ``rhs``. + constexpr interval& operator|=(const interval& rhs) { + // If rhs is an empty interval, then taking the union with it does nothing + if (not static_cast(rhs)) return *this; + + // If lhs is an empty interval, then the union is just rhs + if (not static_cast(*this)) return *this = rhs; + + if (rhs.infimum < infimum) infimum = rhs.infimum; + if (supremum < rhs.supremum) supremum = rhs.supremum; + + return *this; + } + + /// Interection of two intervals. + friend constexpr interval operator&(interval lhs, const interval& rhs) { + lhs &= rhs; + return lhs; + } + + /// Union of two intervals + friend constexpr interval operator|(interval lhs, const interval& rhs) { + lhs |= rhs; + return lhs; + } + + /// The maximum expressible interval + static consteval interval all() { + using limits = std::numeric_limits; + if constexpr (limits::has_infinity) { + return interval(-limits::infinity(), limits::infinity()); + } else { + return interval(limits::lowest(), limits::max()); + } + } + + /// Test whether `x` is a value in the interval. + constexpr bool contains(const T& x) const { return infimum <= x and x <= supremum; } + // dev note: we could support other type combinations in the future + + /// All expressible non-negative values. + static consteval interval nonnegative() { + using limits = std::numeric_limits; + if constexpr (limits::has_infinity) { + return interval(0, limits::infinity()); + } else { + return interval(0, limits::max()); + } + } + + T infimum = 1; + T supremum = 0; +}; + +// Intervals are printable +template +std::ostream& operator<<(std::ostream& os, const interval& in); + +} // namespace dwave::optimization diff --git a/dwave/optimization/src/interval.cpp b/dwave/optimization/src/interval.cpp new file mode 100644 index 000000000..07c22f191 --- /dev/null +++ b/dwave/optimization/src/interval.cpp @@ -0,0 +1,49 @@ +// Copyright 2026 D-Wave +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "dwave-optimization/interval.hpp" + +#include + +namespace dwave::optimization { + +template +std::ostream& operator<<(std::ostream& os, const interval& in) { + if (not static_cast(in)) return os << "[]"; + + os << "["; + + // Not all compilers print all dtypes. So coerce them into a smaller set of + // possible types + if constexpr (std::integral) { + os << static_cast(in.infimum) << ", " << static_cast(in.supremum); + } else if constexpr (std::floating_point) { + os << static_cast(in.infimum) << ", " << static_cast(in.supremum); + } else { + assert(false and "unexpected dtype"); + } + + os << "]"; + return os; +} + +template std::ostream& operator<<(std::ostream&, const interval&); +template std::ostream& operator<<(std::ostream&, const interval&); +template std::ostream& operator<<(std::ostream&, const interval&); +template std::ostream& operator<<(std::ostream&, const interval&); +template std::ostream& operator<<(std::ostream&, const interval&); +template std::ostream& operator<<(std::ostream&, const interval&); +template std::ostream& operator<<(std::ostream&, const interval&); + +} // namespace dwave::optimization diff --git a/meson.build b/meson.build index 1e64bf33c..f93631a1b 100644 --- a/meson.build +++ b/meson.build @@ -55,6 +55,7 @@ dwave_optimization_src = [ 'dwave/optimization/src/array.cpp', 'dwave/optimization/src/fraction.cpp', 'dwave/optimization/src/graph.cpp', + 'dwave/optimization/src/interval.cpp', 'dwave/optimization/src/simplex.cpp', ] diff --git a/releasenotes/notes/interval-ad3c6f2bd14d3747.yaml b/releasenotes/notes/interval-ad3c6f2bd14d3747.yaml new file mode 100644 index 000000000..84d14e4de --- /dev/null +++ b/releasenotes/notes/interval-ad3c6f2bd14d3747.yaml @@ -0,0 +1,3 @@ +--- +features: + - Add a simple C++ ``dwave::optimization::interval`` class. diff --git a/tests/cpp/meson.build b/tests/cpp/meson.build index 3f1421faa..7b9041fe2 100644 --- a/tests/cpp/meson.build +++ b/tests/cpp/meson.build @@ -36,6 +36,7 @@ tests_all = executable( 'test_functional.cpp', 'test_functional_.cpp', 'test_graph.cpp', + 'test_interval.cpp', 'test_iterators.cpp', 'test_simplex.cpp', 'test_type_list.cpp', diff --git a/tests/cpp/test_interval.cpp b/tests/cpp/test_interval.cpp new file mode 100644 index 000000000..5fb5c5b89 --- /dev/null +++ b/tests/cpp/test_interval.cpp @@ -0,0 +1,146 @@ +// Copyright 2026 D-Wave +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include +#include + +#include "dwave-optimization/interval.hpp" + +namespace dwave::optimization { + +TEMPLATE_LIST_TEST_CASE("interval", "", DTypes) { + SECTION("::operator bool()") { + STATIC_REQUIRE(not interval()); + STATIC_REQUIRE(interval(0, 1)); + STATIC_REQUIRE(interval(0, 0)); + STATIC_REQUIRE(interval(1, 1)); + + // We only allow explicit conversion to bool + STATIC_REQUIRE(not std::convertible_to, bool>); + } + + SECTION("::operator==") { + STATIC_REQUIRE(interval(0, 1) == interval(0, 1)); + + STATIC_REQUIRE(interval(0, 0) != interval(0, 1)); + STATIC_REQUIRE(interval(0, 1) != interval()); + + STATIC_REQUIRE(interval(10, -10) == interval(1, 0)); // null always equals null + } + + SECTION("::operator<= (i.e., subset)") { + STATIC_REQUIRE(interval(0, 0) <= interval(0, 1)); + STATIC_REQUIRE(interval(0, 1) <= interval(0, 1)); // equality allowed + } + + SECTION("::operator-") { + if constexpr (std::same_as) { + STATIC_REQUIRE(interval(0, 1) == -interval(0, 1)); + STATIC_REQUIRE(interval(0, 0) == -interval(1, 1)); + STATIC_REQUIRE(interval(1, 1) == -interval(0, 0)); + } else { + STATIC_REQUIRE(interval(1, 10) == -interval(-10, -1)); + } + + if constexpr (std::integral) { + STATIC_REQUIRE(-interval::all() == interval::all()); + } + } + + SECTION("::operator&=/::operator& (i.e., intersection)") { + STATIC_REQUIRE(not static_cast(interval(0, 1) & interval())); + STATIC_REQUIRE(not static_cast(interval() & interval(0, 1))); + + if constexpr (not std::same_as) { + STATIC_REQUIRE( + (interval(0, 5) & interval(2, 3)) == interval(2, 3) + ); + STATIC_REQUIRE( + (interval(0, 5) & interval(3, 10)) == interval(3, 5) + ); + } + } + + SECTION("::operator|=/::operator| (i.e., union)") { + STATIC_REQUIRE( + (interval(0, 1) | interval()) == interval(0, 1) + ); + STATIC_REQUIRE( + (interval() | interval(0, 1)) == interval(0, 1) + ); + + if constexpr (not std::same_as) { + STATIC_REQUIRE( + (interval(0, 5) | interval(2, 3)) == interval(0, 5) + ); + STATIC_REQUIRE( + (interval(0, 5) | interval(3, 10)) == interval(0, 10) + ); + } + } + + SECTION("::contains") { + STATIC_REQUIRE(not interval().contains(0)); + STATIC_REQUIRE(interval(0, 1).contains(0)); + STATIC_REQUIRE(interval(1, 1).contains(1)); + } + + SECTION("printing") { + SECTION("integral") { + std::stringstream ss; + ss << interval(0, 1); + CHECK(ss.str() == "[0, 1]"); + } + + SECTION("floating") { + if constexpr (std::floating_point) { + std::stringstream ss; + ss << interval(.5, 1.5); + CHECK(ss.str() == "[0.5, 1.5]"); + } + } + } + + SECTION("structured binding") { + SECTION("const reference") { + auto in = interval(0, 1); + const auto& [inf, sup] = in; + STATIC_REQUIRE(std::same_as); + STATIC_REQUIRE(std::same_as); + CHECK(inf == 0); + CHECK(sup == 1); + } + + SECTION("rvalue") { + auto in = interval(0, 1); + auto [inf, sup] = in; + STATIC_REQUIRE(std::same_as); + STATIC_REQUIRE(std::same_as); + CHECK(inf == 0); + CHECK(sup == 1); + } + + SECTION("const rvalue") { + const auto in = interval(0, 1); + + auto&& [inf, sup] = in; + STATIC_REQUIRE(std::same_as); + STATIC_REQUIRE(std::same_as); + CHECK(inf == 0); + CHECK(sup == 1); + } + } +} + +} // namespace dwave::optimization From 1d53fd18a8c4b48e212a91d0b049c3b9fafcf6da Mon Sep 17 00:00:00 2001 From: Alexander Condello Date: Thu, 27 Aug 2026 13:25:53 -0700 Subject: [PATCH 2/2] Rework unaryops to carry additional information and UnaryOpNode to use them --- .../include/dwave-optimization/functional.hpp | 276 ++++++++++-- .../dwave-optimization/nodes/unaryop.hpp | 298 +++++++++++-- dwave/optimization/src/nodes/unaryop.cpp | 293 +------------ .../unaryop-rework-23877cf721e8f8a6.yaml | 8 + tests/cpp/nodes/test_unaryop.cpp | 188 ++++---- tests/cpp/test_functional.cpp | 402 +++++++++++++++++- 6 files changed, 1009 insertions(+), 456 deletions(-) create mode 100644 releasenotes/notes/unaryop-rework-23877cf721e8f8a6.yaml diff --git a/dwave/optimization/include/dwave-optimization/functional.hpp b/dwave/optimization/include/dwave-optimization/functional.hpp index e050a7c11..ab87f5efb 100644 --- a/dwave/optimization/include/dwave-optimization/functional.hpp +++ b/dwave/optimization/include/dwave-optimization/functional.hpp @@ -15,45 +15,185 @@ #pragma once #include +#include #include #include #include +#include +#include + +#include "dwave-optimization/interval.hpp" +#include "dwave-optimization/typing.hpp" namespace dwave::optimization::functional { -template -struct abs { - static constexpr T operator()(const T& x) { return std::abs(x); } +enum class Monotonicity { Decreasing = -1, None = 0, Increasing = 1 }; + +template +struct UnaryOpMixin { + template + requires(UnaryOp::monotonic != Monotonicity::None) + static auto operator()(const interval& domain) { + using return_type = interval; + + // op(empty domain) -> empty domain + if (not static_cast(domain)) return return_type(); + + assert( + domain <= UnaryOp::template domain and + "input domain must be a subset of the func's domain" + ); + + // We don't worry about outward rounding here because this overload is meant + // to reflect the behavior of the scalar overload, not necessarily to be + // mathematically correct. + // We *do* assume that UnaryOp (e.g., std::exp()) is monotonic, which + // is not always true, but I think it's an OK assumption for our purposes. + if constexpr (UnaryOp::monotonic == Monotonicity::Increasing) { + return return_type( + UnaryOp::operator()(domain.infimum), UnaryOp::operator()(domain.supremum) + ); + } else if constexpr (UnaryOp::monotonic == Monotonicity::Decreasing) { + return return_type( + UnaryOp::operator()(domain.supremum), UnaryOp::operator()(domain.infimum) + ); + } else { + assert(false and "unexpected monotonicity"); + std::unreachable(); + } + } + + template + static constexpr interval domain = interval::all(); }; -template -struct cos { - static auto operator()(const T& num) { return std::cos(num); } +struct absolute : UnaryOpMixin { + template + static T operator()(const T& x) { + // Unlike NumPy/std, we define std::abs(INT_MIN) to equal INT_MAX under the reasoning + // that it's more important to us to preserve the sign than to preseve the correct value. + if constexpr (std::integral) { + if (x == std::numeric_limits::lowest()) return std::numeric_limits::max(); + } + + // std::abs() is not defined for int8 or int16 so we static_cast to avoid widening. + return static_cast(std::abs(x)); + } + static bool operator()(const bool& x) { return x; } + + template + static interval operator()(const interval& domain) { + if (not static_cast(domain)) return {}; // op(empty domain) -> empty domain + + assert(domain.infimum <= domain.supremum); // implied by non-empty + + // If the domain is non-negative, then absolute is identity + if (0 <= domain.infimum) return domain; + + // If the domain is negative, then absolute is just the inverse + if (domain.supremum < 0) return -domain; + + // Otherwise, the domain straddles 0 + + // Handle the -INT_MIN case. Again we treat abs(-INT_MIN) as INT_MAX under the reasoning + // that [INT_MIN, ...] is probably intended to mean unbounded. + if constexpr (std::integral) { + if (domain.infimum == std::numeric_limits::lowest()) { + return interval(0, std::numeric_limits::max()); + } + } + + return interval( + 0, -domain.infimum < domain.supremum ? domain.supremum : -domain.infimum + ); + } + static interval operator()(const interval& domain) { return domain; } + + static constexpr Monotonicity monotonic = Monotonicity::None; }; -template -struct exp { - static constexpr auto operator()(const T& x) { return std::exp(x); } +struct cos : UnaryOpMixin { + static auto operator()(const DType auto& x) { return std::cos(x); } + + template + static interval operator()(const interval& domain) { + if (not static_cast(domain)) return {}; // op(empty domain) -> empty domain + + // It is possible to be a lot more specific than this by checking whether + // our domain spans a full period or not, but I think this is of dubious + // benefit to the user so for now we just return [-1, +1] + return {-1, +1}; + } + + static constexpr Monotonicity monotonic = Monotonicity::None; }; -template -struct expit { - static constexpr double operator()(const T& x) { return 1.0 / (1.0 + std::exp(-1. * x)); } +struct exp : UnaryOpMixin { + static auto operator()(const DType auto& x) { return std::exp(x); } + using UnaryOpMixin::operator(); + + static constexpr Monotonicity monotonic = Monotonicity::Increasing; }; -template -struct log { - static constexpr auto operator()(const T& x) { return std::log(x); } +struct expit : UnaryOpMixin { + template + static auto operator()(const T& x) { + return 1 / (1 + std::exp(-x)); + } + using UnaryOpMixin::operator(); + + static constexpr Monotonicity monotonic = Monotonicity::Increasing; }; -template -struct logical { - static constexpr bool operator()(const T& x) { return x; } +struct log : UnaryOpMixin { + template + static auto operator()(const T& x) { + assert(domain.contains(x) and "x must be non-negative"); + return std::log(x); + } + using UnaryOpMixin::operator(); + + template + static constexpr interval domain = interval::nonnegative(); + + static constexpr Monotonicity monotonic = Monotonicity::Increasing; +}; + +struct logical : UnaryOpMixin { + static bool operator()(const DType auto& x) { return x; } + + static interval operator()(const interval& domain) { return domain; } + template + static interval operator()(const interval& domain) { + if (not static_cast(domain)) return {}; // op(empty domain) -> empty domain + + if (domain.infimum == 0 and domain.supremum == 0) return interval(false, false); + if (domain.infimum <= 0 and domain.supremum >= 0) return interval(false, true); + return interval(true, true); + } + + static constexpr Monotonicity monotonic = Monotonicity::None; +}; + +struct logical_not : UnaryOpMixin { + static bool operator()(const DType auto& x) { return not x; } + + static interval operator()(const interval& domain) { + if (not static_cast(domain)) return {}; // op(empty domain) -> empty domain + return interval(not domain.supremum, not domain.infimum); + } + template + static interval operator()(const interval& domain) { + // Call the more specific interval overload + return operator()(logical{}(domain)); + } + + static constexpr Monotonicity monotonic = Monotonicity::None; }; template struct logical_xor { - static constexpr bool operator()(const T& x, const T& y) { + static bool operator()(const T& x, const T& y) { return static_cast(x) != static_cast(y); } }; @@ -90,9 +230,28 @@ struct modulus { } }; -template -struct rint { - static constexpr auto operator()(const T& x) { return std::rint(x); } +struct negate : UnaryOpMixin { + template + requires(DType and not std::same_as) // not defined for bool + static auto operator()(const T& x) { + // We define -INT_MIN to equal INT_MAX under the reasoning that it's more + // important to us to preserve the sign than to preseve the correct value. + if constexpr (std::integral) { + if (x == std::numeric_limits::lowest()) return std::numeric_limits::max(); + } + + return static_cast(-x); // so it doesn't widen e.g., int8_t->int + } + using UnaryOpMixin::operator(); + + static constexpr Monotonicity monotonic = Monotonicity::Decreasing; +}; + +struct rint : UnaryOpMixin { + static auto operator()(const DType auto& x) { return std::rint(x); } + using UnaryOpMixin::operator(); + + static constexpr Monotonicity monotonic = Monotonicity::Increasing; }; template @@ -103,24 +262,73 @@ struct safe_divides { } }; -template -struct sin { - static auto operator()(const T& num) { return std::sin(num); } +struct sin : UnaryOpMixin { + static auto operator()(const DType auto& x) { return std::sin(x); } + + template + static interval operator()(const interval& domain) { + if (not static_cast(domain)) return {}; // op(empty domain) -> empty domain + + // It is possible to be a lot more specific than this by checking whether + // our domain spans a full period or not, but I think this is of dubious + // benefit to the user so for now we just return [-1, +1] + return {-1, +1}; + } + + static constexpr Monotonicity monotonic = Monotonicity::None; }; -template -struct square { - static constexpr T operator()(const T& x) { return x * x; } +struct square : UnaryOpMixin { + template + static T operator()(const T& x) { + return x * x; + } + static bool operator()(const bool& x) { return x; } + + template + static interval operator()(const interval& domain) { + if (not static_cast(domain)) return {}; // op(empty domain) -> empty domain + + assert(domain.infimum <= domain.supremum); // implied by non-empty + + square op{}; + T inf_squared = op(domain.infimum); + T sup_squared = op(domain.supremum); + + // Non-negative domain: square is increasing + if (0 <= domain.infimum) return interval(inf_squared, sup_squared); + + // Non-positive domain: square is decreasing + if (domain.supremum <= 0) return interval(sup_squared, inf_squared); + + // Otherwise the domain straddles 0: minimum is 0, maximum is the larger squared endpoint. + + return interval(0, inf_squared < sup_squared ? sup_squared : inf_squared); + } + static interval operator()(const interval& domain) { return domain; } + + static constexpr Monotonicity monotonic = Monotonicity::None; }; -template -struct square_root { - static constexpr auto operator()(const T& x) { return std::sqrt(x); } +struct square_root : UnaryOpMixin { + template + static auto operator()(const T& x) { + assert(domain.contains(x) and "x must be non-negative"); + return std::sqrt(x); + } + using UnaryOpMixin::operator(); + + template + static constexpr interval domain = interval::nonnegative(); + + static constexpr Monotonicity monotonic = Monotonicity::Increasing; }; -template -struct tanh { - static auto operator()(const T& num) { return std::tanh(num); } +struct tanh : UnaryOpMixin { + static auto operator()(const DType auto& num) { return std::tanh(num); } + using UnaryOpMixin::operator(); + + static constexpr Monotonicity monotonic = Monotonicity::Increasing; }; } // namespace dwave::optimization::functional diff --git a/dwave/optimization/include/dwave-optimization/nodes/unaryop.hpp b/dwave/optimization/include/dwave-optimization/nodes/unaryop.hpp index 7e2795339..a60fc9d79 100644 --- a/dwave/optimization/include/dwave-optimization/nodes/unaryop.hpp +++ b/dwave/optimization/include/dwave-optimization/nodes/unaryop.hpp @@ -15,51 +15,82 @@ #pragma once #include -#include -#include +#include +#include +#include #include +#include #include +#include #include "dwave-optimization/array.hpp" #include "dwave-optimization/functional.hpp" #include "dwave-optimization/graph.hpp" +#include "dwave-optimization/interval.hpp" #include "dwave-optimization/state.hpp" namespace dwave::optimization { +/// Apply a unary function to the predecessor node. template class UnaryOpNode : public ArrayOutputMixin { public: - explicit UnaryOpNode(ArrayNode* node_ptr); + template + explicit UnaryOpNode(ArrayNode* node_ptr, Args&&... args) : + ArrayOutputMixin(node_ptr->shape()), + op(std::forward(args)...), + array_ptr_(node_ptr), + domain_(calculate_domain_(node_ptr)), + integral_(calculate_integral_(node_ptr)), + sizeinfo_(array_ptr_->sizeinfo()) { + add_predecessor_(node_ptr); + } - double const* buff(const State& state) const override; - std::span diff(const State& state) const override; + /// @copydoc Array::buff() + double const* buff(const State& state) const override { + return data_ptr_(state)->buffer.data(); + } - bool equal_to(const Node& rhs) const override; - bool equal_to(const UnaryOpNode& rhs) const; + /// @copydoc Node::commit() + void commit(State& state) const override { + auto* state_ptr = data_ptr_(state); + state_ptr->diff.clear(); + state_ptr->previous_size = state_ptr->buffer.size(); + } - /// @copydoc Array::integral() - bool integral() const override; + /// @copydoc Array::diff() + std::span diff(const State& state) const override { + return data_ptr_(state)->diff; + } - /// @copydoc Array::min() - double min() const override; + /// @copydoc Node::equal_to() + bool equal_to(const Node& rhs) const override { + const UnaryOpNode* rhs_ptr = dynamic_cast(&rhs); + if (rhs_ptr == nullptr) return false; // not same type so not equal + return this->equal_to(*rhs_ptr); // use the equal_to(const UnaryOpNode&) overload + } + bool equal_to(const UnaryOpNode& rhs) const { + // If we're the same type, then we just need to check we have the same predecessor + return this->array_ptr_ == rhs.array_ptr_; + } - /// @copydoc Array::max() - double max() const override; + /// @copydoc Node::initialize_state() + void initialize_state(State& state) const override { + emplace_data_ptr_( + state, array_ptr_->view(state) | std::views::transform(op), array_ptr_->shape(state) + ); + } - using ArrayOutputMixin::shape; - std::span shape(const State& state) const override; - using ArrayOutputMixin::size; - ssize_t size(const State& state) const override; - ssize_t size_diff(const State& state) const override; - SizeInfo sizeinfo() const override; + /// @copydoc Array::integral() + bool integral() const override { return integral_; } - void commit(State& state) const override; - void revert(State& state) const override; - void initialize_state(State& state) const override; - void propagate(State& state) const override; + /// @copydoc Array::max() + double max() const override { return domain_.supremum; } - // The predecessor of the operation, as an Array*. + /// @copydoc Array::min() + double min() const override { return domain_.infimum; } + + /// The predecessor of the operation std::span operands() { assert(predecessors().size() == 1); return std::span(&array_ptr_, 1); @@ -69,8 +100,170 @@ class UnaryOpNode : public ArrayOutputMixin { return std::span(&array_ptr_, 1); } + /// @copydoc Node::propagate() + void propagate(State& state) const override { + const auto array_diff = array_ptr_->diff(state); + + // If there are no updates the handle, return early. + if (array_diff.empty()) return; + + auto* state_ptr = data_ptr_(state); + auto& buffer = state_ptr->buffer; + auto& diff = state_ptr->diff; + + for (const auto& update : array_diff) { + assert(0 <= update.index); + + if (update.removed()) { + assert(static_cast(update.index) + 1 == buffer.size()); + + diff.emplace_back(Update::removal(update.index, buffer[update.index])); + buffer.pop_back(); + } else if (update.placed()) { + assert(static_cast(update.index) == buffer.size()); + + buffer.emplace_back(op(update.value)); + diff.emplace_back(Update::placement(update.index, buffer.back())); + } else { + assert(static_cast(update.index) < buffer.size()); + + double& old = buffer[update.index]; + double value = op(update.value); + + if (old == value) continue; // no change to update + + diff.emplace_back(update.index, old, value); + old = value; + } + } + + if (ndim()) state_ptr->shape[0] = array_ptr_->shape(state)[0]; + + if (not diff.empty()) Node::propagate(state); + } + + /// @copydoc Node::revert() + void revert(State& state) const override { + auto* state_ptr = data_ptr_(state); + std::vector& buffer = state_ptr->buffer; + const ssize_t size = state_ptr->previous_size; + std::vector& diff = state_ptr->diff; + ssize_t* shape = state_ptr->shape.get(); + + const ssize_t propagated_size = buffer.size(); + + buffer.resize(size); + for (const auto& [index, old, _] : diff | std::views::reverse) { + assert(0 <= index); + if (size <= index) continue; + buffer[index] = old; + } + diff.clear(); + + // Adjust our shape back to what it should be, while avoiding reading our + // predecessor. + if (size != propagated_size) { + assert(0 < this->ndim()); // if our size changed this must be true + shape[0] = buffer.size(); + + // we could cache the divisor, but that feels like overkill in the context + // of a revert + for (ssize_t i = 1, ndim = this->ndim(); i < ndim; ++i) shape[0] /= shape[i]; + } + } + + /// @copydoc Array::shape() + std::span shape(const State& state) const override { + return std::span(data_ptr_(state)->shape.get(), ndim()); + } + using ArrayOutputMixin::shape; + + /// @copydoc Array::size() + ssize_t size(const State& state) const override { + return data_ptr_(state)->buffer.size(); + } + using ArrayOutputMixin::size; + + /// @copydoc Array::size_diff() + ssize_t size_diff(const State& state) const override { + const auto* state_ptr = data_ptr_(state); + return state_ptr->buffer.size() - state_ptr->previous_size; + } + + /// @copydoc Array::sizeinfo() + SizeInfo sizeinfo() const override { return this->sizeinfo_; } + private: - void replace_predecessor_(ssize_t index, Node* node_ptr) override; + // It would be more convenient to use ArrayStateData from _state.hpp + // but we don't currently have that in a public header. + // If that changes in the future, or we add default state support + // https://github.com/dwavesystems/dwave-optimization/issues/629 then + // we should use that. + struct UnaryOpState_ : NodeStateData { + UnaryOpState_() = delete; + + template + UnaryOpState_(R&& values, std::span shape) : + buffer(std::ranges::to>(std::forward(values))), + previous_size(buffer.size()), + shape(shape.size() ? std::make_unique(shape.size()) : nullptr), + diff() { + for (ssize_t i = 0, ndim = shape.size(); i < ndim; ++i) { + this->shape[i] = shape[i]; + } + } + + std::vector buffer; + ssize_t previous_size; + + std::unique_ptr shape; + + std::vector diff; + }; + + // Calculate the min/max of the op applied to array_ptr + static interval calculate_domain_(Array* array_ptr) { + interval array_domain = interval(array_ptr->min(), array_ptr->max()); + + if constexpr (std::same_as) { + if (not(array_domain <= UnaryOp::template domain)) { + throw std::invalid_argument("SquareRoot's predecessors must be non-negative"); + } + } + if constexpr (std::same_as) { + if (not(array_domain <= UnaryOp::template domain)) { + throw std::invalid_argument("Log's predecessors must be non-negative"); + } + } + + // All the other ufuncs should have unbounded domains but as a sanity check... + assert(array_domain <= UnaryOp::template domain); + + // The range of our UnaryOp, i.e. the output min/max of our UnaryOpNode + return UnaryOp{}(array_domain); + } + + // Determine whether the op applied to array_ptr will always result in an + // integral output. + static bool calculate_integral_(const Array* array_ptr) { + // rint() actually always returns a floating point. But it is an intergral + // floating point. This is a place where having proper dtypes would be + // very very nice. + if constexpr (std::same_as) return true; + + // Otherwise, we ask our op what an integer input would result in. + if (array_ptr->integral()) { + return std::integral; + } else { + return std::integral; + } + } + + void replace_predecessor_(ssize_t index, Node* node_ptr) override { + Node::replace_predecessor_(index, node_ptr); + array_ptr_ = dynamic_cast(node_ptr); + assert(array_ptr_ != nullptr); + } UnaryOp op; @@ -78,22 +271,49 @@ class UnaryOpNode : public ArrayOutputMixin { // predecessors(), but this is more performant ArrayNode* array_ptr_; - const ValuesInfo values_info_; + interval domain_; // the range of possible output values + bool integral_; // whether the node will always output integral values + const SizeInfo sizeinfo_; }; -using AbsoluteNode = UnaryOpNode>; -using CosNode = UnaryOpNode>; -using ExpitNode = UnaryOpNode>; -using ExpNode = UnaryOpNode>; -using LogNode = UnaryOpNode>; -using LogicalNode = UnaryOpNode>; -using NegativeNode = UnaryOpNode>; -using NotNode = UnaryOpNode>; -using RintNode = UnaryOpNode>; -using SinNode = UnaryOpNode>; -using SquareNode = UnaryOpNode>; -using SquareRootNode = UnaryOpNode>; -using TanhNode = UnaryOpNode>; +using AbsoluteNode = UnaryOpNode; +extern template class UnaryOpNode; + +using CosNode = UnaryOpNode; +extern template class UnaryOpNode; + +using ExpNode = UnaryOpNode; +extern template class UnaryOpNode; + +using ExpitNode = UnaryOpNode; +extern template class UnaryOpNode; + +using LogNode = UnaryOpNode; +extern template class UnaryOpNode; + +using LogicalNode = UnaryOpNode; +extern template class UnaryOpNode; + +using NegativeNode = UnaryOpNode; +extern template class UnaryOpNode; + +using NotNode = UnaryOpNode; +extern template class UnaryOpNode; + +using RintNode = UnaryOpNode; +extern template class UnaryOpNode; + +using SinNode = UnaryOpNode; +extern template class UnaryOpNode; + +using SquareNode = UnaryOpNode; +extern template class UnaryOpNode; + +using SquareRootNode = UnaryOpNode; +extern template class UnaryOpNode; + +using TanhNode = UnaryOpNode; +extern template class UnaryOpNode; } // namespace dwave::optimization diff --git a/dwave/optimization/src/nodes/unaryop.cpp b/dwave/optimization/src/nodes/unaryop.cpp index d871be910..f259f1b0e 100644 --- a/dwave/optimization/src/nodes/unaryop.cpp +++ b/dwave/optimization/src/nodes/unaryop.cpp @@ -14,287 +14,20 @@ #include "dwave-optimization/nodes/unaryop.hpp" -#include - -#include "_state.hpp" - namespace dwave::optimization { -template -std::pair calculate_values_minmax(const Array* array_ptr) { - // Do some checks to make sure the resulting domain/range will be valid - if constexpr (std::is_same>::value) { - if (array_ptr->min() < 0) { - throw std::invalid_argument("SquareRoot's predecessors cannot take a negative value"); - } - } else if constexpr (std::is_same>::value) { - if (array_ptr->min() <= 0) { - throw std::invalid_argument("Log's predecessors cannot take a negative or zero value"); - } - } - - // If the output of the operation is boolean, then don't bother caching the result. - using result_type = typename std::invoke_result::type; - if constexpr (std::same_as) { - return {false, true}; - } - - // Likewise for sin/cos/tanh the minmax is -1/+1. We could tighten it if the domain - // of our predecessor is smaller than 2pi, but let's keep it simple for now - if constexpr ( - std::same_as> || - std::same_as> || - std::same_as> - ) { - return {-1, +1}; - } - - // Otherwise the min and max depend on the predecessor - - auto low = array_ptr->min(); - auto high = array_ptr->max(); - assert(low <= high); - - if constexpr (std::same_as>) { - if (low >= 0 && high >= 0) { - return std::make_pair(low, high); - } else if (low >= 0) { - assert(false && "min > max"); - std::unreachable(); - } else if (high >= 0) { - return std::pair(0.0, std::max(-low, high)); - } else { - return std::make_pair(-high, -low); - } - } - if constexpr (std::same_as>) { - return std::make_pair(std::exp(low), std::exp(high)); - } - if constexpr (std::same_as>) { - double expit_low = 1.0 / (1.0 + std::exp(-low)); - double expit_high = 1.0 / (1.0 + std::exp(-high)); - return std::make_pair(expit_low, expit_high); - } - if constexpr (std::same_as>) { - assert(low > 0); // checked by constructor - return std::make_pair(std::log(low), std::log(high)); - } - if constexpr (std::same_as>) { - return std::make_pair(std::rint(low), std::rint(high)); - } - if constexpr (std::same_as>) { - const auto highest = std::numeric_limits::max(); - return std::make_pair( - std::min({low * low, high * high, highest}), - std::min( - std::max({low * low, high * high}), - highest - ) - ); // prevent inf - } - if constexpr (std::same_as>) { - assert(low >= 0); // checked by constructor - return std::make_pair(std::sqrt(low), std::sqrt(high)); - } - if constexpr (std::same_as>) { - return std::make_pair(-high, -low); - } - - assert(false && "not implemeted yet"); - std::unreachable(); -} - -template -bool calculate_integral(const Array*) { - using result_type = typename std::invoke_result::type; - return std::is_integral::value; -} - -template <> -bool calculate_integral>(const Array* array_ptr) { - return array_ptr->integral(); -} - -template <> -bool calculate_integral>(const Array*) { - return false; -} - -template <> -bool calculate_integral>(const Array*) { - return false; -} - -template <> -bool calculate_integral>(const Array*) { - return false; -} - -template <> -bool calculate_integral>(const Array*) { - return false; -} - -template <> -bool calculate_integral>(const Array* array_ptr) { - return array_ptr->integral(); -} - -template <> -bool calculate_integral>(const Array*) { - return true; -} - -template <> -bool calculate_integral>(const Array*) { - return false; -} - -template <> -bool calculate_integral>(const Array* array_ptr) { - return array_ptr->integral(); -} - -template <> -bool calculate_integral>(const Array*) { - return false; -} - -template -UnaryOpNode::UnaryOpNode(ArrayNode* node_ptr) : - ArrayOutputMixin(node_ptr->shape()), - array_ptr_(node_ptr), - values_info_( - calculate_values_minmax(array_ptr_), - calculate_integral(array_ptr_) - ), - sizeinfo_(array_ptr_->sizeinfo()) { - add_predecessor_(node_ptr); -} - -template -void UnaryOpNode::commit(State& state) const { - data_ptr_(state)->commit(); -} - -template -double const* UnaryOpNode::buff(const State& state) const { - return data_ptr_(state)->buff(); -} - -template -std::span UnaryOpNode::diff(const State& state) const { - return data_ptr_(state)->diff(); -} - -template -bool UnaryOpNode::equal_to(const Node& rhs) const { - const UnaryOpNode* rhs_ptr = dynamic_cast(&rhs); - if (rhs_ptr == nullptr) return false; // not same type so not equal - return this->equal_to(*rhs_ptr); // use the equal_to(const UnaryOpNode&) overload -} - -template -bool UnaryOpNode::equal_to(const UnaryOpNode& rhs) const { - // If we're the same type, then we just need to check we have the same predecessor - return this->array_ptr_ == rhs.array_ptr_; -} - -template -void UnaryOpNode::initialize_state(State& state) const { - std::vector values; - values.reserve(array_ptr_->size(state)); - for (const double val : array_ptr_->view(state)) { - values.emplace_back(op(val)); - } - - emplace_data_ptr_(state, std::move(values)); -} - -template -bool UnaryOpNode::integral() const { - return values_info_.integral; -} - -template -double UnaryOpNode::min() const { - return this->values_info_.min; -} - -template -double UnaryOpNode::max() const { - return this->values_info_.max; -} - -template -void UnaryOpNode::propagate(State& state) const { - const auto diff = array_ptr_->diff(state); - // If there are no updates the handle, return early. - if (diff.empty()) return; - - auto node_data = data_ptr_(state); - - for (const auto& update : diff) { - const auto& [idx, _, value] = update; - - if (update.placed()) { - assert(idx == static_cast(node_data->size())); - node_data->emplace_back(op(value)); - } else if (update.removed()) { - assert(idx == static_cast(node_data->size()) - 1); - node_data->pop_back(); - } else { - node_data->set(idx, op(value)); - } - } - - if (node_data->diff().size()) Node::propagate(state); -} - -template -void UnaryOpNode::replace_predecessor_(ssize_t previous_index, Node* node_ptr) { - Node::replace_predecessor_(previous_index, node_ptr); - array_ptr_ = dynamic_cast(node_ptr); - assert(array_ptr_ != nullptr); -} - -template -void UnaryOpNode::revert(State& state) const { - data_ptr_(state)->revert(); -} - -template -std::span UnaryOpNode::shape(const State& state) const { - return array_ptr_->shape(state); -} - -template -ssize_t UnaryOpNode::size(const State& state) const { - return data_ptr_(state)->size(); -} - -template -ssize_t UnaryOpNode::size_diff(const State& state) const { - return data_ptr_(state)->size_diff(); -} - -template -SizeInfo UnaryOpNode::sizeinfo() const { - return this->sizeinfo_; -} - -template class UnaryOpNode>; -template class UnaryOpNode>; -template class UnaryOpNode>; -template class UnaryOpNode>; -template class UnaryOpNode>; -template class UnaryOpNode>; -template class UnaryOpNode>; -template class UnaryOpNode>; -template class UnaryOpNode>; -template class UnaryOpNode>; -template class UnaryOpNode>; -template class UnaryOpNode>; -template class UnaryOpNode>; +template class UnaryOpNode; +template class UnaryOpNode; +template class UnaryOpNode; +template class UnaryOpNode; +template class UnaryOpNode; +template class UnaryOpNode; +template class UnaryOpNode; +template class UnaryOpNode; +template class UnaryOpNode; +template class UnaryOpNode; +template class UnaryOpNode; +template class UnaryOpNode; +template class UnaryOpNode; } // namespace dwave::optimization diff --git a/releasenotes/notes/unaryop-rework-23877cf721e8f8a6.yaml b/releasenotes/notes/unaryop-rework-23877cf721e8f8a6.yaml new file mode 100644 index 000000000..3827288b8 --- /dev/null +++ b/releasenotes/notes/unaryop-rework-23877cf721e8f8a6.yaml @@ -0,0 +1,8 @@ +--- +features: + - Make ``UnaryOpNode`` header-only. + - | + Rework C++ unary ops in ``dwave::optimization::functional`` to provide additional + information about the functions. +upgrade: + - Rename C++ ``abs()`` function to ``absolute()``. diff --git a/tests/cpp/nodes/test_unaryop.cpp b/tests/cpp/nodes/test_unaryop.cpp index 8a960d279..79d6f6b6a 100644 --- a/tests/cpp/nodes/test_unaryop.cpp +++ b/tests/cpp/nodes/test_unaryop.cpp @@ -33,26 +33,26 @@ namespace dwave::optimization { TEMPLATE_TEST_CASE( "UnaryOpNode", "", - functional::abs, - functional::cos, - functional::exp, - functional::expit, - functional::logical, - functional::rint, - functional::sin, - functional::square, - functional::tanh, - std::negate, - std::logical_not + functional::absolute, + functional::cos, + functional::exp, + functional::expit, + functional::logical, + functional::logical_not, + functional::negate, + functional::rint, + functional::sin, + functional::square, + functional::tanh ) { auto graph = Graph(); auto func = TestType(); GIVEN("A constant scalar input") { - auto a_ptr = graph.emplace_node(-5); + auto* a_ptr = graph.emplace_node(-5); - auto p_ptr = graph.emplace_node>(a_ptr); + auto* p_ptr = graph.emplace_node>(a_ptr); THEN("The shape is also a scalar") { CHECK(p_ptr->ndim() == 0); @@ -77,9 +77,9 @@ TEMPLATE_TEST_CASE( } GIVEN("A dynamic array input") { - auto a_ptr = graph.emplace_node(5, 0, 5); + auto* a_ptr = graph.emplace_node(5, 0, 5); - auto p_ptr = graph.emplace_node>(a_ptr); + auto* p_ptr = graph.emplace_node>(a_ptr); graph.emplace_node(p_ptr); @@ -139,13 +139,15 @@ TEMPLATE_TEST_CASE( } GIVEN("A 0-d integer decision input") { - auto a_ptr = graph.emplace_node( + auto* a_ptr = graph.emplace_node( std::span{}, -100, 100 ); // Scalar output - auto p_ptr = graph.emplace_node>(a_ptr); + auto* p_ptr = graph.emplace_node>(a_ptr); + + graph.emplace_node(p_ptr); THEN("The integer is the operand") { CHECK(p_ptr->operands().size() == p_ptr->predecessors().size()); @@ -157,48 +159,42 @@ TEMPLATE_TEST_CASE( a_ptr->initialize_state(state, {-5}); graph.initialize_state(state); - THEN("The output has the value and shape we expect") { - CHECK(p_ptr->size(state) == 1); - CHECK(p_ptr->shape(state).size() == 0); - CHECK(p_ptr->view(state)[0] == func(-5)); - } + CHECK(p_ptr->size(state) == 1); + CHECK(p_ptr->shape(state).size() == 0); + CHECK_THAT(p_ptr->view(state), RangeEquals({func(-5)})); AND_WHEN("We change the integer's state and propagate") { a_ptr->set_value(state, 0, 17); - a_ptr->propagate(state); - p_ptr->propagate(state); + graph.propagate(state); + + CHECK_THAT(p_ptr->view(state), RangeEquals({func(17)})); - THEN("The output is what we expect") { CHECK(p_ptr->view(state)[0] == func(17)); } AND_WHEN("We commit") { - a_ptr->commit(state); - p_ptr->commit(state); + graph.commit(state); - THEN("The value hasn't changed") { CHECK(p_ptr->view(state)[0] == func(17)); } - THEN("The diffs are cleared") { CHECK(p_ptr->diff(state).size() == 0); } + CHECK_THAT(p_ptr->view(state), RangeEquals({func(17)})); + CHECK(p_ptr->diff(state).empty()); } AND_WHEN("We revert") { - a_ptr->revert(state); - p_ptr->revert(state); + graph.revert(state); - THEN("The value reverts to the previous") { - CHECK(p_ptr->view(state)[0] == func(-5)); - } - THEN("The diffs are cleared") { CHECK(p_ptr->diff(state).size() == 0); } + CHECK_THAT(p_ptr->view(state), RangeEquals({func(-5)})); + CHECK(p_ptr->diff(state).empty()); } } } } GIVEN("A 3-d integer decision input") { - auto a_ptr = graph.emplace_node( + auto* a_ptr = graph.emplace_node( std::span({2, 3, 2}), -100, 100 ); // Scalar output - auto p_ptr = graph.emplace_node>(a_ptr); + auto* p_ptr = graph.emplace_node>(a_ptr); THEN("The integers node is the operand") { CHECK(p_ptr->operands().size() == p_ptr->predecessors().size()); @@ -281,8 +277,8 @@ TEST_CASE("UnaryOpNode - AbsoluteNode") { auto graph = Graph(); GIVEN("An integer variable with domain [-3, 2]") { - auto i_ptr = graph.emplace_node(std::vector{}, -3, 2); - auto abs_ptr = graph.emplace_node(i_ptr); + auto* i_ptr = graph.emplace_node(std::vector{}, -3, 2); + auto* abs_ptr = graph.emplace_node(i_ptr); THEN("It has the min/max/integrality we expect") { CHECK(abs_ptr->min() == 0); @@ -292,8 +288,8 @@ TEST_CASE("UnaryOpNode - AbsoluteNode") { } GIVEN("An integer variable with domain [-2, 4]") { - auto i_ptr = graph.emplace_node(std::vector{}, -2, 4); - auto abs_ptr = graph.emplace_node(i_ptr); + auto* i_ptr = graph.emplace_node(std::vector{}, -2, 4); + auto* abs_ptr = graph.emplace_node(i_ptr); THEN("It has the min/max we expect") { CHECK(abs_ptr->min() == 0); @@ -308,14 +304,14 @@ TEST_CASE("UnaryOpNode - CosNode, SinNode, and TanhNode") { auto graph = Graph(); GIVEN("x with min/max of -100/+100, y = cos(x), z = sin(x), a = tanh(x)") { - auto x = graph.emplace_node(std::vector{-100, 1.5, +100}); - auto y = graph.emplace_node(x); - auto z = graph.emplace_node(x); - auto a = graph.emplace_node(x); - - CHECK(!y->integral()); - CHECK(!z->integral()); - CHECK(!a->integral()); + auto* x = graph.emplace_node(std::vector{-100, 1.5, +100}); + auto* y = graph.emplace_node(x); + auto* z = graph.emplace_node(x); + auto* a = graph.emplace_node(x); + + CHECK(not y->integral()); + CHECK(not z->integral()); + CHECK(not a->integral()); CHECK(y->min() == -1); CHECK(y->max() == +1); CHECK(z->min() == -1); @@ -329,8 +325,8 @@ TEST_CASE("UnaryOpNode - ExpitNode") { auto graph = Graph(); GIVEN("An arbitrary number") { double c = 3.0; - auto c_ptr = graph.emplace_node(c); - auto expit_ptr = graph.emplace_node(c_ptr); + auto* c_ptr = graph.emplace_node(c); + auto* expit_ptr = graph.emplace_node(c_ptr); auto state = graph.initialize_state(); CHECK(expit_ptr->min() == 1.0 / (1.0 + std::exp(-c))); CHECK(expit_ptr->max() == 1.0 / (1.0 + std::exp(-c))); @@ -338,16 +334,16 @@ TEST_CASE("UnaryOpNode - ExpitNode") { GIVEN("A negative number") { double c = -4.0; - auto c_ptr = graph.emplace_node(c); - auto expit_ptr = graph.emplace_node(c_ptr); + auto* c_ptr = graph.emplace_node(c); + auto* expit_ptr = graph.emplace_node(c_ptr); auto state = graph.initialize_state(); CHECK(expit_ptr->min() == 1.0 / (1.0 + std::exp(-c))); CHECK(expit_ptr->max() == 1.0 / (1.0 + std::exp(-c))); } GIVEN("A constant 1d array of doubles") { - auto c_ptr = graph.emplace_node(std::vector{-6.0, -0.3, 0.0, 1.2, 5.6}); - auto expit_ptr = graph.emplace_node(c_ptr); + auto* c_ptr = graph.emplace_node(std::vector{-6.0, -0.3, 0.0, 1.2, 5.6}); + auto* expit_ptr = graph.emplace_node(c_ptr); THEN("The min/max are expected") { THEN("expit(x) is not integral") { CHECK_FALSE(expit_ptr->integral()); } @@ -357,8 +353,8 @@ TEST_CASE("UnaryOpNode - ExpitNode") { } GIVEN("An integer with max domain") { - auto i_ptr = graph.emplace_node(std::vector{}); - auto expit_ptr = graph.emplace_node(i_ptr); + auto* i_ptr = graph.emplace_node(std::vector{}); + auto* expit_ptr = graph.emplace_node(i_ptr); graph.emplace_node(expit_ptr); THEN("The min/max are expected") { @@ -372,8 +368,8 @@ TEST_CASE("UnaryOpNode - ExpNode") { auto graph = Graph(); GIVEN("An arbitrary number") { double c = 3.0; - auto c_ptr = graph.emplace_node(c); - auto exp_ptr = graph.emplace_node(c_ptr); + auto* c_ptr = graph.emplace_node(c); + auto* exp_ptr = graph.emplace_node(c_ptr); auto state = graph.initialize_state(); CHECK(exp_ptr->min() == std::exp(c)); CHECK(exp_ptr->max() == std::exp(c)); @@ -381,16 +377,16 @@ TEST_CASE("UnaryOpNode - ExpNode") { GIVEN("A negative number") { double c = -4.0; - auto c_ptr = graph.emplace_node(c); - auto exp_ptr = graph.emplace_node(c_ptr); + auto* c_ptr = graph.emplace_node(c); + auto* exp_ptr = graph.emplace_node(c_ptr); auto state = graph.initialize_state(); CHECK(exp_ptr->min() == std::exp(c)); CHECK(exp_ptr->max() == std::exp(c)); } GIVEN("A constant 1d array of doubles") { - auto c_ptr = graph.emplace_node(std::vector{-6.0, -0.3, 0.0, 1.2, 5.6}); - auto exp_ptr = graph.emplace_node(c_ptr); + auto* c_ptr = graph.emplace_node(std::vector{-6.0, -0.3, 0.0, 1.2, 5.6}); + auto* exp_ptr = graph.emplace_node(c_ptr); THEN("The min/max are expected") { THEN("exp(x) is not integral") { CHECK_FALSE(exp_ptr->integral()); } @@ -400,8 +396,8 @@ TEST_CASE("UnaryOpNode - ExpNode") { } GIVEN("An integer with max domain") { - auto i_ptr = graph.emplace_node(std::vector{}); - auto exp_ptr = graph.emplace_node(i_ptr); + auto* i_ptr = graph.emplace_node(std::vector{}); + auto* exp_ptr = graph.emplace_node(i_ptr); graph.emplace_node(exp_ptr); THEN("The min is expected") { @@ -414,8 +410,8 @@ TEST_CASE("UnaryOpNode - ExpNode") { TEST_CASE("UnaryOpNode - LogNode") { auto graph = Graph(); GIVEN("A constant 1d array of doubles") { - auto c_ptr = graph.emplace_node(std::vector{6.0, 0.3, 1.0, 1.2, 5.6}); - auto log_ptr = graph.emplace_node(c_ptr); + auto* c_ptr = graph.emplace_node(std::vector{6.0, 0.3, 1.0, 1.2, 5.6}); + auto* log_ptr = graph.emplace_node(c_ptr); THEN("The min/max are expected") { THEN("log(x) is not integral") { CHECK_FALSE(log_ptr->integral()); } @@ -425,15 +421,15 @@ TEST_CASE("UnaryOpNode - LogNode") { } GIVEN("An arbitrary number") { double c = 10.0; - auto c_ptr = graph.emplace_node(c); - auto log_ptr = graph.emplace_node(c_ptr); + auto* c_ptr = graph.emplace_node(c); + auto* log_ptr = graph.emplace_node(c_ptr); auto state = graph.initialize_state(); CHECK(log_ptr->min() == std::log(c)); CHECK(log_ptr->max() == std::log(c)); } GIVEN("A negative number") { double c = -10.0; - auto c_ptr = graph.emplace_node(c); + auto* c_ptr = graph.emplace_node(c); CHECK_THROWS(graph.emplace_node(c_ptr)); } } @@ -441,8 +437,8 @@ TEST_CASE("UnaryOpNode - LogNode") { TEST_CASE("UnaryOpNode - LogicalNode") { auto graph = Graph(); GIVEN("A constant of mixed doubles and a negation of it") { - auto c_ptr = graph.emplace_node(std::vector{-2., -1., 0., 1., 2., -.5, .5}); - auto logical_ptr = graph.emplace_node(c_ptr); + auto* c_ptr = graph.emplace_node(std::vector{-2., -1., 0., 1., 2., -.5, .5}); + auto* logical_ptr = graph.emplace_node(c_ptr); THEN("NotNode is logical") { CHECK(logical_ptr->integral()); @@ -463,8 +459,8 @@ TEST_CASE("UnaryOpNode - LogicalNode") { TEST_CASE("UnaryOpNode - NegativeNode") { auto graph = Graph(); GIVEN("An integer array and an asymmetric domain") { - auto i_ptr = graph.emplace_node(std::vector{5}, -3, 8); - auto ni_ptr = graph.emplace_node(i_ptr); + auto* i_ptr = graph.emplace_node(std::vector{5}, -3, 8); + auto* ni_ptr = graph.emplace_node(i_ptr); THEN("Negative has the min/max we expect") { CHECK(i_ptr->min() == -3); @@ -480,8 +476,8 @@ TEST_CASE("UnaryOpNode - NegativeNode") { TEST_CASE("UnaryOpNode - NotNode") { auto graph = Graph(); GIVEN("A constant of mixed doubles and a negation of it") { - auto c_ptr = graph.emplace_node(std::vector{-2., -1., 0., 1., 2., -.5, .5}); - auto nc_ptr = graph.emplace_node(c_ptr); + auto* c_ptr = graph.emplace_node(std::vector{-2., -1., 0., 1., 2., -.5, .5}); + auto* nc_ptr = graph.emplace_node(c_ptr); THEN("NotNode is logical") { CHECK(nc_ptr->integral()); @@ -503,8 +499,8 @@ TEST_CASE("UnaryOpNode - RintNode") { auto graph = Graph(); GIVEN("An arbitrary number") { double c = 10.3; - auto c_ptr = graph.emplace_node(c); - auto rint_ptr = graph.emplace_node(c_ptr); + auto* c_ptr = graph.emplace_node(c); + auto* rint_ptr = graph.emplace_node(c_ptr); auto state = graph.initialize_state(); CHECK(rint_ptr->min() == std::rint(c)); CHECK(rint_ptr->max() == std::rint(c)); @@ -512,16 +508,16 @@ TEST_CASE("UnaryOpNode - RintNode") { GIVEN("A negative number") { double c = -10.5; - auto c_ptr = graph.emplace_node(c); - auto rint_ptr = graph.emplace_node(c_ptr); + auto* c_ptr = graph.emplace_node(c); + auto* rint_ptr = graph.emplace_node(c_ptr); auto state = graph.initialize_state(); CHECK(rint_ptr->min() == std::rint(c)); CHECK(rint_ptr->max() == std::rint(c)); } GIVEN("A constant 1d array of doubles") { - auto c_ptr = graph.emplace_node(std::vector{-3.8, 0.3, 1.2, 5.6}); - auto rint_ptr = graph.emplace_node(c_ptr); + auto* c_ptr = graph.emplace_node(std::vector{-3.8, 0.3, 1.2, 5.6}); + auto* rint_ptr = graph.emplace_node(c_ptr); THEN("The min/max are expected") { CHECK(rint_ptr->integral()); @@ -532,18 +528,18 @@ TEST_CASE("UnaryOpNode - RintNode") { WHEN("We access a constant 1d array using integers from a RintNode") { double c0 = 0.3; - auto c0_ptr = graph.emplace_node(c0); - auto rint0_ptr = graph.emplace_node(c0_ptr); + auto* c0_ptr = graph.emplace_node(c0); + auto* rint0_ptr = graph.emplace_node(c0_ptr); double c3 = 2.8; - auto c3_ptr = graph.emplace_node(c3); - auto rint3_ptr = graph.emplace_node(c3_ptr); + auto* c3_ptr = graph.emplace_node(c3); + auto* rint3_ptr = graph.emplace_node(c3_ptr); ; - auto arr_ptr = graph.emplace_node(std::vector{0, 10, 20, 30}); + auto* arr_ptr = graph.emplace_node(std::vector{0, 10, 20, 30}); - auto a0_ptr = graph.emplace_node(arr_ptr, rint0_ptr); - auto a3_ptr = graph.emplace_node(arr_ptr, rint3_ptr); + auto* a0_ptr = graph.emplace_node(arr_ptr, rint0_ptr); + auto* a3_ptr = graph.emplace_node(arr_ptr, rint3_ptr); auto state = graph.initialize_state(); @@ -557,8 +553,8 @@ TEST_CASE("UnaryOpNode - RintNode") { TEST_CASE("UnaryOpNode - SquareNode") { auto graph = Graph(); GIVEN("An integer with max domain") { - auto i_ptr = graph.emplace_node(std::vector{}); - auto square_ptr = graph.emplace_node(i_ptr); + auto* i_ptr = graph.emplace_node(std::vector{}); + auto* square_ptr = graph.emplace_node(i_ptr); THEN("The min/max are expected") { CHECK(square_ptr->min() == 0); @@ -574,8 +570,8 @@ TEST_CASE("UnaryOpNode - SquareNode") { TEST_CASE("UnaryOpNode - SquareRootNode") { auto graph = Graph(); GIVEN("An integer with max domain") { - auto i_ptr = graph.emplace_node(std::vector{}); - auto square_root_ptr = graph.emplace_node(i_ptr); + auto* i_ptr = graph.emplace_node(std::vector{}); + auto* square_root_ptr = graph.emplace_node(i_ptr); graph.emplace_node(square_root_ptr); THEN("The min/max are expected") { @@ -587,15 +583,15 @@ TEST_CASE("UnaryOpNode - SquareRootNode") { } GIVEN("An arbitrary number") { double c = 10.0; - auto c_ptr = graph.emplace_node(c); - auto square_root_ptr = graph.emplace_node(c_ptr); + auto* c_ptr = graph.emplace_node(c); + auto* square_root_ptr = graph.emplace_node(c_ptr); auto state = graph.initialize_state(); CHECK(square_root_ptr->min() == std::sqrt(c)); CHECK(square_root_ptr->max() == std::sqrt(c)); } GIVEN("A negative number") { double c = -10.0; - auto c_ptr = graph.emplace_node(c); + auto* c_ptr = graph.emplace_node(c); CHECK_THROWS(graph.emplace_node(c_ptr)); } } diff --git a/tests/cpp/test_functional.cpp b/tests/cpp/test_functional.cpp index 931de059c..f0d80f99d 100644 --- a/tests/cpp/test_functional.cpp +++ b/tests/cpp/test_functional.cpp @@ -12,24 +12,412 @@ // See the License for the specific language governing permissions and // limitations under the License. +#include +#include +#include + #include #include #include "dwave-optimization/functional.hpp" +#include "dwave-optimization/interval.hpp" +#include "dwave-optimization/typing.hpp" namespace dwave::optimization::functional { -TEMPLATE_TEST_CASE("modulus", "", double, int) { +TEMPLATE_LIST_TEST_CASE("absolute", "", DTypes) { + constexpr absolute op{}; + + SECTION("absolute(scalar)") { + CHECK(op(TestType(0)) == 0); + CHECK(op(TestType(1)) == 1); + + if constexpr (std::same_as) { + CHECK(op(true) == 1); // abs(bool) is identity + } else if constexpr (std::integral) { + CHECK(op(TestType(-1)) == 1); + CHECK(op(TestType(-10)) == 10); + CHECK(op(TestType(3)) == 3); + // We define abs(lowest) == max (see functional.hpp) + CHECK( + op(std::numeric_limits::lowest()) == std::numeric_limits::max() + ); + } else { // floating + CHECK(op(TestType(-1.5)) == 1.5); + CHECK(op(TestType(1.5)) == 1.5); + } + } + + SECTION("absolute(interval)") { + CHECK(not op(interval())); // op(empty) -> empty + + CHECK(op(interval(0, 0)) == interval(0, 0)); + CHECK(op(interval(0, 1)) == interval(0, 1)); + + if constexpr (not std::same_as) { + CHECK(op(interval(0, 5)) == interval(0, 5)); + CHECK(op(interval(-7, -4)) == interval(4, 7)); + CHECK(op(interval(-3, 1)) == interval(0, 3)); + CHECK(op(interval(-1, 3)) == interval(0, 3)); + CHECK(op(interval(-5, 5)) == interval(0, 5)); + } + if constexpr (std::floating_point) { + CHECK(op(interval(0.5, 5.2)) == interval(0.5, 5.2)); + CHECK(op(interval(-5.2, -0.5)) == interval(0.5, 5.2)); + } + } +} + +TEMPLATE_LIST_TEST_CASE("cos", "", DTypes) { + constexpr cos op{}; + + SECTION("cos(scalar)") { + CHECK(op(TestType(0)) == 1); // cos(0) == 1 exactly + if constexpr (not std::same_as) { + CHECK(op(TestType(1)) == std::cos(TestType(1))); + CHECK(op(TestType(3)) == std::cos(TestType(3))); + } + } + + SECTION("cos(interval)") { + CHECK(not op(interval())); // op(empty) -> empty + CHECK(op(interval(0, 0)) == interval(-1, +1)); + } +} + +TEMPLATE_LIST_TEST_CASE("exp", "", DTypes) { + constexpr exp op{}; + + SECTION("exp(scalar)") { + CHECK(op(TestType(0)) == 1); // exp(0) == 1 exactly + if constexpr (not std::same_as) { + CHECK(op(TestType(1)) == std::exp(TestType(1))); + CHECK(op(TestType(-2)) == std::exp(TestType(-2))); + } + } + + SECTION("exp(interval)") { + CHECK(not op(interval())); // op(empty) -> empty + CHECK(op(interval(0, 1)) == interval(op(TestType(0)), op(TestType(1)))); + if constexpr (not std::same_as) { + CHECK(op(interval(-2, 3)) == interval(op(TestType(-2)), op(TestType(3)))); + } + } + + SECTION("exp domain is unrestricted") { + CHECK(exp::domain == interval::all()); + } +} + +TEMPLATE_LIST_TEST_CASE("expit", "", DTypes) { + constexpr expit op{}; + + SECTION("expit(scalar)") { + CHECK(op(TestType(0)) == 0.5); // 1 / (1 + 1) + if constexpr (std::floating_point) { + // no NaN at the extremes + CHECK(op(TestType(-1000)) == 0); + CHECK(op(TestType(1000)) == 1); + } + } + + SECTION("expit(interval)") { + CHECK(not op(interval())); // op(empty) -> empty + CHECK(op(interval(0, 1)) == interval(op(TestType(0)), op(TestType(1)))); + if constexpr (not std::same_as) { + CHECK(op(interval(-2, 3)) == interval(op(TestType(-2)), op(TestType(3)))); + } + } +} + +TEMPLATE_LIST_TEST_CASE("log", "", DTypes) { + constexpr log op{}; + + SECTION("log(scalar)") { + CHECK(op(TestType(1)) == 0); // log(1) == 0 exactly + if constexpr (std::same_as) { + } else if constexpr (std::integral) { + CHECK(op(TestType(2)) == std::log(TestType(2))); + CHECK(op(TestType(10)) == std::log(TestType(10))); + } else { // floating + CHECK(op(TestType(2.5)) == std::log(TestType(2.5))); + CHECK(op(TestType(0.5)) == std::log(TestType(0.5))); + } + } + + SECTION("log(interval)") { + CHECK(not op(interval())); // op(empty) -> empty + if constexpr (std::same_as) { + CHECK(op(interval(1, 1)) == interval(op(TestType(1)), op(TestType(1)))); + } else { + CHECK(op(interval(1, 4)) == interval(op(TestType(1)), op(TestType(4)))); + CHECK(op(interval(2, 10)) == interval(op(TestType(2)), op(TestType(10)))); + } + } + + SECTION("log domain is non-negative") { + CHECK(log::domain == interval::nonnegative()); + } +} + +TEMPLATE_LIST_TEST_CASE("logical", "", DTypes) { + constexpr logical op{}; + + SECTION("logical()") { + CHECK(op(TestType(0)) == 0); + + if constexpr (std::same_as) { + CHECK(op(true) == 1); + } else if constexpr (std::integral) { + CHECK(op(TestType(-1)) == 1); + CHECK(op(TestType(1)) == 1); + CHECK(op(TestType(3)) == 1); + } else { // floating + CHECK(op(TestType(-.000001)) == 1); + CHECK(op(TestType(.000001)) == 1); + } + } + + SECTION("logical()") { + CHECK(not op(interval())); // op(null) -> null + + CHECK(op(interval(0, 0)) == interval(false, false)); + CHECK(op(interval(1, 1)) == interval(true, true)); + CHECK(op(interval(0, 1)) == interval(false, true)); + + if constexpr (std::same_as) { + // already covered + } else if constexpr (std::integral) { + CHECK(op(interval(0, 5)) == interval(false, true)); + CHECK(op(interval(1, 5)) == interval(true, true)); + + CHECK(op(interval(-3, 5)) == interval(false, true)); + + CHECK(op(interval(-3, 0)) == interval(false, true)); + CHECK(op(interval(-3, -1)) == interval(true, true)); + } else { // floating + CHECK(op(interval(0, .00001)) == interval(false, true)); + CHECK(op(interval(.000001, 5.5)) == interval(true, true)); + + CHECK(op(interval(-3.4, 13.2)) == interval(false, true)); + + CHECK(op(interval(-.00000001, 0)) == interval(false, true)); + CHECK(op(interval(-3.3, -.01)) == interval(true, true)); + } + } +} + +TEMPLATE_LIST_TEST_CASE("logical_not", "", DTypes) { + constexpr logical_not op{}; + + SECTION("logical_not()") { + CHECK(op(TestType(0)) == 1); + + if constexpr (std::same_as) { + CHECK(op(true) == 0); + } else if constexpr (std::integral) { + CHECK(op(TestType(-1)) == 0); + CHECK(op(TestType(1)) == 0); + CHECK(op(TestType(3)) == 0); + } else { // floating + CHECK(op(TestType(-.000001)) == 0); + CHECK(op(TestType(.000001)) == 0); + } + } + + SECTION("logical_not()") { + CHECK(not op(interval())); // op(null) -> null + + CHECK(op(interval(0, 0)) == interval(true, true)); + CHECK(op(interval(1, 1)) == interval(false, false)); + CHECK(op(interval(0, 1)) == interval(false, true)); + + if constexpr (std::same_as) { + // already covered + } else if constexpr (std::integral) { + CHECK(op(interval(0, 5)) == interval(false, true)); + CHECK(op(interval(1, 5)) == interval(false, false)); + + CHECK(op(interval(-3, 5)) == interval(false, true)); + + CHECK(op(interval(-3, 0)) == interval(false, true)); + CHECK(op(interval(-3, -1)) == interval(false, false)); + } else { // floating + CHECK(op(interval(0, .00001)) == interval(false, true)); + CHECK(op(interval(.000001, 5.5)) == interval(false, false)); + + CHECK(op(interval(-3.4, 13.2)) == interval(false, true)); + + CHECK(op(interval(-.00000001, 0)) == interval(false, true)); + CHECK(op(interval(-3.3, -.01)) == interval(false, false)); + } + } +} + +TEMPLATE_LIST_TEST_CASE("modulus", "", DTypes) { constexpr modulus op; - // test for consistency with NumPy CHECK(op(1, 0) == 0); CHECK(op(0, 1) == 0); - CHECK(op(-1, 0) == 0); - CHECK(op(0, -1) == 0); - CHECK(op(-1, -10) == -1); - CHECK(op(-1, 10) == 9); - CHECK(op(1, -10) == -9); + + if constexpr (not std::same_as) { + CHECK(op(-1, 0) == 0); + CHECK(op(0, -1) == 0); + + CHECK(op(-1, -10) == -1); + CHECK(op(-1, 10) == 9); + CHECK(op(1, -10) == -9); + } +} + +TEMPLATE_LIST_TEST_CASE("negate", "", DTypes) { + constexpr negate op{}; + if constexpr (not std::same_as) { + SECTION("negate(scalar)") { + CHECK(op(TestType(0)) == 0); + + if constexpr (std::integral) { + CHECK(op(TestType(3)) == -3); + CHECK(op(TestType(-3)) == 3); + } else { // floating + CHECK(op(TestType(1.5)) == -1.5); + CHECK(op(TestType(-1.5)) == 1.5); + } + } + + SECTION("negate(interval)") { + CHECK(not op(interval())); // op(empty) -> empty + + CHECK(op(interval(0, 1)) == interval(op(TestType(1)), op(TestType(0)))); + CHECK(op(interval(-2, 3)) == interval(op(TestType(3)), op(TestType(-2)))); + CHECK(op(interval(-5, -1)) == interval(op(TestType(-1)), op(TestType(-5)))); + } + } +} + +TEMPLATE_LIST_TEST_CASE("rint", "", DTypes) { + constexpr rint op{}; + + SECTION("rint(scalar)") { + CHECK(op(TestType(0)) == 0); + if constexpr (std::same_as) { + CHECK(op(true) == 1); + } else if constexpr (std::integral) { + CHECK(op(TestType(3)) == 3); + CHECK(op(TestType(-4)) == -4); + } else { // floating: rounds half to even + CHECK(op(TestType(2.5)) == 2); + CHECK(op(TestType(3.5)) == 4); + CHECK(op(TestType(-2.5)) == -2); + CHECK(op(TestType(2.4)) == std::rint(TestType(2.4))); + } + } + + SECTION("rint(interval)") { + CHECK(not op(interval())); // op(empty) -> empty + CHECK(op(interval(0, 1)) == interval(op(TestType(0)), op(TestType(1)))); + if constexpr (not std::same_as) { + CHECK(op(interval(-3, 4)) == interval(op(TestType(-3)), op(TestType(4)))); + } + } +} + +TEMPLATE_LIST_TEST_CASE("sin", "", DTypes) { + constexpr sin op{}; + + SECTION("sin(scalar)") { + CHECK(op(TestType(0)) == 0); // sin(0) == 0 exactly + if constexpr (not std::same_as) { + CHECK(op(TestType(1)) == std::sin(TestType(1))); + CHECK(op(TestType(2)) == std::sin(TestType(2))); + } + } + + SECTION("sin(interval)") { + CHECK(not op(interval())); // op(empty) -> empty + CHECK(op(interval(0, 0)) == interval(-1, +1)); + } +} + +TEMPLATE_LIST_TEST_CASE("square", "", DTypes) { + constexpr square op{}; + + SECTION("square(scalar)") { + CHECK(op(TestType(0)) == 0); + CHECK(op(TestType(1)) == 1); + if constexpr (std::same_as) { + // square(bool) is identity + } else if constexpr (std::integral) { + CHECK(op(TestType(3)) == 9); + CHECK(op(TestType(-3)) == 9); + CHECK(op(TestType(4)) == 16); + } else { // floating + CHECK(op(TestType(2.5)) == 6.25); + CHECK(op(TestType(-1.5)) == 2.25); + } + } + + SECTION("square(interval)") { + CHECK(not op(interval())); // op(empty) -> empty + + if constexpr (std::same_as) { + CHECK(op(interval(0, 1)) == interval(0, 1)); + } else { + CHECK(op(interval(2, 3)) == interval(op(TestType(2)), op(TestType(3)))); + CHECK(op(interval(-3, -2)) == interval(op(TestType(-2)), op(TestType(-3)))); + CHECK(op(interval(-3, 2)) == interval(TestType(0), op(TestType(-3)))); + CHECK(op(interval(-2, 3)) == interval(TestType(0), op(TestType(3)))); + } + } +} + +TEMPLATE_LIST_TEST_CASE("square_root", "", DTypes) { + constexpr square_root op{}; + + SECTION("square_root(scalar)") { + CHECK(op(TestType(0)) == 0); + CHECK(op(TestType(1)) == 1); + if constexpr (not std::same_as) { + CHECK(op(TestType(4)) == 2); + CHECK(op(TestType(9)) == 3); + if constexpr (std::floating_point) { + CHECK(op(TestType(2.0)) == std::sqrt(TestType(2.0))); + } + } + } + + SECTION("square_root(interval)") { + CHECK(not op(interval())); // op(empty) -> empty + CHECK(op(interval(0, 1)) == interval(op(TestType(0)), op(TestType(1)))); + if constexpr (not std::same_as) { + CHECK(op(interval(0, 4)) == interval(op(TestType(0)), op(TestType(4)))); + CHECK(op(interval(1, 9)) == interval(op(TestType(1)), op(TestType(9)))); + } + } + + SECTION("square_root domain is non-negative") { + CHECK(square_root::domain == interval::nonnegative()); + } +} + +TEMPLATE_LIST_TEST_CASE("tanh", "", DTypes) { + constexpr tanh op{}; + + SECTION("tanh(scalar)") { + CHECK(op(TestType(0)) == 0); // tanh(0) == 0 exactly + if constexpr (not std::same_as) { + CHECK(op(TestType(1)) == std::tanh(TestType(1))); + CHECK(op(TestType(-2)) == std::tanh(TestType(-2))); + } + } + + SECTION("tanh(interval)") { + CHECK(not op(interval())); // op(empty) -> empty + CHECK(op(interval(0, 1)) == interval(op(TestType(0)), op(TestType(1)))); + if constexpr (not std::same_as) { + CHECK(op(interval(-2, 3)) == interval(op(TestType(-2)), op(TestType(3)))); + } + } } } // namespace dwave::optimization::functional