Introducing a deterministic decision on whether to compute the analytic centre in MIP

cmake/sources-python.cmake

@@ -311,6 +311,7 @@ set(highs_headers_python
     src/lp_data/HighsRanging.h
     src/lp_data/HighsSolution.h
     src/lp_data/HighsSolutionDebug.h
+    src/lp_data/HighsSolutionStats.h
     src/lp_data/HighsSolve.h
     src/lp_data/HighsStatus.h
     src/lp_data/HStruct.h

cmake/sources.cmake

@@ -315,6 +315,7 @@ set(highs_headers
     lp_data/HighsRanging.h
     lp_data/HighsSolution.h
     lp_data/HighsSolutionDebug.h
+    lp_data/HighsSolutionStats.h
     lp_data/HighsSolve.h
     lp_data/HighsStatus.h
     lp_data/HStruct.h

src/Highs.h

@@ -1219,6 +1219,9 @@ class Highs {
   const HighsSimplexStats& getSimplexStats() const {
     return ekk_instance_.getSimplexStats();
   }
+  void passSimplexStats(const HighsSimplexStats simplex_stats) {
+    return ekk_instance_.passSimplexStats(simplex_stats);
+  }
   void reportSimplexStats(FILE* file) const {
     ekk_instance_.reportSimplexStats(file);
   }

src/ipm/IpxWrapper.cpp

@@ -23,14 +23,15 @@ using std::min;
 
 HighsStatus solveLpIpx(HighsLpSolverObject& solver_object) {
   return solveLpIpx(solver_object.options_, solver_object.timer_,
-                    solver_object.lp_, solver_object.basis_,
-                    solver_object.solution_, solver_object.model_status_,
-                    solver_object.highs_info_, solver_object.callback_);
+                    solver_object.lp_, solver_object.ekk_instance_,
+                    solver_object.basis_, solver_object.solution_,
+                    solver_object.model_status_, solver_object.highs_info_,
+                    solver_object.callback_);
 }
 
 HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer,
-                       const HighsLp& lp, HighsBasis& highs_basis,
-                       HighsSolution& highs_solution,
+                       const HighsLp& lp, const HEkk& ekk_instance,
+                       HighsBasis& highs_basis, HighsSolution& highs_solution,
                        HighsModelStatus& model_status, HighsInfo& highs_info,
                        HighsCallback& callback) {
   // Use IPX to try to solve the LP
@@ -121,6 +122,9 @@ HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer,
   parameters.time_limit = options.time_limit - timer.readRunHighsClock();
   parameters.ipm_maxiter =
       options.ipm_iteration_limit - highs_info.ipm_iteration_count;
+  parameters.cr1_maxiter = options.cr1_iteration_limit;
+  parameters.cr2_maxiter = options.cr2_iteration_limit;
+  parameters.kkt_logging = options.kkt_logging;
   // Determine if crossover is to be run or not
   //
   // When doing analytic centring calculations, crossover must not be
@@ -146,6 +150,8 @@ HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer,
   parameters.max_centring_steps = options.max_centring_steps;
   parameters.centring_ratio_tolerance = options.centring_ratio_tolerance;
 
+  parameters.simplex_stats = ekk_instance.getSimplexStats();
+
   // Set the internal IPX parameters
   lps.SetParameters(parameters);
 
@@ -183,6 +189,8 @@ HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer,
   if (report_solve_data) reportSolveData(options.log_options, ipx_info);
   highs_info.ipm_iteration_count += (HighsInt)ipx_info.iter;
   highs_info.crossover_iteration_count += (HighsInt)ipx_info.updates_crossover;
+  highs_info.max_cr_iteration_count1 = ipx_info.kkt_iter_max1;
+  highs_info.max_cr_iteration_count2 = ipx_info.kkt_iter_max2;
 
   // If not solved...
   if (solve_status != IPX_STATUS_solved) {
@@ -960,6 +968,10 @@ void reportSolveData(const HighsLogOptions& log_options,
               (int)ipx_info.kktiter1);
   highsLogDev(log_options, HighsLogType::kInfo, "    KKT iter 2 = %d\n",
               (int)ipx_info.kktiter2);
+  highsLogDev(log_options, HighsLogType::kInfo, "    KKT iter max 1 = %d\n",
+              (int)ipx_info.kkt_iter_max1);
+  highsLogDev(log_options, HighsLogType::kInfo, "    KKT iter max 2 = %d\n",
+              (int)ipx_info.kkt_iter_max2);
   highsLogDev(log_options, HighsLogType::kInfo, "    Basis repairs = %d\n",
               (int)ipx_info.basis_repairs);
   highsLogDev(log_options, HighsLogType::kInfo, "    Updates start     = %d\n",

src/ipm/IpxWrapper.h

@@ -25,8 +25,8 @@
 HighsStatus solveLpIpx(HighsLpSolverObject& solver_object);
 
 HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer,
-                       const HighsLp& lp, HighsBasis& highs_basis,
-                       HighsSolution& highs_solution,
+                       const HighsLp& lp, const HEkk& ekk_instance,
+                       HighsBasis& highs_basis, HighsSolution& highs_solution,
                        HighsModelStatus& model_status, HighsInfo& highs_info,
                        HighsCallback& callback);
 

src/ipm/ipx/basis.cc

@@ -449,6 +449,13 @@ double Basis::time_update() const {
     return time_update_;
 }
 
+double Basis::current_fill() const {
+    if (fill_factors_.empty())
+        return 0.0;
+    Int num_factors = fill_factors_.size();
+    return fill_factors_[num_factors-1];
+}
+
 double Basis::mean_fill() const {
     if (fill_factors_.empty())
         return 0.0;

src/ipm/ipx/basis.h

@@ -219,6 +219,7 @@ class Basis {
     double time_ftran() const;        // time FTRAN, including partial
     double time_btran() const;        // time BTRAN, including partial
     double time_update() const;       // time LU update
+    double current_fill() const;      // Current LU fill factors
     double mean_fill() const;         // geom. mean of LU fill factors
     double max_fill() const;          // max LU fill factor
 

src/ipm/ipx/conjugate_residuals.cc

@@ -12,6 +12,9 @@ ConjugateResiduals::ConjugateResiduals(const Control& control) :
 void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs,
                               double tol, const double* resscale, Int maxiter,
                               Vector& lhs) {
+
+    const bool cr_logging = control_.kkt_logging();
+
     const Int m = rhs.size();
     Vector residual(m);  // rhs - C*lhs
     Vector step(m);      // update to lhs
@@ -23,9 +26,12 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs,
     errflag_ = 0;
     iter_ = 0;
     time_ = 0.0;
+    Int use_maxiter = maxiter;
     if (maxiter < 0)
-        maxiter = m+100;
+        use_maxiter = m+100;
 
+    if (cr_logging) printf("\nCR(C) maxiter = (entry = %d, use = %d)\n", int(maxiter), int(use_maxiter)); fflush(stdout);
+    maxiter = use_maxiter;
     // Initialize residual, step and Cstep.
     if (Infnorm(lhs) == 0.0) {
         residual = rhs;         // saves a matrix-vector op
@@ -37,9 +43,10 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs,
     step = residual;
     Cstep = Cresidual;
 
+    double resnorm;
     while (true) {
         // Termination check.
-        double resnorm = 0.0;
+        resnorm = 0.0;
         if (resscale)
             for (Int i = 0; i < m; i++)
                 resnorm = std::max(resnorm, std::abs(resscale[i]*residual[i]));
@@ -48,6 +55,7 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs,
         if (resnorm <= tol)
             break;
         if (iter_ == maxiter) {
+ 	  if (cr_logging) printf("CR(C) reached maxiter!\n\n"); fflush(stdout);
             control_.Debug(3)
                 << " CR method not converged in " << maxiter << " iterations."
                 << " residual = " << sci2(resnorm) << ','
@@ -59,7 +67,6 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs,
             errflag_ = IPX_ERROR_cr_matrix_not_posdef;
             break;
         }
-
         // Update lhs, residual and Cresidual.
         const double denom = Dot(Cstep,Cstep);
         const double alpha = cdot/denom;
@@ -79,15 +86,18 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs,
         cdot = cdotnew;
 
         iter_++;
-        if ((errflag_ = control_.InterruptCheck()) != 0)
-            break;
+        if ((errflag_ = control_.InterruptCheck()) != 0) break;
     }
+    if (errflag_ != IPX_ERROR_cr_iter_limit)
+      if (cr_logging) printf("CR(C)    iter = %d; resnorm = %g <= %g = tol? %s\n\n", int(iter_), resnorm, tol, resnorm <= tol ? "T" : "F"); fflush(stdout);
     time_ = timer.Elapsed();
 }
 
 void ConjugateResiduals::Solve(LinearOperator& C, LinearOperator& P,
                               const Vector& rhs, double tol,
                               const double* resscale, Int maxiter, Vector& lhs){
+    const bool cr_logging = control_.kkt_logging();
+
     const Int m = rhs.size();
     Vector residual(m);   // rhs - C*lhs
     Vector sresidual(m);  // preconditioned residual
@@ -109,9 +119,12 @@ void ConjugateResiduals::Solve(LinearOperator& C, LinearOperator& P,
     errflag_ = 0;
     iter_ = 0;
     time_ = 0.0;
+    Int use_maxiter = maxiter;
     if (maxiter < 0)
-        maxiter = m+100;
+        use_maxiter = m+100;
 
+    if (cr_logging) printf("\nCR(C,P) maxiter = (entry = %d, use = %d)\n", int(maxiter), int(use_maxiter)); fflush(stdout);
+    maxiter = use_maxiter;
     // Initialize residual, sresidual, step and Cstep.
     if (Infnorm(lhs) == 0.0) {
         residual = rhs;         // saves a matrix-vector op
@@ -124,9 +137,10 @@ void ConjugateResiduals::Solve(LinearOperator& C, LinearOperator& P,
     step = sresidual;
     Cstep = Csresidual;
 
+    double resnorm;
     while (true) {
         // Termination check.
-        double resnorm = 0.0;
+        resnorm = 0.0;
         if (resscale)
             for (Int i = 0; i < m; i++)
                 resnorm = std::max(resnorm, std::abs(resscale[i]*residual[i]));
@@ -135,7 +149,8 @@ void ConjugateResiduals::Solve(LinearOperator& C, LinearOperator& P,
         if (resnorm <= tol)
             break;
         if (iter_ == maxiter) {
-            control_.Debug(3)
+ 	  if (cr_logging) printf("CR(C,P) reached maxiter!\n\n"); fflush(stdout);
+           control_.Debug(3)
                 << " PCR method not converged in " << maxiter << " iterations."
                 << " residual = " << sci2(resnorm) << ','
                 << " tolerance = " << sci2(tol) << '\n';
@@ -207,6 +222,8 @@ void ConjugateResiduals::Solve(LinearOperator& C, LinearOperator& P,
         if ((errflag_ = control_.InterruptCheck()) != 0)
             break;
     }
+    if (errflag_ != IPX_ERROR_cr_iter_limit)
+      if (cr_logging) printf("CR(C,P)    iter = %d; resnorm = %g <= %g = tol? %s\n\n", int(iter_), resnorm, tol, resnorm <= tol ? "T" : "F"); fflush(stdout);
     time_ = timer.Elapsed();
 }
 

src/ipm/ipx/control.cc

@@ -13,8 +13,10 @@ Control::Control() {
 Int Control::InterruptCheck(const Int ipm_iteration_count) const {
     HighsTaskExecutor::getThisWorkerDeque()->checkInterrupt();
     if (parameters_.time_limit >= 0.0 &&
-        parameters_.time_limit < timer_.Elapsed())
+        parameters_.time_limit < timer_.Elapsed()) {
+      printf("Control::InterruptCheck Reached time limit of %g\n", parameters_.time_limit);
         return IPX_ERROR_time_interrupt;
+    }
     // The pointer callback_ should not be null, since that indicates
     // that it's not been set
     assert(callback_);

src/ipm/ipx/control.h

@@ -68,7 +68,10 @@ class Control {
     double ipm_optimality_tol() const { return parameters_.ipm_optimality_tol; }
     double ipm_drop_primal() const { return parameters_.ipm_drop_primal; }
     double ipm_drop_dual() const { return parameters_.ipm_drop_dual; }
+    bool kkt_logging() const { return parameters_.kkt_logging; }
     double kkt_tol() const { return parameters_.kkt_tol; }
+    ipxint cr1_maxiter() const { return parameters_.cr1_maxiter; }
+    ipxint cr2_maxiter() const { return parameters_.cr2_maxiter; }
     ipxint crash_basis() const { return parameters_.crash_basis; }
     double dependency_tol() const { return parameters_.dependency_tol; }
     double volume_tol() const { return parameters_.volume_tol; }

src/ipm/ipx/info.cc

@@ -53,6 +53,8 @@ std::ostream& operator<<(std::ostream& os, const Info& info) {
     dump(os, "iter", info.iter);
     dump(os, "kktiter1", info.kktiter1);
     dump(os, "kktiter2", info.kktiter2);
+    dump(os, "kkt_iter_max1", info.kkt_iter_max1);
+    dump(os, "kkt_iter_max2", info.kkt_iter_max2);
     dump(os, "basis_repairs", info.basis_repairs);
     dump(os, "updates_start", info.updates_start);
     dump(os, "updates_ipm", info.updates_ipm);

src/ipm/ipx/ipm.cc

@@ -96,6 +96,14 @@ void IPM::Driver(KKTSolver* kkt, Iterate* iterate, Info* info) {
         if (info->errflag)
             break;
         MakeStep(step);
+	//
+	std::stringstream h_logging_stream;
+	h_logging_stream << "IPM::Driver fill factor = " <<
+	  kkt_->current_fill() << "; iter (sum = " <<
+	  kkt_->iterSum() << ", max = " <<
+	  kkt_->iterMax() << ")\n";
+        control_.hLog(h_logging_stream);
+	//
         info->iter++;
         PrintOutput();
     }
@@ -856,7 +864,7 @@ void IPM::PrintOutput() {
     control_.Debug()
       << "  " << Fixed(step_primal_, 4, 2) << " " << Fixed(step_dual_, 4, 2)
       << "  " << Format(kkt_->basis_changes(), 7)
-      << " "  << Format(kkt_->iter(), 7);
+      << " "  << Format(kkt_->iterSum(), 7);
     control_.Debug()
       << "  " << Format(info_->dual_dropped, 7)
       << " "  << Format(info_->primal_dropped, 7);

src/ipm/ipx/ipx_info.h

@@ -58,6 +58,8 @@ struct ipx_info {
     ipxint iter;                /* # interior point iterations */
     ipxint kktiter1;            /* # linear solver iterations before switch */
     ipxint kktiter2;            /* # linear solver iterations after switch */
+    ipxint kkt_iter_max1;       /* # max linear solver iterations before switch */
+    ipxint kkt_iter_max2;       /* # max linear solver iterations after switch */
     ipxint basis_repairs;       /* # basis repairs after crash, < 0 discarded */
     ipxint updates_start;       /* # basis updates for starting basis */
     ipxint updates_ipm;         /* # basis updates in IPM */

src/ipm/ipx/ipx_internal.h

@@ -31,7 +31,10 @@ struct Parameters : public ipx_parameters {
     ipm_optimality_tol = 1e-8;
     ipm_drop_primal = 1e-9;
     ipm_drop_dual = 1e-9;
+    kkt_logging = false;
     kkt_tol = 0.3;
+    cr1_maxiter = -1;
+    cr2_maxiter = -1;
     crash_basis = 1;
     dependency_tol = 1e-6;
     volume_tol = 2.0;

src/ipm/ipx/ipx_parameters.h

@@ -2,6 +2,7 @@
 #define IPX_PARAMETERS_H_
 
 #include "io/HighsIO.h"
+#include "lp_data/HighsSolutionStats.h"
 #include "ipm/ipx/ipx_config.h"
 #include <iostream>
 
@@ -28,6 +29,9 @@ struct ipx_parameters {
     double ipm_drop_dual;
 
     /* Linear solver */
+    bool kkt_logging;
+    ipxint cr1_maxiter;
+    ipxint cr2_maxiter;
     double kkt_tol;
 
     /* Basis construction in IPM */
@@ -65,7 +69,9 @@ struct ipx_parameters {
     /* HiGHS logging parameters */
     bool highs_logging;
     const HighsLogOptions* log_options;
-
+
+    /* Simplex solution stats */
+    HighsSimplexStats simplex_stats;
 };
 
 #ifdef __cplusplus

src/ipm/ipx/kkt_solver.cc

@@ -16,8 +16,10 @@ void KKTSolver::Solve(const Vector& a, const Vector& b, double tol,
     info->time_kkt_solve += timer.Elapsed();
 }
 
-Int KKTSolver::iter() const { return _iter(); }
+Int KKTSolver::iterSum() const { return _iterSum(); }
+Int KKTSolver::iterMax() const { return _iterMax(); }
 Int KKTSolver::basis_changes() const { return _basis_changes(); }
+double KKTSolver::current_fill() const {return _current_fill(); }
 const Basis* KKTSolver::basis() const { return _basis(); }
 
 }  // namespace ipx

src/ipm/ipx/kkt_solver.h

@@ -46,12 +46,19 @@ class KKTSolver {
     // If an iterative method is used, returns the # iterations in all Solve()
     // calls since the last call to Factorize(). A direct solver returns the #
     // iterative refinement steps.
-    Int iter() const;
+    Int iterSum() const;
+
+    // If an iterative method is used, returns the max # iterations in
+    // _all_ Solve() calls.
+    Int iterMax() const;
 
     // If a basis matrix is maintained, returns the # basis changes in the last
     // call to Factorize(). Otherwise returns 0.
     Int basis_changes() const;
 
+    //
+    double current_fill() const;
+
     // If a basis matrix is maintained, returns a pointer to it.
     // Otherwise returns NULL.
     const Basis* basis() const;
@@ -60,8 +67,10 @@ class KKTSolver {
     virtual void _Factorize(Iterate* iterate, Info* info) = 0;
     virtual void _Solve(const Vector& a, const Vector& b, double tol,
                          Vector& x, Vector& y, Info* info) = 0;
-    virtual Int _iter() const = 0;
+    virtual Int _iterSum() const = 0;
+    virtual Int _iterMax() const = 0;
     virtual Int _basis_changes() const { return 0; }
+    virtual double _current_fill() const { return 0.0; }
     virtual const Basis* _basis() const { return nullptr; }
 };