Merge pull request #2064 from ERGO-Code/simplex-stats

Simplex stats

check/TestLpSolvers.cpp

@@ -29,6 +29,7 @@ void testDualObjective(const std::string model) {
       std::max(1.0, std::fabs(primal_objective));
   REQUIRE(relative_primal_dual_gap < 1e-12);
 }
+
 void testSolver(Highs& highs, const std::string solver,
                 IterationCount& default_iteration_count,
                 const HighsInt int_simplex_strategy = 0) {
@@ -641,3 +642,40 @@ TEST_CASE("standard-form-lp", "[highs_lp_solver]") {
         "maximizing\n");
   testStandardForm(highs.getLp());
 }
+
+TEST_CASE("simplex-stats", "[highs_lp_solver]") {
+  HighsStatus return_status;
+
+  Highs h;
+  const HighsSimplexStats& simplex_stats = h.getSimplexStats();
+  h.setOptionValue("output_flag", dev_run);
+  std::string model_file =
+      std::string(HIGHS_DIR) + "/check/instances/adlittle.mps";
+  REQUIRE(h.readModel(model_file) == HighsStatus::kOk);
+
+  REQUIRE(h.run() == HighsStatus::kOk);
+  REQUIRE(simplex_stats.valid);
+  REQUIRE(simplex_stats.iteration_count == 0);
+  REQUIRE(simplex_stats.num_invert == 1);
+  REQUIRE(simplex_stats.last_invert_num_el > 0);
+  REQUIRE(simplex_stats.last_factored_basis_num_el > 0);
+  REQUIRE(simplex_stats.col_aq_density == 0);
+  REQUIRE(simplex_stats.row_ep_density == 0);
+  REQUIRE(simplex_stats.row_ap_density == 0);
+  REQUIRE(simplex_stats.row_DSE_density == 0);
+  if (dev_run) h.reportSimplexStats(stdout);
+
+  h.clearSolver();
+  h.setOptionValue("presolve", kHighsOffString);
+  REQUIRE(h.run() == HighsStatus::kOk);
+  REQUIRE(simplex_stats.valid);
+  REQUIRE(simplex_stats.iteration_count > 0);
+  REQUIRE(simplex_stats.num_invert > 0);
+  REQUIRE(simplex_stats.last_invert_num_el > 0);
+  REQUIRE(simplex_stats.last_factored_basis_num_el > 0);
+  REQUIRE(simplex_stats.col_aq_density > 0);
+  REQUIRE(simplex_stats.row_ep_density > 0);
+  REQUIRE(simplex_stats.row_ap_density > 0);
+  REQUIRE(simplex_stats.row_DSE_density > 0);
+  if (dev_run) h.reportSimplexStats(stdout);
+}

src/Highs.h

