diff --git a/dwave/optimization/include/dwave-optimization/cp/core/base_node.hpp b/dwave/optimization/include/dwave-optimization/cp/core/base_node.hpp index cbd05d8b7..2dc479742 100644 --- a/dwave/optimization/include/dwave-optimization/cp/core/base_node.hpp +++ b/dwave/optimization/include/dwave-optimization/cp/core/base_node.hpp @@ -14,12 +14,70 @@ #pragma once +#include + #include "dwave-optimization/array.hpp" #include "dwave-optimization/cp/core/cpvar.hpp" #include "dwave-optimization/cp/core/engine.hpp" #include "dwave-optimization/cp/core/interval_array.hpp" #include "dwave-optimization/cp/core/propagator.hpp" +#include "dwave-optimization/cp/propagators/binaryop.hpp" #include "dwave-optimization/cp/state/cpvar_state.hpp" #include "dwave-optimization/graph.hpp" +#include "dwave-optimization/nodes.hpp" + +namespace dwave::optimization::cp { + +// enum class CPConversionSatus { OK, Failed }; + +// CPConversionSatus parse_graph(const dwave::optimization::Graph& graph, CPModel& model) { +// CPConversionSatus status = CPConversionSatus::OK; + +// // We need to have a map from node to related CPVar. +// std::unordered_map node_map; + +// for (const auto& uptr : graph.nodes()) { +// const Node* ptr = uptr.get(); + +// if (auto iptr = dynamic_cast(ptr); iptr) { +// // Integer node +// CPVar* i = model.emplace_variable(model, ptr, ptr->topological_index()); +// node_map.insert({iptr, i}); + +// } else if (auto aptr = dynamic_cast(ptr); aptr) { +// // Add node +// CPVar* out = model.emplace_variable(model, ptr, ptr->topological_index()); +// node_map.insert({aptr, out}); + +// auto lhs_it = node_map.find(dynamic_cast(aptr->operands()[0])); +// auto rhs_it = node_map.find(dynamic_cast(aptr->operands()[1])); + +// if (lhs_it == node_map.end() or rhs_it == node_map.end()){ +// throw std::runtime_error("Accessing CP variables before their definition!"); +// } +// model.emplace_propagator(model.num_propagators(), lhs_it->second, +// rhs_it->second, out); +// } else if (auto aptr = dynamic_cast(ptr); aptr) { +// // Add node +// CPVar* out = model.emplace_variable(model, ptr, ptr->topological_index()); +// node_map.insert({aptr, out}); + +// auto lhs_it = node_map.find(dynamic_cast(aptr->operands()[0])); +// auto rhs_it = node_map.find(dynamic_cast(aptr->operands()[1])); + +// if (lhs_it == node_map.end() or rhs_it == node_map.end()){ +// throw std::runtime_error("Accessing CP variables before their definition!"); +// } +// model.emplace_propagator(model.num_propagators(), +// lhs_it->second, rhs_it->second, out); + +// } else { +// status = CPConversionSatus::Failed; +// break; +// } +// } + +// return status; +// } -namespace dwave::optimization::cp {} // namespace dwave::optimization::cp +} // namespace dwave::optimization::cp diff --git a/dwave/optimization/include/dwave-optimization/cp/core/cpvar.hpp b/dwave/optimization/include/dwave-optimization/cp/core/cpvar.hpp index bbd6c0d2c..115bc8bab 100644 --- a/dwave/optimization/include/dwave-optimization/cp/core/cpvar.hpp +++ b/dwave/optimization/include/dwave-optimization/cp/core/cpvar.hpp @@ -128,8 +128,8 @@ class CPVar { CPStatus remove_above(CPVarsState& state, double value, int index) const; CPStatus remove_below(CPVarsState& state, double value, int index) const; CPStatus assign(CPVarsState& state, double value, int index) const; - CPStatus set_min_size(CPVarsState& state, int index) const; - CPStatus set_max_size(CPVarsState& state, int index) const; + CPStatus set_min_size(CPVarsState& state, int new_min_size) const; + CPStatus set_max_size(CPVarsState& state, int new_max_size) const; // actions for the propagation engine void schedule_all(CPState& state, const std::vector& advisors, ssize_t index) const; @@ -176,7 +176,8 @@ class CPVar { void change_max(ssize_t i) override { var_->schedule_all(state_, var_->on_bounds, i); } void change_array_size(ssize_t i) override { - var_->schedule_all(state_, var_->on_array_size_change, i); + var_->schedule_all(state_, var_->on_bounds, i); + var_->schedule_all(state_, var_->on_domain, i); } private: diff --git a/dwave/optimization/include/dwave-optimization/cp/propagators/reduce.hpp b/dwave/optimization/include/dwave-optimization/cp/propagators/reduce.hpp index f97454e9a..130edefa9 100644 --- a/dwave/optimization/include/dwave-optimization/cp/propagators/reduce.hpp +++ b/dwave/optimization/include/dwave-optimization/cp/propagators/reduce.hpp @@ -32,4 +32,19 @@ class ReducePropagator : public Propagator { }; using SumPropagator = ReducePropagator>; + +template +class DynamicReducePropagator : public Propagator { + public: + DynamicReducePropagator(ssize_t index, CPVar* in, CPVar* out); + void initialize_state(CPState& state) const override; + CPStatus propagate(CPPropagatorsState& p_state, CPVarsState& v_state) const override; + + private: + // The variables entering the binary op + CPVar *in_, *out_; +}; + +using DynamicSumPropagator = DynamicReducePropagator>; + } // namespace dwave::optimization::cp diff --git a/dwave/optimization/src/cp/propagators/reduce.cpp b/dwave/optimization/src/cp/propagators/reduce.cpp index f6cb92817..c4595ffad 100644 --- a/dwave/optimization/src/cp/propagators/reduce.cpp +++ b/dwave/optimization/src/cp/propagators/reduce.cpp @@ -14,6 +14,8 @@ #include "dwave-optimization/cp/propagators/reduce.hpp" +#include + namespace dwave::optimization::cp { class SumPropagatorData : public PropagatorData { public: @@ -49,6 +51,8 @@ class SumPropagatorData : public PropagatorData { std::vector min, max; }; +// ----- ReducePropagator ----- + template ReducePropagator::ReducePropagator(ssize_t index, CPVar* in, CPVar* out) : Propagator(index), in_(in), out_(out) { @@ -81,7 +85,7 @@ CPStatus ReducePropagator>::propagate(CPPropagatorsState& p_st CPVarsState& v_state) const { auto data = data_ptr(p_state); auto n = in_->max_size(); - + CPStatus status = CPStatus::OK; std::deque& indices_to_process = data->indices_to_process(); @@ -144,4 +148,223 @@ CPStatus ReducePropagator>::propagate(CPPropagatorsState& p_st template class ReducePropagator>; +// ----- DynamicReducePropagator ----- + +template +DynamicReducePropagator::DynamicReducePropagator(ssize_t index, CPVar* in, CPVar* out) + : Propagator(index), in_(in), out_(out) { + // Default to bound consistency + Advisor in_adv(this, 0, std::make_unique()); + Advisor out_adv(this, 1, std::make_unique()); + + in_->propagate_on_bounds_change(std::move(in_adv)); + out_->propagate_on_bounds_change(std::move(out_adv)); +} + +template +void DynamicReducePropagator::initialize_state(CPState& state) const { + CPPropagatorsState& p_state = state.get_propagators_state(); + assert(propagator_index_ >= 0); + assert(propagator_index_ < static_cast(p_state.size())); + p_state[propagator_index_] = std::make_unique(state.get_state_manager(), 1); +} + +/// Compute the upper and lower bounds derived from the definitely present variables +std::pair compute_present_bounds(const CPVar* in, const CPVar* out, + const CPVarsState& state, ssize_t start_index, + double lb_start, double ub_start) { + // assert(start_index < in->min_size(state)); + + double lb_present = lb_start; + double ub_present = ub_start; + for (ssize_t j = start_index; j < in->min_size(state); ++j) { + assert(in->is_active(state, j)); + lb_present += in->min(state, j); + ub_present += in->max(state, j); + } + assert(lb_present <= ub_present); + return {lb_present, ub_present}; +} + +/// TODO: can make it more C++ :P +void compute_bounds_for_size(const CPVar* in, const CPVar* out, const CPVarsState& state, + std::vector& lb_optional, std::vector& ub_optional, + std::vector& lb_acc, std::vector& ub_acc, + double lb_present, double ub_present) { + // the next two variables used to accumulate the values + double acc_lb = lb_present; + double acc_ub = ub_present; + + lb_optional.push_back(lb_present); + ub_optional.push_back(ub_present); + lb_acc.push_back(lb_present); + ub_acc.push_back(ub_present); + + for (ssize_t j = in->min_size(state); j < in->max_size(state); ++j) { + assert(not in->is_active(state, j)); + assert(in->maybe_active(state, j)); + + lb_optional.push_back(acc_lb + std::min(0., in->min(state, j))); + ub_optional.push_back(acc_ub + std::max(0., in->max(state, j))); + + acc_lb += in->min(state, j); + acc_ub += in->max(state, j); + + // accumulate the lower + lb_acc.push_back(acc_lb); + ub_acc.push_back(acc_ub); + + assert(acc_lb <= acc_ub); + } + assert(ub_acc.size() == lb_acc.size()); + assert(ub_acc.size() == lb_optional.size()); + assert(ub_acc.size() == ub_optional.size()); +} + +template <> +CPStatus DynamicReducePropagator>::propagate(CPPropagatorsState& p_state, + CPVarsState& v_state) const { + auto data = data_ptr(p_state); + + CPStatus status = CPStatus::OK; + + std::deque& indices_to_process = data->indices_to_process(); + + while (data->num_indices_to_process() > 0) { + ssize_t i = indices_to_process.front(); + indices_to_process.pop_front(); + data->set_scheduled(false, i); + assert(i == 0); + + /// ==== Forward propagation ==== + + /// Compute the output interval induced by the present variables + double lb_present, ub_present; + std::tie(lb_present, ub_present) = compute_present_bounds(in_, out_, v_state, 0, 0, 0); + + // The following vectors store the lower and upper bounds calculations for each size of the + // vector. + // They are used for + // 1. Binding the out_ min and max (fwd) + // 2. Binding in_ min size and max size (bck). + // 3. Binding in_ min and max for every index (bck). + std::vector lb_optional, lb_acc; + std::vector ub_optional, ub_acc; + + // Fill the vectors above + compute_bounds_for_size(in_, out_, v_state, lb_optional, ub_optional, lb_acc, ub_acc, + lb_present, ub_present); + + // Prune out_ min value + status = out_->remove_below(v_state, + *std::min_element(lb_optional.begin(), lb_optional.end()), i); + + if (status == CPStatus::Inconsistency) return status; + + // Prune out_ max value + status = out_->remove_above(v_state, + *std::max_element(ub_optional.begin(), ub_optional.end()), i); + + if (status == CPStatus::Inconsistency) return status; + + /// ==== End of forward propagation ==== + + /// ==== Backward propagation ==== + + // Attempt to prune the size of in_. + // We start from min_size < max_size and we loop until we find a valid size. If no valid + // size found, we'll set some inconsistent values, and we return an inconsistency + ssize_t min_size_new = in_->max_size(v_state) + 1; + ssize_t max_size_new = in_->min_size(v_state) - 1; + ssize_t initial_min_size = in_->min_size(v_state); + + double out_min = out_->min(v_state, 0); + double out_max = out_->max(v_state, 0); + + for (ssize_t j = 0; j < static_cast(lb_optional.size()); ++j) { + // Update the min size if the domains overlap + if (not(ub_optional[j] < out_min or out_max < lb_optional[j])) { + min_size_new = std::min(in_->min_size(v_state) + j, min_size_new); + } + + // Update the max size if the domains overlap + if (not(ub_acc[j] < out_min or out_max < lb_acc[j])) { + max_size_new = std::max(in_->min_size(v_state) + j, max_size_new); + } + } + + assert(min_size_new >= in_->min_size(v_state)); + bool changed_min_size = (min_size_new > in_->min_size(v_state)) ? true : false; + + // After this we set the min and max size + status = in_->set_min_size(v_state, min_size_new); + if (status == CPStatus::Inconsistency) return status; + + status = in_->set_max_size(v_state, max_size_new); + if (status == CPStatus::Inconsistency) return status; + + // Now prune the values for each entry + // They might have changed, so need to recompute. + if (changed_min_size) { + // update the bounds from present variables + std::tie(lb_present, ub_present) = compute_present_bounds( + in_, out_, v_state, initial_min_size, lb_present, ub_present); + + // also update the bounds from optional variables, used to prune the + lb_optional.clear(); + ub_optional.clear(); + lb_acc.clear(); + ub_acc.clear(); + compute_bounds_for_size(in_, out_, v_state, lb_optional, ub_optional, lb_acc, ub_acc, + lb_present, ub_present); + } + + // Prune the present variables + // first compute the widest possible sum that there can be gven present and optional + // variables + { + double sum_max = *std::max_element(ub_optional.begin(), ub_optional.end()); + double sum_min = *std::min_element(lb_optional.begin(), lb_optional.end()); + for (ssize_t j = 0; j < in_->min_size(v_state); ++j) { + assert(in_->is_active(v_state, j)); + status = in_->remove_above( + v_state, out_->max(v_state, 0) - (sum_min - in_->min(v_state, j)), j); + if (status == CPStatus::Inconsistency) return status; + + status = in_->remove_below( + v_state, out_->min(v_state, 0) - (sum_max - in_->max(v_state, j)), j); + if (status == CPStatus::Inconsistency) return status; + } + } + + // Prune the optional variables, they receive the bounds from as if they were present and + // the following optional will contribute with either present or not, whichever yields + // bigger gap + { + for (ssize_t j = 0; j < in_->max_size(v_state) - in_->min_size(v_state); ++j) { + assert(not in_->is_active(v_state, in_->min_size(v_state) + j)); + assert(in_->maybe_active(v_state, in_->min_size(v_state) + j)); + + auto smin_it = std::min_element(lb_optional.begin() + j + 1, lb_optional.end()); + auto smax_it = std::max_element(ub_optional.begin() + j + 1, ub_optional.end()); + + double sum_min = (smin_it != lb_optional.end()) ? *smin_it : lb_optional.back(); + double sum_max = (smax_it != ub_optional.end()) ? *smax_it : ub_optional.back(); + + status = in_->remove_above( + v_state, out_->max(v_state, 0) - (sum_min - in_->min(v_state, j)), j); + if (status == CPStatus::Inconsistency) return status; + + status = in_->remove_below( + v_state, out_->min(v_state, 0) - (sum_max - in_->max(v_state, j)), j); + if (status == CPStatus::Inconsistency) return status; + } + } + } + + return status; +} + +template class DynamicReducePropagator>; + } // namespace dwave::optimization::cp diff --git a/tests/cpp/cp/propagators/test_reduce.cpp b/tests/cpp/cp/propagators/test_reduce.cpp index 72412f39c..4cddfd191 100644 --- a/tests/cpp/cp/propagators/test_reduce.cpp +++ b/tests/cpp/cp/propagators/test_reduce.cpp @@ -124,4 +124,129 @@ TEST_CASE("SumPropagator") { } } } + +TEST_CASE("DynamicSumPropagator") { + GIVEN("A Graph with a list node and its sum") { + dwave::optimization::Graph graph; + + // Add an integer node + ListNode* x = graph.emplace_node(5, 0, 5); + std::vector axes_s; + std::optional initial_s; + initial_s.emplace(0); + auto y = graph.emplace_node(x, axes_s, initial_s); + + // Lock the graph + graph.topological_sort(); + + // Construct the CP corresponding model + AND_GIVEN("A CP Model and fix-point engine") { + CPModel model; + CPEngine engine; + CPStatus status; + + // Add the variables to the model (manually) + CPVar* cp_x = model.emplace_variable(model, x, x->topological_index()); + CPVar* cp_y = model.emplace_variable(model, y, y->topological_index()); + + REQUIRE(cp_x->min_size() == 0); + REQUIRE(cp_x->max_size() == 5); + REQUIRE(cp_y->min_size() == 1); + REQUIRE(cp_y->max_size() == 1); + + // Add the propagator for the add node + Propagator* p_add = model.emplace_propagator( + model.num_propagators(), cp_x, cp_y); + + REQUIRE(cp_x->on_bounds.size() == 1); + REQUIRE(cp_y->on_bounds.size() == 1); + + WHEN("We initialize a state") { + CPState state = model.initialize_state(); + CPVarsState& v_state = state.get_variables_state(); + CPPropagatorsState& p_state = state.get_propagators_state(); + + REQUIRE(v_state.size() == 2); + REQUIRE(p_state.size() == 1); + + cp_x->initialize_state(state); + cp_y->initialize_state(state); + p_add->initialize_state(state); + + THEN("The interval of y is correctly set") { + REQUIRE(cp_y->min(v_state, 0) == 0); + REQUIRE(cp_y->max(v_state, 0) == 20); + REQUIRE(cp_y->size(v_state, 0) == 21); + } + + AND_WHEN("We remove x[3]") { + status = cp_x->remove_below(v_state, 6, 3); + REQUIRE(status == CPStatus::OK); + REQUIRE_FALSE(cp_x->is_active(v_state, 3)); + REQUIRE_FALSE(cp_x->maybe_active(v_state, 3)); + status = engine.fix_point(state); + + THEN("The sum output variable is correctly set to [0, 12]") { + REQUIRE(status == CPStatus::OK); + REQUIRE_FALSE(cp_x->is_active(v_state, 4)); + REQUIRE_FALSE(cp_x->maybe_active(v_state, 4)); + REQUIRE(cp_x->min_size(v_state) == 0); + REQUIRE(cp_x->max_size(v_state) == 3); + REQUIRE(cp_y->min(v_state, 0) == 0); + REQUIRE(cp_y->max(v_state, 0) == 12); + } + + AND_WHEN("We make x[0] in [2, 3]") { + status = cp_x->remove_above(v_state, 3, 0); + status = cp_x->remove_below(v_state, 2, 0); + REQUIRE(status == CPStatus::OK); + REQUIRE(cp_x->min(v_state, 0) == 2); + REQUIRE(cp_x->max(v_state, 0) == 3); + status = engine.fix_point(state); + THEN("The sum output is correcly set to [0, 11]") { + REQUIRE(status == CPStatus::OK); + REQUIRE(cp_y->min(v_state, 0) == 0); + REQUIRE(cp_y->max(v_state, 0) == 11); + } + + AND_WHEN("We prune the output to [5, 11]") { + status = cp_y->remove_below(v_state, 5, 0); + REQUIRE(status == CPStatus::OK); + REQUIRE(cp_y->min(v_state, 0) == 5); + status = engine.fix_point(state); + + THEN("The min size increases and the components are unchanged") { + REQUIRE(status == CPStatus::OK); + REQUIRE(cp_x->min_size(v_state) == 2); + REQUIRE(cp_x->max_size(v_state) == 3); + REQUIRE(cp_x->min(v_state, 0) == 2); + REQUIRE(cp_x->max(v_state, 0) == 3); + REQUIRE(cp_x->min(v_state, 1) == 0); + REQUIRE(cp_x->max(v_state, 1) == 4); + REQUIRE(cp_x->min(v_state, 2) == 0); + REQUIRE(cp_x->max(v_state, 2) == 4); + } + + AND_WHEN("We reduce the max size of the input to 2") { + status = cp_x->set_max_size(v_state, 2); + REQUIRE(status == CPStatus::OK); + REQUIRE(cp_x->max_size(v_state) == 2); + status = engine.fix_point(state); + THEN("The input and output are correctly set to x[1] in [2, 4] and " + "y in [5, 7]") { + REQUIRE(status == CPStatus::OK); + REQUIRE(cp_x->min(v_state, 1) == 2); + REQUIRE(cp_x->max(v_state, 1) == 4); + REQUIRE(cp_y->min(v_state, 0) == 5); + REQUIRE(cp_y->max(v_state, 0) == 7); + } + } + } + } + } + } + } + } +} + } // namespace dwave::optimization::cp \ No newline at end of file