Keep small terms big domains lp aggregator

check/TestMipSolver.cpp

@@ -2678,6 +2678,23 @@ TEST_CASE("issue-3461", "[highs_test_mip_solver]") {
   highs.resetThreadScheduler(true);
 }
 
+TEST_CASE("issue-3478", "[highs_test_mip_solver]") {
+  // model tolls.lp in issue 3478. HighsLpAggregator dropped coefficients with
+  // |val| <= small_matrix_value without accounting for their contribution
+  // |val| * bound to the aggregated row, which led to invalid cuts and the
+  // wrong optimal objective -6016.354106
+  std::string filename =
+      std::string(HIGHS_DIR) + "/check/instances/issue-3478.lp";
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+  highs.setOptionValue("mip_rel_gap", 0);
+  REQUIRE(highs.readModel(filename) == HighsStatus::kOk);
+  const HighsModelStatus require_model_status = HighsModelStatus::kOptimal;
+  const double optimal_objective = -6254.734457;
+  solve(highs, kHighsOnString, require_model_status, optimal_objective);
+  highs.resetThreadScheduler(true);
+}
+
 TEST_CASE("conflict-pool-sync-empty-conflict", "[highs_test_mip_solver]") {
   // an empty conflict, which can be derived when a proof is violated by the
   // global bounds alone, is stored at the end of the conflict entries and

check/instances/issue-3478.lp

@@ -0,0 +1,121 @@
+Minimize
+ obj: - 2421 x0 - 1.2105 v0 + 0.0001 z0_0 + 0.0002 z0_8 + 0.0002 z0_5 + 0.0003 z0_2 + 0.0003 z0_3 + 0.0001 z0_1 + 0.0001 z0_9 + 0.0003 z0_7 + 0.0002 z0_10 + 0.0003 z0_6 + 0.0002 z0_11 + 0.0001 z0_4
+ - 2662 x1 - 1.331 v1 + 0.0003 z1_6 + 0.0003 z1_3 + 0.0002 z1_11 + 0.0001 z1_0 + 0.0001 z1_9 + 0.0003 z1_2 + 0.0003 z1_7 + 0.0001 z1_1 + 0.0002 z1_5 + 0.0001 z1_4 + 0.0002 z1_8 + 0.0002 z1_10
+ - 3017 x2 - 1.5085 v2 + 0.0002 z2_10 + 0.0002 z2_11 + 0.0002 z2_8 + 0.0001 z2_9 + 0.0001 z2_1 + 0.0003 z2_6 + 0.0002 z2_5 + 0.0001 z2_0 + 0.0001 z2_4 + 0.0003 z2_2 + 0.0003 z2_7 + 0.0003 z2_3
+ - 3569 x3 - 1.7845 v3 + 0.0001 z3_4 + 0.0001 z3_0 + 0.0002 z3_11 + 0.0003 z3_2 + 0.0002 z3_5 + 0.0003 z3_7 + 0.0002 z3_10 + 0.0001 z3_1 + 0.0002 z3_8 + 0.0003 z3_6 + 0.0001 z3_9 + 0.0003 z3_3
+ - 221 x4 - 0.1105 v4 + 0.0001 z4_9 + 0.0001 z4_0 + 0.0002 z4_5 + 0.0001 z4_4 + 0.0001 z4_1 + 0.0003 z4_6
+ - 272 x5 - 0.136 v5 + 0.0002 z5_11 + 0.0002 z5_5 + 0.0003 z5_2 + 0.0003 z5_7 + 0.0001 z5_1
+Subject To
+ l0: w0 - 1e+06 x0 = 0
+ c0: w0 - 0.0008 z0_0 - 0.0008 z0_8 - 0.0008 z0_5 - 0.0008 z0_2 - 0.0008 z0_3 - 0.0008 z0_1 - 0.0008 z0_9 - 0.0008 z0_7 - 0.0008 z0_10 - 0.0008 z0_6 - 0.0008 z0_11 - 0.0008 z0_4 + 1000 v0 <= 984401
+ l1: w1 - 1e+06 x1 = 0
+ c1: w1 - 0.0008 z1_6 - 0.0008 z1_3 - 0.0008 z1_11 - 0.0008 z1_0 - 0.0008 z1_9 - 0.0008 z1_2 - 0.0008 z1_7 - 0.0008 z1_1 - 0.0008 z1_5 - 0.0008 z1_4 - 0.0008 z1_8 - 0.0008 z1_10 + 1000 v1 <= 990503
+ l2: w2 - 1e+06 x2 = 0
+ c2: w2 - 0.0008 z2_10 - 0.0008 z2_11 - 0.0008 z2_8 - 0.0008 z2_9 - 0.0008 z2_1 - 0.0008 z2_6 - 0.0008 z2_5 - 0.0008 z2_0 - 0.0008 z2_4 - 0.0008 z2_2 - 0.0008 z2_7 - 0.0008 z2_3 + 1000 v2 <= 986979
+ l3: w3 - 1e+06 x3 = 0
+ c3: w3 - 0.0008 z3_4 - 0.0008 z3_0 - 0.0008 z3_11 - 0.0008 z3_2 - 0.0008 z3_5 - 0.0008 z3_7 - 0.0008 z3_10 - 0.0008 z3_1 - 0.0008 z3_8 - 0.0008 z3_6 - 0.0008 z3_9 - 0.0008 z3_3 + 1000 v3 <= 992642
+ l4: w4 - 1e+06 x4 = 0
+ c4: w4 - 0.0008 z4_9 - 0.0008 z4_0 - 0.0008 z4_5 - 0.0008 z4_4 - 0.0008 z4_1 - 0.0008 z4_6 + 1000 v4 <= 998344
+ l5: w5 - 1e+06 x5 = 0
+ c5: w5 - 0.0008 z5_11 - 0.0008 z5_5 - 0.0008 z5_2 - 0.0008 z5_7 - 0.0008 z5_1 + 1000 v5 <= 998770
+ road0: + z0_0 + z1_0 + z2_0 + z3_0 + z4_0 <= 1000000
+ road1: + z0_1 + z1_1 + z2_1 + z3_1 + z4_1 + z5_1 <= 1000000
+ road2: + z0_2 + z1_2 + z2_2 + z3_2 + z5_2 <= 1000000
+ road3: + z0_3 + z1_3 + z2_3 + z3_3 <= 1000000
+ road4: + z0_4 + z1_4 + z2_4 + z3_4 + z4_4 <= 1000000
+ road5: + z0_5 + z1_5 + z2_5 + z3_5 + z4_5 + z5_5 <= 1000000
+ road6: + z0_6 + z1_6 + z2_6 + z3_6 + z4_6 <= 1000000
+ road7: + z0_7 + z1_7 + z2_7 + z3_7 + z5_7 <= 1000000
+ road8: + z0_8 + z1_8 + z2_8 + z3_8 <= 1000000
+ road9: + z0_9 + z1_9 + z2_9 + z3_9 + z4_9 <= 1000000
+ road10: + z0_10 + z1_10 + z2_10 + z3_10 <= 1000000
+ road11: + z0_11 + z1_11 + z2_11 + z3_11 + z5_11 <= 1000000
+ budget0: + 43 x0 + 16 x1 + 14 x2 + 14 x3 + 38 x4 + 46 x5 <= 85
+ budget1: + 47 x0 + 25 x1 + 22 x2 + 19 x3 + 39 x4 + 50 x5 <= 101
+ budget2: + 51 x0 + 54 x1 + 52 x2 + 60 x3 + 28 x4 + 44 x5 <= 144
+Bounds
+ 0 <= x0 <= 1
+ 0 <= w0 <= 1000000
+ 0 <= v0 <= 1500
+ 0 <= z0_0 <= 1000000
+ 0 <= z0_8 <= 1000000
+ 0 <= z0_5 <= 1000000
+ 0 <= z0_2 <= 1000000
+ 0 <= z0_3 <= 1000000
+ 0 <= z0_1 <= 1000000
+ 0 <= z0_9 <= 1000000
+ 0 <= z0_7 <= 1000000
+ 0 <= z0_10 <= 1000000
+ 0 <= z0_6 <= 1000000
+ 0 <= z0_11 <= 1000000
+ 0 <= z0_4 <= 1000000
+ 0 <= x1 <= 1
+ 0 <= w1 <= 1000000
+ 0 <= v1 <= 1500
+ 0 <= z1_6 <= 1000000
+ 0 <= z1_3 <= 1000000
+ 0 <= z1_11 <= 1000000
+ 0 <= z1_0 <= 1000000
+ 0 <= z1_9 <= 1000000
+ 0 <= z1_2 <= 1000000
+ 0 <= z1_7 <= 1000000
+ 0 <= z1_1 <= 1000000
+ 0 <= z1_5 <= 1000000
+ 0 <= z1_4 <= 1000000
+ 0 <= z1_8 <= 1000000
+ 0 <= z1_10 <= 1000000
+ 0 <= x2 <= 1
+ 0 <= w2 <= 1000000
+ 0 <= v2 <= 1500
+ 0 <= z2_10 <= 1000000
+ 0 <= z2_11 <= 1000000
+ 0 <= z2_8 <= 1000000
+ 0 <= z2_9 <= 1000000
+ 0 <= z2_1 <= 1000000
+ 0 <= z2_6 <= 1000000
+ 0 <= z2_5 <= 1000000
+ 0 <= z2_0 <= 1000000
+ 0 <= z2_4 <= 1000000
+ 0 <= z2_2 <= 1000000
+ 0 <= z2_7 <= 1000000
+ 0 <= z2_3 <= 1000000
+ 0 <= x3 <= 1
+ 0 <= w3 <= 1000000
+ 0 <= v3 <= 1500
+ 0 <= z3_4 <= 1000000
+ 0 <= z3_0 <= 1000000
+ 0 <= z3_11 <= 1000000
+ 0 <= z3_2 <= 1000000
+ 0 <= z3_5 <= 1000000
+ 0 <= z3_7 <= 1000000
+ 0 <= z3_10 <= 1000000
+ 0 <= z3_1 <= 1000000
+ 0 <= z3_8 <= 1000000
+ 0 <= z3_6 <= 1000000
+ 0 <= z3_9 <= 1000000
+ 0 <= z3_3 <= 1000000
+ 0 <= x4 <= 1
+ 0 <= w4 <= 1000000
+ 0 <= v4 <= 1500
+ 0 <= z4_9 <= 1000000
+ 0 <= z4_0 <= 1000000
+ 0 <= z4_5 <= 1000000
+ 0 <= z4_4 <= 1000000
+ 0 <= z4_1 <= 1000000
+ 0 <= z4_6 <= 1000000
+ 0 <= x5 <= 1
+ 0 <= w5 <= 1000000
+ 0 <= v5 <= 1500
+ 0 <= z5_11 <= 1000000
+ 0 <= z5_5 <= 1000000
+ 0 <= z5_2 <= 1000000
+ 0 <= z5_7 <= 1000000
+ 0 <= z5_1 <= 1000000
+Binaries
+ x0
+ x1
+ x2
+ x3
+ x4
+ x5
+End

highs/mip/HighsDomain.cpp

@@ -68,13 +68,6 @@ static inline double boundRange(double upper_bound, double lower_bound,
                       : tolerance);
 }
 
