Adding cost scaling without breaking unit tests!

jajhall · jajhall · commit 28245d6098b7 · 2025-10-28T14:09:22.000Z
diff --git a/check/TestBasisSolves.cpp b/check/TestBasisSolves.cpp
@@ -617,19 +617,47 @@ TEST_CASE("Kappa", "[highs_basis_solves]") {
 }
 
 TEST_CASE("scaling-kappa", "[highs_basis_solves]") {
-  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("simplex_scale_strategy", kSimplexScaleStrategyMaxValue);
+  highs.setOptionValue("log_dev_level", kHighsLogDevLevelInfo);
+  const bool model_test = false;
+  if (model_test) {
+    std::string model;
+    // model = "chip";
+    //  model = "avgas";
+    model = "adlittle";
+    std::string filename =
+      std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
+    // Read the LP given by filename
+    highs.readModel(filename);
+    highs.run();
 
-  // Read the LP given by filename
-  highs.readModel(filename);
-
-  highs.setOptionValue("simplex_scale_strategy", kSimplexScaleStrategyMaxValue0157);
+  /*
+  highs.clearSolver();
+  highs.setOptionValue("simplex_scale_strategy", kSimplexScaleStrategyMaxValueMatrixAndCost);
   highs.run();
+  */
+  }
+  const bool lp_test = false;
+  if (lp_test) {
+    const double cost = 1e6;
+    HighsLp lp;
+    lp.num_col_ = 2;
+    lp.num_row_ = 2;
+    lp.offset_ = 1e-4;
+    lp.col_cost_ = {-4.1 * cost, -7.3 * cost};
+    lp.col_lower_ = {0, 0};
+    lp.col_upper_ = {10, 10};
+    lp.row_lower_ = {-kHighsInf, -2};
+    lp.row_upper_ = {6, kHighsInf};
+    lp.a_matrix_.start_ = {0, 2, 4};
+    lp.a_matrix_.index_ = {0, 1, 0, 1};
+    lp.a_matrix_.value_ = {0.95, 0.98, 0.97, -1.01};
+    highs.passModel(lp);
+    highs.setOptionValue("presolve", kHighsOffString);
+    highs.setOptionValue("simplex_scale_strategy",
+			 kSimplexScaleStrategyMaxValueMatrixAndCost);
+    highs.run();
+  }
 }
