Now getting through to postsolve after eliminating a slack in presolve for adlittle

Author: jajhall

src/lp_data/Highs.cpp

@@ -1252,7 +1252,7 @@ HighsStatus Highs::run() {
     bool have_optimal_solution = false;
     // ToDo Put solution of presolved problem in a separate method
 
-    if (!this->options_.presolve_remove_slacks) {
+    //    if (!this->options_.presolve_remove_slacks) {
       HighsLp& reduced_lp = presolve_.getReducedProblem();
       HighsInt num_double_slack = 0;
       HighsInt num_slack = 0;
@@ -1271,19 +1271,20 @@ HighsStatus Highs::run() {
 	  continue;
 	}
 	if (reduced_lp.row_lower_[iRow] != reduced_lp.row_upper_[iRow]) continue;
+	num_slack++;
+	printf("Column %d is slack\n", int(iCol));
 	double coeff = std::fabs(reduced_lp.a_matrix_.value_[reduced_lp.a_matrix_.start_[iCol]]);
 	if (coeff == 1.0) num_unit_coeff_slack++;
 	min_slack_coeff = std::min(coeff,min_slack_coeff);
 	max_slack_coeff = std::max(coeff,max_slack_coeff);
 	found_slack[iRow] = true;
-	num_slack++;
 	if (reduced_lp.col_cost_[iCol] == 0) num_zero_cost_slack++;
       }
       printf("grepSlack,model,col,slack,unit coeff,zero_cost,double,min coeff, max_coeff\n");
       printf("grepSlack,%s,%d, %d, %d, %d, %d, %g, %g\n", this->model_.lp_.model_name_.c_str(),
 	     int(reduced_lp.num_col_), int(num_slack), int(num_unit_coeff_slack), int(num_zero_cost_slack), int(num_double_slack),
 	     min_slack_coeff, max_slack_coeff);
-    }
+      //    }
 
     switch (model_presolve_status_) {
       case HighsPresolveStatus::kNotPresolved: {

src/presolve/HPresolve.cpp

@@ -4360,6 +4360,11 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) {
       break;
     }
 
+    // Now consider removing slacks
+    if (options->presolve_remove_slacks) {
+      HPRESOLVE_CHECKED_CALL(removeSlacks(postsolve_stack));
+    }
+
     report();
   } else {
     highsLogUser(options->log_options, HighsLogType::kInfo,
@@ -4375,6 +4380,62 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) {
   return Result::kOk;
 }
 
+HPresolve::Result HPresolve::removeSlacks(HighsPostsolveStack& postsolve_stack) {
+  // singletonColumns data structure appears not to be retained
+  // throughout presolve
+  //
+  bool unit_coeff_only = true;
+  for (HighsInt iCol = 0; iCol != model->num_col_; ++iCol) {
+    if (colDeleted[iCol]) continue;
+    if (colsize[iCol] != 1) continue;
+    if (model->integrality_[iCol] == HighsVarType::kInteger) continue;
+    HighsInt coliter = colhead[iCol];
+    HighsInt iRow = Arow[coliter];
+    assert(Acol[coliter] == iCol);
+    assert(!rowDeleted[iRow]);
+    if (model->row_lower_[iRow] != model->row_upper_[iRow]) continue;
+    double lower = model->col_lower_[iCol];
+    double upper = model->col_upper_[iCol];
+    double cost = model->col_cost_[iCol];
+    double rhs = model->row_lower_[iRow];
+    double coeff = Avalue[coliter];
+    printf("Col %d is continuous and is singleton in equality row %d with cost %g, bounds [%g, %g], coeff %g and RHS = %g\n", int(iCol), int(iRow), cost, lower, upper, coeff, rhs);
+    if (unit_coeff_only && std::fabs(coeff) != 1.0) continue;
+    assert(coeff);
+    // Slack is s = (rhs - a^Tx)/coeff
+    //
+    if (coeff > 0) {
+      // Constraint bounds become [rhs - coeff * lower, rhs - coeff *
+      // upper]
+      model->col_lower_[iCol] = rhs - coeff * lower;
+      model->col_upper_[iCol] = rhs - coeff * upper;
+    } else {
+      // Constraint bounds become [rhs - coeff * upper, rhs - coeff *
+      // lower]
+      model->col_lower_[iCol] = rhs - coeff * upper;
+      model->col_upper_[iCol] = rhs - coeff * lower;
+    }
+    if (cost) {
+      // Cost is (cost * rhs / coeff) + (col_cost - (cost/coeff) row_values)^Tx
+      double multiplier = cost / coeff;
+      for (const HighsSliceNonzero& nonzero : getRowVector(iRow)) {
+	HighsInt local_iCol = nonzero.index();
+	double local_value = nonzero.value();
+	model->col_cost_[local_iCol] -= multiplier * local_value;
+      }
+      model->offset_ += multiplier * rhs;
+    }
+    // 
+    postsolve_stack.slackColSubstitution(iRow, iCol, rhs, cost, lower, upper, coeff,
+                                        getColumnVector(iCol));
+    markColDeleted(iCol);
+
+    unlink(coliter);
+    
+  }
+  return Result::kOk;
+}
+
 HPresolve::Result HPresolve::checkLimits(HighsPostsolveStack& postsolve_stack) {
   size_t numreductions = postsolve_stack.numReductions();
 

src/presolve/HPresolve.h

@@ -269,6 +269,8 @@ class HPresolve {
 
   Result presolve(HighsPostsolveStack& postsolve_stack);
 
+  Result removeSlacks(HighsPostsolveStack& postsolve_stack);
+  
   Result checkLimits(HighsPostsolveStack& postsolve_stack);
 
   void storeCurrentProblemSize();

src/presolve/HighsPostsolveStack.cpp

@@ -1355,6 +1355,11 @@ void HighsPostsolveStack::SlackColSubstitution::undo(
     const HighsOptions& options, const std::vector<Nonzero>& rowValues,
     const std::vector<Nonzero>& colValues, HighsSolution& solution,
     HighsBasis& basis) {
+
+    assert(111==222);
+    
+
+
   // a (removed) cut may have been used in this reduction.
   bool isModelRow = static_cast<size_t>(row) < solution.row_value.size();
 
@@ -1396,7 +1401,7 @@ void HighsPostsolveStack::SlackColSubstitution::undo(
 
   basis.col_status[col] = HighsBasisStatus::kBasic;
   if (isModelRow)
-    basis.row_status[row] = computeRowStatus(solution.row_dual[row], rowType);
+    basis.row_status[row] = computeRowStatus(solution.row_dual[row], RowType::kEq);
 }
 
 }  // namespace presolve

src/presolve/HighsPostsolveStack.h

@@ -222,9 +222,11 @@ class HighsPostsolveStack {
   struct SlackColSubstitution {
     double rhs;
     double colCost;
+    double colLower;
+    double colUpper;
+    double colCoeff;
     HighsInt row;
     HighsInt col;
-    RowType rowType;
 
     void undo(const HighsOptions& options,
               const std::vector<Nonzero>& rowValues,
@@ -337,21 +339,16 @@ class HighsPostsolveStack {
     reductionAdded(ReductionType::kFreeColSubstitution);
   }
 
-  template <typename RowStorageFormat, typename ColStorageFormat>
+  template <typename ColStorageFormat>
   void slackColSubstitution(HighsInt row, HighsInt col, double rhs,
-                           double colCost, RowType rowType,
-                           const HighsMatrixSlice<RowStorageFormat>& rowVec,
-                           const HighsMatrixSlice<ColStorageFormat>& colVec) {
-    rowValues.clear();
-    for (const HighsSliceNonzero& rowVal : rowVec)
-      rowValues.emplace_back(origColIndex[rowVal.index()], rowVal.value());
-
+			    double colCost, double colLower, double colUpper, double colCoeff, 
+			    const HighsMatrixSlice<ColStorageFormat>& colVec) {
     colValues.clear();
     for (const HighsSliceNonzero& colVal : colVec)
       colValues.emplace_back(origRowIndex[colVal.index()], colVal.value());
 
-    reductionValues.push(SlackColSubstitution{rhs, colCost, origRowIndex[row],
-                                             origColIndex[col], rowType});
+    reductionValues.push(SlackColSubstitution{rhs, colCost, colLower, colUpper, colCoeff, origRowIndex[row],
+                                             origColIndex[col]});
     reductionValues.push(rowValues);
     reductionValues.push(colValues);
     reductionAdded(ReductionType::kSlackColSubstitution);