Add dynamic reduce sum propagator

dwave/optimization/include/dwave-optimization/cp/core/base_node.hpp

@@ -14,12 +14,70 @@
 
 #pragma once
 
+#include <unordered_map>
+
 #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<const Node*, CPVar*> node_map;
+
+//     for (const auto& uptr : graph.nodes()) {
+//         const Node* ptr = uptr.get();
+
+//         if (auto iptr = dynamic_cast<const IntegerNode*>(ptr); iptr) {
+//             // Integer node
+//             CPVar* i = model.emplace_variable<CPVar>(model, ptr, ptr->topological_index());
+//             node_map.insert({iptr, i});
+
+//         } else if (auto aptr = dynamic_cast<const AddNode*>(ptr); aptr) {
+//             // Add node
+//             CPVar* out = model.emplace_variable<CPVar>(model, ptr, ptr->topological_index());
+//             node_map.insert({aptr, out});
+
+//             auto lhs_it = node_map.find(dynamic_cast<const ArrayNode*>(aptr->operands()[0]));
+//             auto rhs_it = node_map.find(dynamic_cast<const ArrayNode*>(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<AddPropagator>(model.num_propagators(), lhs_it->second,
+//             rhs_it->second, out);
+//         } else if (auto aptr = dynamic_cast<const LessEqualNode*>(ptr); aptr) {
+//             // Add node
+//             CPVar* out = model.emplace_variable<CPVar>(model, ptr, ptr->topological_index());
+//             node_map.insert({aptr, out});
+
+//             auto lhs_it = node_map.find(dynamic_cast<const ArrayNode*>(aptr->operands()[0]));
+//             auto rhs_it = node_map.find(dynamic_cast<const ArrayNode*>(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<LessEqualPropagator>(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

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<Advisor>& 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:

dwave/optimization/include/dwave-optimization/cp/propagators/reduce.hpp

@@ -32,4 +32,19 @@ class ReducePropagator : public Propagator {
 };
 
 using SumPropagator = ReducePropagator<std::plus<double>>;
+
+template <class BinaryOp>
+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<std::plus<double>>;
+
 }  // namespace dwave::optimization::cp

dwave/optimization/src/cp/propagators/reduce.cpp

@@ -14,6 +14,8 @@
 
 #include "dwave-optimization/cp/propagators/reduce.hpp"
 
+#include <algorithm>
+
 namespace dwave::optimization::cp {
 class SumPropagatorData : public PropagatorData {
  public:
@@ -49,6 +51,8 @@ class SumPropagatorData : public PropagatorData {
     std::vector<double> min, max;
 };
 
+// ----- ReducePropagator -----
+
 template <class BinaryOp>
 ReducePropagator<BinaryOp>::ReducePropagator(ssize_t index, CPVar* in, CPVar* out)
         : Propagator(index), in_(in), out_(out) {
@@ -81,7 +85,7 @@ CPStatus ReducePropagator<std::plus<double>>::propagate(CPPropagatorsState& p_st
                                                         CPVarsState& v_state) const {
     auto data = data_ptr<SumPropagatorData>(p_state);
     auto n = in_->max_size();
-    
+
     CPStatus status = CPStatus::OK;
 
     std::deque<ssize_t>& indices_to_process = data->indices_to_process();
@@ -144,4 +148,223 @@ CPStatus ReducePropagator<std::plus<double>>::propagate(CPPropagatorsState& p_st
 
 template class ReducePropagator<std::plus<double>>;
 
+// ----- DynamicReducePropagator -----
+
+template <class BinaryOp>
+DynamicReducePropagator<BinaryOp>::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<ReduceTransform>());
+    Advisor out_adv(this, 1, std::make_unique<ElementWiseTransform>());
+
+    in_->propagate_on_bounds_change(std::move(in_adv));
+    out_->propagate_on_bounds_change(std::move(out_adv));
+}
+
+template <class BinaryOp>
+void DynamicReducePropagator<BinaryOp>::initialize_state(CPState& state) const {
+    CPPropagatorsState& p_state = state.get_propagators_state();
+    assert(propagator_index_ >= 0);
+    assert(propagator_index_ < static_cast<ssize_t>(p_state.size()));
+    p_state[propagator_index_] = std::make_unique<PropagatorData>(state.get_state_manager(), 1);
+}
+
+/// Compute the upper and lower bounds derived from the definitely present variables
+std::pair<double, double> 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<double>& lb_optional, std::vector<double>& ub_optional,
+                             std::vector<double>& lb_acc, std::vector<double>& 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<std::plus<double>>::propagate(CPPropagatorsState& p_state,
+                                                               CPVarsState& v_state) const {
+    auto data = data_ptr<SumPropagatorData>(p_state);
+
+    CPStatus status = CPStatus::OK;
+
+    std::deque<ssize_t>& 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<double> lb_optional, lb_acc;
+        std::vector<double> 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<ssize_t>(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<std::plus<double>>;
+
 }  // namespace dwave::optimization::cp