From 919774eae74283ce5986ad6a072b8c53de508aa8 Mon Sep 17 00:00:00 2001 From: JAJHall Date: Tue, 5 Nov 2024 01:52:49 +0000 Subject: [PATCH 1/6] Added trivial heuristics --- src/mip/HighsMipSolver.cpp | 4 + src/mip/HighsMipSolverData.cpp | 188 ++++++++++++++++++++++++++++++--- src/mip/HighsMipSolverData.h | 7 ++ src/mip/MipTimer.h | 16 ++- 4 files changed, 195 insertions(+), 20 deletions(-) diff --git a/src/mip/HighsMipSolver.cpp b/src/mip/HighsMipSolver.cpp index 4c46348b14..e339ee0863 100644 --- a/src/mip/HighsMipSolver.cpp +++ b/src/mip/HighsMipSolver.cpp @@ -175,6 +175,10 @@ void HighsMipSolver::run() { cleanupSolve(); return; } + // Apply the trivial heuristics + analysis_.mipTimerStart(kMipClockTrivialHeuristics); + mipdata_->trivialHeuristics(); + analysis_.mipTimerStop(kMipClockTrivialHeuristics); if (analysis_.analyse_mip_time & !submip) highsLogUser(options_mip_->log_options, HighsLogType::kInfo, "MIP-Timing: %11.2g - starting evaluate root node\n", diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp index 4cd73949ab..224ea05aa2 100644 --- a/src/mip/HighsMipSolverData.cpp +++ b/src/mip/HighsMipSolverData.cpp @@ -60,6 +60,18 @@ std::string HighsMipSolverData::solutionSourceToString( } else if (solution_source == kSolutionSourceUnbounded) { if (code) return "U"; return "Unbounded"; + } else if (solution_source == kSolutionSourceTrivialZ) { + if (code) return "z"; + return "Trivial zero"; + } else if (solution_source == kSolutionSourceTrivialL) { + if (code) return "l"; + return "Trivial lower"; + } else if (solution_source == kSolutionSourceTrivialU) { + if (code) return "u"; + return "Trivial upper"; + } else if (solution_source == kSolutionSourceTrivialP) { + if (code) return "p"; + return "Trivial point"; // } else if (solution_source == kSolutionSourceOpt1) { // if (code) return "1"; // return "1-opt"; @@ -136,6 +148,148 @@ bool HighsMipSolverData::trySolution(const std::vector& solution, return addIncumbent(solution, double(obj), solution_source); } +bool HighsMipSolverData::solutionRowFeasible( + const std::vector& solution) const { + for (HighsInt i = 0; i != mipsolver.model_->num_row_; ++i) { + HighsCDouble c_double_rowactivity = HighsCDouble(0.0); + + HighsInt start = ARstart_[i]; + HighsInt end = ARstart_[i + 1]; + + for (HighsInt j = start; j != end; ++j) + c_double_rowactivity += HighsCDouble(solution[ARindex_[j]] * ARvalue_[j]); + + double rowactivity = double(c_double_rowactivity); + if (rowactivity > mipsolver.rowUpper(i) + feastol) return false; + if (rowactivity < mipsolver.rowLower(i) - feastol) return false; + } + return true; +} + +void HighsMipSolverData::trivialHeuristics() { + // printf("\nHighsMipSolverData::trivialHeuristics() Number of continuous + // columns is %d\n", + // int(continuous_cols.size())); + if (continuous_cols.size() > 0) return; + const HighsInt num_try_heuristic = 4; + const std::vector heuristic_source = { + kSolutionSourceTrivialZ, kSolutionSourceTrivialL, kSolutionSourceTrivialU, + kSolutionSourceTrivialP}; + + const std::vector& col_lower = mipsolver.model_->col_lower_; + const std::vector& col_upper = mipsolver.model_->col_upper_; + const std::vector& row_lower = mipsolver.model_->row_lower_; + const std::vector& row_upper = mipsolver.model_->row_upper_; + const HighsSparseMatrix& matrix = mipsolver.model_->a_matrix_; + // Determine the following properties, according to which some + // trivial heuristics are duplicated or fail immediately + bool all_integer_lower_non_positive = true; + bool all_integer_lower_zero = true; + bool all_integer_lower_finite = true; + bool all_integer_upper_finite = true; + for (HighsInt integer_col = 0; integer_col < numintegercols; integer_col++) { + HighsInt iCol = integer_cols[integer_col]; + if (col_lower[iCol] > 0) all_integer_lower_non_positive = false; + if (col_lower[iCol]) all_integer_lower_zero = false; + if (col_lower[iCol] <= -kHighsInf) all_integer_lower_finite = false; + if (col_upper[iCol] >= kHighsInf) all_integer_upper_finite = false; + // Only continue if one of the properties still holds + if (!(all_integer_lower_non_positive || all_integer_lower_zero || + all_integer_upper_finite)) + break; + } + const bool all_integer_boxed = + all_integer_lower_finite && all_integer_upper_finite; + // printf( + // "Trying trivial heuristics\n" + // " all_integer_lower_non_positive = %d\n" + // " all_integer_lower_zero = %d\n" + // " all_integer_upper_finite = %d\n" + // " all_integer_boxed = %d\n", + // all_integer_lower_non_positive, all_integer_lower_zero, + // all_integer_upper_finite, all_integer_boxed); + const double feasibility_tolerance = + mipsolver.options_mip_->mip_feasibility_tolerance; + // Loop through the trivial heuristics + std::vector solution(mipsolver.model_->num_col_); + for (HighsInt try_heuristic = 0; try_heuristic < num_try_heuristic; + try_heuristic++) { + if (try_heuristic == 0) { + // First heuristic is to see whether all-zero for integer + // variables is feasible + // + // If there is a positive lower bound then the heuristic fails + if (!all_integer_lower_non_positive) continue; + // Determine whether a zero row activity is feasible + bool heuristic_failed = false; + for (HighsInt iRow = 0; iRow < mipsolver.model_->num_row_; iRow++) { + if (row_lower[iRow] > feasibility_tolerance || + row_upper[iRow] < -feasibility_tolerance) { + heuristic_failed = true; + break; + } + } + if (heuristic_failed) continue; + solution.assign(mipsolver.model_->num_col_, 0); + } else if (try_heuristic == 1) { + // Second heuristic is to see whether all-lower for integer + // variables (if distinct from all-zero) is feasible + if (all_integer_lower_zero) continue; + // Trivially feasible for columns + if (!solutionRowFeasible(col_lower)) continue; + solution = col_lower; + } else if (try_heuristic == 2) { + // Third heuristic is to see whether all-upper for integer + // variables is feasible + // + // If there is an infinite upper bound then the heuristic fails + if (!all_integer_upper_finite) continue; + // Trivially feasible for columns + if (!solutionRowFeasible(col_upper)) continue; + solution = col_upper; + } else if (try_heuristic == 3) { + // Fourth heuristic is to see whether the "lock point" is feasible + if (!all_integer_boxed) continue; + for (HighsInt integer_col = 0; integer_col < numintegercols; + integer_col++) { + HighsInt iCol = integer_cols[integer_col]; + HighsInt num_positive_values = 0; + HighsInt num_negative_values = 0; + for (HighsInt iEl = matrix.start_[iCol]; iEl < matrix.start_[iCol + 1]; + iEl++) { + if (matrix.value_[iEl] > 0) + num_positive_values++; + else + num_negative_values++; + } + solution[iCol] = num_positive_values > num_negative_values + ? col_lower[iCol] + : col_upper[iCol]; + } + // Trivially feasible for columns + if (!solutionRowFeasible(solution)) continue; + } + + HighsCDouble cdouble_obj = 0.0; + for (HighsInt iCol = 0; iCol < mipsolver.model_->num_col_; iCol++) + cdouble_obj += mipsolver.colCost(iCol) * solution[iCol]; + double obj = double(cdouble_obj); + const double save_upper_bound = upper_bound; + const bool new_incumbent = + addIncumbent(solution, obj, heuristic_source[try_heuristic]); + const bool lc_report = true; + if (lc_report) { + printf("Trivial heuristic %d has succeeded: objective = %g", + int(try_heuristic), obj); + if (new_incumbent) { + printf("; upper bound from %g to %g\n", save_upper_bound, upper_bound); + } else { + printf("\n"); + } + } + } +} + void HighsMipSolverData::startAnalyticCenterComputation( const highs::parallel::TaskGroup& taskGroup) { taskGroup.spawn([&]() { @@ -1312,9 +1466,12 @@ void HighsMipSolverData::printSolutionSourceKey() { // not a solution source, but used to force the last logging line to // be printed const int last_enum = kSolutionSourceCount - 1; - int half_list = last_enum / 2; + int third_list = 5; + int two_third_list = 10; + std::vector limits = {third_list, two_third_list, last_enum}; + ss.str(std::string()); - for (int k = 0; k < half_list; k++) { + for (int k = 0; k < limits[0]; k++) { if (k == 0) { ss << "\nSrc: "; } else { @@ -1323,21 +1480,22 @@ void HighsMipSolverData::printSolutionSourceKey() { ss << solutionSourceToString(k) << " => " << solutionSourceToString(k, false); } - highsLogUser(mipsolver.options_mip_->log_options, HighsLogType::kInfo, "%s\n", - ss.str().c_str()); - - ss.str(std::string()); - for (int k = half_list; k < last_enum; k++) { - if (k == half_list) { - ss << " "; - } else { - ss << "; "; + highsLogUser(mipsolver.options_mip_->log_options, HighsLogType::kInfo, + "%s;\n", ss.str().c_str()); + for (int line = 0; line < 2; line++) { + ss.str(std::string()); + for (int k = limits[line]; k < limits[line + 1]; k++) { + if (k == limits[line]) { + ss << " "; + } else { + ss << "; "; + } + ss << solutionSourceToString(k) << " => " + << solutionSourceToString(k, false); } - ss << solutionSourceToString(k) << " => " - << solutionSourceToString(k, false); + highsLogUser(mipsolver.options_mip_->log_options, HighsLogType::kInfo, + "%s%s\n", ss.str().c_str(), line == 0 ? ";" : ""); } - highsLogUser(mipsolver.options_mip_->log_options, HighsLogType::kInfo, "%s\n", - ss.str().c_str()); } void HighsMipSolverData::printDisplayLine(const int solution_source) { diff --git a/src/mip/HighsMipSolverData.h b/src/mip/HighsMipSolverData.h index c50ff5f1a2..6a9665f996 100644 --- a/src/mip/HighsMipSolverData.h +++ b/src/mip/HighsMipSolverData.h @@ -47,6 +47,10 @@ enum MipSolutionSource : int { kSolutionSourceSolveLp, kSolutionSourceEvaluateNode, kSolutionSourceUnbounded, + kSolutionSourceTrivialZ, + kSolutionSourceTrivialL, + kSolutionSourceTrivialU, + kSolutionSourceTrivialP, // kSolutionSourceOpt1, // kSolutionSourceOpt2, kSolutionSourceCleanup, @@ -199,6 +203,9 @@ struct HighsMipSolverData { domain.addConflictPool(conflictPool); } + bool solutionRowFeasible(const std::vector& solution) const; + void trivialHeuristics(); + void startAnalyticCenterComputation( const highs::parallel::TaskGroup& taskGroup); void finishAnalyticCenterComputation( diff --git a/src/mip/MipTimer.h b/src/mip/MipTimer.h index 6289a1aca0..bc6534cd93 100644 --- a/src/mip/MipTimer.h +++ b/src/mip/MipTimer.h @@ -24,6 +24,7 @@ enum iClockMip { kMipClockInit, kMipClockRunPresolve, kMipClockRunSetup, + kMipClockTrivialHeuristics, kMipClockEvaluateRootNode, kMipClockPerformAging0, kMipClockSearch, @@ -95,6 +96,8 @@ class MipTimer { clock[kMipClockInit] = timer_pointer->clock_def("Initialise"); clock[kMipClockRunPresolve] = timer_pointer->clock_def("Run presolve"); clock[kMipClockRunSetup] = timer_pointer->clock_def("Run setup"); + clock[kMipClockTrivialHeuristics] = + timer_pointer->clock_def("Trivial heuristics"); clock[kMipClockEvaluateRootNode] = timer_pointer->clock_def("Evaluate root node"); clock[kMipClockPerformAging0] = timer_pointer->clock_def("Perform aging 0"); @@ -249,11 +252,14 @@ class MipTimer { }; void reportMipLevel1Clock(const HighsTimerClock& mip_timer_clock) { - const std::vector mip_clock_list{ - kMipClockInit, kMipClockRunPresolve, - kMipClockRunSetup, kMipClockEvaluateRootNode, - kMipClockPerformAging0, kMipClockSearch, - kMipClockPostsolve}; + const std::vector mip_clock_list{kMipClockInit, + kMipClockRunPresolve, + kMipClockRunSetup, + kMipClockTrivialHeuristics, + kMipClockEvaluateRootNode, + kMipClockPerformAging0, + kMipClockSearch, + kMipClockPostsolve}; reportMipClockList("MipLevl1", mip_clock_list, mip_timer_clock, kMipClockTotal, tolerance_percent_report); }; From 5e6ec3e8951f80e1acd4f8e0b30e67b7ef3660db Mon Sep 17 00:00:00 2001 From: JAJHall Date: Tue, 5 Nov 2024 02:55:46 +0000 Subject: [PATCH 2/6] Corrected TEST_CASE("IP-infeasible-unbounded" now that kInfeasible and kUnbounded are distinguished when presolve is off --- check/TestMipSolver.cpp | 13 +++++++++++-- src/mip/HighsMipSolver.cpp | 9 ++++++++- src/mip/HighsMipSolverData.cpp | 13 ++++++++++--- src/mip/HighsMipSolverData.h | 2 +- 4 files changed, 30 insertions(+), 7 deletions(-) diff --git a/check/TestMipSolver.cpp b/check/TestMipSolver.cpp index eeba016a62..0249725857 100644 --- a/check/TestMipSolver.cpp +++ b/check/TestMipSolver.cpp @@ -677,8 +677,17 @@ TEST_CASE("IP-infeasible-unbounded", "[highs_test_mip_solver]") { // Solve highs.passModel(lp); highs.run(); - REQUIRE(highs.getModelStatus() == - HighsModelStatus::kUnboundedOrInfeasible); + HighsModelStatus required_model_status; + if (k == 0) { + if (l == 0) { + required_model_status = HighsModelStatus::kInfeasible; + } else { + required_model_status = HighsModelStatus::kUnbounded; + } + } else { + required_model_status = HighsModelStatus::kUnboundedOrInfeasible; + } + REQUIRE(highs.getModelStatus() == required_model_status); } highs.setOptionValue("presolve", kHighsOnString); } diff --git a/src/mip/HighsMipSolver.cpp b/src/mip/HighsMipSolver.cpp index e339ee0863..8482442e35 100644 --- a/src/mip/HighsMipSolver.cpp +++ b/src/mip/HighsMipSolver.cpp @@ -177,7 +177,14 @@ void HighsMipSolver::run() { } // Apply the trivial heuristics analysis_.mipTimerStart(kMipClockTrivialHeuristics); - mipdata_->trivialHeuristics(); + HighsModelStatus model_status = mipdata_->trivialHeuristics(); + if (modelstatus_ == HighsModelStatus::kNotset && + model_status == HighsModelStatus::kInfeasible) { + // trivialHeuristics can spot trivial infeasibility, so act on it + modelstatus_ = model_status; + cleanupSolve(); + return; + } analysis_.mipTimerStop(kMipClockTrivialHeuristics); if (analysis_.analyse_mip_time & !submip) highsLogUser(options_mip_->log_options, HighsLogType::kInfo, diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp index 224ea05aa2..aadcf4ab5f 100644 --- a/src/mip/HighsMipSolverData.cpp +++ b/src/mip/HighsMipSolverData.cpp @@ -166,7 +166,7 @@ bool HighsMipSolverData::solutionRowFeasible( return true; } -void HighsMipSolverData::trivialHeuristics() { +HighsModelStatus HighsMipSolverData::trivialHeuristics() { // printf("\nHighsMipSolverData::trivialHeuristics() Number of continuous // columns is %d\n", // int(continuous_cols.size())); @@ -176,8 +176,8 @@ void HighsMipSolverData::trivialHeuristics() { kSolutionSourceTrivialZ, kSolutionSourceTrivialL, kSolutionSourceTrivialU, kSolutionSourceTrivialP}; - const std::vector& col_lower = mipsolver.model_->col_lower_; - const std::vector& col_upper = mipsolver.model_->col_upper_; + std::vector col_lower = mipsolver.model_->col_lower_; + std::vector col_upper = mipsolver.model_->col_upper_; const std::vector& row_lower = mipsolver.model_->row_lower_; const std::vector& row_upper = mipsolver.model_->row_upper_; const HighsSparseMatrix& matrix = mipsolver.model_->a_matrix_; @@ -189,6 +189,12 @@ void HighsMipSolverData::trivialHeuristics() { bool all_integer_upper_finite = true; for (HighsInt integer_col = 0; integer_col < numintegercols; integer_col++) { HighsInt iCol = integer_cols[integer_col]; + // Round bounds in to nearest integer + col_lower[iCol] = std::ceil(col_lower[iCol]); + col_upper[iCol] = std::floor(col_upper[iCol]); + // If bounds are inconsistent then MIP is infeasible + if (col_lower[iCol] > col_upper[iCol]) return HighsModelStatus::kInfeasible; + if (col_lower[iCol] > 0) all_integer_lower_non_positive = false; if (col_lower[iCol]) all_integer_lower_zero = false; if (col_lower[iCol] <= -kHighsInf) all_integer_lower_finite = false; @@ -288,6 +294,7 @@ void HighsMipSolverData::trivialHeuristics() { } } } + return HighsModelStatus::kNotset; } void HighsMipSolverData::startAnalyticCenterComputation( diff --git a/src/mip/HighsMipSolverData.h b/src/mip/HighsMipSolverData.h index 6a9665f996..bf7522b086 100644 --- a/src/mip/HighsMipSolverData.h +++ b/src/mip/HighsMipSolverData.h @@ -204,7 +204,7 @@ struct HighsMipSolverData { } bool solutionRowFeasible(const std::vector& solution) const; - void trivialHeuristics(); + HighsModelStatus trivialHeuristics(); void startAnalyticCenterComputation( const highs::parallel::TaskGroup& taskGroup); From a2ac9583d6057e09ea5c42beada1a3f311ad616e Mon Sep 17 00:00:00 2001 From: JAJHall Date: Tue, 5 Nov 2024 09:36:24 +0000 Subject: [PATCH 3/6] Added return value to HighsMipSolverData::trivialHeuristics() --- src/mip/HighsMipSolverData.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp index aadcf4ab5f..07b599254c 100644 --- a/src/mip/HighsMipSolverData.cpp +++ b/src/mip/HighsMipSolverData.cpp @@ -170,7 +170,7 @@ HighsModelStatus HighsMipSolverData::trivialHeuristics() { // printf("\nHighsMipSolverData::trivialHeuristics() Number of continuous // columns is %d\n", // int(continuous_cols.size())); - if (continuous_cols.size() > 0) return; + if (continuous_cols.size() > 0) return HighsModelStatus::kOk; const HighsInt num_try_heuristic = 4; const std::vector heuristic_source = { kSolutionSourceTrivialZ, kSolutionSourceTrivialL, kSolutionSourceTrivialU, From 21ba8d41f6f5815ec364ad39ab2159baffcadbbd Mon Sep 17 00:00:00 2001 From: JAJHall Date: Tue, 5 Nov 2024 09:40:30 +0000 Subject: [PATCH 4/6] Need to return a HighsModelStatus! --- src/mip/HighsMipSolverData.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp index 07b599254c..7762d7495c 100644 --- a/src/mip/HighsMipSolverData.cpp +++ b/src/mip/HighsMipSolverData.cpp @@ -170,7 +170,7 @@ HighsModelStatus HighsMipSolverData::trivialHeuristics() { // printf("\nHighsMipSolverData::trivialHeuristics() Number of continuous // columns is %d\n", // int(continuous_cols.size())); - if (continuous_cols.size() > 0) return HighsModelStatus::kOk; + if (continuous_cols.size() > 0) return HighsModelStatus::kNotset; const HighsInt num_try_heuristic = 4; const std::vector heuristic_source = { kSolutionSourceTrivialZ, kSolutionSourceTrivialL, kSolutionSourceTrivialU, From cc9975f88326a5f1d603e22df6a1313e585084d9 Mon Sep 17 00:00:00 2001 From: JAJHall Date: Tue, 5 Nov 2024 09:52:00 +0000 Subject: [PATCH 5/6] Need to sort out getColsIntegrality --- src/highs_bindings.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/highs_bindings.cpp b/src/highs_bindings.cpp index 575968b997..08f13fd341 100644 --- a/src/highs_bindings.cpp +++ b/src/highs_bindings.cpp @@ -1042,6 +1042,9 @@ PYBIND11_MODULE(_core, m) { .def("getRowName", &highs_getRowName) .def("getRowByName", &highs_getRowByName) + // .def("getColIntegrality", &Highs::getColIntegrality) + // .def("getColsIntegrality", &Highs::getColsIntegrality) + .def("writeModel", &Highs::writeModel) .def("writePresolvedModel", &Highs::writePresolvedModel) .def("crossover", &Highs::crossover) From fb7b1c7c608fd6b5801aee880772499a3fd8c182 Mon Sep 17 00:00:00 2001 From: Julian Hall Date: Tue, 5 Nov 2024 14:55:20 +0000 Subject: [PATCH 6/6] Added getColIntegrality to highspy --- src/highs_bindings.cpp | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/src/highs_bindings.cpp b/src/highs_bindings.cpp index 08f13fd341..f2543ba1f5 100644 --- a/src/highs_bindings.cpp +++ b/src/highs_bindings.cpp @@ -537,6 +537,14 @@ highs_getColsEntries(Highs* h, HighsInt num_set_entries, py::cast(value)); } +std::tuple +highs_getColIntegrality(Highs* h, HighsInt col) { + HighsInt col_ = static_cast(col); + HighsVarType integrality; + HighsStatus status = h->getColIntegrality(col_, integrality); + return std::make_tuple(status, integrality); +} + std::tuple, dense_array_t, HighsInt> highs_getRows(Highs* h, HighsInt num_set_entries, @@ -1029,11 +1037,13 @@ PYBIND11_MODULE(_core, m) { .def("getCol", &highs_getCol) .def("getColEntries", &highs_getColEntries) + .def("getColIntegrality", &highs_getColIntegrality) .def("getRow", &highs_getRow) .def("getRowEntries", &highs_getRowEntries) .def("getCols", &highs_getCols) .def("getColsEntries", &highs_getColsEntries) + .def("getRows", &highs_getRows) .def("getRowsEntries", &highs_getRowsEntries) @@ -1042,9 +1052,6 @@ PYBIND11_MODULE(_core, m) { .def("getRowName", &highs_getRowName) .def("getRowByName", &highs_getRowByName) - // .def("getColIntegrality", &Highs::getColIntegrality) - // .def("getColsIntegrality", &Highs::getColsIntegrality) - .def("writeModel", &Highs::writeModel) .def("writePresolvedModel", &Highs::writePresolvedModel) .def("crossover", &Highs::crossover)