Merge pull request #2691 from ERGO-Code/fix-2653-local

Fix spurious primal infeasibilities

Author: jajhall

highs/lp_data/HighsSolution.cpp

@@ -549,28 +549,24 @@ void getVariableKktFailures(const double primal_feasibility_tolerance,
                             const HighsVarType integrality,
                             double& primal_infeasibility,
                             double& dual_infeasibility, uint8_t& at_status,
-                            uint8_t& mid_status) {
-  const HighsInt feasibility_tolerance_mu = 0.0;
-  // @primal_infeasibility calculation
-  primal_infeasibility = 0;
-  if (value < lower - primal_feasibility_tolerance * feasibility_tolerance_mu) {
-    // Below lower
-    primal_infeasibility = lower - value;
-  } else if (value >
-             upper + primal_feasibility_tolerance * feasibility_tolerance_mu) {
-    // Above upper
-    primal_infeasibility = value - upper;
-  }
+                            uint8_t& mid_status, const HighsInt index) {
+  // Return the primal residual (ie infeasibility with zero tolerance)
+  // as the primal infeasibility, ensuring (cf #2653) that it doesn't
+  // exceed the primal feasibility tolerance if the standard primal
+  // infeasibility (ie infeasibility exceeding the tolerance) is zero
+  auto infeasibility_residual =
+      infeasibility(lower, value, upper, primal_feasibility_tolerance);
+  primal_infeasibility = infeasibility_residual.second;
   // Determine whether this value is close to a bound
   at_status = kHighsSolutionNo;
-  double residual = std::fabs(lower - value);
-  if (residual * residual <= primal_feasibility_tolerance) {
+  double bound_residual = std::fabs(lower - value);
+  if (bound_residual * bound_residual <= primal_feasibility_tolerance) {
     // Close to lower bound
     at_status = kHighsSolutionLo;
   } else {
     // Not close to lower bound: maybe close to upper bound
-    residual = std::fabs(value - upper);
-    if (residual * residual <= primal_feasibility_tolerance)
+    bound_residual = std::fabs(value - upper);
+    if (bound_residual * bound_residual <= primal_feasibility_tolerance)
       at_status = kHighsSolutionUp;
   }
   // Look for dual sign errors that exceed the tolerance. For boxed

highs/lp_data/HighsSolution.h

@@ -90,7 +90,7 @@ void getVariableKktFailures(const double primal_feasibility_tolerance,
                             const HighsVarType integrality,
                             double& primal_infeasibility,
                             double& dual_infeasibility, uint8_t& at_status,
-                            uint8_t& mid_status);
+                            uint8_t& mid_status, const HighsInt index = 0);
 
 void getPrimalDualGlpsolErrors(const HighsOptions& options, const HighsLp& lp,
                                const std::vector<double>& gradient,

highs/util/HighsUtils.h

@@ -214,4 +214,58 @@ inline T fractionality(T input, T* intval = nullptr) {
   if (intval != nullptr) *intval = val;
   return abs(input - val);
 }
+
+inline std::pair<double, double> infeasibility(const double lower,
+                                               const double value,
+                                               const double upper,
+                                               const double tolerance) {
+  using std::fabs;
+  using std::min;
+  double residual = 0;
+  double infeasibility = 0;
+  // Determine the infeasibility exceeding the tolerance used in
+  // computing the number of infeasibilities in a basic solution -
+  // which defines its feasibility
+  //
+  // @primal_infeasibility calculation
+  if (value < lower - tolerance) infeasibility = lower - value;
+  if (value > upper + tolerance) infeasibility = value - upper;
+  // Determine the residual used in computing the sum of
+  // infeasibilities and max infeasibility - which are just for
+  // reporting
+  if (tolerance > 0) {
+    if (value < lower) residual = lower - value;
+    if (value > upper) residual = value - upper;
+  } else {
+    residual = infeasibility;
+  }
+  // Now, if the bound defining the residual is large, it's possible
+  // for the infeasibility to be zero, but the residual to exceed the
+  // tolerance due to numerical rounding
+  //
+  // Case in point is #2653 where row 1 has l = 157345 and the
+  // activity gives zero infeasibility, but a residual of
+  // 1.00000761449e-06, exceeding the tolerance of 1e-6 by delta =
+  // 7.61449e-12.
+  //
+  // Now delta / l = 4.83937e-17, which is less than machine precision
+  // (2.22045e-16), so the reliable residual value should ignore this
+  // delta.
+  //
+  // In general, it might be possible to subtract delta from the
+  // residual conditional on (something like)
+  //
+  // delta < 1e1 * max(1.0, fabs(bound_value)) * kHighsMacheps
+  //
+  // to give a reasonable value but, in practice, when infeasibility is
+  // 0, it would seem fine to set
+  //
+  // residual = min(residual, tolerance)
+  //
+  // so that values of maximum infeasibility defined by residual
+  // doesn't exceed the tolerance
+  //
+  if (infeasibility == 0) residual = min(residual, tolerance);
+  return std::make_pair(infeasibility, residual);
+}
 #endif  // UTIL_HIGHSUTILS_H_