diff --git a/dimod/include/dimod/abc.h b/dimod/include/dimod/abc.h index 93e30a8af..3ea79e5cd 100644 --- a/dimod/include/dimod/abc.h +++ b/dimod/include/dimod/abc.h @@ -248,6 +248,9 @@ class QuadraticModelBase { template void fix_variable(index_type v, T assignment); + /// Check whether u and v have an interaction + bool has_interaction(index_type u, index_type v) const; + /// Test whether two quadratic models are equal. template bool is_equal(const QuadraticModelBase& other) const; @@ -696,6 +699,24 @@ void QuadraticModelBase::fix_variable(index_type v, T ass remove_variable(v); } +template +bool QuadraticModelBase::has_interaction(index_type u, index_type v) const { + assert(0 <= u && static_cast(u) < num_variables()); + assert(0 <= v && static_cast(v) < num_variables()); + + if (!adj_ptr_) { + return false; + } + + const auto& n = (*adj_ptr_)[u]; + auto it = std::lower_bound(n.cbegin(), n.cend(), v); + if (it == n.cend() || it->v != v) { + return false; + } + + return true; +} + template template bool QuadraticModelBase::is_equal( diff --git a/dimod/include/dimod/expression.h b/dimod/include/dimod/expression.h index db1896de0..0ed6f0ae2 100644 --- a/dimod/include/dimod/expression.h +++ b/dimod/include/dimod/expression.h @@ -209,6 +209,9 @@ class Expression : public abc::QuadraticModelBase { template void fix_variable(index_type v, T assignment); + /// Check whether u and v have an interaction + bool has_interaction(index_type u, index_type v) const; + bool has_variable(index_type v) const; bool is_disjoint(const Expression& other) const; @@ -410,6 +413,19 @@ void Expression::fix_variable(index_type v, T assignment) throw std::logic_error("not implemented - fix_variable"); } +template +bool Expression::has_interaction(index_type u, index_type v) const { + auto uit = indices_.find(u); + auto vit = indices_.find(v); + if (uit == indices_.end() || vit == indices_.end()) { + assert(u >= 0 && static_cast(u) < parent_->num_variables()); + assert(v >= 0 && static_cast(v) < parent_->num_variables()); + return 0; + } + + return base_type::has_interaction(uit->second, vit->second); +} + template bool Expression::has_variable(index_type v) const { return indices_.count(v); @@ -543,6 +559,16 @@ void Expression::reindex_variables(index_type v) { assert(indices_.size() == variables_.size()); } +template +bool Expression::remove_interaction(index_type u, index_type v) { + auto uit = indices_.find(u); + auto vit = indices_.find(v); + if (uit == indices_.end() || vit == indices_.end()) { + return false; + } + return base_type::remove_interaction(uit->second, vit->second); +} + template void Expression::remove_variable(index_type v) { throw std::logic_error("not implemented - remove_variable"); diff --git a/dimod/include/dimod/presolve.h b/dimod/include/dimod/presolve.h index 3996d16ff..c49c14339 100644 --- a/dimod/include/dimod/presolve.h +++ b/dimod/include/dimod/presolve.h @@ -15,6 +15,7 @@ #pragma once #include +#include #include #include @@ -31,6 +32,8 @@ class PostSolver { using size_type = std::size_t; using assignment_type = Assignment; + void add_variable(index_type v); + template std::vector apply(std::vector reduced) const; @@ -40,7 +43,7 @@ class PostSolver { private: // we want to track what changes were made - enum TransformKind { FIX, SUBSTITUTE }; + enum TransformKind { FIX, SUBSTITUTE, ADD }; // todo: we could get fancy with pointers and templates to save a bit of // space here @@ -58,6 +61,12 @@ class PostSolver { std::vector transforms_; }; +template +void PostSolver::add_variable(index_type v) { + transforms_.emplace_back(TransformKind::ADD); + transforms_.back().v = v; +} + template template std::vector PostSolver::apply( @@ -72,6 +81,9 @@ std::vector PostSolver::apply( sample[it->v] *= it->multiplier; sample[it->v] += it->offset; break; + case TransformKind::ADD: + sample.erase(sample.begin() + it->v); + break; } } return sample; @@ -87,8 +99,8 @@ void PostSolver::fix_variable(index_type template void PostSolver::substitute_variable(index_type v, - bias_type multiplier, - bias_type offset) { + bias_type multiplier, + bias_type offset) { assert(multiplier); // cannot undo when it's 0 transforms_.emplace_back(TransformKind::SUBSTITUTE); transforms_.back().v = v; @@ -120,11 +132,46 @@ class PreSolver { private: model_type model_; PostSolver postsolver_; -}; + void substitute_self_loops_expr(Expression& expression, + std::unordered_map& mapping) { + size_type num_variables = expression.num_variables(); + for (size_type i = 0; i < num_variables; ++i) { + index_type v = expression.variables()[i]; + + if (!expression.has_interaction(v, v)) continue; // no self loop + + auto out = mapping.emplace(v, model_.num_variables()); + + if (out.second) { + // we haven't seen this variable before + model_.add_variable(model_.vartype(v), model_.lower_bound(v), + model_.upper_bound(v)); + postsolver_.add_variable(out.first->second); + } + + assert(static_cast(out.first->second) < model_.num_variables()); + + // now set the bias between v and the new variable + expression.add_quadratic(v, out.first->second, expression.quadratic(v, v)); + expression.remove_interaction(v, v); + } + } + + // todo: break into separate presolver + void substitute_self_loops() { + std::unordered_map mapping; + + substitute_self_loops_expr(model_.objective, mapping); + + for (size_type c = 0; c < model_.num_constraints(); ++c) { + substitute_self_loops_expr(model_.constraint_ref(c), mapping); + } + } +}; template -PreSolver::PreSolver(): model_(), postsolver_() {} +PreSolver::PreSolver() : model_(), postsolver_() {} template PreSolver::PreSolver(model_type model) @@ -161,6 +208,9 @@ void PreSolver::apply() { } } + // *-- remove self-loops + substitute_self_loops(); + // Trivial techniques ----------------------------------------------------- bool changes = true; @@ -194,8 +244,8 @@ void PreSolver::apply() { // todo: test if negative if (constraint.sense() == Sense::EQ) { - model_.set_lower_bound(v, std::max(rhs, model_.lower_bound(v))); - model_.set_upper_bound(v, std::min(rhs, model_.upper_bound(v))); + model_.set_lower_bound(v, std::max(rhs, model_.lower_bound(v))); + model_.set_upper_bound(v, std::min(rhs, model_.upper_bound(v))); } else if ((constraint.sense() == Sense::LE) != (a < 0)) { model_.set_upper_bound(v, std::min(rhs, model_.upper_bound(v))); } else { @@ -253,7 +303,8 @@ void PreSolver::load_default_presolvers( } template -const ConstrainedQuadraticModel& PreSolver::model() const { +const ConstrainedQuadraticModel& +PreSolver::model() const { return model_; } diff --git a/testscpp/tests/test_presolve.cpp b/testscpp/tests/test_presolve.cpp index e81f583fb..e08a4adae 100644 --- a/testscpp/tests/test_presolve.cpp +++ b/testscpp/tests/test_presolve.cpp @@ -16,7 +16,6 @@ #include "dimod/constrained_quadratic_model.h" #include "dimod/presolve.h" - namespace dimod { SCENARIO("constrained quadratic models can be presolved") { @@ -127,74 +126,42 @@ SCENARIO("constrained quadratic models can be presolved") { } } - // WHEN("passed through the identity presolver") { - // auto presolver = presolve::PreSolver(cqm); - // presolver.apply(); - - // THEN("nothing has changed") { - // // CHECK(presolver.result.status == unchanged) - - // const auto& newcqm = presolver.model(); - - // REQUIRE(newcqm.num_variables() == 4); - // REQUIRE(newcqm.num_constraints() == 4); - - // CHECK(newcqm.vartype(v0) == Vartype::INTEGER); - // CHECK(newcqm.lower_bound(v0) == 0); - // CHECK(newcqm.upper_bound(v0) == - // vartype_limits::default_max()); - // CHECK(newcqm.constraint_ref(c0).linear(v0) == 1); - // CHECK(newcqm.constraint_ref(c0).sense() == Sense::LE); - // CHECK(newcqm.constraint_ref(c0).rhs() == 5); - - // CHECK(newcqm.vartype(v1) == Vartype::INTEGER); - // CHECK(newcqm.lower_bound(v1) == 5); - // CHECK(newcqm.upper_bound(v1) == 5); - - // CHECK(newcqm.vartype(v2) == Vartype::INTEGER); - // CHECK(newcqm.lower_bound(v2) == 0); - // CHECK(newcqm.upper_bound(v2) == - // vartype_limits::default_max()); - // CHECK(newcqm.constraint_ref(c2).linear(v2) == 1); - // CHECK(newcqm.constraint_ref(c2).sense() == Sense::EQ); - // CHECK(newcqm.constraint_ref(c2).rhs() == 7); - - // CHECK(newcqm.vartype(v3) == Vartype::INTEGER); - // CHECK(newcqm.lower_bound(v3) == 0); - // CHECK(newcqm.upper_bound(v3) == - // vartype_limits::default_max()); - // CHECK(newcqm.constraint_ref(c3a).linear(v3) == 1); - // CHECK(newcqm.constraint_ref(c3a).sense() == Sense::LE); - // CHECK(newcqm.constraint_ref(c3a).rhs() == 5.5); - // CHECK(newcqm.constraint_ref(c3b).linear(v3) == 1); - // CHECK(newcqm.constraint_ref(c3b).sense() == Sense::GE); - // CHECK(newcqm.constraint_ref(c3b).rhs() == 4.5); - // } - // } - - // WHEN("passed through the trivial presolver") { - // auto presolver = presolve::PreSolver(cqm); - // presolver.add_presolver>(); - // presolver.apply(); - - // THEN("several constraints/variables are removed") { - // const auto& newcqm = presolver.model(); - // const auto& postsolver = presolver.postsolver(); - - // CHECK(newcqm.num_constraints() == 0); - - // CHECK(newcqm.num_variables() == 1); - - // // see if we restore the original problem - // auto reduced = std::vector{3}; - - // auto original = postsolver.apply(reduced); - - // CHECK(original == std::vector{3, 5, 7, 5}); - // } - // } - // } -} + GIVEN("a CQM with self-loops") { + auto cqm = ConstrainedQuadraticModel(); + cqm.add_variables(Vartype::INTEGER, 5); + + cqm.objective.set_quadratic(0, 0, 1.5); + cqm.objective.set_quadratic(3, 3, 3.5); + + cqm.add_constraint(); + cqm.constraint_ref(0).add_quadratic(4, 4, 5); + cqm.constraint_ref(0).add_quadratic(3, 3, 6); + + WHEN("presolving is applied") { + auto presolver = presolve::PreSolver(std::move(cqm)); + presolver.load_default_presolvers(); + presolver.apply(); + THEN("the self-loops are removed") { + auto& model = presolver.model(); + + REQUIRE(model.num_variables() == 8); + REQUIRE(model.num_constraints() == 1); + CHECK(model.objective.num_interactions() == 2); + CHECK(model.objective.quadratic(0, 5) == 1.5); + CHECK(model.objective.quadratic(3, 6) == 3.5); + + CHECK(model.constraint_ref(0).num_interactions() == 2); + CHECK(model.constraint_ref(0).quadratic(4, 7) == 5); + CHECK(model.constraint_ref(0).quadratic(3, 6) == 6); + } + + AND_WHEN("we then undo the transformation") { + auto original = presolver.postsolver().apply(std::vector{1, 2, 3, 4, 5, 6, 7, 8}); + CHECK(original == std::vector{1, 2, 3, 4, 5}); + } + } + } +} } // namespace dimod