Standard form of avgas formed correctly

ERGO-Code · Oct 25, 2024 · 9a14344 · 9a14344
1 parent 1d17f81
commit 9a14344
Show file tree

Hide file tree

Showing 7 changed files with 276 additions and 55 deletions.
diff --git a/check/TestLpSolvers.cpp b/check/TestLpSolvers.cpp
@@ -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,64 @@ TEST_CASE("blending-lp-ipm", "[highs_lp_solver]") {
     printf("Sum   dual infeasibilities = %g\n", info.sum_dual_infeasibilities);
   }
 }
+
+TEST_CASE("standard-form-lp", "[highs_lp_solver]") {
+  std::string model;
+  std::string model_file;
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+  model = "avgas";
+  model_file = std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
+  highs.readModel(model_file);
+  highs.run();
+  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.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);
+  if (dev_run) highs.writeModel("");
+  REQUIRE(highs.run() == HighsStatus::kOk);
+  REQUIRE(highs.getModelStatus() == HighsModelStatus::kOptimal);
+  double objective_function_value = 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("");
+  }
+}
diff --git a/src/Highs.h b/src/Highs.h
@@ -404,7 +404,12 @@ class Highs {
    * @brief Identify and the standard form of the HighsLp instance in
    * HiGHS
    */
-  HighsStatus getStandardFormLp(HighsLp& standard_form_lp);
+  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
@@ -1384,6 +1389,8 @@ 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_;
@@ -1439,10 +1446,10 @@ 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.

diff --git a/src/lp_data/Highs.cpp b/src/lp_data/Highs.cpp
@@ -1638,10 +1638,34 @@ HighsStatus Highs::run() {
   return returnFromRun(return_status, undo_mods);
 }
 
-HighsStatus Highs::getStandardFormLp(HighsLp& standard_form_lp) {
-  HighsStatus status = formStandardFormLp();
-  if (status != HighsStatus::kOk) return status;
-  standard_form_lp = this->standard_form_lp_;
+HighsStatus Highs::getStandardFormLp(HighsInt& num_col, HighsInt& num_row,
+                                     HighsInt& num_nz, double offset,
+                                     double* cost, double* rhs, HighsInt* start,
+                                     HighsInt* index, double* value) {
+  if (!this->standard_form_valid_) {
+    HighsStatus status = formStandardFormLp();
+    if (status != HighsStatus::kOk) return status;
+  }
+  num_col = this->standard_form_cost_.size();
+  num_row = this->standard_form_rhs_.size();
+  num_nz = this->standard_form_matrix_.start_[num_col];
+  offset = this->standard_form_offset_;
+  for (HighsInt iCol = 0; iCol < num_col; iCol++) {
+    if (cost) cost[iCol] = this->standard_form_cost_[iCol];
+    if (start) start[iCol] = this->standard_form_matrix_.start_[iCol];
+    if (index || value) {
+      for (HighsInt iEl = this->standard_form_matrix_.start_[iCol];
+           iEl < this->standard_form_matrix_.start_[iCol + 1]; iEl++) {
+        if (index) index[iEl] = this->standard_form_matrix_.index_[iEl];
+        if (value) value[iEl] = this->standard_form_matrix_.value_[iEl];
+      }
+    }
+  }
+  if (start) start[num_col] = this->standard_form_matrix_.start_[num_col];
+  if (rhs) {
+    for (HighsInt iRow = 0; iRow < num_row; iRow++)
+      rhs[iRow] = this->standard_form_rhs_[iRow];
+  }
   return HighsStatus::kOk;
 }
 
@@ -2401,7 +2425,7 @@ HighsStatus Highs::changeObjectiveSense(const ObjSense sense) {
     model_.lp_.sense_ = sense;
     // Nontrivial change
     clearPresolve();
-  clearStandardFormLp();
+    clearStandardFormLp();
     invalidateModelStatusSolutionAndInfo();
   }
   return returnFromHighs(HighsStatus::kOk);
@@ -3451,6 +3475,8 @@ void Highs::clearPresolve() {
 }
 
 void Highs::clearStandardFormLp() {
+  standard_form_valid_ = false;
+  standard_form_offset_ = 0;
   standard_form_cost_.clear();
   standard_form_rhs_.clear();
   standard_form_matrix_.clear();