Fix Expression::fix_variables() for quadratic expressions

Closes dwavesystems#1351

dimod/include/dimod/constrained_quadratic_model.h

@@ -545,9 +545,8 @@ template <class bias_type, class index_type>
 void ConstrainedQuadraticModel<bias_type, index_type>::fix_variables_expr(
         const Expression<bias_type, index_type>& src, Expression<bias_type, index_type>& dst,
         const std::vector<index_type>& old_to_new, const std::vector<bias_type>& assignments) {
-    // We'll want to access the expressions by index for speed
+    // We'll want to access the source expression by index for speed
     const abc::QuadraticModelBase<bias_type, index_type>& isrc = src;
-    abc::QuadraticModelBase<bias_type, index_type>& idst = dst;
 
     // offset
     dst.add_offset(src.offset());
@@ -565,23 +564,27 @@ void ConstrainedQuadraticModel<bias_type, index_type>::fix_variables_expr(
         }
     }
 
-    // quadratic, and it's safe to do everything by index!
-    for (auto it = isrc.cbegin_quadratic(); it != isrc.cend_quadratic(); ++it) {
-        auto u = src.variables()[it->u];
-        auto v = src.variables()[it->v];
+    // quadratic
+    for (auto it = isrc.cbegin_quadratic(), end = isrc.cend_quadratic(); it != end; ++it) {
+        const index_type u = src.variables()[it->u];  // variable u in the source
+        const index_type v = src.variables()[it->v];  // variable v in the source
+        const bias_type bias = it->bias;  // bias in the source
 
-        if (old_to_new[u] < 0 && old_to_new[v] < 0) {
+        const index_type new_u = old_to_new[u];  // variable u in the destination
+        const index_type new_v = old_to_new[v];  // variable v in the destination
+
+        if (new_u < 0 && new_v < 0) {
             // both fixed, becomes offset
-            idst.add_offset(assignments[u] * assignments[v] * it->bias);
-        } else if (old_to_new[u] < 0) {
+            dst.add_offset(assignments[u] * assignments[v] * bias);
+        } else if (new_u < 0) {
             // u fixed, v unfixed
-            idst.add_linear(old_to_new[it->v], assignments[u] * it->bias);
-        } else if (old_to_new[v] < 0) {
+            dst.add_linear(new_v, assignments[u] * bias);
+        } else if (new_v < 0) {
             // u unfixed, v fixed
-            idst.add_linear(old_to_new[it->u], assignments[v] * it->bias);
+            dst.add_linear(new_u, assignments[v] * bias);
         } else {
             // neither fixed
-            idst.add_quadratic_back(old_to_new[it->u], old_to_new[it->v], it->bias);
+            dst.add_quadratic_back(new_u, new_v, bias);
         }
     }
 }

dimod/include/dimod/expression.h

@@ -366,7 +366,7 @@ void Expression<bias_type, index_type>::add_quadratic(ItRow row_iterator, ItCol
 
 template <class bias_type, class index_type>
 void Expression<bias_type, index_type>::add_quadratic_back(index_type u, index_type v, bias_type bias) {
-    throw std::logic_error("not implemented - add_quadratic_back");
+    base_type::add_quadratic_back(enforce_variable(u), enforce_variable(v), bias);
 }
 
 template <class bias_type, class index_type>

releasenotes/notes/fix-quadratic-cqm-fix-variables-43f21519357c228f.yaml

@@ -0,0 +1,7 @@
+---
+fixes:
+  - |
+    Fix ``dimod::ConstrainedQuadraticModel::fix_variables()`` C++ method to work
+    correctly with quadratic objectives and constraints.
+    This fixes a bug introduced in dimod 0.12.5.
+    See `#1351 <https://github.com/dwavesystems/dimod/issues/1351>`_.

testscpp/tests/test_constrained_quadratic_model.cpp

@@ -871,6 +871,40 @@ TEST_CASE("Test CQM.constraint_weak_ptr()") {
     }
 }
 
+TEST_CASE("Test ConstrainedQuadraticModel::fix_variables()") {
+    GIVEN("A CQM with a quadratic constraint") {
+        auto cqm = ConstrainedQuadraticModel<double>();
+
+        cqm.add_variables(Vartype::BINARY, 5);
+
+        auto c = cqm.add_linear_constraint({3, 1, 4, 0}, {1, 2, 3, 4}, Sense::LE, 5);
+        cqm.constraint_ref(c).add_quadratic(3, 1, -15);
+        cqm.constraint_ref(c).add_quadratic(4, 0, 14);
+        cqm.constraint_ref(c).add_quadratic(4, 3, 17);
+
+        WHEN("we fix some variables") {
+            std::vector<int> variables{1, 2, 3};
+            std::vector<double> values{1, 1, 1};
+
+            auto fixed = cqm.fix_variables(variables.begin(), variables.end(), values.begin());
+            // relabels x0 -> x0
+            // relabels x4 -> x1
+
+            THEN("the offset, linear, and quadratic biases are updated as expected") {
+                REQUIRE(fixed.num_constraints() == 1);
+
+                CHECK(fixed.constraint_ref(c).num_variables() == 2);
+                CHECK(fixed.constraint_ref(c).num_interactions() == 1);
+
+                CHECK(fixed.constraint_ref(c).offset() == -12);
+                CHECK(fixed.constraint_ref(c).linear(1) == 20);
+                CHECK(fixed.constraint_ref(c).linear(0) == 4);
+                CHECK(fixed.constraint_ref(c).quadratic(1, 0) == 14);
+            }
+        }
+    }
+}
+
 TEST_CASE("Test Expression::add_quadratic()") {
     GIVEN("A CQM with two variables with vartypes") {
         auto cqm = dimod::ConstrainedQuadraticModel<double>();
@@ -901,6 +935,33 @@ TEST_CASE("Test Expression::add_quadratic()") {
     }
 }
 
+TEST_CASE("Test Expression::add_quadratic_back()") {
+    GIVEN("A CQM with three variables and one constraints") {
+        auto cqm = dimod::ConstrainedQuadraticModel<double>();
+        auto i = cqm.add_variable(Vartype::INTEGER);
+        auto x = cqm.add_variable(Vartype::BINARY);
+        auto y = cqm.add_variable(Vartype::BINARY);
+
+        auto c0 = cqm.add_linear_constraint({i, x, y}, {1, 1, 1}, Sense::LE, 1);
+
+        WHEN("add_quadratic_back() is used on an empty constraint") {
+            cqm.constraint_ref(c0).add_quadratic_back(i, x, 1.5);
+
+            THEN("we can read the value out as expected") {
+                CHECK(cqm.constraint_ref(c0).quadratic(i, x) == 1.5);
+            }
+
+            AND_WHEN("we add another quadratic") {
+                cqm.constraint_ref(c0).add_quadratic_back(y, x, 2.5);
+
+                THEN("we can read the value out as expected") {
+                    CHECK(cqm.constraint_ref(c0).quadratic(y, x) == 2.5);
+                }
+            }
+        }
+    }
+}
+
 TEST_CASE("Test Expression::remove_variables()") {
     GIVEN("A CQM with five variables, an objective, and one constraint") {
         auto cqm = ConstrainedQuadraticModel<double>();