Remove continuous singletons from double-sided rows

highs/presolve/HPresolve.cpp

@@ -3403,7 +3403,10 @@ HPresolve::Result HPresolve::singletonCol(HighsPostsolveStack& postsolve_stack,
     return checkLimits(postsolve_stack);
   }
 
-  // todo: check for zero cost singleton and remove
+  // Remove if col is double-sided finite slack
+  HPRESOLVE_CHECKED_CALL(
+      redundantSingletonColDoubleSidedSlack(postsolve_stack, col));
+
   return Result::kOk;
 }
 
@@ -5147,6 +5150,55 @@ HPresolve::Result HPresolve::singletonColStuffing(
   return Result::kOk;
 }
 
+HPresolve::Result HPresolve::redundantSingletonColDoubleSidedSlack(
+    HighsPostsolveStack& postsolve_stack, HighsInt col) {
+  // 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 ||
+      model->col_lower_[col] == -kHighsInf ||
+      model->col_upper_[col] == kHighsInf) {
+    return Result::kOk;
+  }
+  assert(!colDeleted[col]);
+
+  HighsInt nzPos = colhead[col];
+  HighsInt row = Arow[nzPos];
+  assert(!rowDeleted[row]);
+  if (!isRanged(row) || isEquation(row)) return Result::kOk;
+  double coef = Avalue[nzPos];
+
+  if (std::abs(coef) == kHighsInf) return Result::kOk;
+
+  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.zeroObjSingletonContinuousCol(
+      row, col, model->row_lower_[row], model->row_upper_[row], lb, ub, coef,
+      getStoredRow());
+
+  model->row_lower_[row] = newRowLower;
+  model->row_upper_[row] = newRowUpper;
+
+  // Delete the singleton column
+  markColDeleted(col);
+  unlink(nzPos);
+
+  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

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

highs/presolve/HighsPostsolveStack.cpp

@@ -1365,4 +1365,71 @@ void HighsPostsolveStack::SlackColSubstitution::undo(
   }
 }
 
+void HighsPostsolveStack::ZeroObjSingletonContinuousCol::undo(
+    const HighsOptions& options, const std::vector<Nonzero>& rowValues,
+    HighsSolution& solution, HighsBasis& basis) {
+  // a (removed) cut may have been used in this reduction.
+  bool isModelRow = static_cast<size_t>(row) < solution.row_value.size();
+
+  assert(origRowLower != -kHighsInf && origRowUpper != kHighsInf);
+
+  // Get activity of row without removed singleton
+  HighsCDouble act = 0;
+  for (const auto& rowVal : rowValues) {
+    if (rowVal.index != col) {
+      act += rowVal.value * solution.col_value[rowVal.index];
+    }
+  }
+
+  // Determine domain of potential values
+  double col_lb = lb;
+  double col_ub = ub;
+  if (coef > 0) {
+    col_lb = std::max(col_lb, static_cast<double>((origRowLower - act) / coef));
+    col_ub = std::min(col_ub, static_cast<double>((origRowUpper - act) / coef));
+  } else {
+    col_lb = std::max(col_lb, static_cast<double>((origRowUpper - act) / coef));
+    col_ub = std::min(col_ub, static_cast<double>((origRowLower - act) / coef));
+  }
+
+  // Find a suitable value within the allowed interval
+  if (basis.valid && isModelRow) {
+    if (basis.row_status[row] == HighsBasisStatus::kLower) {
+      solution.col_value[col] = coef > 0 ? col_lb : col_ub;
+    } else if (basis.row_status[row] == HighsBasisStatus::kUpper) {
+      solution.col_value[col] = coef > 0 ? col_ub : col_lb;
+    } else if (0 >= col_lb && 0 <= col_ub) {
+      solution.col_value[col] = 0;
+    } else {
+      solution.col_value[col] = 0.5 * (col_lb + col_ub);
+    }
+  } else {
+    if (0 >= col_lb && 0 <= col_ub) {
+      solution.col_value[col] = 0;
+    } else {
+      solution.col_value[col] = 0.5 * (col_lb + col_ub);
+    }
+  }
+
+  if (isModelRow) {
+    solution.row_value[row] =
+        static_cast<double>(act + (solution.col_value[col] * coef));
+  }
+
+  if (!solution.dual_valid) return;
+
+  solution.col_dual[col] = (isModelRow) ? -coef * solution.row_dual[row] : 0;
+
+  if (!basis.valid) return;
+
+  if (solution.col_value[col] - options.dual_feasibility_tolerance <= lb) {
+    basis.col_status[col] = HighsBasisStatus::kLower;
+  } else if (solution.col_value[col] + options.dual_feasibility_tolerance >=
+             ub) {
+    basis.col_status[col] = HighsBasisStatus::kUpper;
+  } else {
+    basis.col_status[col] = HighsBasisStatus::kBasic;
+  }
+}
+
 }  // namespace presolve