-bool HighsDomain::termIsZero(double val, double lb, double ub, double epsilon) {
-  const double absVal = std::abs(val);
-  if (absVal <= std::numeric_limits<double>::min()) return true;
-  return absVal <= epsilon &&
-         absVal * std::max(std::abs(lb), std::abs(ub)) <= epsilon;
-}
-
 HighsDomain::HighsDomain(HighsMipSolver& mipsolver) : mipsolver(&mipsolver) {
   col_lower_ = mipsolver.model_->col_lower_;
   col_upper_ = mipsolver.model_->col_upper_;

highs/mip/HighsDomain.h

@@ -557,8 +557,6 @@ class HighsDomain {
     return *this;
   }
 
-  static bool termIsZero(double val, double lb, double ub, double epsilon);
-
   void computeMinActivity(HighsInt start, HighsInt end, const HighsInt* ARindex,
                           const double* ARvalue, HighsInt& ninfmin,
                           HighsCDouble& activitymin) const;

highs/mip/HighsLpAggregator.cpp

@@ -9,6 +9,7 @@
 #include "mip/HighsLpAggregator.h"
 
 #include "mip/HighsLpRelaxation.h"
+#include "mip/HighsMipSolverData.h"
 
 HighsLpAggregator::HighsLpAggregator(const HighsLpRelaxation& lprelaxation)
     : lprelaxation(lprelaxation) {
@@ -32,15 +33,17 @@ void HighsLpAggregator::getCurrentAggregation(std::vector<HighsInt>& inds,
                                               bool negate) {
   const double droptol =
       lprelaxation.getMipSolver().options_mip_->small_matrix_value;
+  const HighsDomain& domain = lprelaxation.getMipSolver().mipdata_->getDomain();
   const HighsInt numCol = lprelaxation.numCols();
-  vectorsum.cleanup([droptol, numCol](HighsInt col, double val) {
+  vectorsum.cleanup([&](const HighsInt col, const double val) {
     // only drop values for columns and not for slack variables of rows as the
     // former are the only ones which might be subject to numerical error. The
     // values for row slack variables are exact copies of the weights used
     // (assuming that separators only add each row once). Also all the
     // separators will only allow rows that do have meaningful coefficient
     // contributions with the used weight
-    return col < numCol && std::fabs(val) <= droptol;
+    return col < numCol && termIsZero(val, domain.col_lower_[col],
+                                      domain.col_upper_[col], droptol);
   });
 
   inds = vectorsum.getNonzeros();

highs/mip/HighsLpRelaxation.cpp

@@ -876,9 +876,8 @@ bool HighsLpRelaxation::computeDualProof(const HighsDomain& globaldomain,
 
     double val = static_cast<double>(exactVal);
 
-    if (HighsDomain::termIsZero(val, globaldomain.col_lower_[i],
-                                globaldomain.col_upper_[i],
-                                mipsolver.mipdata_->epsilon))
+    if (termIsZero(val, globaldomain.col_lower_[i], globaldomain.col_upper_[i],
+                   mipsolver.mipdata_->epsilon))
       continue;
 
     bool removeValue = std::fabs(val) <= mipsolver.mipdata_->feastol;
@@ -994,9 +993,8 @@ void HighsLpRelaxation::storeDualInfProof() {
     const HighsCDouble exactVal = row_ap.getValue<HighsCDouble>(i);
     double val = static_cast<double>(exactVal);
 
-    if (HighsDomain::termIsZero(val, globaldomain.col_lower_[i],
-                                globaldomain.col_upper_[i],
-                                mipsolver.mipdata_->epsilon))
+    if (termIsZero(val, globaldomain.col_lower_[i], globaldomain.col_upper_[i],
+                   mipsolver.mipdata_->epsilon))
       continue;
 
     bool removeValue = std::abs(val) <= mipsolver.mipdata_->feastol;

highs/mip/HighsTransformedLp.cpp

@@ -374,8 +374,7 @@ bool HighsTransformedLp::transform(std::vector<double>& vals,
 
     double maxError = 0.0;
     auto IsZero = [&](HighsInt col, double val) {
-      if (HighsDomain::termIsZero(val, getLb(col), getUb(col),
-                                  mip.mipdata_->epsilon))
+      if (termIsZero(val, getLb(col), getUb(col), mip.mipdata_->epsilon))
         return true;
       return false;
     };
@@ -577,9 +576,8 @@ bool HighsTransformedLp::untransform(std::vector<double>& vals,
     bool abort = false;
     auto IsZero = [&](HighsInt col, double val) {
       assert(col < mip.numCol());
-      if (HighsDomain::termIsZero(val, globaldom_.col_lower_[col],
-                                  globaldom_.col_upper_[col],
-                                  mip.mipdata_->epsilon))
+      if (termIsZero(val, globaldom_.col_lower_[col],
+                     globaldom_.col_upper_[col], mip.mipdata_->epsilon))
         return true;
       double absval = std::abs(val);
 

highs/util/HighsUtils.cpp

@@ -1231,6 +1231,14 @@ double nearestPowerOfTwoScale(const double value) {
   return scale;
 }
 
+bool termIsZero(const double val, const double lb, const double ub,
+                const double epsilon) {
+  const double absVal = std::abs(val);
+  if (absVal <= std::numeric_limits<double>::min()) return true;
+  return absVal <= epsilon &&
+         absVal * std::max(std::abs(lb), std::abs(ub)) <= epsilon;
+}
+
 void highsAssert(const bool assert_condition, const std::string& message) {
   if (assert_condition) return;
   printf("Failing highsAssert(\"%s\")\n", message.c_str());

highs/util/HighsUtils.h

@@ -198,6 +198,11 @@ bool computeScatterDataRegressionError(HighsScatterData& scatter_data,
 
 double nearestPowerOfTwoScale(const double value);
 
+// returns true if the term val * x with x in [lb, ub] is negligible, i.e. if
+// |val| and |val| * max(|lb|, |ub|) do not exceed epsilon
+bool termIsZero(const double val, const double lb, const double ub,
+                const double epsilon);
+
 // If assert_condition is false then, if NDEBUG is defined message is
 // printed and abort() is called, otherwise assert is called
 void highsAssert(const bool assert_condition, const std::string& message = "");