Merge branch 'latest' of https://github.com/ERGO-Code/HiGHS into dualBoundTightening

Author: fwesselm

FEATURES.md

@@ -12,3 +12,4 @@ The irreducible infeasibility system (IIS) facility now detects infeasibility du
 
 Prompted by [#2463](https://github.com/ERGO-Code/HiGHS/issues/2463), the HiGHS solution and basis files now match data to any column and row names in the model, only assuming that the data are aligned with column and row indices if there are no names in the model. This requires a new version (v2) of the HiGHS basis file. Basis files from v1 are still read, but deprecated. Now, when writing out a model, basis or solution, column and row names are added to the model - previously they were created temporarily and inconsistentyly on the fly. If the model has existing names, then distinctive names are created to replace any blank names, but names with spaces or duplicate names yield an error status return.
 
+As per [#2487](https://github.com/ERGO-Code/HiGHS/issues/2487), trivial heuristics now run before feasibility jump (FJ), and FJ will use any existing incumbent. FJ will clip any finite variable values in the incumbent to lower and upper bounds, and falls back to the existing logic (lower bound if finite, else upper bound if finite, else 0) for any infinite values in the incumbent.

highs/ipm/ipx/ipm.cc

@@ -821,13 +821,15 @@ void IPM::PrintHeader() {
   std::stringstream h_logging_stream;
   h_logging_stream.str(std::string());
   h_logging_stream
-    << (kTerminationLogging ? "\n" : "")
-    << " "  << Format("Iter", 4)
-    << "  " << Format("P.res", 8) << " " << Format("D.res", 8)
-    << "  " << Format("P.obj", 15) << " " << Format("D.obj", 15)
-    << "  " << Format("mu", 8);
+    << " " << Format("Iter", 4)
+    << "  " << Format("primal obj", 15)
+    << "  " << Format("dual obj", 15)
+    << "  " << Format("pinf", 9)
+    << "  " << Format("dinf", 9)
+    << "  " << Format("gap", 8);
+    //  h_logging_stream << "  " << Format("mu", 8);
   if (!control_.timelessLog())
-    h_logging_stream << "  " << Format("Time", 7);
+    h_logging_stream << "  " << Format("time", 7);
   control_.hLog(h_logging_stream);
   control_.Debug()
     << "  " << Format("stepsizes", 9)
@@ -841,17 +843,29 @@ void IPM::PrintHeader() {
 void IPM::PrintOutput() {
     const bool ipm_optimal = iterate_->feasible() && iterate_->optimal();
 
-    if (kTerminationLogging) PrintHeader();
+    double logging_pobj = iterate_->pobjective_after_postproc();
+    double logging_dobj = iterate_->dobjective_after_postproc();
+    double logging_presidual = iterate_->presidual();
+    double logging_dresidual = iterate_->dresidual();
+
+    // Now logging relative primal and dual infeasibility, and also
+    // the relative primal dual objective gap
+    logging_presidual /= iterate_->bounds_measure_;
+    logging_dresidual /= iterate_->costs_measure_;
+    double logging_gap = std::abs(logging_pobj - logging_dobj) /
+      (1.0+0.5 *std::fabs(logging_pobj + logging_dobj));
+
     std::stringstream h_logging_stream;
     h_logging_stream.str(std::string());
     h_logging_stream
       << " "  << Format(info_->iter, 3)
       << (ipm_optimal ? "*" : " ")
-      << "  " << Scientific(iterate_->presidual(), 8, 2)
-      << " "  << Scientific(iterate_->dresidual(), 8, 2)
-      << "  " << Scientific(iterate_->pobjective_after_postproc(), 15, 8)
-      << " "  << Scientific(iterate_->dobjective_after_postproc(), 15, 8)
-      << "  " << Scientific(iterate_->mu(), 8, 2);
+      << "  " << Scientific(logging_pobj, 15, 8)
+      << "  " << Scientific(logging_dobj, 15, 8)
+      << "  " << Scientific(logging_presidual, 9, 2)
+      << "  " << Scientific(logging_dresidual, 9, 2)
+      << "  " << Scientific(logging_gap, 8, 2);
+    //    h_logging_stream << "  " << Scientific(iterate_->mu(), 8, 2);
     if (!control_.timelessLog())
       h_logging_stream << "  " << Fixed(control_.Elapsed(), 6, 0) << "s";
     control_.hLog(h_logging_stream);

highs/ipm/ipx/iterate.cc

@@ -55,6 +55,9 @@ Iterate::Iterate(const Model& model) : model_(model) {
         }
     }
     assert_consistency();
+    this->bounds_measure_ = 1.0 + model_.norm_bounds();
+    this->costs_measure_ = 1.0 + model_.norm_c();
+    
 }
 
 void Iterate::Initialize(const Vector& x, const Vector& xl, const Vector& xu,
@@ -219,40 +222,21 @@ double Iterate::mu_max() const { Evaluate(); return mu_max_; }
 
 bool Iterate::feasible() const {
     Evaluate();
-    const double bounds_measure = 1.0 + model_.norm_bounds();
-    const double costs_measure = 1.0 + model_.norm_c();
-    const double rel_presidual = presidual_ / bounds_measure;
-    const double rel_dresidual = dresidual_ / costs_measure;
-    const bool primal_feasible = presidual_ <= feasibility_tol_ * (bounds_measure);
-    const bool dual_feasible = dresidual_ <= feasibility_tol_ * (costs_measure);
+    const bool primal_feasible = presidual_ <= feasibility_tol_ * bounds_measure_;
+    const bool dual_feasible = dresidual_ <= feasibility_tol_ * costs_measure_;
     const bool is_feasible = primal_feasible && dual_feasible;
-    if (kTerminationLogging) {
-      printf("\nIterate::feasible presidual_ = %11.4g; bounds_measure = %11.4g; "
-	     "rel_presidual = %11.4g; feasibility_tol = %11.4g: primal_feasible = %d\n",
-	     presidual_, bounds_measure, rel_presidual, feasibility_tol_, primal_feasible);
-      printf("Iterate::feasible dresidual_ = %11.4g;  costs_measure = %11.4g; "
-	     "rel_dresidual = %11.4g; feasibility_tol = %11.4g:   dual_feasible = %d\n",
-	     dresidual_, costs_measure, rel_dresidual, feasibility_tol_, dual_feasible);
-    }
     return is_feasible;
 }
 
 bool Iterate::optimal() const {
     Evaluate();
     double pobj = pobjective_after_postproc();
     double dobj = dobjective_after_postproc();
-    double obj = 0.5 * (pobj + dobj);
+    double ave_obj = 0.5 * (pobj + dobj);
     double gap = pobj - dobj;
     const double abs_gap = std::abs(gap);
-    const double obj_measure = 1.0+std::abs(obj);
+    const double obj_measure = 1.0+std::abs(ave_obj);
     const bool is_optimal = abs_gap <= optimality_tol_ * obj_measure;
-    if (kTerminationLogging) {
-      const double rel_gap = abs_gap / obj_measure; 
-      printf("Iterate::optimal     abs_gap = %11.4g;"
-	     " obj_measure    = %11.4g; rel_gap       = %11.4g;"
-	     "  optimality_tol = %11.4g:   optimal       = %d\n",
-	     abs_gap, obj_measure, rel_gap, optimality_tol_, is_optimal);
-    }
     return is_optimal;
 }
 

highs/ipm/ipx/iterate.h

@@ -199,6 +199,9 @@ class Iterate {
     // The method can only be called after Postprocess().
     void DropToComplementarity(Vector& x, Vector& y, Vector& z) const;
 
+    double bounds_measure_;
+    double costs_measure_;
+
 private:
     // A (primal or dual) variable that is required to be positive in the IPM is
     // not moved closer to zero than kBarrierMin.

highs/ipm/ipx/utils.h

@@ -4,8 +4,6 @@
 #include <vector>
 #include "ipm/ipx/ipx_internal.h"
 
-const bool kTerminationLogging = false;
-
 namespace ipx {
 
 bool AllFinite(const Vector& x);

highs/mip/HighsFeasibilityJump.cpp

@@ -27,6 +27,8 @@ HighsModelStatus HighsMipSolverData::feasibilityJump() {
   std::vector<double> col_value(model->num_col_, 0.0);
   double objective_function_value;
 
+  const bool use_incumbent = !incumbent.empty();
+
   // Configure Feasibility Jump and pass it the problem
   auto solver = external_feasibilityjump::FeasibilityJumpSolver(
       log_options,
@@ -66,10 +68,14 @@ HighsModelStatus HighsMipSolverData::feasibilityJump() {
                   sense_multiplier * model->col_cost_[col]);
 
     double initial_assignment = 0.0;
-    if (std::isfinite(lower)) {
-      initial_assignment = lower;
-    } else if (std::isfinite(upper)) {
-      initial_assignment = upper;
+    if (use_incumbent && std::isfinite(incumbent[col])) {
+      initial_assignment = std::max(lower, std::min(upper, incumbent[col]));
+    } else {
+      if (std::isfinite(lower)) {
+        initial_assignment = lower;
+      } else if (std::isfinite(upper)) {
+        initial_assignment = upper;
+      }
     }
     col_value[col] = initial_assignment;
   }

highs/mip/HighsMipSolver.cpp

@@ -148,8 +148,19 @@ void HighsMipSolver::run() {
       mipdata_->queryExternalSolution(
           solution_objective_, kExternalMipSolutionQueryOriginAfterSetup);
 
+    // Apply the trivial heuristics
+    analysis_.mipTimerStart(kMipClockTrivialHeuristics);
+    HighsModelStatus returned_model_status = mipdata_->trivialHeuristics();
+    analysis_.mipTimerStop(kMipClockTrivialHeuristics);
+    if (modelstatus_ == HighsModelStatus::kNotset &&
+        returned_model_status == HighsModelStatus::kInfeasible) {
+      // trivialHeuristics can spot trivial infeasibility, so act on it
+      modelstatus_ = returned_model_status;
+      cleanupSolve();
+      return;
+    }
+    // Apply the feasibility jump heuristic (if enabled)
     if (options_mip_->mip_heuristic_run_feasibility_jump) {
-      // Apply the feasibility jump before evaluating the root node
       analysis_.mipTimerStart(kMipClockFeasibilityJump);
       HighsModelStatus returned_model_status = mipdata_->feasibilityJump();
       analysis_.mipTimerStop(kMipClockFeasibilityJump);
@@ -161,17 +172,7 @@ void HighsMipSolver::run() {
         return;
       }
     }
-    // Apply the trivial heuristics
-    analysis_.mipTimerStart(kMipClockTrivialHeuristics);
-    HighsModelStatus returned_model_status = mipdata_->trivialHeuristics();
-    analysis_.mipTimerStop(kMipClockTrivialHeuristics);
-    if (modelstatus_ == HighsModelStatus::kNotset &&
-        returned_model_status == HighsModelStatus::kInfeasible) {
-      // trivialHeuristics can spot trivial infeasibility, so act on it
-      modelstatus_ = returned_model_status;
-      cleanupSolve();
-      return;
-    }
+    // End of pre-root-node heuristics
     if (analysis_.analyse_mip_time && !submip)
       if (analysis_.analyse_mip_time & !submip)
         highsLogUser(options_mip_->log_options, HighsLogType::kInfo,

highs/pdlp/cupdlp/cupdlp_solver.c

@@ -843,15 +843,14 @@ void PDHG_Print_Iter(CUPDLPwork *pdhg) {
   else
     cupdlp_snprintf(timeString, 8, "%6ds", (cupdlp_int)timers->dSolvingTime);
 
-  // cupdlp_printf("%9d  %+15.8e  %+15.8e  %+8.2e  %8.2e  %10.2e  %8.2e %7s
-  // [L]\n",
-  //               timers->nIter, resobj->dPrimalObj, resobj->dDualObj,
-  //               resobj->dDualityGap, resobj->dComplementarity,
-  //               resobj->dPrimalFeas, resobj->dDualFeas, timeString);
-
+  CUPDLPscaling *scaling = pdhg->scaling;
+  const double gap_log_value = resobj->dRelObjGap;
+  const double primal_feas_log_value = resobj->dPrimalFeas / (1.0 + scaling->dNormRhs);
+  const double dual_feas_log_value = resobj->dDualFeas / (1.0 + scaling->dNormCost);
   cupdlp_printf("%9d  %+15.8e  %+15.8e  %+8.2e  %10.2e  %8.2e %7s [L]\n",
                 timers->nIter, resobj->dPrimalObj, resobj->dDualObj,
-                resobj->dDualityGap, resobj->dPrimalFeas, resobj->dDualFeas,
+		//                resobj->dDualityGap, resobj->dPrimalFeas, resobj->dDualFeas,
+		gap_log_value, primal_feas_log_value, dual_feas_log_value,
                 timeString);
 }
 
@@ -865,16 +864,15 @@ void PDHG_Print_Iter_Average(CUPDLPwork *pdhg) {
   else
     cupdlp_snprintf(timeString, 8, "%6ds", (cupdlp_int)timers->dSolvingTime);
 
-  // cupdlp_printf("%9d  %+15.8e  %+15.8e  %+8.2e  %8.2e  %10.2e  %8.2e %7s
-  // [A]\n",
-  //               timers->nIter, resobj->dPrimalObjAverage,
-  //               resobj->dDualObjAverage, resobj->dDualityGapAverage,
-  //               resobj->dComplementarityAverage, resobj->dPrimalFeasAverage,
-  //               resobj->dDualFeasAverage, timeString);
+  CUPDLPscaling *scaling = pdhg->scaling;
+  const double gap_log_value = resobj->dRelObjGapAverage;
+  const double primal_feas_log_value = resobj->dPrimalFeasAverage / (1.0 + scaling->dNormRhs);
+  const double dual_feas_log_value = resobj->dDualFeasAverage / (1.0 + scaling->dNormCost);
   cupdlp_printf("%9d  %+15.8e  %+15.8e  %+8.2e  %10.2e  %8.2e %7s [A]\n",
                 timers->nIter, resobj->dPrimalObjAverage,
-                resobj->dDualObjAverage, resobj->dDualityGapAverage,
-                resobj->dPrimalFeasAverage, resobj->dDualFeasAverage,
+                resobj->dDualObjAverage,
+		//		resobj->dDualityGapAverage, resobj->dPrimalFeasAverage, resobj->dDualFeasAverage,
+		gap_log_value, primal_feas_log_value, dual_feas_log_value,
                 timeString);
 }