diff --git a/check/TestLpValidation.cpp b/check/TestLpValidation.cpp
@@ -108,48 +108,6 @@ TEST_CASE("LP-dimension-validation", "[highs_data]") {
   lp.a_matrix_.index_[1] = 1;
   REQUIRE(highs.passModel(lp) == HighsStatus::kOk);
 
-  /*
-  if (dev_run) printf("Give valid scale_.strategy\n");
-  lp.scale_.strategy = kSimplexScaleStrategyOff;
-  REQUIRE(highs.passModel(lp) == HighsStatus::kError);
-
-  if (dev_run) printf("Give valid scale_.num_col\n");
-  lp.scale_.num_col = 0;
-  REQUIRE(highs.passModel(lp) == HighsStatus::kError);
-
-  if (dev_run) printf("Give valid scale_.num_row\n");
-  lp.scale_.num_row = 0;
-  REQUIRE(highs.passModel(lp) == HighsStatus::kError);
-
-  if (dev_run) printf("Give valid scale_.col.size()\n");
-  lp.scale_.col.resize(0);
-  REQUIRE(highs.passModel(lp) == HighsStatus::kError);
-
-  if (dev_run) printf("Give valid scale_.row.size()\n");
-  lp.scale_.row.resize(0);
-  REQUIRE(highs.passModel(lp) == HighsStatus::kOk);
-
-  if (dev_run) printf("Set scale_.strategy =
-  kSimplexScaleStrategyMaxValue015\n"); lp.scale_.strategy =
-  kSimplexScaleStrategyMaxValue015; REQUIRE(highs.passModel(lp) ==
-  HighsStatus::kError);
-
-  if (dev_run) printf("Give valid scale_.num_col\n");
-  lp.scale_.num_col = true_num_col;
-  REQUIRE(highs.passModel(lp) == HighsStatus::kError);
-
-  if (dev_run) printf("Give valid scale_.num_row\n");
-  lp.scale_.num_row = true_num_row;
-  REQUIRE(highs.passModel(lp) == HighsStatus::kError);
-
-  if (dev_run) printf("Give valid scale_.col.size()\n");
-  lp.scale_.col.resize(true_num_col);
-  REQUIRE(highs.passModel(lp) == HighsStatus::kError);
-
-  if (dev_run) printf("Give valid scale_.row.size()\n");
-  lp.scale_.row.resize(true_num_row);
-  REQUIRE(highs.passModel(lp) == HighsStatus::kOk);
-  */
 }
 
 TEST_CASE("LP-validation", "[highs_data]") {
diff --git a/highs/lp_data/HConst.h b/highs/lp_data/HConst.h
@@ -17,7 +17,8 @@
 
 #include "util/HighsInt.h"
 
-const bool kSimplexScaleDev = false;
+const bool kSimplexScaleDev = true;
+const bool kSimplexScaleDevReport = kSimplexScaleDev && false;
 
 const std::string kHighsCopyrightStatement =
     "Copyright (c) 2025 HiGHS under MIT licence terms";
@@ -55,9 +56,8 @@ enum SimplexScaleStrategy {
   kSimplexScaleStrategyEquilibration,                   // 2
   kSimplexScaleStrategyForcedEquilibration,             // 3
   kSimplexScaleStrategyMaxValue,                        // 4
-  kSimplexScaleStrategyMaxValue015 = kSimplexScaleStrategyMaxValue,
-  kSimplexScaleStrategyMaxValue0157 = kSimplexScaleStrategyMaxValue,
-  kSimplexScaleStrategyMax = kSimplexScaleStrategyMaxValue
+  kSimplexScaleStrategyMaxValueMatrixAndCost,           // 5
+  kSimplexScaleStrategyMax = kSimplexScaleStrategyMaxValueMatrixAndCost
 };
 
 enum HighsDebugLevel {
diff --git a/highs/lp_data/HighsLp.cpp b/highs/lp_data/HighsLp.cpp
@@ -281,10 +281,12 @@ void HighsLp::unapplyScale() {
   assert(scale.has_scaling);
   // Unapply the scaling to the bounds, costs and matrix, and record
   // that it has been unapplied
+  assert(scale.cost);
+  const double unscale_cost_mu = 1.0 / scale.cost;
   for (HighsInt iCol = 0; iCol < this->num_col_; iCol++) {
     this->col_lower_[iCol] *= scale.col[iCol];
     this->col_upper_[iCol] *= scale.col[iCol];
-    this->col_cost_[iCol] /= scale.col[iCol];
+    this->col_cost_[iCol] /= (unscale_cost_mu * scale.col[iCol]);
   }
   for (HighsInt iRow = 0; iRow < this->num_row_; iRow++) {
     this->row_lower_[iRow] /= scale.row[iRow];
diff --git a/highs/simplex/HApp.h b/highs/simplex/HApp.h
@@ -155,7 +155,7 @@ 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)
+  if (kSimplexScaleDevReport)
     incumbent_lp.scale_.print("grepSimplexScaling",
                               "After scaling, " + incumbent_lp.model_name_ +
 			      "," + incumbent_lp.origin_name_);
diff --git a/highs/simplex/HEkk.cpp b/highs/simplex/HEkk.cpp
@@ -3613,10 +3613,11 @@ void HEkk::testBasisCondition(const HighsLp& lp,
   exact = false;
   double approx_kappa = this->computeBasisCondition(lp, exact);
   approx_kappa_tt += timer_->read();
-  highsLogUser(options_->log_options, HighsLogType::kInfo,
-               "getKappa,%s,%g,%g,%g,%g,%s,%s\n", lp.model_name_.c_str(),
-               exact_kappa, exact_kappa_tt, approx_kappa, approx_kappa_tt,
-               message.c_str(), lp.origin_name_.c_str());
+  if (kSimplexScaleDevReport) 
+    highsLogUser(options_->log_options, HighsLogType::kInfo,
+		 "getKappa,%s,%g,%g,%g,%g,%s,%s\n", lp.model_name_.c_str(),
+		 exact_kappa, exact_kappa_tt, approx_kappa, approx_kappa_tt,
+		 message.c_str(), lp.origin_name_.c_str());
 }
 
 double HEkk::computeBasisCondition(const HighsLp& lp, const bool exact) const {
diff --git a/highs/simplex/HSimplex.cpp b/highs/simplex/HSimplex.cpp
@@ -484,16 +484,18 @@ void simplexScaleLp(const HighsOptions& options, HighsLp& lp,
 	max_cost = std::max(abs_cost, max_cost);
       }
     }
-    printf("grepSimplexRangeTxt HighsScale: costs in [%g, %g]; matrix in [%g, %g]; %s: %s\n",
-	   min_cost, max_cost,
-	   original_matrix_min_value, original_matrix_max_value,
-	   lp.model_name_.c_str(),
-	   lp.origin_name_.c_str());
-    printf("grepSimplexRangeCsv,%g,%g,%g,%g,%s,%s\n",
-	   min_cost, max_cost,
-	   original_matrix_min_value, original_matrix_max_value,
-	   lp.model_name_.c_str(),
-	   lp.origin_name_.c_str());
+    if (kSimplexScaleDevReport) {
+      printf("grepSimplexRangeTxt HighsScale: costs in [%g, %g]; matrix in [%g, %g]; %s: %s\n",
+	     min_cost, max_cost,
+	     original_matrix_min_value, original_matrix_max_value,
+	     lp.model_name_.c_str(),
+	     lp.origin_name_.c_str());
+      printf("grepSimplexRangeCsv,%g,%g,%g,%g,%s,%s\n",
+	     min_cost, max_cost,
+	     original_matrix_min_value, original_matrix_max_value,
+	     lp.model_name_.c_str(),
+	     lp.origin_name_.c_str());
+    }
   }
   // Possibly force scaling, otherwise base the decision on the range
   // of values in the matrix, values that will be used later for
@@ -504,21 +506,20 @@ void simplexScaleLp(const HighsOptions& options, HighsLp& lp,
                                  no_scaling_original_matrix_min_value) &&
                                     (original_matrix_max_value <=
                                      no_scaling_original_matrix_max_value);
+  HighsScale& scale = lp.scale_;
   bool scaled_matrix = false;
   if (no_scaling) {
-    // No matrix scaling, but possible cost scaling
-    if (options.highs_analysis_level)
-      highsLogDev(options.log_options, HighsLogType::kInfo,
-                  "Scaling: Matrix has [min, max] values of [%g, %g] within "
-                  "[%g, %g] so no scaling performed\n",
-                  original_matrix_min_value, original_matrix_max_value,
-                  no_scaling_original_matrix_min_value,
-                  no_scaling_original_matrix_max_value);
+    // No matrix scaling
+    highsLogDev(options.log_options, HighsLogType::kInfo,
+		"Scaling: Matrix has [min, max] values of [%g, %g] within "
+		"[%g, %g] so no scaling performed\n",
+		original_matrix_min_value, original_matrix_max_value,
+		no_scaling_original_matrix_min_value,
+		no_scaling_original_matrix_max_value);
   } else {
     // Try scaling, so assign unit factors - partly because initial
     // factors may be assumed by the scaling method, but also because
     // scaling factors may not be computed for empty rows/columns
-    HighsScale& scale = lp.scale_;
     scale.col.assign(numCol, 1);
     scale.row.assign(numRow, 1);
     const bool equilibration_scaling =
@@ -549,17 +550,18 @@ void simplexScaleLp(const HighsOptions& options, HighsLp& lp,
       scale.cost = 1.0;
       lp.is_scaled_ = true;
     }
-    bool scaled_costs = false;
-    if (kSimplexScaleDev &&
-	use_scale_strategy == kSimplexScaleStrategyMaxValue0157) {
-      // Consider cost scaling
-      simplexScaleCost(options, lp);
-      scaled_costs = scale.cost != 1.0;
-      lp.is_scaled_ = true;
-    }
-    // If neither costs nor matrix are scaled, then clear the scaling
-    if (!scaled_costs && !scaled_matrix) lp.clearScaling();
   }
+  bool scaled_costs = false;
+  if (kSimplexScaleDev &&
+      use_scale_strategy == kSimplexScaleStrategyMaxValueMatrixAndCost) {
+    // Consider cost scaling
+    simplexScaleCost(options, lp);
+    scaled_costs = scale.cost != 1.0;
+    scale.has_scaling = true;
+    lp.is_scaled_ = true;
+  }
+  // If neither costs nor matrix are scaled, then clear the scaling
+  if (!scaled_costs && !scaled_matrix) lp.clearScaling();
   // Record the scaling strategy used
   lp.scale_.strategy = use_scale_strategy;
 }
@@ -898,26 +900,7 @@ bool maxValueScaleMatrix(const HighsOptions& options, HighsLp& lp,
   vector<HighsInt>& Aindex = lp.a_matrix_.index_;
   vector<double>& Avalue = lp.a_matrix_.value_;
 
-  assert(options.simplex_scale_strategy == kSimplexScaleStrategyMaxValue);
-  assert(kSimplexScaleStrategyMaxValue015 == kSimplexScaleStrategyMaxValue);
-  assert(kSimplexScaleStrategyMaxValue0157 == kSimplexScaleStrategyMaxValue);
-
-  // The 015(7) values refer to bit settings in FICO's scaling options.
-  // Specifically
-  //
-  // 0: Row scaling
-  //
-  // 1: Column scaling
-  //
-  // 5: Scale by maximum element
-  //
-  // 7: Scale objective function for the simplex method
-  //
-  // Note that 7 is not yet implemented, so
-  // kSimplexScaleStrategyMaxValue015 and
-  // kSimplexScaleStrategyMaxValue0157 are equivalent. However, cost
-  // scaling could be well worth adding, now that the unscaled problem
-  // can be solved using scaled NLA
+  assert(options.simplex_scale_strategy >= kSimplexScaleStrategyMaxValue);
 
   const double log2 = log(2.0);
   const double max_allow_scale = pow(2.0, options.allowed_matrix_scale_factor);
@@ -1093,40 +1076,38 @@ void simplexScaleCost(const HighsOptions& options, HighsLp& lp) {
   // Scale the costs by no less than minAlwCostScale
   double max_allowed_cost_scale = pow(2.0, options.allowed_cost_scale_factor);
   double max_nonzero_cost = 0;
-  for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) {
-    if (lp.col_cost_[iCol]) {
+  for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++)
+    if (lp.col_cost_[iCol]) 
       max_nonzero_cost = max(fabs(lp.col_cost_[iCol]), max_nonzero_cost);
-    }
-  }
-  // Scaling the costs up effectively increases the dual tolerance to
-  // which the problem is solved - so, if the max cost is small the
-  // scaling factor pushes it up by a power of 2 so it's close to 1
-  // Scaling the costs down effectively decreases the dual tolerance
-  // to which the problem is solved - so this can't be done too much
-  cost_scale = 1;
-  const double ln2 = log(2.0);
-  // Scale if the max cost is positive and outside the range [1/16, 16]
-  if ((max_nonzero_cost > 0) &&
-      ((max_nonzero_cost < (1.0 / 16)) || (max_nonzero_cost > 16))) {
+  // Scaling the costs up effectively means that the problem is solved
+  // to a smaller dual tolerance, which seems pointless. HiGHS will
+  // warn about excessively small costs, and users can set
+  // user_objective_scale to a positive value.
+  //
+  // Scaling the costs down effectively means that the problem is
+  // solved to a larger dual tolerance, which may require clean-up
+  // iterations when the scaling is removed after solving the scaled
+  // problem
+  cost_scale = 1.0;
+  // Scale if the max cost is greater than 16
+  if (max_nonzero_cost > 16) {
     cost_scale = max_nonzero_cost;
-    cost_scale = pow(2.0, floor(log(cost_scale) / ln2 + 0.5));
+    cost_scale = pow(2.0, floor(log(cost_scale) / log(2.0) + 0.5));
     cost_scale = min(cost_scale, max_allowed_cost_scale);
   }
-  if (cost_scale == 1) {
-    highsLogUser(options.log_options, HighsLogType::kInfo,
-                 "LP cost vector not scaled down: max cost is %g\n",
-                 max_nonzero_cost);
-    return;
-  }
-  // Scale the costs (and record of max_nonzero_cost) by cost_scale, being at
-  // most max_allowed_cost_scale
-  for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) {
-    lp.col_cost_[iCol] /= cost_scale;
+  if (cost_scale == 1.0) {
+    highsLogDev(options.log_options, HighsLogType::kInfo,
+		"LP cost vector not scaled down: max cost is %g\n",
+		max_nonzero_cost);
+  } else {
+    // Scale the costs (and record of max_nonzero_cost) by cost_scale, being at
+    // most max_allowed_cost_scale
+    for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++)
+      lp.col_cost_[iCol] /= cost_scale;
+    highsLogDev(options.log_options, HighsLogType::kInfo,
+		"LP cost vector scaled down by %g: max scaled cost is %g\n", cost_scale,
+		max_nonzero_cost/cost_scale);
   }
-  max_nonzero_cost /= cost_scale;
-  highsLogUser(options.log_options, HighsLogType::kInfo,
-               "LP cost vector scaled down by %g: max cost is %g\n", cost_scale,
-               max_nonzero_cost);
 }
 
 void simplexUnscaleCost(HighsLp& lp) {
