Merge pull request ERGO-Code#2268 from ERGO-Code/fix-2267

Fix 2267

Author: jajhall

FEATURES.md

@@ -2,7 +2,8 @@
 
 ## Code changes
 
-Fixed incorrect assertion in `HighsMipSolver::solutionFeasible()` (fixing #2204)
+Fixed incorrect assertion in `HighsMipSolver::solutionFeasible()` (fixing [#2204](https://github.com/ERGO-Code/HiGHS/issues/2204)))
 
-getColIntegrality now returns `HighsVarType::kContinuous` when `model_.lp_.integrality_` is empty (fixing #2261)
+getColIntegrality now returns `HighsVarType::kContinuous` when `model_.lp_.integrality_` is empty (fixing [#2261](https://github.com/ERGO-Code/HiGHS/issues/2261))
 
+Now ensuring that when solving a scaled LP with useful but unvalidated basis, it does not lose its scaling after validation, since the scaling factors will be applied to the solution (fixing [#2267](https://github.com/ERGO-Code/HiGHS/issues/2267))

check/TestCAPI.c

@@ -2091,6 +2091,38 @@ iteration_count1); assertLogical("Dual", logic);
   Highs_destroy(highs);
 }
 */
+void testDeleteRowResolveWithBasis() {
+  // Created to expose bug in #2267, but also illustrates case where
+  // scaling was performed on the original problem - due to the 6 in
+  // row 1, but would not be performed on the problem after row 1 has
+  // been removed
+  void* highs = Highs_create();
+  Highs_setBoolOptionValue(highs, "output_flag", dev_run);
+  HighsInt ret;
+  double INF = Highs_getInfinity(highs);
+  ret = Highs_addCol(highs, 0.0, 2.0, 2.0, 0, NULL, NULL);
+  ret = Highs_addCol(highs, 0.0, -INF, INF, 0, NULL, NULL);
+  ret = Highs_addCol(highs, 0.0, -INF, INF, 0, NULL, NULL);
+  HighsInt index_1[2] = {0, 2};
+  double value_1[2] = {2.0, -1.0};
+  ret = Highs_addRow(highs, 0.0, 0.0, 2, index_1, value_1);
+  HighsInt index_2[1] = {1};
+  double value_2[1] = {6.0};
+  ret = Highs_addRow(highs, 10.0, INF, 1, index_2, value_2);
+  Highs_run(highs);
+  double col_value[3] = {0.0, 0.0, 0.0};
+  Highs_getSolution(highs, col_value, NULL, NULL, NULL);
+  assertDoubleValuesEqual("col_value[0]", col_value[0], 2.0);
+  ret = Highs_deleteRowsByRange(highs, 1, 1);
+  assert(ret == 0);
+  ret = Highs_run(highs);
+  assert(ret == 0);
+  ret = Highs_getSolution(highs, col_value, NULL, NULL, NULL);
+  assert(ret == 0);
+  assertDoubleValuesEqual("col_value[0]", col_value[0], 2.0);
+  Highs_destroy(highs);
+}
+
 int main() {
   minimalApiIllegalLp();
   testCallback();
@@ -2113,6 +2145,8 @@ int main() {
   testMultiObjective();
   testQpIndefiniteFailure();
   testDualRayTwice();
+
+  testDeleteRowResolveWithBasis();
   return 0;
 }
 //  testSetSolution();

src/interfaces/highs_c_api.cpp

@@ -744,11 +744,11 @@ HighsInt Highs_setCallback(void* highs, HighsCCallbackType user_callback,
   return static_cast<int>(status);
 }
 
-HighsInt Highs_startCallback(void* highs, const int callback_type) {
+HighsInt Highs_startCallback(void* highs, const HighsInt callback_type) {
   return (HighsInt)((Highs*)highs)->startCallback(callback_type);
 }
 
-HighsInt Highs_stopCallback(void* highs, const int callback_type) {
+HighsInt Highs_stopCallback(void* highs, const HighsInt callback_type) {
   return (HighsInt)((Highs*)highs)->stopCallback(callback_type);
 }
 
@@ -1291,9 +1291,9 @@ HighsInt Highs_getPresolvedLp(const void* highs, const HighsInt a_format,
                               a_start, a_index, a_value, integrality);
 }
 
-HighsInt Highs_crossover(void* highs, const int num_col, const int num_row,
-                         const double* col_value, const double* col_dual,
-                         const double* row_dual) {
+HighsInt Highs_crossover(void* highs, const HighsInt num_col,
+                         const HighsInt num_row, const double* col_value,
+                         const double* col_dual, const double* row_dual) {
   HighsSolution solution;
   if (col_value) {
     solution.value_valid = true;

src/interfaces/highs_c_api.h

@@ -770,7 +770,7 @@ HighsInt Highs_getStringOptionValue(const void* highs, const char* option,
  *
  * @param highs     A pointer to the Highs instance.
  * @param option    The name of the option.
- * @param type      An int in which the corresponding `kHighsOptionType`
+ * @param type      A HighsInt in which the corresponding `kHighsOptionType`
  *                  constant should be placed.
  *
  * @returns A `kHighsStatus` constant indicating whether the call succeeded.
@@ -842,7 +842,7 @@ HighsInt Highs_getBoolOptionValues(const void* highs, const char* option,
                                    HighsInt* current_value,
                                    HighsInt* default_value);
 /**
- * Get the current and default values of an int option
+ * Get the current and default values of a HighsInt option
  *
  * @param highs         A pointer to the Highs instance.
  * @param current_value A pointer to the current value of the option.
@@ -924,7 +924,7 @@ HighsInt Highs_getInt64InfoValue(const void* highs, const char* info,
  *
  * @param highs     A pointer to the Highs instance.
  * @param info      The name of the info item.
- * @param type      An int in which the corresponding `kHighsOptionType`
+ * @param type      A HighsInt in which the corresponding `kHighsOptionType`
  *                  constant is stored.
  *
  * @returns A `kHighsStatus` constant indicating whether the call succeeded.
@@ -983,7 +983,7 @@ HighsInt Highs_getModelStatus(const void* highs);
  * LP) gets it if it does not and dual_ray_value is not nullptr.
  *
  * @param highs             A pointer to the Highs instance.
- * @param has_dual_ray      A pointer to an int to store 1 if a dual ray
+ * @param has_dual_ray      A pointer to a HighsInt to store 1 if a dual ray
  *                          currently exists.
  * @param dual_ray_value    An array of length [num_row] filled with the
  *                          unbounded ray.
@@ -1001,9 +1001,10 @@ HighsInt Highs_getDualRay(const void* highs, HighsInt* has_dual_ray,
  *
  * @param highs                                   A pointer to the Highs
  *                                                instance.
- * @param has_dual_unboundedness_direction        A pointer to an int to store 1
- *                                                if the dual unboundedness
- *                                                direction exists.
+ * @param has_dual_unboundedness_direction        A pointer to a HighsInt to
+ *                                                store 1 if the dual
+ *                                                unboundedness direction
+ *                                                exists.
  * @param dual_unboundedness_direction_value      An array of length [num_col]
  *                                                filled with the unboundedness
  *                                                direction.
@@ -1018,7 +1019,7 @@ HighsInt Highs_getDualUnboundednessDirection(
  * LP) gets it if it does not and primal_ray_value is not nullptr.
  *
  * @param highs             A pointer to the Highs instance.
- * @param has_primal_ray    A pointer to an int to store 1 if the primal ray
+ * @param has_primal_ray    A pointer to a HighsInt to store 1 if the primal ray
  *                          exists.
  * @param primal_ray_value  An array of length [num_col] filled with the
  *                          unbounded ray.
@@ -1279,7 +1280,7 @@ HighsInt Highs_setCallback(void* highs, HighsCCallbackType user_callback,
  *
  * @returns A `kHighsStatus` constant indicating whether the call succeeded.
  */
-HighsInt Highs_startCallback(void* highs, const int callback_type);
+HighsInt Highs_startCallback(void* highs, const HighsInt callback_type);
 
 /**
  * Stop callback of given type
@@ -1289,7 +1290,7 @@ HighsInt Highs_startCallback(void* highs, const int callback_type);
  *
  * @returns A `kHighsStatus` constant indicating whether the call succeeded.
  */
-HighsInt Highs_stopCallback(void* highs, const int callback_type);
+HighsInt Highs_stopCallback(void* highs, const HighsInt callback_type);
 
 /**
  * Return the cumulative wall-clock time spent in `Highs_run`.
@@ -1754,15 +1755,15 @@ HighsInt Highs_getObjectiveOffset(const void* highs, double* offset);
  * @param highs         A pointer to the Highs instance.
  * @param from_col      The first column for which to query data for.
  * @param to_col        The last column (inclusive) for which to query data for.
- * @param num_col       An integer populated with the number of columns got from
+ * @param num_col       A HighsInt populated with the number of columns got from
  *                      the model (this should equal `to_col - from_col + 1`).
  * @param costs         An array of size [to_col - from_col + 1] for the column
  *                      cost coefficients.
  * @param lower         An array of size [to_col - from_col + 1] for the column
  *                      lower bounds.
  * @param upper         An array of size [to_col - from_col + 1] for the column
  *                      upper bounds.
- * @param num_nz        An integer to be populated with the number of non-zero
+ * @param num_nz        A HighsInt to be populated with the number of non-zero
  *                      elements in the constraint matrix.
  * @param matrix_start  An array of size [to_col - from_col + 1] with the start
  *                      indices of each column in `matrix_index` and
@@ -1831,13 +1832,13 @@ HighsInt Highs_getColsByMask(const void* highs, const HighsInt* mask,
  * @param highs         A pointer to the Highs instance.
  * @param from_row      The first row for which to query data for.
  * @param to_row        The last row (inclusive) for which to query data for.
- * @param num_row       An integer to be populated with the number of rows got
+ * @param num_row       A HighsInt to be populated with the number of rows got
  *                      from the model.
  * @param lower         An array of size [to_row - from_row + 1] for the row
  *                      lower bounds.
  * @param upper         An array of size [to_row - from_row + 1] for the row
  *                      upper bounds.
- * @param num_nz        An integer to be populated with the number of non-zero
+ * @param num_nz        A HighsInt to be populated with the number of non-zero
  *                      elements in the constraint matrix.
  * @param matrix_start  An array of size [to_row - from_row + 1] with the start
  *                      indices of each row in `matrix_index` and
@@ -1940,7 +1941,7 @@ HighsInt Highs_getColByName(const void* highs, const char* name, HighsInt* col);
  * Get the integrality of a column.
  *
  * @param col          The index of the column to query.
- * @param integrality  An integer in which the integrality of the column should
+ * @param integrality  A HighsInt in which the integrality of the column should
  *                     be placed. The integer is one of the `kHighsVarTypeXXX`
  *                     constants.
  *
@@ -1994,7 +1995,7 @@ HighsInt Highs_deleteColsByMask(void* highs, HighsInt* mask);
  *
  * @returns A `kHighsStatus` constant indicating whether the call succeeded.
  */
-HighsInt Highs_deleteRowsByRange(void* highs, const int from_row,
+HighsInt Highs_deleteRowsByRange(void* highs, const HighsInt from_row,
                                  const HighsInt to_row);
 
 /**
@@ -2222,9 +2223,9 @@ HighsInt Highs_getPresolvedLp(const void* highs, const HighsInt a_format,
  *
  * @returns A `kHighsStatus` constant indicating whether the call succeeded.
  */
-HighsInt Highs_crossover(void* highs, const int num_col, const int num_row,
-                         const double* col_value, const double* col_dual,
-                         const double* row_dual);
+HighsInt Highs_crossover(void* highs, const HighsInt num_col,
+                         const HighsInt num_row, const double* col_value,
+                         const double* col_dual, const double* row_dual);
 
 /**
  * Compute the ranging information for all costs and bounds. For

src/lp_data/HighsLpUtils.cpp

@@ -1439,6 +1439,7 @@ HighsStatus applyScalingToLpRow(HighsLp& lp, const HighsInt row,
 }
 
 void unscaleSolution(HighsSolution& solution, const HighsScale& scale) {
+  assert(scale.has_scaling);
   assert(solution.col_value.size() == static_cast<size_t>(scale.num_col));
   assert(solution.col_dual.size() == static_cast<size_t>(scale.num_col));
   assert(solution.row_value.size() == static_cast<size_t>(scale.num_row));

src/lp_data/HighsSolution.cpp

@@ -1357,8 +1357,9 @@ HighsStatus formSimplexLpBasisAndFactor(HighsLpSolverObject& solver_object,
   HEkk& ekk_instance = solver_object.ekk_instance_;
   HighsSimplexStatus& ekk_status = ekk_instance.status_;
   lp.ensureColwise();
+  const bool passed_scaled = lp.is_scaled_;
   // Consider scaling the LP
-  const bool new_scaling = considerScaling(options, lp);
+  if (!passed_scaled) considerScaling(options, lp);
   const bool check_basis = basis.alien || (!basis.valid && basis.useful);
   if (check_basis) {
     // The basis needs to be checked for rank deficiency, and possibly
@@ -1370,7 +1371,11 @@ HighsStatus formSimplexLpBasisAndFactor(HighsLpSolverObject& solver_object,
     assert(!only_from_known_basis);
     accommodateAlienBasis(solver_object);
     basis.alien = false;
-    lp.unapplyScale();
+    // Unapply any scaling used only for factorization to check and
+    // complete the basis
+    if (!passed_scaled) lp.unapplyScale();
+    // Check that any scaling the LP arrived with has not been removed
+    assert(lp.is_scaled_ == passed_scaled);
     return HighsStatus::kOk;
   }
   // Move the HighsLpSolverObject's LP to EKK

src/simplex/HApp.h

@@ -134,6 +134,7 @@ inline HighsStatus solveLpSimplex(HighsLpSolverObject& solver_object) {
   // then move to EKK
   considerScaling(options, incumbent_lp);
   //
+  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
     // that is not validated
@@ -149,6 +150,8 @@ inline HighsStatus solveLpSimplex(HighsLpSolverObject& solver_object) {
     refineBasis(incumbent_lp, solution, basis);
     basis.valid = true;
   }
+  // Check that scaling has not been removed
+  assert(incumbent_lp.is_scaled_ == was_scaled);
   // Move the LP to EKK, updating other EKK pointers and any simplex
   // NLA pointers, since they may have moved if the LP has been
   // modified
@@ -221,6 +224,9 @@ inline HighsStatus solveLpSimplex(HighsLpSolverObject& solver_object) {
             kSimplexUnscaledSolutionStrategyNone ||
         options.simplex_unscaled_solution_strategy ==
             kSimplexUnscaledSolutionStrategyRefine) {
+      // Check that the incumbent LP has scaling and is scaled
+      assert(incumbent_lp.scale_.has_scaling);
+      assert(incumbent_lp.is_scaled_);
       //
       // Solve the scaled LP!
       //