highs/presolve/HighsPostsolveStack.h

@@ -227,6 +227,20 @@ class HighsPostsolveStack {
               HighsBasis& basis);
   };
 
+  struct ZeroObjSingletonContinuousCol {
+    double origRowLower;
+    double origRowUpper;
+    double lb;
+    double ub;
+    double coef;
+    HighsInt col;
+    HighsInt row;
+
+    void undo(const HighsOptions& options,
+              const std::vector<Nonzero>& rowValues, HighsSolution& solution,
+              HighsBasis& basis);
+  };
+
   /// tags for reduction
   enum class ReductionType : uint8_t {
     kLinearTransform,
@@ -243,6 +257,7 @@ class HighsPostsolveStack {
     kDuplicateRow,
     kDuplicateColumn,
     kSlackColSubstitution,
+    kZeroObjSingletonContinuousCol,
   };
 
   HighsDataStack reductionValues;
@@ -345,6 +360,22 @@ class HighsPostsolveStack {
     reductionAdded(ReductionType::kSlackColSubstitution);
   }
 
+  template <typename RowStorageFormat>
+  void zeroObjSingletonContinuousCol(
+      HighsInt row, HighsInt col, double origRowLower, double origRowUpper,
+      double lb, double ub, double coef,
+      const HighsMatrixSlice<RowStorageFormat>& rowVec) {
+    rowValues.clear();
+    for (const HighsSliceNonzero& rowVal : rowVec)
+      rowValues.emplace_back(origColIndex[rowVal.index()], rowVal.value());
+
+    reductionValues.push(
+        ZeroObjSingletonContinuousCol{origRowLower, origRowUpper, lb, ub, coef,
+                                      origColIndex[col], origRowIndex[row]});
+    reductionValues.push(rowValues);
+    reductionAdded(ReductionType::kZeroObjSingletonContinuousCol);
+  }
+
   template <typename ColStorageFormat>
   void doubletonEquation(HighsInt row, HighsInt colSubst, HighsInt col,
                          double coefSubst, double coef, double rhs,
@@ -741,6 +772,13 @@ class HighsPostsolveStack {
           reduction.undo(options, rowValues, solution, basis);
           break;
         }
+        case ReductionType::kZeroObjSingletonContinuousCol: {
+          ZeroObjSingletonContinuousCol reduction;
+          reductionValues.pop(rowValues);
+          reductionValues.pop(reduction);
+          reduction.undo(options, rowValues, solution, basis);
+          break;
+        }
         default:
           printf("Reduction case %d not handled\n",
                  int(reductions[i - 1].first));
@@ -923,6 +961,13 @@ class HighsPostsolveStack {
           reduction.undo(options, rowValues, solution, basis);
           break;
         }
+        case ReductionType::kZeroObjSingletonContinuousCol: {
+          ZeroObjSingletonContinuousCol reduction;
+          reductionValues.pop(rowValues);
+          reductionValues.pop(reduction);
+          reduction.undo(options, rowValues, solution, basis);
+          break;
+        }
       }
     }
 #ifdef DEBUG_EXTRA