update gscip code

This commit is contained in:
Laurent Perron
2023-12-08 14:48:57 +01:00
parent a87803a881
commit 5f52c87b50
4 changed files with 218 additions and 109 deletions

View File

@@ -937,8 +937,9 @@ absl::StatusOr<GScipResult> GScip::Solve(
"GScip is in an error state due to a previous GScip::Solve()");
}
interrupt_event_handler_.set_interrupter(interrupter);
const absl::Cleanup interrupt_cleanup(
[this]() { interrupt_event_handler_.set_interrupter(nullptr); });
const absl::Cleanup interrupt_cleanup = [this]() {
interrupt_event_handler_.set_interrupter(nullptr);
};
// A four step process:
// 1. Apply parameters.

View File

@@ -72,9 +72,11 @@
#include "absl/strings/string_view.h"
#include "absl/synchronization/mutex.h"
#include "absl/types/span.h"
#include "ortools/base/status_macros.h"
#include "ortools/gscip/gscip.pb.h"
#include "ortools/gscip/gscip_event_handler.h"
#include "ortools/gscip/gscip_message_handler.h" // IWYU pragma: export
#include "ortools/linear_solver/scip_helper_macros.h"
#include "scip/scip.h"
#include "scip/scip_prob.h"
#include "scip/type_cons.h"
@@ -166,6 +168,7 @@ class GScip {
private:
std::atomic<bool> interrupted_{false};
};
// Create a new GScip (the constructor is private). The default objective
// direction is minimization.
static absl::StatusOr<std::unique_ptr<GScip>> Create(
@@ -413,6 +416,23 @@ class GScip {
// CHECK fails if status is OK.
void InterruptSolveFromCallbackOnCallbackError(absl::Status error_status);
// Internal use only. Users should instead call
// GScipConstraintHandler::AddCallbackConstraint().
template <typename ConsHandler, typename ConsData>
inline absl::StatusOr<SCIP_CONS*> AddConstraintForHandler(
ConsHandler* handler, ConsData* data, const std::string& name = "",
const GScipConstraintOptions& options = DefaultGScipConstraintOptions());
// Internal use only.
//
// Replaces +/- inf by +/- ScipInf(), fails when |d| is in [ScipInf(), inf).
absl::StatusOr<double> ScipInfClamp(double d);
// Internal use only.
//
// Returns +/- inf if |d| >= ScipInf(), otherwise returns d.
double ScipInfUnclamp(double d);
private:
// Event handler that it used to call SCIPinterruptSolve() is a safe manner.
//
@@ -463,12 +483,6 @@ class GScip {
absl::string_view legacy_params);
absl::Status FreeTransform();
// Replaces +/- inf by +/- ScipInf(), fails when |d| is in [ScipInf(), inf).
absl::StatusOr<double> ScipInfClamp(double d);
// Returns +/- inf if |d| >= ScipInf(), otherwise returns d.
double ScipInfUnclamp(double d);
// Returns an error if |d| >= ScipInf().
absl::Status CheckScipFinite(double d);
@@ -654,6 +668,28 @@ struct GScipConstraintOptions {
bool keep_alive = true;
};
////////////////////////////////////////////////////////////////////////////////
// Template function implementations
////////////////////////////////////////////////////////////////////////////////
template <typename ConsHandler, typename ConsData>
absl::StatusOr<SCIP_CONS*> GScip::AddConstraintForHandler(
ConsHandler* handler, ConsData* data, const std::string& name,
const GScipConstraintOptions& options) {
SCIP_CONS* constraint = nullptr;
RETURN_IF_SCIP_ERROR(SCIPcreateCons(
scip_, &constraint, name.data(), handler, data, options.initial,
options.separate, options.enforce, options.check, options.propagate,
options.local, options.modifiable, options.dynamic, options.removable,
options.sticking_at_node));
if (constraint == nullptr) {
return absl::InternalError("SCIP failed to create constraint");
}
RETURN_IF_SCIP_ERROR(SCIPaddCons(scip_, constraint));
RETURN_IF_ERROR(MaybeKeepConstraintAlive(constraint, options));
return constraint;
}
} // namespace operations_research
#endif // OR_TOOLS_GSCIP_GSCIP_H_

View File

@@ -20,7 +20,6 @@
#include "absl/log/check.h"
#include "absl/status/status.h"
#include "absl/strings/string_view.h"
#include "absl/types/span.h"
#include "ortools/base/logging.h"
#include "ortools/base/status_macros.h"
@@ -35,7 +34,6 @@
#include "scip/scip_cut.h"
#include "scip/scip_general.h"
#include "scip/scip_lp.h"
#include "scip/scip_prob.h"
#include "scip/scip_sol.h"
#include "scip/scip_solvingstats.h"
#include "scip/scip_tree.h"
@@ -66,22 +64,6 @@ using GScipHandler =
namespace operations_research {
namespace {
// Enum with supported user-implementable callback functions in the SCIP
// constraint handler. Non-user-implementable functions are not included here
// (e.g. CONSFREE). Same order as in type_cons.h.
enum class ConstraintHandlerCallbackType {
kSepaLp, // CONSSEPALP
kSepaSol, // CONSSEPASOL
kEnfoLp, // CONSENFOLP
// Unsupported: kEnfoRelax, // CONSENFORELAX
kEnfoPs, // CONSENFOPS
kConsCheck, // CONSCHECK
// Unsupported: kConsProp, // CONSPROP
// Unsupported: kConsPresol, // CONSPRESOL
// Unsupported: kConsResProp, // CONSRESPROP
kConsLock // CONSLOCK
};
// Returns the "do nothing" callback result. Used to handle the edge case where
// Enfo* callbacks need to return kFeasible instead of kDidNotRun.
GScipCallbackResult DidNotRunCallbackResult(
@@ -94,55 +76,6 @@ GScipCallbackResult DidNotRunCallbackResult(
return GScipCallbackResult::kDidNotRun;
}
// In callbacks, SCIP requires the first SCIP_RESULT in a priority list be
// returned when multiple results are applicable. This is a unified order of the
// priorities extracted from type_cons.h. The higher the result, the higher
// priority it is.
int ConstraintHandlerResultPriority(
const GScipCallbackResult result,
const ConstraintHandlerCallbackType callback_type) {
// In type_cons.h, callback results are consistently ordered across all
// constraint handler callback methods except that SCIP_SOLVELP (kSolveLp)
// takes higher priority than SCIP_BRANCHED (kBranched) in CONSENFOLP, and the
// reverse is true for CONSENFORELAX and CONSENFOLP.
switch (result) {
case GScipCallbackResult::kUnbounded:
return 14;
case GScipCallbackResult::kCutOff:
return 13;
case GScipCallbackResult::kSuccess:
return 12;
case GScipCallbackResult::kConstraintAdded:
return 11;
case GScipCallbackResult::kReducedDomain:
return 10;
case GScipCallbackResult::kSeparated:
return 9;
case GScipCallbackResult::kBranched:
return callback_type == ConstraintHandlerCallbackType::kEnfoLp ? 7 : 8;
case GScipCallbackResult::kSolveLp:
return callback_type == ConstraintHandlerCallbackType::kEnfoLp ? 8 : 7;
case GScipCallbackResult::kInfeasible:
return 6;
case GScipCallbackResult::kFeasible:
return 5;
case GScipCallbackResult::kNewRound:
return 4;
case GScipCallbackResult::kDidNotFind:
return 3;
case GScipCallbackResult::kDidNotRun:
return 2;
case GScipCallbackResult::kDelayed:
return 1;
case GScipCallbackResult::kDelayNode:
return 0;
default:
// kConstraintChanged, kFoundSolution, and kSuspend are not used in
// constraint handlers.
return -1;
}
}
constexpr GScipCallbackResult kMinPriority = GScipCallbackResult::kDelayNode;
// Calls the callback function over a span of constraints, returning the highest
@@ -207,7 +140,8 @@ absl::StatusOr<GScipCallbackResult> ApplyCallback(
return result;
}
GScipCallbackStats GetCallbackStats(SCIP* scip) {
GScipCallbackStats GetCallbackStats(GScip* gscip) {
SCIP* scip = gscip->scip();
GScipCallbackStats stats;
stats.num_processed_nodes = SCIPgetNNodes(scip);
stats.num_processed_nodes_total = SCIPgetNTotalNodes(scip);
@@ -223,6 +157,38 @@ GScipCallbackStats GetCallbackStats(SCIP* scip) {
default:
stats.current_node_id = stats.num_processed_nodes;
}
stats.primal_bound = gscip->ScipInfUnclamp(SCIPgetPrimalbound(scip));
stats.dual_bound = gscip->ScipInfUnclamp(SCIPgetDualbound(scip));
const SCIP_STAGE stage = SCIPgetStage(scip);
switch (stage) {
case SCIP_STAGE_PRESOLVED:
case SCIP_STAGE_SOLVING:
case SCIP_STAGE_SOLVED:
stats.primal_simplex_iterations = SCIPgetNPrimalLPIterations(scip);
stats.dual_simplex_iterations = SCIPgetNDualLPIterations(scip);
stats.num_nodes_left = SCIPgetNNodesLeft(scip);
break;
default:
break;
}
// SCIP counts the focus node (the current node) as explored, but to be
// consistent with gurobi, we want to count it as open instead. In particular,
// for callbacks at the root, we want num_processed_nodes=0
if (stats.num_processed_nodes > 0) {
stats.num_processed_nodes--;
stats.num_processed_nodes_total--;
stats.num_nodes_left++;
}
switch (stage) {
case SCIP_STAGE_SOLVING:
case SCIP_STAGE_SOLVED:
case SCIP_STAGE_EXITSOLVE:
stats.num_cuts_in_lp = SCIPgetNPoolCuts(scip);
break;
default:
break;
}
stats.num_solutions_found = SCIPgetNLimSolsFound(scip);
return stats;
}
@@ -245,6 +211,62 @@ GScipConstraintOptions CallbackLazyConstraintOptions(const bool local,
} // namespace
int ConstraintHandlerResultPriority(
const GScipCallbackResult result,
const ConstraintHandlerCallbackType callback_type) {
// In type_cons.h, callback results are consistently ordered across all
// constraint handler callback methods except that SCIP_SOLVELP (kSolveLp)
// takes higher priority than SCIP_BRANCHED (kBranched) in CONSENFOLP, and the
// reverse is true for CONSENFORELAX and CONSENFOLP.
switch (result) {
case GScipCallbackResult::kUnbounded:
return 14;
case GScipCallbackResult::kCutOff:
return 13;
case GScipCallbackResult::kSuccess:
return 12;
case GScipCallbackResult::kConstraintAdded:
return 11;
case GScipCallbackResult::kReducedDomain:
return 10;
case GScipCallbackResult::kSeparated:
return 9;
case GScipCallbackResult::kBranched:
return callback_type == ConstraintHandlerCallbackType::kEnfoLp ? 7 : 8;
case GScipCallbackResult::kSolveLp:
return callback_type == ConstraintHandlerCallbackType::kEnfoLp ? 8 : 7;
case GScipCallbackResult::kInfeasible:
return 6;
case GScipCallbackResult::kFeasible:
return 5;
case GScipCallbackResult::kNewRound:
return 4;
case GScipCallbackResult::kDidNotFind:
return 3;
case GScipCallbackResult::kDidNotRun:
return 2;
case GScipCallbackResult::kDelayed:
return 1;
case GScipCallbackResult::kDelayNode:
return 0;
default:
// kConstraintChanged, kFoundSolution, and kSuspend are not used in
// constraint handlers.
return -1;
}
}
GScipCallbackResult MergeConstraintHandlerResults(
const GScipCallbackResult result1, const GScipCallbackResult result2,
const ConstraintHandlerCallbackType callback_type) {
const int priority1 = ConstraintHandlerResultPriority(result1, callback_type);
const int priority2 = ConstraintHandlerResultPriority(result2, callback_type);
if (priority2 > priority1) {
return result2;
}
return result1;
}
double GScipConstraintHandlerContext::VariableValue(SCIP_VAR* variable) const {
return SCIPgetSolVal(gscip_->scip(), current_solution_, variable);
}
@@ -355,7 +377,7 @@ static SCIP_DECL_CONSENFOLP(EnforceLpC) {
SCIP_CONSHDLRDATA* scip_handler_data = SCIPconshdlrGetData(conshdlr);
operations_research::GScip* gscip = scip_handler_data->gscip;
operations_research::GScipCallbackStats stats =
operations_research::GetCallbackStats(gscip->scip());
operations_research::GetCallbackStats(gscip);
operations_research::GScipConstraintHandlerContext context(gscip, &stats,
conshdlr, nullptr);
const bool solution_known_infeasible = static_cast<bool>(solinfeasible);
@@ -380,7 +402,7 @@ static SCIP_DECL_CONSENFOPS(EnforcePseudoSolutionC) {
SCIP_CONSHDLRDATA* scip_handler_data = SCIPconshdlrGetData(conshdlr);
operations_research::GScip* gscip = scip_handler_data->gscip;
operations_research::GScipCallbackStats stats =
operations_research::GetCallbackStats(gscip->scip());
operations_research::GetCallbackStats(gscip);
operations_research::GScipConstraintHandlerContext context(gscip, &stats,
conshdlr, nullptr);
const bool solution_known_infeasible = static_cast<bool>(solinfeasible);
@@ -408,7 +430,7 @@ static SCIP_DECL_CONSCHECK(CheckFeasibilityC) {
SCIP_CONSHDLRDATA* scip_handler_data = SCIPconshdlrGetData(conshdlr);
operations_research::GScip* gscip = scip_handler_data->gscip;
operations_research::GScipCallbackStats stats =
operations_research::GetCallbackStats(gscip->scip());
operations_research::GetCallbackStats(gscip);
operations_research::GScipConstraintHandlerContext context(gscip, &stats,
conshdlr, sol);
const bool check_integrality = static_cast<bool>(checkintegrality);
@@ -437,7 +459,7 @@ static SCIP_DECL_CONSSEPALP(SeparateLpC) {
SCIP_CONSHDLRDATA* scip_handler_data = SCIPconshdlrGetData(conshdlr);
operations_research::GScip* gscip = scip_handler_data->gscip;
operations_research::GScipCallbackStats stats =
operations_research::GetCallbackStats(gscip->scip());
operations_research::GetCallbackStats(gscip);
operations_research::GScipConstraintHandlerContext context(gscip, &stats,
conshdlr, nullptr);
GScipHandler* gscip_handler = scip_handler_data->gscip_handler.get();
@@ -460,7 +482,7 @@ static SCIP_DECL_CONSSEPASOL(SeparatePrimalSolutionC) {
SCIP_CONSHDLRDATA* scip_handler_data = SCIPconshdlrGetData(conshdlr);
operations_research::GScip* gscip = scip_handler_data->gscip;
operations_research::GScipCallbackStats stats =
operations_research::GetCallbackStats(gscip->scip());
operations_research::GetCallbackStats(gscip);
operations_research::GScipConstraintHandlerContext context(gscip, &stats,
conshdlr, sol);
GScipHandler* gscip_handler = scip_handler_data->gscip_handler.get();
@@ -549,10 +571,10 @@ absl::Status RegisterConstraintHandler(
return absl::OkStatus();
}
absl::Status AddCallbackConstraint(GScip* gscip, absl::string_view handler_name,
absl::string_view constraint_name,
void* constraint_data,
const GScipConstraintOptions& options) {
absl::StatusOr<SCIP_CONS*> AddCallbackConstraint(
GScip* gscip, const std::string& handler_name,
const std::string& constraint_name, void* constraint_data,
const GScipConstraintOptions& options) {
if (constraint_data == nullptr) {
return absl::InvalidArgumentError(
"Constraint data missing when adding a constraint handler callback");
@@ -567,18 +589,9 @@ absl::Status AddCallbackConstraint(GScip* gscip, absl::string_view handler_name,
}
SCIP_CONSDATA* consdata = new SCIP_CONSDATA;
consdata->data = constraint_data;
SCIP_CONS* constraint = nullptr;
RETURN_IF_SCIP_ERROR(SCIPcreateCons(
scip, &constraint, constraint_name.data(), conshdlr, consdata,
options.initial, options.separate, options.enforce, options.check,
options.propagate, options.local, options.modifiable, options.dynamic,
options.removable, options.sticking_at_node));
if (constraint == nullptr) {
return absl::InternalError("SCIP failed to create constraint");
}
RETURN_IF_SCIP_ERROR(SCIPaddCons(scip, constraint));
RETURN_IF_SCIP_ERROR(SCIPreleaseCons(scip, &constraint));
return absl::OkStatus();
return gscip->AddConstraintForHandler(conshdlr, consdata, constraint_name,
options);
}
} // namespace internal

View File

@@ -21,13 +21,11 @@
#include <cstdint>
#include <memory>
#include <string>
#include <string_view>
#include <utility>
#include <vector>
#include "absl/status/status.h"
#include "absl/status/statusor.h"
#include "absl/strings/string_view.h"
#include "ortools/gscip/gscip.h"
#include "ortools/gscip/gscip_callback_result.h"
#include "scip/type_cons.h"
@@ -140,10 +138,35 @@ struct GScipCallbackStats {
// nodes before a restart), including the focus node. See SCIPgetNNodes().
int64_t num_processed_nodes = 0;
// The number of open nodes left in the current run. See SCIPgetNNodesLeft().
int64_t num_nodes_left = 0;
// The total number of processed nodes in all runs, including the focus node.
// If the solver restarts > 1 time, will be larger than
// num_processed_nodes, otherwise is equal. See SCIPgetNTotalNodes().
int64_t num_processed_nodes_total = 0;
// The global primal bound on the problem. See SCIPgetPrimalBound().
double primal_bound = 0.0;
// The global dual bound on the problem. See SCIPgetDualBound().
double dual_bound = 0.0;
// The number of pivots taken by the primal simplex method. See
// SCIPgetNPrimalLPIterations().
int64_t primal_simplex_iterations = 0;
// The number of pivots taken by the dual simplex method. See
// SCIPgetNDualLPIterations().
int64_t dual_simplex_iterations = 0;
// The number of feasible solutions found that were at least as good as the
// cutoff limit. See SCIPgetNLimSolsFound().
int32_t num_solutions_found = 0;
// The number of rows in the current LP that were added as cuts. See
// SCIPgetNPoolCuts().
int32_t num_cuts_in_lp = 0;
};
// Interface to the callback context and underlying problem. Supports adding
@@ -281,10 +304,10 @@ class GScipConstraintHandler {
// Adds a callback constraint to the model. That is, it attaches to the
// constraint handler a constraint for the given constraint data.
absl::Status AddCallbackConstraint(GScip* gscip,
std::string_view constraint_name,
const ConstraintData* constraint_data,
const GScipConstraintOptions& options);
absl::StatusOr<SCIP_CONS*> AddCallbackConstraint(
GScip* gscip, const std::string& constraint_name,
const ConstraintData* constraint_data,
const GScipConstraintOptions& options = DefaultGScipConstraintOptions());
// Callback function called at SCIP's CONSENFOLP. Must check if an LP solution
// at a node is feasible, and if not, resolve the infeasibility if possible by
@@ -418,6 +441,41 @@ class GScipConstraintHandler {
GScipConstraintHandlerProperties properties_;
};
////////////////////////////////////////////////////////////////////////////////
// Helpers for implementing callbacks
////////////////////////////////////////////////////////////////////////////////
// Enum with supported user-implementable callback functions in the SCIP
// constraint handler. Non-user-implementable functions are not included here
// (e.g. CONSFREE). Same order as in type_cons.h.
enum class ConstraintHandlerCallbackType {
kSepaLp, // CONSSEPALP
kSepaSol, // CONSSEPASOL
kEnfoLp, // CONSENFOLP
// Unsupported: kEnfoRelax, // CONSENFORELAX
kEnfoPs, // CONSENFOPS
kConsCheck, // CONSCHECK
// Unsupported: kConsProp, // CONSPROP
// Unsupported: kConsPresol, // CONSPRESOL
// Unsupported: kConsResProp, // CONSRESPROP
kConsLock // CONSLOCK
};
// In callbacks, SCIP requires the first SCIP_RESULT in a priority list be
// returned when multiple results are applicable. This is a unified order of the
// priorities extracted from type_cons.h. The higher the result, the higher
// priority it is.
//
// Larger number indicates higher priority in what should be returned.
int ConstraintHandlerResultPriority(
GScipCallbackResult result, ConstraintHandlerCallbackType callback_type);
// Returns the GScipCallbackResult with larger priority, see
// ConstraintHandlerResultPriority().
GScipCallbackResult MergeConstraintHandlerResults(
GScipCallbackResult result1, GScipCallbackResult result2,
ConstraintHandlerCallbackType callback_type);
namespace internal {
// These classes implement a void pointer version of GScipConstraintHandler that
@@ -495,10 +553,10 @@ absl::Status RegisterConstraintHandler(
GScip* gscip,
std::unique_ptr<UntypedGScipConstraintHandler> constraint_handler);
absl::Status AddCallbackConstraint(GScip* gscip, std::string_view handler_name,
std::string_view constraint_name,
void* constraint_data,
const GScipConstraintOptions& options);
absl::StatusOr<SCIP_CONS*> AddCallbackConstraint(
GScip* gscip, const std::string& handler_name,
const std::string& constraint_name, void* constraint_data,
const GScipConstraintOptions& options);
} // namespace internal
@@ -513,8 +571,9 @@ absl::Status GScipConstraintHandler<ConstraintData>::Register(GScip* gscip) {
}
template <typename ConstraintData>
absl::Status GScipConstraintHandler<ConstraintData>::AddCallbackConstraint(
GScip* gscip, std::string_view constraint_name,
absl::StatusOr<SCIP_CONS*>
GScipConstraintHandler<ConstraintData>::AddCallbackConstraint(
GScip* gscip, const std::string& constraint_name,
const ConstraintData* constraint_data,
const GScipConstraintOptions& options) {
return internal::AddCallbackConstraint(