diff --git a/dwave/optimization/include/dwave-optimization/nodes/numbers.hpp b/dwave/optimization/include/dwave-optimization/nodes/numbers.hpp index 59c93db1..67dcc66c 100644 --- a/dwave/optimization/include/dwave-optimization/nodes/numbers.hpp +++ b/dwave/optimization/include/dwave-optimization/nodes/numbers.hpp @@ -28,6 +28,55 @@ namespace dwave::optimization { /// A contiguous block of numbers. class NumberNode : public ArrayOutputMixin, public DecisionNode { public: + /// Stateless sum constraint information. + /// + /// A sum constraint constrains the sum of values within slices of the array. + /// The slices are defined along `axis` when `axis` has a value. If + /// `axis == std::nullopt`, the constraint is applied to the entire array, + /// which is treated as a flat array with a single slice. + /// + /// Constraints may be defined either: + /// - for ALL slices (the `operators` and `bounds` vectors have length 1), or + /// - PER slice (their lengths equal the number of slices along `axis`). + /// + /// Each slice sum is constrained by an `Operator` and a corresponding `bound`. + struct SumConstraint { + public: + /// Allowable operators. + enum class Operator { Equal, LessEqual, GreaterEqual }; + + /// To reduce the # of `IntegerNode` and `BinaryNode` constructors, we + /// allow only one constructor. + SumConstraint(std::optional axis, std::vector operators, + std::vector bounds); + + /// Return the axis along which slices are defined. + /// If `std::nullopt`, the sum constraint applies to the entire array. + std::optional axis() const { return axis_; }; + + /// Obtain the bound associated with a given slice. + double get_bound(const ssize_t slice) const; + + /// Obtain the operator associated with a given slice. + Operator get_operator(const ssize_t slice) const; + + /// The number of bounds. + ssize_t num_bounds() const { return bounds_.size(); }; + + /// The number of operators. + ssize_t num_operators() const { return operators_.size(); }; + + private: + /// Axis along which slices are defined (`std::nullopt` = whole array). + std::optional axis_ = std::nullopt; + /// Operator for ALL slices (vector has length one) or operators PER + /// slice (length of vector is equal to the number of slices). + std::vector operators_; + /// Bound for ALL slices (vector has length one) or bounds PER slice + /// (length of vector is equal to the number of slices). + std::vector bounds_; + }; + NumberNode() = delete; // Overloads needed by the Array ABC ************************************** @@ -68,6 +117,11 @@ class NumberNode : public ArrayOutputMixin, public DecisionNode { // Initialize the state of the node randomly template void initialize_state(State& state, Generator& rng) const { + // Currently do not support random node initialization with sum constraints. + if (sum_constraints_.size() > 0) { + throw std::invalid_argument("Cannot randomly initialize_state with sum constraints."); + } + std::vector values; const ssize_t size = this->size(); values.reserve(size); @@ -86,6 +140,9 @@ class NumberNode : public ArrayOutputMixin, public DecisionNode { return initialize_state(state, std::move(values)); } + /// @copydoc Node::propagate() + void propagate(State& state) const override; + // NumberNode methods ***************************************************** // In the given state, swap the value of index i with the value of index j. @@ -106,9 +163,19 @@ class NumberNode : public ArrayOutputMixin, public DecisionNode { // in a given index. void clip_and_set_value(State& state, ssize_t index, double value) const; + /// Return the stateless sum constraints. + const std::vector& sum_constraints() const; + + /// If the node is subject to sum constraints, we track the state + /// dependent sum of the values within each slice per constraint. The + /// returned vector is indexed in the same ordering as the constraints + /// given by `sum_constraints()`. + const std::vector>& sum_constraints_lhs(const State& state) const; + protected: explicit NumberNode(std::span shape, std::vector lower_bound, - std::vector upper_bound); + std::vector upper_bound, + std::vector sum_constraints = {}); // Return truth statement: 'value is valid in a given index'. virtual bool is_valid(ssize_t index, double value) const = 0; @@ -116,11 +183,23 @@ class NumberNode : public ArrayOutputMixin, public DecisionNode { // Default value in a given index. virtual double default_value(ssize_t index) const = 0; + /// Update the relevant sum constraints running sums (`lhs`) given that the + /// value stored at `index` is changed by `value_change` in a given state. + void update_sum_constraints_lhs(State& state, const ssize_t index, + const double value_change) const; + + /// Statelss global minimum and maximum of the values stored in NumberNode. double min_; double max_; + /// Stateless index-wise upper and lower bounds. std::vector lower_bounds_; std::vector upper_bounds_; + + /// Stateless sum constraints. + std::vector sum_constraints_; + /// Indicator variable that all sum constraint operators are "==". + bool sum_constraints_all_equals_; }; /// A contiguous block of integer numbers. @@ -134,33 +213,45 @@ class IntegerNode : public NumberNode { // Default to a single scalar integer with default bounds IntegerNode() : IntegerNode({}) {} - // Create an integer array with the user-defined bounds. - // Defaulting to the specified default bounds. + // Create an integer array with the user-defined index-wise bounds and sum + // constraints. Index-wise bounds default to the specified default bounds. + // By default, there are no sum constraints. IntegerNode(std::span shape, std::optional> lower_bound = std::nullopt, - std::optional> upper_bound = std::nullopt); + std::optional> upper_bound = std::nullopt, + std::vector sum_constraints = {}); IntegerNode(std::initializer_list shape, std::optional> lower_bound = std::nullopt, - std::optional> upper_bound = std::nullopt); + std::optional> upper_bound = std::nullopt, + std::vector sum_constraints = {}); IntegerNode(ssize_t size, std::optional> lower_bound = std::nullopt, - std::optional> upper_bound = std::nullopt); + std::optional> upper_bound = std::nullopt, + std::vector sum_constraints = {}); IntegerNode(std::span shape, double lower_bound, - std::optional> upper_bound = std::nullopt); + std::optional> upper_bound = std::nullopt, + std::vector sum_constraints = {}); IntegerNode(std::initializer_list shape, double lower_bound, - std::optional> upper_bound = std::nullopt); + std::optional> upper_bound = std::nullopt, + std::vector sum_constraints = {}); IntegerNode(ssize_t size, double lower_bound, - std::optional> upper_bound = std::nullopt); + std::optional> upper_bound = std::nullopt, + std::vector sum_constraints = {}); IntegerNode(std::span shape, std::optional> lower_bound, - double upper_bound); + double upper_bound, std::vector sum_constraints = {}); IntegerNode(std::initializer_list shape, - std::optional> lower_bound, double upper_bound); - IntegerNode(ssize_t size, std::optional> lower_bound, double upper_bound); - - IntegerNode(std::span shape, double lower_bound, double upper_bound); - IntegerNode(std::initializer_list shape, double lower_bound, double upper_bound); - IntegerNode(ssize_t size, double lower_bound, double upper_bound); + std::optional> lower_bound, double upper_bound, + std::vector sum_constraints = {}); + IntegerNode(ssize_t size, std::optional> lower_bound, double upper_bound, + std::vector sum_constraints = {}); + + IntegerNode(std::span shape, double lower_bound, double upper_bound, + std::vector sum_constraints = {}); + IntegerNode(std::initializer_list shape, double lower_bound, double upper_bound, + std::vector sum_constraints = {}); + IntegerNode(ssize_t size, double lower_bound, double upper_bound, + std::vector sum_constraints = {}); // Overloads needed by the Node ABC *************************************** @@ -190,33 +281,44 @@ class BinaryNode : public IntegerNode { /// A binary scalar variable with lower_bound = 0.0 and upper_bound = 1.0 BinaryNode() : BinaryNode({}) {} - // Create a binary array with the user-defined bounds. - // Defaulting to lower_bound = 0.0 and upper_bound = 1.0 + // Create a binary array with the user-defined index-wise bounds and sum + // constraints. Index-wise bounds default to lower_bound = 0.0 and + // upper_bound = 1.0. By default, there are no sum constraints. BinaryNode(std::span shape, std::optional> lower_bound = std::nullopt, - std::optional> upper_bound = std::nullopt); + std::optional> upper_bound = std::nullopt, + std::vector sum_constraints = {}); BinaryNode(std::initializer_list shape, std::optional> lower_bound = std::nullopt, - std::optional> upper_bound = std::nullopt); + std::optional> upper_bound = std::nullopt, + std::vector sum_constraints = {}); BinaryNode(ssize_t size, std::optional> lower_bound = std::nullopt, - std::optional> upper_bound = std::nullopt); + std::optional> upper_bound = std::nullopt, + std::vector sum_constraints = {}); BinaryNode(std::span shape, double lower_bound, - std::optional> upper_bound = std::nullopt); + std::optional> upper_bound = std::nullopt, + std::vector sum_constraints = {}); BinaryNode(std::initializer_list shape, double lower_bound, - std::optional> upper_bound = std::nullopt); + std::optional> upper_bound = std::nullopt, + std::vector sum_constraints = {}); BinaryNode(ssize_t size, double lower_bound, - std::optional> upper_bound = std::nullopt); + std::optional> upper_bound = std::nullopt, + std::vector sum_constraints = {}); BinaryNode(std::span shape, std::optional> lower_bound, - double upper_bound); + double upper_bound, std::vector sum_constraints = {}); BinaryNode(std::initializer_list shape, std::optional> lower_bound, - double upper_bound); - BinaryNode(ssize_t size, std::optional> lower_bound, double upper_bound); - - BinaryNode(std::span shape, double lower_bound, double upper_bound); - BinaryNode(std::initializer_list shape, double lower_bound, double upper_bound); - BinaryNode(ssize_t size, double lower_bound, double upper_bound); + double upper_bound, std::vector sum_constraints = {}); + BinaryNode(ssize_t size, std::optional> lower_bound, double upper_bound, + std::vector sum_constraints = {}); + + BinaryNode(std::span shape, double lower_bound, double upper_bound, + std::vector sum_constraints = {}); + BinaryNode(std::initializer_list shape, double lower_bound, double upper_bound, + std::vector sum_constraints = {}); + BinaryNode(ssize_t size, double lower_bound, double upper_bound, + std::vector sum_constraints = {}); // Flip the value (0 -> 1 or 1 -> 0) at index i in the given state. void flip(State& state, ssize_t i) const; diff --git a/dwave/optimization/libcpp/nodes/numbers.pxd b/dwave/optimization/libcpp/nodes/numbers.pxd index 0f08a25b..da37c98e 100644 --- a/dwave/optimization/libcpp/nodes/numbers.pxd +++ b/dwave/optimization/libcpp/nodes/numbers.pxd @@ -12,6 +12,7 @@ # See the License for the specific language governing permissions and # limitations under the License. +from libcpp.optional cimport optional from libcpp.vector cimport vector from dwave.optimization.libcpp.graph cimport ArrayNode @@ -19,16 +20,34 @@ from dwave.optimization.libcpp.state cimport State cdef extern from "dwave-optimization/nodes/numbers.hpp" namespace "dwave::optimization" nogil: - cdef cppclass IntegerNode(ArrayNode): - void initialize_state(State&, vector[double]) except+ - double lower_bound(Py_ssize_t index) - double upper_bound(Py_ssize_t index) - double lower_bound() except+ - double upper_bound() except+ - cdef cppclass BinaryNode(ArrayNode): + cdef cppclass NumberNode(ArrayNode): + struct SumConstraint: + # It appears Cython automatically assumes all (standard) enums are "public". + # Because of this, we use this very explict override. + enum class Operator "dwave::optimization::NumberNode::SumConstraint::Operator": + Equal + LessEqual + GreaterEqual + + SumConstraint(optional[Py_ssize_t] axis, vector[Operator] operators, + vector[double] bounds) + + optional[Py_ssize_t] axis() + double get_bound(Py_ssize_t slice) + Operator get_operator(Py_ssize_t slice) + Py_ssize_t num_bounds() + Py_ssize_t num_operators() + void initialize_state(State&, vector[double]) except+ double lower_bound(Py_ssize_t index) double upper_bound(Py_ssize_t index) double lower_bound() except+ double upper_bound() except+ + const vector[SumConstraint] sum_constraints() + + cdef cppclass IntegerNode(NumberNode): + pass + + cdef cppclass BinaryNode(IntegerNode): + pass diff --git a/dwave/optimization/model.py b/dwave/optimization/model.py index ea0f92f2..ae3fbbd3 100644 --- a/dwave/optimization/model.py +++ b/dwave/optimization/model.py @@ -165,7 +165,10 @@ def objective(self, value: ArraySymbol): def binary(self, shape: None | _ShapeLike = None, lower_bound: None | np.typing.ArrayLike = None, - upper_bound: None | np.typing.ArrayLike = None) -> BinaryVariable: + upper_bound: None | np.typing.ArrayLike = None, + subject_to: None | list[tuple[str, float]] = None, + axes_subject_to: None | list[tuple[int, str | list[str], float | list[float]]] = None + ) -> BinaryVariable: r"""Create a binary symbol as a decision variable. Args: @@ -178,6 +181,29 @@ def binary(self, shape: None | _ShapeLike = None, scalar (one bound for all variables) or an array (one bound for each variable). Non-boolean values are rounded down to the domain [0,1]. If None, the default value of 1 is used. + subject_to (optional): Constraint on the sum of the values in the + array. Must be an array of tuples where each tuple has the form: + (operator, bound). + - operator (str): The constraint operator ("<=", "==", or ">="). + - bound (float): The constraint bound. + If provided, the sum of values within the array must satisfy + the corresponding operator–bound pair. + Note 1: At most one sum constraint may be provided. + Note 2: If provided, axes_subject_to must None. + axes_subject_to (optional): Constraint on the sum of the values in + each slice along a fixed axis in the array. Must be an array of + tuples where each tuple has the form: (axis, operator(s), bound(s)). + - axis (int): The axis that the constraint is applied to. + - operator(s) (str | array[str]): The constraint operator(s) + ("<=", "==", or ">="). A single operator applies to all slice + along the axis; an array specifies one operator per slice. + - bound(s) (float | array[float]): The constraint bound. A + single value applies to all slices; an array specifies one + bound per slice. + If provided, the sum of values within each slice along the + specified axis must satisfy the corresponding operator–bound pair. + Note 1: At most one sum constraint may be provided. + Note 2: If provided, subject_to must None. Returns: A binary symbol. @@ -215,15 +241,43 @@ def binary(self, shape: None | _ShapeLike = None, >>> np.all([1, 0] == b.upper_bound()) True + This example adds a :math:`(2x3)`-sized binary symbol with + index-wise lower bounds and a sum constraint along axis 1. Let + x_i (int i : 0 <= i <= 2) denote the sum of the values within + slice i along axis 1. For each state defined for this symbol: + (x_0 <= 0), (x_1 == 2), and (x_2 >= 1). + + >>> from dwave.optimization.model import Model + >>> import numpy as np + >>> model = Model() + >>> b = model.binary([2, 3], lower_bound=[[0, 1, 1], [0, 1, 0]], + ... axes_subject_to=[(1, ["<=", "==", ">="], [0, 2, 1])]) + >>> np.all(b.sum_constraints() == [(1, ["<=", "==", ">="], [0, 2, 1])]) + True + + This example adds a :math:`6`-sized binary symbol such that + the sum of the values within the array is equal to 2. + + >>> from dwave.optimization.model import Model + >>> import numpy as np + >>> model = Model() + >>> b = model.binary(6, subject_to=[("==", 2)]) + >>> np.all(b.sum_constraints() == [(["=="], [2])]) + True + See Also: :class:`~dwave.optimization.symbols.numbers.BinaryVariable`: equivalent symbol. .. versionchanged:: 0.6.7 - Beginning in version 0.6.7, user-defined bounds and index-wise - bounds are supported. + Beginning in version 0.6.7, user-defined index-wise bounds are + supported. + + .. versionchanged:: 0.6.12 + Beginning in version 0.6.12, user-defined sum constraints are + supported. """ from dwave.optimization.symbols import BinaryVariable # avoid circular import - return BinaryVariable(self, shape, lower_bound, upper_bound) + return BinaryVariable(self, shape, lower_bound, upper_bound, subject_to, axes_subject_to) def constant(self, array_like: numpy.typing.ArrayLike) -> Constant: r"""Create a constant symbol. @@ -478,7 +532,9 @@ def integer( shape: None | _ShapeLike = None, lower_bound: None | numpy.typing.ArrayLike = None, upper_bound: None | numpy.typing.ArrayLike = None, - ) -> IntegerVariable: + subject_to: None | list[tuple[str, float]] = None, + axes_subject_to: None | list[tuple[int, str | list[str], float | list[float]]] = None + ) -> IntegerVariable: r"""Create an integer symbol as a decision variable. Args: @@ -491,7 +547,29 @@ def integer( scalar (one bound for all variables) or an array (one bound for each variable). Non-integer values are down up. If None, the default value is used. - + subject_to (optional): Constraint on the sum of the values in the + array. Must be an array of tuples where each tuple has the form: + (operator, bound). + - operator (str): The constraint operator ("<=", "==", or ">="). + - bound (float): The constraint bound. + If provided, the sum of values within the array must satisfy + the corresponding operator–bound pair. + Note 1: At most one sum constraint may be provided. + Note 2: If provided, axes_subject_to must None. + axes_subject_to (optional): Constraint on the sum of the values in + each slice along a fixed axis in the array. Must be an array of + tuples where each tuple has the form: (axis, operator(s), bound(s)). + - axis (int): The axis that the constraint is applied to. + - operator(s) (str | array[str]): The constraint operator(s) + ("<=", "==", or ">="). A single operator applies to all slice + along the axis; an array specifies one operator per slice. + - bound(s) (float | array[float]): The constraint bound. A + single value applies to all slices; an array specifies one + bound per slice. + If provided, the sum of values within each slice along the + specified axis must satisfy the corresponding operator–bound pair. + Note 1: At most one sum constraint may be provided. + Note 2: If provided, subject_to must None. Returns: An integer symbol. @@ -529,15 +607,44 @@ def integer( >>> np.all([1, 2] == i.upper_bound()) True + This example adds a :math:`(2x3)`-sized integer symbol with general + lower and upper bounds and a sum constraint along axis 1. Let x_i + (int i : 0 <= i <= 2) denote the sum of the values within + slice i along axis 1. For each state defined for this symbol: + (x_0 <= 2), (x_1 <= 4), and (x_2 <= 5). + + >>> from dwave.optimization.model import Model + >>> import numpy as np + >>> model = Model() + >>> i = model.integer([2, 3], lower_bound=1, upper_bound=3, + ... axes_subject_to=[(1, "<=", [2, 4, 5])]) + >>> np.all(i.sum_constraints() == [(1, ["<="], [2, 4, 5])]) + True + + This example adds a :math:`6`-sized integer symbol such that + the sum of the values within the array is less than or equal + to 20. + + >>> from dwave.optimization.model import Model + >>> import numpy as np + >>> model = Model() + >>> i = model.integer(6, subject_to=[("<=", 20)]) + >>> np.all(i.sum_constraints() == [(["<="], [20])]) + True + See Also: :class:`~dwave.optimization.symbols.numbers.IntegerVariable`: equivalent symbol. .. versionchanged:: 0.6.7 Beginning in version 0.6.7, user-defined index-wise bounds are supported. + + .. versionchanged:: 0.6.12 + Beginning in version 0.6.12, user-defined sum constraints are + supported. """ from dwave.optimization.symbols import IntegerVariable # avoid circular import - return IntegerVariable(self, shape, lower_bound, upper_bound) + return IntegerVariable(self, shape, lower_bound, upper_bound, subject_to, axes_subject_to) def list(self, n: int, diff --git a/dwave/optimization/src/nodes/numbers.cpp b/dwave/optimization/src/nodes/numbers.cpp index 5ad26c99..e89a8338 100644 --- a/dwave/optimization/src/nodes/numbers.cpp +++ b/dwave/optimization/src/nodes/numbers.cpp @@ -15,72 +15,436 @@ #include "dwave-optimization/nodes/numbers.hpp" #include +#include +#include #include +#include #include #include +#include #include "_state.hpp" +#include "dwave-optimization/array.hpp" +#include "dwave-optimization/common.hpp" namespace dwave::optimization { -// Base class to be used as interfaces. +NumberNode::SumConstraint::SumConstraint(std::optional axis, + std::vector operators, + std::vector bounds) + : axis_(axis), operators_(std::move(operators)), bounds_(std::move(bounds)) { + const size_t num_operators = operators_.size(); + const size_t num_bounds = bounds_.size(); + + if ((num_operators == 0) || (num_bounds == 0)) { + throw std::invalid_argument("`operators` and `bounds` must have non-zero size."); + } + + if (!axis_.has_value() && (num_operators != 1 || num_bounds != 1)) { + throw std::invalid_argument( + "If `axis` is undefined, `operators` and `bounds` must have size 1."); + } + + // If `operators` and `bounds` are both defined PER slice along `axis`, + // they must have the same size. + if ((num_operators > 1) && (num_bounds > 1) && (num_bounds != num_operators)) { + assert(axis.has_value()); + throw std::invalid_argument( + "`operators` and `bounds` should have same size if neither has size 1."); + } +} + +double NumberNode::SumConstraint::get_bound(const ssize_t slice) const { + assert(0 <= slice); + if (bounds_.size() == 1) return bounds_[0]; + assert(slice < static_cast(bounds_.size())); + return bounds_[slice]; +} + +NumberNode::SumConstraint::Operator NumberNode::SumConstraint::get_operator( + const ssize_t slice) const { + assert(0 <= slice); + if (operators_.size() == 1) return operators_[0]; + assert(slice < static_cast(operators_.size())); + return operators_[slice]; +} + +/// State dependant data attached to NumberNode +struct NumberNodeStateData : public ArrayNodeStateData { + // User does not provide sum constraints. + NumberNodeStateData(std::vector input) : ArrayNodeStateData(std::move(input)) {} + // User provides sum constraints. + NumberNodeStateData(std::vector input, + std::vector> sum_constraints_lhs) + : ArrayNodeStateData(std::move(input)), + sum_constraints_lhs(std::move(sum_constraints_lhs)), + prior_sum_constraints_lhs(this->sum_constraints_lhs) {} + + std::unique_ptr copy() const override { + return std::make_unique(*this); + } + + /// For each sum constraint, track the sum of the values within each slice. + /// `sum_constraints_lhs[i][j]` is the sum of the values within the `j`th slice + /// along the `axis`* defined by the `i`th sum constraint. + /// + /// (*) If `axis == std::nullopt`, the constraint is applied to the entire + /// array, which is treated as a flat array with a single slice. + std::vector> sum_constraints_lhs; + // Store a copy for NumberNode::revert() and commit() + std::vector> prior_sum_constraints_lhs; +}; double const* NumberNode::buff(const State& state) const noexcept { - return data_ptr(state)->buff(); + return data_ptr(state)->buff(); } std::span NumberNode::diff(const State& state) const noexcept { - return data_ptr(state)->diff(); + return data_ptr(state)->diff(); } double NumberNode::min() const { return min_; } double NumberNode::max() const { return max_; } +/// Given a NumberNode and an assignment of its variables (`number_data`), +/// compute and return a vector containing the sum of values for each slice +/// along the specified `axis`*. +/// +/// (*) If `axis == std::nullopt`, the constraint is applied to the entire +/// array, which is treated as a flat array with a single slice. +std::vector> get_sum_constraints_lhs(const NumberNode* node, + const std::vector& number_data) { + std::span node_shape = node->shape(); + const auto& sum_constraints = node->sum_constraints(); + const ssize_t num_sum_constraints = static_cast(sum_constraints.size()); + assert(num_sum_constraints <= static_cast(node_shape.size())); + assert(std::accumulate(node_shape.begin(), node_shape.end(), 1, std::multiplies()) == + static_cast(number_data.size())); + + // For each sum constraint, initialize the sum of the values contained in + // each of its slice to 0. + std::vector> sum_constraints_lhs; + sum_constraints_lhs.reserve(num_sum_constraints); + for (const NumberNode::SumConstraint& constraint : sum_constraints) { + const std::optional axis = constraint.axis(); + // Handle the case where the sum constraint applies to the entire array. + if (!axis.has_value()) { + // Array is treated as a flat array with a single axis. + sum_constraints_lhs.emplace_back(1, 0.0); + continue; + } + assert(axis.has_value()); + assert(0 <= *axis && *axis < static_cast(node_shape.size())); + // Emplace an all zeros vector of size equal to the number of slice + // along the given constrained axis. + sum_constraints_lhs.emplace_back(node_shape[*axis], 0.0); + } + + // Define a BufferIterator for `number_data` given the shape and strides of + // NumberNode and iterate over it. + for (BufferIterator it(number_data.data(), node_shape, node->strides()); + it != std::default_sentinel; ++it) { + // Increment the sum of the appropriate slice per sum constraint. + for (ssize_t i = 0; i < num_sum_constraints; ++i) { + const std::optional axis = sum_constraints[i].axis(); + // Handle the case where the sum constraint applies to the entire array. + if (!axis.has_value()) { + assert(sum_constraints_lhs[i].size() == 1); + sum_constraints_lhs[i].front() += *it; + continue; + } + assert(axis.has_value()); + assert(0 <= *axis && *axis < static_cast(it.location().size())); + const ssize_t slice = it.location()[*axis]; + assert(0 <= slice); + assert(slice < static_cast(sum_constraints_lhs[i].size())); + sum_constraints_lhs[i][slice] += *it; + } + } + + return sum_constraints_lhs; +} + +/// Determine whether the sum constraints are satisfied. +bool satisfies_sum_constraint(const std::vector& sum_constraints, + const std::vector>& sum_constraints_lhs) { + assert(sum_constraints.size() == sum_constraints_lhs.size()); + // Iterate over each sum constraint. + for (ssize_t i = 0, stop_i = static_cast(sum_constraints.size()); i < stop_i; ++i) { + const auto& constraint = sum_constraints[i]; + const auto& lhs = sum_constraints_lhs[i]; + + // Return `false` if any slice does not satisfy the constraint. + for (ssize_t slice = 0, stop_slice = static_cast(lhs.size()); slice < stop_slice; + ++slice) { + switch (constraint.get_operator(slice)) { + case NumberNode::SumConstraint::Operator::Equal: + if (lhs[slice] != constraint.get_bound(slice)) return false; + break; + case NumberNode::SumConstraint::Operator::LessEqual: + if (lhs[slice] > constraint.get_bound(slice)) return false; + break; + case NumberNode::SumConstraint::Operator::GreaterEqual: + if (lhs[slice] < constraint.get_bound(slice)) return false; + break; + default: + assert(false && "Unexpected operator type."); + unreachable(); + } + } + } + return true; +} + void NumberNode::initialize_state(State& state, std::vector&& number_data) const { if (number_data.size() != static_cast(this->size())) { throw std::invalid_argument("Size of data provided does not match node size"); } + for (ssize_t index = 0, stop = this->size(); index < stop; ++index) { if (!is_valid(index, number_data[index])) { throw std::invalid_argument("Invalid data provided for node"); } } - emplace_data_ptr(state, std::move(number_data)); + if (sum_constraints_.size() == 0) { // No sum constraints to consider. + emplace_data_ptr(state, std::move(number_data)); + } else { + // Given the assignment to NumberNode `number_data`, compute the sum + // of the values within each slice per sum constraint. + auto sum_constraints_lhs = get_sum_constraints_lhs(this, number_data); + + if (!satisfies_sum_constraint(sum_constraints_, sum_constraints_lhs)) { + throw std::invalid_argument("Initialized values do not satisfy sum constraint(s)."); + } + + emplace_data_ptr(state, std::move(number_data), + std::move(sum_constraints_lhs)); + } +} + +/// Given a `span` (used for strides or shape data), reorder the values +/// of the span such that the given `axis` is moved to the 0th index. +std::vector shift_axis_data(const std::span span, const ssize_t axis) { + const ssize_t ndim = span.size(); + std::vector output; + output.reserve(ndim); + output.emplace_back(span[axis]); + + for (ssize_t i = 0; i < ndim; ++i) { + if (i != axis) output.emplace_back(span[i]); + } + return output; +} + +/// Reverse the operation defined by `shift_axis_data()`. +std::vector undo_shift_axis_data(const std::span span, const ssize_t axis) { + const ssize_t ndim = span.size(); + std::vector output; + output.reserve(ndim); + + ssize_t i_span = 1; + for (ssize_t i = 0; i < ndim; ++i) { + if (i == axis) + output.emplace_back(span[0]); + else + output.emplace_back(span[i_span++]); + } + return output; +} + +/// Given a `lhs`, operator (`op`), and a `bound`, determine the non-negative amount +/// `delta` needed to be added to `lhs` to satisfy the constraint: (lhs+delta) op bound. +/// e.g. Given (lhs, op, bound) := (10, ==, 12), delta = 2 +/// e.g. Given (lhs, op, bound) := (10, <=, 12), delta = 0 +/// e.g. Given (lhs, op, bound) := (10, >=, 12), delta = 2 +/// Throws an error if `delta` is negative (corresponding with an infeasible sum constraint) +double sum_constraint_delta(const double lhs, const NumberNode::SumConstraint::Operator op, + const double bound) { + switch (op) { + case NumberNode::SumConstraint::Operator::Equal: + if (lhs > bound) throw std::invalid_argument("Infeasible sum constraint."); + // If error was not thrown, return amount needed to satisfy constraint. + return bound - lhs; + case NumberNode::SumConstraint::Operator::LessEqual: + if (lhs > bound) throw std::invalid_argument("Infeasible sum constraint."); + // If error was not thrown, sum satisfies constraint. + return 0.0; + case NumberNode::SumConstraint::Operator::GreaterEqual: + // If sum is less than bound, return the amount needed to equal it. + // Otherwise, sum satisfies constraint. + return (lhs < bound) ? (bound - lhs) : 0.0; + default: + assert(false && "Unexpected operator type."); + unreachable(); + } +} + +/// Given a NumberNode and exactly one sum constraint, assign values to +/// `values` (in-place) to satisfy the constraint. This method +/// A) Initially sets `values[i] = lower_bound(i)` for all i. +/// B) Incremements the values within each slice until they satisfy +/// the constraint (should this be possible). +void construct_state_given_exactly_one_sum_constraint(const NumberNode* node, + std::vector& values) { + const std::span node_shape = node->shape(); + const ssize_t ndim = node_shape.size(); + + // 1) Initialize all elements to their lower bounds. + for (ssize_t i = 0, stop = node->size(); i < stop; ++i) { + values.push_back(node->lower_bound(i)); + } + // 2) Determine the slice sums for the sum constraint. To improve performance, + // compute sum during previous loop. + assert(node->sum_constraints().size() == 1); + const std::vector lhs = get_sum_constraints_lhs(node, values).front(); + // Obtain the stateless sum constraint information. + const NumberNode::SumConstraint& constraint = node->sum_constraints().front(); + const std::optional axis = constraint.axis(); + + // Handle the case where the constraint applies to the entire array. + if (!axis.has_value()) { + assert(lhs.size() == 1); + // Determine the amount needed to adjust the values within the array. + double delta = sum_constraint_delta(lhs.front(), constraint.get_operator(0), + constraint.get_bound(0)); + if (delta == 0) return; // Bound is satisfied for entire array. + + for (ssize_t i = 0, stop = node->size(); i < stop; ++i) { + // Determine allowable amount we can increment the value in at `i`. + const double inc = std::min(delta, node->upper_bound(i) - values[i]); + + if (inc > 0) { // Apply the increment to both `values` and `delta`. + values[i] += inc; + delta -= inc; + if (delta == 0) break; // Bound is satisfied for entire array. + } + } + + if (delta != 0) throw std::invalid_argument("Infeasible sum constraint."); + return; + } + + assert(axis.has_value() && 0 <= *axis && *axis < ndim); + // We need a way to iterate over each slice along the constrainted axis and + // adjust its values until they satisfy the constraint. We do this by + // defining an iterator of `values` that traverses each slice one after + // another. This is equivalent to adjusting the node's shape and strides + // such that the data for the constrained axis is moved to position 0. + const std::vector buff_shape = shift_axis_data(node_shape, *axis); + const std::vector buff_strides = shift_axis_data(node->strides(), *axis); + // Define an iterator for `values` corresponding with the beginning of + // slice 0 along the constrained axis. + const BufferIterator slice_0_it(values.data(), ndim, buff_shape.data(), + buff_strides.data()); + // Determine the size of each slice along the constrained axis. + const ssize_t slice_size = std::accumulate(buff_shape.begin() + 1, buff_shape.end(), 1.0, + std::multiplies()); + + // 3) Iterate over each slice and adjust its values until they satisfy the + // sum constraint. + for (ssize_t slice = 0, stop = node_shape[*axis]; slice < stop; ++slice) { + // Determine the amount needed to adjust the values within the slice. + double delta = sum_constraint_delta(lhs[slice], constraint.get_operator(slice), + constraint.get_bound(slice)); + if (delta == 0) continue; // Sum constraint is satisfied for slice. + assert(delta >= 0); // Should only increment. + + // Determine how much we need to offset `slice_0_it` to get to the + // first index in the given `slice`. + const ssize_t offset = slice * slice_size; + // Iterate over all indices in the given slice. + for (auto slice_it = slice_0_it + offset, slice_end_it = slice_it + slice_size; + slice_it != slice_end_it; ++slice_it) { + assert(slice_it.location()[0] == slice); // We should be in the right slice. + // Determine the "true" index of `slice_it` given the node shape. + ssize_t index = + ravel_multi_index(undo_shift_axis_data(slice_it.location(), *axis), node_shape); + // Sanity check that we can correctly reverse the conversion. + assert(std::ranges::equal(shift_axis_data(unravel_index(index, node_shape), *axis), + slice_it.location())); + assert(0 <= index && index < static_cast(values.size())); + // Determine allowable amount we can increment the value in at `index`. + const double inc = std::min(delta, node->upper_bound(index) - *slice_it); + + if (inc > 0) { // Apply the increment to both `it` and `delta`. + *slice_it += inc; + delta -= inc; + if (delta == 0) break; // Sum constraint is satisfied for slice. + } + } + + if (delta != 0) throw std::invalid_argument("Infeasible sum constraint."); + } } void NumberNode::initialize_state(State& state) const { std::vector values; values.reserve(this->size()); - for (ssize_t i = 0, stop = this->size(); i < stop; ++i) { - values.push_back(default_value(i)); + + if (sum_constraints_.size() == 0) { + // No sum constraint to consider, initialize by default. + for (ssize_t i = 0, stop = this->size(); i < stop; ++i) { + values.push_back(default_value(i)); + } + initialize_state(state, std::move(values)); + } else if (sum_constraints_.size() == 1) { + construct_state_given_exactly_one_sum_constraint(this, values); + initialize_state(state, std::move(values)); + } else { + assert(false && "Multiple sum constraints not yet supported."); + unreachable(); + } +} + +void NumberNode::propagate(State& state) const { + // Should only propagate states that obey the sum constraint(s). + assert(satisfies_sum_constraint(sum_constraints_, sum_constraints_lhs(state))); + // Technically vestigial but will keep it for forms sake. + for (const auto& sv : successors()) { + sv->update(state, sv.index); } - initialize_state(state, std::move(values)); } void NumberNode::commit(State& state) const noexcept { - data_ptr(state)->commit(); + auto node_data = data_ptr(state); + // Manually store a copy of sum_constraints_lhs. + node_data->prior_sum_constraints_lhs = node_data->sum_constraints_lhs; + node_data->commit(); } void NumberNode::revert(State& state) const noexcept { - data_ptr(state)->revert(); + auto node_data = data_ptr(state); + // Manually reset sum_constraints_lhs. + node_data->sum_constraints_lhs = node_data->prior_sum_constraints_lhs; + node_data->revert(); } void NumberNode::exchange(State& state, ssize_t i, ssize_t j) const { - auto ptr = data_ptr(state); + auto ptr = data_ptr(state); // We expect the exchange to obey the index-wise bounds. assert(lower_bound(i) <= ptr->get(j)); assert(upper_bound(i) >= ptr->get(j)); assert(lower_bound(j) <= ptr->get(i)); assert(upper_bound(j) >= ptr->get(i)); - // Assert that i and j are valid indices occurs in ptr->exchange(). - // Exchange occurs IFF (i != j) and (buffer[i] != buffer[j]). - ptr->exchange(i, j); + // assert() that i and j are valid indices occurs in ptr->exchange(). + // State change occurs IFF (i != j) and (buffer[i] != buffer[j]). + if (ptr->exchange(i, j)) { + // If change occurred and sum constraint exist, update running sums. + // Nothing to update if all sum constraints are Equals. + if (!sum_constraints_all_equals_ && sum_constraints_.size() > 0) { + const double difference = ptr->get(i) - ptr->get(j); + // Index i changed from (what is now) ptr->get(j) to ptr->get(i) + update_sum_constraints_lhs(state, i, difference); + // Index j changed from (what is now) ptr->get(i) to ptr->get(j) + update_sum_constraints_lhs(state, j, -difference); + } + } } double NumberNode::get_value(State& state, ssize_t i) const { - return data_ptr(state)->get(i); + return data_ptr(state)->get(i); } double NumberNode::lower_bound(ssize_t index) const { @@ -118,10 +482,24 @@ double NumberNode::upper_bound() const { } void NumberNode::clip_and_set_value(State& state, ssize_t index, double value) const { + auto ptr = data_ptr(state); value = std::clamp(value, lower_bound(index), upper_bound(index)); - // Assert that i is a valid index occurs in data_ptr->set(). - // Set occurs IFF `value` != buffer[i] . - data_ptr(state)->set(index, value); + // assert() that i is a valid index occurs in ptr->set(). + // State change occurs IFF `value` != buffer[index]. + if (ptr->set(index, value)) { + // If change occurred and sum constraint exist, update running sums. + if (sum_constraints_.size() > 0) { + update_sum_constraints_lhs(state, index, value - diff(state).back().old); + } + } +} + +const std::vector& NumberNode::sum_constraints() const { + return sum_constraints_; +} + +const std::vector>& NumberNode::sum_constraints_lhs(const State& state) const { + return data_ptr(state)->sum_constraints_lhs; } template @@ -136,6 +514,18 @@ double get_extreme_index_wise_bound(const std::vector& bound) { return *it; } +bool all_sum_constraint_operators_are_equals( + std::vector& sum_constraints) { + for (const NumberNode::SumConstraint& constraint : sum_constraints) { + for (ssize_t i = 0, stop = constraint.num_operators(); i < stop; ++i) { + if (constraint.get_operator(i) != NumberNode::SumConstraint::Operator::Equal) + return false; + } + } + // Vacuously true if there are no sum constraints. + return true; +} + void check_index_wise_bounds(const NumberNode& node, const std::vector& lower_bounds_, const std::vector& upper_bounds_) { bool index_wise_bound = false; @@ -164,13 +554,80 @@ void check_index_wise_bounds(const NumberNode& node, const std::vector& } } +/// Check the user defined sum constraint(s). +void check_sum_constraints(const NumberNode* node) { + const std::vector& sum_constraints = node->sum_constraints(); + if (sum_constraints.size() == 0) return; // No sum constraints to check. + + const std::span shape = node->shape(); + // Used to assess if an axis is subject to multiple constraints. + std::vector constrained_axis(shape.size(), false); + // Used to assess if array is subject to multiple constraints. + bool constrained_array = false; + + for (const NumberNode::SumConstraint& constraint : sum_constraints) { + const std::optional axis = constraint.axis(); + const ssize_t num_operators = static_cast(constraint.num_operators()); + const ssize_t num_bounds = static_cast(constraint.num_bounds()); + + // Handle the case where the constraint applies to the entire array. + if (!axis.has_value()) { + // Checked in SumConstraint constructor + assert(num_operators == 1 && num_bounds == 1); + + if (constrained_array) + throw std::invalid_argument( + "Cannot define multiple sum constraints for the entire number array."); + constrained_array = true; + continue; + } + + assert(axis.has_value()); + + if (*axis < 0 || *axis >= static_cast(shape.size())) { + throw std::invalid_argument("Invalid constrained axis given number array shape."); + } + + if ((num_operators > 1) && (num_operators != shape[*axis])) { + throw std::invalid_argument("Invalid number of operators given number array shape."); + } + + if ((num_bounds > 1) && (num_bounds != shape[*axis])) { + throw std::invalid_argument("Invalid number of bounds given number array shape."); + } + + // Checked in SumConstraint constructor + assert(num_operators == num_bounds || num_operators == 1 || num_bounds == 1); + + if (constrained_axis[*axis]) { + throw std::invalid_argument( + "Cannot define multiple sum constraints for a single axis."); + } + constrained_axis[*axis] = true; + } + + // *Currently*, we only support one sum constraint. + if (sum_constraints.size() > 1) { + throw std::invalid_argument("Can define at most one sum constraint per number array."); + } + + // There are fasters ways to check whether the sum constraints are feasible. + // For now, fully attempt to construct a state and throw if impossible. + std::vector values; + values.reserve(node->size()); + construct_state_given_exactly_one_sum_constraint(node, values); +} + +// Base class to be used as interfaces. NumberNode::NumberNode(std::span shape, std::vector lower_bound, - std::vector upper_bound) + std::vector upper_bound, std::vector sum_constraints) : ArrayOutputMixin(shape), min_(get_extreme_index_wise_bound(lower_bound)), max_(get_extreme_index_wise_bound(upper_bound)), lower_bounds_(std::move(lower_bound)), - upper_bounds_(std::move(upper_bound)) { + upper_bounds_(std::move(upper_bound)), + sum_constraints_(std::move(sum_constraints)), + sum_constraints_all_equals_(all_sum_constraint_operators_are_equals(sum_constraints_)) { if ((shape.size() > 0) && (shape[0] < 0)) { throw std::invalid_argument("Number array cannot have dynamic size."); } @@ -180,18 +637,71 @@ NumberNode::NumberNode(std::span shape, std::vector lower } check_index_wise_bounds(*this, lower_bounds_, upper_bounds_); + check_sum_constraints(this); +} + +void NumberNode::update_sum_constraints_lhs(State& state, const ssize_t index, + const double value_change) const { + const auto& sum_constraints = this->sum_constraints(); + assert(value_change != 0); // Should not call when no change occurs. + assert(sum_constraints.size() != 0); // Should only call where applicable. + + // Get multidimensional indices for `index` so we can identify the slices + // `index` lies on per sum constraint. + const std::vector multi_index = unravel_index(index, this->shape()); + assert(sum_constraints.size() <= multi_index.size()); + // Get the slice sums for all sum constraints. + auto& sum_constraints_lhs = data_ptr(state)->sum_constraints_lhs; + assert(sum_constraints.size() == sum_constraints_lhs.size()); + + // For each sum constraint. + for (ssize_t i = 0, stop = static_cast(sum_constraints.size()); i < stop; ++i) { + const std::optional axis = sum_constraints[i].axis(); + + // Handle the case where the constraint applies to the entire array. + if (!axis.has_value()) { + assert(sum_constraints_lhs[i].size() == 1); + sum_constraints_lhs[i].front() += value_change; + continue; + } + + assert(axis.has_value() && 0 <= *axis && *axis < static_cast(multi_index.size())); + // Get the slice along the constrained axis the `value_change` occurs in. + const ssize_t slice = multi_index[*axis]; + assert(0 <= slice && slice < static_cast(sum_constraints_lhs[i].size())); + sum_constraints_lhs[i][slice] += value_change; // Offset slice sum. + } } // Integer Node *************************************************************** +/// Check the user defined sum constraint for IntegerNode. +void check_sum_constraint_integrality( + const std::vector& sum_constraints) { + if (sum_constraints.size() == 0) return; // No sum constraints to check. + + for (const NumberNode::SumConstraint& constraint : sum_constraints) { + for (ssize_t slice = 0, stop = constraint.num_bounds(); slice < stop; ++slice) { + const double bound = constraint.get_bound(slice); + if (bound != std::floor(bound)) { + throw std::invalid_argument( + "Sum constraint(s) for integral arrays must be integral."); + } + } + } +} + IntegerNode::IntegerNode(std::span shape, std::optional> lower_bound, - std::optional> upper_bound) - : NumberNode(shape, - lower_bound.has_value() ? std::move(*lower_bound) - : std::vector{default_lower_bound}, - upper_bound.has_value() ? std::move(*upper_bound) - : std::vector{default_upper_bound}) { + std::optional> upper_bound, + std::vector sum_constraints) + : NumberNode( + shape, + lower_bound.has_value() ? std::move(*lower_bound) + : std::vector{default_lower_bound}, + upper_bound.has_value() ? std::move(*upper_bound) + : std::vector{default_upper_bound}, + (check_sum_constraint_integrality(sum_constraints), std::move(sum_constraints))) { if (min_ < minimum_lower_bound || max_ > maximum_upper_bound) { throw std::invalid_argument("range provided for integers exceeds supported range"); } @@ -199,40 +709,59 @@ IntegerNode::IntegerNode(std::span shape, IntegerNode::IntegerNode(std::initializer_list shape, std::optional> lower_bound, - std::optional> upper_bound) - : IntegerNode(std::span(shape), std::move(lower_bound), std::move(upper_bound)) {} + std::optional> upper_bound, + std::vector sum_constraints) + : IntegerNode(std::span(shape), std::move(lower_bound), std::move(upper_bound), + std::move(sum_constraints)) {} IntegerNode::IntegerNode(ssize_t size, std::optional> lower_bound, - std::optional> upper_bound) - : IntegerNode({size}, std::move(lower_bound), std::move(upper_bound)) {} + std::optional> upper_bound, + std::vector sum_constraints) + : IntegerNode({size}, std::move(lower_bound), std::move(upper_bound), + std::move(sum_constraints)) {} IntegerNode::IntegerNode(std::span shape, double lower_bound, - std::optional> upper_bound) - : IntegerNode(shape, std::vector{lower_bound}, std::move(upper_bound)) {} + std::optional> upper_bound, + std::vector sum_constraints) + : IntegerNode(shape, std::vector{lower_bound}, std::move(upper_bound), + std::move(sum_constraints)) {} IntegerNode::IntegerNode(std::initializer_list shape, double lower_bound, - std::optional> upper_bound) - : IntegerNode(std::span(shape), std::vector{lower_bound}, std::move(upper_bound)) {} + std::optional> upper_bound, + std::vector sum_constraints) + : IntegerNode(std::span(shape), std::vector{lower_bound}, std::move(upper_bound), + std::move(sum_constraints)) {} IntegerNode::IntegerNode(ssize_t size, double lower_bound, - std::optional> upper_bound) - : IntegerNode({size}, std::vector{lower_bound}, std::move(upper_bound)) {} + std::optional> upper_bound, + std::vector sum_constraints) + : IntegerNode({size}, std::vector{lower_bound}, std::move(upper_bound), + std::move(sum_constraints)) {} IntegerNode::IntegerNode(std::span shape, - std::optional> lower_bound, double upper_bound) - : IntegerNode(shape, std::move(lower_bound), std::vector{upper_bound}) {} + std::optional> lower_bound, double upper_bound, + std::vector sum_constraints) + : IntegerNode(shape, std::move(lower_bound), std::vector{upper_bound}, + std::move(sum_constraints)) {} IntegerNode::IntegerNode(std::initializer_list shape, - std::optional> lower_bound, double upper_bound) - : IntegerNode(std::span(shape), std::move(lower_bound), std::vector{upper_bound}) {} + std::optional> lower_bound, double upper_bound, + std::vector sum_constraints) + : IntegerNode(std::span(shape), std::move(lower_bound), std::vector{upper_bound}, + std::move(sum_constraints)) {} IntegerNode::IntegerNode(ssize_t size, std::optional> lower_bound, - double upper_bound) - : IntegerNode({size}, std::move(lower_bound), std::vector{upper_bound}) {} - -IntegerNode::IntegerNode(std::span shape, double lower_bound, double upper_bound) - : IntegerNode(shape, std::vector{lower_bound}, std::vector{upper_bound}) {} + double upper_bound, std::vector sum_constraints) + : IntegerNode({size}, std::move(lower_bound), std::vector{upper_bound}, + std::move(sum_constraints)) {} + +IntegerNode::IntegerNode(std::span shape, double lower_bound, double upper_bound, + std::vector sum_constraints) + : IntegerNode(shape, std::vector{lower_bound}, std::vector{upper_bound}, + std::move(sum_constraints)) {} IntegerNode::IntegerNode(std::initializer_list shape, double lower_bound, - double upper_bound) + double upper_bound, std::vector sum_constraints) : IntegerNode(std::span(shape), std::vector{lower_bound}, - std::vector{upper_bound}) {} -IntegerNode::IntegerNode(ssize_t size, double lower_bound, double upper_bound) - : IntegerNode({size}, std::vector{lower_bound}, std::vector{upper_bound}) {} + std::vector{upper_bound}, std::move(sum_constraints)) {} +IntegerNode::IntegerNode(ssize_t size, double lower_bound, double upper_bound, + std::vector sum_constraints) + : IntegerNode({size}, std::vector{lower_bound}, std::vector{upper_bound}, + std::move(sum_constraints)) {} bool IntegerNode::integral() const { return true; } @@ -242,13 +771,19 @@ bool IntegerNode::is_valid(ssize_t index, double value) const { } void IntegerNode::set_value(State& state, ssize_t index, double value) const { + auto ptr = data_ptr(state); // We expect `value` to obey the index-wise bounds and to be an integer. assert(lower_bound(index) <= value); assert(upper_bound(index) >= value); assert(value == std::round(value)); - // Assert that i is a valid index occurs in data_ptr->set(). - // Set occurs IFF `value` != buffer[i] . - data_ptr(state)->set(index, value); + // assert() that i is a valid index occurs in ptr->set(). + // State change occurs IFF `value` != buffer[index]. + if (ptr->set(index, value)) { + // If change occurred and sum constraint exist, update running sums. + if (sum_constraints_.size() > 0) { + update_sum_constraints_lhs(state, index, value - diff(state).back().old); + } + } } double IntegerNode::default_value(ssize_t index) const { @@ -287,69 +822,111 @@ std::vector limit_bound_to_bool_domain(std::optional BinaryNode::BinaryNode(std::span shape, std::optional> lower_bound, - std::optional> upper_bound) + std::optional> upper_bound, + std::vector sum_constraints) : IntegerNode(shape, limit_bound_to_bool_domain(lower_bound), - limit_bound_to_bool_domain(upper_bound)) {} + limit_bound_to_bool_domain(upper_bound), std::move(sum_constraints)) {} BinaryNode::BinaryNode(std::initializer_list shape, std::optional> lower_bound, - std::optional> upper_bound) - : BinaryNode(std::span(shape), std::move(lower_bound), std::move(upper_bound)) {} + std::optional> upper_bound, + std::vector sum_constraints) + : BinaryNode(std::span(shape), std::move(lower_bound), std::move(upper_bound), + std::move(sum_constraints)) {} BinaryNode::BinaryNode(ssize_t size, std::optional> lower_bound, - std::optional> upper_bound) - : BinaryNode({size}, std::move(lower_bound), std::move(upper_bound)) {} + std::optional> upper_bound, + std::vector sum_constraints) + : BinaryNode({size}, std::move(lower_bound), std::move(upper_bound), + std::move(sum_constraints)) {} BinaryNode::BinaryNode(std::span shape, double lower_bound, - std::optional> upper_bound) - : BinaryNode(shape, std::vector{lower_bound}, std::move(upper_bound)) {} + std::optional> upper_bound, + std::vector sum_constraints) + : BinaryNode(shape, std::vector{lower_bound}, std::move(upper_bound), + std::move(sum_constraints)) {} BinaryNode::BinaryNode(std::initializer_list shape, double lower_bound, - std::optional> upper_bound) - : BinaryNode(std::span(shape), std::vector{lower_bound}, std::move(upper_bound)) {} + std::optional> upper_bound, + std::vector sum_constraints) + : BinaryNode(std::span(shape), std::vector{lower_bound}, std::move(upper_bound), + std::move(sum_constraints)) {} BinaryNode::BinaryNode(ssize_t size, double lower_bound, - std::optional> upper_bound) - : BinaryNode({size}, std::vector{lower_bound}, std::move(upper_bound)) {} + std::optional> upper_bound, + std::vector sum_constraints) + : BinaryNode({size}, std::vector{lower_bound}, std::move(upper_bound), + std::move(sum_constraints)) {} BinaryNode::BinaryNode(std::span shape, - std::optional> lower_bound, double upper_bound) - : BinaryNode(shape, std::move(lower_bound), std::vector{upper_bound}) {} + std::optional> lower_bound, double upper_bound, + std::vector sum_constraints) + : BinaryNode(shape, std::move(lower_bound), std::vector{upper_bound}, + std::move(sum_constraints)) {} BinaryNode::BinaryNode(std::initializer_list shape, - std::optional> lower_bound, double upper_bound) - : BinaryNode(std::span(shape), std::move(lower_bound), std::vector{upper_bound}) {} + std::optional> lower_bound, double upper_bound, + std::vector sum_constraints) + : BinaryNode(std::span(shape), std::move(lower_bound), std::vector{upper_bound}, + std::move(sum_constraints)) {} BinaryNode::BinaryNode(ssize_t size, std::optional> lower_bound, - double upper_bound) - : BinaryNode({size}, std::move(lower_bound), std::vector{upper_bound}) {} - -BinaryNode::BinaryNode(std::span shape, double lower_bound, double upper_bound) - : BinaryNode(shape, std::vector{lower_bound}, std::vector{upper_bound}) {} -BinaryNode::BinaryNode(std::initializer_list shape, double lower_bound, double upper_bound) + double upper_bound, std::vector sum_constraints) + : BinaryNode({size}, std::move(lower_bound), std::vector{upper_bound}, + std::move(sum_constraints)) {} + +BinaryNode::BinaryNode(std::span shape, double lower_bound, double upper_bound, + std::vector sum_constraints) + : BinaryNode(shape, std::vector{lower_bound}, std::vector{upper_bound}, + std::move(sum_constraints)) {} +BinaryNode::BinaryNode(std::initializer_list shape, double lower_bound, double upper_bound, + std::vector sum_constraints) : BinaryNode(std::span(shape), std::vector{lower_bound}, - std::vector{upper_bound}) {} -BinaryNode::BinaryNode(ssize_t size, double lower_bound, double upper_bound) - : BinaryNode({size}, std::vector{lower_bound}, std::vector{upper_bound}) {} + std::vector{upper_bound}, std::move(sum_constraints)) {} +BinaryNode::BinaryNode(ssize_t size, double lower_bound, double upper_bound, + std::vector sum_constraints) + : BinaryNode({size}, std::vector{lower_bound}, std::vector{upper_bound}, + std::move(sum_constraints)) {} void BinaryNode::flip(State& state, ssize_t i) const { - auto ptr = data_ptr(state); + auto ptr = data_ptr(state); // Variable should not be fixed. assert(lower_bound(i) != upper_bound(i)); - // Assert that i is a valid index occurs in ptr->set(). - // Set occurs IFF `value` != buffer[i] . - ptr->set(i, !ptr->get(i)); + // assert() that i is a valid index occurs in ptr->set(). + // State change occurs IFF `value` != buffer[i]. + if (ptr->set(i, !ptr->get(i))) { + // If change occurred and sum constraint exist, update running sums. + if (sum_constraints_.size() > 0) { + // If value changed from 0 -> 1, update by 1. + // If value changed from 1 -> 0, update by -1. + update_sum_constraints_lhs(state, i, (ptr->get(i) == 1) ? 1 : -1); + } + } } void BinaryNode::set(State& state, ssize_t i) const { + auto ptr = data_ptr(state); // We expect the set to obey the index-wise bounds. assert(upper_bound(i) == 1.0); - // Assert that i is a valid index occurs in data_ptr->set(). - // Set occurs IFF `value` != buffer[i] . - data_ptr(state)->set(i, 1.0); + // assert() that i is a valid index occurs in ptr->set(). + // State change occurs IFF `value` != buffer[i]. + if (ptr->set(i, 1.0)) { + // If change occurred and sum constraint exist, update running sums. + if (sum_constraints_.size() > 0) { + // If value changed from 0 -> 1, update by 1. + update_sum_constraints_lhs(state, i, 1.0); + } + } } void BinaryNode::unset(State& state, ssize_t i) const { + auto ptr = data_ptr(state); // We expect the set to obey the index-wise bounds. assert(lower_bound(i) == 0.0); - // Assert that i is a valid index occurs in data_ptr->set(). - // Set occurs IFF `value` != buffer[i] . - data_ptr(state)->set(i, 0.0); + // assert() that i is a valid index occurs in ptr->set(). + // State change occurs IFF `value` != buffer[i]. + if (ptr->set(i, 0.0)) { + // If change occurred and sum constraint exist, update running sums. + if (sum_constraints_.size() > 0) { + // If value changed from 1 -> 0, update by -1. + update_sum_constraints_lhs(state, i, -1.0); + } + } } } // namespace dwave::optimization diff --git a/dwave/optimization/symbols/numbers.pyx b/dwave/optimization/symbols/numbers.pyx index 0f98f530..577c4f8b 100644 --- a/dwave/optimization/symbols/numbers.pyx +++ b/dwave/optimization/symbols/numbers.pyx @@ -16,6 +16,7 @@ import json +import collections.abc import numpy as np from cython.operator cimport typeid @@ -27,25 +28,119 @@ from dwave.optimization._model cimport _Graph, _register, ArraySymbol, Symbol from dwave.optimization._utilities cimport as_cppshape from dwave.optimization.libcpp cimport dynamic_cast_ptr from dwave.optimization.libcpp.nodes.numbers cimport ( + NumberNode, BinaryNode, IntegerNode, ) from dwave.optimization.states cimport States +# Convert the str operators "==", "<=", ">=" into their corresponding +# C++ objects. +cdef NumberNode.SumConstraint.Operator _parse_python_operator(str op) except *: + if op == "==": + return NumberNode.SumConstraint.Operator.Equal + elif op == "<=": + return NumberNode.SumConstraint.Operator.LessEqual + elif op == ">=": + return NumberNode.SumConstraint.Operator.GreaterEqual + else: + raise TypeError(f"Invalid sum constraint operator: {op!r}") + + +# Convert the user-defined sum constraints for NumberNode into the +# corresponding C++ objects passed to NumberNode. +cdef vector[NumberNode.SumConstraint] _convert_python_sum_constraints( + subject_to: None | list[tuple[str, float]], + axes_subject_to: None | list[tuple[int, str | list[str], float | list[float]]]) except *: + cdef vector[NumberNode.SumConstraint] output + cdef optional[Py_ssize_t] cpp_axis = nullopt + cdef vector[NumberNode.SumConstraint.Operator] cpp_ops + cdef vector[double] cpp_bounds + cdef double[:] mem + + if subject_to is not None: + for constraint in subject_to: + if not isinstance(constraint, tuple) or len(constraint) != 2: + raise TypeError("A sum constraint on an entire number array must be" + " a tuple with two elements: `operator` and `bound`") + + py_ops, py_bounds = constraint + cpp_axis = nullopt + if not isinstance(py_ops, str): + raise TypeError("Sum constraint operator on entire number array should be a str.") + + cpp_ops.resize(1) + cpp_bounds.resize(1) + cpp_ops[0] = _parse_python_operator(py_ops) + cpp_bounds[0] = py_bounds + output.push_back(NumberNode.SumConstraint(cpp_axis, move(cpp_ops), move(cpp_bounds))) + + if axes_subject_to is not None: + for axis_constraint in axes_subject_to: + if not isinstance(axis_constraint, tuple) or len(axis_constraint) != 3: + raise TypeError("Each axis sum constraint must be a tuple with " + "three elements: axis, operator(s), bound(s)") + + axis, py_ops, py_bounds = axis_constraint + if not isinstance(axis, int): + raise TypeError("Constrained axis must be an int or None.") + cpp_axis = axis + + if isinstance(py_ops, str): + cpp_ops.resize(1) + # One operator defined for all slices. + cpp_ops[0] = _parse_python_operator(py_ops) + elif isinstance(py_ops, collections.abc.Iterable): + # Operator defined per slice. + cpp_ops.reserve(len(py_ops)) + for op in py_ops: + cpp_ops.push_back(_parse_python_operator(op)) + else: + raise TypeError("Axis sum constraint operator(s) should be str or an" + " iterable of str(s).") + + bound_array = np.asarray_chkfinite(py_bounds, dtype=np.double) + if (bound_array.ndim <= 1): + mem = bound_array.ravel() + cpp_bounds.reserve(mem.shape[0]) + for i in range(mem.shape[0]): + cpp_bounds.push_back(mem[i]) + else: + raise TypeError("Axis sum constraint bound(s) should be scalar or 1D-array.") + + output.push_back(NumberNode.SumConstraint(cpp_axis, move(cpp_ops), move(cpp_bounds))) + + return output + +# Convert the C++ operators into their corresponding str +cdef str _parse_cpp_operators(NumberNode.SumConstraint.Operator op): + if op == NumberNode.SumConstraint.Operator.Equal: + return "==" + elif op == NumberNode.SumConstraint.Operator.LessEqual: + return "<=" + elif op == NumberNode.SumConstraint.Operator.GreaterEqual: + return ">=" + else: + raise TypeError(f"Invalid sum constraint operator: {op!r}") + + cdef class BinaryVariable(ArraySymbol): """Binary decision-variable symbol. See also: :meth:`~dwave.optimization.model.Model.binary`: equivalent method. """ - def __init__(self, _Graph model, shape=None, lower_bound=None, upper_bound=None): + def __init__(self, _Graph model, shape=None, lower_bound=None, upper_bound=None, + subject_to: None | list[tuple[str, float]] = None, + axes_subject_to: None | list[tuple[int, str | list[str], float | list[float]]] = None): cdef vector[Py_ssize_t] cppshape = as_cppshape( tuple() if shape is None else shape ) cdef optional[vector[double]] cpplower_bound = nullopt cdef optional[vector[double]] cppupper_bound = nullopt + cdef vector[BinaryNode.SumConstraint] cpp_sum_constraints = _convert_python_sum_constraints(subject_to, axes_subject_to) cdef const double[:] mem if lower_bound is not None: @@ -75,7 +170,7 @@ cdef class BinaryVariable(ArraySymbol): raise ValueError("upper bound should be None, scalar, or the same shape") self.ptr = model._graph.emplace_node[BinaryNode]( - cppshape, cpplower_bound, cppupper_bound + cppshape, cpplower_bound, cppupper_bound, cpp_sum_constraints ) self.initialize_arraynode(model, self.ptr) @@ -116,10 +211,31 @@ cdef class BinaryVariable(ArraySymbol): with zf.open(info, "r") as f: upper_bound = np.load(f, allow_pickle=False) + # needs to be compatible with older versions + try: + info = zf.getinfo(directory + "sum_constraints.json") + except KeyError: + subject_to = None + axes_subject_to = None + else: + with zf.open(info, "r") as f: + subject_to = [] + axes_subject_to = [] + # Note that import is a list of lists, not a list of tuples. + # Hence we convert to tuple. We could also support lists. + for item in json.load(f): + if len(item) == 2: + # Inconvenient but `subject_to` expects scalars, not lists + subject_to.append((item[0][0], item[1][0])) + else: + axes_subject_to.append(tuple(item)) + return BinaryVariable(model, shape=shape_info["shape"], lower_bound=lower_bound, upper_bound=upper_bound, + subject_to=subject_to, + axes_subject_to=axes_subject_to ) def _into_zipfile(self, zf, directory): @@ -143,6 +259,32 @@ cdef class BinaryVariable(ArraySymbol): with zf.open(directory + "upper_bound.npy", mode="w", force_zip64=True) as f: np.save(f, upper_bound, allow_pickle=False) + sum_constraints = self.sum_constraints() + if len(sum_constraints) > 0: + # Using json here converts the tuples to lists + zf.writestr(directory + "sum_constraints.json", encoder.encode(sum_constraints)) + + def sum_constraints(self): + """Sum constraints of Binary symbol as a list of tuples where each tuple + is of the form: ([operator], [bound]) or (axis, [operator(s)], [bound(s)]).""" + cdef vector[NumberNode.SumConstraint] sum_constraints = self.ptr.sum_constraints() + cdef optional[Py_ssize_t] axis + + output = [] + for i in range(sum_constraints.size()): + constraint = &sum_constraints[i] + axis = constraint.axis() + py_ops = [_parse_cpp_operators(constraint.get_operator(j)) for j in + range(constraint.num_operators())] + py_bounds = [constraint.get_bound(j) for j in range(constraint.num_bounds())] + # axis may be nullopt + if axis.has_value(): + output.append((axis.value(), py_ops, py_bounds)) + else: + output.append((py_ops, py_bounds)) + + return output + def lower_bound(self): """Lower bound(s) of Binary symbol.""" try: @@ -212,13 +354,16 @@ cdef class IntegerVariable(ArraySymbol): See Also: :meth:`~dwave.optimization.model.Model.integer`: equivalent method. """ - def __init__(self, _Graph model, shape=None, lower_bound=None, upper_bound=None): + def __init__(self, _Graph model, shape=None, lower_bound=None, upper_bound=None, + subject_to: None | list[tuple[str, float]] = None, + axes_subject_to: None | list[tuple[int, str | list[str], float | list[float]]] = None): cdef vector[Py_ssize_t] cppshape = as_cppshape( tuple() if shape is None else shape ) cdef optional[vector[double]] cpplower_bound = nullopt cdef optional[vector[double]] cppupper_bound = nullopt + cdef vector[IntegerNode.SumConstraint] cpp_sum_constraints = _convert_python_sum_constraints(subject_to, axes_subject_to) cdef const double[:] mem if lower_bound is not None: @@ -248,7 +393,7 @@ cdef class IntegerVariable(ArraySymbol): raise ValueError("upper bound should be None, scalar, or the same shape") self.ptr = model._graph.emplace_node[IntegerNode]( - cppshape, cpplower_bound, cppupper_bound + cppshape, cpplower_bound, cppupper_bound, cpp_sum_constraints ) self.initialize_arraynode(model, self.ptr) @@ -289,10 +434,31 @@ cdef class IntegerVariable(ArraySymbol): with zf.open(info, "r") as f: upper_bound = np.load(f, allow_pickle=False) + # needs to be compatible with older versions + try: + info = zf.getinfo(directory + "sum_constraints.json") + except KeyError: + subject_to = None + axes_subject_to = None + else: + with zf.open(info, "r") as f: + subject_to = [] + axes_subject_to = [] + # Note that import is a list of lists, not a list of tuples. + # Hence we convert to tuple. We could also support lists. + for item in json.load(f): + if len(item) == 2: + # Inconvenient but `subject_to` expects scalars, not lists + subject_to.append((item[0][0], item[1][0])) + else: + axes_subject_to.append(tuple(item)) + return IntegerVariable(model, shape=shape_info["shape"], lower_bound=lower_bound, upper_bound=upper_bound, + subject_to=subject_to, + axes_subject_to=axes_subject_to ) def _into_zipfile(self, zf, directory): @@ -322,6 +488,32 @@ cdef class IntegerVariable(ArraySymbol): with zf.open(directory + "upper_bound.npy", mode="w", force_zip64=True) as f: np.save(f, upper_bound, allow_pickle=False) + sum_constraints = self.sum_constraints() + if len(sum_constraints) > 0: + # Using json here converts the tuples to lists + zf.writestr(directory + "sum_constraints.json", encoder.encode(sum_constraints)) + + def sum_constraints(self): + """Sum constraints of Integer symbol as a list of tuples where each tuple + is of the form: ([operator], [bound]) or (axis, [operator(s)], [bound(s)]).""" + cdef vector[NumberNode.SumConstraint] sum_constraints = self.ptr.sum_constraints() + cdef optional[Py_ssize_t] axis + + output = [] + for i in range(sum_constraints.size()): + constraint = &sum_constraints[i] + axis = constraint.axis() + py_ops = [_parse_cpp_operators(constraint.get_operator(j)) for j in + range(constraint.num_operators())] + py_bounds = [constraint.get_bound(j) for j in range(constraint.num_bounds())] + # axis may be nullopt + if axis.has_value(): + output.append((axis.value(), py_ops, py_bounds)) + else: + output.append((py_ops, py_bounds)) + + return output + def lower_bound(self): """Lower bound(s) of Integer symbol.""" try: diff --git a/releasenotes/notes/numbernode_axis_wise_bounds-594110e581c1115f.yaml b/releasenotes/notes/numbernode_axis_wise_bounds-594110e581c1115f.yaml new file mode 100644 index 00000000..d547c98f --- /dev/null +++ b/releasenotes/notes/numbernode_axis_wise_bounds-594110e581c1115f.yaml @@ -0,0 +1,5 @@ +--- +features: + - | + Axis-wise bounds added to NumberNode. Available to both IntegerNode and + BinaryNode. See #216` `_. diff --git a/tests/cpp/nodes/test_numbers.cpp b/tests/cpp/nodes/test_numbers.cpp index 0761c74e..fc683c58 100644 --- a/tests/cpp/nodes/test_numbers.cpp +++ b/tests/cpp/nodes/test_numbers.cpp @@ -18,6 +18,7 @@ #include "catch2/catch_test_macros.hpp" #include "catch2/matchers/catch_matchers.hpp" #include "catch2/matchers/catch_matchers_all.hpp" +#include "catch2/matchers/catch_matchers_range_equals.hpp" #include "dwave-optimization/graph.hpp" #include "dwave-optimization/nodes/numbers.hpp" @@ -25,6 +26,106 @@ using Catch::Matchers::RangeEquals; namespace dwave::optimization { +using SumConstraint = NumberNode::SumConstraint; +using Operator = NumberNode::SumConstraint::Operator; +using NumberNode::SumConstraint::Operator::Equal; +using NumberNode::SumConstraint::Operator::GreaterEqual; +using NumberNode::SumConstraint::Operator::LessEqual; + +TEST_CASE("SumConstraint") { + GIVEN("SumConstraint(axis = nullopt, operators = {}, bounds = {1.0})") { + REQUIRE_THROWS_WITH(SumConstraint(std::nullopt, {}, {1.0}), + "`operators` and `bounds` must have non-zero size."); + } + + GIVEN("SumConstraint(axis = nullopt, operators = {<=}, bounds = {})") { + REQUIRE_THROWS_WITH(SumConstraint(std::nullopt, {LessEqual}, {}), + "`operators` and `bounds` must have non-zero size."); + } + + GIVEN("SumConstraint(axis = nullopt, operators = {<=, ==}, bounds = {1.0})") { + REQUIRE_THROWS_WITH(SumConstraint(std::nullopt, {LessEqual, Equal}, {1.0}), + "If `axis` is undefined, `operators` and `bounds` must have size 1."); + } + + GIVEN("SumConstraint(axis = nullopt, operators = {<=}, bounds = {1.0, 2.0})") { + REQUIRE_THROWS_WITH(SumConstraint(std::nullopt, {LessEqual}, {1.0, 2.0}), + "If `axis` is undefined, `operators` and `bounds` must have size 1."); + } + + GIVEN("SumConstraint(axis = 0, operators = {}, bounds = {1.0})") { + REQUIRE_THROWS_WITH(SumConstraint(0, {}, {1.0}), + "`operators` and `bounds` must have non-zero size."); + } + + GIVEN("SumConstraint(axis = 0, operators = {<=}, bounds = {})") { + REQUIRE_THROWS_WITH(SumConstraint(0, {LessEqual}, {}), + "`operators` and `bounds` must have non-zero size."); + } + + GIVEN("SumConstraint(axis = 1, operators = {<=, ==, ==}, bounds = {2.0, 1.0})") { + REQUIRE_THROWS_WITH( + SumConstraint(1, {LessEqual, Equal, Equal}, {2.0, 1.0}), + "`operators` and `bounds` should have same size if neither has size 1."); + } + + GIVEN("SumConstraint(axis = nullopt, operators = {==}, bounds = {1.0})") { + SumConstraint sum_constraint(std::nullopt, {Equal}, {1.0}); + + THEN("The sum constraint info is correct") { + CHECK(sum_constraint.axis() == std::nullopt); + CHECK(sum_constraint.num_bounds() == 1); + CHECK(sum_constraint.get_bound(0) == 1.0); + CHECK(sum_constraint.num_operators() == 1); + CHECK(sum_constraint.get_operator(0) == Equal); + } + } + + GIVEN("SumConstraint(axis = 2, operators = {==, <=, >=}, bounds = {1.0})") { + SumConstraint sum_constraint(2, {Equal, LessEqual, GreaterEqual}, {1.0}); + + THEN("The sum constraint info is correct") { + CHECK(sum_constraint.axis() == 2); + CHECK(sum_constraint.num_bounds() == 1); + CHECK(sum_constraint.get_bound(0) == 1.0); + CHECK(sum_constraint.num_operators() == 3); + CHECK(sum_constraint.get_operator(0) == Equal); + CHECK(sum_constraint.get_operator(1) == LessEqual); + CHECK(sum_constraint.get_operator(2) == GreaterEqual); + } + } + + GIVEN("SumConstraint(axis = 2, operators = {==}, bounds = {1.0, 2.0, 3.0})") { + SumConstraint sum_constraint(2, {Equal}, {1.0, 2.0, 3.0}); + + THEN("The sum constraint info is correct") { + CHECK(sum_constraint.axis() == 2); + CHECK(sum_constraint.num_bounds() == 3); + CHECK(sum_constraint.get_bound(0) == 1.0); + CHECK(sum_constraint.get_bound(1) == 2.0); + CHECK(sum_constraint.get_bound(2) == 3.0); + CHECK(sum_constraint.num_operators() == 1); + CHECK(sum_constraint.get_operator(0) == Equal); + } + } + + GIVEN("SumConstraint(axis = 2, operators = {==, <=, >=}, bounds = {1.0, 2.0, 3.0})") { + SumConstraint sum_constraint(2, {Equal, LessEqual, GreaterEqual}, {1.0, 2.0, 3.0}); + + THEN("The sum constraint info is correct") { + CHECK(sum_constraint.axis() == 2); + CHECK(sum_constraint.num_bounds() == 3); + CHECK(sum_constraint.get_bound(0) == 1.0); + CHECK(sum_constraint.get_bound(1) == 2.0); + CHECK(sum_constraint.get_bound(2) == 3.0); + CHECK(sum_constraint.num_operators() == 3); + CHECK(sum_constraint.get_operator(0) == Equal); + CHECK(sum_constraint.get_operator(1) == LessEqual); + CHECK(sum_constraint.get_operator(2) == GreaterEqual); + } + } +} + TEST_CASE("BinaryNode") { auto graph = Graph(); @@ -439,6 +540,752 @@ TEST_CASE("BinaryNode") { REQUIRE_THROWS_WITH(graph.emplace_node(std::initializer_list{-1, 2}), "Number array cannot have dynamic size."); } + + // *********************** Sum Constraint tests ************************* + GIVEN("(2x3)-BinaryNode with a sum constraint on the invalid axis -1") { + std::vector sum_constraints{{-1, {Equal}, {1.0}}}; + + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{2, 3}, std::nullopt, + std::nullopt, sum_constraints), + "Invalid constrained axis given number array shape."); + } + + GIVEN("(2x3)-BinaryNode with a sum constraint on the invalid axis 2") { + std::vector sum_constraints{{2, {Equal}, {1.0}}}; + + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{2, 3}, std::nullopt, + std::nullopt, sum_constraints), + "Invalid constrained axis given number array shape."); + } + + GIVEN("(2x3)-BinaryNode with a sum constraint on axis: 1 with too many operators.") { + std::vector sum_constraints{{1, {LessEqual, Equal, Equal, Equal}, {1.0}}}; + + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{2, 3}, std::nullopt, + std::nullopt, sum_constraints), + "Invalid number of operators given number array shape."); + } + + GIVEN("(2x3)-BinaryNode with a sum constraint on axis: 1 with too few operators.") { + std::vector sum_constraints{{1, {LessEqual, Equal}, {1.0}}}; + + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{2, 3}, std::nullopt, + std::nullopt, sum_constraints), + "Invalid number of operators given number array shape."); + } + + GIVEN("(2x3)-BinaryNode with a sum constraint on axis: 1 with too many bounds.") { + std::vector sum_constraints{{1, {Equal}, {1.0, 2.0, 3.0, 4.0}}}; + + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{2, 3}, std::nullopt, + std::nullopt, sum_constraints), + "Invalid number of bounds given number array shape."); + } + + GIVEN("(2x3)-BinaryNode with a sum constraint on axis: 1 with too few bounds.") { + std::vector sum_constraints{{1, {LessEqual}, {1.0, 2.0}}}; + + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{2, 3}, std::nullopt, + std::nullopt, sum_constraints), + "Invalid number of bounds given number array shape."); + } + + GIVEN("(6)-BinaryNode with duplicate sum constraints over the entire array") { + SumConstraint sum_constraint{std::nullopt, {Equal}, {1.0}}; + std::vector sum_constraints{sum_constraint, sum_constraint}; + + REQUIRE_THROWS_WITH( + graph.emplace_node(6, std::nullopt, std::nullopt, sum_constraints), + "Cannot define multiple sum constraints for the entire number array."); + } + + GIVEN("(2x3)-BinaryNode with duplicate sum constraints on axis: 1") { + SumConstraint sum_constraint{1, {Equal}, {1.0}}; + std::vector sum_constraints{sum_constraint, sum_constraint}; + + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{2, 3}, std::nullopt, + std::nullopt, sum_constraints), + "Cannot define multiple sum constraints for a single axis."); + } + + GIVEN("(2x3)-BinaryNode with sum constraints on axis: 0 and the entire array.") { + SumConstraint sum_constraint{std::nullopt, {LessEqual}, {1.0}}; + SumConstraint sum_constraint_1{1, {LessEqual}, {1.0}}; + std::vector sum_constraints{sum_constraint, sum_constraint_1}; + + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{2, 3}, std::nullopt, + std::nullopt, sum_constraints), + "Can define at most one sum constraint per number array."); + } + + GIVEN("(2x3)-BinaryNode with sum constraints on axes: 0 and 1") { + SumConstraint sum_constraint_0{0, {LessEqual}, {1.0}}; + SumConstraint sum_constraint_1{1, {LessEqual}, {1.0}}; + std::vector sum_constraints{sum_constraint_0, sum_constraint_1}; + + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{2, 3}, std::nullopt, + std::nullopt, sum_constraints), + "Can define at most one sum constraint per number array."); + } + + GIVEN("(2x3x4)-BinaryNode with a non-integral sum constraint") { + std::vector sum_constraints{{1, {Equal}, {0.1}}}; + + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{2, 3, 4}, + std::nullopt, std::nullopt, sum_constraints), + "Sum constraint(s) for integral arrays must be integral."); + } + + GIVEN("(6)-BinaryNode with an infeasible sum constraint over the entire array.") { + auto graph = Graph(); + std::vector sum_constraints{{std::nullopt, {Equal}, {7.0}}}; + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{6}, std::nullopt, + std::nullopt, sum_constraints), + "Infeasible sum constraint."); + } + + GIVEN("(3x2)-BinaryNode with an infeasible sum constraint over the entire array.") { + auto graph = Graph(); + std::vector sum_constraints{{std::nullopt, {GreaterEqual}, {7.0}}}; + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{3, 2}, std::nullopt, + std::nullopt, sum_constraints), + "Infeasible sum constraint."); + } + + GIVEN("(2x2x2)-BinaryNode with an infeasible sum constraint over the entire array.") { + auto graph = Graph(); + std::vector sum_constraints{{std::nullopt, {LessEqual}, {-1.0}}}; + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{2, 2, 2}, + std::nullopt, std::nullopt, sum_constraints), + "Infeasible sum constraint."); + } + + GIVEN("(3x2x2)-BinaryNode with an infeasible sum constraint on axis: 0") { + auto graph = Graph(); + std::vector sum_constraints{ + {0, {Equal, LessEqual, GreaterEqual}, {5.0, 2.0, 3.0}}}; + // Each slice along axis 0 has size 4. There is no feasible assignment + // to the values in slice 0 (along axis 0) that results in a sum equal + // to 5. + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{3, 2, 2}, + std::nullopt, std::nullopt, sum_constraints), + "Infeasible sum constraint."); + } + + GIVEN("(3x2x2)-BinaryNode with an infeasible sum constraint on axis: 1") { + auto graph = Graph(); + std::vector sum_constraints{{1, {Equal, GreaterEqual}, {5.0, 7.0}}}; + // Each slice along axis 1 has size 6. There is no feasible assignment + // to the values in slice 1 (along axis 1) that results in a sum + // greater than or equal to 7. + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{3, 2, 2}, + std::nullopt, std::nullopt, sum_constraints), + "Infeasible sum constraint."); + } + + GIVEN("(3x2x2)-BinaryNode with an infeasible sum constraint on axis: 2") { + auto graph = Graph(); + std::vector sum_constraints{{2, {Equal, LessEqual}, {5.0, -1.0}}}; + // Each slice along axis 2 has size 6. There is no feasible assignment + // to the values in slice 1 (along axis 2) that results in a sum less + // than or equal to -1. + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{3, 2, 2}, + std::nullopt, std::nullopt, sum_constraints), + "Infeasible sum constraint."); + } + + GIVEN("(6)-BinaryNode with a feasible sum constraint over the entire array.") { + auto graph = Graph(); + std::vector lower_bounds{0, 0, 1, 0, 0, 1}; + std::vector upper_bounds{0, 1, 1, 1, 1, 1}; + std::vector sum_constraints{{std::nullopt, {Equal}, {3.0}}}; + auto bnode_ptr = + graph.emplace_node(6, lower_bounds, upper_bounds, sum_constraints); + + THEN("Sum constraint is correct") { + CHECK(bnode_ptr->sum_constraints().size() == 1); + SumConstraint bnode_sum_constraint = bnode_ptr->sum_constraints()[0]; + CHECK(bnode_sum_constraint.axis() == std::nullopt); + CHECK(bnode_sum_constraint.num_bounds() == 1); + CHECK(bnode_sum_constraint.get_bound(0) == 3.0); + CHECK(bnode_sum_constraint.num_operators() == 1); + CHECK(bnode_sum_constraint.get_operator(0) == Equal); + } + + WHEN("We create a state by initialize_state()") { + auto state = graph.initialize_state(); + graph.initialize_state(state); + std::vector expected_init{0, 1, 1, 0, 0, 1}; + auto sum_constraints_lhs = bnode_ptr->sum_constraints_lhs(state); + + THEN("Sum constraint sums and state are correct") { + CHECK(bnode_ptr->sum_constraints_lhs(state).size() == 1); + CHECK(bnode_ptr->sum_constraints_lhs(state).data()[0].size() == 1); + CHECK_THAT(bnode_ptr->sum_constraints_lhs(state)[0], RangeEquals({3})); + CHECK_THAT(bnode_ptr->view(state), RangeEquals(expected_init)); + } + } + } + + GIVEN("(3x2x2)-BinaryNode with a feasible sum constraint on axis: 0") { + auto graph = Graph(); + std::vector lower_bounds{0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}; + std::vector upper_bounds{0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1}; + std::vector sum_constraints{ + {0, {Equal, LessEqual, GreaterEqual}, {1.0, 2.0, 3.0}}}; + auto bnode_ptr = + graph.emplace_node(std::initializer_list{3, 2, 2}, + lower_bounds, upper_bounds, sum_constraints); + + THEN("Sum constraint is correct") { + CHECK(bnode_ptr->sum_constraints().size() == 1); + SumConstraint bnode_sum_constraint = bnode_ptr->sum_constraints()[0]; + CHECK(bnode_sum_constraint.axis() == 0); + CHECK(bnode_sum_constraint.num_bounds() == 3); + CHECK(bnode_sum_constraint.get_bound(0) == 1.0); + CHECK(bnode_sum_constraint.get_bound(1) == 2.0); + CHECK(bnode_sum_constraint.get_bound(2) == 3.0); + CHECK(bnode_sum_constraint.num_operators() == 3); + CHECK(bnode_sum_constraint.get_operator(0) == Equal); + CHECK(bnode_sum_constraint.get_operator(1) == LessEqual); + CHECK(bnode_sum_constraint.get_operator(2) == GreaterEqual); + } + + WHEN("We create a state by initialize_state()") { + auto state = graph.initialize_state(); + graph.initialize_state(state); + // import numpy as np + // a = np.asarray([i for i in range(3*2*2)]).reshape(3, 2, 2) + // print(a[0, :, :].flatten()) + // >>> [0 1 2 3] + // print(a[1, :, :].flatten()) + // >>> [4 5 6 7] + // print(a[2, :, :].flatten()) + // >>> [ 8 9 10 11] + // + // Cannonically least state that satisfies the index-wise bounds + // and sum constraints. + // slice 0 slice 1 slice 2 + // 0, 0 0, 0 1, 1 + // 1, 0 0, 0 0, 1 + std::vector expected_init{0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1}; + auto sum_constraints_lhs = bnode_ptr->sum_constraints_lhs(state); + + THEN("Sum constraint sums and state are correct") { + CHECK(bnode_ptr->sum_constraints_lhs(state).size() == 1); + CHECK(bnode_ptr->sum_constraints_lhs(state).data()[0].size() == 3); + CHECK_THAT(bnode_ptr->sum_constraints_lhs(state)[0], RangeEquals({1, 0, 3})); + CHECK_THAT(bnode_ptr->view(state), RangeEquals(expected_init)); + } + } + } + + GIVEN("(3x2x2)-BinaryNode with a feasible sum constraint on axis: 1") { + auto graph = Graph(); + std::vector lower_bounds{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}; + std::vector upper_bounds{0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1}; + std::vector sum_constraints{{1, {LessEqual, GreaterEqual}, {1.0, 5.0}}}; + auto bnode_ptr = + graph.emplace_node(std::initializer_list{3, 2, 2}, + lower_bounds, upper_bounds, sum_constraints); + + THEN("Sum constraint is correct") { + CHECK(bnode_ptr->sum_constraints().size() == 1); + SumConstraint bnode_sum_constraint = bnode_ptr->sum_constraints()[0]; + CHECK(bnode_sum_constraint.axis() == 1); + CHECK(bnode_sum_constraint.num_bounds() == 2); + CHECK(bnode_sum_constraint.get_bound(0) == 1.0); + CHECK(bnode_sum_constraint.get_bound(1) == 5.0); + CHECK(bnode_sum_constraint.num_operators() == 2); + CHECK(bnode_sum_constraint.get_operator(0) == LessEqual); + CHECK(bnode_sum_constraint.get_operator(1) == GreaterEqual); + } + + WHEN("We create a state by initialize_state()") { + auto state = graph.initialize_state(); + graph.initialize_state(state); + // import numpy as np + // a = np.asarray([i for i in range(3*2*2)]).reshape(3, 2, 2) + // print(a[:, 0, :].flatten()) + // >>> [0 1 4 5 8 9] + // print(a[:, 1, :].flatten()) + // >>> [ 2 3 6 7 10 11] + // + // Cannonically least state that satisfies sum constraint + // slice 0 slice 1 + // 0, 0 1, 1 + // 0, 0 1, 1 + // 0, 0 0, 1 + std::vector expected_init{0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1}; + + auto sum_constraints_lhs = bnode_ptr->sum_constraints_lhs(state); + + THEN("Sum constraint sums and state are correct") { + CHECK(bnode_ptr->sum_constraints_lhs(state).size() == 1); + CHECK(bnode_ptr->sum_constraints_lhs(state).data()[0].size() == 2); + CHECK_THAT(bnode_ptr->sum_constraints_lhs(state)[0], RangeEquals({0, 5})); + CHECK_THAT(bnode_ptr->view(state), RangeEquals(expected_init)); + } + } + } + + GIVEN("(3x2x2)-BinaryNode with a feasible sum constraint on axis: 2") { + auto graph = Graph(); + std::vector lower_bounds{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0}; + std::vector upper_bounds{0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; + std::vector sum_constraints{{2, {Equal, GreaterEqual}, {3.0, 6.0}}}; + auto bnode_ptr = + graph.emplace_node(std::initializer_list{3, 2, 2}, + lower_bounds, upper_bounds, sum_constraints); + + THEN("Sum constraint is correct") { + CHECK(bnode_ptr->sum_constraints().size() == 1); + SumConstraint bnode_sum_constraint = bnode_ptr->sum_constraints()[0]; + CHECK(bnode_sum_constraint.axis() == 2); + CHECK(bnode_sum_constraint.num_bounds() == 2); + CHECK(bnode_sum_constraint.get_bound(0) == 3.0); + CHECK(bnode_sum_constraint.get_bound(1) == 6.0); + CHECK(bnode_sum_constraint.num_operators() == 2); + CHECK(bnode_sum_constraint.get_operator(0) == Equal); + CHECK(bnode_sum_constraint.get_operator(1) == GreaterEqual); + } + + WHEN("We create a state by initialize_state()") { + auto state = graph.initialize_state(); + graph.initialize_state(state); + // import numpy as np + // a = np.asarray([i for i in range(3*2*2)]).reshape(3, 2, 2) + // print(a[:, :, 0].flatten()) + // >>> [ 0 2 4 6 8 10] + // print(a[:, :, 1].flatten()) + // >>> [ 1 3 5 7 9 11] + // + // Cannonically least state that satisfies the index-wise bounds + // and sum constraint. + // slice 0 slice 1 + // 0, 1 1, 1 + // 1, 0 1, 1 + // 0, 1 1, 1 + std::vector expected_init{0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1}; + auto sum_constraints_lhs = bnode_ptr->sum_constraints_lhs(state); + + THEN("Sum constraint sums and state are correct") { + CHECK(bnode_ptr->sum_constraints_lhs(state).size() == 1); + CHECK(bnode_ptr->sum_constraints_lhs(state).data()[0].size() == 2); + CHECK_THAT(bnode_ptr->sum_constraints_lhs(state)[0], RangeEquals({3, 6})); + CHECK_THAT(bnode_ptr->view(state), RangeEquals(expected_init)); + } + } + } + + GIVEN("(2)-BinaryNode with a sum constraint over the entire array") { + auto graph = Graph(); + std::vector sum_constraints{{std::nullopt, {Equal}, {1}}}; + auto bnode_ptr = + graph.emplace_node(2, std::nullopt, std::nullopt, sum_constraints); + + WHEN("We initialize an invalid states") { + auto state = graph.empty_state(); + std::vector init_values{0, 0}; + CHECK_THROWS_WITH(bnode_ptr->initialize_state(state, init_values), + "Initialized values do not satisfy sum constraint(s)."); + } + + WHEN("We initialize an invalid states") { + auto state = graph.empty_state(); + std::vector init_values{1, 1}; + CHECK_THROWS_WITH(bnode_ptr->initialize_state(state, init_values), + "Initialized values do not satisfy sum constraint(s)."); + } + } + + GIVEN("(2)-BinaryNode with a sum constraint over the entire array") { + auto graph = Graph(); + std::vector sum_constraints{{std::nullopt, {GreaterEqual}, {1}}}; + auto bnode_ptr = + graph.emplace_node(2, std::nullopt, std::nullopt, sum_constraints); + + WHEN("We initialize an invalid states") { + auto state = graph.empty_state(); + std::vector init_values{0, 0}; + CHECK_THROWS_WITH(bnode_ptr->initialize_state(state, init_values), + "Initialized values do not satisfy sum constraint(s)."); + } + } + + GIVEN("(2)-BinaryNode with a sum constraint over the entire array") { + auto graph = Graph(); + std::vector sum_constraints{{std::nullopt, {LessEqual}, {1}}}; + auto bnode_ptr = + graph.emplace_node(2, std::nullopt, std::nullopt, sum_constraints); + + WHEN("We initialize an invalid states") { + auto state = graph.empty_state(); + std::vector init_values{1, 1}; + CHECK_THROWS_WITH(bnode_ptr->initialize_state(state, init_values), + "Initialized values do not satisfy sum constraint(s)."); + } + } + + GIVEN("(2x2x2)-BinaryNode with a sum constraint over the entire array") { + auto graph = Graph(); + std::vector sum_constraints{{std::nullopt, {LessEqual}, {5}}}; + auto bnode_ptr = + graph.emplace_node(std::initializer_list{2, 2, 2}, + std::nullopt, std::nullopt, sum_constraints); + + THEN("Sum constraint is correct") { + CHECK(bnode_ptr->sum_constraints().size() == 1); + SumConstraint bnode_sum_constraint = bnode_ptr->sum_constraints()[0]; + CHECK(bnode_sum_constraint.axis() == std::nullopt); + CHECK(bnode_sum_constraint.num_bounds() == 1); + CHECK(bnode_sum_constraint.get_bound(0) == 5.0); + CHECK(bnode_sum_constraint.num_operators() == 1); + CHECK(bnode_sum_constraint.get_operator(0) == LessEqual); + } + + WHEN("We create a state using a random number generator") { + auto state = graph.empty_state(); + auto rng = std::default_random_engine(42); + CHECK_THROWS_WITH(bnode_ptr->initialize_state(state, rng), + "Cannot randomly initialize_state with sum constraints."); + } + + WHEN("We initialize a valid state") { + auto state = graph.empty_state(); + std::vector init_values{0, 0, 0, 1, 1, 0, 0, 0}; + bnode_ptr->initialize_state(state, init_values); + graph.initialize_state(state); + + auto sum_constraints_lhs = bnode_ptr->sum_constraints_lhs(state); + + THEN("Sum constraint sums and state are correct") { + CHECK(bnode_ptr->sum_constraints_lhs(state).size() == 1); + CHECK(bnode_ptr->sum_constraints_lhs(state).data()[0].size() == 1); + CHECK_THAT(bnode_ptr->sum_constraints_lhs(state)[0], RangeEquals({2.0})); + CHECK_THAT(bnode_ptr->view(state), RangeEquals(init_values)); + } + + THEN("We exchange() some values") { + bnode_ptr->exchange(state, 0, 1); // Does nothing. + bnode_ptr->exchange(state, 2, 3); + std::swap(init_values[0], init_values[1]); + std::swap(init_values[2], init_values[3]); + // state is now: [0, 0, 1, 0, 1, 0, 0, 0] + + THEN("Sum constraint sums and state updated correctly") { + CHECK_THAT(bnode_ptr->sum_constraints_lhs(state)[0], RangeEquals({2.0})); + CHECK(bnode_ptr->diff(state).size() == 2); // 2 updates per exchange + CHECK_THAT(bnode_ptr->view(state), RangeEquals(init_values)); + } + + AND_WHEN("We revert") { + graph.revert(state); + + THEN("Sum constraint sums reverted correctly") { + CHECK_THAT(bnode_ptr->sum_constraints_lhs(state)[0], RangeEquals({2.0})); + CHECK(bnode_ptr->diff(state).size() == 0); + } + } + } + } + } + + GIVEN("(3x2x2)-BinaryNode with a sum constraint on axis: 0") { + auto graph = Graph(); + std::vector sum_constraints{ + {0, {Equal, LessEqual, GreaterEqual}, {1.0, 2.0, 3.0}}}; + auto bnode_ptr = + graph.emplace_node(std::initializer_list{3, 2, 2}, + std::nullopt, std::nullopt, sum_constraints); + + THEN("Sum constraint is correct") { + CHECK(bnode_ptr->sum_constraints().size() == 1); + SumConstraint bnode_sum_constraint = bnode_ptr->sum_constraints()[0]; + CHECK(bnode_sum_constraint.axis() == 0); + CHECK(bnode_sum_constraint.num_bounds() == 3); + CHECK(bnode_sum_constraint.get_bound(0) == 1.0); + CHECK(bnode_sum_constraint.get_bound(1) == 2.0); + CHECK(bnode_sum_constraint.get_bound(2) == 3.0); + CHECK(bnode_sum_constraint.num_operators() == 3); + CHECK(bnode_sum_constraint.get_operator(0) == Equal); + CHECK(bnode_sum_constraint.get_operator(1) == LessEqual); + CHECK(bnode_sum_constraint.get_operator(2) == GreaterEqual); + } + + WHEN("We create a state using a random number generator") { + auto state = graph.empty_state(); + auto rng = std::default_random_engine(42); + CHECK_THROWS_WITH(bnode_ptr->initialize_state(state, rng), + "Cannot randomly initialize_state with sum constraints."); + } + + WHEN("We initialize three invalid states") { + auto state = graph.empty_state(); + // This state violates slice 0 along axis 0 + std::vector init_values{1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1}; + // import numpy as np + // a = np.asarray([1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1]) + // a = a.reshape(3, 2, 2) + // a.sum(axis=(1, 2)) + // >>> array([2, 2, 4]) + CHECK_THROWS_WITH(bnode_ptr->initialize_state(state, init_values), + "Initialized values do not satisfy sum constraint(s)."); + + state = graph.empty_state(); + // This state violates the slice 1 along axis 0 + init_values = {0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1}; + // import numpy as np + // a = np.asarray([0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1]) + // a = a.reshape(3, 2, 2) + // a.sum(axis=(1, 2)) + // >>> array([1, 3, 4]) + CHECK_THROWS_WITH(bnode_ptr->initialize_state(state, init_values), + "Initialized values do not satisfy sum constraint(s)."); + + state = graph.empty_state(); + // This state violates the slice 2 along axis 0 + init_values = {0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0}; + // import numpy as np + // a = np.asarray([0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0]) + // a = a.reshape(3, 2, 2) + // a.sum(axis=(1, 2)) + // >>> array([1, 2, 2]) + CHECK_THROWS_WITH(bnode_ptr->initialize_state(state, init_values), + "Initialized values do not satisfy sum constraint(s)."); + } + + WHEN("We initialize a valid state") { + auto state = graph.empty_state(); + std::vector init_values{0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1}; + bnode_ptr->initialize_state(state, init_values); + graph.initialize_state(state); + + auto sum_constraints_lhs = bnode_ptr->sum_constraints_lhs(state); + + THEN("Sum constraint sums and state are correct") { + // **Python Code 1** + // import numpy as np + // a = np.asarray([0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1]) + // a = a.reshape(3, 2, 2) + // a.sum(axis=(1, 2)) + // >>> array([1, 2, 4]) + CHECK(bnode_ptr->sum_constraints_lhs(state).size() == 1); + CHECK(bnode_ptr->sum_constraints_lhs(state).data()[0].size() == 3); + CHECK_THAT(bnode_ptr->sum_constraints_lhs(state)[0], RangeEquals({1, 2, 4})); + CHECK_THAT(bnode_ptr->view(state), RangeEquals(init_values)); + } + + THEN("We exchange() some values") { + bnode_ptr->exchange(state, 0, 3); // Does nothing. + bnode_ptr->exchange(state, 1, 6); // Does nothing. + bnode_ptr->exchange(state, 1, 3); + std::swap(init_values[0], init_values[3]); + std::swap(init_values[1], init_values[6]); + std::swap(init_values[1], init_values[3]); + // state is now: [0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1] + + THEN("Sum constraint sums and state updated correctly") { + // Cont. w/ Python code at **Python Code 1** + // a[np.unravel_index(1, a.shape)] = 0 + // a[np.unravel_index(3, a.shape)] = 1 + // a.sum(axis=(1, 2)) + // >>> array([1, 2, 4]) + CHECK_THAT(bnode_ptr->sum_constraints_lhs(state)[0], RangeEquals({1, 2, 4})); + CHECK(bnode_ptr->diff(state).size() == 2); // 2 updates per exchange + CHECK_THAT(bnode_ptr->view(state), RangeEquals(init_values)); + } + + AND_WHEN("We revert") { + graph.revert(state); + + THEN("Sum constraint sums reverted correctly") { + CHECK_THAT(bnode_ptr->sum_constraints_lhs(state)[0], + RangeEquals({1, 2, 4})); + CHECK(bnode_ptr->diff(state).size() == 0); + } + } + } + + THEN("We clip_and_set_value() some values") { + bnode_ptr->clip_and_set_value(state, 5, -1); // Does nothing. + bnode_ptr->clip_and_set_value(state, 7, -1); + bnode_ptr->clip_and_set_value(state, 9, 1); // Does nothing. + bnode_ptr->clip_and_set_value(state, 11, 0); + bnode_ptr->clip_and_set_value(state, 11, 1); + bnode_ptr->clip_and_set_value(state, 10, 0); + init_values[5] = 0; + init_values[7] = 0; + init_values[9] = 1; + init_values[11] = 1; + init_values[10] = 0; + // state is now: [0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1] + + THEN("Sum constraint sums and state updated correctly") { + // Cont. w/ Python code at **Python Code 1** + // a[np.unravel_index(5, a.shape)] = 0 + // a[np.unravel_index(7, a.shape)] = 0 + // a[np.unravel_index(9, a.shape)] = 1 + // a[np.unravel_index(11, a.shape)] = 1 + // a[np.unravel_index(10, a.shape)] = 0 + // a.sum(axis=(1, 2)) + // >>> array([1, 1, 3]) + CHECK_THAT(bnode_ptr->sum_constraints_lhs(state)[0], RangeEquals({1, 1, 3})); + CHECK(bnode_ptr->diff(state).size() == 4); + CHECK_THAT(bnode_ptr->view(state), RangeEquals(init_values)); + } + + AND_WHEN("We revert") { + graph.revert(state); + + THEN("Sum constraint sums reverted correctly") { + CHECK_THAT(bnode_ptr->sum_constraints_lhs(state)[0], + RangeEquals({1, 2, 4})); + CHECK(bnode_ptr->diff(state).size() == 0); + } + } + } + + THEN("We set_value() some values") { + bnode_ptr->set_value(state, 0, 0); // Does nothing. + bnode_ptr->set_value(state, 6, 0); + bnode_ptr->set_value(state, 7, 0); + bnode_ptr->set_value(state, 4, 1); + bnode_ptr->set_value(state, 10, 1); // Does nothing. + bnode_ptr->set_value(state, 11, 0); + init_values[0] = 0; + init_values[6] = 0; + init_values[7] = 0; + init_values[4] = 1; + init_values[10] = 1; + init_values[11] = 0; + // state is now: [0, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0] + + THEN("Sum constraint sums and state updated correctly") { + // Cont. w/ Python code at **Python Code 1** + // a[np.unravel_index(0, a.shape)] = 0 + // a[np.unravel_index(6, a.shape)] = 0 + // a[np.unravel_index(7, a.shape)] = 0 + // a[np.unravel_index(4, a.shape)] = 1 + // a[np.unravel_index(10, a.shape)] = 1 + // a[np.unravel_index(11, a.shape)] = 0 + // a.sum(axis=(1, 2)) + // >>> array([1, 1, 3]) + CHECK_THAT(bnode_ptr->sum_constraints_lhs(state)[0], RangeEquals({1, 1, 3})); + CHECK(bnode_ptr->diff(state).size() == 4); + CHECK_THAT(bnode_ptr->view(state), RangeEquals(init_values)); + } + + AND_WHEN("We revert") { + graph.revert(state); + + THEN("Sum constraint sums reverted correctly") { + CHECK_THAT(bnode_ptr->sum_constraints_lhs(state)[0], + RangeEquals({1, 2, 4})); + CHECK(bnode_ptr->diff(state).size() == 0); + } + } + } + + THEN("We flip() some values") { + bnode_ptr->flip(state, 6); // 1 -> 0 + bnode_ptr->flip(state, 4); // 0 -> 1 + bnode_ptr->flip(state, 11); // 1 -> 0 + init_values[6] = !init_values[6]; + init_values[4] = !init_values[4]; + init_values[11] = !init_values[11]; + // state is now: [0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0] + + THEN("Sum constraint sums and state updated correctly") { + // Cont. w/ Python code at **Python Code 1** + // a[np.unravel_index(6, a.shape)] = 0 + // a[np.unravel_index(4, a.shape)] = 1 + // a[np.unravel_index(11, a.shape)] = 0 + // a.sum(axis=(1, 2)) + // >>> array([1, 2, 3]) + CHECK_THAT(bnode_ptr->sum_constraints_lhs(state)[0], RangeEquals({1, 2, 3})); + CHECK(bnode_ptr->diff(state).size() == 3); + CHECK_THAT(bnode_ptr->view(state), RangeEquals(init_values)); + } + + AND_WHEN("We revert") { + graph.revert(state); + + THEN("Sum constraint sums reverted correctly") { + CHECK_THAT(bnode_ptr->sum_constraints_lhs(state)[0], + RangeEquals({1, 2, 4})); + CHECK(bnode_ptr->diff(state).size() == 0); + } + } + } + + THEN("We unset() some values") { + bnode_ptr->unset(state, 0); // Does nothing. + bnode_ptr->unset(state, 6); + bnode_ptr->unset(state, 11); + init_values[0] = 0; + init_values[6] = 0; + init_values[11] = 0; + // state is now: [0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0] + + THEN("Sum constraint sums and state updated correctly") { + // Cont. w/ Python code at **Python Code 1** + // a[np.unravel_index(0, a.shape)] = 0 + // a[np.unravel_index(6, a.shape)] = 0 + // a[np.unravel_index(11, a.shape)] = 0 + // a.sum(axis=(1, 2)) + // >>> array([1, 1, 3]) + CHECK_THAT(bnode_ptr->sum_constraints_lhs(state)[0], RangeEquals({1, 1, 3})); + CHECK(bnode_ptr->diff(state).size() == 2); + CHECK_THAT(bnode_ptr->view(state), RangeEquals(init_values)); + } + + AND_WHEN("We commit and set() some values") { + graph.commit(state); + + bnode_ptr->set(state, 10); // Does nothing. + bnode_ptr->set(state, 11); + init_values[10] = 1; + init_values[11] = 1; + // state is now: [0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1] + + THEN("sum constraint sums updated correctly") { + CHECK_THAT(bnode_ptr->sum_constraints_lhs(state)[0], + RangeEquals({1, 1, 4})); + CHECK(bnode_ptr->diff(state).size() == 1); + CHECK_THAT(bnode_ptr->view(state), RangeEquals(init_values)); + } + + AND_WHEN("We revert") { + graph.revert(state); + + THEN("Sum constraint sums reverted correctly") { + CHECK_THAT(bnode_ptr->sum_constraints_lhs(state)[0], + RangeEquals({1, 1, 3})); + CHECK(bnode_ptr->diff(state).size() == 0); + } + } + } + } + } + } + // *********************** Sum Constraint tests ************************* } TEST_CASE("IntegerNode") { @@ -459,7 +1306,8 @@ TEST_CASE("IntegerNode") { } } - GIVEN("Double precision numbers, which may fall outside integer range or are not integral") { + GIVEN("Double precision numbers, which may fall outside integer range or are not " + "integral") { IntegerNode inode({1}); THEN("The state is not deterministic") { CHECK(!inode.deterministic_state()); } @@ -736,6 +1584,680 @@ TEST_CASE("IntegerNode") { REQUIRE_THROWS_WITH(graph.emplace_node(std::initializer_list{-1, 3}), "Number array cannot have dynamic size."); } + + // *********************** Sum Constraint tests ************************* + GIVEN("(2x3)-IntegerNode with a sum constraint on the invalid axis -2") { + std::vector sum_constraints{{-2, {Equal}, {20.0}}}; + + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{2, 3}, std::nullopt, + std::nullopt, sum_constraints), + "Invalid constrained axis given number array shape."); + } + + GIVEN("(2x3x4)-IntegerNode with a sum constraint on the invalid axis 3") { + std::vector sum_constraints{{3, {Equal}, {10.0}}}; + + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{2, 3, 4}, + std::nullopt, std::nullopt, sum_constraints), + "Invalid constrained axis given number array shape."); + } + + GIVEN("(2x3x4)-IntegerNode with a sum constraint on axis: 1 with too many operators.") { + std::vector sum_constraints{{1, {LessEqual, Equal, Equal, Equal}, {-10.0}}}; + + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{2, 3, 4}, + std::nullopt, std::nullopt, sum_constraints), + "Invalid number of operators given number array shape."); + } + + GIVEN("(2x3x4)-IntegerNode with a sum constraint on axis: 1 with too few operators.") { + std::vector sum_constraints{{1, {LessEqual, Equal}, {-11.0}}}; + + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{2, 3, 4}, + std::nullopt, std::nullopt, sum_constraints), + "Invalid number of operators given number array shape."); + } + + GIVEN("(2x3x4)-IntegerNode with a sum constraint on axis: 1 with too many bounds.") { + std::vector sum_constraints{{1, {LessEqual}, {-10.0, 20.0, 30.0, 40.0}}}; + + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{2, 3, 4}, + std::nullopt, std::nullopt, sum_constraints), + "Invalid number of bounds given number array shape."); + } + + GIVEN("(2x3x4)-IntegerNode with a sum constraint on axis: 1 with too few bounds.") { + std::vector sum_constraints{{1, {LessEqual}, {111.0, -223.0}}}; + + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{2, 3, 4}, + std::nullopt, std::nullopt, sum_constraints), + "Invalid number of bounds given number array shape."); + } + + GIVEN("(6)-IntegerNode with duplicate sum constraints over the entire array") { + SumConstraint sum_constraint{std::nullopt, {Equal}, {10.0}}; + std::vector sum_constraints{sum_constraint, sum_constraint}; + + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{6}, std::nullopt, + std::nullopt, sum_constraints), + "Cannot define multiple sum constraints for the entire number array."); + } + + GIVEN("(2x3x4)-IntegerNode with duplicate sum constraints on axis: 1") { + std::vector sum_constraints{{1, {Equal}, {100.0}}, {1, {Equal}, {100.0}}}; + + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{2, 3, 4}, + std::nullopt, std::nullopt, sum_constraints), + "Cannot define multiple sum constraints for a single axis."); + } + + GIVEN("(2x3x4)-IntegerNode with sum constraints on axis: 1 and the entire array.") { + std::vector sum_constraints{{std::nullopt, {Equal}, {100.0}}, + {1, {Equal}, {100.0}}}; + + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{2, 3, 4}, + std::nullopt, std::nullopt, sum_constraints), + "Can define at most one sum constraint per number array."); + } + + GIVEN("(2x3x4)-IntegerNode with sum constraints on axes: 0 and 1") { + std::vector sum_constraints{{0, {Equal}, {100.0}}, {1, {Equal}, {100.0}}}; + + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{2, 3, 4}, + std::nullopt, std::nullopt, sum_constraints), + "Can define at most one sum constraint per number array."); + } + + GIVEN("(2x3x4)-IntegerNode with a non-integral sum constraint") { + std::vector sum_constraints{{1, {LessEqual}, {11.0, 12.0001, 0.0}}}; + + REQUIRE_THROWS_WITH( + graph.emplace_node(std::initializer_list{2, 3, 4}, + std::nullopt, std::nullopt, sum_constraints), + "Sum constraint(s) for integral arrays must be integral."); + } + + GIVEN("(6)-IntegerNode with an infeasible sum constraint over the entire array.") { + auto graph = Graph(); + std::vector sum_constraints{{std::nullopt, {Equal}, {-7.0}}}; + REQUIRE_THROWS_WITH(graph.emplace_node(6, -1, 8, sum_constraints), + "Infeasible sum constraint."); + } + + GIVEN("(6)-IntegerNode with an infeasible sum constraint over the entire array.") { + auto graph = Graph(); + std::vector sum_constraints{{std::nullopt, {LessEqual}, {-7.0}}}; + REQUIRE_THROWS_WITH(graph.emplace_node(6, -1, 8, sum_constraints), + "Infeasible sum constraint."); + } + + GIVEN("(6)-IntegerNode with an infeasible sum constraint over the entire array.") { + auto graph = Graph(); + std::vector sum_constraints{{std::nullopt, {GreaterEqual}, {13}}}; + REQUIRE_THROWS_WITH(graph.emplace_node(6, -1, 2, sum_constraints), + "Infeasible sum constraint."); + } + + GIVEN("(2x3x2)-IntegerNode with an infeasible sum constraint on axis: 0") { + auto graph = Graph(); + std::vector sum_constraints{{0, {Equal, LessEqual}, {5.0, -31.0}}}; + // Each slice along axis 0 has size 6. There is no feasible assignment + // to the values in slice 1 (along axis 0) that results in a sum less + // than or equal to -5*6 - 1 = -31. + REQUIRE_THROWS_WITH(graph.emplace_node(std::initializer_list{2, 3, 2}, + -5, 8, sum_constraints), + "Infeasible sum constraint."); + } + + GIVEN("(2x3x2)-IntegerNode with an infeasible sum constraint on axis: 1") { + auto graph = Graph(); + std::vector sum_constraints{ + {1, {GreaterEqual, Equal, Equal}, {33.0, 0.0, 0.0}}}; + // Each slice along axis 1 has size 4. There is no feasible assignment + // to the values in slice 0 (along axis 1) that results in a sum + // greater than or equal to 4*8 + 1 = 33. + REQUIRE_THROWS_WITH(graph.emplace_node(std::initializer_list{2, 3, 2}, + -5, 8, sum_constraints), + "Infeasible sum constraint."); + } + + GIVEN("(2x3x2)-IntegerNode with an infeasible sum constraint on axis: 2") { + auto graph = Graph(); + std::vector sum_constraints{{2, {GreaterEqual, Equal}, {-1.0, 49.0}}}; + // Each slice along axis 2 has size 6. There is no feasible assignment + // to the values in slice 1 (along axis 2) that results in a sum or + // equal to 6*8 + 1 = 49 + REQUIRE_THROWS_WITH(graph.emplace_node(std::initializer_list{2, 3, 2}, + -5, 8, sum_constraints), + "Infeasible sum constraint."); + } + + GIVEN("(2x2x2)-IntegerNode with a feasible sum constraint over the entire array ") { + auto graph = Graph(); + std::vector sum_constraints{{std::nullopt, {GreaterEqual}, {40}}}; + auto inode_ptr = graph.emplace_node(std::initializer_list{2, 2, 2}, + -5, 8, sum_constraints); + + THEN("Sum constraint is correct") { + CHECK(inode_ptr->sum_constraints().size() == 1); + SumConstraint inode_sum_constraint = inode_ptr->sum_constraints()[0]; + CHECK(inode_sum_constraint.axis() == std::nullopt); + CHECK(inode_sum_constraint.num_bounds() == 1); + CHECK(inode_sum_constraint.get_bound(0) == 40.0); + CHECK(inode_sum_constraint.num_operators() == 1); + CHECK(inode_sum_constraint.get_operator(0) == GreaterEqual); + } + + WHEN("We create a state by initialize_state()") { + auto state = graph.initialize_state(); + graph.initialize_state(state); + std::vector expected_init{8, 8, 8, 8, 8, 8, -3, -5}; + auto sum_constraints_lhs = inode_ptr->sum_constraints_lhs(state); + + THEN("Sum constraint sums and state are correct") { + CHECK(inode_ptr->sum_constraints_lhs(state).size() == 1); + CHECK(inode_ptr->sum_constraints_lhs(state).data()[0].size() == 1); + CHECK_THAT(inode_ptr->sum_constraints_lhs(state)[0], RangeEquals({40})); + CHECK_THAT(inode_ptr->view(state), RangeEquals(expected_init)); + } + } + } + + GIVEN("(2x3x2)-IntegerNode with a feasible sum constraint on axis: 0") { + auto graph = Graph(); + std::vector sum_constraints{{0, {Equal, GreaterEqual}, {-21.0, 9.0}}}; + auto inode_ptr = graph.emplace_node(std::initializer_list{2, 3, 2}, + -5, 8, sum_constraints); + + THEN("Sum constraint is correct") { + CHECK(inode_ptr->sum_constraints().size() == 1); + SumConstraint inode_sum_constraint = inode_ptr->sum_constraints()[0]; + CHECK(inode_sum_constraint.axis() == 0); + CHECK(inode_sum_constraint.num_bounds() == 2); + CHECK(inode_sum_constraint.get_bound(0) == -21.0); + CHECK(inode_sum_constraint.get_bound(1) == 9.0); + CHECK(inode_sum_constraint.num_operators() == 2); + CHECK(inode_sum_constraint.get_operator(0) == Equal); + CHECK(inode_sum_constraint.get_operator(1) == GreaterEqual); + } + + WHEN("We create a state by initialize_state()") { + auto state = graph.initialize_state(); + graph.initialize_state(state); + // import numpy as np + // a = np.asarray([i for i in range(2*3*2)]).reshape(2, 3, 2) + // print(a[0, :, :].flatten()) + // >>> [0 1 2 3 4 5] + // print(a[1, :, :].flatten()) + // >>> [ 6 7 8 9 10 11] + // + // The method `construct_state_given_exactly_one_sum_constraint()` + // will construct a state as follows: + // [-5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5] + // repair slice 0 + // [4, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5] + // repair slice 1 + // [4, -5, -5, -5, -5, -5, 8, 8, 8, -5, -5, -5] + std::vector expected_init{4, -5, -5, -5, -5, -5, 8, 8, 8, -5, -5, -5}; + auto sum_constraints_lhs = inode_ptr->sum_constraints_lhs(state); + + THEN("Sum constraint sums and state are correct") { + CHECK(inode_ptr->sum_constraints_lhs(state).size() == 1); + CHECK(inode_ptr->sum_constraints_lhs(state).data()[0].size() == 2); + CHECK_THAT(inode_ptr->sum_constraints_lhs(state)[0], RangeEquals({-21.0, 9.0})); + CHECK_THAT(inode_ptr->view(state), RangeEquals(expected_init)); + } + } + } + + GIVEN("(2x3x2)-IntegerNode with a feasible sum constraint on axis: 1") { + auto graph = Graph(); + std::vector sum_constraints{ + {1, {Equal, GreaterEqual, LessEqual}, {0.0, -2.0, 0.0}}}; + auto inode_ptr = graph.emplace_node(std::initializer_list{2, 3, 2}, + -5, 8, sum_constraints); + + THEN("Sum constraint is correct") { + CHECK(inode_ptr->sum_constraints().size() == 1); + SumConstraint inode_sum_constraint = inode_ptr->sum_constraints()[0]; + CHECK(inode_sum_constraint.axis() == 1); + CHECK(inode_sum_constraint.num_bounds() == 3); + CHECK(inode_sum_constraint.get_bound(0) == 0.0); + CHECK(inode_sum_constraint.get_bound(1) == -2.0); + CHECK(inode_sum_constraint.get_bound(2) == 0.0); + CHECK(inode_sum_constraint.num_operators() == 3); + CHECK(inode_sum_constraint.get_operator(0) == Equal); + CHECK(inode_sum_constraint.get_operator(1) == GreaterEqual); + CHECK(inode_sum_constraint.get_operator(2) == LessEqual); + } + + WHEN("We create a state by initialize_state()") { + auto state = graph.initialize_state(); + graph.initialize_state(state); + // import numpy as np + // a = np.asarray([i for i in range(2*3*2)]).reshape(2, 3, 2) + // print(a[:, 0, :].flatten()) + // >>> [0 1 6 7] + // print(a[:, 1, :].flatten()) + // >>> [2 3 8 9] + // print(a[:, 2, :].flatten()) + // >>> [ 4 5 10 11] + // + // The method `construct_state_given_exactly_one_sum_constraint()` + // will construct a state as follows: + // [-5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5] + // repair slice 0 w/ [8, 2, -5, -5] + // [8, 2, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5] + // repair slice 1 w/ [8, 0, -5, -5] + // [8, 2, 8, 0, -5, -5, -5, -5, -5, -5, -5, -5] + // no need to repair slice 2 + std::vector expected_init{8, 2, 8, 0, -5, -5, -5, -5, -5, -5, -5, -5}; + auto sum_constraints_lhs = inode_ptr->sum_constraints_lhs(state); + + THEN("Sum constraint sums and state are correct") { + CHECK(inode_ptr->sum_constraints_lhs(state).size() == 1); + CHECK(inode_ptr->sum_constraints_lhs(state).data()[0].size() == 3); + CHECK_THAT(inode_ptr->sum_constraints_lhs(state)[0], + RangeEquals({0.0, -2.0, -20.0})); + CHECK_THAT(inode_ptr->view(state), RangeEquals(expected_init)); + } + } + } + + GIVEN("(2x3x2)-IntegerNode with a feasible sum constraint on axis: 2") { + auto graph = Graph(); + std::vector sum_constraints{{2, {Equal, GreaterEqual}, {23.0, 14.0}}}; + auto inode_ptr = graph.emplace_node(std::initializer_list{2, 3, 2}, + -5, 8, sum_constraints); + + THEN("Sum constraint is correct") { + CHECK(inode_ptr->sum_constraints().size() == 1); + SumConstraint inode_sum_constraint = inode_ptr->sum_constraints()[0]; + CHECK(inode_sum_constraint.axis() == 2); + CHECK(inode_sum_constraint.num_bounds() == 2); + CHECK(inode_sum_constraint.get_bound(0) == 23.0); + CHECK(inode_sum_constraint.get_bound(1) == 14.0); + CHECK(inode_sum_constraint.num_operators() == 2); + CHECK(inode_sum_constraint.get_operator(0) == Equal); + CHECK(inode_sum_constraint.get_operator(1) == GreaterEqual); + } + + WHEN("We create a state by initialize_state()") { + auto state = graph.initialize_state(); + graph.initialize_state(state); + // import numpy as np + // a = np.asarray([i for i in range(2*3*2)]).reshape(2, 3, 2) + // print(a[:, :, 0].flatten()) + // >>> [ 0 2 4 6 8 10] + // print(a[:, :, 1].flatten()) + // >>> [ 1 3 5 7 9 11] + // + // The method `construct_state_given_exactly_one_sum_constraint()` + // will construct a state as follows: + // [-5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5] + // repair slice 0 w/ [8, 8, 8, 8, -4, -5] + // [8, -5, 8, -5, 8, -5, 8, -5, -4, -5, -5, -5] + // repair slice 0 w/ [8, 8, 8, 0, -5, -5] + // [8, 8, 8, 8, 8, 8, 8, 0, -4, -5, -5, -5] + std::vector expected_init{8, 8, 8, 8, 8, 8, 8, 0, -4, -5, -5, -5}; + auto sum_constraints_lhs = inode_ptr->sum_constraints_lhs(state); + + THEN("Sum constraint sums and state are correct") { + CHECK(inode_ptr->sum_constraints_lhs(state).size() == 1); + CHECK(inode_ptr->sum_constraints_lhs(state).data()[0].size() == 2); + CHECK_THAT(inode_ptr->sum_constraints_lhs(state)[0], RangeEquals({23.0, 14.0})); + CHECK_THAT(inode_ptr->view(state), RangeEquals(expected_init)); + } + } + } + + GIVEN("(2)-IntegerNode with a sum constraint over the entire array") { + auto graph = Graph(); + std::vector sum_constraints{{std::nullopt, {Equal}, {15}}}; + auto inode_ptr = graph.emplace_node(2, -5, 8, sum_constraints); + + WHEN("We initialize two invalid states") { + auto state = graph.empty_state(); + std::vector init_values{0.0, 0.0}; + CHECK_THROWS_WITH(inode_ptr->initialize_state(state, init_values), + "Initialized values do not satisfy sum constraint(s)."); + + state = graph.empty_state(); + init_values = {8.0, 8.0}; + CHECK_THROWS_WITH(inode_ptr->initialize_state(state, init_values), + "Initialized values do not satisfy sum constraint(s)."); + } + } + + GIVEN("(2)-IntegerNode with a sum constraint over the entire array") { + auto graph = Graph(); + std::vector sum_constraints{{std::nullopt, {LessEqual}, {10}}}; + auto inode_ptr = graph.emplace_node(2, -5, 8, sum_constraints); + + WHEN("We initialize an invalid states") { + auto state = graph.empty_state(); + std::vector init_values{8.0, 7.0}; + CHECK_THROWS_WITH(inode_ptr->initialize_state(state, init_values), + "Initialized values do not satisfy sum constraint(s)."); + } + } + + GIVEN("(2)-IntegerNode with a sum constraint over the entire array") { + auto graph = Graph(); + std::vector sum_constraints{{std::nullopt, {GreaterEqual}, {10}}}; + auto inode_ptr = graph.emplace_node(2, -5, 8, sum_constraints); + + WHEN("We initialize an invalid states") { + auto state = graph.empty_state(); + std::vector init_values{-5.0, -4.0}; + CHECK_THROWS_WITH(inode_ptr->initialize_state(state, init_values), + "Initialized values do not satisfy sum constraint(s)."); + } + } + + GIVEN("(2x2)-BinaryNode with a sum constraint over the entire array") { + auto graph = Graph(); + std::vector sum_constraints{{std::nullopt, {GreaterEqual}, {5.0}}}; + auto inode_ptr = graph.emplace_node(std::initializer_list{2, 2}, -5, + 8, sum_constraints); + + THEN("Sum constraint is correct") { + CHECK(inode_ptr->sum_constraints().size() == 1); + SumConstraint inode_sum_constraint = inode_ptr->sum_constraints()[0]; + CHECK(inode_sum_constraint.axis() == std::nullopt); + CHECK(inode_sum_constraint.num_bounds() == 1); + CHECK(inode_sum_constraint.get_bound(0) == 5.0); + CHECK(inode_sum_constraint.num_operators() == 1); + CHECK(inode_sum_constraint.get_operator(0) == GreaterEqual); + } + + WHEN("We create a state using a random number generator") { + auto state = graph.empty_state(); + auto rng = std::default_random_engine(42); + CHECK_THROWS_WITH(inode_ptr->initialize_state(state, rng), + "Cannot randomly initialize_state with sum constraints."); + } + + WHEN("We initialize a valid state") { + auto state = graph.empty_state(); + std::vector init_values{1.0, -1.0, 0.0, 5.0}; + inode_ptr->initialize_state(state, init_values); + graph.initialize_state(state); + + auto sum_constraints_lhs = inode_ptr->sum_constraints_lhs(state); + + THEN("Sum constraint sums and state are correct") { + CHECK(inode_ptr->sum_constraints_lhs(state).size() == 1); + CHECK(inode_ptr->sum_constraints_lhs(state).data()[0].size() == 1); + CHECK_THAT(inode_ptr->sum_constraints_lhs(state)[0], RangeEquals({5.0})); + CHECK_THAT(inode_ptr->view(state), RangeEquals(init_values)); + } + + THEN("We set_value() some values") { + inode_ptr->set_value(state, 0, 1); // Does nothing. + inode_ptr->set_value(state, 2, 3); + init_values[0] = 1; + init_values[2] = 3; + // state is now: [1.0, -1.0, 3.0, 5.0] + + THEN("Sum constraint sums and state updated correctly") { + CHECK_THAT(inode_ptr->sum_constraints_lhs(state)[0], RangeEquals({8.0})); + CHECK(inode_ptr->diff(state).size() == 1); + CHECK_THAT(inode_ptr->view(state), RangeEquals(init_values)); + } + + AND_WHEN("We revert") { + graph.revert(state); + + THEN("Sum constraint sums reverted correctly") { + CHECK_THAT(sum_constraints_lhs[0], RangeEquals({5.0})); + CHECK(inode_ptr->diff(state).size() == 0); + } + } + } + } + } + + GIVEN("(2x3x2)-IntegerNode with index-wise bounds and a sum constraint on axis: 1") { + auto graph = Graph(); + std::vector sum_constraints{ + {1, {Equal, LessEqual, GreaterEqual}, {11.0, 2.0, 5.0}}}; + auto inode_ptr = graph.emplace_node(std::initializer_list{2, 3, 2}, + -5, 8, sum_constraints); + + THEN("Sum constraint is correct") { + CHECK(inode_ptr->sum_constraints().size() == 1); + SumConstraint inode_sum_constraint = inode_ptr->sum_constraints()[0]; + CHECK(inode_sum_constraint.axis() == 1); + CHECK(inode_sum_constraint.num_bounds() == 3); + CHECK(inode_sum_constraint.get_bound(0) == 11.0); + CHECK(inode_sum_constraint.get_bound(1) == 2.0); + CHECK(inode_sum_constraint.get_bound(2) == 5.0); + CHECK(inode_sum_constraint.num_operators() == 3); + CHECK(inode_sum_constraint.get_operator(0) == Equal); + CHECK(inode_sum_constraint.get_operator(1) == LessEqual); + CHECK(inode_sum_constraint.get_operator(2) == GreaterEqual); + } + + WHEN("We create a state using a random number generator") { + auto state = graph.empty_state(); + auto rng = std::default_random_engine(42); + CHECK_THROWS_WITH(inode_ptr->initialize_state(state, rng), + "Cannot randomly initialize_state with sum constraints."); + } + + WHEN("We initialize three invalid states") { + auto state = graph.empty_state(); + // This state violates the slice 0 along axis 1 + std::vector init_values{5, 6, 0, 0, 3, 1, 4, 0, 2, 0, 0, 3}; + // import numpy as np + // a = np.asarray([5, 6, 0, 0, 3, 1, 4, 0, 2, 0, 0, 3]) + // a = a.reshape(2, 3, 2) + // a.sum(axis=(0, 2)) + // >>> array([15, 2, 7]) + CHECK_THROWS_WITH(inode_ptr->initialize_state(state, init_values), + "Initialized values do not satisfy sum constraint(s)."); + + state = graph.empty_state(); + // This state violates the slice 1 along axis 1 + init_values = {5, 2, 0, 0, 3, 1, 4, 0, 2, 1, 0, 3}; + // import numpy as np + // a = np.asarray([5, 2, 0, 0, 3, 1, 4, 0, 2, 1, 0, 3]) + // a = a.reshape(2, 3, 2) + // a.sum(axis=(0, 2)) + // >>> array([11, 3, 7]) + CHECK_THROWS_WITH(inode_ptr->initialize_state(state, init_values), + "Initialized values do not satisfy sum constraint(s)."); + + state = graph.empty_state(); + // This state violates the slice 2 along axis 1 + init_values = {5, 2, 0, 0, 3, 1, 4, 0, 1, 0, 0, 0}; + // import numpy as np + // a = np.asarray([5, 2, 0, 0, 3, 1, 4, 0, 1, 0, 0, 0]) + // a = a.reshape(2, 3, 2) + // a.sum(axis=(0, 2)) + // >>> array([11, 1, 4]) + CHECK_THROWS_WITH(inode_ptr->initialize_state(state, init_values), + "Initialized values do not satisfy sum constraint(s)."); + } + + WHEN("We initialize a valid state") { + auto state = graph.empty_state(); + std::vector init_values{5, 2, 0, 0, 3, 1, 4, 0, 2, 0, 0, 3}; + inode_ptr->initialize_state(state, init_values); + graph.initialize_state(state); + + auto sum_constraints_lhs = inode_ptr->sum_constraints_lhs(state); + + THEN("Sum constraint sums and state are correct") { + // **Python Code 2** + // import numpy as np + // a = np.asarray([5, 2, 0, 0, 3, 1, 4, 0, 2, 0, 0, 3]) + // a = a.reshape(2, 3, 2) + // a.sum(axis=(0, 2)) + // >>> array([11, 2, 7]) + CHECK(inode_ptr->sum_constraints_lhs(state).size() == 1); + CHECK(inode_ptr->sum_constraints_lhs(state).data()[0].size() == 3); + CHECK_THAT(inode_ptr->sum_constraints_lhs(state)[0], RangeEquals({11, 2, 7})); + CHECK_THAT(inode_ptr->view(state), RangeEquals(init_values)); + } + + THEN("We exchange() some values") { + inode_ptr->exchange(state, 2, 3); // Does nothing. + inode_ptr->exchange(state, 1, 8); // Does nothing. + inode_ptr->exchange(state, 8, 10); + inode_ptr->exchange(state, 0, 1); + std::swap(init_values[2], init_values[3]); + std::swap(init_values[1], init_values[8]); + std::swap(init_values[8], init_values[10]); + std::swap(init_values[0], init_values[1]); + // state is now: [2, 5, 0, 0, 3, 1, 4, 0, 0, 0, 2, 3] + + THEN("Sum constraint sums and state updated correctly") { + // Cont. w/ Python code at **Python Code 2** + // a[np.unravel_index(8, a.shape)] = 0 + // a[np.unravel_index(10, a.shape)] = 2 + // a[np.unravel_index(0, a.shape)] = 2 + // a[np.unravel_index(1, a.shape)] = 5 + // a.sum(axis=(0, 2)) + // >>> array([11, 0, 9]) + CHECK_THAT(inode_ptr->sum_constraints_lhs(state)[0], RangeEquals({11, 0, 9})); + CHECK(inode_ptr->diff(state).size() == 4); // 2 updates per exchange + CHECK_THAT(inode_ptr->view(state), RangeEquals(init_values)); + } + + AND_WHEN("We revert") { + graph.revert(state); + + THEN("Sum constraint sums reverted correctly") { + CHECK_THAT(inode_ptr->sum_constraints_lhs(state)[0], + RangeEquals({11, 2, 7})); + CHECK(inode_ptr->diff(state).size() == 0); + } + } + } + + THEN("We clip_and_set_value() some values") { + inode_ptr->clip_and_set_value(state, 0, 5); // Does nothing. + inode_ptr->clip_and_set_value(state, 8, -300); + inode_ptr->clip_and_set_value(state, 10, 100); + init_values[8] = -5; + init_values[10] = 8; + // state is now: [5, 2, 0, 0, 3, 1, 4, 0, -5, 0, 8, 3] + + THEN("Sum constraint sums and state updated correctly") { + // Cont. w/ Python code at **Python Code 2** + // a[np.unravel_index(8, a.shape)] = -5 + // a[np.unravel_index(10, a.shape)] = 8 + // a.sum(axis=(0, 2)) + // >>> array([11, -5, 15]) + CHECK_THAT(inode_ptr->sum_constraints_lhs(state)[0], RangeEquals({11, -5, 15})); + CHECK(inode_ptr->diff(state).size() == 2); + CHECK_THAT(inode_ptr->view(state), RangeEquals(init_values)); + } + + AND_WHEN("We revert") { + graph.revert(state); + + THEN("Sum constraint sums reverted correctly") { + CHECK_THAT(inode_ptr->sum_constraints_lhs(state)[0], + RangeEquals({11, 2, 7})); + CHECK(inode_ptr->diff(state).size() == 0); + } + } + } + + THEN("We set_value() some values") { + inode_ptr->set_value(state, 0, 5); // Does nothing. + inode_ptr->set_value(state, 8, 0); + inode_ptr->set_value(state, 9, 1); + inode_ptr->set_value(state, 10, 5); + inode_ptr->set_value(state, 11, 0); + init_values[0] = 5; + init_values[8] = 0; + init_values[9] = 1; + init_values[10] = 5; + init_values[11] = 0; + // state is now: [5, 2, 0, 0, 3, 1, 4, 0, 0, 1, 5, 0] + + THEN("Sum constraint sums and state updated correctly") { + // Cont. w/ Python code at **Python Code 2** + // a[np.unravel_index(0, a.shape)] = 5 + // a[np.unravel_index(8, a.shape)] = 0 + // a[np.unravel_index(9, a.shape)] = 1 + // a[np.unravel_index(10, a.shape)] = 5 + // a[np.unravel_index(11, a.shape)] = 0 + // a.sum(axis=(0, 2)) + // >>> array([11, 1, 9]) + CHECK_THAT(inode_ptr->sum_constraints_lhs(state)[0], RangeEquals({11, 1, 9})); + CHECK(inode_ptr->diff(state).size() == 4); + CHECK_THAT(inode_ptr->view(state), RangeEquals(init_values)); + } + + AND_WHEN("We revert") { + graph.revert(state); + + THEN("Sum constraint sums reverted correctly") { + CHECK_THAT(sum_constraints_lhs[0], RangeEquals({11, 2, 7})); + CHECK(inode_ptr->diff(state).size() == 0); + } + } + } + } + } + + GIVEN("(2x3)-IntegerNode and a sum constraint on axis: 0 with operator `==`") { + auto graph = Graph(); + std::vector sum_constraints{{0, {Equal}, {1.0}}}; + auto inode_ptr = graph.emplace_node( + std::initializer_list{2, 3}, std::nullopt, std::nullopt, sum_constraints); + + THEN("Sum constraint is correct") { + CHECK(inode_ptr->sum_constraints().size() == 1); + SumConstraint inode_sum_constraint = inode_ptr->sum_constraints()[0]; + CHECK(inode_sum_constraint.axis() == 0); + CHECK(inode_sum_constraint.num_bounds() == 1); + CHECK(inode_sum_constraint.get_bound(0) == 1.0); + CHECK(inode_sum_constraint.num_operators() == 1); + CHECK(inode_sum_constraint.get_operator(0) == Equal); + } + + WHEN("We initialize a valid state by construction") { + auto state = graph.empty_state(); + graph.initialize_state(state); + + auto sum_constraints_lhs = inode_ptr->sum_constraints_lhs(state); + + THEN("Sum constraint sums and state are correct") { + CHECK(inode_ptr->sum_constraints_lhs(state).size() == 1); + CHECK(inode_ptr->sum_constraints_lhs(state).data()[0].size() == 2); + CHECK_THAT(inode_ptr->sum_constraints_lhs(state)[0], RangeEquals({1.0, 1.0})); + CHECK_THAT(inode_ptr->view(state), RangeEquals({1, 0, 0, 1, 0, 0})); + } + + THEN("We exchange() some values") { + inode_ptr->exchange(state, 0, 1); + inode_ptr->exchange(state, 3, 4); + + THEN("Sum constraint sums and state updated correctly") { + CHECK_THAT(inode_ptr->sum_constraints_lhs(state)[0], RangeEquals({1.0, 1.0})); + CHECK(inode_ptr->diff(state).size() == 4); // 2 updates per exchange + CHECK_THAT(inode_ptr->view(state), RangeEquals({0, 1, 0, 0, 1, 0})); + } + } + } + } + // *********************** Sum Constraint tests ************************* } } // namespace dwave::optimization diff --git a/tests/test_symbols.py b/tests/test_symbols.py index 9220e288..25041fc1 100644 --- a/tests/test_symbols.py +++ b/tests/test_symbols.py @@ -711,7 +711,7 @@ def test(self): model.binary([10]) - def test_bounds(self): + def test_index_wise_bounds(self): model = Model() x = model.binary(lower_bound=0, upper_bound=1) self.assertEqual(x.lower_bound(), 0) @@ -730,6 +730,67 @@ def test_bounds(self): with self.assertRaises(ValueError): model.binary((2, 3), upper_bound=np.arange(6)) + def test_sum_constraints(self): + model = Model() + + # stores correct sum constraint + x = model.binary((2, 3), subject_to=[("==", 1)]) + self.assertEqual(x.sum_constraints(), [(["=="], [1])]) + x = model.binary((2, 3), subject_to=[("==", 1)]) + self.assertEqual(x.sum_constraints(), [(["=="], [1])]) + x = model.binary((2, 3), axes_subject_to=[(0, ["<=", "=="], [1, 2])]) + self.assertEqual(x.sum_constraints(), [(0, ["<=", "=="], [1, 2])]) + x = model.binary((2, 3), axes_subject_to=[(1, "<=", [1, 2, 1])]) + self.assertEqual(x.sum_constraints(), [(1, ["<="], [1, 2, 1])]) + x = model.binary((2, 3), axes_subject_to=[(0, ["<=", "=="], 1)]) + self.assertEqual(x.sum_constraints(), [(0, ["<=", "=="], [1])]) + x = model.binary((2, 3), axes_subject_to=[(0, "<=", 1)]) + self.assertEqual(x.sum_constraints(), [(0, ["<="], [1])]) + x = model.binary((2, 3), axes_subject_to=[(0, ["<=", "=="], np.asarray([1, 2]))]) + self.assertEqual(x.sum_constraints(), [(0, ["<=", "=="], [1, 2])]) + + # infeasible sum constraint + with self.assertRaises(ValueError): + model.binary((2, 3), upper_bound=[0, 1, 0, 0, 1, 0], subject_to=[(">=", 3)]) + with self.assertRaises(ValueError): + model.binary((2, 3), lower_bound=[0, 1, 0, 0, 1, 0], axes_subject_to=[(0, "==", 0)]) + with self.assertRaises(ValueError): + model.binary((2, 3), lower_bound=[0, 1, 0, 0, 1, 0], axes_subject_to=[(0, "<=", 0)]) + with self.assertRaises(ValueError): + model.binary((2, 3), upper_bound=[0, 1, 0, 0, 1, 0], axes_subject_to=[(0, ">=", 2)]) + + # incorrect number of operators and or bounds + with self.assertRaises(TypeError): + model.binary((2, 3), subject_to=[("==", [0, 0, 0])]) + with self.assertRaises(TypeError): + model.binary((2, 3), subject_to=[(["==", "<=", "=="], [0])]) + with self.assertRaises(ValueError): + model.binary((2, 3), axes_subject_to=[(0, "==", [0, 0, 0])]) + with self.assertRaises(ValueError): + model.binary((2, 3), axes_subject_to=[(0, ["==", "<=", "=="], [0, 0])]) + + # check bad argument format + with self.assertRaises(TypeError): + model.binary((2, 3), axes_subject_to=[(1.1, "<=", [0, 0, 0])]) + with self.assertRaises(TypeError): + model.binary((2, 3), axes_subject_to=[(1, 4, [0, 0, 0])]) + with self.assertRaises(TypeError): + model.binary((2, 3), axes_subject_to=[(1, ["!="], [0, 0, 0])]) + with self.assertRaises(TypeError): + model.binary((2, 3), axes_subject_to=[(1, ["=="], [[0, 0, 0]])]) + with self.assertRaises(TypeError): + model.binary((2, 3), axes_subject_to=[(1, [["<="]], [0, 0, 0])]) + with self.assertRaises(TypeError): + model.binary((2, 3), subject_to=[([["<="]], [0, 0, 0])]) + + # invalid number of sum constraints + with self.assertRaises(ValueError): + model.binary((2, 3), subject_to=[("==", 1), ("<=", 0)]) + with self.assertRaises(ValueError): + model.binary((2, 3), subject_to=[("==", 1)], axes_subject_to=[(1, "<=", [1, 1, 1])]) + with self.assertRaises(ValueError): + model.binary((2, 3), axes_subject_to=[(0, "==", 1), (1, "<=", [1, 1, 1])]) + def test_no_shape(self): model = Model() x = model.binary() @@ -765,6 +826,10 @@ def test_serialization(self): model.binary(), model.binary(3, lower_bound=1), model.binary(2, upper_bound=[0,1]), + model.binary(6, subject_to=[("<=", 2)]), + model.binary((2, 3), subject_to=[("<=", 2)]), + model.binary((2, 3), axes_subject_to=[(1, "<=", [0, 1, 2])]), + model.binary((2, 3), axes_subject_to=[(0, ["<=", "=="], 1)]), ] model.lock() @@ -776,6 +841,7 @@ def test_serialization(self): for i in range(old.size()): self.assertTrue(np.all(old.lower_bound() == new.lower_bound())) self.assertTrue(np.all(old.upper_bound() == new.upper_bound())) + self.assertEqual(old.sum_constraints(), new.sum_constraints()) def test_set_state(self): with self.subTest("array-like"): @@ -798,7 +864,7 @@ def test_set_state(self): with np.testing.assert_raises(ValueError): x.set_state(0, 2) - with self.subTest("Simple bounds test"): + with self.subTest("Simple index-wise bounds test"): model = Model() model.states.resize(1) x = model.binary(2, lower_bound=[-1, 0.9], upper_bound=[1.1, 1.2]) @@ -808,6 +874,32 @@ def test_set_state(self): with np.testing.assert_raises(ValueError): x.set_state(1, 0) + with self.subTest("Simple sum constraint test"): + model = Model() + model.states.resize(1) + + x = model.binary((2, 2), subject_to=[("<=", 1)]) + x.set_state(0, [0, 1, 0, 0]) + # Do not satisfy sum constraint + with np.testing.assert_raises(ValueError): + x.set_state(0, [1, 1, 0, 1]) + + x = model.binary((2, 3), axes_subject_to=[(0, "==", 1)]) + x.set_state(0, [0, 1, 0, 1, 0, 0]) + # Do not satisfy sum constraint + with np.testing.assert_raises(ValueError): + x.set_state(0, [1, 1, 0, 1, 0, 0]) + with np.testing.assert_raises(ValueError): + x.set_state(0, [0, 1, 0, 0, 0, 0]) + + x = model.binary((2, 2), axes_subject_to=[(1, ["<=", ">="], [0, 2])]) + x.set_state(0, [0, 1, 0, 1]) + # Do not satisfy sum constraint + with np.testing.assert_raises(ValueError): + x.set_state(0, [1, 1, 0, 1]) + with np.testing.assert_raises(ValueError): + x.set_state(0, [0, 0, 0, 1]) + with self.subTest("invalid state index"): model = Model() x = model.binary(5) @@ -1830,7 +1922,7 @@ def test_no_shape(self): model.states.resize(1) self.assertEqual(x.state(0).shape, tuple()) - def test_bounds(self): + def test_index_wise_bounds(self): model = Model() x = model.integer(lower_bound=4, upper_bound=5) self.assertEqual(x.lower_bound(), 4) @@ -1849,6 +1941,67 @@ def test_bounds(self): with self.assertRaises(ValueError): model.integer((2, 3), upper_bound=np.arange(6)) + def test_sum_constraints(self): + model = Model() + + # stores correct sum constraints + x = model.integer((2, 3), subject_to=[("==", 2)]) + self.assertEqual(x.sum_constraints(), [(["=="], [2])]) + x = model.integer((2, 3), subject_to=[("==", 2)]) + self.assertEqual(x.sum_constraints(), [(["=="], [2])]) + x = model.integer((2, 3), axes_subject_to=[(0, ["<=", "=="], [1, 2])]) + self.assertEqual(x.sum_constraints(), [(0, ["<=", "=="], [1, 2])]) + x = model.integer((2, 3), axes_subject_to=[(1, "<=", [1, 2, 1])]) + self.assertEqual(x.sum_constraints(), [(1, ["<="], [1, 2, 1])]) + x = model.integer((2, 3), axes_subject_to=[(0, ["<=", "=="], 1)]) + self.assertEqual(x.sum_constraints(), [(0, ["<=", "=="], [1])]) + x = model.integer((2, 3), axes_subject_to=[(0, "<=", 1)]) + self.assertEqual(x.sum_constraints(), [(0, ["<="], [1])]) + x = model.integer((2, 3), axes_subject_to=[(0, ["<=", "=="], np.asarray([1, 2]))]) + self.assertEqual(x.sum_constraints(), [(0, ["<=", "=="], [1, 2])]) + + # infeasible sum constraint + with self.assertRaises(ValueError): + model.integer((2, 2), upper_bound=2, subject_to=[(">=", 9)]) + with self.assertRaises(ValueError): + model.integer((2, 3), axes_subject_to=[(0, "==", -1)]) + with self.assertRaises(ValueError): + model.integer((2, 3), lower_bound=0, axes_subject_to=[(0, "<=", -1)]) + with self.assertRaises(ValueError): + model.integer((2, 3), upper_bound=2, axes_subject_to=[(0, ">=", 7)]) + + # incorrect number of operators and or bounds + with self.assertRaises(TypeError): + model.integer((2, 3), subject_to=[("==", [10, 20, 30])]) + with self.assertRaises(TypeError): + model.integer((2, 3), subject_to=[(["==", "<=", "=="], 10)]) + with self.assertRaises(ValueError): + model.integer((2, 3), axes_subject_to=[(0, "==", [10, 20, 30])]) + with self.assertRaises(ValueError): + model.integer((2, 3), axes_subject_to=[(0, ["==", "<=", "=="], [10, 20])]) + + # bad argument format + with self.assertRaises(TypeError): + model.integer((2, 3), subject_to=[([["=="]], 0)]) + with self.assertRaises(TypeError): + model.integer((2, 3), axes_subject_to=[(1.1, "<=", [0, 0, 0])]) + with self.assertRaises(TypeError): + model.integer((2, 3), axes_subject_to=[(1, 4, [0, 0, 0])]) + with self.assertRaises(TypeError): + model.integer((2, 3), axes_subject_to=[(1, ["!="], [0, 0, 0])]) + with self.assertRaises(TypeError): + model.integer((2, 3), axes_subject_to=[(1, ["=="], [[0, 0, 0]])]) + with self.assertRaises(TypeError): + model.integer((2, 3), axes_subject_to=[(1, [["=="]], [0, 0, 0])]) + + # invalid number of sum constraints + with self.assertRaises(ValueError): + model.integer((2, 3), subject_to=[("==", 1), ("<=", 0)]) + with self.assertRaises(ValueError): + model.integer((2, 3), subject_to=[("==", 1)], axes_subject_to=[(1, "<=", [1, 1, 1])]) + with self.assertRaises(ValueError): + model.integer((2, 3), axes_subject_to=[(0, "==", 1), (1, "<=", [1, 1, 1])]) + # Todo: we can generalize many of these tests for all decisions that can have # their state set @@ -1869,6 +2022,10 @@ def test_serialization(self): model.integer(upper_bound=105), model.integer(15, lower_bound=4, upper_bound=6), model.integer(2, lower_bound=[1, 2], upper_bound=[3, 4]), + model.integer(6, subject_to=[("<=", 2)]), + model.integer((2, 3), subject_to=[("<=", 2)]), + model.integer((2, 3), axes_subject_to=[(1, "<=", [0, 1, 2])]), + model.integer((2, 3), axes_subject_to=[(0, ["<=", ">="], 2)]), ] model.lock() @@ -1880,6 +2037,7 @@ def test_serialization(self): for i in range(old.size()): self.assertTrue(np.all(old.lower_bound() == new.lower_bound())) self.assertTrue(np.all(old.upper_bound() == new.upper_bound())) + self.assertEqual(old.sum_constraints(), new.sum_constraints()) def test_set_state(self): with self.subTest("Simple positive integer"): @@ -1904,7 +2062,7 @@ def test_set_state(self): with np.testing.assert_raises(ValueError): x.set_state(0, -1234) - with self.subTest("Simple bounds test"): + with self.subTest("Simple index-wise bounds test"): model = Model() model.states.resize(1) x = model.integer(1, lower_bound=-1, upper_bound=1) @@ -1915,6 +2073,34 @@ def test_set_state(self): with np.testing.assert_raises(ValueError): x.set_state(0, -2) + with self.subTest("Simple sum constraint test"): + model = Model() + model.states.resize(1) + + x = model.integer((2, 2), subject_to=[("==", 2)]) + x.set_state(0, [0, 2, 0, 0]) + # Do not satisfy sum constraint + with np.testing.assert_raises(ValueError): + x.set_state(0, [1, 1, 1, 1]) + with np.testing.assert_raises(ValueError): + x.set_state(0, [0, 0, 0, 1]) + + x = model.integer((2, 3), axes_subject_to=[(0, "==", 3)]) + x.set_state(0, [0, 3, 0, 1, 1, 1]) + # Do not satisfy sum constraint + with np.testing.assert_raises(ValueError): + x.set_state(0, [0, 3, 1, 1, 1, 1]) + with np.testing.assert_raises(ValueError): + x.set_state(0, [0, 3, 0, 1, 1, 0]) + + x = model.integer((2, 2), axes_subject_to=[(1, ["<=", ">="], [2, 6])]) + x.set_state(0, [1, 6, 1, 10]) + # Do not satisfy sum constraint + with np.testing.assert_raises(ValueError): + x.set_state(0, [1, 2, 1, 1]) + with np.testing.assert_raises(ValueError): + x.set_state(0, [1, 6, 2, 10]) + with self.subTest("array-like"): model = Model() model.states.resize(1)