Added double computeBasisCondition() { return computeBasisCondition(t…

…his->lp_); } and now using passed HighsLp; Formatted
ERGO-Code · Sep 20, 2024 · c6028d8 · c6028d8
1 parent 091ab68
commit c6028d8
Show file tree

Hide file tree

Showing 6 changed files with 45 additions and 26 deletions.
diff --git a/check/TestBasisSolves.cpp b/check/TestBasisSolves.cpp
@@ -566,31 +566,48 @@ TEST_CASE("Basis-solves", "[highs_basis_solves]") {
 }
 
 TEST_CASE("Kappa", "[highs_basis_solves]") {
-  std::string filename;
-  filename = std::string(HIGHS_DIR) + "/check/instances/chip.mps";
-  // filename = std::string(HIGHS_DIR) + "/check/instances/avgas.mps";
-  // filename = std::string(HIGHS_DIR) + "/check/instances/adlittle.mps";
-  //  filename = std::string(HIGHS_DIR) + "/check/instances/25fv47.mps";
+  // chip optimal basis matrix is B=[1, 2; 1, 4] with
+  //
+  // ||B||_1=6; ||B^{-1}||_1=5/2 so kappa_1 = 15
+  //
+  // ||B||_2=4.56; ||B^{-1}||_inf=2.355 so kappa_inf = 10.9
+  //
+  // ||B||_inf=6; ||B^{-1}||_inf=5/2 so kappa_inf = 15
+  //
+  double chip_kappa = 15;
+
+  std::string model;
+  model = "chip";
+  // model = "avgas";
+  // model = "adlittle";
+  std::string filename =
+      std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
 
   Highs highs;
-  //highs.setOptionValue("output_flag", dev_run);
+  highs.setOptionValue("output_flag", dev_run);
 
- // Read the LP given by filename
+  // Read the LP given by filename
   highs.readModel(filename);
 
   double kappa;
   REQUIRE(highs.getKappa(kappa) == HighsStatus::kError);
 
   highs.run();
-  
+
   REQUIRE(highs.getKappa(kappa) == HighsStatus::kOk);
+  if (dev_run)
+    printf("highs.getKappa for %s yields %g\n", model.c_str(), kappa);
+
+  if (model == "chip") REQUIRE(std::fabs(kappa - chip_kappa) < 1e-4);
 
   highs.clearModel();
-
-  filename = std::string(HIGHS_DIR) + "/check/instances/flugpl.mps";
-
-  highs.run();
-
-  REQUIRE(highs.getKappa(kappa) == HighsStatus::kError);
 
+  const bool test_mip = true;
+  if (test_mip) {
+    model = "flugpl";
+    filename = std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
+    highs.readModel(filename);
+    highs.run();
+    REQUIRE(highs.getKappa(kappa) == HighsStatus::kError);
+  }
 }
diff --git a/src/lp_data/Highs.cpp b/src/lp_data/Highs.cpp
@@ -1948,12 +1948,13 @@ HighsStatus Highs::getReducedColumn(const HighsInt col, double* col_vector,
 }
 
 HighsStatus Highs::getKappa(double& kappa) {
-  printf("Highs::getKappa basis_.valid = %d, ekk_instance_.status_.has_invert = %d\n",
-	 int(basis_.valid),
-	 int(ekk_instance_.status_.has_invert));
+  printf(
+      "Highs::getKappa basis_.valid = %d, ekk_instance_.status_.has_invert = "
+      "%d\n",
+      int(basis_.valid), int(ekk_instance_.status_.has_invert));
   if (!ekk_instance_.status_.has_invert)
     return invertRequirementError("getBasisInverseRow");
-  kappa = ekk_instance_.computeBasisCondition();
+  kappa = ekk_instance_.computeBasisCondition(this->model_.lp_);
   return HighsStatus::kOk;
 }
 

diff --git a/src/lp_data/HighsOptions.h b/src/lp_data/HighsOptions.h
@@ -555,7 +555,7 @@ struct HighsOptionsStruct {
 #endif
         mip_improving_solution_save(false),
         mip_improving_solution_report_sparse(false),
-        mip_improving_solution_file("") {};
+        mip_improving_solution_file(""){};
 };
 
 // For now, but later change so HiGHS properties are string based so that new

diff --git a/src/simplex/HEkk.cpp b/src/simplex/HEkk.cpp
@@ -3634,18 +3634,18 @@ HighsStatus HEkk::returnFromSolve(const HighsStatus return_status) {
   return return_status;
 }
 
-double HEkk::computeBasisCondition() {
-  HighsInt solver_num_row = lp_.num_row_;
-  HighsInt solver_num_col = lp_.num_col_;
+double HEkk::computeBasisCondition(const HighsLp& lp) {
+  HighsInt solver_num_row = lp.num_row_;
+  HighsInt solver_num_col = lp.num_col_;
   vector<double> bs_cond_x;
   vector<double> bs_cond_y;
   vector<double> bs_cond_z;
   vector<double> bs_cond_w;
   HVector row_ep;
   row_ep.setup(solver_num_row);
 
-  const HighsInt* Astart = lp_.a_matrix_.start_.data();
-  const double* Avalue = lp_.a_matrix_.value_.data();
+  const HighsInt* Astart = lp.a_matrix_.start_.data();
+  const double* Avalue = lp.a_matrix_.value_.data();
   // Compute the Hager condition number estimate for the basis matrix
   const double expected_density = 1;
   bs_cond_x.resize(solver_num_row);

diff --git a/src/simplex/HEkk.h b/src/simplex/HEkk.h
@@ -134,7 +134,8 @@ class HEkk {
   HighsBasis getHighsBasis(HighsLp& use_lp) const;
 
   const SimplexBasis& getSimplexBasis() { return basis_; }
-  double computeBasisCondition();
+  double computeBasisCondition(const HighsLp& lp);
+  double computeBasisCondition() { return computeBasisCondition(this->lp_); }
 
   HighsStatus initialiseSimplexLpBasisAndFactor(
       const bool only_from_known_basis = false);

diff --git a/src/util/HFactor.h b/src/util/HFactor.h
@@ -132,7 +132,7 @@ class HFactor {
         build_timer_(nullptr),
         nwork(0),
         u_merit_x(0),
-        u_total_x(0) {};
+        u_total_x(0){};
 
   /**
    * @brief Copy problem size and pointers of constraint matrix, and set