Adding code for testing scaling and condition estimates

Author: jajhall

check/TestBasisSolves.cpp

@@ -615,3 +615,20 @@ TEST_CASE("Kappa", "[highs_basis_solves]") {
 
   highs.resetGlobalScheduler(true);
 }
+
+TEST_CASE("scaling-kappa", "[highs_basis_solves]") {
+  std::string model;
+  // model = "chip";
+  //  model = "avgas";
+  model = "afiro";
+  std::string filename =
+      std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
+
+  Highs highs;
+  //  highs.setOptionValue("output_flag", dev_run);
+
+  // Read the LP given by filename
+  highs.readModel(filename);
+
+  highs.run();
+}

highs/Highs.h

@@ -690,11 +690,9 @@ class Highs {
 
   /**
    * @brief Get the condition number of the current basis matrix,
-   * possibly computing it exactly and reporting the error in the
-   * approximate condition number
+   * possibly computing it exactly
    */
-  HighsStatus getKappa(double& kappa, const bool exact = false,
-                       const bool report = false) const;
+  HighsStatus getKappa(double& kappa, const bool exact = false) const;
 
   /**
    * @brief Get the number of columns in the incumbent model

highs/interfaces/highs_c_api.h

@@ -784,7 +784,8 @@ HighsInt Highs_getDoubleOptionValue(const void* highs, const char* option,
  * @param highs     A pointer to the Highs instance.
  * @param option    The name of the option.
  * @param value     A pointer to allocated memory to store the current value of
- *                  the option. This must have length `kHighsMaximumStringLength`.
+ *                  the option. This must have length
+ * `kHighsMaximumStringLength`.
  *
  * @returns A `kHighsStatus` constant indicating whether the call succeeded.
  */

highs/lp_data/HConst.h

@@ -17,6 +17,8 @@
 
 #include "util/HighsInt.h"
 
+const bool kSimplexScaleDev = true;
+
 const std::string kHighsCopyrightStatement =
     "Copyright (c) 2025 HiGHS under MIT licence terms";
 

highs/lp_data/HStruct.h

@@ -100,6 +100,8 @@ struct HighsScale {
   double cost;
   std::vector<double> col;
   std::vector<double> row;
+  void print(const std::string& prefix = "",
+             const std::string& message = "") const;
 };
 
 struct HighsLpMods {

highs/lp_data/Highs.cpp

@@ -2338,11 +2338,10 @@ HighsStatus Highs::getReducedColumn(const HighsInt col, double* col_vector,
   return HighsStatus::kOk;
 }
 
-HighsStatus Highs::getKappa(double& kappa, const bool exact,
-                            const bool report) const {
+HighsStatus Highs::getKappa(double& kappa, const bool exact) const {
   if (!ekk_instance_.status_.has_invert)
-    return invertRequirementError("getBasisInverseRow");
-  kappa = ekk_instance_.computeBasisCondition(this->model_.lp_, exact, report);
+    return invertRequirementError("getKappa");
+  kappa = ekk_instance_.computeBasisCondition(this->model_.lp_, exact);
   return HighsStatus::kOk;
 }
 

highs/lp_data/HighsLpUtils.cpp

@@ -3077,3 +3077,29 @@ bool HighsUserScaleData::scaleWarning(std::string& message) const {
   message = ss.str();
   return true;
 }
+
+void HighsScale::print(const std::string& prefix,
+                       const std::string& message) const {
+  double min_row_scale = kHighsInf;
+  double max_row_scale = -kHighsInf;
+  double min_col_scale = kHighsInf;
+  double max_col_scale = -kHighsInf;
+  HighsInt num_row = this->row.size();
+  HighsInt num_col = this->col.size();
+  for (HighsInt iCol = 0; iCol < num_col; iCol++) {
+    min_col_scale = std::min(this->col[iCol], min_col_scale);
+    max_col_scale = std::max(this->col[iCol], max_col_scale);
+  }
+  for (HighsInt iRow = 0; iRow < num_row; iRow++) {
+    min_row_scale = std::min(this->row[iRow], min_row_scale);
+    max_row_scale = std::max(this->row[iRow], max_row_scale);
+  }
+  printf(
+      "%sTxt HighsScale: Cost scale = %g; col scale in [%g, %g]; row scale in "
+      "[%g, %g]; %s\n",
+      prefix.c_str(), this->cost, min_col_scale, max_col_scale, min_row_scale,
+      max_row_scale, message.c_str());
+  printf("%sCsv,%s,%g,%g,%g,%g,%g\n", prefix.c_str(), message.c_str(),
+         this->cost, min_col_scale, max_col_scale, min_row_scale,
+         max_row_scale);
+}

highs/simplex/HApp.h

@@ -154,7 +154,9 @@ inline HighsStatus solveLpSimplex(HighsLpSolverObject& solver_object) {
   // considering computing scaling factors if there are none - and
   // then move to EKK
   considerSimplexScaling(options, incumbent_lp);
-  //
+  if (kSimplexScaleDev)
+    incumbent_lp.scale_.print("grepSimplexScaling",
+                              "After scaling, " + incumbent_lp.model_name_);
   const bool was_scaled = incumbent_lp.is_scaled_;
   if (!status.has_basis && !basis.valid && basis.useful) {
     // There is no simplex basis, but there is a useful HiGHS basis
@@ -218,10 +220,10 @@ inline HighsStatus solveLpSimplex(HighsLpSolverObject& solver_object) {
   bool solve_unscaled_lp = false;
   bool solved_unscaled_lp = false;
   if (!incumbent_lp.scale_.has_scaling) {
+    solve_unscaled_lp = true;
     //
     // Solve the unscaled LP with unscaled NLA
     //
-    solve_unscaled_lp = true;
     return_status = ekk_instance.solve();
     solved_unscaled_lp = true;
     ekk_instance.unpermute();
@@ -447,6 +449,7 @@ inline HighsStatus solveLpSimplex(HighsLpSolverObject& solver_object) {
     assert(basis.valid);
     highs_info.basis_validity = kBasisValidityValid;
   }
+  if (kSimplexScaleDev) ekk_instance.testBasisCondition();
   // Move the incumbent LP back from Ekk
   incumbent_lp = std::move(ekk_lp);
   incumbent_lp.is_moved_ = false;

highs/simplex/HEkk.cpp

@@ -3603,8 +3603,22 @@ HighsStatus HEkk::returnFromSolve(const HighsStatus return_status) {
   return return_status;
 }
 
-double HEkk::computeBasisCondition(const HighsLp& lp, const bool exact,
-                                   const bool report) const {
+void HEkk::testBasisCondition(const std::string& message) const {
+  double exact_kappa_tt = -timer_->read();
+  bool exact = true;
+  double exact_kappa = this->computeBasisCondition(this->lp_, exact);
+  exact_kappa_tt += timer_->read();
+  double approx_kappa_tt = -timer_->read();
+  exact = false;
+  double approx_kappa = this->computeBasisCondition(this->lp_, exact);
+  approx_kappa_tt += timer_->read();
+  highsLogUser(options_->log_options, HighsLogType::kInfo,
+               "getKappa,%s,%g,%g,%g,%g,%s,%s\n", this->lp_.model_name_.c_str(),
+               exact_kappa, exact_kappa_tt, approx_kappa, approx_kappa_tt,
+               message.c_str(), this->lp_.origin_name_.c_str());
+}
+
+double HEkk::computeBasisCondition(const HighsLp& lp, const bool exact) const {
   HighsInt solver_num_row = lp.num_row_;
   HighsInt solver_num_col = lp.num_col_;
   vector<double> bs_cond_x;
@@ -3616,9 +3630,10 @@ double HEkk::computeBasisCondition(const HighsLp& lp, const bool exact,
 
   const HighsInt* Astart = lp.a_matrix_.start_.data();
   const double* Avalue = lp.a_matrix_.value_.data();
-  double exact_norm_Binv = 0;
+  double norm_Binv = 0;
   if (exact) {
-    // Compute the exact norm of B^{-1}
+    // Compute the exact 1-norm of B^{-1}
+    norm_Binv = 0;
     for (HighsInt r_n = 0; r_n < solver_num_row; r_n++) {
       row_ep.clear();
       row_ep.index[row_ep.count] = r_n;
@@ -3630,79 +3645,80 @@ double HEkk::computeBasisCondition(const HighsLp& lp, const bool exact,
       double c_norm = 0.0;
       for (HighsInt iX = 0; iX < row_ep.count; iX++)
         c_norm += std::fabs(row_ep.array[row_ep.index[iX]]);
-      exact_norm_Binv = std::max(c_norm, exact_norm_Binv);
-    }
-  }
-  // Compute the Hager condition number estimate for the basis matrix
-  const double expected_density = 1;
-  bs_cond_x.resize(solver_num_row);
-  bs_cond_y.resize(solver_num_row);
-  bs_cond_z.resize(solver_num_row);
-  bs_cond_w.resize(solver_num_row);
-  // x = ones(n,1)/n;
-  // y = A\x;
-  double mu = 1.0 / solver_num_row;
-  double norm_Binv = 0;
-  for (HighsInt r_n = 0; r_n < solver_num_row; r_n++) bs_cond_x[r_n] = mu;
-  row_ep.clear();
-  for (HighsInt r_n = 0; r_n < solver_num_row; r_n++) {
-    double value = bs_cond_x[r_n];
-    if (value) {
-      row_ep.index[row_ep.count] = r_n;
-      row_ep.array[r_n] = value;
-      row_ep.count++;
-    }
-  }
-  for (HighsInt ps_n = 1; ps_n <= 5; ps_n++) {
-    row_ep.packFlag = false;
-    simplex_nla_.ftran(row_ep, expected_density);
-
-    // zeta = sign(y);
-    for (HighsInt r_n = 0; r_n < solver_num_row; r_n++) {
-      bs_cond_y[r_n] = row_ep.array[r_n];
-      if (bs_cond_y[r_n] > 0)
-        bs_cond_w[r_n] = 1.0;
-      else if (bs_cond_y[r_n] < 0)
-        bs_cond_w[r_n] = -1.0;
-      else
-        bs_cond_w[r_n] = 0.0;
+      norm_Binv = std::max(c_norm, norm_Binv);
     }
-    // z=A'\zeta;
+  } else {
+    // Compute the Hager 1-norm condition number estimate for the basis matrix
+    const double expected_density = 1;
+    bs_cond_x.resize(solver_num_row);
+    bs_cond_y.resize(solver_num_row);
+    bs_cond_z.resize(solver_num_row);
+    bs_cond_w.resize(solver_num_row);
+    // x = ones(n,1)/n;
+    // y = A\x;
+    double mu = 1.0 / solver_num_row;
+    norm_Binv = 0;
+    for (HighsInt r_n = 0; r_n < solver_num_row; r_n++) bs_cond_x[r_n] = mu;
     row_ep.clear();
     for (HighsInt r_n = 0; r_n < solver_num_row; r_n++) {
-      double value = bs_cond_w[r_n];
+      double value = bs_cond_x[r_n];
       if (value) {
         row_ep.index[row_ep.count] = r_n;
         row_ep.array[r_n] = value;
         row_ep.count++;
       }
     }
-    row_ep.packFlag = false;
-    simplex_nla_.btran(row_ep, expected_density);
-    double norm_z = 0.0;
-    double ztx = 0.0;
-    norm_Binv = 0.0;
-    HighsInt argmax_z = -1;
-    for (HighsInt r_n = 0; r_n < solver_num_row; r_n++) {
-      bs_cond_z[r_n] = row_ep.array[r_n];
-      double abs_z_v = fabs(bs_cond_z[r_n]);
-      if (abs_z_v > norm_z) {
-        norm_z = abs_z_v;
-        argmax_z = r_n;
+    for (HighsInt ps_n = 1; ps_n <= 5; ps_n++) {
+      row_ep.packFlag = false;
+      simplex_nla_.ftran(row_ep, expected_density);
+      // zeta = sign(y);
+      for (HighsInt r_n = 0; r_n < solver_num_row; r_n++) {
+        bs_cond_y[r_n] = row_ep.array[r_n];
+        if (bs_cond_y[r_n] > 0)
+          bs_cond_w[r_n] = 1.0;
+        else if (bs_cond_y[r_n] < 0)
+          bs_cond_w[r_n] = -1.0;
+        else
+          bs_cond_w[r_n] = 0.0;
+      }
+      // z=A'\zeta;
+      row_ep.clear();
+      for (HighsInt r_n = 0; r_n < solver_num_row; r_n++) {
+        double value = bs_cond_w[r_n];
+        if (value) {
+          row_ep.index[row_ep.count] = r_n;
+          row_ep.array[r_n] = value;
+          row_ep.count++;
+        }
       }
-      ztx += bs_cond_z[r_n] * bs_cond_x[r_n];
-      norm_Binv += fabs(bs_cond_y[r_n]);
+      row_ep.packFlag = false;
+      simplex_nla_.btran(row_ep, expected_density);
+      double norm_z = 0.0;
+      double ztx = 0.0;
+      norm_Binv = 0.0;
+      HighsInt argmax_z = -1;
+      for (HighsInt r_n = 0; r_n < solver_num_row; r_n++) {
+        bs_cond_z[r_n] = row_ep.array[r_n];
+        double abs_z_v = fabs(bs_cond_z[r_n]);
+        if (abs_z_v > norm_z) {
+          norm_z = abs_z_v;
+          argmax_z = r_n;
+        }
+        ztx += bs_cond_z[r_n] * bs_cond_x[r_n];
+        norm_Binv += fabs(bs_cond_y[r_n]);
+      }
+      if (norm_z <= ztx) break;
+      // x = zeros(n,1);
+      // x(fd_i) = 1;
+      for (HighsInt r_n = 0; r_n < solver_num_row; r_n++) bs_cond_x[r_n] = 0.0;
+      row_ep.clear();
+      row_ep.count = 1;
+      row_ep.index[0] = argmax_z;
+      row_ep.array[argmax_z] = 1.0;
+      bs_cond_x[argmax_z] = 1.0;
     }
-    if (norm_z <= ztx) break;
-    // x = zeros(n,1);
-    // x(fd_i) = 1;
-    for (HighsInt r_n = 0; r_n < solver_num_row; r_n++) bs_cond_x[r_n] = 0.0;
-    row_ep.clear();
-    row_ep.count = 1;
-    row_ep.index[0] = argmax_z;
-    row_ep.array[argmax_z] = 1.0;
-    bs_cond_x[argmax_z] = 1.0;
   }
+  // Compute the exact 1-norm of B
   double norm_B = 0.0;
   for (HighsInt r_n = 0; r_n < solver_num_row; r_n++) {
     HighsInt vr_n = basis_.basicIndex_[r_n];
@@ -3714,19 +3730,7 @@ double HEkk::computeBasisCondition(const HighsLp& lp, const bool exact,
       c_norm += 1.0;
     norm_B = max(c_norm, norm_B);
   }
-  double cond_B = norm_Binv * norm_B;
-  double exact_cond_B = exact_norm_Binv * norm_B;
-  if (exact) {
-    assert(exact_norm_Binv > 0);
-    if (report)
-      highsLogUser(
-          options_->log_options, HighsLogType::kInfo,
-          "HEkk::computeBasisCondition: grep_kappa model,||B||_1,approx "
-          "||B^{-1}||_1,approx_kappa,||B^{-1}||_1,kappa = ,%s,%g,%g,%g,%g,%g\n",
-          lp.model_name_.c_str(), norm_B, norm_Binv, cond_B, exact_norm_Binv,
-          exact_cond_B);
-    return exact_cond_B;
-  }
+  double cond_B = norm_B * norm_Binv;
   return cond_B;
 }
 

highs/simplex/HEkk.h

@@ -127,11 +127,10 @@ class HEkk {
   HighsBasis getHighsBasis(HighsLp& use_lp) const;
 
   const SimplexBasis& getSimplexBasis() { return basis_; }
-  double computeBasisCondition(const HighsLp& lp, const bool exact = false,
-                               const bool report = false) const;
-  double computeBasisCondition() const {
-    return computeBasisCondition(this->lp_, false, false);
-  }
+
+  void testBasisCondition(const std::string& message = "") const;
+  double computeBasisCondition(const HighsLp& lp,
+                               const bool exact = false) const;
 
   HighsStatus initialiseSimplexLpBasisAndFactor(
       const bool only_from_known_basis = false);