Extracted effective-costs changes from fix-initial-sweep-interrupt

Author: jajhall

check/TestMipSolver.cpp

@@ -2033,5 +2033,6 @@ TEST_CASE("dominated-column-double-fixing", "[highs_test_mip_solver]") {
   highs.setOptionValue("output_flag", dev_run);
   REQUIRE(highs.passModel(lp) == HighsStatus::kOk);
   solve(highs, kHighsOffString, HighsModelStatus::kInfeasible);
+  highs.setOptionValue("output_flag", dev_run);
   solve(highs, kHighsOnString, HighsModelStatus::kInfeasible);
 }

check/TestPresolve.cpp

@@ -1241,3 +1241,56 @@ TEST_CASE("issue-3140", "[highs_test_presolve]") {
 
   highs.resetGlobalScheduler(true);
 }
+
+TEST_CASE("presolve-light-no-crossover", "[highs_test_presolve]") {
+  Highs h;
+  h.setOptionValue("output_flag", dev_run);
+  std::string source = "instances";  //"miplib2017";
+  std::string model = "adlittle";    //"germanrr";  //"neos-2746589-doon";//
+  std::string type = "mps";
+  std::string filename =
+      (source == "instances" ? std::string(HIGHS_DIR) + "/check/instances/"
+                             : "/srv/" + source) +
+      "/" + model + "." + type + (source == "instances" ? "" : ".gz");
+
+  REQUIRE(h.readModel(filename) == HighsStatus::kOk);
+  const HighsLp& lp = h.getLp();
+
+  REQUIRE(h.setOptionValue("presolve_light", kHighsOnString) ==
+          HighsStatus::kOk);
+  const std::string solver = kIpxString;  // kHipoString;  //
+  REQUIRE(h.setOptionValue("solver", solver) == HighsStatus::kOk);
+  REQUIRE(h.setOptionValue("run_crossover", kHighsOffString) ==
+          HighsStatus::kOk);
+  REQUIRE(h.setOptionValue("solve_relaxation", true) == HighsStatus::kOk);
+  h.run();
+
+  REQUIRE(h.getModelStatus() == HighsModelStatus::kOptimal);
+
+  h.resetGlobalScheduler(true);
+}
+
+/*
+TEST_CASE("test-fuzzing", "[highs_test_presolve]") {
+  Highs h;
+  //  h.setOptionValue("output_flag", dev_run);
+  h.setOptionValue("presolve_rule_logging", true);
+  h.setOptionValue("log_dev_level", 1);
+
+  const std::string model = "issue-008";
+  std::string model_file = std::string(HIGHS_DIR) + "/build/OscarFuzzing/" +
+                           model + "/" + model + ".mps";
+
+  REQUIRE(h.readModel(model_file) == HighsStatus::kOk);
+
+  std::string options_file =
+      std::string(HIGHS_DIR) + "/build/OscarFuzzing/" + model + "/options.txt";
+  REQUIRE(h.readOptions(options_file) == HighsStatus::kOk);
+
+  h.writeOptions("", true);
+
+  h.run();
+
+  h.resetGlobalScheduler(true);
+}
+*/

check/TestPresolveRules.cpp

