Merge branch 'latest' of https://github.com/fwesselm/HiGHS into aggregatedImplications

Author: fwesselm

check/TestPresolve.cpp

@@ -1276,14 +1276,100 @@ TEST_CASE("test-non-stop-initial-sweep", "[highs_test_presolve]") {
   h.resetGlobalScheduler(true);
 }
 
+TEST_CASE("test-duplicate-row", "[highs_test_presolve]") {
+  HighsLp lp;
+  lp.num_col_ = 3;
+  lp.num_row_ = 2;
+  lp.col_cost_ = {2, -1, 1};
+  lp.col_lower_ = {0, 0, 0};
+  lp.col_upper_ = {kHighsInf, 2, kHighsInf};
+  lp.a_matrix_.start_ = {0, 2, 4, 6};
+  lp.a_matrix_.index_ = {0, 1, 0, 1, 0, 1};
+  Highs h;
+  h.setOptionValue("output_flag", dev_run);
+  h.setOptionValue("presolve_rule_logging", true);
+  h.setOptionValue(kSolverString, kIpxString);
+
+  lp.row_lower_ = {-kHighsInf, -4};
+  lp.row_upper_ = {0, 0};
+  for (HighsInt bound_flip = 0; bound_flip < 2; bound_flip++) {
+    lp.a_matrix_.value_ = {1, -4, -1, 4, 3, -12};
+    for (HighsInt lhs_flip = 0; lhs_flip < 2; lhs_flip++) {
+      // Base model has parallel rows
+      //
+      // r0:         x -  y +  3z <= 0
+      // r1: -4 <= -4x + 4y - 12z <= 0
+      //
+      // After x is fixed at 0 (dominated column) these parallel rows
+      // are deduced to be the doubleton equation
+      //
+      // -y + 3z = 0
+      //
+      // with r1 being removed - presumably because the elimination
+      // multiplier (-0.25) is then less than 1 in magnitude.
+      //
+      // This allows z to be substituted and y is then fixed at 2,
+      // reducing the problem to empty
+      //
+      // In substitution, the row dual for the equation is 1/3
+      //
+      // * r0 is made basic because its lower bound was tightened, and
+      // the sign of the dual means that it can't be nonbasic so it
+      // can't be nonbasic
+      //
+      // * r1 is made nonbasic and the dual is now correctly scaled by
+      // _multiplying_ by (-0.25) - since the row values are larger,
+      // the dual must be reduced in magnitude - and should be viewed
+      // as being at its upper bound
+      //
+      // However, in DuplicateRow::undo, computeRowDualAndStatus was
+      // previously only passed "tightened", and no indication of
+      // whether that was at its lower or upper bound. Hence it made
+      // r1 nonbasic at its lower bound when it tightens the lower
+      // bound of r0, due to the negative sign of the scale factor.
+      //
+      // However, it's the upper bound on r1 that tightens the lower
+      // bound on r0. By passing the direction sign -1 (+1) if r0 is
+      // tightened at its lower (upper) bound, r1 is now set to be
+      // nonbasic at the correct bound.
+      //
+      // The spurious lower bound on r1 is necessary to expose the
+      // consequences of the error since the simplex solver only needs
+      // the basic/nonbasic status to set values of variables to
+      // bound, unless they are ranged, in which case the
+      // HighsBasisStatus being lower or upper is used.
+      //
+      // The multiple passes ensure code coverage in
+      // DuplicateRow::undo - all four cases are passed to
+      // computeRowDualAndStatus - and test the correctness of both
+      // primal-dual and basis postsolve.
+      //
+      h.setOptionValue("run_crossover", kHighsOnString);
+      for (HighsInt k = 0; k < 2; k++) {
+        REQUIRE(h.passModel(lp) == HighsStatus::kOk);
+        REQUIRE(h.run() == HighsStatus::kOk);
+
+        REQUIRE(h.getModelStatus() == HighsModelStatus::kOptimal);
+        h.setOptionValue("run_crossover", kHighsOffString);
+      }
+      lp.a_matrix_.value_ = {-1, 4, 1, -4, -3, 12};
+    }
+    lp.row_lower_ = {0, 0};
+    lp.row_upper_ = {kHighsInf, 4};
+  }
+
+  h.resetGlobalScheduler(true);
+}
+
 /*
 TEST_CASE("test-fuzzing", "[highs_test_presolve]") {
   Highs h;
   //  h.setOptionValue("output_flag", dev_run);
   //  if (dev_run) {
   printf("\n====================\nWithout presolve\n====================\n");
 
-  const std::string model = "issue-009";
+  const std::string model = "issue-002";
+  const bool reduces_to_empty = true;
   std::string model_file = std::string(HIGHS_DIR) + "/build/OscarFuzzing/" +
                            model + "/" + model + ".mps";
 
@@ -1303,23 +1389,23 @@ TEST_CASE("test-fuzzing", "[highs_test_presolve]") {
   std::string options_file =
       std::string(HIGHS_DIR) + "/build/OscarFuzzing/" + model + "/options.txt";
   REQUIRE(h.readOptions(options_file) == HighsStatus::kOk);
-  //  REQUIRE(h.setOptionValue("presolve_rule_off", 1 <<
-kPresolveRuleColStuffing) == HighsStatus::kOk); HighsOptions options =
-h.getOptions();
+  HighsOptions options = h.getOptions();
 
-  printf("\n====================\nPresolved LP\n====================\n");
-  h.presolve();
+  if (!reduces_to_empty) {
+    printf("\n====================\nPresolved LP\n====================\n");
+    h.presolve();
 
-  HighsLp lp = h.getPresolvedLp();
+    HighsLp lp = h.getPresolvedLp();
 
-  h.clear();
-  h.passModel(lp);
+    h.clear();
+    h.passModel(lp);
 
-  h.setOptionValue(kPresolveString, kHighsOffString);
+    h.setOptionValue(kPresolveString, kHighsOffString);
 
-  h.run();
-  h.writeSolution("", 1);
-  h.clear();
+    h.run();
+    h.writeSolution("", 1);
+    h.clear();
+  }
 
   printf(
       "\n====================\nPresolve no crossover\n====================\n");
@@ -1330,6 +1416,8 @@ h.getOptions();
   h.setOptionValue("log_dev_level", 1);
   h.writeOptions("", true);
 
+  //  h.setOptionValue(kSolverString, kIpxString);
+
   h.run();
   h.writeSolution("", 1);
 

highs/lp_data/HighsOptions.h

@@ -1500,7 +1500,7 @@ class HighsOptions : public HighsOptionsStruct {
 
     record_bool = new OptionRecordBool(
         "use_implied_bounds_from_presolve",
-        "Use relaxed implied bounds from presolve", advanced,
+        "Use relaxed implied bounds from presolve: redundant option!", advanced,
         &use_implied_bounds_from_presolve, false);
     records.push_back(record_bool);
 

highs/mip/feasibilityjump.hh

@@ -145,7 +145,7 @@ class Problem {
     // If we are relaxing continuous variables, an equality needs to be split
     // into Gte and Lte.
     if (relax_continuous > 0 && sense == RowType::Equal)
-      if (std::any_of(rowVarIdxs, rowVarIdxs + numCoeffs, [&](double varIdx) {
+      if (std::any_of(rowVarIdxs, rowVarIdxs + numCoeffs, [&](int varIdx) {
             return vars[varIdx].vartype == VarType::Continuous;
           })) {
         addConstraint(RowType::Gte, rhs, numCoeffs, rowVarIdxs, rowCoeffs,
@@ -350,7 +350,6 @@ class JumpMove {
       else {
         constraintBounds.emplace_back(-std::numeric_limits<double>::infinity(),
                                       constraint.rhs);
-        constraintBounds.emplace_back(constraint.rhs, constraint.rhs);
         constraintBounds.emplace_back(constraint.rhs,
                                       std::numeric_limits<double>::infinity());
       }
@@ -359,9 +358,17 @@ class JumpMove {
         double residualIncumbent =
             constraint.incumbentLhs - cell.coeff * varIncumbentValue;
 
+        double boundForValidLower, boundForValidUpper;
+        if (cell.coeff < 0.) {
+          boundForValidLower = bound.second;
+          boundForValidUpper = bound.first;
+        } else {
+          boundForValidLower = bound.first;
+          boundForValidUpper = bound.second;
+        }
         std::pair<double, double> validRange = {
-            ((1.0 / cell.coeff) * (bound.first - residualIncumbent)),
-            ((1.0 / cell.coeff) * (bound.second - residualIncumbent)),
+            ((1.0 / cell.coeff) * (boundForValidLower - residualIncumbent)),
+            ((1.0 / cell.coeff) * (boundForValidUpper - residualIncumbent)),
         };
 
         if (problem.vars[varIdx].vartype == VarType::Integer)
@@ -370,8 +377,6 @@ class JumpMove {
               std::floor(validRange.second + equalityTolerance),
           };
 
-        if (validRange.first > validRange.second) continue;
-
         if (validRange.first > currentValue) {
           currentSlope -= constraint.weight;
           if (validRange.first < upper)
@@ -415,6 +420,9 @@ class JumpMove {
       if (eq(bestValue, varIncumbentValue, equalityTolerance) ||
           (!eq(currentValue, varIncumbentValue, equalityTolerance) &&
            currentScore < bestScore)) {
+        // Different to everywhere else (!)
+        // Here the score represents infeasibility (not improvement)
+        // Hence the minimization
         bestScore = currentScore;
         bestValue = currentValue;
       }
@@ -549,7 +557,7 @@ class FeasibilityJumpSolver {
       if (problem.vars.size() == 0) break;
 
       uint32_t var = selectVariable();
-      if (var == UINT_MAX){
+      if (var == UINT_MAX) {
         break;
       }
       doVariableMove(var);
@@ -653,7 +661,8 @@ class FeasibilityJumpSolver {
       dt += problem.vars.size();
       for (size_t varIdx = 0; varIdx < problem.vars.size(); varIdx += 1)
         forEachMove(varIdx, [&](Move& move) {
-          move.score += weightUpdateIncrement *
+          // -= to align objective minimization with score maximization
+          move.score -= weightUpdateIncrement *
                         problem.vars[varIdx].objectiveCoeff *
                         (move.value - problem.incumbentAssignment[varIdx]);
         });
@@ -754,7 +763,8 @@ class FeasibilityJumpSolver {
 
     forEachMove(varIdx, [&](Move& move) {
       move.score = 0.0;
-      move.score += objectiveWeight * problem.vars[varIdx].objectiveCoeff *
+      // -= to align objective minimization with score maximization
+      move.score -= objectiveWeight * problem.vars[varIdx].objectiveCoeff *
                     (move.value - problem.incumbentAssignment[varIdx]);
 
       for (auto& cell : problem.vars[varIdx].coeffs) {