Merge pull request #2014 from ERGO-Code/fix-1946

Compute primal-dual integral

src/highs_bindings.cpp

@@ -832,7 +832,9 @@ PYBIND11_MODULE(_core, m) {
       .def_readwrite("max_complementarity_violation",
                      &HighsInfo::max_complementarity_violation)
       .def_readwrite("sum_complementarity_violations",
-                     &HighsInfo::sum_complementarity_violations);
+                     &HighsInfo::sum_complementarity_violations)
+      .def_readwrite("primal_dual_integral",
+                     &HighsInfo::primal_dual_integral);
   py::class_<HighsOptions>(m, "HighsOptions")
       .def(py::init<>())
       .def_readwrite("presolve", &HighsOptions::presolve)

src/interfaces/highs_c_api.cpp

@@ -529,6 +529,17 @@ HighsInt Highs_getDualRay(const void* highs, HighsInt* has_dual_ray,
   return retcode;
 }
 
+HighsInt Highs_getDualUnboundednessDirection(
+    const void* highs, HighsInt* has_dual_unboundedness_direction,
+    double* dual_unboundedness_direction_value) {
+  bool v;
+  HighsInt retcode = (HighsInt)((Highs*)highs)
+                         ->getDualUnboundednessDirection(
+                             v, dual_unboundedness_direction_value);
+  *has_dual_unboundedness_direction = (HighsInt)v;
+  return retcode;
+}
+
 HighsInt Highs_getPrimalRay(const void* highs, HighsInt* has_primal_ray,
                             double* primal_ray_value) {
   bool v;

src/interfaces/highs_c_api.h

@@ -915,11 +915,13 @@ HighsInt Highs_getBasis(const void* highs, HighsInt* col_status,
 HighsInt Highs_getModelStatus(const void* highs);
 
 /**
- * Get an unbounded dual ray that is a certificate of primal infeasibility.
+ * Indicates whether a dual ray that is a certificate of primal
+ * infeasibility currently exists, and (at the expense of solving an
+ * 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 the dual ray
- *                          exists.
+ * @param has_dual_ray      A pointer to an int to store 1 if a dual ray
+ *                          currently exists.
  * @param dual_ray_value    An array of length [num_row] filled with the
  *                          unbounded ray.
  *
@@ -929,7 +931,28 @@ HighsInt Highs_getDualRay(const void* highs, HighsInt* has_dual_ray,
                           double* dual_ray_value);
 
 /**
- * Get an unbounded primal ray that is a certificate of dual infeasibility.
+ * Indicates whether a dual unboundedness direction (corresponding to a
+ * certificate of primal infeasibility) exists, and (at the expense of
+ * solving an LP) gets it if it does not and
+ * dual_unboundedness_direction is not nullptr
+ *
+ * @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 dual_unboundedness_direction_value      An array of length [num_col]
+ *                                                filled with the unboundedness
+ *                                                direction.
+ */
+HighsInt getDualUnboundednessDirection(
+    const void* highs, HighsInt* has_dual_unboundedness_direction,
+    double* dual_unboundedness_direction_value);
+
+/**
+ * Indicates whether a primal ray that is a certificate of primal
+ * unboundedness currently exists, and (at the expense of solving an
+ * 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

src/lp_data/Highs.cpp

@@ -3938,6 +3938,7 @@ HighsStatus Highs::callSolveMip() {
   info_.mip_node_count = solver.node_count_;
   info_.mip_dual_bound = solver.dual_bound_;
   info_.mip_gap = solver.gap_;
+  info_.primal_dual_integral = solver.primal_dual_integral_;
   // Get the number of LP iterations, avoiding overflow if the int64_t
   // value is too large
   int64_t mip_total_lp_iterations = solver.total_lp_iterations_;
@@ -4502,6 +4503,7 @@ HighsStatus Highs::returnFromHighs(HighsStatus highs_return_status) {
 }
 
 void Highs::reportSolvedLpQpStats() {
+  if (!options_.output_flag) return;
   HighsLogOptions& log_options = options_.log_options;
   if (this->model_.lp_.model_name_.length())
     highsLogUser(log_options, HighsLogType::kInfo, "Model name          : %s\n",
@@ -4541,11 +4543,7 @@ void Highs::reportSolvedLpQpStats() {
         std::fabs(info_.objective_function_value - dual_objective_value) /
         std::max(1.0, std::fabs(info_.objective_function_value));
     highsLogUser(log_options, HighsLogType::kInfo,
-                 "Highs::reportSolvedLpQpStats Objective for %s: primal = "
-                 "%17.10e; dual = %17.10e; rel gap = %17.10e\n",
-                 this->model_.lp_.model_name_.c_str(),
-                 info_.objective_function_value, dual_objective_value,
-                 relative_primal_dual_gap);
+                 "Relative P-D gap    : %17.10e\n", relative_primal_dual_gap);
   }
   double run_time = timer_.readRunHighsClock();
   highsLogUser(log_options, HighsLogType::kInfo,

src/lp_data/HighsInfo.cpp

@@ -40,6 +40,7 @@ void HighsInfo::invalidate() {
   sum_dual_infeasibilities = kHighsIllegalInfeasibilityMeasure;
   max_complementarity_violation = kHighsIllegalComplementarityViolation;
   sum_complementarity_violations = kHighsIllegalComplementarityViolation;
+  primal_dual_integral = -kHighsInf;
 }
 
 static std::string infoEntryTypeToString(const HighsInfoType type) {

src/lp_data/HighsInfo.h

@@ -115,6 +115,7 @@ struct HighsInfoStruct {
   double sum_dual_infeasibilities;
   double max_complementarity_violation;
   double sum_complementarity_violations;
+  double primal_dual_integral;
 };
 
 class HighsInfo : public HighsInfoStruct {
@@ -273,6 +274,11 @@ class HighsInfo : public HighsInfoStruct {
         "sum_complementarity_violations", "Sum of complementarity violations",
         advanced, &sum_complementarity_violations, 0);
     records.push_back(record_double);
+
+    record_double =
+        new InfoRecordDouble("primal_dual_integral", "Primal-dual integral",
+                             advanced, &primal_dual_integral, 0);
+    records.push_back(record_double);
   }
 
  public:

src/mip/HighsMipSolver.cpp

@@ -225,8 +225,16 @@ void HighsMipSolver::run() {
   search.setLpRelaxation(&mipdata_->lp);
   sepa.setLpRelaxation(&mipdata_->lp);
 
+  double prev_lower_bound = mipdata_->lower_bound;
+
   mipdata_->lower_bound = mipdata_->nodequeue.getBestLowerBound();
 
+  bool bound_change = mipdata_->lower_bound != prev_lower_bound;
+  if (!submip && bound_change)
+    mipdata_->updatePrimalDualIntegral(prev_lower_bound, mipdata_->lower_bound,
+                                       mipdata_->upper_bound,
+                                       mipdata_->upper_bound);
+
   mipdata_->printDisplayLine();
   search.installNode(mipdata_->nodequeue.popBestBoundNode());
   int64_t numStallNodes = 0;
@@ -348,8 +356,16 @@ void HighsMipSolver::run() {
     search.flushStatistics();
 
     if (limit_reached) {
+      double prev_lower_bound = mipdata_->lower_bound;
+
       mipdata_->lower_bound = std::min(mipdata_->upper_bound,
                                        mipdata_->nodequeue.getBestLowerBound());
+
+      bool bound_change = mipdata_->lower_bound != prev_lower_bound;
+      if (!submip && bound_change)
+        mipdata_->updatePrimalDualIntegral(
+            prev_lower_bound, mipdata_->lower_bound, mipdata_->upper_bound,
+            mipdata_->upper_bound);
       mipdata_->printDisplayLine();
       break;
     }
@@ -371,13 +387,29 @@ void HighsMipSolver::run() {
     if (mipdata_->domain.infeasible()) {
       mipdata_->nodequeue.clear();
       mipdata_->pruned_treeweight = 1.0;
+
+      double prev_lower_bound = mipdata_->lower_bound;
+
       mipdata_->lower_bound = std::min(kHighsInf, mipdata_->upper_bound);
+
+      bool bound_change = mipdata_->lower_bound != prev_lower_bound;
+      if (!submip && bound_change)
+        mipdata_->updatePrimalDualIntegral(
+            prev_lower_bound, mipdata_->lower_bound, mipdata_->upper_bound,
+            mipdata_->upper_bound);
       mipdata_->printDisplayLine();
       break;
     }
 
+    double prev_lower_bound = mipdata_->lower_bound;
+
     mipdata_->lower_bound = std::min(mipdata_->upper_bound,
                                      mipdata_->nodequeue.getBestLowerBound());
+    bool bound_change = mipdata_->lower_bound != prev_lower_bound;
+    if (!submip && bound_change)
+      mipdata_->updatePrimalDualIntegral(
+          prev_lower_bound, mipdata_->lower_bound, mipdata_->upper_bound,
+          mipdata_->upper_bound);
     mipdata_->printDisplayLine();
     if (mipdata_->nodequeue.empty()) break;
 
@@ -538,7 +570,16 @@ void HighsMipSolver::run() {
         if (mipdata_->domain.infeasible()) {
           mipdata_->nodequeue.clear();
           mipdata_->pruned_treeweight = 1.0;
+
+          double prev_lower_bound = mipdata_->lower_bound;
+
           mipdata_->lower_bound = std::min(kHighsInf, mipdata_->upper_bound);
+
+          bool bound_change = mipdata_->lower_bound != prev_lower_bound;
+          if (!submip && bound_change)
+            mipdata_->updatePrimalDualIntegral(
+                prev_lower_bound, mipdata_->lower_bound, mipdata_->upper_bound,
+                mipdata_->upper_bound);
           break;
         }
 
@@ -547,9 +588,16 @@ void HighsMipSolver::run() {
           break;
         }
 
+        double prev_lower_bound = mipdata_->lower_bound;
+
         mipdata_->lower_bound = std::min(
             mipdata_->upper_bound, mipdata_->nodequeue.getBestLowerBound());
 
+        bool bound_change = mipdata_->lower_bound != prev_lower_bound;
+        if (!submip && bound_change)
+          mipdata_->updatePrimalDualIntegral(
+              prev_lower_bound, mipdata_->lower_bound, mipdata_->upper_bound,
+              mipdata_->upper_bound);
         mipdata_->printDisplayLine();
 
         if (!mipdata_->domain.getChangedCols().empty()) {
@@ -579,7 +627,16 @@ void HighsMipSolver::run() {
         search.openNodesToQueue(mipdata_->nodequeue);
         mipdata_->nodequeue.clear();
         mipdata_->pruned_treeweight = 1.0;
+
+        double prev_lower_bound = mipdata_->lower_bound;
+
         mipdata_->lower_bound = std::min(kHighsInf, mipdata_->upper_bound);
+
+        bool bound_change = mipdata_->lower_bound != prev_lower_bound;
+        if (!submip && bound_change)
+          mipdata_->updatePrimalDualIntegral(
+              prev_lower_bound, mipdata_->lower_bound, mipdata_->upper_bound,
+              mipdata_->upper_bound);
         break;
       }
 
@@ -609,12 +666,20 @@ void HighsMipSolver::run() {
 }
 
 void HighsMipSolver::cleanupSolve() {
+  // Force a final logging line
   mipdata_->printDisplayLine(kSolutionSourceCleanup);
   // Stop the solve clock - which won't be running if presolve
   // determines the model status
   if (analysis_.mipTimerRunning(kMipClockSolve))
     analysis_.mipTimerStop(kMipClockSolve);
+
+  // Need to complete the calculation of P-D integral, checking for NO
+  // gap change
+  mipdata_->updatePrimalDualIntegral(
+      mipdata_->lower_bound, mipdata_->lower_bound, mipdata_->upper_bound,
+      mipdata_->upper_bound, false);
   analysis_.mipTimerStart(kMipClockPostsolve);
+
   bool havesolution = solution_objective_ != kHighsInf;
   bool feasible;
   if (havesolution)
@@ -639,6 +704,7 @@ void HighsMipSolver::cleanupSolve() {
   node_count_ = mipdata_->num_nodes;
   total_lp_iterations_ = mipdata_->total_lp_iterations;
   dual_bound_ = std::min(dual_bound_, primal_bound_);
+  primal_dual_integral_ = mipdata_->primal_dual_integral.value;
 
   // adjust objective sense in case of maximization problem
   if (orig_model_->sense_ == ObjSense::kMaximize) {
@@ -716,9 +782,11 @@ void HighsMipSolver::cleanupSolve() {
                "  Primal bound      %.12g\n"
                "  Dual bound        %.12g\n"
                "  Gap               %s\n"
+               "  P-D integral      %.12g\n"
                "  Solution status   %s\n",
                utilModelStatusToString(modelstatus_).c_str(), primal_bound_,
-               dual_bound_, gapString.data(), solutionstatus.c_str());
+               dual_bound_, gapString.data(),
+               mipdata_->primal_dual_integral.value, solutionstatus.c_str());
   if (solutionstatus != "-")
     highsLogUser(options_mip_->log_options, HighsLogType::kInfo,
                  "                    %.12g (objective)\n"

src/mip/HighsMipSolver.h

@@ -41,6 +41,7 @@ class HighsMipSolver {
   double gap_;
   int64_t node_count_;
   int64_t total_lp_iterations_;
+  double primal_dual_integral_;
 
   FILE* improving_solution_file_;
   std::vector<HighsObjectiveSolution> saved_objective_and_solution_;