WIP: scipy callback concept for IPX

CMakeLists.txt

@@ -7,11 +7,11 @@ cmake_minimum_required(VERSION 3.15)
 # set preference for clang compiler and intel compiler over gcc and other compilers
 include(Platform/${CMAKE_SYSTEM_NAME}-Determine-C OPTIONAL)
 include(Platform/${CMAKE_SYSTEM_NAME}-C OPTIONAL)
-set(CMAKE_C_COMPILER_NAMES clang icc cc ${CMAKE_C_COMPILER_NAMES})
+set(CMAKE_C_COMPILER_NAMES cc clang icc ${CMAKE_C_COMPILER_NAMES})
 
 include(Platform/${CMAKE_SYSTEM_NAME}-Determine-CXX OPTIONAL)
 include(Platform/${CMAKE_SYSTEM_NAME}-CXX OPTIONAL)
-set(CMAKE_CXX_COMPILER_NAMES clang++ icpc c++ ${CMAKE_CXX_COMPILER_NAMES})
+set(CMAKE_CXX_COMPILER_NAMES c++ clang++ icpc ${CMAKE_CXX_COMPILER_NAMES})
 
 # set default build type before project call, as it otherwise seems to fail for some plattforms
 if(NOT CMAKE_BUILD_TYPE)

src/Highs.h

@@ -18,6 +18,7 @@
 
 #include <sstream>
 
+#include "HighsClbk.h"
 #include "lp_data/HighsLpUtils.h"
 #include "lp_data/HighsRanging.h"
 #include "lp_data/HighsSolutionDebug.h"
