Skip to content

Commit

Permalink
Implement self-loop removal
Browse files Browse the repository at this point in the history
  • Loading branch information
@arcondello
arcondello committed Oct 18, 2022
1 parent 712c432 commit 9cdbdae
Show file tree
Hide file tree
Showing 4 changed files with 142 additions and 77 deletions.
21 changes: 21 additions & 0 deletions dimod/include/dimod/abc.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -248,6 +248,9 @@ class QuadraticModelBase {
template <class T>
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 <class B, class I>
bool is_equal(const QuadraticModelBase<B, I>& other) const;
Expand Down Expand Up @@ -696,6 +699,24 @@ void QuadraticModelBase<bias_type, index_type>::fix_variable(index_type v, T ass
remove_variable(v);
}

template <class bias_type, class index_type>
bool QuadraticModelBase<bias_type, index_type>::has_interaction(index_type u, index_type v) const {
assert(0 <= u && static_cast<size_type>(u) < num_variables());
assert(0 <= v && static_cast<size_type>(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 <class bias_type, class index_type>
template <class B, class I>
bool QuadraticModelBase<bias_type, index_type>::is_equal(
Expand Down
26 changes: 26 additions & 0 deletions dimod/include/dimod/expression.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -209,6 +209,9 @@ class Expression : public abc::QuadraticModelBase<Bias, Index> {
template <class T>
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;
Expand Down Expand Up @@ -410,6 +413,19 @@ void Expression<bias_type, index_type>::fix_variable(index_type v, T assignment)
throw std::logic_error("not implemented - fix_variable");
}

template <class bias_type, class index_type>
bool Expression<bias_type, index_type>::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<size_type>(u) < parent_->num_variables());
assert(v >= 0 && static_cast<size_type>(v) < parent_->num_variables());
return 0;
}

return base_type::has_interaction(uit->second, vit->second);
}

template <class bias_type, class index_type>
bool Expression<bias_type, index_type>::has_variable(index_type v) const {
return indices_.count(v);
Expand Down Expand Up @@ -543,6 +559,16 @@ void Expression<bias_type, index_type>::reindex_variables(index_type v) {
assert(indices_.size() == variables_.size());
}

template <class bias_type, class index_type>
bool Expression<bias_type, index_type>::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 <class bias_type, class index_type>
void Expression<bias_type, index_type>::remove_variable(index_type v) {
throw std::logic_error("not implemented - remove_variable");
Expand Down
67 changes: 59 additions & 8 deletions dimod/include/dimod/presolve.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
#pragma once

#include <algorithm>
#include <unordered_map>
#include <utility>
#include <vector>

Expand All @@ -31,6 +32,8 @@ class PostSolver {
using size_type = std::size_t;
using assignment_type = Assignment;

void add_variable(index_type v);

template <class T>
std::vector<T> apply(std::vector<T> reduced) const;

Expand All @@ -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
Expand All @@ -58,6 +61,12 @@ class PostSolver {
std::vector<Transform> transforms_;
};

template <class bias_type, class index_type, class assignment_type>
void PostSolver<bias_type, index_type, assignment_type>::add_variable(index_type v) {
transforms_.emplace_back(TransformKind::ADD);
transforms_.back().v = v;
}

template <class bias_type, class index_type, class assignment_type>
template <class T>
std::vector<T> PostSolver<bias_type, index_type, assignment_type>::apply(
Expand All @@ -72,6 +81,9 @@ std::vector<T> PostSolver<bias_type, index_type, assignment_type>::apply(
sample[it->v] *= it->multiplier;
sample[it->v] += it->offset;
break;
case TransformKind::ADD:
sample.erase(sample.begin() + it->v);
break;
}
}
return sample;
Expand All @@ -87,8 +99,8 @@ void PostSolver<bias_type, index_type, assignment_type>::fix_variable(index_type

template <class bias_type, class index_type, class assignment_type>
void PostSolver<bias_type, index_type, assignment_type>::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;
Expand Down Expand Up @@ -120,11 +132,46 @@ class PreSolver {
private:
model_type model_;
PostSolver<bias_type, index_type, assignment_type> postsolver_;
};

void substitute_self_loops_expr(Expression<bias_type, index_type>& expression,
std::unordered_map<index_type, index_type>& 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<size_type>(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<index_type, index_type> 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 <class bias_type, class index_type, class assignment_type>
PreSolver<bias_type, index_type, assignment_type>::PreSolver(): model_(), postsolver_() {}
PreSolver<bias_type, index_type, assignment_type>::PreSolver() : model_(), postsolver_() {}

template <class bias_type, class index_type, class assignment_type>
PreSolver<bias_type, index_type, assignment_type>::PreSolver(model_type model)
Expand Down Expand Up @@ -161,6 +208,9 @@ void PreSolver<bias_type, index_type, assignment_type>::apply() {
}
}

// *-- remove self-loops
substitute_self_loops();

// Trivial techniques -----------------------------------------------------

bool changes = true;
Expand Down Expand Up @@ -194,8 +244,8 @@ void PreSolver<bias_type, index_type, assignment_type>::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 {
Expand Down Expand Up @@ -253,7 +303,8 @@ void PreSolver<bias_type, index_type, assignment_type>::load_default_presolvers(
}

template <class bias_type, class index_type, class assignment_type>
const ConstrainedQuadraticModel<bias_type, index_type>& PreSolver<bias_type, index_type, assignment_type>::model() const {
const ConstrainedQuadraticModel<bias_type, index_type>&
PreSolver<bias_type, index_type, assignment_type>::model() const {
return model_;
}

Expand Down
105 changes: 36 additions & 69 deletions testscpp/tests/test_presolve.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,6 @@
#include "dimod/constrained_quadratic_model.h"
#include "dimod/presolve.h"


namespace dimod {

SCENARIO("constrained quadratic models can be presolved") {
Expand Down Expand Up @@ -127,74 +126,42 @@ SCENARIO("constrained quadratic models can be presolved") {
}
}

// WHEN("passed through the identity presolver") {
// auto presolver = presolve::PreSolver<double>(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<double, Vartype::INTEGER>::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<double, Vartype::INTEGER>::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<double, Vartype::INTEGER>::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<double>(cqm);
// presolver.add_presolver<presolve::techniques::TrivialPresolver<double>>();
// 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<int>{3};

// auto original = postsolver.apply(reduced);

// CHECK(original == std::vector<int>{3, 5, 7, 5});
// }
// }
// }
}
GIVEN("a CQM with self-loops") {
auto cqm = ConstrainedQuadraticModel<double>();
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<double>(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<int>{1, 2, 3, 4, 5, 6, 7, 8});
CHECK(original == std::vector<int>{1, 2, 3, 4, 5});
}
}
}
}

} // namespace dimod

0 comments on commit 9cdbdae

Please sign in to comment.