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

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.

README.md

@@ -46,7 +46,7 @@ linear optimization problems of the form
 
 $$ \min \quad \dfrac{1}{2}x^TQx + c^Tx \qquad \textrm{s.t.}~ \quad L \leq Ax \leq U; \quad l \leq x \leq u $$
 
-where Q must be positive semi-definite and, if Q is zero, there may be a requirement that some of the variables take integer values. Thus HiGHS can solve linear programming (LP) problems, convex quadratic programming (QP) problems, and mixed integer programming (MIP) problems. It is mainly written in C++, but also has some C. It has been developed and tested on various Linux, MacOS and Windows installations. No third-party dependencies are required.
+where $Q$ must be positive semi-definite and, if $Q$ is zero, there may be a requirement that some of the variables take integer values. Thus HiGHS can solve linear programming (LP) problems, convex quadratic programming (QP) problems, and mixed integer programming (MIP) problems. It is mainly written in C++, but also has some C. It has been developed and tested on various Linux, MacOS and Windows installations. No third-party dependencies are required.
 
 HiGHS has primal and dual revised simplex solvers, originally written by Qi Huangfu and further developed by Julian Hall. It also has an interior point solver for LP written by Lukas Schork, an active set solver for QP written by Michael Feldmeier, and a MIP solver written by Leona Gottwald. Other features have been added by Julian Hall and Ivet Galabova, who manages the software engineering of HiGHS and interfaces to C, C#, FORTRAN, Julia and Python.
 

check/TestBasis.cpp

@@ -315,3 +315,43 @@ void testBasisRestart(Highs& highs, const std::string& basis_file,
 
   REQUIRE(info.simplex_iteration_count == 0);
 }
+
+TEST_CASE("Basis-read", "[highs_basis_data]") {
+  // Duplicates test_read_basis in test_highspy.py
+  const std::string test_name = Catch::getResultCapture().getCurrentTestName();
+
+  HighsLp lp;
+  lp.num_col_ = 2;
+  lp.num_row_ = 2;
+  lp.col_cost_ = {0, 1};
+  lp.col_lower_.assign(lp.num_col_, -kHighsInf);
+  lp.col_upper_.assign(lp.num_col_, kHighsInf);
+  lp.row_lower_ = {2, 0};
+  lp.row_upper_.assign(lp.num_row_, kHighsInf);
+  lp.a_matrix_.start_ = {0, 2, 4};
+  lp.a_matrix_.index_ = {0, 1, 0, 1};
+  lp.a_matrix_.value_ = {-1, 1, 1, 1};
+
+  HighsBasisStatus status_before = HighsBasisStatus::kNonbasic;
+  HighsBasisStatus status_after = HighsBasisStatus::kBasic;
+  Highs h1;
+  const HighsBasis& basis1 = h1.getBasis();
+  h1.passModel(lp);
+  REQUIRE(basis1.col_status[0] == status_before);
+  h1.run();
+  REQUIRE(basis1.col_status[0] == status_after);
+
+  Highs h2;
+  const HighsBasis& basis2 = h2.getBasis();
+  h2.passModel(lp);
+  REQUIRE(basis2.col_status[0] == status_before);
+
+  const std::string basis_file = test_name + ".bas";
+  h1.writeBasis(basis_file);
+  h2.readBasis(basis_file);
+  REQUIRE(basis2.col_status[0] == status_after);
+
+  std::remove(basis_file.c_str());
+  h1.resetGlobalScheduler(true);
+  h2.resetGlobalScheduler(true);
+}

check/TestIpm.cpp

@@ -200,3 +200,19 @@ TEST_CASE("test-2087", "[highs_ipm]") {
 
   h.resetGlobalScheduler(true);
 }
