Merge branch 'hi-pdlp-jh' into hi-pdlp

Author: Yanyu000

FEATURES.md

@@ -14,3 +14,4 @@ Prompted by [#2463](https://github.com/ERGO-Code/HiGHS/issues/2463), the HiGHS s
 
 Only for LPs is there a choice of solver. Previously, when setting the `solver` option to anything other than "choose", any incumbent model was solved as an LP, using that LP solver. This has caused confusiuon for users, and is unnecessary now that there is the `solve_relaxation` option. Now, if the incumbent model is a QP or MIP, it is solved as such (unless `solve_relaxation` is true for a MIP), and the value of the `solver` option only determines what solver is used to solve an LP. If the value of `solver` is "choose", then HiGHS will use what it expects to be the best solver for the problem; if value of `solver` is "ipm", then HiGHS will use what it expects to be the better IPM solver (of HiPO and IPX) for the problem; if value of `solver` is "hipo", then HiGHS will use the HiPO IPM solver (if available in the build); if value of `solver` is "ipx", then HiGHS will use the IPX IPM solver; if value of `solver` is "pdlp", then HiGHS will use the PDLP first-order solver. The option `mip_lp_solver` has been introduced to define which LP solver is used when solving LPs in the MIP solver for which an advanced basis is not known - typically the "root node" LP. Note that The PDLP solver cannot be used to solve such LPs, since it does not yield a basic solution. If an interior point solver fails to obtain a basic solution, the simplex solver will then be used. The option `mip_ipm_solver` has been introduced to define which IPM solver is used when solving LPs in the MIP solver for which IPM is mandatory - typically the analytic centre calculation. When LPs are to be solved by an IPM solver, the HiPO solver is used (if available in the build) unless IPX has been specified explicitly. 
 
+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

@@ -208,3 +208,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);
+}

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

@@ -822,13 +822,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)
@@ -842,17 +844,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/lp_solver.cc

@@ -522,8 +522,8 @@ void LpSolver::BuildStartingBasis() {
         info_.status_ipm = IPX_STATUS_debug;
         return;
     }
-    basis_.reset(new Basis(control_, model_));
     control_.hLog(" Constructing starting basis...\n");
+    basis_.reset(new Basis(control_, model_));
     StartingBasis(iterate_.get(), basis_.get(), &info_);
     if (info_.errflag == IPX_ERROR_user_interrupt) {
         info_.errflag = 0;