Fourier-Motzkin presolve

check/TestCAPI.c

@@ -726,17 +726,17 @@ void testNames() {
       printf("Row    %" HIGHSINT_FORMAT " has name %s\n", iRow, name);
   }
 
-  // Check extraction of names for the presolved LP, in which the
-  // first row is removed
+  // Check extraction of names for the presolved LP
   Highs_presolve(highs);
   if (dev_run) Highs_writePresolvedModel(highs, "");
 
   HighsInt presolved_num_col = Highs_getPresolvedNumCol(highs);
   HighsInt presolved_num_row = Highs_getPresolvedNumRow(highs);
-  assert(presolved_num_col == num_col);
-  assert(presolved_num_row == num_row - 1);
+  // Fourier-Motzkin presolve reduction may add columns/rows
+  //assert(presolved_num_col == num_col);
+  //assert(presolved_num_row == num_row-1);
 
-  char presolved_name[5];
+  char presolved_name[512];
 
   return_status = Highs_getPresolvedColName(highs, -1, presolved_name);
   assert(return_status == kHighsStatusError);
@@ -1511,7 +1511,10 @@ void passPresolveGetLp() {
     double* presolved_row_upper =
         (double*)malloc(sizeof(double) * presolved_num_row);
     HighsInt* presolved_a_start =
-        (HighsInt*)malloc(sizeof(HighsInt) * (presolved_num_col + 1));
+        (HighsInt*)malloc(sizeof(HighsInt) *
+                          (presolved_a_format == kHighsMatrixFormatColwise
+                               ? presolved_num_col + 1
+                               : presolved_num_row + 1));
     HighsInt* presolved_a_index =
         (HighsInt*)malloc(sizeof(HighsInt) * presolved_num_nz);
     double* presolved_a_value =
@@ -1538,9 +1541,9 @@ void passPresolveGetLp() {
     assert(return_status == kHighsStatusOk);
     return_status = Highs_run(local_highs);
 
-    double* col_value = (double*)malloc(sizeof(double) * num_col);
-    double* col_dual = (double*)malloc(sizeof(double) * num_col);
-    double* row_dual = (double*)malloc(sizeof(double) * num_row);
+    double* col_value = (double*)malloc(sizeof(double) * presolved_num_col);
+    double* col_dual = (double*)malloc(sizeof(double) * presolved_num_col);
+    double* row_dual = (double*)malloc(sizeof(double) * presolved_num_row);
 
     return_status =
         Highs_getSolution(local_highs, col_value, col_dual, NULL, row_dual);

check/TestPresolve.cpp

@@ -71,6 +71,8 @@ TEST_CASE("postsolve-no-basis", "[highs_test_presolve]") {
           "Col      Primal  Col      Primal\n");
     for (HighsInt iCol = 0; iCol < presolved_lp.num_col_; iCol++) {
       HighsInt original_iCol = original_col_indices[iCol];
+      // Skip columns added by presolve (e.g. FME objective reformulation)
+      if (original_iCol >= highs.getNumCol()) continue;
       if (dev_run)
         printf("%3d %11.5g  %3d %11.5g\n", int(iCol), solution.col_value[iCol],
                int(original_iCol), postsolve_solution.col_value[original_iCol]);
@@ -136,11 +138,14 @@ TEST_CASE("presolve", "[highs_test_presolve]") {
   // Have to set matrix dimensions to match presolved_model.lp_
   lp.setMatrixDimensions();
   highs.passModel(lp);
+  // Disable Fourier-Motzkin so this LP is not reduced
+  highs.setOptionValue("presolve_rule_off", 1 << kPresolveRuleFourierMotzkin);
   REQUIRE(highs.presolve() == HighsStatus::kOk);
   REQUIRE(lp.equalButForNames(presolved_model.lp_));
   REQUIRE(highs.getModelPresolveStatus() == HighsPresolveStatus::kNotReduced);
   REQUIRE(highs.getModelStatus() == HighsModelStatus::kNotset);
   REQUIRE(!presolved_model.isEmpty());
+  highs.setOptionValue("presolve_rule_off", 0);
 
   special_lps.primalDualInfeasible1Lp(lp, require_model_status);
   highs.passModel(lp);

check/TestPresolveRules.cpp

@@ -4,7 +4,7 @@
 #include "Highs.h"
 #include "catch.hpp"
 
-const bool dev_run = false;
+const bool dev_run = true;
 
 void presolveOffOn(const std::string& message, const HighsLp& lp, Highs& h,
                    const HighsInt require_presolved_model_num_col = 0,
@@ -79,6 +79,67 @@ TEST_CASE("test-col-stuffing", "[highs_test_presolve_rules]") {
   h.resetGlobalScheduler(true);
 }
 
+TEST_CASE("test-fourier-motzkin", "[highs_test_presolve_rules]") {
+  Highs h;
+  h.setOptionValue("output_flag", dev_run);
+  h.setOptionValue("presolve_rule_test", kPresolveRuleFourierMotzkin);
+  h.setOptionValue("presolve_rule_logging", true);
+  h.setOptionValue("log_dev_level", 1);
+
+  const bool lp0 = true;
+  const bool lp1 = true;  // Makes eliminations marginal, and leaves x2=0
+  const bool lp2 = true;
+
+  // From "A novel linear optimization presolve technique based on
+  // Fourier-Motzkin elimination", Zhang, Ploskas and Sahinidis,
+  // Mathematical Programming Computation (2026) 18:345–378
+  HighsLp lp;
+
+  lp.num_col_ = 4;
+  lp.num_row_ = 3;
+
+  lp.col_cost_.assign(lp.num_col_, 0);
+  lp.col_lower_.assign(lp.num_col_, 0);
+  lp.col_upper_.assign(lp.num_col_, kHighsInf);
+  lp.col_upper_[0] = 40.0;
+
+  lp.row_lower_.assign(lp.num_row_, -kHighsInf);
+  lp.row_upper_ = {-30, 50, 40};
+  lp.a_matrix_.format_ = MatrixFormat::kRowwise;
+  lp.a_matrix_.start_ = {0, 3, 6, 9};
+  lp.a_matrix_.index_ = {0, 1, 3, 1, 2, 3, 1, 2, 3};
+  lp.a_matrix_.value_ = {-1, 1, -1, 2, 1, 2, 3, -1, 3};
+
+  if (lp0) {
+    REQUIRE(h.passModel(lp) == HighsStatus::kOk);
+    presolveOffOn("FM example from paper", lp, h);
+  }
+
+  lp.col_upper_[0] = 5.0;
+  lp.row_upper_ = {-30, 75, 50};
+
+  if (lp1) {
+    REQUIRE(h.passModel(lp) == HighsStatus::kOk);
+    presolveOffOn("FM example from paper - tightened", lp, h);
+  }
+
+  lp.col_cost_ = {1, 2, 3, 4};
+
+  REQUIRE(h.passModel(lp) == HighsStatus::kOk);
+
+  if (lp2) {
+    HighsInt require_presolved_model_num_col = 1;
+    HighsInt require_presolved_model_num_row = 6;
+    HighsInt require_presolved_model_num_nz = 6;
+    presolveOffOn("FM example from paper - tightened and with costs", lp, h,
+                  require_presolved_model_num_col,
+                  require_presolved_model_num_row,
+                  require_presolved_model_num_nz);
+  }
+
+  h.resetGlobalScheduler(true);
+}
+
 void presolveOffOn(const std::string& message, const HighsLp& lp, Highs& h,
                    const HighsInt require_presolved_model_num_col,
                    const HighsInt require_presolved_model_num_row,
@@ -89,7 +150,7 @@ void presolveOffOn(const std::string& message, const HighsLp& lp, Highs& h,
   // solvers cannot be tested
   const bool reduce_to_empty = require_presolved_model_num_col == 0 &&
                                require_presolved_model_num_row == 0;
-  const HighsInt to_k = reduce_to_empty ? 2 : 5;
+  const HighsInt to_k = reduce_to_empty ? 2 : 4;
   for (int k = 0; k < to_k; k++) {
     std::string solver = kSimplexString;
     std::string run_crossover = kHighsOnString;
@@ -108,9 +169,6 @@ void presolveOffOn(const std::string& message, const HighsLp& lp, Highs& h,
         solver = kIpmString;
         run_crossover = kHighsOffString;
         basis_postsolve = false;
-      } else {
-        solver = kHiPdlpString;
-        basis_postsolve = false;
       }
     }
     std::string presolve = presolve_on ? kHighsOnString : kHighsOffString;
@@ -138,7 +196,7 @@ void presolveOffOn(const std::string& message, const HighsLp& lp, Highs& h,
       REQUIRE(h.getModelStatus() == HighsModelStatus::kOptimal);
       REQUIRE(h.getInfo().num_primal_infeasibilities == 0);
       REQUIRE(h.getInfo().num_dual_infeasibilities == 0);
-      REQUIRE(h.getInfo().simplex_iteration_count == 0);
+      if (reduce_to_empty) REQUIRE(h.getInfo().simplex_iteration_count == 0);
       // Ensure that any basis postsolve is correct
       if (basis_postsolve)
         REQUIRE(run_data.num_simplex_iterations_after_postsolve == 0);

check/TestSemiVariables.cpp

@@ -335,6 +335,9 @@ TEST_CASE("3015", "[highs_test_semi_variables]") {
   double optimal_objective_value = -1407973.679417;
   Highs highs;
   highs.setOptionValue("output_flag", dev_run);
+  // Disable Fourier-Motzkin presolve so that the semi-variable
+  // infeasibility is still triggered with default mip_feasibility_tolerance
+  highs.setOptionValue("presolve_rule_off", 1 << kPresolveRuleFourierMotzkin);
   highs.readModel(filename);
   HighsStatus status = highs.run();
   REQUIRE(status == HighsStatus::kError);

highs/Highs.h

@@ -535,15 +535,15 @@ class Highs {
    * the presolved model
    */
   const HighsInt* getPresolveOrigColsIndex() const {
-    return presolve_.data_.postSolveStack.getOrigColsIndex();
+    return presolve_.data_.postSolveStack.getOrigColIndex().data();
   }
 
   /**
    * @brief Return a const pointer to the original row indices for the
    * presolved model
    */
   const HighsInt* getPresolveOrigRowsIndex() const {
-    return presolve_.data_.postSolveStack.getOrigRowsIndex();
+    return presolve_.data_.postSolveStack.getOrigRowIndex().data();
   }
 
   /**

highs/lp_data/HConst.h

@@ -282,7 +282,8 @@ enum PresolveRuleType : int {
   kPresolveRuleDualFixing,
   kPresolveRuleColStuffing,
   kPresolveRuleInitialSweep,
-  kPresolveRuleMax = kPresolveRuleInitialSweep,
+  kPresolveRuleFourierMotzkin,
+  kPresolveRuleMax = kPresolveRuleFourierMotzkin,
   kPresolveRuleLastAllowOff = kPresolveRuleMax,
   kPresolveRuleCount
 };

highs/lp_data/HStruct.h

@@ -30,9 +30,6 @@ struct HighsSolution {
   void clear();
   void print(const std::string& prefix = "",
              const std::string& message = "") const;
-  bool isModelRow(HighsInt row) const {
-    return static_cast<size_t>(row) < row_value.size();
-  }
 };
 
 struct HighsObjectiveSolution {

highs/lp_data/Highs.cpp

@@ -4002,7 +4002,7 @@ HighsPostsolveStatus Highs::runPostsolve() {
   const HighsInt report_3040_col = -21792;
   presolve_.data_.postSolveStack.undo(
       options_, presolve_.data_.recovered_solution_,
-      presolve_.data_.recovered_basis_, report_3040_col);
+      presolve_.data_.recovered_basis_, 0, report_3040_col);
   // Compute the row activities
   assert(model_.lp_.a_matrix_.isColwise());
   calculateRowValuesQuad(model_.lp_, presolve_.data_.recovered_solution_);

highs/lp_data/HighsLp.cpp

@@ -226,6 +226,7 @@ void HighsLp::clear() {
   this->is_moved_ = false;
   this->cost_row_location_ = -1;
   this->has_infinite_cost_ = false;
+  this->fme_obj_col_ = -1;
   this->mods_.clear();
 }
 

highs/lp_data/HighsLp.h

@@ -55,6 +55,7 @@ class HighsLp {
   bool is_moved_;
   HighsInt cost_row_location_;
   bool has_infinite_cost_;
+  HighsInt fme_obj_col_ = -1;
   HighsLpMods mods_;
 
   bool operator==(const HighsLp& lp) const;

highs/lp_data/HighsModelUtils.cpp

@@ -1521,6 +1521,8 @@ std::string utilPresolveRuleTypeToString(const HighsInt rule_type) {
     return "Col stuffing";
   } else if (rule_type == kPresolveRuleInitialSweep) {
     return "Initial sweep";
+  } else if (rule_type == kPresolveRuleFourierMotzkin) {
+    return "Fourier-Motzkin";
   }
   assert(1 == 0);
   return "????";

highs/lp_data/HighsOptions.h

@@ -461,6 +461,7 @@ struct HighsOptionsStruct {
   HighsInt presolve_substitution_maxfillin;
   HighsInt presolve_rule_off;
   HighsInt presolve_rule_test;
+  HighsInt presolve_fm_level;
   bool presolve_rule_logging;
   bool presolve_remove_slacks;
   bool no_unnecessary_rebuild_refactor;
@@ -636,6 +637,7 @@ struct HighsOptionsStruct {
         presolve_substitution_maxfillin(0),
         presolve_rule_off(0),
         presolve_rule_test(0),
+        presolve_fm_level(0),
         presolve_rule_logging(false),
         presolve_remove_slacks(false),
         no_unnecessary_rebuild_refactor(false),
@@ -1684,6 +1686,11 @@ class HighsOptions : public HighsOptionsStruct {
         &presolve_rule_test, 0, 0, kPresolveRuleMax);
     records.push_back(record_int);
 
+    record_int = new OptionRecordInt("presolve_fm_level",
+                                     "Fourier-Motzkin elimination level",
+                                     advanced, &presolve_fm_level, 0, 1, 1);
+    records.push_back(record_int);
+
     record_bool = new OptionRecordBool(
         "presolve_rule_logging", "Log effectiveness of presolve rules for LP",
         advanced, &presolve_rule_logging, false);

highs/mip/HighsCliqueTable.cpp

@@ -1283,9 +1283,10 @@ void HighsCliqueTable::extractCliques(HighsMipSolver& mipsolver,
     HighsInt start = mipsolver.mipdata_->ARstart_[i];
     HighsInt end = mipsolver.mipdata_->ARstart_[i + 1];
 
-    if (mipsolver.mipdata_->postSolveStack.getOrigRowIndex(i) >=
-        mipsolver.orig_model_->num_row_)
-      break;
+    if (mipsolver.mipdata_->postSolveStack.isCutRow(i)) {
+      if (!mipsolver.mipdata_->postSolveStack.hasAppendedRows()) break;
+      continue;
+    }
 
     // catch set packing and partitioning constraints that already have the form
     // of a clique without transformations and add those cliques with the rows

highs/mip/HighsCliqueTable.h

@@ -164,11 +164,7 @@ class HighsCliqueTable {
   int64_t numNeighbourhoodQueries;
 
   HighsCliqueTable(HighsInt ncols) {
-    invertedHashList.resize(2 * static_cast<size_t>(ncols));
-    invertedHashListSizeTwo.resize(2 * static_cast<size_t>(ncols));
-    numcliquesvar.resize(2 * static_cast<size_t>(ncols), 0);
-    colsubstituted.resize(ncols);
-    colDeleted.resize(ncols, false);
+    resize(static_cast<size_t>(ncols));
     nfixings = 0;
     numNeighbourhoodQueries = 0;
     numEntries = 0;
@@ -178,6 +174,14 @@ class HighsCliqueTable {
     allowParallel = true;
   }
 
+  void resize(size_t ncols) {
+    invertedHashList.resize(2 * ncols);
+    invertedHashListSizeTwo.resize(2 * ncols);
+    numcliquesvar.resize(2 * ncols, 0);
+    colsubstituted.resize(ncols);
+    colDeleted.resize(ncols, false);
+  }
+
   void setPresolveFlag(bool inPresolve) { this->inPresolve = inPresolve; }
 
   bool getPresolveFlag() const { return inPresolve; }

highs/mip/HighsImplications.h

@@ -58,15 +58,10 @@ class HighsImplications {
   std::vector<HighsSubstitution> substitutions;
   std::vector<HighsBool> colsubstituted;
   HighsImplications(const HighsMipSolver& mipsolver) : mipsolver(mipsolver) {
-    HighsInt numcol = mipsolver.numCol();
-    implications.resize(2 * static_cast<size_t>(numcol));
-    colsubstituted.resize(numcol);
-    vubs.resize(numcol);
-    vlbs.resize(numcol);
     nextCleanupCall = mipsolver.numNonzero();
     numImplications = 0;
     numVarBounds = 0;
-    maxVarBounds = calcMaxVarBounds(numcol);
+    resize(mipsolver.numCol());
   }
 
   std::function<void(HighsInt, HighsInt, HighsInt, double)>
@@ -78,22 +73,24 @@ class HighsImplications {
     implications.clear();
     implications.shrink_to_fit();
 
-    HighsInt numcol = mipsolver.numCol();
-    implications.resize(2 * static_cast<size_t>(numcol));
-    colsubstituted.resize(numcol);
     numImplications = 0;
     vubs.clear();
     vubs.shrink_to_fit();
-    vubs.resize(numcol);
     vlbs.clear();
     vlbs.shrink_to_fit();
-    vlbs.resize(numcol);
+    resize(mipsolver.numCol());
     numVarBounds = 0;
-    maxVarBounds = calcMaxVarBounds(numcol);
-
     nextCleanupCall = mipsolver.numNonzero();
   }
 
+  void resize(HighsInt ncols) {
+    implications.resize(2 * static_cast<size_t>(ncols));
+    colsubstituted.resize(ncols);
+    vubs.resize(ncols);
+    vlbs.resize(ncols);
+    maxVarBounds = calcMaxVarBounds(ncols);
+  }
+
   constexpr static int64_t calcMaxVarBounds(HighsInt numcol) {
     return int64_t{5000000} + 10 * static_cast<int64_t>(numcol);
   };