@@ -26,6 +26,10 @@ TEST_CASE("test-col-stuffing", "[highs_test_presolve_rules]") {
   Highs h;
   h.setOptionValue("output_flag", dev_run);
   h.setOptionValue("presolve_rule_test", kPresolveRuleColStuffing);
+  REQUIRE(h.setOptionValue("presolve_rule_logging", true) == HighsStatus::kOk);
+  // Initial sweep doesn't yield reductions, but switch it off for clarity
+  REQUIRE(h.setOptionValue("presolve_rule_off",
+                           1 << kPresolveRuleInitialSweep) == HighsStatus::kOk);
   const bool lp0 = true;
   const bool lp1 = true;
   const bool lp1a = true;
@@ -90,6 +94,67 @@ TEST_CASE("test-col-stuffing", "[highs_test_presolve_rules]") {
   h.resetGlobalScheduler(true);
 }
 
+/*
+TEST_CASE("test-weakly-dominated-col-upper", "[highs_test_presolve_rules]") {
+  Highs h;
+  h.setOptionValue("output_flag", dev_run);
+  REQUIRE(h.setOptionValue("presolve_rule_logging", true) == HighsStatus::kOk);
+  // LP is
+  //
+  // min -y, subject to x+y <= 0, x >= 0; 0 <= x <= 1, y free
+  //
+  // Optimal solution is x = 1; y = -1, with x nonbasic with dual -1, and
+  HighsLp lp;
+  lp.num_col_ = 2;
+  lp.num_row_ = 2;
+  lp.col_lower_ = {-kHighsInf, -kHighsInf};
+  lp.col_upper_ = {1,  kHighsInf};
+  lp.row_lower_ = {-kHighsInf,         1};
+  lp.row_upper_ = {         0, kHighsInf};
+  lp.a_matrix_.format_ = MatrixFormat::kRowwise;
+  lp.a_matrix_.start_ = {0, 2, 3};
+  lp.a_matrix_.index_ = {0, 1, 0};
+  lp.a_matrix_.value_ = {1, 1, 1};
+
+  bool maximize_first = true;
+  std::string sense_string = "";
+  std::string test_string = "";
+
+  for (HighsInt k = 0; k < 2; k++) {
+    // Passes are minimize c^Tx and maximize -c^Tx according to
+    // maximize_first
+    if (maximize_first) {
+      lp.sense_ = ObjSense::kMaximize;
+      sense_string = "maximize";
+      lp.col_cost_ = {0, 1};
+    } else {
+      lp.sense_ = ObjSense::kMinimize;
+      sense_string = "minimize";
+      lp.col_cost_ = {0, -1};
+    }
+    //  REQUIRE(h.setOptionValue("presolve_rule_test", 0) == HighsStatus::kOk);
+    //  test_string = "vanilla-presolve-" + sense_string;
+    //  presolveOffOn(test_string, lp, h);
+
+    REQUIRE(h.setOptionValue("presolve_rule_test",
+kPresolveRuleWeaklyDominatedColUpper) == HighsStatus::kOk);
+
+    // test_string = "initial-sweep+test-weakly-dominated-col-upper-" +
+sense_string;
+    // presolveOffOn(test_string, lp, h, 1, 1, 1);
+
+    REQUIRE(h.setOptionValue("presolve_rule_off", 1 <<
+kPresolveRuleInitialSweep) == HighsStatus::kOk);
+
+    test_string = "test-weakly-dominated-col-upper-" + sense_string;
+    presolveOffOn(test_string, lp, h, 1, 2, 1);
+
+    maximize_first = !maximize_first;
+  }
+  h.resetGlobalScheduler(true);
+}
+*/
+
 TEST_CASE("test-parallel-rows-cut-ordering", "[highs_test_presolve_rules]") {
   // Rows 0 and 1 are parallel (both [1, 1]). Row 0 is marked as a
   // cut. detectParallelRowsAndCols must remove the cut row (0) and
@@ -113,6 +178,7 @@ TEST_CASE("test-parallel-rows-cut-ordering", "[highs_test_presolve_rules]") {
 
   HighsOptions options;
   options.presolve_rule_test = kPresolveRuleParallelRowsAndCols;
+  options.presolve_rule_off = 1 << kPresolveRuleInitialSweep;
   options.output_flag = dev_run;
 
   HighsTimer timer;
@@ -137,6 +203,117 @@ TEST_CASE("test-parallel-rows-cut-ordering", "[highs_test_presolve_rules]") {
   REQUIRE(!postsolve_stack.isCutRow(0));
 }
 
+TEST_CASE("test-effective-costs", "[highs_test_presolve]") {
+  // Debugging ZeroObjSingletonContinuousCol for germanrr highlighted
+  // the deficiency in computing the active_cost_norm when the
+  // objective is f = z, with z = c^Tx and z free. In
+  // HighsSolution.cpp is the method getEffectiveCosts that
+  // substitutes all free column singletons into the objective to get
+  // the "effective costs".
+  Highs h;
+  h.setOptionValue("output_flag", dev_run);
+  bool test_all = true;
+  bool test_lp0 = test_all;
+  bool test_lp1 = test_all;
+  bool test_lp2 = test_all;
+
+  if (test_lp0) {
+    HighsLp lp;
+    // First LP is
+    //
+    // min 4z
+    //
+    // -1 <=    x + y - 2z <= 1
+    //
+    // -1 <= 201x + y      <= 1
+    //
+    // 0 <= x <= 1, y, z free
+    //
+    // where the bounds on the two constraints and non-unit
+    // coefficients of z in the objective and first contraint give
+    // code coverage
+    //
+    // Aiming to minimize 4z, and bound is given by 2z >= x + y - 1,
+    // so substitute z = (x+y-1)/2 into the objective to give
+    //
+    // min 2x + 2y - 2
+    //
+    // y is then minimized with bound is given by y >= -201x - 1, so
+    // substitute y = -201x - 1 into the objective to give
+    //
+    // min 2x +(-402x-2) - 2 = -400x - 4
+    //
+    // This function is minimized when x = 1 to give y = -202 and z =
+    // -101 with objective -404
+    //
+    // The optimal dual values are -400 for x, -2 for row 0 and 2 for
+    // row 1. However, although this example tests code coverage on
+    // identifying free column singletons and a double free column
+    // singleton identified in getEffectiveCosts, the dual of -400 for
+    // the only nonbasic column means that there are no active costs,
+    // so active_cost_norm is zero (hence absolute and relative dual
+    // infeasibility measures are identical).
+    lp.model_name_ = "LP0";
+    lp.num_col_ = 3;
+    lp.num_row_ = 2;
+    lp.col_cost_ = {0, 0, 4};
+    lp.col_lower_ = {0, -kHighsInf, -kHighsInf};
+    lp.col_upper_ = {1, kHighsInf, kHighsInf};
+    lp.a_matrix_.format_ = MatrixFormat::kRowwise;
+    lp.a_matrix_.start_ = {0, 3, 5};
+    lp.a_matrix_.index_ = {0, 1, 2, 0, 1};
+    lp.a_matrix_.value_ = {1, 1, -2, 201, 1};
+    lp.row_lower_ = {-1, -1};
+    lp.row_upper_ = {1, 1};
+    h.passModel(lp);
+    h.setOptionValue("log_dev_level", 1);
+    h.setOptionValue("presolve_rule_logging", kHighsOnString);
+    h.run();
+    REQUIRE(h.getInfo().active_cost_norm == 0);
+  }
+  if (test_lp1) {
+    HighsLp lp;
+    // Here's a simpler example that reflects the behaviour observed
+    // with germanrr, where the cost row of the matrix introduced many
+    // large costs. Hence the presolved model had a large value for
+    // active_cost_norm but, after postsolve, the model had
+    // active_cost_norm = 1.
+
+    double cost = 1e5;
+    double eps = 1e-4;
+    lp.model_name_ = "LP1";
+    lp.num_col_ = 3;
+    lp.num_row_ = 2;
+    lp.col_cost_ = {0, 0, 1};
+    lp.col_lower_ = {0, 0, -kHighsInf};
+    lp.col_upper_ = {1, 1, kHighsInf};
+    lp.a_matrix_.format_ = MatrixFormat::kRowwise;
+    lp.a_matrix_.start_ = {0, 3, 5};
+    lp.a_matrix_.index_ = {0, 1, 2, 0, 1};
+    lp.a_matrix_.value_ = {cost, cost - eps, 1, 1, 1, 1};
+    lp.row_lower_ = {0, 1};
+    lp.row_upper_ = {0, 1};
+    h.passModel(lp);
+
+    h.run();
+    REQUIRE(h.getInfo().active_cost_norm == cost);
+  }
+  if (test_lp2) {
+    // Finally gas11 has 61 free column singletons: 55 in the first
+    // pass, and 6 in the second.
+    const std::string model = "gas11";
+    std::string model_file =
+        std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
+    REQUIRE(h.readModel(model_file) == HighsStatus::kWarning);
+    REQUIRE(h.setOptionValue(kPresolveString, kHighsOffString) ==
+            HighsStatus::kOk);
+    HighsStatus return_status = h.run();
+    REQUIRE(return_status == HighsStatus::kOk);
+    double active_cost_norm = 2.000000001e+7;
+    REQUIRE(std::fabs(h.getInfo().active_cost_norm - active_cost_norm) <= 1e-8);
+  }
+}
+
 TEST_CASE("test-fourier-motzkin", "[highs_test_presolve_rules]") {
   Highs h;
   h.setOptionValue("output_flag", dev_run);

cmake/sources.cmake

@@ -417,6 +417,7 @@ set(highs_sources
     presolve/HighsSymmetry.cpp
     presolve/HPresolve.cpp
     presolve/HPresolveAnalysis.cpp
+    presolve/HPresolveDebug.cpp
     presolve/HPresolveTest.cpp
     presolve/ICrash.cpp
     presolve/ICrashUtil.cpp

highs/lp_data/HighsInfo.cpp

@@ -53,6 +53,7 @@ void HighsInfo::invalidatePrimalKkt() {
   num_complementarity_violations = kHighsIllegalComplementarityCount;
   max_complementarity_violation = kHighsIllegalComplementarityViolation;
   primal_dual_objective_error = kHighsIllegalComplementarityViolation;
+  active_bound_norm = kHighsIllegalDoubleMeasure;
 }
 
 void HighsInfo::invalidateDualKkt() {
@@ -69,6 +70,7 @@ void HighsInfo::invalidateDualKkt() {
   num_complementarity_violations = kHighsIllegalComplementarityCount;
   max_complementarity_violation = kHighsIllegalComplementarityViolation;
   primal_dual_objective_error = kHighsIllegalComplementarityViolation;
+  active_cost_norm = kHighsIllegalDoubleMeasure;
 }
 
 bool HighsInfo::equal(const HighsInfo& info_) const {
@@ -149,6 +151,8 @@ bool HighsInfo::equal(const HighsInfo& info_) const {
     return false;
   if (info_.primal_dual_objective_error != this->primal_dual_objective_error)
     return false;
+  if (info_.active_cost_norm != this->active_cost_norm) return false;
+  if (info_.active_bound_norm != this->active_bound_norm) return false;
   return true;
 }
 

highs/lp_data/HighsInfo.h

@@ -129,6 +129,8 @@ struct HighsInfoStruct {
   double max_complementarity_violation;
   double primal_dual_objective_error;
   double primal_dual_integral;
+  double active_cost_norm;
+  double active_bound_norm;
 };
 
 class HighsInfo : public HighsInfoStruct {
@@ -386,6 +388,15 @@ class HighsInfo : public HighsInfoStruct {
         new InfoRecordDouble("primal_dual_integral", "Primal-dual integral",
                              advanced, &primal_dual_integral, 0);
     records.push_back(record_double);
+
+    record_double = new InfoRecordDouble("active_cost_norm", "Active cost norm",
+                                         advanced, &active_cost_norm, 0);
+    records.push_back(record_double);
+
+    record_double =
+        new InfoRecordDouble("active_bound_norm", "Active bound norm", advanced,
+                             &active_bound_norm, 0);
+    records.push_back(record_double);
   }
 
  public: