diff --git a/check/TestCAPI.c b/check/TestCAPI.c
index cad521f9c0a..6c9b3f5f6c5 100644
--- a/check/TestCAPI.c
+++ b/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);
diff --git a/check/TestPresolve.cpp b/check/TestPresolve.cpp
index 1485df41f59..624fb9ed440 100644
--- a/check/TestPresolve.cpp
+++ b/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);
diff --git a/check/TestPresolveRules.cpp b/check/TestPresolveRules.cpp
index 6a1717b72fe..2e49aefe515 100644
--- a/check/TestPresolveRules.cpp
+++ b/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);
diff --git a/check/TestSemiVariables.cpp b/check/TestSemiVariables.cpp
index 89cf9679fc2..0443d7eee4f 100644
--- a/check/TestSemiVariables.cpp
+++ b/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);
diff --git a/highs/Highs.h b/highs/Highs.h
index d777745db44..f6c72607b71 100644
--- a/highs/Highs.h
+++ b/highs/Highs.h
@@ -535,7 +535,7 @@ class Highs {
    * the presolved model
    */
   const HighsInt* getPresolveOrigColsIndex() const {
-    return presolve_.data_.postSolveStack.getOrigColsIndex();
+    return presolve_.data_.postSolveStack.getOrigColIndex().data();
   }
 
   /**
@@ -543,7 +543,7 @@ class Highs {
    * presolved model
    */
   const HighsInt* getPresolveOrigRowsIndex() const {
-    return presolve_.data_.postSolveStack.getOrigRowsIndex();
+    return presolve_.data_.postSolveStack.getOrigRowIndex().data();
   }
 
   /**
diff --git a/highs/lp_data/HConst.h b/highs/lp_data/HConst.h
index e150f7838a3..4449a734593 100644
--- a/highs/lp_data/HConst.h
+++ b/highs/lp_data/HConst.h
@@ -282,7 +282,8 @@ enum PresolveRuleType : int {
   kPresolveRuleDualFixing,
   kPresolveRuleColStuffing,
   kPresolveRuleInitialSweep,
-  kPresolveRuleMax = kPresolveRuleInitialSweep,
+  kPresolveRuleFourierMotzkin,
+  kPresolveRuleMax = kPresolveRuleFourierMotzkin,
   kPresolveRuleLastAllowOff = kPresolveRuleMax,
   kPresolveRuleCount
 };
diff --git a/highs/lp_data/HStruct.h b/highs/lp_data/HStruct.h
index 0f0081c5a0e..bd6bd74ceca 100644
--- a/highs/lp_data/HStruct.h
+++ b/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 {
diff --git a/highs/lp_data/Highs.cpp b/highs/lp_data/Highs.cpp
index fd912246ce2..b6010f9edd8 100644
--- a/highs/lp_data/Highs.cpp
+++ b/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_);
diff --git a/highs/lp_data/HighsLp.cpp b/highs/lp_data/HighsLp.cpp
index ceb4de3e51c..aef964ac2c0 100644
--- a/highs/lp_data/HighsLp.cpp
+++ b/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();
 }
 
diff --git a/highs/lp_data/HighsLp.h b/highs/lp_data/HighsLp.h
index 77236f406f2..a9200075d04 100644
--- a/highs/lp_data/HighsLp.h
+++ b/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;
diff --git a/highs/lp_data/HighsModelUtils.cpp b/highs/lp_data/HighsModelUtils.cpp
index 1d7f0b9a057..fc2f0d9d1b2 100644
--- a/highs/lp_data/HighsModelUtils.cpp
+++ b/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 "????";
diff --git a/highs/lp_data/HighsOptions.h b/highs/lp_data/HighsOptions.h
index 063918be135..10c21337587 100644
--- a/highs/lp_data/HighsOptions.h
+++ b/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);
diff --git a/highs/mip/HighsCliqueTable.cpp b/highs/mip/HighsCliqueTable.cpp
index 40a90902dc2..f44e1b32fdc 100644
--- a/highs/mip/HighsCliqueTable.cpp
+++ b/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
diff --git a/highs/mip/HighsCliqueTable.h b/highs/mip/HighsCliqueTable.h
index afaaedaae06..c9f27e81c9e 100644
--- a/highs/mip/HighsCliqueTable.h
+++ b/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; }
diff --git a/highs/mip/HighsImplications.h b/highs/mip/HighsImplications.h
index 50fa737395e..d574685de84 100644
--- a/highs/mip/HighsImplications.h
+++ b/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);
   };
diff --git a/highs/mip/HighsMipSolverData.cpp b/highs/mip/HighsMipSolverData.cpp
index 5dfd43feb2f..54ac726d0a6 100644
--- a/highs/mip/HighsMipSolverData.cpp
+++ b/highs/mip/HighsMipSolverData.cpp
@@ -1345,7 +1345,7 @@ void HighsMipSolverData::performRestart() {
   HighsInt numLpRows = getLp().getLp().num_row_;
   HighsInt numModelRows = mipsolver.numRow();
   HighsInt numCuts = numLpRows - numModelRows;
-  if (numCuts > 0) postSolveStack.appendCutsToModel(numCuts);
+  postSolveStack.appendCutsToModel(numCuts);
   auto integrality = std::move(presolvedModel.integrality_);
   double offset = presolvedModel.offset_;
   presolvedModel = getLp().getLp();
@@ -1361,13 +1361,14 @@ void HighsMipSolverData::performRestart() {
     // if we have a basis after solving the root LP, we expand it to the
     // original space so that it can be used for constructing a starting basis
     // for the presolved model after the restart
-    root_basis.col_status.resize(postSolveStack.getOrigNumCol());
-    root_basis.row_status.resize(postSolveStack.getOrigNumRow(),
+    root_basis.col_status.resize(postSolveStack.getNextColIndex());
+    root_basis.row_status.resize(postSolveStack.getNextRowIndex(),
                                  HighsBasisStatus::kBasic);
     root_basis.valid = true;
     root_basis.useful = true;
 
-    for (HighsInt i = 0; i < mipsolver.numCol(); ++i)
+    HighsInt numCol = basis.col_status.size();
+    for (HighsInt i = 0; i < numCol; ++i)
       root_basis.col_status[postSolveStack.getOrigColIndex(i)] =
           basis.col_status[i];
 
@@ -1492,13 +1493,19 @@ void HighsMipSolverData::basisTransfer() {
     firstrootbasis.alien = true;
     firstrootbasis.useful = true;
 
-    for (HighsInt i = 0; i < numRow; ++i) {
+    for (HighsInt i = 0;
+         i < static_cast<HighsInt>(postSolveStack.getOrigRowIndex().size());
+         ++i) {
+      if (!postSolveStack.isOrigRow(i)) break;
       HighsBasisStatus status =
           mipsolver.rootbasis->row_status[postSolveStack.getOrigRowIndex(i)];
       firstrootbasis.row_status[i] = status;
     }
 
-    for (HighsInt i = 0; i < numCol; ++i) {
+    for (HighsInt i = 0;
+         i < static_cast<HighsInt>(postSolveStack.getOrigColIndex().size());
+         ++i) {
+      if (!postSolveStack.isOrigCol(i)) break;
       HighsBasisStatus status =
           mipsolver.rootbasis->col_status[postSolveStack.getOrigColIndex(i)];
       firstrootbasis.col_status[i] = status;
diff --git a/highs/mip/HighsPseudocost.cpp b/highs/mip/HighsPseudocost.cpp
index 73d7eeded3f..6c670c19e79 100644
--- a/highs/mip/HighsPseudocost.cpp
+++ b/highs/mip/HighsPseudocost.cpp
@@ -47,6 +47,7 @@ HighsPseudocost::HighsPseudocost(const HighsMipSolver& mipsolver)
         mipsolver.pscostinit->conflict_avg_score * mipsolver.numCol();
 
     for (HighsInt i = 0; i != mipsolver.numCol(); ++i) {
+      if (!mipsolver.mipdata_->postSolveStack.isOrigCol(i)) continue;
       HighsInt origCol = mipsolver.mipdata_->postSolveStack.getOrigColIndex(i);
 
       pseudocostup[i] = mipsolver.pscostinit->pseudocostup[origCol];
@@ -115,6 +116,7 @@ HighsPseudocostInitialization::HighsPseudocostInitialization(
   conflict_avg_score /= ncols * pscost.conflict_weight;
 
   for (HighsInt i = 0; i != ncols; ++i) {
+    if (!postsolveStack.isOrigCol(i)) continue;
     pseudocostup[postsolveStack.getOrigColIndex(i)] = pscost.pseudocostup[i];
     pseudocostdown[postsolveStack.getOrigColIndex(i)] =
         pscost.pseudocostdown[i];
diff --git a/highs/presolve/HPresolve.cpp b/highs/presolve/HPresolve.cpp
index 1989287b47f..3a3d6624c02 100644
--- a/highs/presolve/HPresolve.cpp
+++ b/highs/presolve/HPresolve.cpp
@@ -499,6 +499,7 @@ void HPresolve::chooseRules() {
     presolve_light_rule_off[kPresolveRuleEnumeration] = true;
     presolve_light_rule_off[kPresolveRuleDualFixing] = true;
     presolve_light_rule_off[kPresolveRuleColStuffing] = true;
+    presolve_light_rule_off[kPresolveRuleFourierMotzkin] = true;
   }
 
   if (!silent && options->log_dev_level) {
@@ -807,8 +808,7 @@ HPresolve::Result HPresolve::updateColImpliedBounds(HighsInt row, HighsInt col,
     // row can be removed and should not be used, e.g., to identify a
     // column as implied free
     bool useImplBound = mipsolver == nullptr ||
-                        mipsolver->mipdata_->postSolveStack.getOrigRowIndex(
-                            row) < mipsolver->orig_model_->num_row_;
+                        !mipsolver->mipdata_->postSolveStack.isCutRow(row);
 
     if (direction * val > 0) {
       // upper bound
@@ -993,6 +993,8 @@ void HPresolve::shrinkProblem(HighsPostsolveStack& postsolve_stack) {
       }
     }
   }
+  if (model->fme_obj_col_ >= 0)
+    model->fme_obj_col_ = newColIndex[model->fme_obj_col_];
   colDeleted.assign(model->num_col_, false);
   model->col_cost_.resize(model->num_col_);
   model->col_lower_.resize(model->num_col_);
@@ -1165,8 +1167,6 @@ void HPresolve::shrinkProblem(HighsPostsolveStack& postsolve_stack) {
 
     mipsolver->mipdata_->debugSolution.shrink(newColIndex);
     numProbes.resize(model->num_col_);
-    // Need to set the constraint matrix dimensions
-    model->setMatrixDimensions();
   }
   // Need to set the constraint matrix dimensions
   model->setMatrixDimensions();
@@ -1461,12 +1461,8 @@ HPresolve::Result HPresolve::dominatedColumns(
         if (!tryToFix) numDomChecksPredBndAnalysis++;
         // check for domination
         if (checkDomination(direction, col, direction_k, k)) {
-          // Re-check the implied bound condition since earlier fixings in
-          // this dominatedColumns call may have changed the model state
-          bool currentBoundImplied =
-              direction > 0 ? isUpperImplied(col) : isLowerImplied(col);
           if (tryToFix &&
-              (currentBoundImplied ||
+              (boundImplied ||
                mipsolver->mipdata_->cliquetable.haveCommonClique(
                    HighsCliqueTable::CliqueVar(col, direction > 0 ? 1 : 0),
                    HighsCliqueTable::CliqueVar(k, direction_k > 0 ? 1 : 0)))) {
@@ -1498,9 +1494,8 @@ HPresolve::Result HPresolve::dominatedColumns(
 
     // lambda for finding a domination relationship in the given row
     auto checkRow = [&](HighsInt row, HighsInt col, HighsInt direction,
-                        double bestVal, bool boundImplied, bool hasCliques) {
+                        double bestVal, bool hasCliques) {
       storeRow(row);
-      bool onlyPredBndAnalysis = !boundImplied && !hasCliques;
       for (const HighsSliceNonzero& nonz : getStoredRow()) {
         // get column index
         HighsInt k = nonz.index();
@@ -1514,8 +1509,12 @@ HPresolve::Result HPresolve::dominatedColumns(
         // check if variables have the same type
         bool sameVarType = varsHaveSameType(col, k);
 
+        // check if bound is implied (computed fresh due to earlier fixings)
+        bool boundImplied =
+            direction > 0 ? isUpperImplied(col) : isLowerImplied(col);
+
         // skip checks if nothing to do
-        if (onlyPredBndAnalysis && !sameVarType) continue;
+        if (!boundImplied && !hasCliques && !sameVarType) continue;
 
         // try to fix variables or strengthen bounds
         // check already known non-zeros in respective columns in advance to
@@ -1550,15 +1549,13 @@ HPresolve::Result HPresolve::dominatedColumns(
     if (bestRowMinus != -1 &&
         (allowPredBndAnalysis || lowerImplied || hasNegCliques))
       HPRESOLVE_CHECKED_CALL(checkRow(bestRowMinus, j, HighsInt{-1},
-                                      ajBestRowMinus, lowerImplied,
-                                      hasNegCliques));
+                                      ajBestRowMinus, hasNegCliques));
 
     // use row 'bestRowPlus'
     if (!colDeleted[j] && bestRowPlus != -1 &&
         (allowPredBndAnalysis || upperImplied || hasPosCliques))
-      HPRESOLVE_CHECKED_CALL(checkRow(bestRowPlus, j, HighsInt{1},
-                                      ajBestRowPlus, upperImplied,
-                                      hasPosCliques));
+      HPRESOLVE_CHECKED_CALL(
+          checkRow(bestRowPlus, j, HighsInt{1}, ajBestRowPlus, hasPosCliques));
 
     // do not use predictive bound analysis if it requires many domination
     // checks and only yields few fixings or improved bounds on average
@@ -2264,6 +2261,101 @@ void HPresolve::addToMatrix(const HighsInt row, const HighsInt col,
   }
 }
 
+bool HPresolve::addToMatrix(
+    HighsPostsolveStack& postsolve_stack, const std::vector<double>& row_lower,
+    const std::vector<double>& row_upper,
+    const std::vector<std::vector<row_entry>>& row_entries) {
+  // update number of rows
+  HighsInt num_rows = static_cast<HighsInt>(row_entries.size());
+  if (num_rows == 0) return true;
+  HighsInt oldNumRows = model->num_row_;
+  model->num_row_ += num_rows;
+  model->a_matrix_.num_row_ += num_rows;
+
+  // resize postsolve vectors
+  postsolve_stack.appendRowsToModel(num_rows);
+
+  // add row bounds
+  model->row_lower_.insert(model->row_lower_.end(), row_lower.begin(),
+                           row_lower.end());
+  model->row_upper_.insert(model->row_upper_.end(), row_upper.begin(),
+                           row_upper.end());
+
+  // initialise row sizes
+  if (!okResize(rowroot, model->num_row_, HighsInt{-1})) return false;
+  if (!okResize(rowsize, model->num_row_, HighsInt{0})) return false;
+  if (!okResize(rowsizeInteger, model->num_row_, HighsInt{0})) return false;
+  if (!okResize(rowsizeImplInt, model->num_row_, HighsInt{0})) return false;
+
+  // initialise row duals
+  if (!okResize(rowDualLower, model->num_row_, -kHighsInf)) return false;
+  if (!okResize(rowDualUpper, model->num_row_, kHighsInf)) return false;
+  for (HighsInt i = oldNumRows; i < model->num_row_; i++) {
+    if (model->row_lower_[i] == -kHighsInf) rowDualUpper[i] = 0;
+    if (model->row_upper_[i] == kHighsInf) rowDualLower[i] = 0;
+  }
+
+  // initialise implied row duals
+  if (!okResize(implRowDualLower, model->num_row_, -kHighsInf)) return false;
+  if (!okResize(implRowDualUpper, model->num_row_, kHighsInf)) return false;
+  if (!okResize(rowDualLowerSource, model->num_row_, HighsInt{-1}))
+    return false;
+  if (!okResize(rowDualUpperSource, model->num_row_, HighsInt{-1}))
+    return false;
+  if (!okResize(colImplSourceByRow, model->num_row_, std::set<HighsInt>{}))
+    return false;
+
+  // initialise flags
+  if (!okResize(changedRowFlag, model->num_row_, uint8_t{0})) return false;
+  if (!okResize(rowDeleted, model->num_row_, uint8_t{0})) return false;
+  if (!okResize(singleEquationChecked, model->num_row_, uint8_t{0}))
+    return false;
+
+  // initialise row names
+  if (!okResize(model->row_names_, model->num_row_, std::string{}))
+    return false;
+
+  // resize vector for equations
+  if (!okResize(eqiters, model->num_row_, equations.end())) return false;
+
+  // resize vectors for implied row bounds
+  impliedRowBounds.setNumSums(model->num_row_);
+
+  // set bound arrays again (pointers may get invalidated by reallocation)
+  impliedDualRowBounds.setBoundArrays(
+      rowDualLower.data(), rowDualUpper.data(), implRowDualLower.data(),
+      implRowDualUpper.data(), rowDualLowerSource.data(),
+      rowDualUpperSource.data());
+
+  for (HighsInt i = 0; i < num_rows; i++) {
+    // new row index
+    HighsInt row = oldNumRows + i;
+
+    // add non-zeros
+    for (const auto& entry : row_entries[i])
+      addToMatrix(row, entry.col, entry.val);
+
+    // add row singleton
+    if (rowsize[row] == 1) singletonRows.push_back(row);
+
+    // add equation
+    if (isEquation(row))
+      eqiters[row] = equations.emplace(rowsize[row], row).first;
+  }
+
+  return true;
+}
+
+bool HPresolve::addToMatrix(HighsPostsolveStack& postsolve_stack,
+                            double row_lower, double row_upper,
+                            std::vector<row_entry> row_entries) {
+  std::vector<double> rl = {row_lower};
+  std::vector<double> ru = {row_upper};
+  std::vector<std::vector<row_entry>> re;
+  re.push_back(std::move(row_entries));
+  return addToMatrix(postsolve_stack, rl, ru, re);
+}
+
 HighsTripletListSlice HPresolve::getColumnVector(HighsInt col) const {
   return HighsTripletListSlice(Arow.data(), Avalue.data(), Anext.data(),
                                colhead[col]);
@@ -2305,6 +2397,7 @@ void HPresolve::markColDeleted(HighsInt col) {
     colDeleted[col] = true;
   }
   ++numDeletedCols;
+  if (col == model->fme_obj_col_) model->fme_obj_col_ = -1;
 }
 
 HPresolve::Result HPresolve::changeColUpper(HighsInt col, double newUpper) {
@@ -3808,11 +3901,11 @@ HPresolve::Result HPresolve::rowPresolve(HighsPostsolveStack& postsolve_stack,
   double origRowUpper = model->row_upper_[row];
   double origRowLower = model->row_lower_[row];
 
+  // Convert to equality constraint and record for dual postsolve
   if (!isEquation(row)) {
     if (isImpliedEquationAtLower(row)) {
-      // Convert to equality constraint (note that currently postsolve will not
-      // know about this conversion)
       model->row_upper_[row] = model->row_lower_[row];
+      postsolve_stack.impliedEquation(row, true, getRowVector(row));
       // Since row upper bound is now finite, lower bound on row dual is
       // -kHighsInf
       changeRowDualLower(row, -kHighsInf);
@@ -3820,9 +3913,8 @@ HPresolve::Result HPresolve::rowPresolve(HighsPostsolveStack& postsolve_stack,
         HPRESOLVE_CHECKED_CALL(
             checkRedundantBounds(rowDualLowerSource[row], row));
     } else if (isImpliedEquationAtUpper(row)) {
-      // Convert to equality constraint (note that currently postsolve will not
-      // know about this conversion)
       model->row_lower_[row] = model->row_upper_[row];
+      postsolve_stack.impliedEquation(row, false, getRowVector(row));
       // Since row lower bound is now finite, upper bound on row dual is
       // kHighsInf
       changeRowDualUpper(row, kHighsInf);
@@ -6143,6 +6235,8 @@ HPresolve::Result HPresolve::initialSweep(
   model->a_matrix_.start_.resize(num_col + 1);
   model->a_matrix_.index_.resize(nnz);
   model->a_matrix_.value_.resize(nnz);
+  if (model->fme_obj_col_ >= 0)
+    model->fme_obj_col_ = newColIndex[model->fme_obj_col_];
   postsolve_stack.compressColIndexMap(newColIndex);
   HPRESOLVE_CHECKED_CALL(checkLimits(postsolve_stack));
 
@@ -6528,6 +6622,8 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) {
         mipsolver != nullptr || !options->lp_presolve_requires_basis_postsolve;
 #endif
     bool tryProbing = mipsolver != nullptr;
+    bool tryFourierMotzkin = true;
+
     HighsInt numCliquesBeforeProbing = -1;
     bool domcolAfterProbingCalled = false;
     bool dependentEquationsCalled = mipsolver != nullptr;
@@ -6564,8 +6660,10 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) {
             applyConflictGraphSubstitutions(postsolve_stack, numDelCol));
       }
 
-      const bool reduced_to_empty = numDeletedCols == model->num_col_ &&
-                                    numDeletedRows == model->num_row_;
+      HighsInt numColsEliminatedFourierMotzkin = 0;
+      if (tryFourierMotzkin && this->allow_rule_[kPresolveRuleFourierMotzkin])
+        HPRESOLVE_CHECKED_CALL(
+            fourierMotzkin(postsolve_stack, numColsEliminatedFourierMotzkin));
 
       if (reducedToEmpty()) break;
 
@@ -6575,7 +6673,10 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) {
         analysis_.presolveTimerStop(kPresolveClockAggregator);
       }
 
-      if (problemSizeReduction() > 0.05) continue;
+      // check if there were reductions
+      bool haveReductions = problemSizeReduction() > 0.05;
+      tryFourierMotzkin = haveReductions || numColsEliminatedFourierMotzkin > 0;
+      if (haveReductions) continue;
 
       if (trySparsify && this->allow_rule_[kPresolveRuleSparsify]) {
         HighsInt numNz = numNonzeros();
@@ -7077,10 +7178,8 @@ HighsModelStatus HPresolve::run(HighsPostsolveStack& postsolve_stack) {
       cutvals.reserve(model->num_col_);
       HighsInt numcuts = 0;
       for (HighsInt i = model->num_row_ - 1; i >= 0; --i) {
-        // check if we already reached the original rows
-        if (postsolve_stack.getOrigRowIndex(i) <
-            mipsolver->orig_model_->num_row_)
-          break;
+        if (postsolve_stack.isOrigRow(i)) break;
+        if (!postsolve_stack.isCutRow(i)) continue;
 
         // row is a cut, remove it from matrix but add to cutpool
         ++numcuts;
@@ -7103,10 +7202,37 @@ HighsModelStatus HPresolve::run(HighsPostsolveStack& postsolve_stack) {
         for (HighsInt j : rowpositions) unlink(j);
       }
 
-      model->num_row_ -= numcuts;
-      model->row_lower_.resize(model->num_row_);
-      model->row_upper_.resize(model->num_row_);
-      model->row_names_.resize(model->num_row_);
+      // Compact deleted cut rows. shrinkProblem must not be used here
+      // because it replaces the cutpool with a new empty one, destroying
+      // the cuts that were just added above.
+      auto compactDeletedRows = [&]() {
+        HighsInt oldNumRow = model->num_row_;
+        std::vector<HighsInt> newRowIndex(oldNumRow);
+        HighsInt newNumRow = 0;
+        for (HighsInt i = 0; i < oldNumRow; ++i) {
+          if (rowDeleted[i])
+            newRowIndex[i] = -1;
+          else
+            newRowIndex[i] = newNumRow++;
+        }
+        model->num_row_ = newNumRow;
+
+        for (HighsInt i = 0; i < oldNumRow; ++i) {
+          if (newRowIndex[i] == -1 || newRowIndex[i] == i) continue;
+          model->row_lower_[newRowIndex[i]] = model->row_lower_[i];
+          model->row_upper_[newRowIndex[i]] = model->row_upper_[i];
+        }
+        model->row_lower_.resize(model->num_row_);
+        model->row_upper_.resize(model->num_row_);
+        model->row_names_.resize(model->num_row_);
+
+        for (size_t i = 0; i < Avalue.size(); ++i) {
+          if (Avalue[i] == 0) continue;
+          assert(newRowIndex[Arow[i]] != -1);
+          Arow[i] = newRowIndex[Arow[i]];
+        }
+      };
+      compactDeletedRows();
     }
   }
 
@@ -7547,6 +7673,852 @@ HPresolve::Result HPresolve::aggregator(HighsPostsolveStack& postsolve_stack) {
   return Result::kOk;
 }
 
+HPresolve::Result HPresolve::fourierMotzkin(
+    HighsPostsolveStack& postsolve_stack, HighsInt& numColsEliminated) {
+  assert(this->allow_rule_[kPresolveRuleFourierMotzkin]);
+  const bool logging_on = analysis_.logging_on_;
+  if (logging_on) analysis_.startPresolveRuleLog(kPresolveRuleFourierMotzkin);
+
+  using FmeRow = HighsPostsolveStack::FmeRowData<HighsTripletTreeSlicePreOrder>;
+  using FmeAncestryEntry = HighsPostsolveStack::FmeAncestryEntry;
+  using FmeBlockStep = HighsPostsolveStack::FmeBlockStep;
+
+  // max. absolute coefficient
+  const double maxCoef = 1e3;
+
+  // max. number of consecutive failures (while trying to build the heap)
+  const HighsInt maxNumFails = 100;
+  // max. size of the heap
+  const HighsInt maxHeapSize = 10000;
+
+  // sentinel row indices for variable bounds and objective row
+  const HighsInt kUpperBoundRow = -2;
+  const HighsInt kLowerBoundRow = -3;
+  const HighsInt kObjectiveRow = -4;
+
+  // structs
+  struct Heap {
+    struct Entry {
+      HighsInt col;
+      int64_t neRed;
+      int64_t mrRed;
+    };
+
+    std::vector<Entry> entries;
+    std::vector<HighsInt> pos;
+
+    bool empty() const { return entries.empty(); }
+    HighsInt size() const { return static_cast<HighsInt>(entries.size()); }
+    HighsInt top() const { return entries[0].col; }
+    bool contains(HighsInt col) const { return pos[col] != -1; }
+
+    void reset(HighsInt numCol, HighsInt reserveSize) {
+      entries.clear();
+      entries.reserve(reserveSize);
+      pos.assign(numCol, -1);
+    }
+
+    void push(HighsInt col, int64_t neRed, int64_t mrRed) {
+      pos[col] = size();
+      entries.push_back({col, neRed, mrRed});
+    }
+
+    void insert(HighsInt col, int64_t neRed, int64_t mrRed) {
+      push(col, neRed, mrRed);
+      siftUp(pos[col]);
+    }
+
+    void remove(HighsInt col) {
+      HighsInt p = pos[col];
+      if (p == -1) return;
+      swap(p, size() - 1);
+      pos[col] = -1;
+      entries.pop_back();
+      siftUp(p);
+      siftDown(p);
+    }
+
+    void update(HighsInt col, int64_t neRed, int64_t mrRed) {
+      HighsInt p = pos[col];
+      if (p == -1) return;
+      entries[p].neRed = neRed;
+      entries[p].mrRed = mrRed;
+      siftUp(p);
+      siftDown(p);
+    }
+
+    void heapify() {
+      for (HighsInt i = size() / 2 - 1; i >= 0; --i) siftDown(i);
+    }
+
+   private:
+    bool better(HighsInt i, HighsInt j) const {
+      if (entries[i].neRed != entries[j].neRed)
+        return entries[i].neRed > entries[j].neRed;
+      return entries[i].mrRed > entries[j].mrRed;
+    }
+
+    void swap(HighsInt i, HighsInt j) {
+      if (i == j) return;
+      std::swap(entries[i], entries[j]);
+      pos[entries[i].col] = i;
+      pos[entries[j].col] = j;
+    }
+
+    void siftUp(HighsInt i) {
+      if (i >= size()) return;
+      while (i > 0) {
+        HighsInt parent = (i - 1) / 2;
+        if (!better(i, parent)) break;
+        swap(i, parent);
+        i = parent;
+      }
+    }
+
+    void siftDown(HighsInt i) {
+      HighsInt n = size();
+      if (i >= n) return;
+      while (true) {
+        HighsInt best = i;
+        HighsInt left = 2 * i + 1;
+        HighsInt right = 2 * i + 2;
+        if (left < n && better(left, best)) best = left;
+        if (right < n && better(right, best)) best = right;
+        if (best == i) break;
+        swap(i, best);
+        i = best;
+      }
+    }
+  };
+
+  struct newRowEntry {
+    HighsInt col;
+    HighsCDouble val;
+  };
+
+  struct newRow {
+    std::vector<newRowEntry> entries;
+    double lower;
+    double upper;
+    HighsInt plusIndex;
+    HighsInt minusIndex;
+    double plusScale;
+    double minusScale;
+  };
+
+  struct NewRowOrigin {
+    HighsInt plusRow;
+    HighsInt minusRow;
+    double plusScale;
+    double minusScale;
+  };
+
+  auto finalise = [&]() {
+    analysis_.logging_on_ = logging_on;
+    if (logging_on) analysis_.stopPresolveRuleLog(kPresolveRuleFourierMotzkin);
+    return checkLimits(postsolve_stack);
+  };
+
+  auto acceptCoef = [&](double val) {
+    double absval = std::abs(val);
+    return absval == 0.0 || (absval >= 1.0 / maxCoef && absval <= maxCoef);
+  };
+
+  auto isCandidate = [&](HighsInt col) {
+    if (colDeleted[col]) return false;
+    if (colsize[col] == 0) return false;
+    if (col == model->fme_obj_col_) return false;
+    if (model->integrality_[col] != HighsVarType::kContinuous) return false;
+    if (!acceptCoef(model->col_cost_[col])) return false;
+    if (options->presolve_fm_level < 1 && model->col_cost_[col] != 0.0)
+      return false;
+    for (const auto& nz : getColumnVector(col))
+      if (isEquation(nz.index()) || !acceptCoef(nz.value())) return false;
+    return true;
+  };
+
+  auto computeCandidates = [&](std::vector<HighsInt>& candidates) {
+    candidates.clear();
+    for (HighsInt col = 0; col < model->num_col_; col++)
+      if (isCandidate(col)) candidates.push_back(col);
+    pdqsort(candidates.begin(), candidates.end(),
+            [&](HighsInt a, HighsInt b) { return colsize[a] < colsize[b]; });
+    return !candidates.empty();
+  };
+
+  auto checkRows = [&](HighsInt col, const std::vector<HighsInt>& objRowCols,
+                       std::vector<HighsInt>& iPlus,
+                       std::vector<HighsInt>& iMinus, int64_t& nePlus,
+                       int64_t& neMinus) {
+    nePlus = 0;
+    neMinus = 0;
+    iPlus.clear();
+    iMinus.clear();
+    for (const auto& nz : getColumnVector(col)) {
+      HighsInt row = nz.index();
+      if (rowDeleted[row]) continue;
+
+      if (isRanged(row)) {
+        iPlus.push_back(row);
+        nePlus += rowsize[row];
+        iMinus.push_back(row);
+        neMinus += rowsize[row];
+      } else {
+        HighsInt direction;
+        if (model->row_lower_[row] == -kHighsInf &&
+            model->row_upper_[row] != kHighsInf)
+          direction = 1;
+        else
+          direction = -1;
+
+        if (direction * nz.value() > 0) {
+          iPlus.push_back(row);
+          nePlus += rowsize[row];
+        } else {
+          iMinus.push_back(row);
+          neMinus += rowsize[row];
+        }
+      }
+    }
+
+    // include finite variable bounds as singleton rows
+    if (model->col_upper_[col] != kHighsInf) {
+      iPlus.push_back(kUpperBoundRow);
+      nePlus += 1;
+    }
+    if (model->col_lower_[col] != -kHighsInf) {
+      iMinus.push_back(kLowerBoundRow);
+      neMinus += 1;
+    }
+
+    // simulate the objective constraint row for candidates with nonzero
+    // cost when the reformulation has not yet been performed
+    if (!objRowCols.empty() && model->col_cost_[col] != 0.0) {
+      int64_t objRowSize = static_cast<int64_t>(objRowCols.size());
+      if (model->col_cost_[col] > 0.0) {
+        iPlus.push_back(kObjectiveRow);
+        nePlus += objRowSize;
+      } else {
+        iMinus.push_back(kObjectiveRow);
+        neMinus += objRowSize;
+      }
+    }
+  };
+
+  auto collectAffectedCols = [&](HighsInt col, const std::vector<HighsInt>& set,
+                                 const std::vector<HighsInt>& objRowCols,
+                                 std::vector<HighsInt>& mark,
+                                 std::vector<HighsInt>& otherMark,
+                                 std::vector<HighsInt>& affectedCols) {
+    for (HighsInt row : set) {
+      if (row == kObjectiveRow) {
+        for (HighsInt k : objRowCols) {
+          if (k == col) continue;
+          if (mark[k] == 0 && otherMark[k] == 0) affectedCols.push_back(k);
+          mark[k]++;
+        }
+      } else {
+        if (row < 0) continue;
+        for (const auto& nz : getRowVector(row)) {
+          HighsInt k = nz.index();
+          if (k == col) continue;
+          if (mark[k] == 0 && otherMark[k] == 0) affectedCols.push_back(k);
+          mark[k]++;
+        }
+      }
+    }
+  };
+
+  auto checkNonZeros = [&](HighsInt col,
+                           const std::vector<HighsInt>& objRowCols,
+                           std::vector<HighsInt>& iPlus,
+                           std::vector<HighsInt>& iMinus,
+                           std::vector<HighsInt>& pPlus,
+                           std::vector<HighsInt>& pMinus,
+                           std::vector<HighsInt>& affectedCols, int64_t& neRed,
+                           int64_t& mrRed) {
+    // initialise
+    neRed = 0;
+    mrRed = 0;
+
+    // check rows
+    int64_t nePlus;
+    int64_t neMinus;
+    checkRows(col, objRowCols, iPlus, iMinus, nePlus, neMinus);
+
+    if (iPlus.size() == 0 || iMinus.size() == 0) {
+      // other presolve reductions may handle this case (e.g., implied free
+      // column substitution)
+      iPlus.clear();
+      iMinus.clear();
+      return false;
+    }
+
+    // take into account other variables present in the rows
+    collectAffectedCols(col, iPlus, objRowCols, pPlus, pMinus, affectedCols);
+    collectAffectedCols(col, iMinus, objRowCols, pMinus, pPlus, affectedCols);
+
+    // compute correction term
+    int64_t correction = 0;
+    for (HighsInt k : affectedCols) {
+      correction += static_cast<int64_t>(pPlus[k]) * pMinus[k];
+      pPlus[k] = 0;
+      pMinus[k] = 0;
+    }
+
+    int64_t mPlus = static_cast<int64_t>(iPlus.size());
+    int64_t mMinus = static_cast<int64_t>(iMinus.size());
+    int64_t neOld = nePlus + neMinus;
+    // note that we subtract the entries for column 'col' since these are
+    // eliminated
+    int64_t neNew =
+        mPlus * (neMinus - mMinus) + mMinus * (nePlus - mPlus) - correction;
+    neRed = neOld - neNew;
+    mrRed = mPlus + mMinus - mPlus * mMinus;
+    return true;
+  };
+
+  auto checkNewRow = [&](const newRow& nr, bool& isRedundant) {
+    HighsCDouble impliedLower = 0;
+    HighsCDouble impliedUpper = 0;
+    bool lowerFinite = true;
+    bool upperFinite = true;
+    isRedundant = false;
+    for (const auto& e : nr.entries) {
+      double lb = model->col_lower_[e.col];
+      double ub = model->col_upper_[e.col];
+      if (e.val > 0) {
+        lowerFinite = lowerFinite && lb != -kHighsInf;
+        if (lowerFinite) impliedLower += e.val * lb;
+        upperFinite = upperFinite && ub != kHighsInf;
+        if (upperFinite) impliedUpper += e.val * ub;
+      } else {
+        lowerFinite = lowerFinite && ub != kHighsInf;
+        if (lowerFinite) impliedLower += e.val * ub;
+        upperFinite = upperFinite && lb != -kHighsInf;
+        if (upperFinite) impliedUpper += e.val * lb;
+      }
+      if (!lowerFinite && !upperFinite) return Result::kOk;
+    }
+
+    double lower = lowerFinite ? static_cast<double>(impliedLower) : -kHighsInf;
+    double upper = upperFinite ? static_cast<double>(impliedUpper) : kHighsInf;
+
+    // check for infeasibility
+    if (lower > nr.upper + primal_feastol || upper < nr.lower - primal_feastol)
+      return Result::kPrimalInfeasible;
+
+    // check for redundancy
+    isRedundant = lower >= nr.lower - primal_feastol &&
+                  upper <= nr.upper + primal_feastol;
+
+    return Result::kOk;
+  };
+
+  auto getRowData = [&](HighsInt row, HighsInt col, HighsInt multiplier,
+                        double& absCoef, HighsInt& direction, double& bound) {
+    if (row < 0) {
+      // artificial lower / upper bound row
+      direction = 1;
+      absCoef = 1.0;
+      bound = multiplier > 0 ? model->col_upper_[col] : -model->col_lower_[col];
+    } else {
+      HighsInt pPos = findNonzero(row, col);
+      assert(pPos != -1);
+      direction = multiplier * Avalue[pPos] > 0 ? HighsInt{1} : HighsInt{-1};
+      absCoef = std::abs(Avalue[pPos]);
+      bound = direction > 0 ? model->row_upper_[row] : -model->row_lower_[row];
+    }
+  };
+
+  auto collectRowEntries = [&](HighsInt row, HighsInt col, double scale,
+                               std::vector<newRowEntry>& newRowEntries,
+                               std::vector<HighsInt>& newRowMark) {
+    if (row < 0) return;
+    for (const auto& nz : getRowVector(row)) {
+      if (nz.index() == col) continue;
+      double val = scale * nz.value();
+      if (newRowMark[nz.index()] == -1) {
+        newRowMark[nz.index()] = static_cast<HighsInt>(newRowEntries.size());
+        newRowEntries.push_back({nz.index(), val});
+      } else {
+        newRowEntries[newRowMark[nz.index()]].val += val;
+      }
+    }
+  };
+
+  auto isReduction = [](int64_t neRed, int64_t mrRed) {
+    return neRed > 0 || (neRed == 0 && mrRed > 0);
+  };
+
+  auto insertOriginals =
+      [&](std::set<HighsInt>& rows,
+          const std::unordered_map<HighsInt, std::set<HighsInt>>& originals,
+          HighsInt row, HighsInt col) {
+        if (row == kUpperBoundRow)
+          rows.insert(-(2 * col + 1));
+        else if (row == kLowerBoundRow)
+          rows.insert(-(2 * col + 2));
+        else {
+          auto it = originals.find(row);
+          if (it != originals.end())
+            rows.insert(it->second.begin(), it->second.end());
+          else
+            rows.insert(row);
+        }
+      };
+
+  auto mergeOriginals =
+      [&](std::set<HighsInt>& rows,
+          const std::unordered_map<HighsInt, std::set<HighsInt>>& originals,
+          HighsInt plusRow, HighsInt minusRow, HighsInt col) {
+        rows.clear();
+        insertOriginals(rows, originals, plusRow, col);
+        insertOriginals(rows, originals, minusRow, col);
+      };
+
+  auto cernikovRedundant =
+      [&](std::set<HighsInt>& rows,
+          const std::unordered_map<HighsInt, std::set<HighsInt>>& originals,
+          HighsInt plusRow, HighsInt minusRow, HighsInt col,
+          HighsInt numColsElim) {
+        mergeOriginals(rows, originals, plusRow, minusRow, col);
+        return static_cast<HighsInt>(rows.size()) > numColsElim + 2;
+      };
+
+  // reformulate objective as a constraint: min c^T x + offset becomes
+  // min z with c^T x - z <= -offset. this allows FME to eliminate
+  // continuous columns with nonzero cost.
+  auto reformulateObjective = [&]() {
+    if (model->fme_obj_col_ != -1) {
+      assert(!colDeleted[model->fme_obj_col_]);
+      return;
+    }
+
+    HighsInt zCol = model->num_col_;
+    model->num_col_++;
+    model->a_matrix_.num_col_++;
+
+    // extend model vectors
+    model->col_cost_.push_back(1.0);
+    model->col_lower_.push_back(-kHighsInf);
+    model->col_upper_.push_back(kHighsInf);
+    model->integrality_.push_back(HighsVarType::kContinuous);
+    model->a_matrix_.start_.push_back(model->a_matrix_.start_.back());
+    if (model->col_names_.size() > 0) model->col_names_.push_back("fme_obj_z");
+
+    // extend presolve vectors
+    colhead.push_back(-1);
+    colsize.push_back(0);
+    colDeleted.push_back(0);
+    implColLower.push_back(-kHighsInf);
+    implColUpper.push_back(kHighsInf);
+    colLowerSource.push_back(-1);
+    colUpperSource.push_back(-1);
+    implRowDualSourceByCol.push_back({});
+    changedColFlag.push_back(1);
+    numProbes.push_back(0);
+
+    // update implied bound structures (pointers may be invalidated by
+    // reallocation of column vectors above)
+    impliedRowBounds.setBoundArrays(
+        model->col_lower_.data(), model->col_upper_.data(), implColLower.data(),
+        implColUpper.data(), colLowerSource.data(), colUpperSource.data());
+    impliedDualRowBounds.setNumSums(model->num_col_);
+
+    // register in postsolve stack
+    postsolve_stack.appendColToModel();
+
+    // build the objective constraint row: c^T x - z <= -offset
+    double offset = model->offset_;
+    std::vector<row_entry> objRow;
+    for (HighsInt j = 0; j < zCol; ++j) {
+      if (!colDeleted[j] && model->col_cost_[j] != 0.0)
+        objRow.push_back({j, model->col_cost_[j]});
+    }
+    objRow.push_back({zCol, -1.0});
+
+    // zero out original costs and offset
+    for (HighsInt j = 0; j < zCol; ++j) model->col_cost_[j] = 0.0;
+    model->offset_ = 0.0;
+
+    // add the constraint row to the matrix
+    addToMatrix(postsolve_stack, -kHighsInf, -offset, objRow);
+
+    // register reduction so getReducedPrimalSolution can compute z
+    std::vector<HighsPostsolveStack::Nonzero> costEntries;
+    for (const auto& entry : objRow)
+      costEntries.emplace_back(entry.col, entry.val);
+    postsolve_stack.fourierMotzkinObjCol(zCol, offset, costEntries);
+
+    model->fme_obj_col_ = zCol;
+
+    shrinkProblem(postsolve_stack);
+  };
+
+  auto collectCandidatesAndBuildHeap =
+      [&](std::vector<HighsInt>& candidates, Heap& heap,
+          std::vector<HighsInt>& iPlus, std::vector<HighsInt>& iMinus,
+          std::vector<HighsInt>& pPlus, std::vector<HighsInt>& pMinus,
+          std::vector<HighsInt>& affectedCols,
+          const std::vector<HighsInt>& objRowCols) {
+        // compute candidates
+        if (!computeCandidates(candidates)) return false;
+        // set up data structures for heap
+        heap.reset(model->num_col_, static_cast<HighsInt>(candidates.size()));
+        pPlus.assign(model->num_col_, 0);
+        pMinus.assign(model->num_col_, 0);
+        iPlus.reserve(model->num_row_);
+        iMinus.reserve(model->num_row_);
+        affectedCols.reserve(model->num_col_);
+        // inspect candidates (with limits)
+        HighsInt numFails = 0;
+        for (HighsInt col : candidates) {
+          int64_t neRed;
+          int64_t mrRed;
+          bool elimCandidate =
+              checkNonZeros(col, objRowCols, iPlus, iMinus, pPlus, pMinus,
+                            affectedCols, neRed, mrRed);
+          affectedCols.clear();
+          if (!elimCandidate || !isReduction(neRed, mrRed)) {
+            // count number of failures
+            if (++numFails > maxNumFails) break;
+            continue;
+          }
+          // add to heap
+          numFails = 0;
+          heap.push(col, neRed, mrRed);
+          if (heap.size() >= maxHeapSize) break;
+        }
+        if (heap.empty()) return false;
+        // heapify
+        heap.heapify();
+        return true;
+      };
+
+  // find index of a row within a list
+  auto findRowIndex = [](HighsInt row,
+                         const std::vector<FmeRow>& rows) -> HighsInt {
+    for (HighsInt i = 0; i < static_cast<HighsInt>(rows.size()); ++i)
+      if (rows[i].row == row) return i;
+    return -1;
+  };
+
+  auto collectRows = [&](const std::vector<HighsInt>& rows) {
+    std::vector<FmeRow> result;
+    for (HighsInt r : rows) {
+      if (r < 0) continue;
+      result.push_back(
+          {r, model->row_lower_[r], model->row_upper_[r], getRowVector(r)});
+    }
+    return result;
+  };
+
+  auto inheritAncestry =
+      [&](std::unordered_map<HighsInt, std::vector<FmeAncestryEntry>>&
+              rowAncestry,
+          HighsInt newModelRow, HighsInt parentRow, HighsInt parentRowIndex,
+          HighsInt stepIndex, double scale, bool isMinus) {
+        if (parentRow < 0) return;
+        auto it = rowAncestry.find(parentRow);
+        if (it != rowAncestry.end()) {
+          for (const auto& a : it->second)
+            rowAncestry[newModelRow].push_back(
+                {a.step, a.parentRowIndex, a.scale * scale, a.isMinus});
+        }
+        if (parentRowIndex >= 0)
+          rowAncestry[newModelRow].push_back(
+              {stepIndex, parentRowIndex, scale, isMinus});
+      };
+
+  auto printLog = [&](HighsInt colsRemoved, HighsInt rowsRemoved,
+                      HighsInt rowsAdded) {
+    highsLogDev(options->log_options, HighsLogType::kInfo,
+                "Fourier-Motzkin (%s objective reformulation) added "
+                "%" HIGHSINT_FORMAT " rows and eliminated %" HIGHSINT_FORMAT
+                " rows and %" HIGHSINT_FORMAT " columns\n",
+                options->presolve_fm_level >= 1 ? "with" : "without", rowsAdded,
+                rowsRemoved, colsRemoved);
+  };
+
+  // workspace vectors
+  std::vector<HighsInt> candidates;
+  std::vector<HighsInt> iPlus;
+  std::vector<HighsInt> iMinus;
+  std::vector<HighsInt> pPlus;
+  std::vector<HighsInt> pMinus;
+  std::vector<HighsInt> affectedCols;
+
+  // indexed max-heap
+  Heap heap;
+
+  // precompute the objective row: columns with nonzero cost
+  // used to simulate the objective constraint in checkRows before
+  // reformulation actually happens
+  std::vector<HighsInt> objRowCols;
+  if (model->fme_obj_col_ == -1 && options->presolve_fm_level >= 1) {
+    for (HighsInt j = 0; j < model->num_col_; ++j) {
+      if (!colDeleted[j] && model->col_cost_[j] != 0.0) objRowCols.push_back(j);
+    }
+  }
+
+  // compute candidates and build initial heap
+  if (!collectCandidatesAndBuildHeap(candidates, heap, iPlus, iMinus, pPlus,
+                                     pMinus, affectedCols, objRowCols))
+    return finalise();
+
+  // vectors for computing new row entries
+  std::vector<newRowEntry> newRowEntries;
+  std::vector<HighsInt> newRowMark(model->num_col_, -1);
+
+  // vector for storing new rows
+  std::vector<newRow> newRows;
+
+  // workspace for filtering new rows
+  std::vector<double> rowLower;
+  std::vector<double> rowUpper;
+  std::vector<std::vector<row_entry>> rowEntries;
+  std::vector<NewRowOrigin> newRowOrigins;
+
+  // vector for saving affected candidates
+  std::vector<HighsInt> saveAffectedCols;
+
+  // counters for numbers of eliminations
+  numColsEliminated = 0;
+  HighsInt numColsEliminatedBlock = 0;
+  HighsInt numRowsEliminated = 0;
+  HighsInt numRowsAdded = 0;
+
+  // FM block data for postsolve
+  std::vector<FmeBlockStep> blockSteps;
+
+  // surviving row to its ancestry (which parent rows it descends from)
+  std::unordered_map<HighsInt, std::vector<FmeAncestryEntry>> rowAncestry;
+
+  // distinct original parent rows for each derived row (Cernikov check)
+  std::unordered_map<HighsInt, std::set<HighsInt>> rowOriginals;
+  std::set<HighsInt> mergedOriginals;
+
+  // main loop: eliminate variables from heap
+  while (!heap.empty()) {
+    HighsInt col = heap.top();
+    heap.remove(col);
+
+    // if this candidate has nonzero cost and objective has not yet been
+    // reformulated, perform the reformulation now and rebuild the heap
+    if (model->fme_obj_col_ == -1 && model->col_cost_[col] != 0.0) {
+      // finalise any in-progress FM block before reformulating, since
+      // reformulateObjective pushes other reductions onto the data stack
+      if (!blockSteps.empty()) {
+        postsolve_stack.fourierMotzkinBlockFinalise(blockSteps, rowAncestry);
+        printLog(numColsEliminatedBlock, numRowsEliminated, numRowsAdded);
+        blockSteps.clear();
+        rowAncestry.clear();
+        rowOriginals.clear();
+        numColsEliminatedBlock = 0;
+        numRowsEliminated = 0;
+        numRowsAdded = 0;
+      }
+      // reformulate objective
+      reformulateObjective();
+      // clear vector for objective and resize marker
+      objRowCols.clear();
+      newRowMark.resize(model->num_col_, -1);
+      // re-compute candidates and re-build heap
+      if (!collectCandidatesAndBuildHeap(candidates, heap, iPlus, iMinus, pPlus,
+                                         pMinus, affectedCols, objRowCols))
+        return finalise();
+      continue;
+    }
+
+    // compute affected columns
+    int64_t neRed;
+    int64_t mrRed;
+    bool elimCandidate = checkNonZeros(col, objRowCols, iPlus, iMinus, pPlus,
+                                       pMinus, affectedCols, neRed, mrRed);
+
+    // heap data should be up-to-date
+    assert(elimCandidate && isReduction(neRed, mrRed));
+
+    HighsInt stepIdx = static_cast<HighsInt>(blockSteps.size());
+
+    // perform elimination: generate new rows
+    newRows.clear();
+    for (HighsInt pRow : iPlus) {
+      double pCoefAbs;
+      double pBound;
+      HighsInt pDirection;
+      getRowData(pRow, col, HighsInt{1}, pCoefAbs, pDirection, pBound);
+
+      for (HighsInt mRow : iMinus) {
+        double mCoefAbs;
+        double mBound;
+        HighsInt mDirection;
+        getRowData(mRow, col, HighsInt{-1}, mCoefAbs, mDirection, mBound);
+
+        // scale factor to preserve violation tolerances (see section 4.3):
+        double s = (pCoefAbs * mCoefAbs) / (pCoefAbs + mCoefAbs);
+        double pScale = s / pCoefAbs;
+        double mScale = s / mCoefAbs;
+
+        // collect row entries
+        collectRowEntries(pRow, col, pDirection * pScale, newRowEntries,
+                          newRowMark);
+        collectRowEntries(mRow, col, mDirection * mScale, newRowEntries,
+                          newRowMark);
+
+        // reset marker before removing near-zeros
+        for (const auto& e : newRowEntries) newRowMark[e.col] = -1;
+
+        // remove near-zero entries
+        newRowEntries.erase(
+            std::remove_if(newRowEntries.begin(), newRowEntries.end(),
+                           [&](const newRowEntry& e) {
+                             return abs(e.val) <= options->small_matrix_value;
+                           }),
+            newRowEntries.end());
+
+        // store new row
+        double new_upper =
+            static_cast<double>(static_cast<HighsCDouble>(pScale) * pBound +
+                                static_cast<HighsCDouble>(mScale) * mBound);
+        newRows.push_back({newRowEntries, -kHighsInf, new_upper, pRow, mRow,
+                           pDirection * pScale, mDirection * mScale});
+
+        // clear vector
+        newRowEntries.clear();
+      }
+    }
+
+    // add new rows, filtering out redundant ones
+    rowLower.clear();
+    rowUpper.clear();
+    rowEntries.clear();
+    newRowOrigins.clear();
+
+    for (const auto& nr : newRows) {
+      bool redundant = false;
+      HPRESOLVE_CHECKED_CALL(checkNewRow(nr, redundant));
+      if (redundant) continue;
+
+      // Cernikov redundancy check
+      if (cernikovRedundant(mergedOriginals, rowOriginals, nr.plusIndex,
+                            nr.minusIndex, col, numColsEliminated))
+        continue;
+
+      std::vector<row_entry> entries;
+      entries.reserve(nr.entries.size());
+      for (const auto& e : nr.entries)
+        entries.push_back({e.col, static_cast<double>(e.val)});
+      rowLower.push_back(nr.lower);
+      rowUpper.push_back(nr.upper);
+      rowEntries.push_back(std::move(entries));
+      newRowOrigins.push_back(
+          {nr.plusIndex, nr.minusIndex, nr.plusScale, nr.minusScale});
+    }
+
+    // serialize row data for postsolve before addToMatrix invalidates slices
+    std::vector<FmeRow> plusRows = collectRows(iPlus);
+    std::vector<FmeRow> minusRows = collectRows(iMinus);
+
+    // push row data for this elimination step onto the postsolve stack
+    postsolve_stack.fourierMotzkinBlockPushStep(col, plusRows, minusRows);
+
+    // save block metadata
+    assert(model->col_cost_[col] == 0.0);
+    blockSteps.push_back({col,
+                          model->col_lower_[col],
+                          model->col_upper_[col],
+                          static_cast<HighsInt>(plusRows.size()),
+                          static_cast<HighsInt>(minusRows.size()),
+                          {}});
+
+    // add new rows to matrix
+    HighsInt firstNewRow = model->num_row_;
+    if (!addToMatrix(postsolve_stack, rowLower, rowUpper, rowEntries))
+      return finalise();
+    numRowsAdded += static_cast<HighsInt>(rowEntries.size());
+
+    // build FmeNewRow data and ancestry for this step
+    auto& stepNewRows = blockSteps.back().newRows;
+    stepNewRows.reserve(newRowOrigins.size());
+    for (HighsInt k = 0; k < static_cast<HighsInt>(newRowOrigins.size()); ++k) {
+      HighsInt newModelRow = firstNewRow + k;
+      const auto& origin = newRowOrigins[k];
+      HighsInt pIdx = findRowIndex(origin.plusRow, plusRows);
+      HighsInt mIdx = findRowIndex(origin.minusRow, minusRows);
+      inheritAncestry(rowAncestry, newModelRow, origin.plusRow, pIdx, stepIdx,
+                      origin.plusScale, false);
+      inheritAncestry(rowAncestry, newModelRow, origin.minusRow, mIdx, stepIdx,
+                      origin.minusScale, true);
+      mergeOriginals(mergedOriginals, rowOriginals, origin.plusRow,
+                     origin.minusRow, col);
+      rowOriginals[newModelRow] = mergedOriginals;
+      stepNewRows.push_back({newModelRow, pIdx, mIdx});
+    }
+
+    // mark column as deleted
+    markColDeleted(col);
+    ++numColsEliminatedBlock;
+    ++numColsEliminated;
+
+    // remove old rows containing col (skip bound rows)
+    for (HighsInt rp : iPlus) {
+      if (rp < 0) continue;
+      rowAncestry.erase(rp);
+      rowOriginals.erase(rp);
+      removeRow(rp);
+      ++numRowsEliminated;
+    }
+    for (HighsInt rm : iMinus) {
+      if (rm < 0) continue;
+      rowAncestry.erase(rm);
+      rowOriginals.erase(rm);
+      if (rowDeleted[rm]) continue;
+      removeRow(rm);
+      ++numRowsEliminated;
+    }
+
+    // update affected candidates in the heap
+    saveAffectedCols.swap(affectedCols);
+    for (HighsInt k : saveAffectedCols) {
+      // check if variable is a candidate
+      bool isCandidateCol = isCandidate(k);
+      // skip variable if it is not on the heap and no candidate
+      if (!heap.contains(k) && !isCandidateCol) continue;
+      // check column non-zeros
+      int64_t ne, mr;
+      bool elimCandidate =
+          isCandidateCol && checkNonZeros(k, objRowCols, iPlus, iMinus, pPlus,
+                                          pMinus, affectedCols, ne, mr);
+      affectedCols.clear();
+      if (!elimCandidate || !isReduction(ne, mr)) {
+        // no candidate or not beneficial -> remove from heap
+        heap.remove(k);
+      } else if (!heap.contains(k)) {
+        // new candidate -> insert into heap
+        heap.insert(k, ne, mr);
+      } else {
+        // update heap
+        heap.update(k, ne, mr);
+      }
+    }
+    saveAffectedCols.clear();
+
+    if (checkLimits(postsolve_stack) != Result::kOk) break;
+  }
+
+  if (numColsEliminatedBlock > 0) {
+    // finalize the FM block
+    postsolve_stack.fourierMotzkinBlockFinalise(blockSteps, rowAncestry);
+
+    // log message
+    printLog(numColsEliminatedBlock, numRowsEliminated, numRowsAdded);
+  }
+
+  return finalise();
+}
+
 void HPresolve::substitute(HighsInt substcol, HighsInt staycol, double offset,
                            double scale) {
   // substitute the column in each row where it occurs
@@ -8988,7 +9960,7 @@ void HPresolve::debug(const HighsLp& lp, const HighsOptions& options) {
 
     sol = reducedsol;
     basis = reducedbasis;
-    postsolve_stack.undoUntil(options, sol, basis, tmp.numReductions());
+    postsolve_stack.undo(options, sol, basis, tmp.numReductions());
 
     HighsBasis temp_basis;
     HighsSolution temp_sol;
@@ -9034,7 +10006,7 @@ void HPresolve::debug(const HighsLp& lp, const HighsOptions& options) {
                          ARstart, ARindex, ARvalue);
     sol = reducedsol;
     basis = reducedbasis;
-    postsolve_stack.undoUntil(options, sol, basis, reductionLim);
+    postsolve_stack.undo(options, sol, basis, reductionLim);
 
     calculateRowValuesQuad(model, sol);
     kktinfo = dev_kkt_check::initInfo();
diff --git a/highs/presolve/HPresolve.h b/highs/presolve/HPresolve.h
index 3c30f894962..e653bceba71 100644
--- a/highs/presolve/HPresolve.h
+++ b/highs/presolve/HPresolve.h
@@ -163,6 +163,11 @@ class HPresolve {
     explicit operator Result() const { return my_result; };
   };
 
+  struct row_entry {
+    HighsInt col;
+    double val;
+  };
+
   HighsPresolveStatus presolve_status_;
   HPresolveAnalysis analysis_;
 
@@ -397,6 +402,14 @@ class HPresolve {
 
   void addToMatrix(const HighsInt row, const HighsInt col, const double val);
 
+  bool addToMatrix(HighsPostsolveStack& postsolve_stack,
+                   const std::vector<double>& row_lower,
+                   const std::vector<double>& row_upper,
+                   const std::vector<std::vector<row_entry>>& row_entries);
+
+  bool addToMatrix(HighsPostsolveStack& postsolve_stack, double row_lower,
+                   double row_upper, std::vector<row_entry> row_entries);
+
   Result prepareProbing(HighsPostsolveStack& postsolve_stack, bool& firstCall);
 
   Result finaliseProbing(HighsPostsolveStack& postsolve_stack, bool firstCall,
@@ -490,6 +503,9 @@ class HPresolve {
 
   Result aggregator(HighsPostsolveStack& postsolve_stack);
 
+  Result fourierMotzkin(HighsPostsolveStack& postsolve_stack,
+                        HighsInt& numColsEliminated);
+
   Result removeRowSingletons(HighsPostsolveStack& postsolve_stack);
 
   Result presolveColSingletons(HighsPostsolveStack& postsolve_stack);
@@ -531,6 +547,7 @@ class HPresolve {
 
   Result presolveRuleTest(HighsPostsolveStack& postsolve_stack);
   Result presolveRuleTestColStuffing(HighsPostsolveStack& postsolve_stack);
+  Result presolveRuleTestFourierMotzkin(HighsPostsolveStack& postsolve_stack);
 
   // Not currently called
   static void debug(const HighsLp& lp, const HighsOptions& options);
diff --git a/highs/presolve/HPresolveAnalysis.h b/highs/presolve/HPresolveAnalysis.h
index 707ccd06035..3a8b189b8b9 100644
--- a/highs/presolve/HPresolveAnalysis.h
+++ b/highs/presolve/HPresolveAnalysis.h
@@ -28,8 +28,6 @@ class HPresolveAnalysis {
   HighsInt original_num_row_;
 
  public:
-  std::vector<HighsBool> allow_rule_;
-
   bool allow_logging_;
   bool logging_on_;
 
diff --git a/highs/presolve/HPresolveTest.cpp b/highs/presolve/HPresolveTest.cpp
index 77a24d1e98d..6e0df138d4f 100644
--- a/highs/presolve/HPresolveTest.cpp
+++ b/highs/presolve/HPresolveTest.cpp
@@ -14,9 +14,12 @@ HPresolve::Result HPresolve::presolveRuleTest(
   assert(options->presolve_rule_test);
   if (options->presolve_rule_test == kPresolveRuleColStuffing) {
     return presolveRuleTestColStuffing(postsolve_stack);
+  } else if (options->presolve_rule_test == kPresolveRuleFourierMotzkin) {
+    return presolveRuleTestFourierMotzkin(postsolve_stack);
   }
   return Result::kOk;
 }
+
 HPresolve::Result HPresolve::presolveRuleTestColStuffing(
     HighsPostsolveStack& postsolve_stack) {
   assert(options->presolve_rule_test == kPresolveRuleColStuffing);
@@ -36,4 +39,23 @@ HPresolve::Result HPresolve::presolveRuleTestColStuffing(
   // Possibly remove the row
   return rowPresolve(postsolve_stack, 0);
 }
+
+HPresolve::Result HPresolve::presolveRuleTestFourierMotzkin(
+    HighsPostsolveStack& postsolve_stack) {
+  assert(options->presolve_rule_test == kPresolveRuleFourierMotzkin);
+  highsLogUser(options->log_options, HighsLogType::kInfo,
+               "HPresolve::presolveRuleTestFourierMotzkin\n");
+
+  HighsInt numColsEliminated;
+  HPresolve::Result result = fourierMotzkin(postsolve_stack, numColsEliminated);
+  if (result != Result::kOk) return result;
+
+  highsLogUser(options->log_options, HighsLogType::kInfo,
+               "HPresolve::presolveRuleTestFourierMotzkin: Removed %d "
+               "rows and %d columns\n",
+               int(numDeletedRows), int(numDeletedCols));
+  // Possibly remove the row
+  //  result = rowPresolve(postsolve_stack, 0);
+  return result;
+}
 }  // namespace presolve
diff --git a/highs/presolve/HighsPostsolveStack.cpp b/highs/presolve/HighsPostsolveStack.cpp
index 4f8168a2232..d5392d11409 100644
--- a/highs/presolve/HighsPostsolveStack.cpp
+++ b/highs/presolve/HighsPostsolveStack.cpp
@@ -21,10 +21,14 @@ void HighsPostsolveStack::initializeIndexMaps(HighsInt numRow,
                                               HighsInt numCol) {
   origNumRow = numRow;
   origNumCol = numCol;
+  nextRowIndex = numRow;
+  nextColIndex = numCol;
 
   origRowIndex.resize(numRow);
   std::iota(origRowIndex.begin(), origRowIndex.end(), 0);
 
+  origRowType.resize(numRow, OrigRowType::kOriginal);
+
   origColIndex.resize(numCol);
   std::iota(origColIndex.begin(), origColIndex.end(), 0);
 
@@ -41,6 +45,7 @@ void HighsPostsolveStack::compressIndexMaps(
 void HighsPostsolveStack::compressRowIndexMap(
     const std::vector<HighsInt>& newRowIndex) {
   compressIndexMap(newRowIndex, this->origRowIndex);
+  compressIndexMap(newRowIndex, this->origRowType);
 }
 
 void HighsPostsolveStack::compressColIndexMap(
@@ -48,22 +53,6 @@ void HighsPostsolveStack::compressColIndexMap(
   compressIndexMap(newColIndex, this->origColIndex);
 }
 
-void HighsPostsolveStack::compressIndexMap(
-    const std::vector<HighsInt>& newIndex, std::vector<HighsInt>& origIndex) {
-  // loop over entries, decrease entry counter for deleted entries
-  // (marked with -1), otherwise store original index at new index
-  // position
-  HighsInt numEn = origIndex.size();
-  for (size_t i = 0; i != newIndex.size(); ++i) {
-    if (newIndex[i] == -1)
-      --numEn;
-    else
-      origIndex[newIndex[i]] = origIndex[i];
-  }
-  // resize original index array to new size
-  origIndex.resize(numEn);
-}
-
 void HighsPostsolveStack::LinearTransform::undo(const HighsOptions& options,
                                                 HighsSolution& solution) const {
   solution.col_value[col] *= scale;
@@ -78,6 +67,24 @@ void HighsPostsolveStack::LinearTransform::transformToPresolvedSpace(
   primalSol[col] /= scale;
 }
 
+void HighsPostsolveStack::FourierMotzkinObjCol::transformToPresolvedSpace(
+    const std::vector<Nonzero>& costEntries,
+    std::vector<double>& primalSol) const {
+  double val = offset;
+  for (const Nonzero& entry : costEntries)
+    val += entry.value * primalSol[entry.index];
+  primalSol[col] = val;
+}
+
+void HighsPostsolveStack::FourierMotzkinObjCol::undo(
+    const std::vector<Nonzero>& costEntries, HighsSolution& solution) const {
+  if (!solution.dual_valid) return;
+  double zDual = solution.col_dual[col];
+  for (const Nonzero& entry : costEntries)
+    solution.col_dual[entry.index] += entry.value * zDual;
+  solution.col_dual[col] = 0.0;
+}
+
 static HighsBasisStatus computeRowStatus(double dual,
                                          HighsPostsolveStack::RowType rowType) {
   if (rowType == HighsPostsolveStack::RowType::kEq)
@@ -89,7 +96,8 @@ static HighsBasisStatus computeRowStatus(double dual,
 }
 
 void HighsPostsolveStack::FreeColSubstitution::undo(
-    const HighsOptions& options, const std::vector<Nonzero>& rowValues,
+    const HighsPostsolveStack& postsolveStack, const HighsOptions& options,
+    const std::vector<Nonzero>& rowValues,
     const std::vector<Nonzero>& colValues, HighsSolution& solution,
     HighsBasis& basis) {
   // compute primal values
@@ -104,7 +112,7 @@ void HighsPostsolveStack::FreeColSubstitution::undo(
 
   assert(colCoef != 0);
   // Row values aren't fully postsolved, so why do this?
-  if (solution.isModelRow(row))
+  if (postsolveStack.isModelRow(row))
     solution.row_value[row] =
         static_cast<double>(rowValue + colCoef * solution.col_value[col]);
   solution.col_value[col] = static_cast<double>((rhs - rowValue) / colCoef);
@@ -113,11 +121,11 @@ void HighsPostsolveStack::FreeColSubstitution::undo(
   if (!solution.dual_valid) return;
 
   // compute the row dual value such that reduced cost of basic column is 0
-  if (solution.isModelRow(row)) {
+  if (postsolveStack.isModelRow(row)) {
     solution.row_dual[row] = 0;
     HighsCDouble dualval = colCost;
     for (const auto& colVal : colValues) {
-      if (solution.isModelRow(colVal.index))
+      if (postsolveStack.isModelRow(colVal.index))
         dualval -= colVal.value * solution.row_dual[colVal.index];
     }
     solution.row_dual[row] = static_cast<double>(dualval / colCoef);
@@ -129,7 +137,7 @@ void HighsPostsolveStack::FreeColSubstitution::undo(
   if (!basis.valid) return;
 
   basis.col_status[col] = HighsBasisStatus::kBasic;
-  if (solution.isModelRow(row))
+  if (postsolveStack.isModelRow(row))
     basis.row_status[row] = computeRowStatus(solution.row_dual[row], rowType);
 }
 
@@ -154,8 +162,9 @@ static HighsBasisStatus computeStatus(double dual,
 }
 
 void HighsPostsolveStack::DoubletonEquation::undo(
-    const HighsOptions& options, const std::vector<Nonzero>& colValues,
-    HighsSolution& solution, HighsBasis& basis) const {
+    const HighsPostsolveStack& postsolveStack, const HighsOptions& options,
+    const std::vector<Nonzero>& colValues, HighsSolution& solution,
+    HighsBasis& basis) const {
   // retrieve the row and column index, the row side and the two
   // coefficients then compute the primal values
   solution.col_value[colSubst] = static_cast<double>(
@@ -179,10 +188,10 @@ void HighsPostsolveStack::DoubletonEquation::undo(
   // multiplier of this row i implicitly increases the dual multiplier of this
   // doubleton equation row with that scale.
   HighsCDouble rowDual = 0.0;
-  if (solution.isModelRow(row)) {
+  if (postsolveStack.isModelRow(row)) {
     solution.row_dual[row] = 0;
     for (const auto& colVal : colValues) {
-      if (solution.isModelRow(colVal.index))
+      if (postsolveStack.isModelRow(colVal.index))
         rowDual -= colVal.value * solution.row_dual[colVal.index];
     }
     rowDual /= coefSubst;
@@ -199,7 +208,7 @@ void HighsPostsolveStack::DoubletonEquation::undo(
     // so alter the dual multiplier of the row to make the dual multiplier of
     // column zero
     double rowDualDelta = solution.col_dual[col] / coef;
-    if (solution.isModelRow(row))
+    if (postsolveStack.isModelRow(row))
       solution.row_dual[row] = static_cast<double>(rowDual + rowDualDelta);
     solution.col_dual[col] = 0.0;
     solution.col_dual[colSubst] = static_cast<double>(
@@ -220,7 +229,7 @@ void HighsPostsolveStack::DoubletonEquation::undo(
     // otherwise make the reduced cost of the substituted column zero and make
     // that column basic
     double rowDualDelta = solution.col_dual[colSubst] / coefSubst;
-    if (solution.isModelRow(row))
+    if (postsolveStack.isModelRow(row))
       solution.row_dual[row] = static_cast<double>(rowDual + rowDualDelta);
     solution.col_dual[colSubst] = 0.0;
     solution.col_dual[col] =
@@ -231,15 +240,17 @@ void HighsPostsolveStack::DoubletonEquation::undo(
 
   if (!basis.valid) return;
 
-  if (solution.isModelRow(row))
+  if (postsolveStack.isModelRow(row))
     basis.row_status[row] = computeRowStatus(solution.row_dual[row], rowType);
 }
 
 void HighsPostsolveStack::EqualityRowAddition::undo(
-    const HighsOptions& options, const std::vector<Nonzero>& eqRowValues,
-    HighsSolution& solution, HighsBasis& basis) const {
+    const HighsPostsolveStack& postsolveStack, const HighsOptions& options,
+    const std::vector<Nonzero>& eqRowValues, HighsSolution& solution,
+    HighsBasis& basis) const {
   // (removed) cuts may have been used in this reduction.
-  if (!solution.isModelRow(row) || !solution.isModelRow(addedEqRow)) return;
+  if (!postsolveStack.isModelRow(row) || !postsolveStack.isModelRow(addedEqRow))
+    return;
 
   // nothing more to do if the row is zero in the dual solution or there is
   // no dual solution
@@ -255,11 +266,12 @@ void HighsPostsolveStack::EqualityRowAddition::undo(
 }
 
 void HighsPostsolveStack::EqualityRowAdditions::undo(
-    const HighsOptions& options, const std::vector<Nonzero>& eqRowValues,
+    const HighsPostsolveStack& postsolveStack, const HighsOptions& options,
+    const std::vector<Nonzero>& eqRowValues,
     const std::vector<Nonzero>& targetRows, HighsSolution& solution,
     HighsBasis& basis) const {
   // a (removed) cut may have been used in this reduction.
-  if (!solution.isModelRow(addedEqRow)) return;
+  if (!postsolveStack.isModelRow(addedEqRow)) return;
 
   // nothing more to do if the row is zero in the dual solution or there is
   // no dual solution
@@ -270,7 +282,7 @@ void HighsPostsolveStack::EqualityRowAdditions::undo(
   // used for adding the equation
   HighsCDouble eqRowDual = solution.row_dual[addedEqRow];
   for (const auto& targetRow : targetRows) {
-    if (solution.isModelRow(targetRow.index))
+    if (postsolveStack.isModelRow(targetRow.index))
       eqRowDual += static_cast<HighsCDouble>(targetRow.value) *
                    solution.row_dual[targetRow.index];
   }
@@ -280,8 +292,9 @@ void HighsPostsolveStack::EqualityRowAdditions::undo(
 }
 
 void HighsPostsolveStack::ForcingColumn::undo(
-    const HighsOptions& options, const std::vector<Nonzero>& colValues,
-    HighsSolution& solution, HighsBasis& basis) const {
+    const HighsPostsolveStack& postsolveStack, const HighsOptions& options,
+    const std::vector<Nonzero>& colValues, HighsSolution& solution,
+    HighsBasis& basis) const {
   HighsInt nonbasicRow = -1;
   HighsBasisStatus nonbasicRowStatus = HighsBasisStatus::kNonbasic;
   double colValFromNonbasicRow = colBound;
@@ -294,7 +307,7 @@ void HighsPostsolveStack::ForcingColumn::undo(
     for (const auto& colVal : colValues) {
       // Row values aren't fully postsolved, so how can this work?
       debug_num_use_row_value++;
-      if (solution.isModelRow(colVal.index)) {
+      if (postsolveStack.isModelRow(colVal.index)) {
         double colValFromRow = solution.row_value[colVal.index] / colVal.value;
         if (direction * colValFromRow > direction * colValFromNonbasicRow) {
           nonbasicRow = colVal.index;
@@ -344,10 +357,11 @@ void HighsPostsolveStack::ForcingColumn::undo(
 }
 
 void HighsPostsolveStack::ForcingColumnRemovedRow::undo(
-    const HighsOptions& options, const std::vector<Nonzero>& rowValues,
-    HighsSolution& solution, HighsBasis& basis) const {
+    const HighsPostsolveStack& postsolveStack, const HighsOptions& options,
+    const std::vector<Nonzero>& rowValues, HighsSolution& solution,
+    HighsBasis& basis) const {
   // a (removed) cut may have been used in this reduction.
-  if (!solution.isModelRow(row)) return;
+  if (!postsolveStack.isModelRow(row)) return;
 
   // we use the row value as storage for the scaled value implied on the
   // column dual
@@ -362,9 +376,9 @@ void HighsPostsolveStack::ForcingColumnRemovedRow::undo(
   if (basis.valid) basis.row_status[row] = HighsBasisStatus::kBasic;
 }
 
-void HighsPostsolveStack::SingletonRow::undo(const HighsOptions& options,
-                                             HighsSolution& solution,
-                                             HighsBasis& basis) const {
+void HighsPostsolveStack::SingletonRow::undo(
+    const HighsPostsolveStack& postsolveStack, const HighsOptions& options,
+    HighsSolution& solution, HighsBasis& basis) const {
   // nothing to do if the row's dual value is zero in the dual
   // solution or there is no dual solution
   if (!solution.dual_valid) return;
@@ -380,7 +394,7 @@ void HighsPostsolveStack::SingletonRow::undo(const HighsOptions& options,
       (!colUpperTightened || colStatus != HighsBasisStatus::kUpper)) {
     // the tightened bound is not used in the basic solution
     // hence we simply make the row basic and give it a dual multiplier of 0
-    if (solution.isModelRow(row)) {
+    if (postsolveStack.isModelRow(row)) {
       if (basis.valid) basis.row_status[row] = HighsBasisStatus::kBasic;
       solution.row_dual[row] = 0;
     }
@@ -389,13 +403,13 @@ void HighsPostsolveStack::SingletonRow::undo(const HighsOptions& options,
 
   // choose the row dual value such that the columns reduced cost becomes
   // zero
-  if (solution.isModelRow(row))
+  if (postsolveStack.isModelRow(row))
     solution.row_dual[row] = solution.col_dual[col] / coef;
   solution.col_dual[col] = 0;
 
   if (!basis.valid) return;
 
-  if (solution.isModelRow(row)) {
+  if (postsolveStack.isModelRow(row)) {
     switch (colStatus) {
       case HighsBasisStatus::kLower:
         assert(colLowerTightened);
@@ -425,10 +439,10 @@ void HighsPostsolveStack::SingletonRow::undo(const HighsOptions& options,
 }
 
 // column fixed to lower or upper bound
-void HighsPostsolveStack::FixedCol::undo(const HighsOptions& options,
-                                         const std::vector<Nonzero>& colValues,
-                                         HighsSolution& solution,
-                                         HighsBasis& basis) const {
+void HighsPostsolveStack::FixedCol::undo(
+    const HighsPostsolveStack& postsolveStack, const HighsOptions& options,
+    const std::vector<Nonzero>& colValues, HighsSolution& solution,
+    HighsBasis& basis) const {
   // set solution value
   solution.col_value[col] = fixValue;
 
@@ -438,7 +452,7 @@ void HighsPostsolveStack::FixedCol::undo(const HighsOptions& options,
 
   HighsCDouble reducedCost = colCost;
   for (const auto& colVal : colValues) {
-    if (solution.isModelRow(colVal.index))
+    if (postsolveStack.isModelRow(colVal.index))
       reducedCost -= colVal.value * solution.row_dual[colVal.index];
   }
 
@@ -454,11 +468,11 @@ void HighsPostsolveStack::FixedCol::undo(const HighsOptions& options,
   }
 }
 
-void HighsPostsolveStack::RedundantRow::undo(const HighsOptions& options,
-                                             HighsSolution& solution,
-                                             HighsBasis& basis) const {
+void HighsPostsolveStack::RedundantRow::undo(
+    const HighsPostsolveStack& postsolveStack, const HighsOptions& options,
+    HighsSolution& solution, HighsBasis& basis) const {
   // a (removed) cut may have been used in this reduction.
-  if (!solution.isModelRow(row)) return;
+  if (!postsolveStack.isModelRow(row)) return;
 
   // set row dual to zero if dual solution requested
   if (!solution.dual_valid) return;
@@ -468,9 +482,23 @@ void HighsPostsolveStack::RedundantRow::undo(const HighsOptions& options,
   if (basis.valid) basis.row_status[row] = HighsBasisStatus::kBasic;
 }
 
+void HighsPostsolveStack::ImpliedEquation::undo(
+    const HighsPostsolveStack& postsolveStack,
+    const std::vector<Nonzero>& rowValues, HighsSolution& solution) const {
+  if (!solution.dual_valid) return;
+  if (!postsolveStack.isModelRow(row)) return;
+  double oldDual = solution.row_dual[row];
+  if (atLower ? (oldDual < 0) : (oldDual > 0)) {
+    solution.row_dual[row] = 0;
+    for (const auto& nz : rowValues)
+      solution.col_dual[nz.index] += nz.value * oldDual;
+  }
+}
+
 void HighsPostsolveStack::ForcingRow::undo(
-    const HighsOptions& options, const std::vector<Nonzero>& rowValues,
-    HighsSolution& solution, HighsBasis& basis) const {
+    const HighsPostsolveStack& postsolveStack, const HighsOptions& options,
+    const std::vector<Nonzero>& rowValues, HighsSolution& solution,
+    HighsBasis& basis) const {
   if (!solution.dual_valid) return;
 
   // compute the row dual multiplier and determine the new basic column
@@ -489,7 +517,7 @@ void HighsPostsolveStack::ForcingRow::undo(
   }
 
   if (basicCol != -1) {
-    if (solution.isModelRow(row))
+    if (postsolveStack.isModelRow(row))
       solution.row_dual[row] = solution.row_dual[row] + dualDelta;
     for (const auto& rowVal : rowValues) {
       solution.col_dual[rowVal.index] = static_cast<double>(
@@ -499,7 +527,7 @@ void HighsPostsolveStack::ForcingRow::undo(
     solution.col_dual[basicCol] = 0;
 
     if (basis.valid) {
-      if (solution.isModelRow(row))
+      if (postsolveStack.isModelRow(row))
         basis.row_status[row] =
             (rowType == RowType::kGeq ? HighsBasisStatus::kLower
                                       : HighsBasisStatus::kUpper);
@@ -509,17 +537,17 @@ void HighsPostsolveStack::ForcingRow::undo(
   }
 }
 
-void HighsPostsolveStack::DuplicateRow::undo(const HighsOptions& options,
-                                             HighsSolution& solution,
-                                             HighsBasis& basis) const {
+void HighsPostsolveStack::DuplicateRow::undo(
+    const HighsPostsolveStack& postsolveStack, const HighsOptions& options,
+    HighsSolution& solution, HighsBasis& basis) const {
   // (removed) cuts may have been used in this reduction.
-  if (!solution.isModelRow(row)) return;
+  if (!postsolveStack.isModelRow(row)) return;
 
   if (!solution.dual_valid) return;
   if (!rowUpperTightened && !rowLowerTightened) {
     // simple case of row2 being redundant, in which case it just gets a
     // dual multiplier of 0 and is made basic
-    if (solution.isModelRow(duplicateRow)) {
+    if (postsolveStack.isModelRow(duplicateRow)) {
       solution.row_dual[duplicateRow] = 0.0;
       if (basis.valid)
         basis.row_status[duplicateRow] = HighsBasisStatus::kBasic;
@@ -534,9 +562,9 @@ void HighsPostsolveStack::DuplicateRow::undo(const HighsOptions& options,
           : computeStatus(solution.row_dual[row], basis.row_status[row],
                           options.dual_feasibility_tolerance);
 
-  auto computeRowDualAndStatus = [&](bool tighened) {
-    if (tighened) {
-      if (solution.isModelRow(duplicateRow)) {
+  auto computeRowDualAndStatus = [&](bool tightened) {
+    if (tightened) {
+      if (postsolveStack.isModelRow(duplicateRow)) {
         solution.row_dual[duplicateRow] =
             solution.row_dual[row] / duplicateRowScale;
         if (basis.valid) {
@@ -548,7 +576,7 @@ void HighsPostsolveStack::DuplicateRow::undo(const HighsOptions& options,
       }
       solution.row_dual[row] = 0.0;
       if (basis.valid) basis.row_status[row] = HighsBasisStatus::kBasic;
-    } else if (solution.isModelRow(duplicateRow)) {
+    } else if (postsolveStack.isModelRow(duplicateRow)) {
       solution.row_dual[duplicateRow] = 0.0;
       if (basis.valid)
         basis.row_status[duplicateRow] = HighsBasisStatus::kBasic;
@@ -562,7 +590,7 @@ void HighsPostsolveStack::DuplicateRow::undo(const HighsOptions& options,
   switch (rowStatus) {
     case HighsBasisStatus::kBasic:
       // if row is basic the parallel row is also basic
-      if (solution.isModelRow(duplicateRow)) {
+      if (postsolveStack.isModelRow(duplicateRow)) {
         solution.row_dual[duplicateRow] = 0.0;
         if (basis.valid)
           basis.row_status[duplicateRow] = HighsBasisStatus::kBasic;
@@ -1332,8 +1360,9 @@ void HighsPostsolveStack::DuplicateColumn::transformToPresolvedSpace(
 }
 
 void HighsPostsolveStack::SlackColSubstitution::undo(
-    const HighsOptions& options, const std::vector<Nonzero>& rowValues,
-    HighsSolution& solution, HighsBasis& basis) {
+    const HighsPostsolveStack& postsolveStack, const HighsOptions& options,
+    const std::vector<Nonzero>& rowValues, HighsSolution& solution,
+    HighsBasis& basis) {
   bool debug_print = false;
   // May have to determine row dual and basis status unless doing
   // primal-only transformation in MIP solver, in which case row may
@@ -1352,7 +1381,7 @@ void HighsPostsolveStack::SlackColSubstitution::undo(
 
   assert(colCoef != 0);
   // Row values aren't fully postsolved, so why do this?
-  if (solution.isModelRow(row))
+  if (postsolveStack.isModelRow(row))
     solution.row_value[row] =
         static_cast<double>(rowValue + colCoef * solution.col_value[col]);
 
@@ -1362,14 +1391,14 @@ void HighsPostsolveStack::SlackColSubstitution::undo(
   if (!solution.dual_valid) return;
 
   // Row retains its dual value, and column has this dual value scaled by coeff
-  if (solution.isModelRow(row))
+  if (postsolveStack.isModelRow(row))
     solution.col_dual[col] = -solution.row_dual[row] / colCoef;
 
   // Set basis status if necessary
   if (!basis.valid) return;
 
   // If row is basic, then slack is basic, otherwise row retains its status
-  if (solution.isModelRow(row)) {
+  if (postsolveStack.isModelRow(row)) {
     HighsBasisStatus save_row_basis_status = basis.row_status[row];
     if (basis.row_status[row] == HighsBasisStatus::kBasic) {
       basis.col_status[col] = HighsBasisStatus::kBasic;
@@ -1399,4 +1428,366 @@ void HighsPostsolveStack::SlackColSubstitution::undo(
   }
 }
 
+std::vector<HighsPostsolveStack::FmeStepData>
+HighsPostsolveStack::popFourierMotzkinBlock(HighsDataStack& stack) {
+  HighsInt numSteps;
+  stack.pop(numSteps);
+
+  std::vector<FmeStepData> steps(numSteps);
+
+  // step headers
+  for (HighsInt s = numSteps - 1; s >= 0; --s) stack.pop(steps[s].header);
+
+  // new row origins
+  for (HighsInt s = numSteps - 1; s >= 0; --s) stack.pop(steps[s].newRows);
+
+  // descendants
+  for (HighsInt s = numSteps - 1; s >= 0; --s) {
+    HighsInt numMinus = steps[s].header.numMinus;
+    steps[s].minusDescendants.resize(numMinus);
+    for (HighsInt m = numMinus - 1; m >= 0; --m)
+      stack.pop(steps[s].minusDescendants[m]);
+    HighsInt numPlus = steps[s].header.numPlus;
+    steps[s].plusDescendants.resize(numPlus);
+    for (HighsInt p = numPlus - 1; p >= 0; --p)
+      stack.pop(steps[s].plusDescendants[p]);
+  }
+
+  // row data
+  for (HighsInt s = numSteps - 1; s >= 0; --s) {
+    // minus row data
+    stack.pop(steps[s].minusHeaders);
+    stack.pop(steps[s].minusCoefs);
+    HighsInt numMinus = static_cast<HighsInt>(steps[s].minusCoefs.size());
+    steps[s].minusEntries.resize(numMinus);
+    for (HighsInt r = numMinus - 1; r >= 0; --r)
+      stack.pop(steps[s].minusEntries[r]);
+
+    // plus row data
+    stack.pop(steps[s].plusHeaders);
+    stack.pop(steps[s].plusCoefs);
+    HighsInt numPlus = static_cast<HighsInt>(steps[s].plusCoefs.size());
+    steps[s].plusEntries.resize(numPlus);
+    for (HighsInt r = numPlus - 1; r >= 0; --r)
+      stack.pop(steps[s].plusEntries[r]);
+  }
+
+  return steps;
+}
+
+void HighsPostsolveStack::undoFourierMotzkinBlock(
+    const std::vector<FmeStepData>& steps, const HighsOptions& options,
+    HighsSolution& solution, HighsBasis& basis) {
+  const double tol = options.primal_feasibility_tolerance;
+  const double dual_tol = options.dual_feasibility_tolerance;
+
+  HighsInt numSteps = static_cast<HighsInt>(steps.size());
+
+  // primal postsolve (Algorithm 3): process in reverse elimination order
+  for (HighsInt s = numSteps - 1; s >= 0; --s) {
+    const auto& step = steps[s];
+    HighsInt col = step.header.col;
+    double lower = step.header.colLower;
+    double upper = step.header.colUpper;
+
+    auto tightenBounds = [&](const std::vector<FmeRowHeader>& headers,
+                             const std::vector<double>& coefs,
+                             const std::vector<std::vector<Nonzero>>& entries,
+                             double& lowerBound, double& upperBound) {
+      for (size_t r = 0; r < headers.size(); ++r) {
+        double aij = coefs[r];
+        HighsCDouble sum = 0.0;
+        for (const auto& nz : entries[r])
+          sum += static_cast<HighsCDouble>(nz.value) *
+                 solution.col_value[nz.index];
+        HighsInt direction = aij > 0 ? HighsInt{1} : HighsInt{-1};
+        double rhs_upper =
+            direction > 0 ? headers[r].rowUpper : headers[r].rowLower;
+        double rhs_lower =
+            direction > 0 ? headers[r].rowLower : headers[r].rowUpper;
+        if (direction * rhs_upper != kHighsInf) {
+          double bound = static_cast<double>((rhs_upper - sum) / aij);
+          upperBound = std::min(upperBound, bound);
+        }
+        if (direction * rhs_lower != -kHighsInf) {
+          double bound = static_cast<double>((rhs_lower - sum) / aij);
+          lowerBound = std::max(lowerBound, bound);
+        }
+      }
+    };
+
+    tightenBounds(step.plusHeaders, step.plusCoefs, step.plusEntries, lower,
+                  upper);
+    tightenBounds(step.minusHeaders, step.minusCoefs, step.minusEntries, lower,
+                  upper);
+
+    if (lower <= tol && upper >= -tol)
+      solution.col_value[col] = 0.0;
+    else if (lower > 0.0)
+      solution.col_value[col] = lower;
+    else
+      solution.col_value[col] = upper;
+  }
+
+  if (!solution.dual_valid) return;
+
+  // dual postsolve (Algorithm 4): process in reverse elimination order
+  for (HighsInt s = numSteps - 1; s >= 0; --s) {
+    const auto& step = steps[s];
+    HighsInt col = step.header.col;
+    HighsInt numPlus = step.header.numPlus;
+    HighsInt numMinus = step.header.numMinus;
+
+    // u_i = Σ_{k ∈ K^j_i} λ_k * scaleFactor
+    auto recoverDual =
+        [&](const std::vector<FmeRowHeader>& headers,
+            const std::vector<std::vector<FmeDescendant>>& descendants) {
+          for (size_t r = 0; r < headers.size(); ++r) {
+            HighsCDouble dual = 0.0;
+            for (const auto& desc : descendants[r])
+              dual += static_cast<HighsCDouble>(solution.row_dual[desc.row]) *
+                      desc.scaleFactor;
+            solution.row_dual[headers[r].row] += static_cast<double>(dual);
+          }
+        };
+    recoverDual(step.plusHeaders, step.plusDescendants);
+    recoverDual(step.minusHeaders, step.minusDescendants);
+
+    // col_dual = -Σ a_{ij} * row_dual[i] (cost is zero after reformulation)
+    HighsCDouble colDual = 0.0;
+    std::vector<bool> visited(solution.row_dual.size(), false);
+    for (HighsInt r = 0; r < numPlus; ++r) {
+      HighsInt row = step.plusHeaders[r].row;
+      colDual -=
+          static_cast<HighsCDouble>(step.plusCoefs[r]) * solution.row_dual[row];
+      visited[row] = true;
+    }
+    for (HighsInt r = 0; r < numMinus; ++r) {
+      HighsInt row = step.minusHeaders[r].row;
+      if (visited[row]) continue;
+      colDual -= static_cast<HighsCDouble>(step.minusCoefs[r]) *
+                 solution.row_dual[row];
+    }
+    solution.col_dual[col] = static_cast<double>(colDual);
+  }
+
+  // basis postsolve: use dual solution to determine basis status
+  if (!basis.valid) return;
+
+  // pre-compute lower and upper slacks for each row
+  auto computeSlacks =
+      [&](HighsInt col, const std::vector<FmeRowHeader>& headers,
+          const std::vector<double>& coefs,
+          const std::vector<std::vector<Nonzero>>& entries,
+          std::vector<double>& lowerSlacks, std::vector<double>& upperSlacks) {
+        HighsInt n = static_cast<HighsInt>(headers.size());
+        lowerSlacks.resize(n);
+        upperSlacks.resize(n);
+        for (HighsInt r = 0; r < n; ++r) {
+          HighsCDouble activity =
+              static_cast<HighsCDouble>(coefs[r]) * solution.col_value[col];
+          for (const auto& nz : entries[r])
+            activity += static_cast<HighsCDouble>(nz.value) *
+                        solution.col_value[nz.index];
+          double act = static_cast<double>(activity);
+          lowerSlacks[r] = headers[r].rowLower != -kHighsInf
+                               ? act - headers[r].rowLower
+                               : kHighsInf;
+          upperSlacks[r] = headers[r].rowUpper != kHighsInf
+                               ? headers[r].rowUpper - act
+                               : kHighsInf;
+        }
+      };
+
+  // row must be basic if it has zero dual and activity strictly
+  // between bounds (complementary slackness)
+  auto rowMustBeBasic = [&](HighsInt row, double lowerSlack,
+                            double upperSlack) {
+    return std::abs(solution.row_dual[row]) <= dual_tol && lowerSlack > tol &&
+           upperSlack > tol;
+  };
+
+  // assign row as basic
+  auto assignBasicRowStatus = [&](HighsInt row, HighsInt& basicAssigned) {
+    if (basis.row_status[row] == HighsBasisStatus::kBasic ||
+        std::abs(solution.row_dual[row]) > dual_tol)
+      return false;
+    basis.row_status[row] = HighsBasisStatus::kBasic;
+    basicAssigned++;
+    return true;
+  };
+
+  // assign row as non-basic
+  auto assignNonBasicRowStatus = [&](HighsInt row, double lowerSlack,
+                                     double upperSlack) {
+    if (solution.row_dual[row] > dual_tol)
+      basis.row_status[row] = HighsBasisStatus::kLower;
+    else if (solution.row_dual[row] < -dual_tol)
+      basis.row_status[row] = HighsBasisStatus::kUpper;
+    else
+      basis.row_status[row] = upperSlack < lowerSlack
+                                  ? HighsBasisStatus::kUpper
+                                  : HighsBasisStatus::kLower;
+  };
+
+  // assign row status
+  auto assignRowStatus = [&](HighsInt row, double lowerSlack, double upperSlack,
+                             HighsInt& basicAssigned,
+                             bool forceNonBasic = false) {
+    if (forceNonBasic || !assignBasicRowStatus(row, basicAssigned))
+      assignNonBasicRowStatus(row, lowerSlack, upperSlack);
+  };
+
+  // collect a single candidate for basic assignment
+  auto collectCandidate =
+      [&](HighsInt row, bool forcedNonBasic, double lowerSlack,
+          double upperSlack, const std::vector<Nonzero>& entries,
+          HighsInt parentIndex, bool isMinus,
+          std::vector<std::tuple<HighsInt, HighsInt, bool>>& candidates) {
+        if (basis.row_status[row] == HighsBasisStatus::kBasic) return;
+        if (forcedNonBasic || std::abs(solution.row_dual[row]) > dual_tol) {
+          assignNonBasicRowStatus(row, lowerSlack, upperSlack);
+          return;
+        }
+        HighsInt nonBasicCount = 0;
+        for (const auto& nz : entries)
+          if (basis.col_status[nz.index] != HighsBasisStatus::kBasic)
+            nonBasicCount++;
+        candidates.emplace_back(nonBasicCount, parentIndex, isMinus);
+      };
+
+  for (HighsInt s = numSteps - 1; s >= 0; --s) {
+    const auto& step = steps[s];
+    HighsInt col = step.header.col;
+    HighsInt numPlus = step.header.numPlus;
+    HighsInt numMinus = step.header.numMinus;
+
+    // compute slacks
+    std::vector<double> plusLowerSlack;
+    std::vector<double> plusUpperSlack;
+    std::vector<double> minusLowerSlack;
+    std::vector<double> minusUpperSlack;
+    computeSlacks(col, step.plusHeaders, step.plusCoefs, step.plusEntries,
+                  plusLowerSlack, plusUpperSlack);
+    computeSlacks(col, step.minusHeaders, step.minusCoefs, step.minusEntries,
+                  minusLowerSlack, minusUpperSlack);
+
+    // non-basic propagation: if a generated row is non-basic (with nonzero
+    // dual), both its parents are forced non-basic. mark them so the greedy
+    // passes skip them. only force if the parent doesn't must-be-basic.
+    std::vector<bool> forcedNonBasicPlus(numPlus, false);
+    std::vector<bool> forcedNonBasicMinus(numMinus, false);
+    for (const auto& nr : step.newRows) {
+      // get indices of parent rows
+      HighsInt p = nr.plusParentIdx;
+      HighsInt m = nr.minusParentIdx;
+      // skip basic rows (zero dual) and degenerate non-basic rows (with zero
+      // dual)
+      if (p < 0 || m < 0 || std::abs(solution.row_dual[nr.row]) <= dual_tol)
+        continue;
+      // mark rows that do not have to be basic
+      if (!rowMustBeBasic(step.plusHeaders[p].row, plusLowerSlack[p],
+                          plusUpperSlack[p]))
+        forcedNonBasicPlus[p] = true;
+      if (!rowMustBeBasic(step.minusHeaders[m].row, minusLowerSlack[m],
+                          minusUpperSlack[m]))
+        forcedNonBasicMinus[m] = true;
+    }
+
+    // mark ranged rows (appearing in both plus and minus sets)
+    std::vector<bool> isMinusRowRanged(numMinus, false);
+    HighsInt numRanged = 0;
+    for (HighsInt m = 0; m < numMinus; ++m)
+      for (HighsInt p = 0; p < numPlus; ++p)
+        if (step.minusHeaders[m].row == step.plusHeaders[p].row) {
+          isMinusRowRanged[m] = true;
+          numRanged++;
+          break;
+        }
+
+    // count number of basic new rows
+    HighsInt numNewRows = static_cast<HighsInt>(step.newRows.size());
+    HighsInt numBasicNewRows = 0;
+    for (const auto& nr : step.newRows)
+      if (basis.row_status[nr.row] == HighsBasisStatus::kBasic)
+        numBasicNewRows++;
+
+    // how many basic variables are needed?
+    HighsInt basicNeeded =
+        (numPlus + numMinus - numRanged - numNewRows) + numBasicNewRows;
+    HighsInt basicAssigned = 0;
+
+    // determine col status
+    bool colMustBeBasic =
+        solution.col_value[col] > step.header.colLower + tol &&
+        solution.col_value[col] < step.header.colUpper - tol;
+    bool colCanBeBasic =
+        colMustBeBasic || std::abs(solution.col_dual[col]) <= dual_tol;
+
+    // pass 1: assign all must-be-basic (col and rows)
+    if (colMustBeBasic) {
+      basis.col_status[col] = HighsBasisStatus::kBasic;
+      basicAssigned++;
+    }
+    for (HighsInt p = 0; p < numPlus; ++p) {
+      if (forcedNonBasicPlus[p]) continue;
+      if (rowMustBeBasic(step.plusHeaders[p].row, plusLowerSlack[p],
+                         plusUpperSlack[p]))
+        assignBasicRowStatus(step.plusHeaders[p].row, basicAssigned);
+    }
+    for (HighsInt m = 0; m < numMinus; ++m) {
+      if (isMinusRowRanged[m] || forcedNonBasicMinus[m]) continue;
+      if (rowMustBeBasic(step.minusHeaders[m].row, minusLowerSlack[m],
+                         minusUpperSlack[m]))
+        assignBasicRowStatus(step.minusHeaders[m].row, basicAssigned);
+    }
+
+    // pass 2: assign can-be-basic col (if not already assigned)
+    if (!colMustBeBasic) {
+      if (colCanBeBasic && basicAssigned < basicNeeded) {
+        basis.col_status[col] = HighsBasisStatus::kBasic;
+        basicAssigned++;
+      } else if (solution.col_value[col] <= step.header.colLower + tol) {
+        basis.col_status[col] = HighsBasisStatus::kLower;
+      } else {
+        basis.col_status[col] = HighsBasisStatus::kUpper;
+      }
+    }
+
+    // pass 3: assign can-be-basic rows, sorted by non-basic support count
+    // to reduce risk of rank deficiency in degenerate cases
+    std::vector<std::tuple<HighsInt, HighsInt, bool>> candidates;
+    for (HighsInt p = 0; p < numPlus; ++p)
+      collectCandidate(step.plusHeaders[p].row, forcedNonBasicPlus[p],
+                       plusLowerSlack[p], plusUpperSlack[p],
+                       step.plusEntries[p], p, false, candidates);
+    for (HighsInt m = 0; m < numMinus; ++m) {
+      if (isMinusRowRanged[m]) continue;
+      collectCandidate(step.minusHeaders[m].row, forcedNonBasicMinus[m],
+                       minusLowerSlack[m], minusUpperSlack[m],
+                       step.minusEntries[m], m, true, candidates);
+    }
+    // sort descending by non-basic support count
+    std::sort(candidates.begin(), candidates.end(),
+              [](const std::tuple<HighsInt, HighsInt, bool>& a,
+                 const std::tuple<HighsInt, HighsInt, bool>& b) {
+                return std::get<0>(a) > std::get<0>(b);
+              });
+    for (const auto& cand : candidates) {
+      HighsInt parentIndex = std::get<1>(cand);
+      if (std::get<2>(cand)) {
+        assignRowStatus(step.minusHeaders[parentIndex].row,
+                        minusLowerSlack[parentIndex],
+                        minusUpperSlack[parentIndex], basicAssigned,
+                        basicAssigned >= basicNeeded);
+      } else {
+        assignRowStatus(step.plusHeaders[parentIndex].row,
+                        plusLowerSlack[parentIndex],
+                        plusUpperSlack[parentIndex], basicAssigned,
+                        basicAssigned >= basicNeeded);
+      }
+    }
+  }
+}
+
 }  // namespace presolve
diff --git a/highs/presolve/HighsPostsolveStack.h b/highs/presolve/HighsPostsolveStack.h
index 60b3da9fdd0..da6e25dbbb5 100644
--- a/highs/presolve/HighsPostsolveStack.h
+++ b/highs/presolve/HighsPostsolveStack.h
@@ -18,6 +18,7 @@
 #include <cmath>
 #include <numeric>
 #include <tuple>
+#include <unordered_map>
 #include <vector>
 
 #include "lp_data/HConst.h"
@@ -46,6 +47,8 @@ class HighsPostsolveStack {
   // the constructor call, and should restore primal/dual solution values, as
   // well as the basis status as appropriate.
  public:
+  enum class OrigRowType : uint8_t { kOriginal, kCut, kAppended };
+
   enum class RowType {
     kGeq,
     kLeq,
@@ -59,6 +62,55 @@ class HighsPostsolveStack {
     Nonzero() = default;
   };
 
+  template <typename RowStorageFormat>
+  struct FmeRowData {
+    HighsInt row;
+    double rowLower;
+    double rowUpper;
+    HighsMatrixSlice<RowStorageFormat> rowVec;
+  };
+
+  struct FmeRowHeader {
+    HighsInt row;
+    double rowLower;
+    double rowUpper;
+  };
+
+  struct FmeStepHeader {
+    double colLower;
+    double colUpper;
+    HighsInt col;
+    HighsInt numPlus;
+    HighsInt numMinus;
+  };
+
+  struct FmeDescendant {
+    HighsInt row;
+    double scaleFactor;
+  };
+
+  struct FmeNewRow {
+    HighsInt row;
+    HighsInt plusParentIdx;   // index into plus parents (-1 if bound row)
+    HighsInt minusParentIdx;  // index into minus parents (-1 if bound row)
+  };
+
+  struct FmeAncestryEntry {
+    HighsInt step;
+    HighsInt parentRowIndex;
+    double scale;
+    bool isMinus;
+  };
+
+  struct FmeBlockStep {
+    HighsInt col;
+    double colLower;
+    double colUpper;
+    HighsInt numPlus;
+    HighsInt numMinus;
+    std::vector<FmeNewRow> newRows;
+  };
+
   size_t debug_prev_numreductions = 0;
   double debug_prev_col_lower = 0;
   double debug_prev_col_upper = 0;
@@ -78,6 +130,17 @@ class HighsPostsolveStack {
     void transformToPresolvedSpace(std::vector<double>& primalSol) const;
   };
 
+  struct FourierMotzkinObjCol {
+    double offset;
+    HighsInt col;
+
+    void transformToPresolvedSpace(const std::vector<Nonzero>& costEntries,
+                                   std::vector<double>& primalSol) const;
+
+    void undo(const std::vector<Nonzero>& costEntries,
+              HighsSolution& solution) const;
+  };
+
   struct FreeColSubstitution {
     double rhs;
     double colCost;
@@ -85,7 +148,8 @@ class HighsPostsolveStack {
     HighsInt col;
     RowType rowType;
 
-    void undo(const HighsOptions& options,
+    void undo(const HighsPostsolveStack& postsolveStack,
+              const HighsOptions& options,
               const std::vector<Nonzero>& rowValues,
               const std::vector<Nonzero>& colValues, HighsSolution& solution,
               HighsBasis& basis);
@@ -105,7 +169,8 @@ class HighsPostsolveStack {
     bool upperTightened;
     RowType rowType;
 
-    void undo(const HighsOptions& options,
+    void undo(const HighsPostsolveStack& postsolveStack,
+              const HighsOptions& options,
               const std::vector<Nonzero>& colValues, HighsSolution& solution,
               HighsBasis& basis) const;
   };
@@ -115,7 +180,8 @@ class HighsPostsolveStack {
     HighsInt addedEqRow;
     double eqRowScale;
 
-    void undo(const HighsOptions& options,
+    void undo(const HighsPostsolveStack& postsolveStack,
+              const HighsOptions& options,
               const std::vector<Nonzero>& eqRowValues, HighsSolution& solution,
               HighsBasis& basis) const;
   };
@@ -123,7 +189,8 @@ class HighsPostsolveStack {
   struct EqualityRowAdditions {
     HighsInt addedEqRow;
 
-    void undo(const HighsOptions& options,
+    void undo(const HighsPostsolveStack& postsolveStack,
+              const HighsOptions& options,
               const std::vector<Nonzero>& eqRowValues,
               const std::vector<Nonzero>& targetRows, HighsSolution& solution,
               HighsBasis& basis) const;
@@ -135,7 +202,8 @@ class HighsPostsolveStack {
     bool colLowerTightened;
     bool colUpperTightened;
 
-    void undo(const HighsOptions& options, HighsSolution& solution,
+    void undo(const HighsPostsolveStack& postsolveStack,
+              const HighsOptions& options, HighsSolution& solution,
               HighsBasis& basis) const;
   };
 
@@ -146,7 +214,8 @@ class HighsPostsolveStack {
     HighsInt col;
     HighsBasisStatus fixType;
 
-    void undo(const HighsOptions& options,
+    void undo(const HighsPostsolveStack& postsolveStack,
+              const HighsOptions& options,
               const std::vector<Nonzero>& colValues, HighsSolution& solution,
               HighsBasis& basis) const;
   };
@@ -154,16 +223,27 @@ class HighsPostsolveStack {
   struct RedundantRow {
     HighsInt row;
 
-    void undo(const HighsOptions& options, HighsSolution& solution,
+    void undo(const HighsPostsolveStack& postsolveStack,
+              const HighsOptions& options, HighsSolution& solution,
               HighsBasis& basis) const;
   };
 
+  struct ImpliedEquation {
+    HighsInt row;
+    bool atLower;
+
+    void undo(const HighsPostsolveStack& postsolveStack,
+              const std::vector<Nonzero>& rowValues,
+              HighsSolution& solution) const;
+  };
+
   struct ForcingRow {
     double side;
     HighsInt row;
     RowType rowType;
 
-    void undo(const HighsOptions& options,
+    void undo(const HighsPostsolveStack& postsolveStack,
+              const HighsOptions& options,
               const std::vector<Nonzero>& rowValues, HighsSolution& solution,
               HighsBasis& basis) const;
   };
@@ -175,7 +255,8 @@ class HighsPostsolveStack {
     bool atInfiniteUpper;
     bool colIntegral;
 
-    void undo(const HighsOptions& options,
+    void undo(const HighsPostsolveStack& postsolveStack,
+              const HighsOptions& options,
               const std::vector<Nonzero>& colValues, HighsSolution& solution,
               HighsBasis& basis) const;
   };
@@ -183,7 +264,8 @@ class HighsPostsolveStack {
   struct ForcingColumnRemovedRow {
     double rhs;
     HighsInt row;
-    void undo(const HighsOptions& options,
+    void undo(const HighsPostsolveStack& postsolveStack,
+              const HighsOptions& options,
               const std::vector<Nonzero>& rowValues, HighsSolution& solution,
               HighsBasis& basis) const;
   };
@@ -195,7 +277,8 @@ class HighsPostsolveStack {
     bool rowLowerTightened;
     bool rowUpperTightened;
 
-    void undo(const HighsOptions& options, HighsSolution& solution,
+    void undo(const HighsPostsolveStack& postsolveStack,
+              const HighsOptions& options, HighsSolution& solution,
               HighsBasis& basis) const;
   };
 
@@ -223,11 +306,31 @@ class HighsPostsolveStack {
     HighsInt row;
     HighsInt col;
 
-    void undo(const HighsOptions& options,
+    void undo(const HighsPostsolveStack& postsolveStack,
+              const HighsOptions& options,
               const std::vector<Nonzero>& rowValues, HighsSolution& solution,
               HighsBasis& basis);
   };
 
+  struct FmeStepData {
+    FmeStepHeader header;
+    std::vector<FmeRowHeader> plusHeaders;
+    std::vector<double> plusCoefs;
+    std::vector<std::vector<Nonzero>> plusEntries;
+    std::vector<std::vector<FmeDescendant>> plusDescendants;
+    std::vector<FmeRowHeader> minusHeaders;
+    std::vector<double> minusCoefs;
+    std::vector<std::vector<Nonzero>> minusEntries;
+    std::vector<std::vector<FmeDescendant>> minusDescendants;
+    std::vector<FmeNewRow> newRows;
+  };
+
+  static std::vector<FmeStepData> popFourierMotzkinBlock(HighsDataStack& stack);
+  static void undoFourierMotzkinBlock(const std::vector<FmeStepData>& steps,
+                                      const HighsOptions& options,
+                                      HighsSolution& solution,
+                                      HighsBasis& basis);
+
   /// tags for reduction
   enum class ReductionType : uint8_t {
     kLinearTransform,
@@ -238,69 +341,112 @@ class HighsPostsolveStack {
     kSingletonRow,
     kFixedCol,
     kRedundantRow,
+    kImpliedEquation,
     kForcingRow,
     kForcingColumn,
     kForcingColumnRemovedRow,
     kDuplicateRow,
     kDuplicateColumn,
     kSlackColSubstitution,
+    kFourierMotzkinBlock,
+    kFourierMotzkinObjCol,
   };
 
   HighsDataStack reductionValues;
   std::vector<std::pair<ReductionType, size_t>> reductions;
   std::vector<HighsInt> origColIndex;
   std::vector<HighsInt> origRowIndex;
+  std::vector<OrigRowType> origRowType;
   std::vector<HighsBool> linearlyTransformable;
 
   std::vector<Nonzero> rowValues;
   std::vector<Nonzero> colValues;
   HighsInt origNumCol = -1;
   HighsInt origNumRow = -1;
+  HighsInt numAppendedRows = 0;
+  HighsInt nextRowIndex = -1;
+  HighsInt nextColIndex = -1;
 
   void reductionAdded(ReductionType type) {
     size_t position = reductionValues.getCurrentDataSize();
     reductions.emplace_back(type, position);
   }
 
+  bool isModelRow(HighsInt row) const { return row < nextRowIndex; }
+
  public:
-  const HighsInt* getOrigRowsIndex() const { return origRowIndex.data(); }
+  const std::vector<HighsInt>& getOrigColIndex() const { return origColIndex; }
+  HighsInt getOrigColIndex(HighsInt col) const { return origColIndex[col]; }
 
-  const HighsInt* getOrigColsIndex() const { return origColIndex.data(); }
+  const std::vector<HighsInt>& getOrigRowIndex() const { return origRowIndex; }
+  HighsInt getOrigRowIndex(HighsInt row) const { return origRowIndex[row]; }
 
-  HighsInt getOrigRowIndex(HighsInt row) const {
-    assert(static_cast<size_t>(row) < origRowIndex.size());
-    return origRowIndex[row];
+  bool isOrigCol(HighsInt col) const { return origColIndex[col] < origNumCol; }
+
+  bool isOrigRow(HighsInt row) const {
+    return origRowType[row] == OrigRowType::kOriginal;
   }
 
-  HighsInt getOrigColIndex(HighsInt col) const {
-    assert(static_cast<size_t>(col) < origColIndex.size());
-    return origColIndex[col];
+  bool isAppendedRow(HighsInt row) const {
+    return origRowType[row] == OrigRowType::kAppended;
   }
 
-  void appendCutsToModel(HighsInt numCuts) {
+  bool isCutRow(HighsInt row) const {
+    return origRowType[row] == OrigRowType::kCut;
+  }
+
+  bool hasAppendedRows() const { return numAppendedRows > 0; }
+
+  void appendToModel(HighsInt& numRows, HighsInt numRowsToAppend,
+                     OrigRowType rowType) {
+    if (numRowsToAppend <= 0) return;
     size_t currNumRow = origRowIndex.size();
-    size_t newNumRow = currNumRow + numCuts;
+    size_t newNumRow = currNumRow + numRowsToAppend;
     origRowIndex.resize(newNumRow);
+    origRowType.resize(newNumRow, rowType);
     for (size_t i = currNumRow; i != newNumRow; ++i)
-      origRowIndex[i] = origNumRow++;
+      origRowIndex[i] = nextRowIndex++;
+    numRows += numRowsToAppend;
+  }
+
+  void appendCutsToModel(HighsInt numCuts) {
+    appendToModel(origNumRow, numCuts, OrigRowType::kCut);
+  }
+
+  void appendRowsToModel(HighsInt numRows) {
+    appendToModel(numAppendedRows, numRows, OrigRowType::kAppended);
   }
 
   void removeCutsFromModel(HighsInt numCuts) {
+    if (numCuts <= 0) return;
     origNumRow -= numCuts;
-
-    size_t origRowIndexSize = origRowIndex.size();
-    for (size_t i = origRowIndex.size(); i > 0; --i) {
-      if (origRowIndex[i - 1] < origNumRow) break;
-      --origRowIndexSize;
+    size_t newSize = 0;
+    for (size_t i = 0; i < origRowIndex.size(); ++i) {
+      if (origRowType[i] != OrigRowType::kCut) {
+        if (i != newSize) {
+          origRowIndex[newSize] = origRowIndex[i];
+          origRowType[newSize] = origRowType[i];
+        }
+        ++newSize;
+      }
     }
-
-    origRowIndex.resize(origRowIndexSize);
+    origRowIndex.resize(newSize);
+    origRowType.resize(newSize);
   }
 
   HighsInt getOrigNumRow() const { return origNumRow; }
 
   HighsInt getOrigNumCol() const { return origNumCol; }
 
+  HighsInt getNextRowIndex() const { return nextRowIndex; }
+
+  HighsInt getNextColIndex() const { return nextColIndex; }
+
+  void appendColToModel() {
+    origColIndex.push_back(nextColIndex++);
+    linearlyTransformable.push_back(false);
+  }
+
   void initializeIndexMaps(HighsInt numRow, HighsInt numCol);
 
   void compressIndexMaps(const std::vector<HighsInt>& newRowIndex,
@@ -308,8 +454,19 @@ class HighsPostsolveStack {
 
   void compressRowIndexMap(const std::vector<HighsInt>& newRowIndex);
   void compressColIndexMap(const std::vector<HighsInt>& newColIndex);
+
+  template <typename T>
   void compressIndexMap(const std::vector<HighsInt>& newIndex,
-                        std::vector<HighsInt>& origIndex);
+                        std::vector<T>& origIndex) {
+    size_t numEn = origIndex.size();
+    for (size_t i = 0; i != newIndex.size(); ++i) {
+      if (newIndex[i] == -1)
+        --numEn;
+      else
+        origIndex[newIndex[i]] = origIndex[i];
+    }
+    origIndex.resize(numEn);
+  }
 
   /// transform a column x by a linear mapping with a new column x'.
   /// I.e. substitute x = scale * x' + constant
@@ -478,6 +635,17 @@ class HighsPostsolveStack {
     reductionAdded(ReductionType::kRedundantRow);
   }
 
+  template <typename RowStorageFormat>
+  void impliedEquation(HighsInt row, bool atLower,
+                       const HighsMatrixSlice<RowStorageFormat>& rowVec) {
+    rowValues.clear();
+    for (const HighsSliceNonzero& rowVal : rowVec)
+      rowValues.emplace_back(origColIndex[rowVal.index()], rowVal.value());
+    reductionValues.push(ImpliedEquation{origRowIndex[row], atLower});
+    reductionValues.push(rowValues);
+    reductionAdded(ReductionType::kImpliedEquation);
+  }
+
   template <typename RowStorageFormat>
   void forcingRow(HighsInt row,
                   const HighsMatrixSlice<RowStorageFormat>& rowVec, double side,
@@ -520,6 +688,127 @@ class HighsPostsolveStack {
     reductionAdded(ReductionType::kForcingColumnRemovedRow);
   }
 
+  template <typename RowStorageFormat>
+  void fourierMotzkinBlockPushStep(
+      HighsInt col, const std::vector<FmeRowData<RowStorageFormat>>& plusRows,
+      const std::vector<FmeRowData<RowStorageFormat>>& minusRows) {
+    // push plus row entries
+    std::vector<FmeRowHeader> plusHeaders;
+    std::vector<double> plusCoefs;
+    plusHeaders.reserve(plusRows.size());
+    plusCoefs.reserve(plusRows.size());
+    for (const auto& rd : plusRows) {
+      std::vector<Nonzero> translated;
+      double coef = 0.0;
+      for (const HighsSliceNonzero& nz : rd.rowVec) {
+        if (nz.index() == col)
+          coef = nz.value();
+        else
+          translated.push_back({origColIndex[nz.index()], nz.value()});
+      }
+      reductionValues.push(translated);
+      plusCoefs.push_back(coef);
+      plusHeaders.push_back({origRowIndex[rd.row], rd.rowLower, rd.rowUpper});
+    }
+    reductionValues.push(plusCoefs);
+    reductionValues.push(plusHeaders);
+
+    // push minus row entries
+    std::vector<FmeRowHeader> minusHeaders;
+    std::vector<double> minusCoefs;
+    minusHeaders.reserve(minusRows.size());
+    minusCoefs.reserve(minusRows.size());
+    for (const auto& rd : minusRows) {
+      std::vector<Nonzero> translated;
+      double coef = 0.0;
+      for (const HighsSliceNonzero& nz : rd.rowVec) {
+        if (nz.index() == col)
+          coef = nz.value();
+        else
+          translated.push_back({origColIndex[nz.index()], nz.value()});
+      }
+      reductionValues.push(translated);
+      minusCoefs.push_back(coef);
+      minusHeaders.push_back({origRowIndex[rd.row], rd.rowLower, rd.rowUpper});
+    }
+    reductionValues.push(minusCoefs);
+    reductionValues.push(minusHeaders);
+  }
+
+  void fourierMotzkinBlockFinalise(
+      const std::vector<FmeBlockStep>& blockSteps,
+      const std::unordered_map<HighsInt, std::vector<FmeAncestryEntry>>&
+          rowAncestry) {
+    HighsInt numSteps = static_cast<HighsInt>(blockSteps.size());
+
+    // build K^j_i mapping from ancestry
+    std::vector<std::vector<std::vector<FmeDescendant>>> plusDescendantsAll(
+        numSteps);
+    std::vector<std::vector<std::vector<FmeDescendant>>> minusDescendantsAll(
+        numSteps);
+    for (HighsInt s = 0; s < numSteps; ++s) {
+      plusDescendantsAll[s].resize(blockSteps[s].numPlus);
+      minusDescendantsAll[s].resize(blockSteps[s].numMinus);
+    }
+    for (const auto& entry : rowAncestry) {
+      HighsInt row = entry.first;
+      HighsInt origRow = origRowIndex[row];
+      for (const auto& a : entry.second) {
+        if (a.isMinus)
+          minusDescendantsAll[a.step][a.parentRowIndex].push_back(
+              {origRow, a.scale});
+        else
+          plusDescendantsAll[a.step][a.parentRowIndex].push_back(
+              {origRow, a.scale});
+      }
+    }
+
+    // push descendants for each step's parents (plus then minus)
+    for (HighsInt s = 0; s < numSteps; ++s) {
+      assert(static_cast<HighsInt>(plusDescendantsAll[s].size()) ==
+             blockSteps[s].numPlus);
+      for (HighsInt p = 0; p < blockSteps[s].numPlus; ++p)
+        reductionValues.push(plusDescendantsAll[s][p]);
+      assert(static_cast<HighsInt>(minusDescendantsAll[s].size()) ==
+             blockSteps[s].numMinus);
+      for (HighsInt m = 0; m < blockSteps[s].numMinus; ++m)
+        reductionValues.push(minusDescendantsAll[s][m]);
+    }
+
+    // push new row origins for each step (translate row to orig space)
+    for (HighsInt s = 0; s < numSteps; ++s) {
+      std::vector<FmeNewRow> translated;
+      translated.reserve(blockSteps[s].newRows.size());
+      for (const auto& nr : blockSteps[s].newRows)
+        translated.push_back(
+            {origRowIndex[nr.row], nr.plusParentIdx, nr.minusParentIdx});
+      reductionValues.push(translated);
+    }
+
+    // push step headers
+    for (HighsInt s = 0; s < numSteps; ++s) {
+      FmeStepHeader header{blockSteps[s].colLower, blockSteps[s].colUpper,
+                           origColIndex[blockSteps[s].col],
+                           blockSteps[s].numPlus, blockSteps[s].numMinus};
+      reductionValues.push(header);
+    }
+
+    reductionValues.push(numSteps);
+    reductionAdded(ReductionType::kFourierMotzkinBlock);
+  }
+
+  void fourierMotzkinObjCol(HighsInt col, double offset,
+                            const std::vector<Nonzero>& costEntries) {
+    reductionValues.push(FourierMotzkinObjCol{offset, origColIndex[col]});
+    std::vector<Nonzero> translatedEntries;
+    translatedEntries.reserve(costEntries.size());
+    for (const Nonzero& entry : costEntries)
+      if (entry.index != col)
+        translatedEntries.emplace_back(origColIndex[entry.index], entry.value);
+    reductionValues.push(translatedEntries);
+    reductionAdded(ReductionType::kFourierMotzkinObjCol);
+  }
+
   void duplicateRow(HighsInt row, bool rowUpperTightened,
                     bool rowLowerTightened, HighsInt duplicateRow,
                     double duplicateRowScale) {
@@ -560,6 +849,7 @@ class HighsPostsolveStack {
   std::vector<double> getReducedPrimalSolution(
       const std::vector<double>& origPrimalSolution) {
     std::vector<double> reducedSolution = origPrimalSolution;
+    reducedSolution.resize(nextColIndex, 0.0);
 
     for (const std::pair<ReductionType, size_t>& primalColTransformation :
          reductions) {
@@ -578,6 +868,15 @@ class HighsPostsolveStack {
           linearTransform.transformToPresolvedSpace(reducedSolution);
           break;
         }
+        case ReductionType::kFourierMotzkinObjCol: {
+          reductionValues.setPosition(primalColTransformation.second);
+          std::vector<Nonzero> costEntries;
+          reductionValues.pop(costEntries);
+          FourierMotzkinObjCol fmObjCol;
+          reductionValues.pop(fmObjCol);
+          fmObjCol.transformToPresolvedSpace(costEntries, reducedSolution);
+          break;
+        }
         default:
           continue;
       }
@@ -632,8 +931,8 @@ class HighsPostsolveStack {
 
   /// undo presolve steps for primal dual solution and basis
   void undo(const HighsOptions& options, HighsSolution& solution,
-            HighsBasis& basis, const HighsInt report_col = -1,
-            const bool thread_safe = false) {
+            HighsBasis& basis, size_t numReductions = 0,
+            const HighsInt report_col = -1, const bool thread_safe = false) {
     HighsDataStack reductionValuesCopy;
     std::vector<Nonzero> colValuesCopy;
     std::vector<Nonzero> rowValuesCopy;
@@ -656,26 +955,26 @@ class HighsPostsolveStack {
     bool perform_basis_postsolve = basis.valid;
 
     // expand solution to original index space
-    assert(origNumCol > 0);
-    undoIterateBackwards(solution.col_value, origColIndex, origNumCol, 0.0);
+    assert(nextColIndex > 0);
+    undoIterateBackwards(solution.col_value, origColIndex, nextColIndex, 0.0);
 
-    assert(origNumRow >= 0);
-    undoIterateBackwards(solution.row_value, origRowIndex, origNumRow, 0.0);
+    assert(nextRowIndex >= 0);
+    undoIterateBackwards(solution.row_value, origRowIndex, nextRowIndex, 0.0);
 
     if (perform_dual_postsolve) {
       // if dual solution is given, expand dual solution and basis to original
       // index space
-      undoIterateBackwards(solution.col_dual, origColIndex, origNumCol, 0.0);
+      undoIterateBackwards(solution.col_dual, origColIndex, nextColIndex, 0.0);
 
-      undoIterateBackwards(solution.row_dual, origRowIndex, origNumRow, 0.0);
+      undoIterateBackwards(solution.row_dual, origRowIndex, nextRowIndex, 0.0);
     }
 
     if (perform_basis_postsolve) {
       // if basis is given, expand basis status values to original index space
-      undoIterateBackwards(basis.col_status, origColIndex, origNumCol,
+      undoIterateBackwards(basis.col_status, origColIndex, nextColIndex,
                            HighsBasisStatus::kNonbasic);
 
-      undoIterateBackwards(basis.row_status, origRowIndex, origNumRow,
+      undoIterateBackwards(basis.row_status, origRowIndex, nextRowIndex,
                            HighsBasisStatus::kNonbasic);
     }
 
@@ -733,7 +1032,11 @@ class HighsPostsolveStack {
       solutionLogging("After solving presolved LP");
     */
     // now undo the changes
-    for (size_t i = reductions.size(); i > 0; --i) {
+    for (size_t i = reductions.size(); i > numReductions; --i) {
+      if (report_col >= 0)
+        printf("Before  reduction %2d (type %2d): col_value[%2d] = %g\n",
+               int(i - 1), int(reductions[i - 1].first), int(report_col),
+               solution.col_value[report_col]);
       /*
       if (i - 1 == check_reduction) {
         printf("Checking reduction %d\n", int(check_reduction));
@@ -752,21 +1055,22 @@ class HighsPostsolveStack {
           reductionValues_.pop(colValues_);
           reductionValues_.pop(rowValues_);
           reductionValues_.pop(reduction);
-          reduction.undo(options, rowValues_, colValues_, solution, basis);
+          reduction.undo(*this, options, rowValues_, colValues_, solution,
+                         basis);
           break;
         }
         case ReductionType::kDoubletonEquation: {
           DoubletonEquation reduction;
           reductionValues_.pop(colValues_);
           reductionValues_.pop(reduction);
-          reduction.undo(options, colValues_, solution, basis);
+          reduction.undo(*this, options, colValues_, solution, basis);
           break;
         }
         case ReductionType::kEqualityRowAddition: {
           EqualityRowAddition reduction;
           reductionValues_.pop(rowValues_);
           reductionValues_.pop(reduction);
-          reduction.undo(options, rowValues_, solution, basis);
+          reduction.undo(*this, options, rowValues_, solution, basis);
           break;
         }
         case ReductionType::kEqualityRowAdditions: {
@@ -774,53 +1078,61 @@ class HighsPostsolveStack {
           reductionValues_.pop(colValues_);
           reductionValues_.pop(rowValues_);
           reductionValues_.pop(reduction);
-          reduction.undo(options, rowValues_, colValues_, solution, basis);
+          reduction.undo(*this, options, rowValues_, colValues_, solution,
+                         basis);
           break;
         }
         case ReductionType::kSingletonRow: {
           SingletonRow reduction;
           reductionValues_.pop(reduction);
-          reduction.undo(options, solution, basis);
+          reduction.undo(*this, options, solution, basis);
           break;
         }
         case ReductionType::kFixedCol: {
           FixedCol reduction;
           reductionValues_.pop(colValues_);
           reductionValues_.pop(reduction);
-          reduction.undo(options, colValues_, solution, basis);
+          reduction.undo(*this, options, colValues_, solution, basis);
           break;
         }
         case ReductionType::kRedundantRow: {
           RedundantRow reduction;
           reductionValues_.pop(reduction);
-          reduction.undo(options, solution, basis);
+          reduction.undo(*this, options, solution, basis);
+          break;
+        }
+        case ReductionType::kImpliedEquation: {
+          ImpliedEquation reduction;
+          reductionValues_.pop(rowValues_);
+          reductionValues_.pop(reduction);
+          reduction.undo(*this, rowValues_, solution);
           break;
         }
         case ReductionType::kForcingRow: {
           ForcingRow reduction;
           reductionValues_.pop(rowValues_);
           reductionValues_.pop(reduction);
-          reduction.undo(options, rowValues_, solution, basis);
+          reduction.undo(*this, options, rowValues_, solution, basis);
           break;
         }
         case ReductionType::kForcingColumn: {
           ForcingColumn reduction;
           reductionValues_.pop(colValues_);
           reductionValues_.pop(reduction);
-          reduction.undo(options, colValues_, solution, basis);
+          reduction.undo(*this, options, colValues_, solution, basis);
           break;
         }
         case ReductionType::kForcingColumnRemovedRow: {
           ForcingColumnRemovedRow reduction;
           reductionValues_.pop(rowValues_);
           reductionValues_.pop(reduction);
-          reduction.undo(options, rowValues_, solution, basis);
+          reduction.undo(*this, options, rowValues_, solution, basis);
           break;
         }
         case ReductionType::kDuplicateRow: {
           DuplicateRow reduction;
           reductionValues_.pop(reduction);
-          reduction.undo(options, solution, basis);
+          reduction.undo(*this, options, solution, basis);
           break;
         }
         case ReductionType::kDuplicateColumn: {
@@ -833,7 +1145,20 @@ class HighsPostsolveStack {
           SlackColSubstitution reduction;
           reductionValues_.pop(rowValues_);
           reductionValues_.pop(reduction);
-          reduction.undo(options, rowValues_, solution, basis);
+          reduction.undo(*this, options, rowValues_, solution, basis);
+          break;
+        }
+        case ReductionType::kFourierMotzkinBlock: {
+          auto steps = popFourierMotzkinBlock(reductionValues_);
+          undoFourierMotzkinBlock(steps, options, solution, basis);
+          break;
+        }
+        case ReductionType::kFourierMotzkinObjCol: {
+          std::vector<Nonzero> costEntries;
+          reductionValues_.pop(costEntries);
+          FourierMotzkinObjCol reduction;
+          reductionValues_.pop(reduction);
+          reduction.undo(costEntries, solution);
           break;
         }
         default:
@@ -845,6 +1170,14 @@ class HighsPostsolveStack {
     }
     //    if (report_col >= 0) reportColLogging(-2);
 
+    solution.col_value.resize(origNumCol);
+    if (perform_dual_postsolve) solution.col_dual.resize(origNumCol);
+    if (perform_basis_postsolve) basis.col_status.resize(origNumCol);
+
+    solution.row_value.resize(origNumRow);
+    if (perform_dual_postsolve) solution.row_dual.resize(origNumRow);
+    if (perform_basis_postsolve) basis.row_status.resize(origNumRow);
+
 #ifdef DEBUG_EXTRA
     // solution should not contain NaN or Inf
     assert(!containsNanOrInf(solution.col_value));
@@ -865,7 +1198,7 @@ class HighsPostsolveStack {
     HighsBasis basis;
     basis.valid = false;
     solution.dual_valid = false;
-    undo(options, solution, basis, report_col, thread_safe);
+    undo(options, solution, basis, 0, report_col, thread_safe);
   }
 
   /*
@@ -881,158 +1214,6 @@ class HighsPostsolveStack {
   }
   */
 
-  // Only used for debugging
-  void undoUntil(const HighsOptions& options, HighsSolution& solution,
-                 HighsBasis& basis, size_t numReductions) {
-    reductionValues.resetPosition();
-
-    // Do these returns ever happen? How is it known that undo has not
-    // been performed?
-    assert(solution.col_value.size() == origColIndex.size());
-    assert(solution.row_value.size() == origRowIndex.size());
-    // This should be a better measure of whether undo can be
-    // performed
-    assert(solution.value_valid);
-    if (solution.col_value.size() != origColIndex.size()) return;
-    if (solution.row_value.size() != origRowIndex.size()) return;
-
-    bool perform_dual_postsolve = solution.dual_valid;
-    assert((solution.col_dual.size() == solution.col_value.size()) ==
-           perform_dual_postsolve);
-    bool perform_basis_postsolve = basis.valid;
-
-    // expand solution to original index space
-    undoIterateBackwards(solution.col_value, origColIndex, origNumCol, 0.0);
-
-    undoIterateBackwards(solution.row_value, origRowIndex, origNumRow, 0.0);
-
-    if (perform_dual_postsolve) {
-      // if dual solution is given, expand dual solution and basis to original
-      // index space
-      undoIterateBackwards(solution.col_dual, origColIndex, origNumCol, 0.0);
-
-      undoIterateBackwards(solution.row_dual, origRowIndex, origNumRow, 0.0);
-    }
-
-    if (perform_basis_postsolve) {
-      // if basis is given, expand basis status values to original index space
-      undoIterateBackwards(basis.col_status, origColIndex, origNumCol,
-                           HighsBasisStatus::kNonbasic);
-
-      undoIterateBackwards(basis.row_status, origRowIndex, origNumRow,
-                           HighsBasisStatus::kNonbasic);
-    }
-
-    // now undo the changes
-    for (size_t i = reductions.size(); i > numReductions; --i) {
-      switch (reductions[i - 1].first) {
-        case ReductionType::kLinearTransform: {
-          LinearTransform reduction;
-          reductionValues.pop(reduction);
-          reduction.undo(options, solution);
-          break;
-        }
-        case ReductionType::kFreeColSubstitution: {
-          FreeColSubstitution reduction;
-          reductionValues.pop(colValues);
-          reductionValues.pop(rowValues);
-          reductionValues.pop(reduction);
-          reduction.undo(options, rowValues, colValues, solution, basis);
-          break;
-        }
-        case ReductionType::kDoubletonEquation: {
-          DoubletonEquation reduction;
-          reductionValues.pop(colValues);
-          reductionValues.pop(reduction);
-          reduction.undo(options, colValues, solution, basis);
-          break;
-        }
-        case ReductionType::kEqualityRowAddition: {
-          EqualityRowAddition reduction;
-          reductionValues.pop(rowValues);
-          reductionValues.pop(reduction);
-          reduction.undo(options, rowValues, solution, basis);
-          break;
-        }
-        case ReductionType::kEqualityRowAdditions: {
-          EqualityRowAdditions reduction;
-          reductionValues.pop(colValues);
-          reductionValues.pop(rowValues);
-          reductionValues.pop(reduction);
-          reduction.undo(options, rowValues, colValues, solution, basis);
-          break;
-        }
-        case ReductionType::kSingletonRow: {
-          SingletonRow reduction;
-          reductionValues.pop(reduction);
-          reduction.undo(options, solution, basis);
-          break;
-        }
-        case ReductionType::kFixedCol: {
-          FixedCol reduction;
-          reductionValues.pop(colValues);
-          reductionValues.pop(reduction);
-          reduction.undo(options, colValues, solution, basis);
-          break;
-        }
-        case ReductionType::kRedundantRow: {
-          RedundantRow reduction;
-          reductionValues.pop(reduction);
-          reduction.undo(options, solution, basis);
-          break;
-        }
-        case ReductionType::kForcingRow: {
-          ForcingRow reduction;
-          reductionValues.pop(rowValues);
-          reductionValues.pop(reduction);
-          reduction.undo(options, rowValues, solution, basis);
-          break;
-        }
-        case ReductionType::kForcingColumn: {
-          ForcingColumn reduction;
-          reductionValues.pop(colValues);
-          reductionValues.pop(reduction);
-          reduction.undo(options, colValues, solution, basis);
-          break;
-        }
-        case ReductionType::kForcingColumnRemovedRow: {
-          ForcingColumnRemovedRow reduction;
-          reductionValues.pop(rowValues);
-          reductionValues.pop(reduction);
-          reduction.undo(options, rowValues, solution, basis);
-          break;
-        }
-        case ReductionType::kDuplicateRow: {
-          DuplicateRow reduction;
-          reductionValues.pop(reduction);
-          reduction.undo(options, solution, basis);
-          break;
-        }
-        case ReductionType::kDuplicateColumn: {
-          DuplicateColumn reduction;
-          reductionValues.pop(reduction);
-          reduction.undo(options, solution, basis);
-          break;
-        }
-        case ReductionType::kSlackColSubstitution: {
-          SlackColSubstitution reduction;
-          reductionValues.pop(rowValues);
-          reductionValues.pop(reduction);
-          reduction.undo(options, rowValues, solution, basis);
-          break;
-        }
-      }
-    }
-#ifdef DEBUG_EXTRA
-    // solution should not contain NaN or Inf
-    assert(!containsNanOrInf(solution.col_value));
-    // row values are not determined by postsolve
-    // assert(!containsNanOrInf(solution.row_value));
-    assert(!containsNanOrInf(solution.col_dual));
-    assert(!containsNanOrInf(solution.row_dual));
-#endif
-  }
-
   size_t numReductions() const { return reductions.size(); }
 };