Merge pull request #1884 from few/highspy-feasibilityRelaxation

highspy: Expose new feasibility relaxation feature

examples/call_highs_from_python.py

@@ -339,3 +339,158 @@ def user_interrupt_callback(
 for icol in range(num_var-2, num_var):
     print(icol, solution.col_value[icol],
           h.basisStatusToString(basis.col_status[icol]))
+
+print("computing IIS for lp-incompatible-bounds")
+"""
+LP has row0 and col2 with inconsistent bounds.
+
+When prioritising rows, row0 and its constituent columns (1, 2) should be found
+When prioritising columns, col2 and its constituent rows (0, 1) should be found
+"""
+# Define the LP
+lp = highspy.HighsLp()
+lp.num_col_ = 3
+lp.num_row_ = 2
+lp.col_cost_ = np.array([0, 0, 0], dtype=np.double)
+lp.col_lower_ = np.array([0, 0, 0], dtype=np.double)
+lp.col_upper_ = np.array([1, 1, -1], dtype=np.double)
+lp.row_lower_ = np.array([1, 0], dtype=np.double)
+lp.row_upper_ = np.array([0, 1], dtype=np.double)
+lp.a_matrix_.format_ = highspy.MatrixFormat.kRowwise
+lp.a_matrix_.start_ = np.array([0, 2, 4])
+lp.a_matrix_.index_ = np.array([1, 2, 0, 2])
+lp.a_matrix_.value_ = np.array([1, 1, 1, 1], dtype=np.double)
+
+h.clear()
+h.passModel(lp)
+h.run()
+assert h.getModelStatus() == highspy.HighsModelStatus.kInfeasible
+
+# Set IIS strategy to row priority and get IIS
+h.setOptionValue("iis_strategy", highspy.IisStrategy.kIisStrategyFromLpRowPriority)
+
+iis = highspy.HighsIis()
+assert h.getIis(iis) == highspy.HighsStatus.kOk
+assert len(iis.col_index) == 0
+assert len(iis.row_index) == 1
+assert iis.row_index[0] == 0
+assert iis.row_bound[0] == highspy.IisBoundStatus.kIisBoundStatusBoxed
+
+# Set IIS strategy to column priority and get IIS
+h.setOptionValue("iis_strategy", highspy.IisStrategy.kIisStrategyFromLpColPriority)
+iis.invalidate()
+assert h.getIis(iis) == highspy.HighsStatus.kOk
+assert len(iis.col_index) == 1
+assert len(iis.row_index) == 0
+assert iis.col_index[0] == 2
+assert iis.col_bound[0] == highspy.IisBoundStatus.kIisBoundStatusBoxed
+
+print("IIS computation completed successfully")
+
+print("computing feasibility relaxation")
+h.clear()
+inf = h.getInfinity()
+
+num_col = 2
+num_row = 3
+num_nz = 6
+a_format = highspy.MatrixFormat.kColwise
+sense = highspy.ObjSense.kMinimize
+offset = 0
+col_cost = np.array([1, -2], dtype=np.double)
+col_lower = np.array([5, -inf], dtype=np.double)
+col_upper = np.array([inf, inf], dtype=np.double)
+row_lower = np.array([2, -inf, -inf], dtype=np.double)
+row_upper = np.array([inf, 1, 20], dtype=np.double)
+a_start = np.array([0, 3])
+a_index = np.array([0, 1, 2, 0, 1, 2])
+a_value = np.array([-1, -3, 20, 21, 2, 1], dtype=np.double)
+integrality = np.array([highspy.HighsVarType.kInteger, highspy.HighsVarType.kInteger])
+
+h.passModel(
+    num_col, num_row, num_nz, a_format, sense, offset,
+    col_cost, col_lower, col_upper,
+    row_lower, row_upper,
+    a_start, a_index, a_value,
+    integrality
+)
+
+assert h.feasibilityRelaxation(1, 1, 1) == highspy.HighsStatus.kOk
+
+info = h.getInfo()
+objective_function_value = info.objective_function_value
+
+solution = h.getSolution()
+
+assert abs(objective_function_value - 5) < 1e-6, f"Expected objective value 5, got {objective_function_value}"
+assert abs(solution.col_value[0] - 1) < 1e-6, f"Expected solution[0] = 1, got {solution.col_value[0]}"
+assert abs(solution.col_value[1] - 1) < 1e-6, f"Expected solution[1] = 1, got {solution.col_value[1]}"
+
+print("Feasibility Relaxation Test Passed")
+
+# Using infeasible LP from AMPL documentation
+h = highspy.Highs()
+lp = highspy.HighsLp()
+lp.num_col_ = 2
+lp.num_row_ = 3
+lp.col_cost_ = np.array([1, -2], dtype=np.double)
+lp.col_lower_ = np.array([5, -h.getInfinity()], dtype=np.double)
+lp.col_upper_ = np.array([h.getInfinity(), h.getInfinity()], dtype=np.double)
+lp.col_names_ = ["X", "Y"]
+lp.row_lower_ = np.array([2, -h.getInfinity(), -h.getInfinity()], dtype=np.double)
+lp.row_upper_ = np.array([h.getInfinity(), 1, 20], dtype=np.double)
+lp.row_names_ = ["R0", "R1", "R2"]
+lp.a_matrix_.start_ = np.array([0, 3, 6])
+lp.a_matrix_.index_ = np.array([0, 1, 2, 0, 1, 2])
+lp.a_matrix_.value_ = np.array([-1, -3, 20, 21, 2, 1], dtype=np.double)
+lp.integrality_ = np.array([highspy.HighsVarType.kInteger, highspy.HighsVarType.kInteger])
+h.setOptionValue("output_flag", False)
+h.passModel(lp)
+
+# Vanilla feasibility relaxation
+print("Vanilla feasibility relaxation")
+h.feasibilityRelaxation(1, 1, 1)
+solution = h.getSolution()
+print(f"Solution: ({solution.col_value[0]}, {solution.col_value[1]})")
+print(f"Slacks: ({solution.row_value[0] - lp.row_lower_[0]}, "
+      f"{lp.row_upper_[1] - solution.row_value[1]}, "
+      f"{lp.row_upper_[2] - solution.row_value[2]})")
+
+# Respect all lower bounds
+print("\nRespect all lower bounds")
+h.feasibilityRelaxation(-1, 1, 1)
+solution = h.getSolution()
+print(f"Solution: ({solution.col_value[0]}, {solution.col_value[1]})")
+print(f"Slacks: ({solution.row_value[0] - lp.row_lower_[0]}, "
+      f"{lp.row_upper_[1] - solution.row_value[1]}, "
+      f"{lp.row_upper_[2] - solution.row_value[2]})")
+
+# Local penalties RHS {1, -1, 10}
+print("\nLocal penalties RHS {1, -1, 10}")
+local_rhs_penalty = np.array([1, -1, 10], dtype=np.double)
+h.feasibilityRelaxation(1, 1, 0, None, None, local_rhs_penalty)
+solution = h.getSolution()
+print(f"Solution: ({solution.col_value[0]}, {solution.col_value[1]})")
+print(f"Slacks: ({solution.row_value[0] - lp.row_lower_[0]}, "
+      f"{lp.row_upper_[1] - solution.row_value[1]}, "
+      f"{lp.row_upper_[2] - solution.row_value[2]})")
+
+# Local penalties RHS {10, 1, 1}
+print("\nLocal penalties RHS {10, 1, 1}")
+local_rhs_penalty = np.array([10, 1, 1], dtype=np.double)
+h.feasibilityRelaxation(1, 1, 0, None, None, local_rhs_penalty)
+solution = h.getSolution()
+print(f"Solution: ({solution.col_value[0]}, {solution.col_value[1]})")
+print(f"Slacks: ({solution.row_value[0] - lp.row_lower_[0]}, "
+      f"{lp.row_upper_[1] - solution.row_value[1]}, "
+      f"{lp.row_upper_[2] - solution.row_value[2]})")
+
+iis = highspy.HighsIis()
+assert h.getIis(iis) == highspy.HighsStatus.kOk
+
+print("\nIIS")
+print("row_index:", iis.row_index)
+print("row_bound:", iis.row_bound)
+
+print("col_index:", iis.col_index)
+print("col_bound:", iis.col_bound)

src/highs_bindings.cpp

@@ -671,6 +671,20 @@ PYBIND11_MODULE(_core, m) {
       .value("kWarning", HighsLogType::kWarning)
       .value("kError", HighsLogType::kError);
       // .export_values();
+  py::enum_<IisStrategy>(m, "IisStrategy")
+      .value("kIisStrategyMin", IisStrategy::kIisStrategyMin)
+      .value("kIisStrategyFromLpRowPriority", IisStrategy::kIisStrategyFromLpRowPriority)
+      .value("kIisStrategyFromLpColPriority", IisStrategy::kIisStrategyFromLpColPriority)
+      .value("kIisStrategyMax", IisStrategy::kIisStrategyMax);
+      // .export_values();
+  py::enum_<IisBoundStatus>(m, "IisBoundStatus")
+      .value("kIisBoundStatusDropped", IisBoundStatus::kIisBoundStatusDropped)
+      .value("kIisBoundStatusNull", IisBoundStatus::kIisBoundStatusNull)
+      .value("kIisBoundStatusFree", IisBoundStatus::kIisBoundStatusFree)
+      .value("kIisBoundStatusLower", IisBoundStatus::kIisBoundStatusLower)
+      .value("kIisBoundStatusUpper", IisBoundStatus::kIisBoundStatusUpper)
+      .value("kIisBoundStatusBoxed", IisBoundStatus::kIisBoundStatusBoxed);
+      // .export_values();
   // Classes
   py::class_<HighsSparseMatrix>(m, "HighsSparseMatrix")
       .def(py::init<>())
@@ -861,6 +875,37 @@ PYBIND11_MODULE(_core, m) {
       .def("postsolve", &highs_postsolve)
       .def("postsolve", &highs_mipPostsolve)
       .def("run", &Highs::run)
+      .def("feasibilityRelaxation", 
+     [](Highs& self, double global_lower_penalty, double global_upper_penalty, double global_rhs_penalty,
+        py::object local_lower_penalty, py::object local_upper_penalty, py::object local_rhs_penalty) {
+         std::vector<double> llp, lup, lrp;
+         const double* llp_ptr = nullptr;
+         const double* lup_ptr = nullptr;
+         const double* lrp_ptr = nullptr;
+
+         if (!local_lower_penalty.is_none()) {
+             llp = local_lower_penalty.cast<std::vector<double>>();
+             llp_ptr = llp.data();
+         }
+         if (!local_upper_penalty.is_none()) {
+             lup = local_upper_penalty.cast<std::vector<double>>();
+             lup_ptr = lup.data();
+         }
+         if (!local_rhs_penalty.is_none()) {
+             lrp = local_rhs_penalty.cast<std::vector<double>>();
+             lrp_ptr = lrp.data();
+         }
+
+         return self.feasibilityRelaxation(global_lower_penalty, global_upper_penalty, global_rhs_penalty,
+                                           llp_ptr, lup_ptr, lrp_ptr);
+     },
+     py::arg("global_lower_penalty"),
+     py::arg("global_upper_penalty"),
+     py::arg("global_rhs_penalty"),
+     py::arg("local_lower_penalty") = py::none(),
+     py::arg("local_upper_penalty") = py::none(),
+     py::arg("local_rhs_penalty") = py::none())
+      .def("getIis", &Highs::getIis)
       .def("presolve", &Highs::presolve)
       .def("writeSolution", &highs_writeSolution)
       .def("readSolution", &Highs::readSolution)
@@ -975,6 +1020,17 @@ PYBIND11_MODULE(_core, m) {
                                    &Highs::startCallback))
       .def("stopCallbackInt", static_cast<HighsStatus (Highs::*)(const int)>(
                                   &Highs::stopCallback));
+
+  py::class_<HighsIis>(m, "HighsIis")
+    .def(py::init<>())
+    .def("invalidate", &HighsIis::invalidate)
+    .def_readwrite("valid", &HighsIis::valid_)
+    .def_readwrite("strategy", &HighsIis::strategy_)
+    .def_readwrite("col_index", &HighsIis::col_index_)
+    .def_readwrite("row_index", &HighsIis::row_index_)
+    .def_readwrite("col_bound", &HighsIis::col_bound_)
+    .def_readwrite("row_bound", &HighsIis::row_bound_)
+    .def_readwrite("info", &HighsIis::info_);
   // structs
   py::class_<HighsSolution>(m, "HighsSolution")
       .def(py::init<>())
@@ -1013,6 +1069,10 @@ PYBIND11_MODULE(_core, m) {
       .def_readwrite("col_bound_dn", &HighsRanging::col_bound_dn)
       .def_readwrite("row_bound_up", &HighsRanging::row_bound_up)
       .def_readwrite("row_bound_dn", &HighsRanging::row_bound_dn);
+  py::class_<HighsIisInfo>(m, "HighsIisInfo")
+    .def(py::init<>())
+    .def_readwrite("simplex_time", &HighsIisInfo::simplex_time)
+    .def_readwrite("simplex_iterations", &HighsIisInfo::simplex_iterations);
   // constants
   m.attr("kHighsInf") = kHighsInf;
   m.attr("kHighsIInf") = kHighsIInf;

src/highspy/__init__.py

@@ -14,13 +14,16 @@
     HighsInfoType, \
     HighsStatus, \
     HighsLogType, \
+    IisStrategy, \
+    IisBoundStatus, \
     HighsSparseMatrix, \
     HighsLp, \
     HighsHessian, \
     HighsModel, \
     HighsInfo, \
     HighsOptions, \
     _Highs, \
+    HighsIis, \
     HighsSolution, \
     HighsObjectiveSolution, \
     HighsBasis, \
@@ -111,6 +114,8 @@
            "HighsInfoType",
            "HighsStatus",
            "HighsLogType",
+           "IisStrategy",
+           "IisBoundStatus",
            "HighsSparseMatrix",
            "HighsLp",
            "HighsHessian",
@@ -119,6 +124,7 @@
            "HighsOptions",
            "_Highs",
            "Highs",
+           "HighsIis",
            "HighsSolution",
            "HighsObjectiveSolution",
            "HighsBasis",