Fix edge-case assertion failures in presolve and domain propagation code

highs/mip/HighsDebugSol.cpp

@@ -68,10 +68,12 @@ void HighsDebugSol::activate() {
 
       HighsCDouble debugsolobj = 0.0;
       for (HighsInt i = 0; i != mipsolver->orig_model_->num_col_; ++i)
-        debugsolobj += mipsolver->orig_model_->col_cost_[i] *
-                       HighsCDouble(debugOrigSolution[i]);
+        debugsolobj +=
+            static_cast<HighsCDouble>(mipsolver->orig_model_->col_cost_[i]) *
+            debugOrigSolution[i];
 
-      debugSolObjective = double(debugsolobj + mipsolver->orig_model_->offset_);
+      debugSolObjective =
+          static_cast<double>(debugsolobj + mipsolver->orig_model_->offset_);
       debugSolActive = true;
       printf("debug sol active\n");
       registerDomain(mipsolver->mipdata_->getDomain());
@@ -101,9 +103,10 @@ void HighsDebugSol::shrink(const std::vector<HighsInt>& newColIndex) {
   debugSolution.resize(mipsolver->model_->num_col_);
   HighsCDouble debugsolobj = 0.0;
   for (HighsInt i = 0; i != mipsolver->model_->num_col_; ++i)
-    debugsolobj += mipsolver->model_->col_cost_[i] * debugSolution[i];
+    debugsolobj += static_cast<HighsCDouble>(mipsolver->model_->col_cost_[i]) *
+                   debugSolution[i];
 
-  debugSolObjective = double(debugsolobj);
+  debugSolObjective = static_cast<double>(debugsolobj);
 
   conflictingBounds.clear();
 }

highs/mip/HighsDomain.h

@@ -537,7 +537,7 @@ class HighsDomain {
 
   double getObjectiveLowerBound() const {
     if (objProp_.isActive() && objProp_.numInfObjLower == 0)
-      return double(objProp_.objectiveLower);
+      return static_cast<double>(objProp_.objectiveLower);
 
     return -kHighsInf;
   }
@@ -561,7 +561,7 @@ class HighsDomain {
     std::vector<HighsDomainChange> reducedstack;
     reducedstack.reserve(domchgstack_.size());
     branchingPositions.reserve(branchPos_.size());
-    for (HighsInt i = 0; i < (HighsInt)domchgstack_.size(); ++i) {
+    for (HighsInt i = 0; i < static_cast<HighsInt>(domchgstack_.size()); ++i) {
       // keep only the tightest bound change for each variable
       if ((domchgstack_[i].boundtype == HighsBoundType::kLower &&
            colLowerPos_[domchgstack_[i].column] != i) ||
@@ -620,11 +620,13 @@ class HighsDomain {
                            double& rhs) const;
 
   double getMinActivity(HighsInt row) const {
-    return activitymininf_[row] == 0 ? double(activitymin_[row]) : -kHighsInf;
+    return activitymininf_[row] == 0 ? static_cast<double>(activitymin_[row])
+                                     : -kHighsInf;
   }
 
   double getMaxActivity(HighsInt row) const {
-    return activitymaxinf_[row] == 0 ? double(activitymax_[row]) : kHighsInf;
+    return activitymaxinf_[row] == 0 ? static_cast<double>(activitymax_[row])
+                                     : kHighsInf;
   }
 
   double getMinCutActivity(const HighsCutPool& cutpool, HighsInt cut) const;
@@ -654,6 +656,8 @@ class HighsDomain {
 
   double feastol() const;
 
+  double epsilon() const;
+
   HighsInt numModelNonzeros() const { return mipsolver->numNonzero(); }
 
   bool inSubmip() const { return mipsolver->submip; }

highs/mip/HighsLpRelaxation.cpp

@@ -446,7 +446,7 @@ double HighsLpRelaxation::computeBestEstimate(const HighsPseudocost& ps) const {
     estimate += double(increase);
   }
 
-  return double(estimate);
+  return static_cast<double>(estimate);
 }
 
 double HighsLpRelaxation::computeLPDegneracy(
@@ -1612,7 +1612,8 @@ HighsLpRelaxation::Status HighsLpRelaxation::resolveLp(HighsDomain* domain) {
           }
 
           for (HighsInt i = 0; i != mipsolver.numCol(); ++i)
-            objsum += roundsol[i] * mipsolver.colCost(i);
+            objsum +=
+                static_cast<HighsCDouble>(roundsol[i]) * mipsolver.colCost(i);
 
           if (!mipsolver.mipdata_->parallelLockActive() || !worker_) {
             mipsolver.mipdata_->addIncumbent(
@@ -1625,7 +1626,8 @@ HighsLpRelaxation::Status HighsLpRelaxation::resolveLp(HighsDomain* domain) {
         }
 
         for (HighsInt i = 0; i != mipsolver.numCol(); ++i)
-          objsum += sol.col_value[i] * mipsolver.colCost(i);
+          objsum += static_cast<HighsCDouble>(sol.col_value[i]) *
+                    mipsolver.colCost(i);
 
         if (fractionalints.empty() && !unscaledPrimalFeasible(status)) {
           std::vector<double> fixSol = sol.col_value;
@@ -1647,7 +1649,7 @@ HighsLpRelaxation::Status HighsLpRelaxation::resolveLp(HighsDomain* domain) {
           }
         }
 
-        objective = double(objsum);
+        objective = static_cast<double>(objsum);
         break;
       }
       case Status::kInfeasible:

highs/mip/HighsMipSolver.cpp

@@ -1365,7 +1365,7 @@ bool HighsMipSolver::solutionFeasible(const HighsLp* lp,
     assert(col_value.size() == static_cast<size_t>(lp->num_col_));
   for (HighsInt i = 0; i != lp->num_col_; ++i) {
     const double value = col_value[i];
-    obj += lp->col_cost_[i] * value;
+    obj += static_cast<HighsCDouble>(lp->col_cost_[i]) * value;
 
     if (lp->integrality_[i] == HighsVarType::kInteger) {
       integrality_violation =

highs/mip/HighsMipSolverData.cpp

@@ -226,7 +226,7 @@ bool HighsMipSolverData::trySolution(const std::vector<double>& solution,
     if (mipsolver.isColInteger(i) && fractionality(solution[i]) > feastol)
       return false;
 
-    obj += mipsolver.colCost(i) * solution[i];
+    obj += static_cast<HighsCDouble>(mipsolver.colCost(i)) * solution[i];
   }
 
   for (HighsInt i = 0; i != mipsolver.numRow(); ++i) {
@@ -389,8 +389,9 @@ HighsModelStatus HighsMipSolverData::trivialHeuristics() {
 
     HighsCDouble cdouble_obj = 0.0;
     for (HighsInt iCol = 0; iCol < mipsolver.numCol(); iCol++)
-      cdouble_obj += mipsolver.colCost(iCol) * solution[iCol];
-    double obj = double(cdouble_obj);
+      cdouble_obj +=
+          static_cast<HighsCDouble>(mipsolver.colCost(iCol)) * solution[iCol];
+    double obj = static_cast<double>(cdouble_obj);
     const double save_upper_bound = upper_bound;
     const bool new_incumbent =
         addIncumbent(solution, obj, heuristic_source[try_heuristic]);
@@ -1145,8 +1146,8 @@ void HighsMipSolverData::runSetup() {
     debugSolution.debugSolObjective = 0;
     HighsCDouble debugsolobj = 0.0;
     for (HighsInt i = 0; i != mipsolver.numCol(); ++i)
-      debugsolobj +=
-          mipsolver.colCost(i) * HighsCDouble(debugSolution.debugSolution[i]);
+      debugsolobj += static_cast<HighsCDouble>(mipsolver.colCost(i)) *
+                     debugSolution.debugSolution[i];
     debugSolution.debugSolObjective = static_cast<double>(debugsolobj);
     debugSolution.registerDomain(getDomain());
     assert(checkSolution(debugSolution.debugSolution));

highs/mip/HighsMipWorker.cpp

@@ -127,7 +127,7 @@ bool HighsMipWorker::trySolution(const std::vector<double>& solution,
         fractionality(solution[i]) > mipdata_.feastol)
       return false;
 
-    obj += mipsolver_.colCost(i) * solution[i];
+    obj += static_cast<HighsCDouble>(mipsolver_.colCost(i)) * solution[i];
   }
 
   for (HighsInt i = 0; i != mipsolver_.model_->num_row_; ++i) {

highs/mip/HighsSearch.cpp

@@ -44,7 +44,7 @@ double HighsSearch::checkSol(const std::vector<double>& sol,
   HighsCDouble objval = 0.0;
   integerfeasible = true;
   for (HighsInt i = 0; i != mipsolver.numCol(); ++i) {
-    objval += sol[i] * mipsolver.colCost(i);
+    objval += static_cast<HighsCDouble>(sol[i]) * mipsolver.colCost(i);
     assert(std::isfinite(sol[i]));
 
     if (!integerfeasible || !mipsolver.isColInteger(i)) continue;

highs/presolve/HPresolve.cpp

@@ -2952,7 +2952,8 @@ void HPresolve::substitute(HighsInt row, HighsInt col, double rhs) {
 
   // substitute column in the objective function
   if (model->col_cost_[col] != 0.0) {
-    HighsCDouble objscale = model->col_cost_[col] * substrowscale;
+    HighsCDouble objscale =
+        static_cast<HighsCDouble>(model->col_cost_[col]) * substrowscale;
     model->offset_ = static_cast<double>(model->offset_ - objscale * rhs);
     assert(std::isfinite(model->offset_));
     for (HighsInt rowiter : rowpositions) {
@@ -2967,7 +2968,7 @@ void HPresolve::substitute(HighsInt row, HighsInt col, double rhs) {
     }
     assert(std::abs(model->col_cost_[col]) <=
            std::max(options->dual_feasibility_tolerance,
-                    kHighsTiny * std::abs(static_cast<double>(objscale))));
+                    static_cast<double>(kHighsTiny * abs(objscale))));
     model->col_cost_[col] = 0.0;
   }
 
@@ -7104,10 +7105,12 @@ HPresolve::Result HPresolve::strengthenInequalities(
       double ub = model->col_upper_[col] - model->col_lower_[col];
       if (weight > 0) {
         comp = 1;
-        maxviolation += model->col_upper_[col] * weight;
+        maxviolation +=
+            static_cast<HighsCDouble>(model->col_upper_[col]) * weight;
       } else {
         comp = -1;
-        maxviolation += model->col_lower_[col] * weight;
+        maxviolation +=
+            static_cast<HighsCDouble>(model->col_lower_[col]) * weight;
         weight = -weight;
       }
 
@@ -7227,16 +7230,15 @@ HPresolve::Result HPresolve::strengthenInequalities(
                               HighsInt direction) {
       for (HighsInt i : indices) {
         assert(Arow[positions[i]] == row);
-        double coefdelta =
-            direction * static_cast<double>(reducedcost[i] - maxviolation);
+        HighsCDouble coefdelta = direction * (reducedcost[i] - maxviolation);
         HighsInt col = Acol[positions[i]];
 
         if (complementation[i] == -1) {
           rhs += coefdelta * model->col_lower_[col];
-          addToMatrix(row, col, coefdelta);
+          addToMatrix(row, col, static_cast<double>(coefdelta));
         } else {
           rhs -= coefdelta * model->col_upper_[col];
-          addToMatrix(row, col, -coefdelta);
+          addToMatrix(row, col, static_cast<double>(-coefdelta));
         }
       }
     };