Merge pull request #1918 from fwesselm/minorStuffAgain

Add function to compute fractional part

src/interfaces/highs_c_api.h

@@ -18,7 +18,7 @@
 // Highs_mipCall or Highs_qpCall, and these methods return the
 // appropriate solution information
 //
-// For sophisticated applications, where esoteric solutiuon
+// For sophisticated applications, where esoteric solution
 // information is needed, or if a sequence of modified models need to
 // be solved, use the Highs_create method to generate a pointer to an
 // instance of the C++ Highs class, and then use any of a large number

src/lp_data/Highs.cpp

@@ -3924,10 +3924,9 @@ HighsStatus Highs::callRunPostsolve(const HighsSolution& solution,
       max_integrality_violation = 0;
       for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) {
         if (lp.integrality_[iCol] == HighsVarType::kInteger) {
-          const double value = this->solution_.col_value[iCol];
-          double intval = std::floor(value + 0.5);
           max_integrality_violation =
-              std::max(fabs(intval - value), max_integrality_violation);
+              std::max(fractionality(this->solution_.col_value[iCol]),
+                       max_integrality_violation);
         }
       }
       highsLogUser(

src/lp_data/HighsInterface.cpp

@@ -1712,7 +1712,7 @@ HighsStatus Highs::elasticityFilterReturn(
 
   if (return_status == HighsStatus::kOk) {
     // Solution is invalidated by deleting rows and columns, but
-    // primal values are correct. Have to recompute row acivities,
+    // primal values are correct. Have to recompute row activities,
     // though
     this->model_.lp_.a_matrix_.productQuad(this->solution_.row_value,
                                            this->solution_.col_value);
@@ -1842,15 +1842,15 @@ HighsStatus Highs::elasticityFilter(
       const double lower_penalty = has_local_lower_penalty
                                        ? local_lower_penalty[iCol]
                                        : global_lower_penalty;
-      // Negative lower penalty and infininte upper bound implies that the
+      // Negative lower penalty and infinite upper bound implies that the
       // bounds cannot be violated
       if (lower_penalty < 0 && upper >= kHighsInf) continue;
 
       // Get the penalty for violating the upper bounds on this column
       const double upper_penalty = has_local_upper_penalty
                                        ? local_upper_penalty[iCol]
                                        : global_upper_penalty;
-      // Infininte upper bound and negative lower penalty implies that the
+      // Infinite upper bound and negative lower penalty implies that the
       // bounds cannot be violated
       if (lower <= -kHighsInf && upper_penalty < 0) continue;
       erow_lower.push_back(lower);

src/lp_data/HighsLpUtils.cpp

@@ -2305,8 +2305,7 @@ void assessColPrimalSolution(const HighsOptions& options, const double primal,
   }
   integer_infeasibility = 0;
   if (type == HighsVarType::kInteger || type == HighsVarType::kSemiInteger) {
-    double nearest_integer = std::floor(primal + 0.5);
-    integer_infeasibility = std::fabs(primal - nearest_integer);
+    integer_infeasibility = fractionality(primal);
   }
   if (col_infeasibility > 0 && (type == HighsVarType::kSemiContinuous ||
                                 type == HighsVarType::kSemiInteger)) {

src/lp_data/HighsSolve.cpp

@@ -136,7 +136,7 @@ HighsStatus solveLp(HighsLpSolverObject& solver_object, const string message) {
       // and duality gap that are within the tolerances supplied by
       // HiGHS, the HiGHS primal and dual feasibility tolerances may
       // not be satisfied since they are absolute, and in PDLP they
-      // are relative. Note that, even when only one PDLP row activit
+      // are relative. Note that, even when only one PDLP row activity
       // fails to satisfy the absolute tolerance, the absolute norm
       // measure reported by PDLP will not necessarily be the same as
       // with HiGHS, since PDLP uses the 2-norm, and HiGHS the

src/mip/HighsLpRelaxation.cpp

@@ -1224,9 +1224,8 @@ HighsLpRelaxation::Status HighsLpRelaxation::resolveLp(HighsDomain* domain) {
           double val = std::max(
               std::min(sol.col_value[i], lpsolver.getLp().col_upper_[i]),
               lpsolver.getLp().col_lower_[i]);
-          double intval = std::floor(val + 0.5);
 
-          if (std::abs(val - intval) > mipsolver.mipdata_->feastol) {
+          if (fractionality(val) > mipsolver.mipdata_->feastol) {
             HighsInt col = i;
             if (roundable && mipsolver.mipdata_->uplocks[col] != 0 &&
                 mipsolver.mipdata_->downlocks[col] != 0)

src/mip/HighsMipSolver.cpp

@@ -63,9 +63,8 @@ HighsMipSolver::HighsMipSolver(HighsCallback& callback,
       obj += orig_model_->col_cost_[i] * value;
 
       if (orig_model_->integrality_[i] == HighsVarType::kInteger) {
-        double intval = std::floor(value + 0.5);
         integrality_violation_ =
-            std::max(fabs(intval - value), integrality_violation_);
+            std::max(fractionality(value), integrality_violation_);
       }
 
       const double lower = orig_model_->col_lower_[i];

src/mip/HighsMipSolverData.cpp

@@ -27,7 +27,7 @@ bool HighsMipSolverData::checkSolution(
     if (solution[i] < mipsolver.model_->col_lower_[i] - feastol) return false;
     if (solution[i] > mipsolver.model_->col_upper_[i] + feastol) return false;
     if (mipsolver.variableType(i) == HighsVarType::kInteger &&
-        std::abs(solution[i] - std::floor(solution[i] + 0.5)) > feastol)
+        fractionality(solution[i]) > feastol)
       return false;
   }
 
@@ -57,7 +57,7 @@ bool HighsMipSolverData::trySolution(const std::vector<double>& solution,
     if (solution[i] < mipsolver.model_->col_lower_[i] - feastol) return false;
     if (solution[i] > mipsolver.model_->col_upper_[i] + feastol) return false;
     if (mipsolver.variableType(i) == HighsVarType::kInteger &&
-        std::abs(solution[i] - std::floor(solution[i] + 0.5)) > feastol)
+        fractionality(solution[i]) > feastol)
       return false;
 
     obj += mipsolver.colCost(i) * solution[i];
@@ -523,14 +523,10 @@ void HighsMipSolverData::runSetup() {
     HighsInt end = ARstart_[i + 1];
     bool integral = true;
     for (HighsInt j = start; j != end; ++j) {
-      if (integral) {
-        if (mipsolver.variableType(ARindex_[j]) == HighsVarType::kContinuous)
-          integral = false;
-        else {
-          double intval = std::floor(ARvalue_[j] + 0.5);
-          if (std::abs(ARvalue_[j] - intval) > epsilon) integral = false;
-        }
-      }
+      integral =
+          integral &&
+          mipsolver.variableType(ARindex_[j]) != HighsVarType::kContinuous &&
+          fractionality(ARvalue_[j]) <= epsilon;
 
       maxabsval = std::max(maxabsval, std::abs(ARvalue_[j]));
     }
@@ -597,8 +593,7 @@ void HighsMipSolverData::runSetup() {
         break;
       case HighsVarType::kInteger:
         if (domain.isFixed(i)) {
-          if (std::abs(domain.col_lower_[i] -
-                       std::floor(domain.col_lower_[i] + 0.5)) > feastol) {
+          if (fractionality(domain.col_lower_[i]) > feastol) {
             // integer variable is fixed to a fractional value -> infeasible
             mipsolver.modelstatus_ = HighsModelStatus::kInfeasible;
             lower_bound = kHighsInf;
@@ -729,8 +724,7 @@ double HighsMipSolverData::transformNewIntegerFeasibleSolution(
         mipsolver.orig_model_->col_cost_[i] * value;
 
     if (mipsolver.orig_model_->integrality_[i] == HighsVarType::kInteger) {
-      double intval = std::floor(value + 0.5);
-      double integrality_infeasibility = std::fabs(intval - value);
+      double integrality_infeasibility = fractionality(value);
       if (integrality_infeasibility >
           mipsolver.options_mip_->mip_feasibility_tolerance) {
         if (debug_report)

src/mip/HighsSearch.cpp

@@ -48,8 +48,7 @@ double HighsSearch::checkSol(const std::vector<double>& sol,
     if (!integerfeasible || mipsolver.variableType(i) != HighsVarType::kInteger)
       continue;
 
-    double intval = std::floor(sol[i] + 0.5);
-    if (std::abs(sol[i] - intval) > mipsolver.mipdata_->feastol) {
+    if (fractionality(sol[i]) > mipsolver.mipdata_->feastol) {
       integerfeasible = false;
     }
   }

src/mip/HighsTableauSeparator.cpp

@@ -65,25 +65,23 @@ void HighsTableauSeparator::separateLpSolution(HighsLpRelaxation& lpRelaxation,
   std::vector<FractionalInteger> fractionalBasisvars;
   fractionalBasisvars.reserve(numRow);
   for (HighsInt i = 0; i < numRow; ++i) {
-    double fractionality;
+    double my_fractionality;
     if (basisinds[i] >= numCol) {
       HighsInt row = basisinds[i] - numCol;
 
       if (!lpRelaxation.isRowIntegral(row)) continue;
 
-      double solval = lpSolution.row_value[row];
-      fractionality = std::fabs(std::round(solval) - solval);
+      my_fractionality = fractionality(lpSolution.row_value[row]);
     } else {
       HighsInt col = basisinds[i];
       if (mip.variableType(col) == HighsVarType::kContinuous) continue;
 
-      double solval = lpSolution.col_value[col];
-      fractionality = std::fabs(std::round(solval) - solval);
+      my_fractionality = fractionality(lpSolution.col_value[col]);
     }
 
-    if (fractionality < 1000 * mip.mipdata_->feastol) continue;
+    if (my_fractionality < 1000 * mip.mipdata_->feastol) continue;
 
-    fractionalBasisvars.emplace_back(i, fractionality);
+    fractionalBasisvars.emplace_back(i, my_fractionality);
   }
 
   if (fractionalBasisvars.empty()) return;

src/presolve/HPresolve.cpp

@@ -173,8 +173,7 @@ bool HPresolve::okSetInput(HighsMipSolver& mipsolver,
 
 bool HPresolve::rowCoefficientsIntegral(HighsInt row, double scale) const {
   for (const HighsSliceNonzero& nz : getRowVector(row)) {
-    double val = nz.value() * scale;
-    if (std::abs(val - std::round(val)) > options->small_matrix_value)
+    if (fractionality(nz.value() * scale) > options->small_matrix_value)
       return false;
   }
 
@@ -252,9 +251,8 @@ bool HPresolve::isImpliedIntegral(HighsInt col) {
       double scale = 1.0 / nz.value();
       if (!rowCoefficientsIntegral(nz.index(), scale)) continue;
 
-      double rhs = model->row_lower_[nz.index()] * scale;
-
-      if (std::abs(rhs - std::round(rhs)) > primal_feastol) {
+      if (fractionality(model->row_lower_[nz.index()] * scale) >
+          primal_feastol) {
         // todo infeasible
       }
 
@@ -323,9 +321,9 @@ bool HPresolve::isImpliedInteger(HighsInt col) {
       // if there is an equation the dual detection does not need to be tried
       runDualDetection = false;
       double scale = 1.0 / nz.value();
-      double rhs = model->row_lower_[nz.index()] * scale;
 
-      if (std::abs(rhs - std::round(rhs)) > primal_feastol) {
+      if (fractionality(model->row_lower_[nz.index()] * scale) >
+          primal_feastol) {
         continue;
       }
 
@@ -338,24 +336,20 @@ bool HPresolve::isImpliedInteger(HighsInt col) {
   if (!runDualDetection) return false;
 
   if ((model->col_lower_[col] != -kHighsInf &&
-       std::abs(std::round(model->col_lower_[col]) - model->col_lower_[col]) >
-           options->small_matrix_value) ||
+       fractionality(model->col_lower_[col]) > options->small_matrix_value) ||
       (model->col_upper_[col] != -kHighsInf &&
-       std::abs(std::round(model->col_upper_[col]) - model->col_upper_[col]) >
-           options->small_matrix_value))
+       fractionality(model->col_upper_[col]) > options->small_matrix_value))
     return false;
 
   for (const HighsSliceNonzero& nz : getColumnVector(col)) {
     double scale = 1.0 / nz.value();
-    if (model->row_upper_[nz.index()] != kHighsInf) {
-      double rhs = model->row_upper_[nz.index()];
-      if (std::abs(rhs - std::round(rhs)) > primal_feastol) return false;
-    }
+    if (model->row_upper_[nz.index()] != kHighsInf &&
+        fractionality(model->row_upper_[nz.index()]) > primal_feastol)
+      return false;
 
-    if (model->row_lower_[nz.index()] != -kHighsInf) {
-      double rhs = model->row_lower_[nz.index()];
-      if (std::abs(rhs - std::round(rhs)) > primal_feastol) return false;
-    }
+    if (model->row_lower_[nz.index()] != -kHighsInf &&
+        fractionality(model->row_lower_[nz.index()]) > primal_feastol)
+      return false;
 
     if (!rowCoefficientsIntegral(nz.index(), scale)) return false;
   }
@@ -3305,7 +3299,7 @@ HPresolve::Result HPresolve::rowPresolve(HighsPostsolveStack& postsolve_stack,
 
           if (intScale != 0.0 && intScale <= 1e3) {
             double rhs = rowUpper * intScale;
-            if (std::abs(rhs - std::round(rhs)) > primal_feastol)
+            if (fractionality(rhs) > primal_feastol)
               return Result::kPrimalInfeasible;
 
             rhs = std::round(rhs);
@@ -3341,8 +3335,7 @@ HPresolve::Result HPresolve::rowPresolve(HighsPostsolveStack& postsolve_stack,
               HighsInt x1Pos = rowpositions[x1Cand];
               HighsInt x1 = Acol[x1Pos];
               double rhs2 = rhs / static_cast<double>(d);
-              if (std::abs(std::round(rhs2) - rhs2) <=
-                  mipsolver->mipdata_->epsilon) {
+              if (fractionality(rhs2) <= mipsolver->mipdata_->epsilon) {
                 // the right hand side is integral, so we can substitute
                 // x1 = d * z
 
@@ -5723,10 +5716,9 @@ HPresolve::Result HPresolve::detectParallelRowsAndCols(
         }
 
         double scaleCand = colMax[duplicateCol].first / colMax[col].first;
-        colScale = std::round(scaleCand);
-        assert(std::abs(colScale) >= 1.0);
-        if (std::abs(colScale - scaleCand) > options->small_matrix_value)
+        if (fractionality(scaleCand, &colScale) > options->small_matrix_value)
           continue;
+        assert(std::abs(colScale) >= 1.0);
 
         // if the scale is larger than 1, duplicate column cannot compensate for
         // all values of scaled col due to integrality as the scaled column

src/presolve/HighsPostsolveStack.cpp

@@ -15,6 +15,7 @@
 #include "lp_data/HConst.h"
 #include "lp_data/HighsOptions.h"
 #include "util/HighsCDouble.h"
+#include "util/HighsUtils.h"
 
 namespace presolve {
 
@@ -728,8 +729,7 @@ void HighsPostsolveStack::DuplicateColumn::undo(const HighsOptions& options,
   } else if (duplicateColIntegral) {
     // Doesn't set basis.col_status[duplicateCol], so assume no basis
     assert(!basis.valid);
-    double roundVal = std::round(solution.col_value[duplicateCol]);
-    if (std::abs(roundVal - solution.col_value[duplicateCol]) >
+    if (fractionality(solution.col_value[duplicateCol]) >
         options.mip_feasibility_tolerance) {
       solution.col_value[duplicateCol] =
           std::floor(solution.col_value[duplicateCol]);
@@ -927,9 +927,7 @@ bool HighsPostsolveStack::DuplicateColumn::okMerge(
   if (x_int) {
     if (y_int) {
       // Scale must be integer and not exceed (x_u-x_l)+1 in magnitude
-      double int_scale = std::floor(scale + 0.5);
-      bool scale_is_int = std::fabs(int_scale - scale) <= tolerance;
-      if (!scale_is_int) {
+      if (fractionality(scale) > tolerance) {
         if (debug_report)
           printf(
               "%sDuplicateColumn::checkMerge: scale must be integer, but is "
@@ -1012,8 +1010,7 @@ void HighsPostsolveStack::DuplicateColumn::undoFix(
   //=============================================================================================
 
   auto isInteger = [&](const double v) {
-    double int_v = std::floor(v + 0.5);
-    return std::fabs(int_v - v) <= mip_feasibility_tolerance;
+    return (fractionality(v) <= mip_feasibility_tolerance);
   };
 
   auto isFeasible = [&](const double l, const double v, const double u) {

src/util/HVectorBase.h

@@ -77,7 +77,7 @@ class HVectorBase {
   vector<Real> packValue;      //!< Packed values
 
   /**
-   * @brief Copy from another HVector structure to this instanc
+   * @brief Copy from another HVector structure to this instance
    */
   template <typename FromReal>
   void copy(const HVectorBase<FromReal>*

src/util/HighsUtils.h

@@ -207,4 +207,14 @@ void highsAssert(const bool assert_condition, const std::string message = "");
 // If pause_condition is true, then keyboard input is required. Allows
 // breakpoints in VScode where optimization might prevent them.
 bool highsPause(const bool pause_condition, const std::string message = "");
+
+// Utility for computing fractional part
+template <typename T>
+inline T fractionality(T input, T* intval = nullptr) {
+  using std::abs;
+  using std::round;
+  T val = round(input);
+  if (intval != nullptr) *intval = val;
+  return abs(input - val);
+}
 #endif  // UTIL_HIGHSUTILS_H_