@@ -1007,6 +1008,8 @@ class Highs {
   HighsStatus setHighsOptionValue(const std::string& option,
                                   const HighsInt value);
 
+  void setClbk(HighsClbk clbk_fun) { clbk_fun_ = clbk_fun; }
+
 #ifdef HIGHSINT64
   HighsStatus setHighsOptionValue(const std::string& option,
                                   const int value  //!< The option value
@@ -1115,6 +1118,8 @@ class Highs {
 
   bool written_log_header = false;
 
+  HighsClbk clbk_fun_ = nullptr;
+
   void exactResizeModel() {
     this->model_.lp_.exactResize();
     this->model_.hessian_.exactResize();

src/ipm/IpxWrapper.cpp

@@ -18,6 +18,7 @@
 
 #include <cassert>
 
+#include "HighsClbk.h"
 #include "lp_data/HighsOptions.h"
 #include "lp_data/HighsSolution.h"
 
@@ -26,7 +27,8 @@ 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.model_status_, solver_object.highs_info_,
+                    solver_object.clbk_fun_);
 }
 
 HighsStatus solveLpIpx(const HighsOptions& options,
@@ -35,7 +37,8 @@ HighsStatus solveLpIpx(const HighsOptions& options,
                        HighsBasis& highs_basis,
 		       HighsSolution& highs_solution,
                        HighsModelStatus& model_status,
-                       HighsInfo& highs_info) {
+                       HighsInfo& highs_info,
+                       HighsClbk clbk_fun) {
   // Use IPX to try to solve the LP
   //
   // Can return HighsModelStatus (HighsStatus) values:
@@ -90,6 +93,7 @@ HighsStatus solveLpIpx(const HighsOptions& options,
   } else if (options.log_dev_level == kHighsLogDevLevelVerbose) {
     parameters.debug = 4;
   }
+  parameters.clbk_fun = clbk_fun;
   // Just test feasibility and optimality tolerances for now
   // ToDo Set more parameters
   parameters.ipm_feasibility_tol = min(options.primal_feasibility_tolerance,

src/ipm/IpxWrapper.h

@@ -29,7 +29,8 @@ HighsStatus solveLpIpx(HighsLpSolverObject& solver_object);
 HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer,
                        const HighsLp& lp, HighsBasis& highs_basis,
                        HighsSolution& highs_solution,
-                       HighsModelStatus& model_status, HighsInfo& highs_info);
+                       HighsModelStatus& model_status, HighsInfo& highs_info,
+                       HighsClbk clbk_t);
 
 void fillInIpxData(const HighsLp& lp, ipx::Int& num_col, ipx::Int& num_row,
                    std::vector<double>& obj, std::vector<double>& col_lb,

src/ipm/ipx/ipm.cc

@@ -39,7 +39,8 @@ void IPM::StartingPoint(KKTSolver* kkt, Iterate* iterate, Info* info) {
     }
 }
 
-void IPM::Driver(KKTSolver* kkt, Iterate* iterate, Info* info) {
+void IPM::Driver(KKTSolver* kkt, Iterate* iterate, Info* info,
+                 HighsClbk clbk_fun) {
     const Model& model = iterate->model();
     const Int m = model.rows();
     const Int n = model.cols();
@@ -75,6 +76,12 @@ void IPM::Driver(KKTSolver* kkt, Iterate* iterate, Info* info) {
             break;
         MakeStep(step);
         info->iter++;
+        if (clbk_fun) {
+            HighsClbkInfo clbk_info {
+                info->iter
+            };
+            clbk_fun(&clbk_info);
+        }
         PrintOutput();
     }
 

src/ipm/ipx/ipm.h

@@ -1,9 +1,10 @@
 #ifndef IPX_IPM_H_
 #define IPX_IPM_H_
 
+#include "HighsClbk.h"
 #include "ipm/ipx/control.h"
-#include "ipm/ipx/kkt_solver.h"
 #include "ipm/ipx/iterate.h"
+#include "ipm/ipx/kkt_solver.h"
 
 namespace ipx {
 
@@ -33,7 +34,8 @@ class IPM {
     // IPX_STATUS_no_progress   if no progress over a number of iterations,
     // IPX_STATUS_time_limit    if interrupted by time limit,
     // IPX_STATUS_failed        if the KKT solver failed with info->errflag.
-    void Driver(KKTSolver* kkt, Iterate* iterate, Info* info);
+    void Driver(KKTSolver* kkt, Iterate* iterate, Info* info,
+                HighsClbk clbk_fun);
 
     Int maxiter() const { return maxiter_; }
     void maxiter(Int i) { maxiter_ = i; }

src/ipm/ipx/ipx_internal.h

@@ -4,6 +4,7 @@
 #include <cstring>
 #include <valarray>
 
+#include "HighsClbk.h"
 #include "ipm/ipx/ipx_config.h"
 #include "ipm/ipx/ipx_info.h"
 #include "ipm/ipx/ipx_parameters.h"
@@ -15,6 +16,7 @@ using Int = ipxint;
 
 struct Info : public ipx_info {
   Info() { std::memset(this, 0, sizeof(Info)); }
+
 };
 
 struct Parameters : public ipx_parameters {
@@ -23,6 +25,7 @@ struct Parameters : public ipx_parameters {
     logfile = nullptr;
     print_interval = 5.0;
     analyse_basis_data = false;
+    HighsClbk clbk_fun = nullptr;
     time_limit = -1.0;
     dualize = -1;
     scale = 1;

src/ipm/ipx/ipx_parameters.h

@@ -1,6 +1,7 @@
 #ifndef IPX_PARAMETERS_H_
 #define IPX_PARAMETERS_H_
 
+#include "HighsClbk.h"
 #include "ipm/ipx/ipx_config.h"
 
 #ifdef __cplusplus
@@ -13,6 +14,7 @@ struct ipx_parameters {
     double print_interval;
     double time_limit;
     bool analyse_basis_data;
+    HighsClbk clbk_fun;
 
     /* Preprocessing */
     ipxint dualize;

src/ipm/ipx/lp_solver.cc

@@ -474,7 +474,7 @@ void LpSolver::RunInitialIPM(IPM& ipm) {
     } else {
         ipm.maxiter(std::min(switchiter, control_.ipm_maxiter()));
     }
-    ipm.Driver(&kkt, iterate_.get(), &info_);
+    ipm.Driver(&kkt, iterate_.get(), &info_, GetParameters().clbk_fun);
     switch (info_.status_ipm) {
     case IPX_STATUS_optimal:
         // If the IPM reached its termination criterion in the initial
@@ -536,7 +536,7 @@ void LpSolver::RunMainIPM(IPM& ipm) {
     KKTSolverBasis kkt(control_, *basis_);
     Timer timer;
     ipm.maxiter(control_.ipm_maxiter());
-    ipm.Driver(&kkt, iterate_.get(), &info_);
+    ipm.Driver(&kkt, iterate_.get(), &info_, GetParameters().clbk_fun);
     info_.time_ipm2 = timer.Elapsed();
 }
 

src/lp_data/Highs.cpp

@@ -1428,7 +1428,8 @@ HighsStatus Highs::getRanging() {
   // Create a HighsLpSolverObject of references to data in the Highs
   // class, and the scaled/unscaled model status
   HighsLpSolverObject solver_object(model_.lp_, basis_, solution_, info_,
-                                    ekk_instance_, options_, timer_);
+                                    ekk_instance_, options_, timer_,
+                                    clbk_fun_);
   solver_object.model_status_ = model_status_;
   return getRangingData(this->ranging_, solver_object);
 }
@@ -1730,7 +1731,8 @@ HighsStatus Highs::setBasis(const HighsBasis& basis,
     HighsBasis modifiable_basis = basis;
     modifiable_basis.was_alien = true;
     HighsLpSolverObject solver_object(model_.lp_, modifiable_basis, solution_,
-                                      info_, ekk_instance_, options_, timer_);
+                                      info_, ekk_instance_, options_, timer_,
+                                      clbk_fun_);
     HighsStatus return_status = formSimplexLpBasisAndFactor(solver_object);
     if (return_status != HighsStatus::kOk) return HighsStatus::kError;
     // Update the HiGHS basis
@@ -2752,7 +2754,7 @@ HighsStatus Highs::callSolveLp(HighsLp& lp, const string message) {
   // Create a HighsLpSolverObject of references to data in the Highs
   // class, and the scaled/unscaled model status
   HighsLpSolverObject solver_object(lp, basis_, solution_, info_, ekk_instance_,
-                                    options_, timer_);
+                                    options_, timer_, clbk_fun_);
 
   // Check that the model is column-wise
   assert(model_.lp_.a_matrix_.isColwise());

src/lp_data/HighsInterface.cpp

@@ -1093,7 +1093,8 @@ HighsStatus Highs::getBasicVariablesInterface(HighsInt* basic_variables) {
     //
     // Create a HighsLpSolverObject
     HighsLpSolverObject solver_object(lp, basis_, solution_, info_,
-                                      ekk_instance_, options_, timer_);
+                                      ekk_instance_, options_, timer_,
+                                      clbk_fun_);
     const bool only_from_known_basis = true;
     return_status = interpretCallStatus(
         options_.log_options,

src/lp_data/HighsLpSolverObject.h

@@ -16,6 +16,7 @@
 #ifndef LP_DATA_HIGHS_LP_SOLVER_OBJECT_H_
 #define LP_DATA_HIGHS_LP_SOLVER_OBJECT_H_
 
+#include "HighsClbk.h"
 #include "lp_data/HighsInfo.h"
 #include "lp_data/HighsOptions.h"
 #include "simplex/HEkk.h"
@@ -24,14 +25,16 @@ class HighsLpSolverObject {
  public:
   HighsLpSolverObject(HighsLp& lp, HighsBasis& basis, HighsSolution& solution,
                       HighsInfo& highs_info, HEkk& ekk_instance,
-                      HighsOptions& options, HighsTimer& timer)
+                      HighsOptions& options, HighsTimer& timer,
+                      HighsClbk clbk_fun)
       : lp_(lp),
         basis_(basis),
         solution_(solution),
         highs_info_(highs_info),
         ekk_instance_(ekk_instance),
         options_(options),
-        timer_(timer) {}
+        timer_(timer),
+        clbk_fun_(clbk_fun) {}
 
   HighsLp& lp_;
   HighsBasis& basis_;
@@ -40,6 +43,7 @@ class HighsLpSolverObject {
   HEkk& ekk_instance_;
   HighsOptions& options_;
   HighsTimer& timer_;
+  HighsClbk clbk_fun_;
 
   HighsModelStatus model_status_ = HighsModelStatus::kNotset;
 };

src/util/HighsClbk.h

@@ -0,0 +1,13 @@
+#ifndef HIGHS_HIGHSCLBK_H
+#define HIGHS_HIGHSCLBK_H
+
+#include "HighsInt.h"
+
+struct HighsClbkInfo {
+    explicit HighsClbkInfo(const HighsInt& iteration)
+      : iteration_(iteration) {}
+    const HighsInt& iteration_;
+};
+typedef void(*HighsClbk)(HighsClbkInfo*);
+
+#endif  // HIGHS_HIGHSCLBK_H