Merge pull request #2004 from ERGO-Code/fix-1999

Construct standard form LP

check/TestLpSolvers.cpp

@@ -2,7 +2,7 @@
 #include "Highs.h"
 #include "catch.hpp"
 
-const bool dev_run = false;
+const bool dev_run = true;  // false;
 
 struct IterationCount {
   HighsInt simplex;
@@ -283,7 +283,7 @@ TEST_CASE("LP-solver", "[highs_lp_solver]") {
   const HighsInfo& info = highs.getInfo();
   REQUIRE(info.num_dual_infeasibilities == 0);
 
-  REQUIRE(info.simplex_iteration_count == 472);  // 476);  // 444);
+  REQUIRE(info.simplex_iteration_count == 472);
 
   HighsModelStatus model_status = highs.getModelStatus();
   REQUIRE(model_status == HighsModelStatus::kOptimal);
@@ -296,7 +296,7 @@ TEST_CASE("LP-solver", "[highs_lp_solver]") {
   return_status = highs.run();
   REQUIRE(return_status == HighsStatus::kOk);
 
-  REQUIRE(info.simplex_iteration_count == 592);  // 621);  // 584);  //
+  REQUIRE(info.simplex_iteration_count == 592);
 }
 
 TEST_CASE("mip-with-lp-solver", "[highs_lp_solver]") {
@@ -470,3 +470,117 @@ TEST_CASE("blending-lp-ipm", "[highs_lp_solver]") {
     printf("Sum   dual infeasibilities = %g\n", info.sum_dual_infeasibilities);
   }
 }
+
+void testStandardForm(const HighsLp& lp) {
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+  HighsInt sense = HighsInt(lp.sense_);
+  highs.passModel(lp);
+  highs.run();
+  //  highs.writeSolution("", kSolutionStylePretty);
+  double required_objective_function_value =
+      highs.getInfo().objective_function_value;
+
+  HighsInt num_col;
+  HighsInt num_row;
+  HighsInt num_nz;
+  double offset;
+  REQUIRE(highs.getStandardFormLp(num_col, num_row, num_nz, offset) ==
+          HighsStatus::kOk);
+
+  std::vector<double> cost(num_col);
+  std::vector<double> rhs(num_row);
+  std::vector<HighsInt> start(num_col + 1);
+  std::vector<HighsInt> index(num_nz);
+  std::vector<double> value(num_nz);
+  REQUIRE(highs.getStandardFormLp(num_col, num_row, num_nz, offset, cost.data(),
+                                  rhs.data(), start.data(), index.data(),
+                                  value.data()) == HighsStatus::kOk);
+
+  HighsLp standard_form_lp;
+  standard_form_lp.num_col_ = num_col;
+  standard_form_lp.num_row_ = num_row;
+  standard_form_lp.offset_ = offset;
+  standard_form_lp.col_cost_ = cost;
+  standard_form_lp.col_lower_.assign(num_col, 0);
+  standard_form_lp.col_upper_.assign(num_col, kHighsInf);
+  standard_form_lp.row_lower_ = rhs;
+  standard_form_lp.row_upper_ = rhs;
+  standard_form_lp.a_matrix_.start_ = start;
+  standard_form_lp.a_matrix_.index_ = index;
+  standard_form_lp.a_matrix_.value_ = value;
+  REQUIRE(highs.passModel(standard_form_lp) == HighsStatus::kOk);
+  // highs.writeModel("");
+  REQUIRE(highs.run() == HighsStatus::kOk);
+  REQUIRE(highs.getModelStatus() == HighsModelStatus::kOptimal);
+  highs.writeSolution("", kSolutionStylePretty);
+  double objective_function_value =
+      sense * highs.getInfo().objective_function_value;
+  double objective_difference =
+      std::fabs(objective_function_value - required_objective_function_value) /
+      std::max(1.0, std::fabs(required_objective_function_value));
+  REQUIRE(objective_difference < 1e-10);
+  const bool look_at_presolved_lp = false;
+  if (look_at_presolved_lp) {
+    // Strange that presolve doesn't convert the constraints
+    //
+    // Ax-s = b; s >= 0 into Ax >= b
+    REQUIRE(highs.passModel(standard_form_lp) == HighsStatus::kOk);
+    REQUIRE(highs.presolve() == HighsStatus::kOk);
+    HighsLp presolved_lp = highs.getPresolvedLp();
+    REQUIRE(highs.passModel(presolved_lp) == HighsStatus::kOk);
+    highs.writeModel("");
+  }
+}
+
+void testStandardFormModel(const std::string model) {
+  const std::string model_file =
+      std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
+  ;
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+  highs.readModel(model_file);
+  HighsLp lp = highs.getLp();
+  testStandardForm(lp);
+}
+
+TEST_CASE("standard-form-mps", "[highs_lp_solver]") {
+  testStandardFormModel("avgas");
+  testStandardFormModel("afiro");
+}
+
+TEST_CASE("standard-form-lp", "[highs_lp_solver]") {
+  HighsLp lp;
+  lp.offset_ = -0.5;
+  lp.num_col_ = 4;
+  lp.num_row_ = 3;
+  lp.col_cost_ = {1, 1, 1, -1};
+  lp.col_lower_ = {1, -kHighsInf, -kHighsInf, -1};
+  lp.col_upper_ = {kHighsInf, kHighsInf, 2, 3};
+  lp.row_lower_ = {0, 1, -kHighsInf};
+  lp.row_upper_ = {4, kHighsInf, 4};
+  lp.a_matrix_.start_ = {0, 2, 4, 6, 8};
+  lp.a_matrix_.index_ = {0, 2, 0, 1, 1, 2, 0, 2};
+  lp.a_matrix_.value_ = {1, 1, 1, 1, 1, 1, 1, 1};
+
+  testStandardForm(lp);
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+
+  std::vector<HighsInt> index;
+  std::vector<double> value;
+  // Add a fixed column and a fixed row, and maximize
+  highs.passModel(lp);
+  index = {0, 1, 2};
+  value = {-1, 1, -1};
+  REQUIRE(highs.addCol(-2.0, 1.0, 1.0, 3, index.data(), value.data()) ==
+          HighsStatus::kOk);
+  index = {0, 1, 2, 3};
+  value = {-2, -1, 1, 3};
+  REQUIRE(highs.addRow(1.0, 1.0, 4, index.data(), value.data()) ==
+          HighsStatus::kOk);
+  REQUIRE(highs.changeObjectiveSense(ObjSense::kMaximize) == HighsStatus::kOk);
+  printf(
+      "\nNow test by adding a fixed column and a fixed row, and maximizing\n");
+  testStandardForm(highs.getLp());
+}

src/Highs.h

@@ -400,6 +400,17 @@ class Highs {
    * Methods for model output
    */
 
+  /**
+   * @brief Identify and the standard form of the HighsLp instance in
+   * HiGHS
+   */
+  HighsStatus getStandardFormLp(HighsInt& num_col, HighsInt& num_row,
+                                HighsInt& num_nz, double& offset,
+                                double* cost = nullptr, double* rhs = nullptr,
+                                HighsInt* start = nullptr,
+                                HighsInt* index = nullptr,
+                                double* value = nullptr);
+
   /**
    * @brief Return a const reference to the presolved HighsLp instance in HiGHS
    */
@@ -1378,6 +1389,12 @@ class Highs {
       HighsPresolveStatus::kNotPresolved;
   HighsModelStatus model_status_ = HighsModelStatus::kNotset;
 
+  bool standard_form_valid_;
+  double standard_form_offset_;
+  std::vector<double> standard_form_cost_;
+  std::vector<double> standard_form_rhs_;
+  HighsSparseMatrix standard_form_matrix_;
+
   HEkk ekk_instance_;
 
   HighsPresolveLog presolve_log_;
@@ -1430,6 +1447,9 @@ class Highs {
   // Clears the presolved model and its status
   void clearPresolve();
   //
+  // Clears the standard form LP
+  void clearStandardFormLp();
+  //
   // Methods to clear solver data for users in Highs class members
   // before (possibly) updating them with data from trying to solve
   // the incumbent model.
@@ -1473,6 +1493,7 @@ class Highs {
   void underDevelopmentLogMessage(const std::string& method_name);
 
   // Interface methods
+  HighsStatus formStandardFormLp();
   HighsStatus basisForSolution();
   HighsStatus addColsInterface(
       HighsInt ext_num_new_col, const double* ext_col_cost,