From 524b33cc3f112d9c497788024fb7d90096af839c Mon Sep 17 00:00:00 2001 From: Corentin Le Molgat Date: Mon, 9 Oct 2023 14:23:57 +0200 Subject: [PATCH 1/5] update gscip from google3 --- ortools/gscip/gscip.cc | 54 +- ortools/gscip/gscip.h | 25 +- ortools/gscip/gscip_callback_result.cc | 61 ++ ortools/gscip/gscip_callback_result.h | 53 ++ ortools/gscip/gscip_constraint_handler.cc | 586 +++++++++++++++++++ ortools/gscip/gscip_constraint_handler.h | 653 ++++++++++++++++++++++ 6 files changed, 1429 insertions(+), 3 deletions(-) create mode 100644 ortools/gscip/gscip_callback_result.cc create mode 100644 ortools/gscip/gscip_callback_result.h create mode 100644 ortools/gscip/gscip_constraint_handler.cc create mode 100644 ortools/gscip/gscip_constraint_handler.h diff --git a/ortools/gscip/gscip.cc b/ortools/gscip/gscip.cc index 76e4fa8127b..e9a2301cf3f 100644 --- a/ortools/gscip/gscip.cc +++ b/ortools/gscip/gscip.cc @@ -17,7 +17,6 @@ #include #include -#include #include #include #include @@ -26,12 +25,14 @@ #include "absl/container/flat_hash_map.h" #include "absl/container/flat_hash_set.h" +#include "absl/log/check.h" #include "absl/memory/memory.h" #include "absl/status/status.h" #include "absl/status/statusor.h" #include "absl/strings/str_cat.h" #include "absl/strings/str_format.h" #include "absl/strings/string_view.h" +#include "absl/synchronization/mutex.h" #include "absl/types/span.h" #include "ortools/base/logging.h" #include "ortools/base/status_macros.h" @@ -43,14 +44,27 @@ #include "ortools/util/status_macros.h" #include "scip/cons_linear.h" #include "scip/cons_quadratic.h" +#include "scip/def.h" #include "scip/scip.h" +#include "scip/scip_cons.h" #include "scip/scip_general.h" +#include "scip/scip_message.h" +#include "scip/scip_numerics.h" #include "scip/scip_param.h" #include "scip/scip_prob.h" +#include "scip/scip_sol.h" +#include "scip/scip_solve.h" #include "scip/scip_solvingstats.h" +#include "scip/scip_var.h" #include "scip/scipdefplugins.h" +#include "lpi/lpi.h" #include "scip/type_cons.h" +#include "scip/type_paramset.h" +#include "scip/type_retcode.h" #include "scip/type_scip.h" +#include "scip/type_set.h" +#include "scip/type_sol.h" +#include "scip/type_stat.h" #include "scip/type_var.h" namespace operations_research { @@ -304,6 +318,18 @@ bool GScip::InterruptSolve() { return SCIPinterruptSolve(scip_) == SCIP_OKAY; } +void GScip::InterruptSolveFromCallback(absl::Status error_status) { + CHECK(!error_status.ok()); + { + const absl::MutexLock lock(&callback_status_mutex_); + if (!callback_status_.ok()) { + return; + } + callback_status_ = std::move(error_status); + } + InterruptSolve(); +} + absl::Status GScip::CleanUp() { if (scip_ != nullptr) { for (SCIP_VAR* variable : variables_) { @@ -832,6 +858,11 @@ absl::StatusOr GScip::SuggestHint( absl::StatusOr GScip::Solve( const GScipParameters& params, absl::string_view legacy_params, const GScipMessageHandler message_handler) { + if (InErrorState()) { + return absl::InvalidArgumentError( + "GScip is in an error state due to a previous GScip::Solve()"); + } + // A four step process: // 1. Apply parameters. // 2. Solve the problem. @@ -919,6 +950,22 @@ absl::StatusOr GScip::Solve( " when writing solve stats.")); } } + absl::Status callback_status; + { + const absl::MutexLock callback_status_lock(&callback_status_mutex_); + callback_status = util::StatusBuilder(callback_status_) + << "error in a callback that interrupted the solve"; + } + if (!callback_status.ok()) { + const absl::Status status = FreeTransform(); + if (status.ok()) { + return callback_status; + } + LOG(ERROR) << "GScip::FreeTransform() failed after interrupting " + "the solve due to an error in a callback: " + << callback_status; + return status; + } // Step 3: Extract solution information. // Some outputs are available unconditionally, and some are only ready if at // least presolve succeeded. @@ -1070,6 +1117,11 @@ absl::Status GScip::CheckScipFinite(double d) { return absl::OkStatus(); } +bool GScip::InErrorState() { + const absl::MutexLock lock(&callback_status_mutex_); + return !callback_status_.ok(); +} + #undef RETURN_ERROR_UNLESS } // namespace operations_research diff --git a/ortools/gscip/gscip.h b/ortools/gscip/gscip.h index 9ed179457d7..99a47ef5ce6 100644 --- a/ortools/gscip/gscip.h +++ b/ortools/gscip/gscip.h @@ -26,7 +26,9 @@ // but < inf result in an error. Changing the underlying SCIP's infinity is // not supported. // * absl::Status and absl::StatusOr are used to propagate SCIP errors (and on -// a best effort basis, also filter out bad input to gSCIP functions). +// a best effort basis, also filter out bad input to gSCIP functions). In +// constraint handlers, we also use absl::Status and absl::StatusOr for +// error propagation and not SCIP_RETCODE. // // A note on error propagation and reliability: // Many methods on SCIP return an error code. Errors can be triggered by @@ -55,11 +57,13 @@ #include #include +#include "absl/base/thread_annotations.h" #include "absl/container/flat_hash_map.h" #include "absl/container/flat_hash_set.h" #include "absl/status/status.h" #include "absl/status/statusor.h" #include "absl/strings/string_view.h" +#include "absl/synchronization/mutex.h" #include "absl/types/span.h" #include "ortools/gscip/gscip.pb.h" #include "ortools/gscip/gscip_message_handler.h" // IWYU pragma: export @@ -123,7 +127,8 @@ enum class GScipHintResult; // A thin wrapper around the SCIP solver that provides C++ bindings that are // idiomatic for Google. Unless callbacks are used, the SCIP stage is always -// PROBLEM. +// PROBLEM. If any GScip function returns an absl::Status error, then the GScip +// object should be considered to be in an error state. class GScip { public: // Create a new GScip (the constructor is private). The default objective @@ -140,6 +145,7 @@ class GScip { // The returned StatusOr will contain an error only if an: // * An underlying function from SCIP fails. // * There is an I/O error with managing SCIP output. + // * A user-defined callback function fails. // The above cases are not mutually exclusive. If the problem is infeasible, // this will be reflected in the value of GScipResult::gscip_output::status. absl::StatusOr Solve( @@ -362,6 +368,19 @@ class GScip { absl::StatusOr DefaultStringParamValue( const std::string& parameter_name); + // Returns true if GScip is in an error state. Currently only checks if a + // callback has failed, but in the future it may check for other failures. + bool InErrorState(); + + // Internal use. Used by constraint handlers to propagate user-returned Status + // errors to GSCIP. Interrupts a solve and passes an error status that + // GScip::Solve() must return after SCIP finishes interrupting the solve. Does + // not do anything if previously called with a (non-OK) status (i.e. the first + // status error is returned by SCIP). + // + // CHECK fails if status is OK. + void InterruptSolveFromCallback(absl::Status error_status); + private: explicit GScip(SCIP* scip); // Releases SCIP memory. @@ -386,6 +405,8 @@ class GScip { SCIP* scip_; absl::flat_hash_set variables_; absl::flat_hash_set constraints_; + absl::Mutex callback_status_mutex_; + absl::Status callback_status_ ABSL_GUARDED_BY(callback_status_mutex_); }; // Advanced features below diff --git a/ortools/gscip/gscip_callback_result.cc b/ortools/gscip/gscip_callback_result.cc new file mode 100644 index 00000000000..15606aa8ec5 --- /dev/null +++ b/ortools/gscip/gscip_callback_result.cc @@ -0,0 +1,61 @@ +// Copyright 2010-2022 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "ortools/gscip/gscip_callback_result.h" + +#include "scip/type_result.h" + +namespace operations_research { + +SCIP_RESULT ConvertGScipCallbackResult(const GScipCallbackResult result) { + switch (result) { + case GScipCallbackResult::kDidNotRun: + return SCIP_DIDNOTRUN; + case GScipCallbackResult::kDelayed: + return SCIP_DELAYED; + case GScipCallbackResult::kDidNotFind: + return SCIP_DIDNOTFIND; + case GScipCallbackResult::kFeasible: + return SCIP_FEASIBLE; + case GScipCallbackResult::kInfeasible: + return SCIP_INFEASIBLE; + case GScipCallbackResult::kUnbounded: + return SCIP_UNBOUNDED; + case GScipCallbackResult::kCutOff: + return SCIP_CUTOFF; + case GScipCallbackResult::kSeparated: + return SCIP_SEPARATED; + case GScipCallbackResult::kNewRound: + return SCIP_NEWROUND; + case GScipCallbackResult::kReducedDomain: + return SCIP_REDUCEDDOM; + case GScipCallbackResult::kConstraintAdded: + return SCIP_CONSADDED; + case GScipCallbackResult::kConstraintChanged: + return SCIP_CONSCHANGED; + case GScipCallbackResult::kBranched: + return SCIP_BRANCHED; + case GScipCallbackResult::kSolveLp: + return SCIP_SOLVELP; + case GScipCallbackResult::kFoundSolution: + return SCIP_FOUNDSOL; + case GScipCallbackResult::kSuspend: + return SCIP_SUSPENDED; + case GScipCallbackResult::kSuccess: + return SCIP_SUCCESS; + case GScipCallbackResult::kDelayNode: + return SCIP_DELAYNODE; + } +} + +} // namespace operations_research diff --git a/ortools/gscip/gscip_callback_result.h b/ortools/gscip/gscip_callback_result.h new file mode 100644 index 00000000000..27676b7ca2e --- /dev/null +++ b/ortools/gscip/gscip_callback_result.h @@ -0,0 +1,53 @@ +// Copyright 2010-2022 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#ifndef OR_TOOLS_GSCIP_GSCIP_CALLBACK_RESULT_H_ +#define OR_TOOLS_GSCIP_GSCIP_CALLBACK_RESULT_H_ + +#include "scip/type_result.h" +#include "scip/type_retcode.h" + +namespace operations_research { + +// Equivalent to type_result.h in SCIP. +enum class GScipCallbackResult { + kDidNotRun = 1, // The method was not executed. + kDelayed = 2, // The method was not executed, should be called again later. + kDidNotFind = 3, // The method was executed, but failed finding anything. + kFeasible = 4, // No infeasibility could be found. + kInfeasible = 5, // An infeasibility was detected. + kUnbounded = 6, // An unboundedness was detected. + kCutOff = 7, // The current node is infeasible and can be cut off. + kSeparated = 8, // The method added a cutting plane. + // The method added a cutting plane and a new separation round should + // immediately start. + kNewRound = 9, + kReducedDomain = 10, // The method reduced the domain of a variable. + kConstraintAdded = 11, // The method added a constraint. + kConstraintChanged = 12, // The method changed a constraint. + kBranched = 13, // The method created a branching. + kSolveLp = 14, // The current node's LP must be solved. + kFoundSolution = 15, // The method found a feasible primal solution. + // The method interrupted its execution, but can continue if needed. + kSuspend = 16, + kSuccess = 17, // The method was successfully executed. + // The processing of the branch-and-bound node should stopped and continued + // later. + kDelayNode = 18 +}; + +SCIP_RESULT ConvertGScipCallbackResult(GScipCallbackResult result); + +} // namespace operations_research + +#endif // OR_TOOLS_GSCIP_GSCIP_CALLBACK_RESULT_H_ diff --git a/ortools/gscip/gscip_constraint_handler.cc b/ortools/gscip/gscip_constraint_handler.cc new file mode 100644 index 00000000000..f5452490416 --- /dev/null +++ b/ortools/gscip/gscip_constraint_handler.cc @@ -0,0 +1,586 @@ +// Copyright 2010-2022 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "ortools/gscip/gscip_constraint_handler.h" + +#include +#include +#include +#include + +#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" +#include "ortools/gscip/gscip.h" +#include "ortools/gscip/gscip_callback_result.h" +#include "ortools/linear_solver/scip_helper_macros.h" +#include "scip/def.h" +#include "scip/pub_cons.h" +#include "scip/pub_message.h" +#include "scip/pub_var.h" +#include "scip/scip_cons.h" +#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" +#include "scip/scip_var.h" +#include "scip/struct_cons.h" +#include "scip/struct_tree.h" // IWYU pragma: keep +#include "scip/type_cons.h" +#include "scip/type_lp.h" +#include "scip/type_retcode.h" +#include "scip/type_scip.h" +#include "scip/type_set.h" +#include "scip/type_tree.h" +#include "scip/type_var.h" + +struct SCIP_ConshdlrData { + std::unique_ptr + gscip_handler; + operations_research::GScip* gscip = nullptr; +}; + +struct SCIP_ConsData { + void* data; +}; + +using GScipHandler = + operations_research::internal::UntypedGScipConstraintHandler; + +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( + const ConstraintHandlerCallbackType callback_type) { + // TODO(user): Add kEnfoRelax when we support it. + if (callback_type == ConstraintHandlerCallbackType::kEnfoLp || + callback_type == ConstraintHandlerCallbackType::kEnfoPs) { + return GScipCallbackResult::kFeasible; + } + 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 +// priority callback result, along with a SCIP return code. +absl::StatusOr ApplyCallback( + absl::Span constraints, + std::function callback_function, + const ConstraintHandlerCallbackType callback_type) { + if (constraints.empty()) { + return DidNotRunCallbackResult(callback_type); + } + GScipCallbackResult callback_result = kMinPriority; + for (SCIP_Cons* cons : constraints) { + if (cons == nullptr) { + return absl::InternalError("Constraint handler has null constraint"); + } + SCIP_CONSDATA* consdata = SCIPconsGetData(cons); + if (consdata == nullptr || consdata->data == nullptr) { + return absl::InternalError("Constraint handler has null constraint data"); + } + const GScipCallbackResult cons_result = callback_function(consdata->data); + if (ConstraintHandlerResultPriority(cons_result, callback_type) > + ConstraintHandlerResultPriority(callback_result, callback_type)) { + callback_result = cons_result; + } + } + return callback_result; +} + +// Calls the callback function over all the constraints of a constraint handler, +// prioritizing the ones SCIP deems more useful. Returns the highest priority +// callback result, along with a SCIP return code. +absl::StatusOr ApplyCallback( + SCIP_Cons** constraints, const int num_useful_constraints, + const int total_num_constraints, + std::function callback_function, + const ConstraintHandlerCallbackType callback_type) { + absl::Span all_constraints = + absl::MakeSpan(constraints, total_num_constraints); + absl::Span useful_constraints = + all_constraints.subspan(0, num_useful_constraints); + ASSIGN_OR_RETURN( + const GScipCallbackResult result, + ApplyCallback(useful_constraints, callback_function, callback_type)); + // The first num_useful_constraints are the ones that are more likely to be + // violated. If no violation was found, we consider the remaining + // constraints. + if (result == GScipCallbackResult::kDidNotFind || + result == GScipCallbackResult::kDidNotRun || + result == GScipCallbackResult::kFeasible) { + absl::Span remaining_constraints = + all_constraints.subspan(num_useful_constraints); + ASSIGN_OR_RETURN( + const GScipCallbackResult remaining_result, + ApplyCallback(remaining_constraints, callback_function, callback_type)); + if (remaining_result != GScipCallbackResult::kDidNotFind && + remaining_result != GScipCallbackResult::kDidNotRun && + remaining_result != GScipCallbackResult::kFeasible) { + return remaining_result; + } + } + return result; +} + +GScipCallbackStats GetCallbackStats(SCIP* scip) { + GScipCallbackStats stats; + stats.num_processed_nodes = SCIPgetNNodes(scip); + stats.num_processed_nodes_total = SCIPgetNTotalNodes(scip); + switch (SCIPgetStage(scip)) { + case SCIP_STAGE_INITPRESOLVE: + case SCIP_STAGE_PRESOLVING: + case SCIP_STAGE_EXITPRESOLVE: + case SCIP_STAGE_SOLVING: { + SCIP_NODE* node = SCIPgetCurrentNode(scip); + stats.current_node_id = node == nullptr ? -1 : node->number; + break; + } + default: + stats.current_node_id = stats.num_processed_nodes; + } + return stats; +} + +GScipConstraintOptions CallbackLazyConstraintOptions(const bool local, + const bool dynamic) { + GScipConstraintOptions result; + result.initial = true; + result.separate = true; + result.enforce = true; + result.check = true; + result.propagate = true; + result.local = local; + result.modifiable = false; + result.dynamic = dynamic; + result.removable = true; + result.sticking_at_node = false; + result.keep_alive = false; + return result; +} + +} // namespace + +double GScipConstraintHandlerContext::VariableValue(SCIP_VAR* variable) const { + return SCIPgetSolVal(gscip_->scip(), current_solution_, variable); +} + +absl::StatusOr GScipConstraintHandlerContext::AddCut( + const GScipLinearRange& range, const std::string& name, + const GScipCutOptions& options) { + SCIP* scip = gscip_->scip(); + SCIP_ROW* row = nullptr; + RETURN_IF_SCIP_ERROR(SCIPcreateEmptyRowConshdlr( + scip, &row, current_handler_, name.data(), range.lower_bound, + range.upper_bound, options.local, options.modifiable, options.removable)); + RETURN_IF_SCIP_ERROR(SCIPcacheRowExtensions(scip, row)); + if (range.coefficients.size() != range.variables.size()) { + return absl::InternalError( + "GScipLinearRange variables and coefficients do not match in size"); + } + for (int i = 0; i < range.variables.size(); ++i) { + RETURN_IF_SCIP_ERROR(SCIPaddVarToRow(scip, row, range.variables.at(i), + range.coefficients.at(i))); + } + RETURN_IF_SCIP_ERROR(SCIPflushRowExtensions(scip, row)); + SCIP_Bool infeasible; + RETURN_IF_SCIP_ERROR(SCIPaddRow(scip, row, options.force_cut, &infeasible)); + RETURN_IF_SCIP_ERROR(SCIPreleaseRow(scip, &row)); + return infeasible ? GScipCallbackResult::kCutOff + : GScipCallbackResult::kSeparated; +} + +absl::Status GScipConstraintHandlerContext::AddLazyLinearConstraint( + const GScipLinearRange& range, const std::string& name, + const GScipLazyConstraintOptions& options) { + return gscip_ + ->AddLinearConstraint( + range, name, + CallbackLazyConstraintOptions(options.local, options.dynamic)) + .status(); +} + +double GScipConstraintHandlerContext::LocalVarLb(SCIP_VAR* var) const { + return SCIPvarGetLbLocal(var); +} +double GScipConstraintHandlerContext::LocalVarUb(SCIP_VAR* var) const { + return SCIPvarGetUbLocal(var); +} +double GScipConstraintHandlerContext::GlobalVarLb(SCIP_VAR* var) const { + return SCIPvarGetLbGlobal(var); +} +double GScipConstraintHandlerContext::GlobalVarUb(SCIP_VAR* var) const { + return SCIPvarGetUbGlobal(var); +} + +absl::Status GScipConstraintHandlerContext::SetLocalVarLb(SCIP_VAR* var, + double value) { + return SCIP_TO_STATUS( + SCIPchgVarLbNode(gscip_->scip(), /*node=*/nullptr, var, value)); +} +absl::Status GScipConstraintHandlerContext::SetLocalVarUb(SCIP_VAR* var, + double value) { + return SCIP_TO_STATUS( + SCIPchgVarUbNode(gscip_->scip(), /*node=*/nullptr, var, value)); +} +absl::Status GScipConstraintHandlerContext::SetGlobalVarLb(SCIP_VAR* var, + double value) { + return SCIP_TO_STATUS(SCIPchgVarLbGlobal(gscip_->scip(), var, value)); +} +absl::Status GScipConstraintHandlerContext::SetGlobalVarUb(SCIP_VAR* var, + double value) { + return SCIP_TO_STATUS(SCIPchgVarUbGlobal(gscip_->scip(), var, value)); +} + +} // namespace operations_research + +// /////////////////////////////////////////////////////////////////////////// +// SCIP callback implementation +// ////////////////////////////////////////////////////////////////////////// + +extern "C" { + +// Destructor of the constraint handler to free user data (called when SCIP is +// exiting). +static SCIP_DECL_CONSFREE(ConstraintHandlerFreeC) { + if (scip == nullptr) { + SCIPerrorMessage("SCIP not found in SCIP_DECL_CONSFREE"); + return SCIP_ERROR; + } + SCIP_CONSHDLRDATA* scip_handler_data = SCIPconshdlrGetData(conshdlr); + if (scip_handler_data == nullptr) { + SCIPerrorMessage("SCIP handler data not found in SCIP_DECL_CONSFREE"); + return SCIP_ERROR; + } + delete scip_handler_data; + SCIPconshdlrSetData(conshdlr, nullptr); + return SCIP_OKAY; +} + +static SCIP_DECL_CONSDELETE(ConstraintDataDeleteC) { + if (consdata == nullptr || *consdata == nullptr) { + SCIPerrorMessage("SCIP constraint data not found in SCIP_DECL_CONSDELETE"); + return SCIP_ERROR; + } + delete *consdata; + cons->consdata = nullptr; + return SCIP_OKAY; +} + +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::GScipConstraintHandlerContext context(gscip, &stats, + conshdlr, nullptr); + const bool solution_known_infeasible = static_cast(solinfeasible); + GScipHandler* gscip_handler = scip_handler_data->gscip_handler.get(); + auto DoEnforceLp = [=](void* constraint_data) { + return gscip_handler->CallEnforceLp(context, constraint_data, + solution_known_infeasible); + }; + const absl::StatusOr gresult = + operations_research::ApplyCallback( + conss, nusefulconss, nconss, DoEnforceLp, + operations_research::ConstraintHandlerCallbackType::kEnfoLp); + if (!gresult.ok()) { + SCIPerrorMessage(gresult.status().ToString().c_str()); + return SCIP_ERROR; + } + *result = operations_research::ConvertGScipCallbackResult(*gresult); + return SCIP_OKAY; +} + +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::GScipConstraintHandlerContext context(gscip, &stats, + conshdlr, nullptr); + const bool solution_known_infeasible = static_cast(solinfeasible); + const bool solution_infeasible_by_objective = + static_cast(objinfeasible); + GScipHandler* gscip_handler = scip_handler_data->gscip_handler.get(); + auto DoEnforcePseudoSolution = [=](void* constraint_data) { + return gscip_handler->CallEnforcePseudoSolution( + context, constraint_data, solution_known_infeasible, + solution_infeasible_by_objective); + }; + const absl::StatusOr gresult = + operations_research::ApplyCallback( + conss, nusefulconss, nconss, DoEnforcePseudoSolution, + operations_research::ConstraintHandlerCallbackType::kEnfoPs); + if (!gresult.ok()) { + SCIPerrorMessage(gresult.status().ToString().c_str()); + return SCIP_ERROR; + } + *result = operations_research::ConvertGScipCallbackResult(*gresult); + return SCIP_OKAY; +} + +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::GScipConstraintHandlerContext context(gscip, &stats, + conshdlr, sol); + const bool check_integrality = static_cast(checkintegrality); + const bool check_lp_rows = static_cast(checklprows); + const bool print_reason = static_cast(printreason); + const bool complete = static_cast(completely); + GScipHandler* gscip_handler = scip_handler_data->gscip_handler.get(); + auto DoCheckIsFeasible = [=](void* constraint_data) { + return gscip_handler->CallCheckIsFeasible(context, constraint_data, + check_integrality, check_lp_rows, + print_reason, complete); + }; + const absl::StatusOr gresult = + operations_research::ApplyCallback( + conss, nconss, nconss, DoCheckIsFeasible, + operations_research::ConstraintHandlerCallbackType::kEnfoPs); + if (!gresult.ok()) { + SCIPerrorMessage(gresult.status().ToString().c_str()); + return SCIP_ERROR; + } + *result = operations_research::ConvertGScipCallbackResult(*gresult); + return SCIP_OKAY; +} + +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::GScipConstraintHandlerContext context(gscip, &stats, + conshdlr, nullptr); + GScipHandler* gscip_handler = scip_handler_data->gscip_handler.get(); + auto DoSeparateLp = [=](void* constraint_data) { + return gscip_handler->CallSeparateLp(context, constraint_data); + }; + const absl::StatusOr gresult = + operations_research::ApplyCallback( + conss, nusefulconss, nconss, DoSeparateLp, + operations_research::ConstraintHandlerCallbackType::kSepaLp); + if (!gresult.ok()) { + SCIPerrorMessage(gresult.status().ToString().c_str()); + return SCIP_ERROR; + } + *result = operations_research::ConvertGScipCallbackResult(*gresult); + return SCIP_OKAY; +} + +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::GScipConstraintHandlerContext context(gscip, &stats, + conshdlr, sol); + GScipHandler* gscip_handler = scip_handler_data->gscip_handler.get(); + auto DoSeparateSolution = [=](void* constraint_data) { + return gscip_handler->CallSeparateSolution(context, constraint_data); + }; + const absl::StatusOr gresult = + operations_research::ApplyCallback( + conss, nusefulconss, nconss, DoSeparateSolution, + operations_research::ConstraintHandlerCallbackType::kSepaSol); + if (!gresult.ok()) { + SCIPerrorMessage(gresult.status().ToString().c_str()); + return SCIP_ERROR; + } + *result = operations_research::ConvertGScipCallbackResult(*gresult); + return SCIP_OKAY; +} + +static SCIP_DECL_CONSLOCK(VariableRoundingLockC) { + SCIP_CONSHDLRDATA* scip_handler_data = SCIPconshdlrGetData(conshdlr); + operations_research::GScip* gscip = scip_handler_data->gscip; + GScipHandler* gscip_handler = scip_handler_data->gscip_handler.get(); + SCIP_CONSDATA* consdata = SCIPconsGetData(cons); + if (consdata == nullptr) { + return SCIP_ERROR; + } + const bool lock_type_is_model = locktype == SCIP_LOCKTYPE_MODEL; + for (const auto [locked_var, lock_direction] : + gscip_handler->RoundingLock(gscip, consdata->data, lock_type_is_model)) { + int lock_down; + int lock_up; + switch (lock_direction) { + case operations_research::RoundingLockDirection::kUp: + lock_down = nlocksneg; + lock_up = nlockspos; + break; + case operations_research::RoundingLockDirection::kDown: + lock_down = nlockspos; + lock_up = nlocksneg; + break; + case operations_research::RoundingLockDirection::kBoth: + lock_down = nlocksneg + nlockspos; + lock_up = nlocksneg + nlockspos; + break; + } + SCIP_CALL( + SCIPaddVarLocksType(scip, locked_var, locktype, lock_down, lock_up)); + } + return SCIP_OKAY; +} +} + +namespace operations_research { + +namespace internal { + +absl::Status RegisterConstraintHandler( + GScip* gscip, + std::unique_ptr constraint_handler) { + SCIP_CONSHDLR* c_scip_handler; + SCIP_CONSHDLRDATA* scip_handler_data = new SCIP_CONSHDLRDATA; + scip_handler_data->gscip_handler = std::move(constraint_handler); + scip_handler_data->gscip = gscip; + SCIP* scip = gscip->scip(); + const GScipConstraintHandlerProperties& properties = + scip_handler_data->gscip_handler->properties(); + + RETURN_IF_SCIP_ERROR(SCIPincludeConshdlrBasic( + scip, &c_scip_handler, properties.name.c_str(), + properties.description.c_str(), properties.enforcement_priority, + properties.feasibility_check_priority, properties.eager_frequency, + properties.needs_constraints, EnforceLpC, EnforcePseudoSolutionC, + CheckFeasibilityC, VariableRoundingLockC, scip_handler_data)); + if (c_scip_handler == nullptr) { + return absl::InternalError("SCIP failed to add constraint handler"); + } + RETURN_IF_SCIP_ERROR(SCIPsetConshdlrSepa( + scip, c_scip_handler, SeparateLpC, SeparatePrimalSolutionC, + properties.separation_frequency, properties.separation_priority, + properties.delay_separation)); + RETURN_IF_SCIP_ERROR( + SCIPsetConshdlrFree(scip, c_scip_handler, ConstraintHandlerFreeC)); + RETURN_IF_SCIP_ERROR( + SCIPsetConshdlrDelete(scip, c_scip_handler, ConstraintDataDeleteC)); + + return absl::OkStatus(); +} + +absl::Status AddCallbackConstraint(GScip* gscip, absl::string_view handler_name, + absl::string_view constraint_name, + void* constraint_data, + const GScipConstraintOptions& options) { + if (constraint_data == nullptr) { + return absl::InvalidArgumentError( + "Constraint data missing when adding a constraint handler callback"); + } + SCIP* scip = gscip->scip(); + SCIP_CONSHDLR* conshdlr = SCIPfindConshdlr(scip, handler_name.data()); + if (conshdlr == nullptr) { + return util::InternalErrorBuilder() + << "Constraint handler " << handler_name + << " not registered with SCIP. Check if you " + "registered the constraint handler before adding constraints."; + } + 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(); +} + +} // namespace internal + +} // namespace operations_research diff --git a/ortools/gscip/gscip_constraint_handler.h b/ortools/gscip/gscip_constraint_handler.h new file mode 100644 index 00000000000..bdaf45180ef --- /dev/null +++ b/ortools/gscip/gscip_constraint_handler.h @@ -0,0 +1,653 @@ +// Copyright 2010-2022 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// Provides a safe C++ interface for SCIP constraint handlers, which are +// described at https://www.scipopt.org/doc/html/CONS.php. For instructions to +// write a constraint handler, see the documentation of GScipConstraintHandler. +// Examples can be found in gscip_constraint_handler_test.cc. +#ifndef OR_TOOLS_GSCIP_GSCIP_CONSTRAINT_HANDLER_H_ +#define OR_TOOLS_GSCIP_GSCIP_CONSTRAINT_HANDLER_H_ + +#include +#include +#include +#include +#include + +#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" +#include "scip/type_sol.h" +#include "scip/type_var.h" + +namespace operations_research { + +// Properties for the constraint handler. It is recommended to set priorities +// and frequencies manually. +struct GScipConstraintHandlerProperties { + // For members below, the corresponding SCIP constraint handler property name + // is provided. See https://www.scipopt.org/doc/html/CONS.php#CONS_PROPERTIES + // for details. + // + // While it is recommended to set your own parameters, we provide a default + // set of parameters that has the following behavior: + // * Enforcement and feasibility checking is done right after enforcing + // integrality, but before any other constraint handlers. This implies that + // it is only performed on integer solutions by default. + // * Obsolete constraints are revisited every 100 nodes (eager frequency). + // This default follows the most common frequency in SCIP's existing + // constraint handlers. + // * If separation is used, it is run before all constraint handlers and at + // every node. Note however that all separators are always run before any + // constraint handler separation. A user may control separation frequency + // either by changing this parameter or implementing a check in the + // callback. + + // Name of the constraint handler. See CONSHDLR_NAME in SCIP + // documentation above. + std::string name; + + // Description of the constraint handler. See CONSHDLR_DESC in SCIP + // documentation above. + std::string description; + + // Determines the order this constraint class is checked at each LP node. If + // negative, the enforcement is only performed on integer solutions. See + // CONSHDLR_ENFOPRIORITY in the SCIP documentation above. Only relevant if + // enforcement callbacks are implemented. + int enforcement_priority = -1; + + // Determines the order this constraint class runs in when testing solution + // feasibility. If negative, the feasibility check is only performed on + // integer solutions. See CONSHDLR_CHECKPRIORITY in the SCIP documentation + // above. Only relevant if check callback is implemented. + int feasibility_check_priority = -1; + + // Determines the order the separation from this constraint handler runs in + // the cut loop. Note that separators are run before constraint handlers. + // See CONSHDLR_SEPAPRIORITY in SCIP documentation above. Only relevant if + // separation callbacks are implemented. + int separation_priority = 3000000; + + // Frequency for separating cuts. See CONSHDLR_SEPAFREQ in the SCIP + // documentation above. Only relevant if separation callbacks are implemented. + int separation_frequency = 1; + + // Determines if this separator be delayed if another separator has already + // found a cut. See CONSHDLR_DELAYSEPA in the SCIP documentation above. + bool delay_separation = false; + + // Frequency for using all instead of only the useful constraints in + // separation, propagation, and enforcement. For example, some constraints may + // be aged out by SCIP if they are not relevant for several iterations. See + // CONSHDLR_EAGERFREQ in SCIP documentation above. + int eager_frequency = 100; + + // Indicates whether the constraint handler can be skipped if no constraints + // from this handler are present in the model. In most cases, this should be + // true. This should only be false for constraints that are not added + // explicitly as a constraint, such as integrality. See CONSHDLR_NEEDSCONS in + // SCIP documentation above. + bool needs_constraints = true; +}; + +// Advanced use only. Indicates that if a variable moves in this direction, it +// can cause a constraint violation. `kBoth` is the safest option and always +// valid, but it is the least flexible for SCIP. +enum class RoundingLockDirection { kUp, kDown, kBoth }; + +// Options passed to SCIP when adding a cut. +struct GScipCutOptions { + // Cut is only valid for the current subtree. + bool local = false; + // Cut is modifiable during node processing (subject to column generation). + bool modifiable = false; + // Cut can be removed from the LP due to aging or cleanup. + bool removable = true; + // Cut is forced to enter the LP. + bool force_cut = false; +}; + +// Options passed to SCIP when adding a lazy constraint. +struct GScipLazyConstraintOptions { + // Cut is only valid for the current subtree. + bool local = false; + // Constraint is subject to aging. + bool dynamic = false; +}; + +struct GScipCallbackStats { + // A unique id within a run, assigned consecutively by order of creation. -1 + // if no nodes have been created yet, or num_processed_nodes if + // search is over. See SCIPgetCurrentNode(). + int64_t current_node_id = 0; + + // The number of processed nodes in the current run (i.e. does not include + // nodes before a restart), including the focus node. See SCIPgetNNodes(). + int64_t num_processed_nodes = 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; +}; + +// Interface to the callback context and underlying problem. Supports adding +// cuts and lazy constraints, and setting bounds. Prefer to use this context to +// query information instead of a raw SCIP pointer. Passed by value. +// TODO(user): Add support for branching. +class GScipConstraintHandlerContext { + public: + // Construct the context for the given handler. Following SCIP convention, if + // SCIP_SOL is nullptr, then the current solution from the LP is used. + GScipConstraintHandlerContext(GScip* gscip, const GScipCallbackStats* stats, + SCIP_CONSHDLR* current_handler, + SCIP_SOL* current_solution) + : gscip_(gscip), + stats_(stats), + current_handler_(current_handler), + current_solution_(current_solution) {} + + GScip* gscip() { return gscip_; } + + // Returns the current solution value of a variable. This may be for a given + // solution (e.g. in CONS_SEPASOL) or the current LP/pseudosolution (e.g. in + // CONS_SEPALP). Equivalent to calling SCIPgetSolVal. + double VariableValue(SCIP_VAR* variable) const; + + // Adds a cut (row) to the SCIP separation storage. + // + // If this is called and succeeds, the callback result must be the one + // returned or a higher priority result. The result returned is either kCutOff + // (SCIP_CUTOFF) if SCIP determined that the cut results in infeasibility + // based on local bounds, or kSeparated (SCIP_SEPARATED) otherwise. + absl::StatusOr AddCut( + const GScipLinearRange& range, const std::string& name, + const GScipCutOptions& options = GScipCutOptions()); + + // Adds a lazy constraint as a SCIP linear constraint. This is similar to + // adding it as a row (and it would be valid to add a lazy constraint with + // AddCut and proper options), but it is treated as a higher-level object and + // may affect other portions of SCIP such as propagation. This is a thin + // wrapper on GScip::AddLinearConstraint() with different defaults. + // + // If this is called and succeeds, the callback result must be kConsAdded + // (equivalent to SCIP_CONSADDED) or a higher priority result. + absl::Status AddLazyLinearConstraint( + const GScipLinearRange& range, const std::string& name, + const GScipLazyConstraintOptions& options = GScipLazyConstraintOptions()); + + // The functions below set variable bounds. If they are used to cut off a + // solution, then the callback result must be kReducedDomain + // (SCIP_REDUCEDDOM) or a higher priority result. + + absl::Status SetLocalVarLb(SCIP_VAR* var, double value); + absl::Status SetLocalVarUb(SCIP_VAR* var, double value); + absl::Status SetGlobalVarLb(SCIP_VAR* var, double value); + absl::Status SetGlobalVarUb(SCIP_VAR* var, double value); + + double LocalVarLb(SCIP_VAR* var) const; + double LocalVarUb(SCIP_VAR* var) const; + double GlobalVarLb(SCIP_VAR* var) const; + double GlobalVarUb(SCIP_VAR* var) const; + + const GScipCallbackStats& stats() const { return *stats_; } + + private: + GScip* gscip_; + const GScipCallbackStats* stats_; + SCIP_CONSHDLR* current_handler_; + SCIP_SOL* current_solution_; +}; + +// Constraint handler class. To implement a constraint handler, the user can +// inherit this class and implement the desired callback functions. The +// templated ConstraintData is the equivalent of SCIP's SCIP_CONSHDLRDATA, and +// can hold the data needed for the constraint. To then use it, the function +// Register must be called once, and AddCallbackConstraint must be called for +// each constraint to be added in this constraint handler. +// +// There is a one-to-one mapping between relevant SCIP callback functions and +// the functions in this class; see SCIP documentation for which types of +// callbacks to use. Make sure to follow SCIP's rules (e.g. if implementing +// enforcement, all enforcement and check callbacks must be implemented). +// +// For examples of usage, see gscip_constraint_handler_test.cc. +// +// Implementation details: +// +// * Default implementations: All callback functions have a default +// implementation that returns "did not run" or "feasible" accordingly. For +// rounding lock, the default implementation locks both directions. +// +// * Status errors: If the user returns an absl::Status error, then the solve is +// interrupted via SCIPinterruptSolve(), and the status error is ultimately +// returned by GScip::Solve() after SCIP completes the interruption. The +// callback function returns SCIP_OKAY to SCIP except for internal errors. We +// try to avoid returning SCIP_ERROR in the middle of a callback since SCIP +// might not stay in a fully clean state (e.g. calling SCIPfree might hit an +// assert). +// +// * Constraint priority: SCIP informs the callback which subset of constraints +// are more likely to be violated. The callback is called on those constraints +// first, and if the highest priority result is kDidNotFind, kDidNotRun, or +// kFeasible, it is called for the remaining ones. +// +// Supported SCIP callback functions: +// * SCIP_DECL_CONSENFOLP +// * SCIP_DECL_CONSENFOPS +// * SCIP_DECL_CONSCHECK +// * SCIP_DECL_CONSLOCK +// * SCIP_DECL_CONSSEPALP +// * SCIP_DECL_CONSSEPASOL +// +// Used, but not customizable: +// * SCIP_DECL_CONSFREE +// * SCIP_DECL_CONSINIT +// * SCIP_DECL_CONSDELETE +template +class GScipConstraintHandler { + public: + // Constructs a constraint handler that will be registered using the given + // properties. It is recommended to set priorities and frequencies manually in + // properties. + explicit GScipConstraintHandler( + const GScipConstraintHandlerProperties& properties) + : properties_(properties) {} + + virtual ~GScipConstraintHandler() = default; + + const GScipConstraintHandlerProperties& properties() const { + return properties_; + } + + // Registers this constraint handler with GScip. If the handler has already + // been registered, returns an error. + absl::Status Register(GScip* gscip); + + // 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); + + // 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 + // branching, reducing variable domains, or separating the solution with a + // cutting plane. If properties_.enforcement_priority < 0, then this only acts + // on integer solutions. + // + // SCIP CONSENFOLP callback arguments: + // * solution_infeasible: solinfeasible in SCIP, indicates if the solution was + // already declared infeasible by a constraint handler. + // + // It is the user's responsibility to return a valid result for CONSENFOLP; + // see SCIP's documentation (e.g. type_cons.h). + virtual absl::StatusOr EnforceLp( + GScipConstraintHandlerContext context, + const ConstraintData& constraint_data, bool solution_infeasible); + + // Callback function called at SCIP's CONSENFOPS. Must check if a + // pseudosolution is feasible, and if not, resolve the infeasibility if + // possible by branching, reducing variable domains, or adding an additional + // constraint. Separating with a cutting plane is not possible since there is + // no corresponding LP (i.e. kSeparated cannot be returned). If + // properties_.enforcement_priority < 0, then this only acts on integer + // solutions. + // + // SCIP CONSENFOPS callback arguments: + // * solution_infeasible: solinfeasible in SCIP, indicates if the solution was + // already declared infeasible by a constraint handler. + // * objective_infeasible: objinfeasible in SCIP, indicates if the solution is + // infeasible due to violating objective bound. + // + // It is the user's responsibility to return a valid result for CONSENFOPS; + // see SCIP's documentation (e.g. type_cons.h). + virtual absl::StatusOr EnforcePseudoSolution( + GScipConstraintHandlerContext context, + const ConstraintData& constraint_data, bool solution_infeasible, + bool objective_infeasible); + + // Callback function called at SCIP's CONSCHECK. Must return true if the + // current solution stored in the context satisfies all constraints of the + // constraint handler, or false otherwise. If + // properties_.feasibility_check_priority < 0, then this only acts on integer + // solutions. + // + // SCIP CONSCHECK callback arguments: + // * check_integrality: checkintegrality in SCIP, indicates if integrality + // must be checked. Used to avoid redundant checks in cases where + // integrality is already checked or implicit. + // * check_lp_rows: checklprows in SCIP, indicates if the constraints + // represented by rows in the current LP must be checked. Used to avoid + // redundant checks in cases where row feasibility is already checked or + // implicit. + // * print_reason: printreason in SCIP, indicates if the reason for the + // violation should be printed. + // * check_completely: completely in SCIP, indicates if all violations should + // be checked. + virtual absl::StatusOr CheckIsFeasible( + GScipConstraintHandlerContext context, + const ConstraintData& constraint_data, bool check_integrality, + bool check_lp_rows, bool print_reason, bool check_completely); + + // Callback function called at SCIP's CONSLOCK. Must return, for each + // variable, whether the constraint may be violated by decreasing and/or + // increasing the variable value. It is always safe to claim that both + // directions can violate the constraint, which is the default implementation, + // but it may affect SCIP's capabilities. + // + // SCIP CONSLOCK callback arguments: + // * lock_type_is_model: if locktype == SCIP_LOCKTYPE_MODEL in SCIP. If true, + // this callback is called for model constraints, otherwise it is called for + // conflict constraints. + // + // It is the user's responsibility to return a valid result for CONSLOCK; see + // SCIP's documentation (e.g. type_cons.h). + virtual std::vector> RoundingLock( + GScip* gscip, const ConstraintData& constraint_data, + bool lock_type_is_model); + + // Callback function called at SCIP's CONSSEPALP. Separates all constraints of + // the constraint handler for LP solutions. + // + // It is the user's responsibility to return a valid result for CONSSEPALP; + // see SCIP's documentation (e.g. type_cons.h). + virtual absl::StatusOr SeparateLp( + GScipConstraintHandlerContext context, + const ConstraintData& constraint_data); + + // Callback function called at SCIP's CONSSEPASOL. Separates all constraints + // of the constraint handler for solutions that do not come from LP (e.g. + // relaxators and primal heuristics). + // + // It is the user's responsibility to return a valid result for CONSSEPASOL; + // see SCIP's documentation (e.g. type_cons.h). + virtual absl::StatusOr SeparateSolution( + GScipConstraintHandlerContext context, + const ConstraintData& constraint_data); + + // The functions below wrap each callback function to manage status. + + GScipCallbackResult CallEnforceLp(GScipConstraintHandlerContext context, + const ConstraintData& constraint_data, + bool solution_infeasible); + + GScipCallbackResult CallEnforcePseudoSolution( + GScipConstraintHandlerContext context, + const ConstraintData& constraint_data, bool solution_infeasible, + bool objective_infeasible); + + GScipCallbackResult CallCheckIsFeasible(GScipConstraintHandlerContext context, + const ConstraintData& constraint_data, + bool check_integrality, + bool check_lp_rows, bool print_reason, + bool check_completely); + + GScipCallbackResult CallSeparateLp(GScipConstraintHandlerContext context, + const ConstraintData& constraint_data); + + GScipCallbackResult CallSeparateSolution( + GScipConstraintHandlerContext context, + const ConstraintData& constraint_data); + + private: + // If the result status is not OK, stores the status for GScip to later return + // to the user, and interrupts the solve. Otherwise, returns the result + // itself. + GScipCallbackResult HandleCallbackStatus( + absl::StatusOr result, + GScipConstraintHandlerContext context, + GScipCallbackResult default_callback_result); + + GScipConstraintHandlerProperties properties_; +}; + +namespace internal { + +// These classes implement a void pointer version of GScipConstraintHandler that +// performs a conversion to the more user-friendly templated implementation with +// ConstraintData. This parent class is needed as an untemplated version that +// SCIP_ConshdlrData can store. +class UntypedGScipConstraintHandler : public GScipConstraintHandler { + public: + explicit UntypedGScipConstraintHandler( + const GScipConstraintHandlerProperties& properties) + : GScipConstraintHandler(properties) {} +}; + +template +class UntypedGScipConstraintHandlerImpl : public UntypedGScipConstraintHandler { + public: + explicit UntypedGScipConstraintHandlerImpl( + GScipConstraintHandler* constraint_handler) + : UntypedGScipConstraintHandler(constraint_handler->properties()), + actual_handler_(constraint_handler) {} + + absl::StatusOr EnforceLp( + GScipConstraintHandlerContext context, void* const& constraint_data, + bool solution_infeasible) override { + return actual_handler_->EnforceLp( + context, *static_cast(constraint_data), + solution_infeasible); + } + + absl::StatusOr EnforcePseudoSolution( + GScipConstraintHandlerContext context, void* const& constraint_data, + bool solution_infeasible, bool objective_infeasible) override { + return actual_handler_->EnforcePseudoSolution( + context, *static_cast(constraint_data), + solution_infeasible, objective_infeasible); + } + + absl::StatusOr CheckIsFeasible(GScipConstraintHandlerContext context, + void* const& constraint_data, + bool check_integrality, + bool check_lp_rows, bool print_reason, + bool completely) override { + return actual_handler_->CheckIsFeasible( + context, *static_cast(constraint_data), + check_integrality, check_lp_rows, print_reason, completely); + } + + std::vector> RoundingLock( + GScip* gscip, void* const& constraint_data, + bool lock_type_is_model) override { + return actual_handler_->RoundingLock( + gscip, *static_cast(constraint_data), + lock_type_is_model); + } + + absl::StatusOr SeparateLp( + GScipConstraintHandlerContext context, + void* const& constraint_data) override { + return actual_handler_->SeparateLp( + context, *static_cast(constraint_data)); + } + + absl::StatusOr SeparateSolution( + GScipConstraintHandlerContext context, + void* const& constraint_data) override { + return actual_handler_->SeparateSolution( + context, *static_cast(constraint_data)); + } + + private: + GScipConstraintHandler* actual_handler_; +}; + +absl::Status RegisterConstraintHandler( + GScip* gscip, + std::unique_ptr constraint_handler); + +absl::Status AddCallbackConstraint(GScip* gscip, std::string_view handler_name, + std::string_view constraint_name, + void* constraint_data, + const GScipConstraintOptions& options); + +} // namespace internal + +// Template implementations. + +template +absl::Status GScipConstraintHandler::Register(GScip* gscip) { + return internal::RegisterConstraintHandler( + gscip, + std::make_unique< + internal::UntypedGScipConstraintHandlerImpl>(this)); +} + +template +absl::Status GScipConstraintHandler::AddCallbackConstraint( + GScip* gscip, std::string_view constraint_name, + const ConstraintData* constraint_data, + const GScipConstraintOptions& options) { + return internal::AddCallbackConstraint( + gscip, properties().name, constraint_name, + static_cast(const_cast(constraint_data)), + options); +} + +// Default callback implementations. + +template +absl::StatusOr +GScipConstraintHandler::EnforceLp( + GScipConstraintHandlerContext /*context*/, + const ConstraintData& /*constraint_data*/, bool /*solution_infeasible*/) { + return GScipCallbackResult::kFeasible; +} + +template +absl::StatusOr +GScipConstraintHandler::EnforcePseudoSolution( + GScipConstraintHandlerContext /*context*/, + const ConstraintData& /*constraint_data*/, bool /*solution_infeasible*/, + bool /*objective_infeasible*/) { + return GScipCallbackResult::kFeasible; +} + +template +absl::StatusOr GScipConstraintHandler::CheckIsFeasible( + GScipConstraintHandlerContext /*context*/, + const ConstraintData& /*constraint_data*/, bool /*check_integrality*/, + bool /*check_lp_rows*/, bool /*print_reason*/, bool /*check_completely*/) { + return true; +} + +template +std::vector> +GScipConstraintHandler::RoundingLock( + GScip* gscip, const ConstraintData& /*constraint_data*/, + bool /*lock_type_is_model*/) { + std::vector> result; + for (SCIP_VAR* var : gscip->variables()) { + result.push_back({var, RoundingLockDirection::kBoth}); + } + return result; +} + +template +absl::StatusOr +GScipConstraintHandler::SeparateLp( + GScipConstraintHandlerContext /*context*/, + const ConstraintData& /*constraint_data*/) { + return GScipCallbackResult::kDidNotRun; +} + +template +absl::StatusOr +GScipConstraintHandler::SeparateSolution( + GScipConstraintHandlerContext /*context*/, + const ConstraintData& /*constraint_data*/) { + return GScipCallbackResult::kDidNotRun; +} + +// Internal functions to handle status. + +template +GScipCallbackResult +GScipConstraintHandler::HandleCallbackStatus( + const absl::StatusOr result, + GScipConstraintHandlerContext context, + const GScipCallbackResult default_callback_result) { + if (!result.ok()) { + context.gscip()->InterruptSolveFromCallback(result.status()); + return default_callback_result; + } + return result.value(); +} + +template +GScipCallbackResult GScipConstraintHandler::CallEnforceLp( + GScipConstraintHandlerContext context, + const ConstraintData& constraint_data, bool solution_infeasible) { + return HandleCallbackStatus( + EnforceLp(context, constraint_data, solution_infeasible), context, + GScipCallbackResult::kFeasible); +} + +template +GScipCallbackResult +GScipConstraintHandler::CallEnforcePseudoSolution( + GScipConstraintHandlerContext context, + const ConstraintData& constraint_data, bool solution_infeasible, + bool objective_infeasible) { + return HandleCallbackStatus( + EnforcePseudoSolution(context, constraint_data, solution_infeasible, + objective_infeasible), + context, GScipCallbackResult::kFeasible); +} + +template +GScipCallbackResult GScipConstraintHandler::CallCheckIsFeasible( + GScipConstraintHandlerContext context, + const ConstraintData& constraint_data, bool check_integrality, + bool check_lp_rows, bool print_reason, bool check_completely) { + const absl::StatusOr result = + CheckIsFeasible(context, constraint_data, check_integrality, + check_lp_rows, print_reason, check_completely); + if (!result.ok()) { + return HandleCallbackStatus(result.status(), context, + GScipCallbackResult::kFeasible); + } else { + return HandleCallbackStatus(*result ? GScipCallbackResult::kFeasible + : GScipCallbackResult::kInfeasible, + context, GScipCallbackResult::kFeasible); + } +} + +template +GScipCallbackResult GScipConstraintHandler::CallSeparateLp( + GScipConstraintHandlerContext context, + const ConstraintData& constraint_data) { + return HandleCallbackStatus(SeparateLp(context, constraint_data), context, + GScipCallbackResult::kDidNotRun); +} + +template +GScipCallbackResult +GScipConstraintHandler::CallSeparateSolution( + GScipConstraintHandlerContext context, + const ConstraintData& constraint_data) { + return HandleCallbackStatus(SeparateSolution(context, constraint_data), + context, GScipCallbackResult::kDidNotRun); +} + +} // namespace operations_research + +#endif // OR_TOOLS_GSCIP_GSCIP_CONSTRAINT_HANDLER_H_ From 753a7ab8740b750186332c7b2ecc671e470392b3 Mon Sep 17 00:00:00 2001 From: Corentin Le Molgat Date: Mon, 9 Oct 2023 14:30:15 +0200 Subject: [PATCH 2/5] update constraint_solver from google3 --- .../constraint_solver/python/pywrapcp_test.py | 2 +- ortools/constraint_solver/routing.h | 7 + ortools/constraint_solver/routing_search.cc | 358 +++++++++++------- ortools/constraint_solver/routing_search.h | 27 +- 4 files changed, 253 insertions(+), 141 deletions(-) diff --git a/ortools/constraint_solver/python/pywrapcp_test.py b/ortools/constraint_solver/python/pywrapcp_test.py index b553e803f58..051afe3b397 100755 --- a/ortools/constraint_solver/python/pywrapcp_test.py +++ b/ortools/constraint_solver/python/pywrapcp_test.py @@ -877,7 +877,7 @@ class SumFilter(pywrapcp.IntVarLocalSearchFilter): """Filter to speed up LS computation.""" def __init__(self, int_vars): - pywrapcp.IntVarLocalSearchFilter.__init__(self, int_vars) + super.__init__(int_vars) self.__sum = 0 def OnSynchronize(self, delta): diff --git a/ortools/constraint_solver/routing.h b/ortools/constraint_solver/routing.h index 7a650cd6d04..e04eb2c21b2 100644 --- a/ortools/constraint_solver/routing.h +++ b/ortools/constraint_solver/routing.h @@ -907,6 +907,13 @@ class RoutingModel { const std::vector >& GetDeliveryIndexPairs(int64_t node_index) const; // clang-format on + /// Returns whether the node is a pickup (resp. delivery). + bool IsPickup(int64_t node_index) const { + return !GetPickupIndexPairs(node_index).empty(); + } + bool IsDelivery(int64_t node_index) const { + return !GetDeliveryIndexPairs(node_index).empty(); + } /// Sets the Pickup and delivery policy of all vehicles. It is equivalent to /// calling SetPickupAndDeliveryPolicyOfVehicle on all vehicles. diff --git a/ortools/constraint_solver/routing_search.cc b/ortools/constraint_solver/routing_search.cc index f09e4078f5b..3d5f472c66b 100644 --- a/ortools/constraint_solver/routing_search.cc +++ b/ortools/constraint_solver/routing_search.cc @@ -24,11 +24,9 @@ #include #include #include -#include #include #include #include -#include #include #include #include @@ -36,6 +34,7 @@ #include #include +#include "absl/algorithm/container.h" #include "absl/base/attributes.h" #include "absl/container/flat_hash_map.h" #include "absl/container/flat_hash_set.h" @@ -532,46 +531,67 @@ CheapestInsertionFilteredHeuristic::CheapestInsertionFilteredHeuristic( evaluator_(std::move(evaluator)), penalty_evaluator_(std::move(penalty_evaluator)) {} +namespace { +void ProcessVehicleStartEndCosts( + const RoutingModel& model, int node, + const std::function& process_cost, + const Bitset64& vehicle_set, bool ignore_start = false, + bool ignore_end = false) { + const auto start_end_cost = [&model, ignore_start, ignore_end](int64_t node, + int v) { + const int64_t start_cost = + ignore_start ? 0 : model.GetArcCostForVehicle(model.Start(v), node, v); + const int64_t end_cost = + ignore_end ? 0 : model.GetArcCostForVehicle(model.End(v), node, v); + return CapAdd(start_cost, end_cost); + }; + + // Iterating over an IntVar domain is faster than calling Contains. + // Therefore we iterate on 'vehicles' only if it's smaller than the domain + // size of the VehicleVar. + const IntVar* const vehicle_var = model.VehicleVar(node); + if (vehicle_var->Size() < vehicle_set.size()) { + std::unique_ptr it(vehicle_var->MakeDomainIterator(false)); + for (const int v : InitAndGetValues(it.get())) { + if (v < 0 || !vehicle_set[v]) { + continue; + } + process_cost(start_end_cost(node, v), v); + } + } else { + for (const int v : vehicle_set) { + if (!vehicle_var->Contains(v)) continue; + process_cost(start_end_cost(node, v), v); + } + } +} +} // namespace + std::vector> CheapestInsertionFilteredHeuristic::ComputeStartEndDistanceForVehicles( const std::vector& vehicles) { // TODO(user): consider checking search limits. - const absl::flat_hash_set vehicle_set(vehicles.begin(), vehicles.end()); + const RoutingModel& model = *this->model(); std::vector> start_end_distances_per_node( - model()->Size()); + model.Size()); - for (int node = 0; node < model()->Size(); node++) { + Bitset64 vehicle_set(model.vehicles()); + for (int v : vehicles) vehicle_set.Set(v); + + for (int node = 0; node < model.Size(); node++) { if (Contains(node)) continue; std::vector& start_end_distances = start_end_distances_per_node[node]; + start_end_distances.reserve( + std::min(model.VehicleVar(node)->Size(), vehicles.size())); + + ProcessVehicleStartEndCosts( + model, node, + [&start_end_distances](int64_t dist, int v) { + start_end_distances.push_back({dist, v}); + }, + vehicle_set); - const auto add_distance = [this, node, &start_end_distances](int vehicle) { - const int64_t start = model()->Start(vehicle); - const int64_t end = model()->End(vehicle); - // We compute the distance of node to the start/end nodes of the route. - const int64_t distance = - CapAdd(model()->GetArcCostForVehicle(start, node, vehicle), - model()->GetArcCostForVehicle(node, end, vehicle)); - start_end_distances.push_back({distance, vehicle}); - }; - // Iterating over an IntVar domain is faster than calling Contains. - // Therefore we iterate on 'vehicles' only if it's smaller than the domain - // size of the VehicleVar. - const IntVar* const vehicle_var = model()->VehicleVar(node); - if (vehicle_var->Size() < vehicles.size()) { - std::unique_ptr it( - vehicle_var->MakeDomainIterator(false)); - for (const int64_t vehicle : InitAndGetValues(it.get())) { - if (vehicle < 0 || !vehicle_set.contains(vehicle)) continue; - add_distance(vehicle); - } - } else { - start_end_distances.reserve(vehicles.size()); - for (const int vehicle : vehicles) { - if (!vehicle_var->Contains(vehicle)) continue; - add_distance(vehicle); - } - } // Sort the distances for the node to all start/ends of available vehicles // in decreasing order. std::sort(start_end_distances.begin(), start_end_distances.end(), @@ -582,7 +602,7 @@ CheapestInsertionFilteredHeuristic::ComputeStartEndDistanceForVehicles( return start_end_distances_per_node; } -void CheapestInsertionFilteredHeuristic::AddSeedToQueue( +void CheapestInsertionFilteredHeuristic::AddSeedNodeToQueue( int node, std::vector* start_end_distances, SeedQueue* sq) { if (start_end_distances->empty()) { return; @@ -596,7 +616,7 @@ void CheapestInsertionFilteredHeuristic::AddSeedToQueue( const uint64_t num_allowed_vehicles = model()->VehicleVar(node)->Size(); const int64_t neg_penalty = CapOpp(model()->UnperformedPenalty(node)); sq->priority_queue.push( - {num_allowed_vehicles, neg_penalty, start_end_value, node}); + {num_allowed_vehicles, neg_penalty, start_end_value, true, node}); start_end_distances->pop_back(); } @@ -608,7 +628,7 @@ void CheapestInsertionFilteredHeuristic::InitializeSeedQueue( for (int node = 0; node < num_nodes; node++) { if (Contains(node)) continue; - AddSeedToQueue(node, &start_end_distances_per_node->at(node), sq); + AddSeedNodeToQueue(node, &start_end_distances_per_node->at(node), sq); } } @@ -791,8 +811,8 @@ bool GlobalCheapestInsertionFilteredHeuristic::BuildSolutionInternal() { } std::map> nodes_by_bucket; for (int node = 0; node < model()->Size(); ++node) { - if (!Contains(node) && model()->GetPickupIndexPairs(node).empty() && - model()->GetDeliveryIndexPairs(node).empty()) { + if (!Contains(node) && !model()->IsPickup(node) && + !model()->IsDelivery(node)) { nodes_by_bucket[GetBucketOfNode(node)].push_back(node); } } @@ -1460,7 +1480,8 @@ int GlobalCheapestInsertionFilteredHeuristic::InsertSeedNode( while (!priority_queue.empty()) { if (StopSearch()) return -1; const Seed& seed = priority_queue.top(); - const int seed_node = seed.node; + const int seed_node = seed.index; + DCHECK(seed.is_node_index); const int seed_vehicle = seed.start_end_value.vehicle; priority_queue.pop(); @@ -1489,7 +1510,7 @@ int GlobalCheapestInsertionFilteredHeuristic::InsertSeedNode( // Either the vehicle is already used, or the Commit() wasn't successful. // In both cases, we insert the next StartEndValue from // start_end_distances_per_node[seed_node] in the priority queue. - AddSeedToQueue(seed_node, &other_start_end_values, sq); + AddSeedNodeToQueue(seed_node, &other_start_end_values, sq); } // No seed node was inserted. return -1; @@ -2222,19 +2243,140 @@ LocalCheapestInsertionFilteredHeuristic:: : CheapestInsertionFilteredHeuristic(model, std::move(stop_search), std::move(evaluator), nullptr, filter_manager), - update_start_end_distances_per_node_(true), pair_insertion_strategy_(pair_insertion_strategy), bin_capacities_(bin_capacities) {} void LocalCheapestInsertionFilteredHeuristic::Initialize() { - // Avoid recomputing if used in a local search operator. - if (update_start_end_distances_per_node_) { - update_start_end_distances_per_node_ = false; - std::vector all_vehicles(model()->vehicles()); - std::iota(std::begin(all_vehicles), std::end(all_vehicles), 0); - start_end_distances_per_node_ = - ComputeStartEndDistanceForVehicles(all_vehicles); + // NOTE(user): Keeping the code in a separate function as opposed to + // inlining here, to allow for future additions to this function. + ComputeInsertionOrder(); +} + +namespace { +// Returns the opposite of the maximum cost between all pickup/delivery nodes of +// the given pair from their "closest" vehicle. +// For a given pickup/delivery, the cost from the closest vehicle is defined as +// min_{v in vehicles}(ArcCost(start[v]->pickup) + ArcCost(delivery->end[v])). +int64_t GetNegMaxDistanceFromVehicles(const RoutingModel& model, + int pair_index) { + const RoutingModel::IndexPair& pickup_deliveries = + model.GetPickupAndDeliveryPairs()[pair_index]; + + Bitset64 vehicle_set(model.vehicles()); + for (int v = 0; v < model.vehicles(); ++v) vehicle_set.Set(v); + + // Precompute the cost from vehicle starts to every pickup in the pair. + std::vector> pickup_costs(model.Size()); + for (int64_t pickup : pickup_deliveries.first) { + std::vector& cost_from_start = pickup_costs[pickup]; + cost_from_start.resize(model.vehicles(), -1); + + ProcessVehicleStartEndCosts( + model, pickup, + [&cost_from_start](int64_t cost, int v) { cost_from_start[v] = cost; }, + vehicle_set, /*ignore_start=*/false, /*ignore_end=*/true); + } + + // Precompute the cost from every delivery in the pair to vehicle ends. + std::vector> delivery_costs(model.Size()); + for (int64_t delivery : pickup_deliveries.second) { + std::vector& cost_to_end = delivery_costs[delivery]; + cost_to_end.resize(model.vehicles(), -1); + + ProcessVehicleStartEndCosts( + model, delivery, + [&cost_to_end](int64_t cost, int v) { cost_to_end[v] = cost; }, + vehicle_set, /*ignore_start=*/true, /*ignore_end=*/false); + } + + int64_t max_pair_distance = 0; + for (int64_t pickup : pickup_deliveries.first) { + const std::vector& cost_from_start = pickup_costs[pickup]; + for (int64_t delivery : pickup_deliveries.second) { + const std::vector& cost_to_end = delivery_costs[delivery]; + int64_t closest_vehicle_distance = std::numeric_limits::max(); + for (int v = 0; v < model.vehicles(); v++) { + if (cost_from_start[v] < 0 || cost_to_end[v] < 0) { + // Vehicle not in the pickup and/or delivery's vehicle var domain. + continue; + } + closest_vehicle_distance = + std::min(closest_vehicle_distance, + CapAdd(cost_from_start[v], cost_to_end[v])); + } + max_pair_distance = std::max(max_pair_distance, closest_vehicle_distance); + } } + return CapOpp(max_pair_distance); +} +} // namespace + +void LocalCheapestInsertionFilteredHeuristic::ComputeInsertionOrder() { + if (!insertion_order_.empty()) return; + + // We consider pairs and single nodes simultaneously, to make sure the most + // critical ones (fewer allowed vehicles and high penalties) get inserted + // first. + // TODO(user): Explore other metrics to evaluate criticality, more directly + // mixing penalties and allowed vehicles (such as a ratio between the two). + + const RoutingModel& model = *this->model(); + insertion_order_.reserve(model.Size() + + model.GetPickupAndDeliveryPairs().size()); + + // Iterating on pickup and delivery pairs + const RoutingModel::IndexPairs& index_pairs = + model.GetPickupAndDeliveryPairs(); + + for (int pair_index = 0; pair_index < index_pairs.size(); ++pair_index) { + uint64_t num_allowed_vehicles = std::numeric_limits::max(); + int64_t pickup_penalty = 0; + for (int64_t pickup : index_pairs[pair_index].first) { + num_allowed_vehicles = + std::min(num_allowed_vehicles, model.VehicleVar(pickup)->Size()); + pickup_penalty = + std::max(pickup_penalty, model.UnperformedPenalty(pickup)); + } + int64_t delivery_penalty = 0; + for (int64_t delivery : index_pairs[pair_index].second) { + num_allowed_vehicles = + std::min(num_allowed_vehicles, model.VehicleVar(delivery)->Size()); + delivery_penalty = + std::max(delivery_penalty, model.UnperformedPenalty(delivery)); + } + insertion_order_.push_back( + {num_allowed_vehicles, + CapOpp(CapAdd(pickup_penalty, delivery_penalty)), + {GetNegMaxDistanceFromVehicles(model, pair_index), 0}, + false, + pair_index}); + } + + Bitset64 vehicle_set(model.vehicles()); + for (int v = 0; v < model.vehicles(); ++v) vehicle_set.Set(v); + + for (int node = 0; node < model.Size(); ++node) { + if (model.IsStart(node) || model.IsEnd(node)) continue; + + int64_t min_distance = std::numeric_limits::max(); + ProcessVehicleStartEndCosts( + model, node, + [&min_distance](int64_t dist, int) { + min_distance = std::min(min_distance, dist); + }, + vehicle_set); + + const uint64_t num_allowed_vehicles = model.VehicleVar(node)->Size(); + const int64_t neg_penalty = CapOpp(model.UnperformedPenalty(node)); + insertion_order_.push_back({num_allowed_vehicles, + neg_penalty, + {CapOpp(min_distance), 0}, + true, + node}); + } + + absl::c_sort(insertion_order_, std::greater()); + absl::c_reverse(insertion_order_); } bool LocalCheapestInsertionFilteredHeuristic::InsertPair( @@ -2450,56 +2592,13 @@ bool LocalCheapestInsertionFilteredHeuristic::BuildSolutionInternal() { // Marking if we've tried inserting a node. visited_.assign(model()->Size(), false); - // Iterating on pickup and delivery pairs + // Multitour needs to know if a node is a pickup, delivery, or single. + // TODO(user): node_is_[pickup/delivery] is not actually required as we + // have access to this information through + // model()->IsPickup/Delivery(), so use this as callback wherever necessary + // instead of these 2 vectors. const RoutingModel::IndexPairs& index_pairs = model()->GetPickupAndDeliveryPairs(); - // Sort pairs according to number of possible vehicles. - struct PairDomainSize { - uint64_t num_allowed_vehicles; - int64_t neg_penalty; - int pair_index; - - bool operator<(const PairDomainSize& other) const { - return std::tie(num_allowed_vehicles, neg_penalty, pair_index) < - std::tie(other.num_allowed_vehicles, other.neg_penalty, - other.pair_index); - } - }; - std::vector pair_domain_sizes; - absl::flat_hash_set pair_nodes; - for (int pair_index = 0; pair_index < index_pairs.size(); ++pair_index) { - bool pickup_is_contained = false; - uint64_t domain_size = std::numeric_limits::max(); - int64_t pickup_penalty = 0; - for (int64_t pickup : index_pairs[pair_index].first) { - domain_size = std::min(domain_size, model()->VehicleVar(pickup)->Size()); - pickup_penalty = - std::max(pickup_penalty, model()->UnperformedPenalty(pickup)); - pickup_is_contained |= Contains(pickup); - } - bool delivery_is_contained = false; - int64_t delivery_penalty = 0; - for (int64_t delivery : index_pairs[pair_index].second) { - domain_size = - std::min(domain_size, model()->VehicleVar(delivery)->Size()); - delivery_penalty = - std::max(delivery_penalty, model()->UnperformedPenalty(delivery)); - delivery_is_contained |= Contains(delivery); - } - if (pickup_is_contained && delivery_is_contained) { - SetIndexPairVisited(index_pairs[pair_index]); - } else if (!pickup_is_contained && !delivery_is_contained) { - pair_domain_sizes.push_back( - {domain_size, CapOpp(CapAdd(pickup_penalty, delivery_penalty)), - pair_index}); - pair_nodes.insert(index_pairs[pair_index].first.begin(), - index_pairs[pair_index].first.end()); - pair_nodes.insert(index_pairs[pair_index].second.begin(), - index_pairs[pair_index].second.end()); - } - } - std::sort(pair_domain_sizes.begin(), pair_domain_sizes.end()); - // Multitour needs to know if a node is a pickup, delivery, or single. std::vector node_is_pickup, node_is_delivery; if (pair_insertion_strategy_ == RoutingSearchParameters::BEST_PICKUP_DELIVERY_PAIR_MULTITOUR) { @@ -2526,37 +2625,34 @@ bool LocalCheapestInsertionFilteredHeuristic::BuildSolutionInternal() { } } - // TODO(user): A priority queue is not actually needed here as we only use - // the "Seed" values to determine the static order of insertion of the nodes. - // This can be done with a simple sorted container, which could also contain - // the data for pairs in the same place to determine the overall order of - // insertion. - SeedQueue sq(/*prioritize_farthest_nodes=*/true); - InitializeSeedQueue(&start_end_distances_per_node_, &sq); - auto& node_queue = sq.priority_queue; + std::vector ignore_pair_index(index_pairs.size(), false); + for (int pair_index = 0; pair_index < index_pairs.size(); ++pair_index) { + bool pickup_contained = false; + for (int64_t pickup : index_pairs[pair_index].first) { + if (Contains(pickup)) { + pickup_contained = true; + break; + } + } + bool delivery_contained = false; + for (int64_t delivery : index_pairs[pair_index].second) { + if (Contains(delivery)) { + delivery_contained = true; + break; + } + } + ignore_pair_index[pair_index] = pickup_contained || delivery_contained; + if (pickup_contained && delivery_contained) { + SetIndexPairVisited(index_pairs[pair_index]); + } + } - // Inserting pairs and single nodes simultaneously, to make sure the most - // critical ones (fewer allowed vehicles and high penalties) get inserted - // first. - // TODO(user): Explore other metrics to evaluate criticality, more directly - // mixing penalties and allowed vehicles (such as a ratio between the two). - int pair_domain_index = 0; - while (pair_domain_index < pair_domain_sizes.size() || !node_queue.empty()) { - bool insert_pair = false; - if (node_queue.empty()) { - insert_pair = true; - } else if (pair_domain_index < pair_domain_sizes.size()) { - const PairDomainSize& pair_domain_size = - pair_domain_sizes[pair_domain_index]; - const Seed& seed_node = node_queue.top(); - insert_pair = - std::tie(pair_domain_size.num_allowed_vehicles, - pair_domain_size.neg_penalty) <= - std::tie(seed_node.num_allowed_vehicles, seed_node.neg_penalty); - } - if (insert_pair) { - const auto index_pair = - index_pairs[pair_domain_sizes[pair_domain_index].pair_index]; + for (const Seed& seed : insertion_order_) { + const int index = seed.index; + if (!seed.is_node_index) { + if (ignore_pair_index[index]) continue; + + const auto& index_pair = index_pairs[index]; switch (pair_insertion_strategy_) { case RoutingSearchParameters::AUTOMATIC: case RoutingSearchParameters::BEST_PICKUP_DELIVERY_PAIR: @@ -2576,23 +2672,21 @@ bool LocalCheapestInsertionFilteredHeuristic::BuildSolutionInternal() { return MakeUnassignedNodesUnperformed() && Evaluate(true).has_value(); } SetIndexPairVisited(index_pair); - pair_domain_index++; } else { - const int node = node_queue.top().node; - node_queue.pop(); - if (Contains(node) || visited_[node] || pair_nodes.contains(node)) { + if (Contains(index) || visited_[index] || model()->IsPickup(index) || + model()->IsDelivery(index)) { continue; } for (const NodeInsertion& insertion : - ComputeEvaluatorSortedPositions(node)) { + ComputeEvaluatorSortedPositions(index)) { if (StopSearch()) { return MakeUnassignedNodesUnperformed() && Evaluate(true).has_value(); } - InsertBetween(node, insertion.insert_after, + InsertBetween(index, insertion.insert_after, Value(insertion.insert_after), insertion.vehicle); if (Evaluate(/*commit=*/true).has_value()) { if (bin_capacities_) { - bin_capacities_->AddItemToBin(node, insertion.vehicle); + bin_capacities_->AddItemToBin(index, insertion.vehicle); } break; } diff --git a/ortools/constraint_solver/routing_search.h b/ortools/constraint_solver/routing_search.h index 7dfc0bb74a5..6bbf0894636 100644 --- a/ortools/constraint_solver/routing_search.h +++ b/ortools/constraint_solver/routing_search.h @@ -40,6 +40,7 @@ #include "ortools/constraint_solver/constraint_solver.h" #include "ortools/constraint_solver/constraint_solveri.h" #include "ortools/constraint_solver/routing.h" +#include "ortools/constraint_solver/routing_parameters.pb.h" #include "ortools/constraint_solver/routing_utils.h" namespace operations_research { @@ -343,13 +344,16 @@ class CheapestInsertionFilteredHeuristic : public RoutingFilteredHeuristic { uint64_t num_allowed_vehicles; int64_t neg_penalty; StartEndValue start_end_value; - int node; + /// Indicates whether this Seed corresponds to a pair or a single node. + /// If false, the 'index' is the pair_index, otherwise it's the node index. + bool is_node_index = true; + int index; bool operator>(const Seed& other) const { return std::tie(num_allowed_vehicles, neg_penalty, start_end_value, - node) > std::tie(other.num_allowed_vehicles, - other.neg_penalty, other.start_end_value, - other.node); + is_node_index, index) > + std::tie(other.num_allowed_vehicles, other.neg_penalty, + other.start_end_value, other.is_node_index, other.index); } }; // clang-format off @@ -382,10 +386,14 @@ class CheapestInsertionFilteredHeuristic : public RoutingFilteredHeuristic { SeedQueue* sq); // clang-format on + /// Creates and returns a Seed corresponding to the given node/pair_index. + Seed CreateSeed(int index, bool is_pair_index, StartEndValue start_end_value); + /// Adds a Seed corresponding to the given 'node' to sq.priority_queue, based /// on the last entry in its 'start_end_distances' (from which it's deleted). - void AddSeedToQueue(int node, std::vector* start_end_distances, - SeedQueue* sq); + void AddSeedNodeToQueue(int node, + std::vector* start_end_distances, + SeedQueue* sq); /// Inserts 'node' just after 'predecessor', and just before 'successor' on /// the route of 'vehicle', resulting in the following subsequence: @@ -1072,6 +1080,10 @@ class LocalCheapestInsertionFilteredHeuristic void Initialize() override; private: + /// Computes the order of insertion of the node/pairs in the model based on + /// the "Seed" values (number of allowed vehicles, penalty, distance from + /// vehicle start/ends), and stores them in insertion_order_. + void ComputeInsertionOrder(); /// Computes the possible insertion positions of 'node' and sorts them /// according to the current cost evaluator. /// 'node' is a variable index corresponding to a node. @@ -1107,8 +1119,7 @@ class LocalCheapestInsertionFilteredHeuristic // Sets all nodes of pair alternatives as visited. void SetIndexPairVisited(const RoutingModel::IndexPair& index_pair); - bool update_start_end_distances_per_node_; - std::vector> start_end_distances_per_node_; + std::vector insertion_order_; const RoutingSearchParameters::PairInsertionStrategy pair_insertion_strategy_; InsertionSequenceContainer insertion_container_; InsertionSequenceGenerator insertion_generator_; From 69874c9b26c021dcdeeea812397f3e8a09c9987b Mon Sep 17 00:00:00 2001 From: Corentin Le Molgat Date: Mon, 9 Oct 2023 14:23:43 +0200 Subject: [PATCH 3/5] fixup from google3 --- ortools/algorithms/dynamic_partition.h | 3 +-- ortools/constraint_solver/python/constraint_solver.i | 6 ------ ortools/linear_solver/python/model_builder_helper.cc | 9 ++------- ortools/lp_data/sparse_vector.h | 4 ++-- 4 files changed, 5 insertions(+), 17 deletions(-) diff --git a/ortools/algorithms/dynamic_partition.h b/ortools/algorithms/dynamic_partition.h index 1d5072659d0..4b62db70392 100644 --- a/ortools/algorithms/dynamic_partition.h +++ b/ortools/algorithms/dynamic_partition.h @@ -34,8 +34,8 @@ #include #include +#include "absl/log/check.h" #include "absl/types/span.h" -#include "ortools/base/logging.h" namespace operations_research { @@ -405,7 +405,6 @@ inline std::vector> SimpleDynamicPartition::GetParts( starts.clear(); return result; } - } // namespace operations_research #endif // OR_TOOLS_ALGORITHMS_DYNAMIC_PARTITION_H_ diff --git a/ortools/constraint_solver/python/constraint_solver.i b/ortools/constraint_solver/python/constraint_solver.i index 36790357f3d..422d7ae2902 100644 --- a/ortools/constraint_solver/python/constraint_solver.i +++ b/ortools/constraint_solver/python/constraint_solver.i @@ -116,10 +116,6 @@ PY_CONVERT_HELPER_PTR(LocalSearchFilter); PY_CONVERT_HELPER_PTR(LocalSearchFilterManager); PY_CONVERT_HELPER_INTEXPR_AND_INTVAR(); -%{ - -%} - // Actual conversions. This also includes the conversion to std::vector. PY_CONVERT(IntVar); PY_CONVERT(IntExpr); @@ -2198,6 +2194,4 @@ class PyConstraint(Constraint): def DebugString(self): return "PyConstraint" - - } // %pythoncode diff --git a/ortools/linear_solver/python/model_builder_helper.cc b/ortools/linear_solver/python/model_builder_helper.cc index 190707af272..4759581fac2 100644 --- a/ortools/linear_solver/python/model_builder_helper.cc +++ b/ortools/linear_solver/python/model_builder_helper.cc @@ -416,13 +416,8 @@ PYBIND11_MODULE(model_builder_helper, m) { .def("enable_output", &ModelSolverHelper::EnableOutput, arg("output")) .def("has_solution", &ModelSolverHelper::has_solution) .def("has_response", &ModelSolverHelper::has_response) - .def("status", - [](const ModelSolverHelper& solver) { - // TODO(user): - // - Return the true enum when pybind11_protobuf is working. - // - Return the response proto - return static_cast(solver.status()); - }) + .def("response", &ModelSolverHelper::response) + .def("status", &ModelSolverHelper::status) .def("status_string", &ModelSolverHelper::status_string) .def("wall_time", &ModelSolverHelper::wall_time) .def("user_time", &ModelSolverHelper::user_time) diff --git a/ortools/lp_data/sparse_vector.h b/ortools/lp_data/sparse_vector.h index 139a0664495..74caf9a2b9d 100644 --- a/ortools/lp_data/sparse_vector.h +++ b/ortools/lp_data/sparse_vector.h @@ -383,8 +383,8 @@ class SparseVector { EntryIndex capacity_; // Pointers to the first elements of the index and coefficient arrays. - Index* index_; - Fractional* coefficient_; + Index* index_ = nullptr; + Fractional* coefficient_ = nullptr; // This is here to speed up the CheckNoDuplicates() methods and is mutable // so we can perform checks on const argument. From 83c39f0826c904e6d0849730e959d5f99620841c Mon Sep 17 00:00:00 2001 From: Corentin Le Molgat Date: Wed, 6 Sep 2023 11:57:50 +0200 Subject: [PATCH 4/5] math_opt: add mathopt_info.cc sample --- ortools/math_opt/samples/mathopt_info.cc | 50 ++++++++++++++++++++++++ 1 file changed, 50 insertions(+) create mode 100644 ortools/math_opt/samples/mathopt_info.cc diff --git a/ortools/math_opt/samples/mathopt_info.cc b/ortools/math_opt/samples/mathopt_info.cc new file mode 100644 index 00000000000..f67bc00db08 --- /dev/null +++ b/ortools/math_opt/samples/mathopt_info.cc @@ -0,0 +1,50 @@ +// Copyright 2010-2022 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// Check MathOpt available solvers. +#include + +#include "absl/strings/str_cat.h" +#include "absl/strings/str_join.h" +#include "ortools/base/init_google.h" + +#include "ortools/math_opt/cpp/math_opt.h" +#include "ortools/math_opt/parameters.pb.h" +#include "ortools/math_opt/core/solver_interface.h" +//#include "ortools/math_opt/cpp/enums.h" + +namespace { + +struct SolverTypeProtoFormatter { + void operator()( + std::string* const out, + const operations_research::math_opt::SolverTypeProto solver_type) { + out->append(EnumToString(EnumFromProto(solver_type).value())); + } +}; + +} // namespace + +int main(int argc, char* argv[]) { + InitGoogle(argv[0], &argc, &argv, /*remove_flags=*/true); + + std::cout << + absl::StrCat( + "MathOpt is configured to support the following solvers: ", + absl::StrJoin( + operations_research::math_opt::AllSolversRegistry::Instance()->RegisteredSolvers(), + ", ", + SolverTypeProtoFormatter())) << std::endl; + + return 0; +} From 457e927b85992eb44fce31e25a1c1c75fd85bdeb Mon Sep 17 00:00:00 2001 From: Corentin Le Molgat Date: Mon, 9 Oct 2023 16:58:51 +0200 Subject: [PATCH 5/5] indent cleanup --- ortools/gscip/gscip.cc | 2 +- tools/export_to_ipynb.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/ortools/gscip/gscip.cc b/ortools/gscip/gscip.cc index e9a2301cf3f..5749281dff6 100644 --- a/ortools/gscip/gscip.cc +++ b/ortools/gscip/gscip.cc @@ -34,6 +34,7 @@ #include "absl/strings/string_view.h" #include "absl/synchronization/mutex.h" #include "absl/types/span.h" +#include "lpi/lpi.h" #include "ortools/base/logging.h" #include "ortools/base/status_macros.h" #include "ortools/gscip/gscip.pb.h" @@ -57,7 +58,6 @@ #include "scip/scip_solvingstats.h" #include "scip/scip_var.h" #include "scip/scipdefplugins.h" -#include "lpi/lpi.h" #include "scip/type_cons.h" #include "scip/type_paramset.h" #include "scip/type_retcode.h" diff --git a/tools/export_to_ipynb.py b/tools/export_to_ipynb.py index 9900d94e58b..1b4897c9d26 100755 --- a/tools/export_to_ipynb.py +++ b/tools/export_to_ipynb.py @@ -127,7 +127,7 @@ print(f'Removing import {c_block.module}.{c_block.names[0].name}...') # rewrite absl flag import elif (isinstance(c_block, ast.ImportFrom) and c_block.module == 'absl' - and c_block.names[0].name =='flags'): + and c_block.names[0].name == 'flags'): print(f'Rewrite import {c_block.module}.{c_block.names[0].name}...') full_text += 'from ortools.sat.colab import flags\n' # Unwrap __main__ function @@ -156,7 +156,7 @@ print('Appending block...') filtered_lines = lines[s:e] for i, line in enumerate(filtered_lines): - filtered_lines[i] = line.replace('DEFINE_', 'define_') + filtered_lines[i] = line.replace('DEFINE_', 'define_') filtered_lines = list( filter(lambda l: not re.search(r'# \[START .*\]$', l), filtered_lines)) filtered_lines = list(