Remove continuous singletons from double-sided rows

check/TestBasis.cpp

@@ -25,9 +25,7 @@ TEST_CASE("Basis-file", "[highs_basis_file]") {
       std::string(HIGHS_DIR) + "/check/instances/avgas.mps";
 
   Highs highs;
-  if (!dev_run) {
-    highs.setOptionValue("output_flag", false);
-  }
+  highs.setOptionValue("output_flag", dev_run);
   assert(model0_file != model1_file);
 
   return_status = highs.readModel(model0_file);

check/TestPresolve.cpp

@@ -192,11 +192,13 @@ void presolveSolvePostsolve(const std::string& model_file,
     if (dev_run)
       printf("Presolve timeout: return status = %d\n", (int)return_status);
   }
+  REQUIRE(model_presolve_status != HighsPresolveStatus::kUnboundedOrInfeasible);
   HighsLp lp = highs0.getPresolvedLp();
   highs1.passModel(lp);
   highs1.setOptionValue("solve_relaxation", solve_relaxation);
   highs1.setOptionValue("presolve", kHighsOffString);
   highs1.run();
+  REQUIRE(highs1.getModelStatus() == HighsModelStatus::kOptimal);
   HighsSolution solution = highs1.getSolution();
   const double objective_value = highs1.getInfo().objective_function_value;
   if (lp.isMip() && !solve_relaxation) {
@@ -910,3 +912,20 @@ TEST_CASE("presolve-issue-2874", "[highs_test_presolve]") {
   REQUIRE(highs.presolve() == HighsStatus::kOk);
   REQUIRE(highs.getModelPresolveStatus() == HighsPresolveStatus::kInfeasible);
 }
+
+TEST_CASE("issue-2962", "[highs_test_presolve]") {
+  // Model that exposes the case where zeroCostSingleton in postsolve
+  // requires the column to inherit the basic status of the row
+  std::string model_file =
+      std::string(HIGHS_DIR) + "/check/instances/p0548.mps";
+  Highs h;
+  h.setOptionValue("output_flag", dev_run);
+  h.setOptionValue("solve_relaxation", true);
+  h.setOptionValue("log_dev_level", 1);
+  h.setOptionValue("presolve_rule_logging", true);
+
+  h.readModel(model_file);
+  REQUIRE(h.run() == HighsStatus::kOk);
+
+  h.resetGlobalScheduler(true);
+}

highs/lp_data/HConst.h

@@ -275,7 +275,9 @@ enum PresolveRuleType : int {
   kPresolveRuleSparsify,
   kPresolveRuleProbing,
   kPresolveRuleEnumeration,
-  kPresolveRuleMax = kPresolveRuleEnumeration,
+  kPresolveRuleDualFixing,
+  kPresolveRuleZeroCostSingleton,
+  kPresolveRuleMax = kPresolveRuleZeroCostSingleton,
   kPresolveRuleLastAllowOff = kPresolveRuleMax,
   kPresolveRuleCount
 };

highs/lp_data/HighsModelUtils.cpp

@@ -1379,21 +1379,27 @@ std::string utilSolutionStatusToString(const HighsInt solution_status) {
 }
 
 // Return a string representation of HighsBasisStatus
-std::string utilBasisStatusToString(const HighsBasisStatus basis_status) {
+std::string utilBasisStatusToString(const HighsBasisStatus basis_status,
+                                    const bool s2) {
   switch (basis_status) {
     case HighsBasisStatus::kLower:
+      if (s2) return "LO";
       return "At lower/fixed bound";
       break;
     case HighsBasisStatus::kBasic:
+      if (s2) return "BS";
       return "Basic";
       break;
     case HighsBasisStatus::kUpper:
+      if (s2) return "UP";
       return "At upper bound";
       break;
     case HighsBasisStatus::kZero:
+      if (s2) return "ZE";
       return "Free at zero";
       break;
     case HighsBasisStatus::kNonbasic:
+      if (s2) return "NB";
       return "Nonbasic";
       break;
     default:
@@ -1515,6 +1521,10 @@ std::string utilPresolveRuleTypeToString(const HighsInt rule_type) {
     return "Probing";
   } else if (rule_type == kPresolveRuleEnumeration) {
     return "Enumeration";
+  } else if (rule_type == kPresolveRuleDualFixing) {
+    return "Dual fixing";
+  } else if (rule_type == kPresolveRuleZeroCostSingleton) {
+    return "Zero cost singleton";
   }
   assert(1 == 0);
   return "????";

highs/lp_data/HighsModelUtils.h

@@ -109,7 +109,8 @@ std::string utilModelStatusToString(const HighsModelStatus model_status);
 
 std::string utilSolutionStatusToString(const HighsInt solution_status);
 
-std::string utilBasisStatusToString(const HighsBasisStatus basis_status);
+std::string utilBasisStatusToString(const HighsBasisStatus basis_status,
+                                    const bool s2 = false);
 
 std::string utilBasisValidityToString(const HighsInt basis_validity);
 

highs/lp_data/HighsSolution.cpp

@@ -2051,8 +2051,8 @@ bool reportKktFailures(const HighsLp& lp, const HighsOptions& options,
   double primal_residual_tolerance = options.primal_residual_tolerance;
   double dual_residual_tolerance = options.dual_residual_tolerance;
   double optimality_tolerance = options.optimality_tolerance;
-  const bool is_mip = lp.isMip();
-  if (is_mip) {
+  const bool solved_as_mip = lp.isMip() && !options.solve_relaxation;
+  if (solved_as_mip) {
     primal_feasibility_tolerance = mip_feasibility_tolerance;
   } else if (options.kkt_tolerance != kDefaultKktTolerance) {
     mip_feasibility_tolerance = options.kkt_tolerance;
@@ -2065,9 +2065,10 @@ bool reportKktFailures(const HighsLp& lp, const HighsOptions& options,
 
   const bool force_report = options.log_dev_level >= kHighsLogDevLevelInfo;
   const bool complementarity_error =
-      !is_mip && info.primal_dual_objective_error > optimality_tolerance;
+      !solved_as_mip && info.primal_dual_objective_error > optimality_tolerance;
   const bool integrality_error =
-      is_mip && info.max_integrality_violation >= mip_feasibility_tolerance;
+      solved_as_mip &&
+      info.max_integrality_violation >= mip_feasibility_tolerance;
   const bool has_kkt_failures =
       integrality_error || info.num_primal_infeasibilities > 0 ||
       info.num_dual_infeasibilities > 0 ||
@@ -2080,7 +2081,7 @@ bool reportKktFailures(const HighsLp& lp, const HighsOptions& options,
 
   highsLogUser(log_options, log_type, "Solution optimality conditions%s%s\n",
                message == "" ? "" : ": ", message == "" ? "" : message.c_str());
-  if (is_mip && info.max_integrality_violation >= 0)
+  if (solved_as_mip && info.max_integrality_violation >= 0)
     highsLogUser(log_options, HighsLogType::kInfo,
                  "    max      %8.3g                                  "
                  "integrality violations"

highs/presolve/HPresolve.cpp

@@ -2144,7 +2144,6 @@ HighsTripletTreeSliceInOrder HPresolve::getSortedRowVector(HighsInt row) const {
 
 void HPresolve::markRowDeleted(HighsInt row) {
   assert(!rowDeleted[row]);
-
   // remove equations from set of equations
   if (isEquation(row) && eqiters[row] != equations.end()) {
     equations.erase(eqiters[row]);
@@ -2154,15 +2153,16 @@ void HPresolve::markRowDeleted(HighsInt row) {
   // prevents row from being added to change vector
   changedRowFlag[row] = true;
   rowDeleted[row] = true;
+  assert(!analysis_.logging_on_);
   ++numDeletedRows;
 }
 
 void HPresolve::markColDeleted(HighsInt col) {
   assert(!colDeleted[col]);
-
   // prevents col from being added to change vector
   changedColFlag[col] = true;
   colDeleted[col] = true;
+  assert(!analysis_.logging_on_);
   ++numDeletedCols;
 }
 
@@ -3369,8 +3369,10 @@ HPresolve::Result HPresolve::singletonCol(HighsPostsolveStack& postsolve_stack,
         static_cast<Result>(convertImpliedInteger(col, row)));
 
   // dual fixing
-  HPRESOLVE_CHECKED_CALL(dualFixing(postsolve_stack, col));
-  if (colDeleted[col]) return Result::kOk;
+  if (analysis_.allow_rule_[kPresolveRuleDualFixing]) {
+    HPRESOLVE_CHECKED_CALL(dualFixing(postsolve_stack, col));
+    if (colDeleted[col]) return Result::kOk;
+  }
 
   // singleton column stuffing
   HPRESOLVE_CHECKED_CALL(singletonColStuffing(postsolve_stack, col));
@@ -3409,7 +3411,11 @@ HPresolve::Result HPresolve::singletonCol(HighsPostsolveStack& postsolve_stack,
     return checkLimits(postsolve_stack);
   }
 
-  // todo: check for zero cost singleton and remove
+  if (analysis_.allow_rule_[kPresolveRuleZeroCostSingleton]) {
+    // Remove if col is double-sided finite slack
+    HPRESOLVE_CHECKED_CALL(zeroCostSingleton(postsolve_stack, col));
+  }
+
   return Result::kOk;
 }
 
@@ -4493,8 +4499,10 @@ HPresolve::Result HPresolve::colPresolve(HighsPostsolveStack& postsolve_stack,
   }
 
   // dual fixing
-  HPRESOLVE_CHECKED_CALL(dualFixing(postsolve_stack, col));
-  if (colDeleted[col]) return Result::kOk;
+  if (analysis_.allow_rule_[kPresolveRuleDualFixing]) {
+    HPRESOLVE_CHECKED_CALL(dualFixing(postsolve_stack, col));
+    if (colDeleted[col]) return Result::kOk;
+  }
 
   // singleton column stuffing
   HPRESOLVE_CHECKED_CALL(singletonColStuffing(postsolve_stack, col));
@@ -4913,13 +4921,19 @@ HPresolve::Result HPresolve::dualFixing(HighsPostsolveStack& postsolve_stack,
   // compute locks
   computeLocks(col, true, lockCallback);
 
-  // check if variable can be fixed
+  // Check if variable is fixed: only possible if there are no up
+  // (down) locks, and the cost forces it up (down) to its bound
   if (numDownLocks == 0 || numUpLocks == 0) {
-    // fix variable
-    if (numDownLocks == 0)
+    // fix variable if cost is driving it to its bound
+    if (numDownLocks == 0 &&
+        model->col_cost_[col] <= options->dual_feasibility_tolerance &&
+        model->col_lower_[col] > -kHighsInf) {
       HPRESOLVE_CHECKED_CALL(fixColToLower(postsolve_stack, col));
-    else
+    } else if (numUpLocks == 0 &&
+               model->col_cost_[col] >= -options->dual_feasibility_tolerance &&
+               model->col_upper_[col] < kHighsInf) {
       HPRESOLVE_CHECKED_CALL(fixColToUpper(postsolve_stack, col));
+    }
   } else {
     bool hasSingleDownLock = numDownLocks == 1 && downLockRow != -1;
     bool hasSingleUpLock = numUpLocks == 1 && upLockRow != -1;
@@ -5206,6 +5220,68 @@ HPresolve::Result HPresolve::singletonColStuffing(
   return Result::kOk;
 }
 
+HPresolve::Result HPresolve::zeroCostSingleton(
+    HighsPostsolveStack& postsolve_stack, HighsInt col) {
+  assert(analysis_.allow_rule_[kPresolveRuleZeroCostSingleton]);
+  // For a double-sided row b_0 <= a^Tx + cs <= b_1,
+  // where s is a singleton continuous column with 0 cost,
+  // relax s out as its value can be determined in postsolve.
+  // The row may now admit additional reductions afterwards.
+  // Dual fixing already handles single-sided row case with fixings.
+  if (model->integrality_[col] != HighsVarType::kContinuous ||
+      model->col_cost_[col] != 0.0 || colsize[col] != 1) {
+    return Result::kOk;
+  }
+  assert(!colDeleted[col]);
+
+  HighsInt nzPos = colhead[col];
+  HighsInt row = Arow[nzPos];
+  assert(!rowDeleted[row]);
+  // Row must be ranged, but why can't it be an equation?
+  const bool was_equation = isEquation(row);
+  if (!isRanged(row)) return Result::kOk;
+
+  double coef = Avalue[nzPos];
+
+  if (std::abs(coef) == kHighsInf) return Result::kOk;
+
+  const bool logging_on = analysis_.logging_on_;
+  if (logging_on)
+    analysis_.startPresolveRuleLog(kPresolveRuleZeroCostSingleton);
+  storeRow(row);
+
+  double lb = model->col_lower_[col];
+  double ub = model->col_upper_[col];
+  double change_from_col_lb = coef * lb;
+  double change_from_col_ub = coef * ub;
+
+  double newRowLower =
+      model->row_lower_[row] - std::max(change_from_col_lb, change_from_col_ub);
+  double newRowUpper =
+      model->row_upper_[row] - std::min(change_from_col_lb, change_from_col_ub);
+
+  postsolve_stack.zeroCostSingleton(row, col, model->row_lower_[row],
+                                    model->row_upper_[row], newRowLower,
+                                    newRowUpper, lb, ub, coef, getStoredRow());
+
+  model->row_lower_[row] = newRowLower;
+  model->row_upper_[row] = newRowUpper;
+  if (was_equation && newRowLower != newRowUpper &&
+      eqiters[row] != equations.end()) {
+    equations.erase(eqiters[row]);
+    eqiters[row] = equations.end();
+  }
+
+  // Delete the singleton column
+  markColDeleted(col);
+  unlink(nzPos);
+
+  analysis_.logging_on_ = logging_on;
+  if (logging_on) analysis_.stopPresolveRuleLog(kPresolveRuleZeroCostSingleton);
+
+  return Result::kOk;
+}
+
 HPresolve::Result HPresolve::enumerateSolutions(
     HighsPostsolveStack& postsolve_stack) {
   // enumerate all solutions for pure binary constraints with a small number of

highs/presolve/HPresolve.h

@@ -421,6 +421,8 @@ class HPresolve {
   Result singletonColStuffing(HighsPostsolveStack& postsolve_stack,
                               HighsInt col);
 
+  Result zeroCostSingleton(HighsPostsolveStack& postsolve_stack, HighsInt col);
+
   Result enumerateSolutions(HighsPostsolveStack& postsolve_stack);
 
   double computeImpliedLowerBound(HighsInt col, HighsInt boundCol = -1,

highs/presolve/HPresolveAnalysis.cpp

@@ -48,7 +48,7 @@ void HPresolveAnalysis::setup(const HighsLp* model_,
           if (!allow ||
               (!options->presolve_rule_off && options_->log_dev_level))
             highsLogUser(options->log_options, HighsLogType::kInfo,
-                         "   Rule %2d (set bit %2d = %5d): %s\n",
+                         "   Rule %2d (set bit %2d = %6d): %s\n",
                          int(rule_type), int(rule_type), int(bit),
                          utilPresolveRuleTypeToString(rule_type).c_str());
       } else if (!allow && !silent) {
@@ -65,6 +65,21 @@ void HPresolveAnalysis::setup(const HighsLp* model_,
   }
   // Allow logging if option is set and model is not a MIP
   allow_logging_ = options_->presolve_rule_logging && !model_->isMip();
+  // NB logging_on_ is also used to determine whether logging has
+  // started to prevent double-accounting. Specifically,
+  //
+  // * If allow_logging_ is false, then startPresolveRuleLog is never called
+  //
+  // * If allow_logging_ is true, then startPresolveRuleLog will be
+  //   called at some point, and within startPresolveRuleLog,
+  //   allow_logging_ is set to false. Hence, it's never called again -
+  //   which would lead to double-accounting. When stopPresolveRuleLog
+  //   is called - and this is done by saving the value of allow_logging_
+  //   before a posible call to startPresolveRuleLog - allow_logging_
+  //   is set to true
+  //
+  // Since all reductions are logged, logging_on_ must be false in all
+  // calls to markRowDeleted and markColDeleted
   logging_on_ = allow_logging_;
   log_rule_type_ = kPresolveRuleIllegal;
   resetNumDeleted();