@@ -1215,6 +1215,14 @@ class Highs {
   static void resetGlobalScheduler(bool blocking = false);
 
   // Start of advanced methods for HiGHS MIP solver
+
+  const HighsSimplexStats& getSimplexStats() const {
+    return ekk_instance_.getSimplexStats();
+  }
+  void reportSimplexStats(FILE* file) const {
+    ekk_instance_.reportSimplexStats(file);
+  }
+
   /**
    * @brief Get the hot start basis data from the most recent simplex
    * solve. Advanced method: for HiGHS MIP solver

src/lp_data/HStruct.h

@@ -154,4 +154,18 @@ struct HighsLinearObjective {
   void clear();
 };
 
+struct HighsSimplexStats {
+  bool valid;
+  HighsInt iteration_count;
+  HighsInt num_invert;
+  HighsInt last_invert_num_el;
+  HighsInt last_factored_basis_num_el;
+  double col_aq_density;
+  double row_ep_density;
+  double row_ap_density;
+  double row_DSE_density;
+  void report(FILE* file, const std::string message = "") const;
+  void initialise(const HighsInt iteration_count_ = 0);
+};
+
 #endif /* LP_DATA_HSTRUCT_H_ */

src/lp_data/HighsSolve.cpp

@@ -17,7 +17,7 @@
 #include "pdlp/CupdlpWrapper.h"
 #include "simplex/HApp.h"
 
-// The method below runs simplex or ipx solver on the lp.
+// The method below runs simplex, ipx or pdlp solver on the lp.
 HighsStatus solveLp(HighsLpSolverObject& solver_object, const string message) {
   HighsStatus return_status = HighsStatus::kOk;
   HighsStatus call_status;

src/simplex/HApp.h

@@ -112,6 +112,9 @@ inline HighsStatus solveLpSimplex(HighsLpSolverObject& solver_object) {
   // return
   resetModelStatusAndHighsInfo(solver_object);
 
+  // Initialise the simplex stats
+  ekk_instance.initialiseSimplexStats();
+
   // Assumes that the LP has a positive number of rows, since
   // unconstrained LPs should be solved in solveLp
   bool positive_num_row = solver_object.lp_.num_row_ > 0;

src/simplex/HEkk.cpp

@@ -296,6 +296,7 @@ void HEkk::invalidate() {
   assert(!this->status_.is_permuted);
   this->status_.initialised_for_solve = false;
   this->invalidateBasisMatrix();
+  this->simplex_stats_.initialise();
 }
 
 void HEkk::invalidateBasisMatrix() {
@@ -2082,6 +2083,7 @@ HighsInt HEkk::computeFactor() {
   // number of updates shouldn't be positive
   info_.update_count = 0;
 
+  simplex_stats_.num_invert++;
   return rank_deficiency;
 }
 
@@ -3500,7 +3502,19 @@ HighsStatus HEkk::returnFromEkkSolve(const HighsStatus return_status) {
   // Note that in timeReporting(1), analysis_.analyse_simplex_time
   // reverts to its value given by options_
   if (analysis_.analyse_simplex_time) analysis_.reportSimplexTimer();
-
+  simplex_stats_.valid = true;
+  // Since HEkk::iteration_count_ includes iteration on presolved LP,
+  // simplex_stats_.iteration_count is initialised to -
+  // HEkk::iteration_count_
+  simplex_stats_.iteration_count += iteration_count_;
+  //  simplex_stats_.num_invert is incremented internally
+  simplex_stats_.last_invert_num_el = simplex_nla_.factor_.invert_num_el;
+  simplex_stats_.last_factored_basis_num_el =
+      simplex_nla_.factor_.basis_matrix_num_el;
+  simplex_stats_.col_aq_density = analysis_.col_aq_density;
+  simplex_stats_.row_ep_density = analysis_.row_ep_density;
+  simplex_stats_.row_ap_density = analysis_.row_ap_density;
+  simplex_stats_.row_DSE_density = analysis_.row_DSE_density;
   return return_status;
 }
 
@@ -4406,3 +4420,30 @@ void HEkk::unitBtranResidual(const HighsInt row_out, const HVector& row_ep,
     residual_norm = max(fabs(residual.array[iRow]), residual_norm);
   }
 }
+
+void HighsSimplexStats::report(FILE* file, std::string message) const {
+  fprintf(file, "\nSimplex stats: %s\n", message.c_str());
+  fprintf(file, "   valid                      = %d\n", this->valid);
+  fprintf(file, "   iteration_count            = %d\n", this->iteration_count);
+  fprintf(file, "   num_invert                 = %d\n", this->num_invert);
+  fprintf(file, "   last_invert_num_el         = %d\n",
+          this->last_invert_num_el);
+  fprintf(file, "   last_factored_basis_num_el = %d\n",
+          this->last_factored_basis_num_el);
+  fprintf(file, "   col_aq_density             = %g\n", this->col_aq_density);
+  fprintf(file, "   row_ep_density             = %g\n", this->row_ep_density);
+  fprintf(file, "   row_ap_density             = %g\n", this->row_ap_density);
+  fprintf(file, "   row_DSE_density            = %g\n", this->row_DSE_density);
+}
+
+void HighsSimplexStats::initialise(const HighsInt iteration_count_) {
+  valid = false;
+  iteration_count = -iteration_count_;
+  num_invert = 0;
+  last_invert_num_el = 0;
+  last_factored_basis_num_el = 0;
+  col_aq_density = 0;
+  row_ep_density = 0;
+  row_ap_density = 0;
+  row_DSE_density = 0;
+}

src/simplex/HEkk.h

@@ -156,6 +156,12 @@ class HEkk {
                            const vector<double>& rowLower,
                            const vector<double>& rowUpper);
 
+  const HighsSimplexStats& getSimplexStats() const { return simplex_stats_; }
+  void initialiseSimplexStats() { simplex_stats_.initialise(iteration_count_); }
+  void reportSimplexStats(FILE* file, const std::string message = "") const {
+    simplex_stats_.report(file, message);
+  }
+
   // Make this private later
   void chooseSimplexStrategyThreads(const HighsOptions& options,
                                     HighsSimplexInfo& info);
@@ -255,6 +261,8 @@ class HEkk {
   std::vector<HighsSimplexBadBasisChangeRecord> bad_basis_change_;
   std::vector<double> primal_phase1_dual_;
 
+  HighsSimplexStats simplex_stats_;
+
  private:
   bool isUnconstrainedLp();
   void initialiseForSolve();

src/simplex/HighsSimplexAnalysis.h

@@ -251,6 +251,7 @@ class HighsSimplexAnalysis {
   void reportFactorTimer();
   void updateInvertFormData(const HFactor& factor);
   void reportInvertFormData();
+  HighsInt numInvert() { return num_invert; }
 
   // Control methods to be moved to HEkkControl
   void dualSteepestEdgeWeightError(const double computed_edge_weight,