+
+TEST_CASE("test-2527", "[highs_ipm]") {
+  std::string filename =
+      std::string(HIGHS_DIR) + "/check/instances/primal1.mps";
+  Highs h;
+  //  h.setOptionValue("output_flag", dev_run);
+  REQUIRE(h.readModel(filename) == HighsStatus::kOk);
+  HighsLp lp = h.getLp();
+  lp.col_cost_.assign(lp.num_col_, 0);
+  REQUIRE(h.passModel(lp) == HighsStatus::kOk);
+  h.setOptionValue("solver", kIpmString);
+  h.setOptionValue("presolve", kHighsOffString);
+  REQUIRE(h.run() == HighsStatus::kOk);
+
+  h.resetGlobalScheduler(true);
+}

check/TestPresolve.cpp

@@ -838,3 +838,27 @@ TEST_CASE("presolve-egout-ac", "[highs_test_presolve]") {
 
   h.resetGlobalScheduler(true);
 }
+
+TEST_CASE("dual-bound-tightening", "[highs_test_presolve]") {
+  std::string model_file =
+      std::string(HIGHS_DIR) + "/check/instances/gesa2.mps";
+
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+  highs.readModel(model_file);
+
+  // complement variables to get code coverage
+  HighsLp lp = highs.getLp();
+  std::transform(lp.a_matrix_.value_.begin(), lp.a_matrix_.value_.end(),
+                 lp.a_matrix_.value_.begin(), [](double v) { return -v; });
+  std::transform(lp.col_cost_.begin(), lp.col_cost_.end(), lp.col_cost_.begin(),
+                 [](double v) { return -v; });
+  std::transform(lp.col_upper_.begin(), lp.col_upper_.end(),
+                 lp.col_upper_.begin(), [](double v) { return -v; });
+  std::transform(lp.col_lower_.begin(), lp.col_lower_.end(),
+                 lp.col_lower_.begin(), [](double v) { return -v; });
+  std::swap(lp.col_lower_, lp.col_upper_);
+
+  highs.passModel(lp);
+  REQUIRE(highs.presolve() == HighsStatus::kOk);
+}

highs/interfaces/highs_c_api.h

@@ -400,7 +400,9 @@ HighsInt Highs_presolve(void* highs);
 HighsInt Highs_run(void* highs);
 
 /**
- * Postsolve a model using a primal (and possibly dual) solution.
+ * Postsolve a model using a primal (and possibly dual) solution. The
+ * postsolved solution can be retrieved later by calling
+ * `Highs_getSolution`.
  *
  * @param highs       A pointer to the Highs instance.
  * @param col_value   An array of length [num_col] with the column solution

highs/ipm/ipx/basis.cc

@@ -119,12 +119,21 @@ Int Basis::Factorize() {
 
     // Build column pointers for passing to LU factorization.
     std::vector<Int> begin(m), end(m);
+    Int basis_num_nz = 0;
     for (Int i = 0; i < m; i++) {
         assert(basis_[i] >= 0);
         begin[i] = AI.begin(basis_[i]);
         end[i] = AI.end(basis_[i]);
+	basis_num_nz += (end[i]-begin[i]);
     }
 
+    std::stringstream h_logging_stream;
+    h_logging_stream.str(std::string());
+    h_logging_stream <<
+      " Start  factorization " << num_factorizations_+1 <<
+      ": nonzeros in basis = " << basis_num_nz << "\n";
+    control_.hIntervalLog(h_logging_stream);
+
     Int err = 0;                // return code
     while (true) {
         Int flag = lu_->Factorize(begin.data(), end.data(), AI.rowidx(),
@@ -151,6 +160,11 @@ Int Basis::Factorize() {
     }
     time_factorize_ += timer.Elapsed();
     factorization_is_fresh_ = true;
+    h_logging_stream.str(std::string());
+    h_logging_stream <<
+      " Finish factorization " << num_factorizations_ <<
+      ": fill factor = " << lu_->fill_factor() << "\n";
+    control_.hIntervalLog(h_logging_stream);
     return err;
 }
 

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.