26#include "absl/cleanup/cleanup.h"
27#include "absl/container/flat_hash_map.h"
28#include "absl/container/flat_hash_set.h"
29#include "absl/log/check.h"
30#include "absl/log/die_if_null.h"
31#include "absl/memory/memory.h"
32#include "absl/status/status.h"
33#include "absl/status/statusor.h"
34#include "absl/strings/str_cat.h"
35#include "absl/strings/str_join.h"
36#include "absl/strings/string_view.h"
37#include "absl/time/clock.h"
38#include "absl/time/time.h"
39#include "absl/types/span.h"
40#include "google/protobuf/repeated_ptr_field.h"
47#include "ortools/gscip/gscip.pb.h"
50#include "ortools/math_opt/callback.pb.h"
57#include "ortools/math_opt/infeasible_subsystem.pb.h"
58#include "ortools/math_opt/model.pb.h"
59#include "ortools/math_opt/model_parameters.pb.h"
60#include "ortools/math_opt/model_update.pb.h"
61#include "ortools/math_opt/parameters.pb.h"
62#include "ortools/math_opt/result.pb.h"
63#include "ortools/math_opt/solution.pb.h"
66#include "ortools/math_opt/sparse_containers.pb.h"
70#include "scip/type_cons.h"
71#include "scip/type_var.h"
78constexpr double kInf = std::numeric_limits<double>::infinity();
88int64_t SafeId(
const VariablesProto& variables,
int index) {
89 if (variables.ids().empty()) {
92 return variables.ids(index);
95const std::string& EmptyString() {
96 static const std::string*
const empty_string =
new std::string;
100const std::string& SafeName(
const VariablesProto& variables,
int index) {
101 if (variables.names().empty()) {
102 return EmptyString();
104 return variables.names(index);
107int64_t SafeId(
const LinearConstraintsProto& linear_constraints,
int index) {
108 if (linear_constraints.ids().empty()) {
111 return linear_constraints.ids(index);
114absl::string_view SafeName(
const LinearConstraintsProto& linear_constraints,
116 if (linear_constraints.names().empty()) {
117 return EmptyString();
119 return linear_constraints.names(index);
122absl::flat_hash_map<int64_t, double> SparseDoubleVectorAsMap(
123 const SparseDoubleVectorProto& vector) {
124 CHECK_EQ(vector.ids_size(), vector.values_size());
125 absl::flat_hash_map<int64_t, double> result;
126 result.reserve(vector.ids_size());
127 for (
int i = 0;
i < vector.ids_size(); ++
i) {
128 result[vector.ids(i)] = vector.values(i);
137inline int FindRowStart(
const SparseDoubleMatrixProto& matrix,
138 const int64_t row_id,
const int scan_start) {
139 int result = scan_start;
140 while (result < matrix.row_ids_size() && matrix.row_ids(result) < row_id) {
146struct LinearConstraintView {
147 int64_t linear_constraint_id;
150 absl::string_view name;
151 absl::Span<const int64_t> variable_ids;
152 absl::Span<const double> coefficients;
161class LinearConstraintIterator {
163 LinearConstraintIterator(
164 const LinearConstraintsProto* linear_constraints,
165 const SparseDoubleMatrixProto* linear_constraint_matrix)
166 : linear_constraints_(ABSL_DIE_IF_NULL(linear_constraints)),
167 linear_constraint_matrix_(ABSL_DIE_IF_NULL(linear_constraint_matrix)) {
169 const int64_t first_constraint = SafeId(*linear_constraints_, 0);
171 FindRowStart(*linear_constraint_matrix_, first_constraint, 0);
172 matrix_end_ = FindRowStart(*linear_constraint_matrix_,
173 first_constraint + 1, matrix_start_);
176 matrix_end_ = matrix_start_;
180 bool IsDone()
const {
185 LinearConstraintView Current()
const {
187 LinearConstraintView result;
188 result.lower_bound = linear_constraints_->lower_bounds(current_con_);
189 result.upper_bound = linear_constraints_->upper_bounds(current_con_);
190 result.name = SafeName(*linear_constraints_, current_con_);
191 result.linear_constraint_id = SafeId(*linear_constraints_, current_con_);
193 const auto vars_begin = linear_constraint_matrix_->column_ids().data();
194 result.variable_ids = absl::MakeConstSpan(vars_begin + matrix_start_,
195 vars_begin + matrix_end_);
196 const auto coefficients_begins =
197 linear_constraint_matrix_->coefficients().data();
198 result.coefficients = absl::MakeConstSpan(
199 coefficients_begins + matrix_start_, coefficients_begins + matrix_end_);
209 matrix_end_ = matrix_start_;
212 const int64_t current_row_id = SafeId(*linear_constraints_, current_con_);
214 FindRowStart(*linear_constraint_matrix_, current_row_id, matrix_end_);
216 matrix_end_ = FindRowStart(*linear_constraint_matrix_, current_row_id + 1,
222 const LinearConstraintsProto*
const linear_constraints_;
224 const SparseDoubleMatrixProto*
const linear_constraint_matrix_;
227 int current_con_ = 0;
240 int matrix_start_ = 0;
244inline GScipVarType GScipVarTypeFromIsInteger(
const bool is_integer) {
245 return is_integer ? GScipVarType::kInteger : GScipVarType::kContinuous;
250 const gtl::linked_hash_map<int64_t, T>& id_map,
251 const absl::flat_hash_map<T, double>& value_map,
252 const SparseVectorFilterProto& filter) {
254 SparseDoubleVectorProto result;
255 for (
const auto [mathopt_id, scip_id] : id_map) {
256 const double value = value_map.at(scip_id);
257 if (predicate.AcceptsAndUpdate(mathopt_id, value)) {
258 result.add_ids(mathopt_id);
259 result.add_values(value);
267absl::Status GScipSolver::AddVariables(
268 const VariablesProto& variables,
269 const absl::flat_hash_map<int64_t, double>& linear_objective_coefficients) {
271 const int64_t
id = SafeId(variables, i);
276 const bool inverted_bounds =
277 variables.lower_bounds(i) > variables.upper_bounds(i);
281 variables.lower_bounds(i),
282 inverted_bounds ? variables.lower_bounds(i)
283 : variables.upper_bounds(i),
285 GScipVarTypeFromIsInteger(variables.integers(i)),
286 SafeName(variables, i)));
287 if (inverted_bounds) {
288 const double ub = variables.upper_bounds(i);
289 if (gscip_->VarType(v) == GScipVarType::kBinary && ub != 0.0 &&
301 return absl::OkStatus();
304absl::StatusOr<bool> GScipSolver::UpdateVariables(
305 const VariableUpdatesProto& variable_updates) {
306 for (
const auto [
id, is_integer] :
MakeView(variable_updates.integers())) {
309 SCIP_VAR*
const var = variables_.at(
id);
310 if (gscip_->VarType(var) == GScipVarType::kBinary) {
316 gscip_->SetVarType(var, GScipVarTypeFromIsInteger(is_integer)));
318 for (
const auto [
id, lb] :
MakeView(variable_updates.lower_bounds())) {
319 SCIP_VAR*
const var = variables_.at(
id);
320 if (gscip_->VarType(var) == GScipVarType::kBinary) {
327 for (
const auto [
id, ub] :
MakeView(variable_updates.upper_bounds())) {
328 SCIP_VAR*
const var = variables_.at(
id);
329 if (gscip_->VarType(var) == GScipVarType::kBinary) {
344absl::Status GScipSolver::AddQuadraticObjectiveTerms(
345 const SparseDoubleMatrixProto& new_qp_terms,
const bool maximize) {
346 const int num_qp_terms = new_qp_terms.row_ids_size();
347 if (num_qp_terms == 0) {
348 return absl::OkStatus();
351 SCIP_VAR*
const qp_auxiliary_variable,
352 gscip_->AddVariable(-
kInf,
kInf, 1.0, GScipVarType::kContinuous));
353 GScipQuadraticRange range{
354 .lower_bound = maximize ? 0.0 : -
kInf,
355 .linear_variables = {qp_auxiliary_variable},
356 .linear_coefficients = {-1.0},
357 .upper_bound = maximize ?
kInf : 0.0,
359 range.quadratic_variables1.reserve(num_qp_terms);
360 range.quadratic_variables2.reserve(num_qp_terms);
361 range.quadratic_coefficients.reserve(num_qp_terms);
362 for (
int i = 0;
i < num_qp_terms; ++
i) {
363 range.quadratic_variables1.push_back(
364 variables_.at(new_qp_terms.row_ids(i)));
365 range.quadratic_variables2.push_back(
366 variables_.at(new_qp_terms.column_ids(i)));
367 range.quadratic_coefficients.push_back(new_qp_terms.coefficients(i));
370 has_quadratic_objective_ =
true;
371 return absl::OkStatus();
374absl::Status GScipSolver::AddLinearConstraints(
375 const LinearConstraintsProto& linear_constraints,
376 const SparseDoubleMatrixProto& linear_constraint_matrix) {
377 for (LinearConstraintIterator lin_con_it(&linear_constraints,
378 &linear_constraint_matrix);
379 !lin_con_it.IsDone(); lin_con_it.Next()) {
380 const LinearConstraintView current = lin_con_it.Current();
382 GScipLinearRange range;
383 range.lower_bound = current.lower_bound;
384 range.upper_bound = current.upper_bound;
385 range.coefficients = std::vector<double>(current.coefficients.begin(),
386 current.coefficients.end());
387 range.variables.reserve(current.variable_ids.size());
388 for (
const int64_t var_id : current.variable_ids) {
389 range.variables.push_back(variables_.at(var_id));
392 SCIP_CONS*
const scip_con,
393 gscip_->AddLinearConstraint(range, std::string(current.name)));
397 return absl::OkStatus();
400absl::Status GScipSolver::UpdateLinearConstraints(
401 const LinearConstraintUpdatesProto linear_constraint_updates,
402 const SparseDoubleMatrixProto& linear_constraint_matrix,
403 const std::optional<int64_t> first_new_var_id,
404 const std::optional<int64_t> first_new_cstr_id) {
405 for (
const auto [
id, lb] :
406 MakeView(linear_constraint_updates.lower_bounds())) {
408 gscip_->SetLinearConstraintLb(linear_constraints_.at(
id), lb));
410 for (
const auto [
id, ub] :
411 MakeView(linear_constraint_updates.upper_bounds())) {
413 gscip_->SetLinearConstraintUb(linear_constraints_.at(
id), ub));
416 linear_constraint_matrix, 0,
417 first_new_cstr_id, 0,
419 for (
const auto& [var_id, value] : var_coeffs) {
421 linear_constraints_.at(lin_con_id), variables_.at(var_id), value));
424 return absl::OkStatus();
427absl::Status GScipSolver::AddQuadraticConstraints(
428 const google::protobuf::Map<int64_t, QuadraticConstraintProto>&
429 quadratic_constraints) {
430 for (
const auto& [
id, constraint] : quadratic_constraints) {
431 GScipQuadraticRange range{
432 .lower_bound = constraint.lower_bound(),
433 .upper_bound = constraint.upper_bound(),
436 const int num_linear_terms = constraint.linear_terms().ids_size();
437 range.linear_variables.reserve(num_linear_terms);
438 range.linear_coefficients.reserve(num_linear_terms);
439 for (
const auto [var_id, coeff] :
MakeView(constraint.linear_terms())) {
440 range.linear_variables.push_back(variables_.at(var_id));
441 range.linear_coefficients.push_back(coeff);
445 const SparseDoubleMatrixProto& quad_terms = constraint.quadratic_terms();
446 const int num_quad_terms = constraint.quadratic_terms().row_ids_size();
447 range.quadratic_variables1.reserve(num_quad_terms);
448 range.quadratic_variables2.reserve(num_quad_terms);
449 range.quadratic_coefficients.reserve(num_quad_terms);
450 for (
int i = 0;
i < num_quad_terms; ++
i) {
451 range.quadratic_variables1.push_back(
452 variables_.at(quad_terms.row_ids(i)));
453 range.quadratic_variables2.push_back(
454 variables_.at(quad_terms.column_ids(i)));
455 range.quadratic_coefficients.push_back(quad_terms.coefficients(i));
459 gscip_->AddQuadraticConstraint(range, constraint.name()));
462 return absl::OkStatus();
465absl::Status GScipSolver::AddIndicatorConstraints(
466 const google::protobuf::Map<int64_t, IndicatorConstraintProto>&
467 indicator_constraints) {
468 for (
const auto& [
id, constraint] : indicator_constraints) {
469 if (!constraint.has_indicator_id()) {
473 SCIP_VAR*
const indicator_var = variables_.at(constraint.indicator_id());
477 GScipIndicatorConstraint data{
478 .indicator_variable = indicator_var,
479 .negate_indicator = constraint.activate_on_zero(),
481 const double lb = constraint.lower_bound();
482 const double ub = constraint.upper_bound();
485 <<
"gSCIP does not support indicator constraints with ranged "
486 "implied constraints; bounds are: "
487 << lb <<
" <= ... <= " << ub;
493 data.upper_bound = ub;
496 data.upper_bound = -lb;
500 const int num_terms = constraint.expression().ids_size();
501 data.variables.reserve(num_terms);
502 data.coefficients.reserve(num_terms);
503 for (
const auto [var_id, coeff] :
MakeView(constraint.expression())) {
504 data.variables.push_back(variables_.at(var_id));
505 data.coefficients.push_back(scaling * coeff);
509 gscip_->AddIndicatorConstraint(data, constraint.name()));
511 std::make_pair(scip_con, constraint.indicator_id()));
513 return absl::OkStatus();
516absl::StatusOr<std::pair<SCIP_VAR*, SCIP_CONS*>>
517GScipSolver::AddSlackVariableEqualToExpression(
518 const LinearExpressionProto& expression) {
521 gscip_->AddVariable(-
kInf,
kInf, 0.0, GScipVarType::kContinuous));
522 GScipLinearRange range{
523 .lower_bound = -expression.offset(),
524 .upper_bound = -expression.offset(),
526 range.variables.push_back(aux_var);
527 range.coefficients.push_back(-1.0);
528 for (
int i = 0;
i < expression.ids_size(); ++
i) {
529 range.variables.push_back(variables_.at(expression.ids(i)));
530 range.coefficients.push_back(expression.coefficients(i));
532 ASSIGN_OR_RETURN(SCIP_CONS * aux_constr, gscip_->AddLinearConstraint(range));
533 return std::make_pair(aux_var, aux_constr);
536absl::Status GScipSolver::AuxiliaryStructureHandler::DeleteStructure(
538 for (SCIP_CONS*
const constraint : constraints) {
541 for (SCIP_VAR*
const variable : variables) {
546 return absl::OkStatus();
549absl::StatusOr<std::pair<GScipSOSData, GScipSolver::AuxiliaryStructureHandler>>
550GScipSolver::ProcessSosProto(
const SosConstraintProto& sos_constraint) {
552 AuxiliaryStructureHandler handler;
553 absl::flat_hash_set<int64_t> reused_variables;
554 for (
const LinearExpressionProto& expr : sos_constraint.expressions()) {
558 if (expr.ids_size() == 1 && expr.coefficients(0) == 1.0 &&
559 expr.offset() == 0.0 && !reused_variables.contains(expr.ids(0))) {
560 data.variables.push_back(variables_.at(expr.ids(0)));
561 reused_variables.insert(expr.ids(0));
564 AddSlackVariableEqualToExpression(expr));
565 handler.variables.push_back(slack_var);
566 handler.constraints.push_back(slack_constr);
567 data.variables.push_back(slack_var);
570 for (
const double weight : sos_constraint.weights()) {
571 data.weights.push_back(weight);
573 return std::make_pair(data, handler);
576absl::Status GScipSolver::AddSos1Constraints(
577 const google::protobuf::Map<int64_t, SosConstraintProto>&
579 for (
const auto& [
id, constraint] : sos1_constraints) {
581 ProcessSosProto(constraint));
583 SCIP_CONS*
const scip_con,
584 gscip_->AddSOS1Constraint(std::move(sos_data), constraint.name()));
586 std::make_pair(scip_con, std::move(slack_handler)));
588 return absl::OkStatus();
591absl::Status GScipSolver::AddSos2Constraints(
592 const google::protobuf::Map<int64_t, SosConstraintProto>&
594 for (
const auto& [
id, constraint] : sos2_constraints) {
596 ProcessSosProto(constraint));
598 gscip_->AddSOS2Constraint(sos_data, constraint.name()));
600 std::make_pair(scip_con, std::move(slack_handler)));
602 return absl::OkStatus();
608 return GScipParameters::OFF;
610 return GScipParameters::FAST;
611 case EMPHASIS_MEDIUM:
612 case EMPHASIS_UNSPECIFIED:
613 return GScipParameters::DEFAULT_META_PARAM_VALUE;
615 case EMPHASIS_VERY_HIGH:
616 return GScipParameters::AGGRESSIVE;
618 LOG(FATAL) <<
"Unsupported MathOpt Emphasis value: "
620 <<
" unknown, error setting gSCIP parameters";
625 const SolveParametersProto& solve_parameters) {
630 GScipParameters result;
631 std::vector<std::string> warnings;
637 if (solve_parameters.has_time_limit()) {
642 if (solve_parameters.has_iteration_limit()) {
643 warnings.push_back(
"parameter iteration_limit not supported for gSCIP.");
645 if (solve_parameters.has_node_limit()) {
646 (*result.mutable_long_params())[
"limits/totalnodes"] =
647 solve_parameters.node_limit();
649 if (solve_parameters.has_cutoff_limit()) {
650 result.set_objective_limit(solve_parameters.cutoff_limit());
652 if (solve_parameters.has_objective_limit()) {
653 warnings.push_back(
"parameter objective_limit not supported for gSCIP.");
655 if (solve_parameters.has_best_bound_limit()) {
656 warnings.push_back(
"parameter best_bound_limit not supported for gSCIP.");
658 if (solve_parameters.has_solution_limit()) {
659 (*result.mutable_int_params())[
"limits/solutions"] =
660 solve_parameters.solution_limit();
668 result.set_silence_output(!solve_parameters.enable_output());
669 if (solve_parameters.has_threads()) {
672 if (solve_parameters.has_random_seed()) {
675 if (solve_parameters.has_absolute_gap_tolerance()) {
676 (*result.mutable_real_params())[
"limits/absgap"] =
677 solve_parameters.absolute_gap_tolerance();
679 if (solve_parameters.has_relative_gap_tolerance()) {
680 (*result.mutable_real_params())[
"limits/gap"] =
681 solve_parameters.relative_gap_tolerance();
683 if (solve_parameters.has_solution_pool_size()) {
684 result.set_num_solutions(solve_parameters.solution_pool_size());
690 (*result.mutable_int_params())[
"limits/maxsol"] =
691 std::max(100, solve_parameters.solution_pool_size());
692 (*result.mutable_int_params())[
"limits/maxorigsol"] =
693 solve_parameters.solution_pool_size();
696 if (solve_parameters.lp_algorithm() != LP_ALGORITHM_UNSPECIFIED) {
698 switch (solve_parameters.lp_algorithm()) {
699 case LP_ALGORITHM_PRIMAL_SIMPLEX:
702 case LP_ALGORITHM_DUAL_SIMPLEX:
705 case LP_ALGORITHM_BARRIER:
709 "parameter lp_algorithm with value BARRIER is not supported for "
710 "gSCIP in ortools.");
713 case LP_ALGORITHM_FIRST_ORDER:
715 "parameter lp_algorithm with value FIRST_ORDER is not supported.");
718 LOG(FATAL) <<
"LPAlgorithm: "
720 <<
" unknown, error setting gSCIP parameters";
722 (*result.mutable_char_params())[
"lp/initalgorithm"] = alg;
724 if (solve_parameters.presolve() != EMPHASIS_UNSPECIFIED) {
727 if (solve_parameters.cuts() != EMPHASIS_UNSPECIFIED) {
730 if (solve_parameters.heuristics() != EMPHASIS_UNSPECIFIED) {
731 result.set_heuristics(
734 if (solve_parameters.scaling() != EMPHASIS_UNSPECIFIED) {
736 switch (solve_parameters.scaling()) {
741 case EMPHASIS_MEDIUM:
745 case EMPHASIS_VERY_HIGH:
749 LOG(FATAL) <<
"Scaling emphasis: "
751 <<
" unknown, error setting gSCIP parameters";
753 (*result.mutable_int_params())[
"lp/scaling"] = scaling_value;
756 result.MergeFrom(solve_parameters.gscip());
758 if (!warnings.empty()) {
759 return absl::InvalidArgumentError(absl::StrJoin(warnings,
"; "));
766std::string JoinDetails(absl::string_view gscip_detail,
767 absl::string_view math_opt_detail) {
768 if (gscip_detail.empty()) {
769 return std::string(math_opt_detail);
771 if (math_opt_detail.empty()) {
772 return std::string(gscip_detail);
774 return absl::StrCat(gscip_detail,
"; ", math_opt_detail);
777double GetBestPrimalObjective(
778 const double gscip_best_objective,
779 const google::protobuf::RepeatedPtrField<SolutionProto>& solutions,
780 const bool is_maximize) {
781 double result = gscip_best_objective;
782 for (
const auto&
solution : solutions) {
783 if (
solution.has_primal_solution() &&
784 solution.primal_solution().feasibility_status() ==
785 SOLUTION_STATUS_FEASIBLE) {
788 ? std::max(result,
solution.primal_solution().objective_value())
795absl::StatusOr<TerminationProto> ConvertTerminationReason(
796 const GScipResult& gscip_result,
const bool is_maximize,
797 const google::protobuf::RepeatedPtrField<SolutionProto>& solutions,
798 const bool had_cutoff) {
799 const bool has_feasible_solution = !solutions.empty();
800 const GScipOutput::Status gscip_status = gscip_result.gscip_output.status();
801 absl::string_view gscip_status_detail =
802 gscip_result.gscip_output.status_detail();
803 const GScipSolvingStats& gscip_stats = gscip_result.gscip_output.stats();
806 const double best_primal_objective = GetBestPrimalObjective(
807 gscip_stats.best_objective(), solutions, is_maximize);
808 const std::optional<double> optional_finite_primal_objective =
809 has_feasible_solution ? std::make_optional(best_primal_objective)
813 const std::optional<double> optional_dual_objective =
814 std::isfinite(gscip_stats.best_bound())
815 ? std::make_optional(gscip_stats.best_bound())
817 switch (gscip_status) {
818 case GScipOutput::USER_INTERRUPT:
820 is_maximize, LIMIT_INTERRUPTED, optional_finite_primal_objective,
821 optional_dual_objective,
822 JoinDetails(gscip_status_detail,
823 "underlying gSCIP status: USER_INTERRUPT"));
824 case GScipOutput::NODE_LIMIT:
826 is_maximize, LIMIT_NODE, optional_finite_primal_objective,
827 optional_dual_objective,
828 JoinDetails(gscip_status_detail,
829 "underlying gSCIP status: NODE_LIMIT"));
830 case GScipOutput::TOTAL_NODE_LIMIT:
832 is_maximize, LIMIT_NODE, optional_finite_primal_objective,
833 optional_dual_objective,
834 JoinDetails(gscip_status_detail,
835 "underlying gSCIP status: TOTAL_NODE_LIMIT"));
836 case GScipOutput::STALL_NODE_LIMIT:
838 is_maximize, LIMIT_SLOW_PROGRESS, optional_finite_primal_objective,
839 optional_dual_objective, gscip_status_detail);
840 case GScipOutput::TIME_LIMIT:
842 is_maximize, LIMIT_TIME, optional_finite_primal_objective,
843 optional_dual_objective, gscip_status_detail);
844 case GScipOutput::MEM_LIMIT:
846 is_maximize, LIMIT_MEMORY, optional_finite_primal_objective,
847 optional_dual_objective, gscip_status_detail);
848 case GScipOutput::SOL_LIMIT:
850 is_maximize, LIMIT_SOLUTION, optional_finite_primal_objective,
851 optional_dual_objective,
852 JoinDetails(gscip_status_detail,
853 "underlying gSCIP status: SOL_LIMIT"));
854 case GScipOutput::BEST_SOL_LIMIT:
856 is_maximize, LIMIT_SOLUTION, optional_finite_primal_objective,
857 optional_dual_objective,
858 JoinDetails(gscip_status_detail,
859 "underlying gSCIP status: BEST_SOL_LIMIT"));
860 case GScipOutput::RESTART_LIMIT:
862 is_maximize, LIMIT_OTHER, optional_finite_primal_objective,
863 optional_dual_objective,
864 JoinDetails(gscip_status_detail,
865 "underlying gSCIP status: RESTART_LIMIT"));
866 case GScipOutput::OPTIMAL:
868 best_primal_objective,
869 gscip_stats.best_bound(),
870 JoinDetails(gscip_status_detail,
"underlying gSCIP status: OPTIMAL"));
871 case GScipOutput::GAP_LIMIT:
873 best_primal_objective,
874 gscip_stats.best_bound(),
875 JoinDetails(gscip_status_detail,
876 "underlying gSCIP status: GAP_LIMIT"));
877 case GScipOutput::INFEASIBLE:
882 const FeasibilityStatusProto dual_feasibility_status =
883 FEASIBILITY_STATUS_FEASIBLE;
885 gscip_status_detail);
887 case GScipOutput::UNBOUNDED: {
888 if (has_feasible_solution) {
891 JoinDetails(gscip_status_detail,
892 "underlying gSCIP status was UNBOUNDED, both primal "
893 "ray and feasible solution are present"));
897 FEASIBILITY_STATUS_INFEASIBLE,
900 "underlying gSCIP status was UNBOUNDED, but only primal ray "
901 "was given, no feasible solution was found"));
904 case GScipOutput::INF_OR_UNBD:
907 FEASIBILITY_STATUS_UNDETERMINED,
908 JoinDetails(gscip_status_detail,
909 "underlying gSCIP status: INF_OR_UNBD"));
910 case GScipOutput::TERMINATE:
912 is_maximize, LIMIT_INTERRUPTED, optional_finite_primal_objective,
913 optional_dual_objective,
914 JoinDetails(gscip_status_detail,
915 "underlying gSCIP status: TERMINATE"));
916 case GScipOutput::INVALID_SOLVER_PARAMETERS:
917 return absl::InvalidArgumentError(gscip_status_detail);
918 case GScipOutput::UNKNOWN:
919 return absl::InternalError(JoinDetails(
920 gscip_status_detail,
"Unexpected GScipOutput.status: UNKNOWN"));
922 return absl::InternalError(JoinDetails(
923 gscip_status_detail, absl::StrCat(
"Missing GScipOutput.status case: ",
930absl::StatusOr<SolveResultProto> GScipSolver::CreateSolveResultProto(
931 GScipResult gscip_result,
const ModelSolveParametersProto& model_parameters,
932 const std::optional<double> cutoff) {
933 SolveResultProto solve_result;
934 const bool is_maximize = gscip_->ObjectiveIsMaximize();
937 const auto meets_cutoff = [cutoff, is_maximize](
const double obj_value) {
938 if (!cutoff.has_value()) {
942 return obj_value >= *cutoff;
944 return obj_value <= *cutoff;
948 CHECK_EQ(gscip_result.solutions.size(), gscip_result.objective_values.size());
949 for (
int i = 0;
i < gscip_result.solutions.size(); ++
i) {
951 if (!meets_cutoff(gscip_result.objective_values[i])) {
954 SolutionProto*
const solution = solve_result.add_solutions();
956 solution->mutable_primal_solution();
957 primal_solution->set_objective_value(gscip_result.objective_values[i]);
961 model_parameters.variable_values_filter());
963 if (!gscip_result.primal_ray.empty()) {
964 *solve_result.add_primal_rays()->mutable_variable_values() =
966 model_parameters.variable_values_filter());
969 ConvertTerminationReason(gscip_result, is_maximize,
970 solve_result.solutions(),
971 cutoff.has_value()));
972 SolveStatsProto*
const common_stats = solve_result.mutable_solve_stats();
973 const GScipSolvingStats& gscip_stats = gscip_result.gscip_output.stats();
975 common_stats->set_node_count(gscip_stats.node_count());
976 common_stats->set_simplex_iterations(gscip_stats.primal_simplex_iterations() +
977 gscip_stats.dual_simplex_iterations());
978 common_stats->set_barrier_iterations(gscip_stats.barrier_iterations());
979 *solve_result.mutable_gscip_output() = std::move(gscip_result.gscip_output);
983GScipSolver::GScipSolver(std::unique_ptr<GScip> gscip)
984 : gscip_(std::move(ABSL_DIE_IF_NULL(gscip))) {}
987 const ModelProto& model,
const InitArgs&) {
989 ASSIGN_OR_RETURN(std::unique_ptr<GScip> gscip, GScip::Create(model.name()));
991 RETURN_IF_ERROR(gscip->SetObjectiveOffset(model.objective().offset()));
992 auto solver = absl::WrapUnique(
new GScipSolver(std::move(gscip)));
993 RETURN_IF_ERROR(solver->constraint_handler_.Register(solver->gscip_.get()));
996 SparseDoubleVectorAsMap(model.objective().linear_coefficients())));
998 model.objective().quadratic_coefficients(),
999 model.objective().maximize()));
1001 model.linear_constraints(), model.linear_constraint_matrix()));
1003 solver->AddQuadraticConstraints(model.quadratic_constraints()));
1005 solver->AddIndicatorConstraints(model.indicator_constraints()));
1006 RETURN_IF_ERROR(solver->AddSos1Constraints(model.sos1_constraints()));
1007 RETURN_IF_ERROR(solver->AddSos2Constraints(model.sos2_constraints()));
1012 const SolveParametersProto& parameters,
1013 const ModelSolveParametersProto& model_parameters,
1015 const CallbackRegistrationProto& callback_registration,
Callback cb,
1018 model_parameters, kGscipSupportedStructures,
"SCIP"));
1019 const absl::Time start = absl::Now();
1021 GScip::Interrupter gscip_interrupter;
1023 interrupter, [&]() { gscip_interrupter.Interrupt(); });
1024 const bool use_interrupter = interrupter !=
nullptr || cb !=
nullptr;
1027 callback_registration,
1028 {CALLBACK_EVENT_MIP_SOLUTION,
1029 CALLBACK_EVENT_MIP_NODE}));
1030 if (constraint_data_ !=
nullptr) {
1031 return absl::InternalError(
1032 "constraint_data_ should always be null at the start of "
1033 "GScipSolver::Solver()");
1035 SCIP_CONS* callback_cons =
nullptr;
1036 if (cb !=
nullptr) {
1039 constraint_data_ = std::make_unique<GScipSolverConstraintData>();
1040 constraint_data_->user_callback = std::move(cb);
1041 constraint_data_->SetWhenRunAndAdds(callback_registration);
1042 constraint_data_->solve_start_time = absl::Now();
1043 constraint_data_->variables = &variables_;
1044 constraint_data_->variable_node_filter =
1045 &callback_registration.mip_node_filter();
1046 constraint_data_->variable_solution_filter =
1047 &callback_registration.mip_solution_filter();
1048 constraint_data_->interrupter = &gscip_interrupter;
1054 constraint_handler_.AddCallbackConstraint(
1055 gscip_.get(),
"mathopt_callback_constraint",
1056 constraint_data_.get()));
1058 const auto cleanup_constraint_data = absl::Cleanup([
this]() {
1059 if (constraint_data_ !=
nullptr) {
1060 *constraint_data_ = {};
1065 auto message_cb_cleanup = absl::MakeCleanup(
1066 [&buffered_message_callback]() { buffered_message_callback.Flush(); });
1070 for (
const SolutionHintProto& hint : model_parameters.solution_hints()) {
1071 absl::flat_hash_map<SCIP_VAR*, double> partial_solution;
1072 for (
const auto [
id, val] :
MakeView(hint.variable_values())) {
1073 partial_solution.insert({variables_.at(
id), val});
1077 for (
const auto [
id, value] :
1078 MakeView(model_parameters.branching_priorities())) {
1079 RETURN_IF_ERROR(gscip_->SetBranchingPriority(variables_.at(
id), value));
1087 if (buffered_message_callback.has_user_message_callback()) {
1088 gscip_msg_cb = [&buffered_message_callback](
1089 const auto,
const absl::string_view message) {
1090 buffered_message_callback.OnMessage(message);
1095 GScipResult gscip_result,
1096 gscip_->Solve(gscip_parameters,
1097 "", std::move(gscip_msg_cb),
1098 use_interrupter ? &gscip_interrupter :
nullptr));
1101 std::move(message_cb_cleanup).Invoke();
1103 if (callback_cons !=
nullptr) {
1105 constraint_data_.reset();
1109 SolveResultProto result,
1110 CreateSolveResultProto(std::move(gscip_result), model_parameters,
1111 parameters.has_cutoff_limit()
1112 ? std::make_optional(parameters.cutoff_limit())
1117 for (
const auto [
id, unused] :
1118 MakeView(model_parameters.branching_priorities())) {
1123 absl::Now() - start, result.mutable_solve_stats()->mutable_solve_time()));
1127absl::flat_hash_set<SCIP_VAR*> GScipSolver::LookupAllVariables(
1128 absl::Span<const int64_t> variable_ids) {
1129 absl::flat_hash_set<SCIP_VAR*> result;
1130 result.reserve(variable_ids.size());
1131 for (
const int64_t var_id : variable_ids) {
1132 result.insert(variables_.
at(var_id));
1138InvertedBounds GScipSolver::ListInvertedBounds()
const {
1140 InvertedBounds inverted_bounds;
1141 for (
const auto& [
id, var] : variables_) {
1142 if (gscip_->Lb(var) > gscip_->Ub(var)) {
1143 inverted_bounds.variables.push_back(
id);
1146 for (
const auto& [
id, cstr] : linear_constraints_) {
1147 if (gscip_->LinearConstraintLb(cstr) > gscip_->LinearConstraintUb(cstr)) {
1148 inverted_bounds.linear_constraints.push_back(
id);
1153 std::sort(inverted_bounds.variables.begin(), inverted_bounds.variables.end());
1154 std::sort(inverted_bounds.linear_constraints.begin(),
1155 inverted_bounds.linear_constraints.end());
1156 return inverted_bounds;
1160 InvalidIndicators invalid_indicators;
1161 for (
const auto& [constraint_id, gscip_data] : indicator_constraints_) {
1162 if (!gscip_data.has_value()) {
1165 const auto [gscip_constraint, indicator_id] = *gscip_data;
1166 SCIP_VAR*
const indicator_var = variables_.at(indicator_id);
1167 if (gscip_->VarType(indicator_var) == GScipVarType::kContinuous ||
1168 gscip_->Lb(indicator_var) < 0.0 || gscip_->Ub(indicator_var) > 1.0) {
1169 invalid_indicators.invalid_indicators.push_back(
1170 {.variable = indicator_id, .constraint = constraint_id});
1173 invalid_indicators.
Sort();
1174 return invalid_indicators;
1179 ->CanSafeBulkDelete(
1180 LookupAllVariables(model_update.deleted_variable_ids()))
1189 if (has_quadratic_objective_ &&
1190 (model_update.objective_updates().has_direction_update() ||
1191 model_update.objective_updates()
1192 .quadratic_coefficients()
1193 .row_ids_size() > 0)) {
1197 for (
const int64_t constraint_id :
1198 model_update.deleted_linear_constraint_ids()) {
1199 SCIP_CONS*
const scip_cons = linear_constraints_.at(constraint_id);
1200 linear_constraints_.erase(constraint_id);
1204 const absl::flat_hash_set<SCIP_VAR*> vars_to_delete =
1205 LookupAllVariables(model_update.deleted_variable_ids());
1206 for (
const int64_t deleted_variable_id :
1207 model_update.deleted_variable_ids()) {
1208 variables_.erase(deleted_variable_id);
1213 const std::optional<int64_t> first_new_var_id =
1215 const std::optional<int64_t> first_new_cstr_id =
1218 if (model_update.objective_updates().has_direction_update()) {
1220 model_update.objective_updates().direction_update()));
1222 if (model_update.objective_updates().has_offset_update()) {
1224 model_update.objective_updates().offset_update()));
1226 if (
const auto response = UpdateVariables(model_update.variable_updates());
1227 !response.ok() || !(*response)) {
1230 const absl::flat_hash_map<int64_t, double> linear_objective_updates =
1231 SparseDoubleVectorAsMap(
1232 model_update.objective_updates().linear_coefficients());
1233 for (
const auto& obj_pair : linear_objective_updates) {
1235 if (!first_new_var_id.has_value() || obj_pair.first < *first_new_var_id) {
1237 gscip_->SetObjCoef(variables_.at(obj_pair.first), obj_pair.second));
1269 AddVariables(model_update.new_variables(), linear_objective_updates));
1272 model_update.objective_updates().quadratic_coefficients(),
1273 gscip_->ObjectiveIsMaximize()));
1277 UpdateLinearConstraints(model_update.linear_constraint_updates(),
1278 model_update.linear_constraint_matrix_updates(),
1280 first_new_cstr_id));
1283 const std::optional first_new_variable_id =
1285 if (first_new_variable_id.has_value()) {
1286 for (
const auto& [lin_con_id, var_coeffs] :
1290 *first_new_variable_id,
1292 for (
const auto& [var_id, value] : var_coeffs) {
1295 linear_constraints_.at(lin_con_id), variables_.at(var_id), value));
1302 AddLinearConstraints(model_update.new_linear_constraints(),
1303 model_update.linear_constraint_matrix_updates()));
1306 for (
const int64_t constraint_id :
1307 model_update.quadratic_constraint_updates().deleted_constraint_ids()) {
1309 quadratic_constraints_.extract(constraint_id).mapped()));
1312 model_update.quadratic_constraint_updates().new_constraints()));
1315 for (
const int64_t constraint_id :
1316 model_update.indicator_constraint_updates().deleted_constraint_ids()) {
1317 const auto gscip_data =
1318 indicator_constraints_.extract(constraint_id).mapped();
1319 if (!gscip_data.has_value()) {
1322 SCIP_CONS*
const gscip_constraint = gscip_data->first;
1323 CHECK_NE(gscip_constraint,
nullptr);
1327 model_update.indicator_constraint_updates().new_constraints()));
1330 for (
const int64_t constraint_id :
1331 model_update.sos1_constraint_updates().deleted_constraint_ids()) {
1332 auto [gscip_constraint, slack_handler] =
1333 sos1_constraints_.extract(constraint_id).mapped();
1338 model_update.sos1_constraint_updates().new_constraints()));
1341 for (
const int64_t constraint_id :
1342 model_update.sos2_constraint_updates().deleted_constraint_ids()) {
1343 auto [gscip_constraint, slack_handler] =
1344 sos2_constraints_.extract(constraint_id).mapped();
1349 model_update.sos2_constraint_updates().new_constraints()));
1354absl::StatusOr<ComputeInfeasibleSubsystemResultProto>
1358 return absl::UnimplementedError(
1359 "SCIP does not provide a method to compute an infeasible subsystem");
#define ASSIGN_OR_RETURN(lhs, rexpr)
#define RETURN_IF_ERROR(expr)
mapped_type & at(const key_arg< K > &key)
absl::StatusOr< SolveResultProto > Solve(const SolveParametersProto ¶meters, const ModelSolveParametersProto &model_parameters, MessageCallback message_cb, const CallbackRegistrationProto &callback_registration, Callback cb, const SolveInterrupter *interrupter) override
absl::StatusOr< bool > Update(const ModelUpdateProto &model_update) override
static absl::StatusOr< GScipParameters > MergeParameters(const SolveParametersProto &solve_parameters)
absl::StatusOr< ComputeInfeasibleSubsystemResultProto > ComputeInfeasibleSubsystem(const SolveParametersProto ¶meters, MessageCallback message_cb, const SolveInterrupter *interrupter) override
static absl::StatusOr< std::unique_ptr< SolverInterface > > New(const ModelProto &model, const InitArgs &init_args)
std::function< void(const std::vector< std::string > &)> MessageCallback
std::function< absl::StatusOr< CallbackResultProto >( const CallbackDataProto &)> Callback
void InsertOrDie(Collection *const collection, const typename Collection::value_type &value)
const MapUtilMappedT< Collection > & FindWithDefault(const Collection &collection, const KeyType &key, const MapUtilMappedT< Collection > &value)
An object oriented wrapper for quadratic constraints in ModelStorage.
SparseDoubleVectorProto FillSparseDoubleVector(absl::Span< const int64_t > ids_in_order, const absl::flat_hash_map< int64_t, IndexType > &id_map, const glop::StrictITIVector< IndexType, glop::Fractional > &values, const SparseVectorFilterProto &filter)
absl::Status ModelIsSupported(const ModelProto &model, const SupportedProblemStructures &support_menu, const absl::string_view solver_name)
std::optional< int64_t > FirstLinearConstraintId(const LinearConstraintsProto &linear_constraints)
absl::Status ModelSolveParametersAreSupported(const ModelSolveParametersProto &model_parameters, const SupportedProblemStructures &support_menu, const absl::string_view solver_name)
int NumMatrixNonzeros(const SparseDoubleMatrixProto &matrix)
int NumVariables(const VariablesProto &variables)
TerminationProto OptimalTerminationProto(const double finite_primal_objective, const double dual_objective, const absl::string_view detail)
TerminationProto LimitTerminationProto(const bool is_maximize, const LimitProto limit, const std::optional< double > optional_finite_primal_objective, const std::optional< double > optional_dual_objective, const absl::string_view detail)
TerminationProto InfeasibleOrUnboundedTerminationProto(bool is_maximize, const FeasibilityStatusProto dual_feasibility_status, const absl::string_view detail)
bool UpdateIsSupported(const ModelUpdateProto &update, const SupportedProblemStructures &support_menu)
absl::Status CheckRegisteredCallbackEvents(const CallbackRegistrationProto ®istration, const absl::flat_hash_set< CallbackEventProto > &supported_events)
SparseVectorView< T > MakeView(absl::Span< const int64_t > ids, const Collection &values)
TerminationProto CutoffTerminationProto(bool is_maximize, const absl::string_view detail)
Calls NoSolutionFoundTerminationProto() with LIMIT_CUTOFF LIMIT.
int NumConstraints(const LinearConstraintsProto &linear_constraints)
GScipParameters::MetaParamValue ConvertMathOptEmphasis(EmphasisProto emphasis)
SparseSubmatrixRowsView SparseSubmatrixByRows(const SparseDoubleMatrixProto &matrix, const int64_t start_row_id, const std::optional< int64_t > end_row_id, const int64_t start_col_id, const std::optional< int64_t > end_col_id)
std::optional< int64_t > FirstVariableId(const VariablesProto &variables)
TerminationProto InfeasibleTerminationProto(bool is_maximize, const FeasibilityStatusProto dual_feasibility_status, const absl::string_view detail)
TerminationProto UnboundedTerminationProto(const bool is_maximize, const absl::string_view detail)
In SWIG mode, we don't want anything besides these top-level includes.
Select next search node to expand Select next item_i to add this new search node to the search Generate a new search node where item_i is not in the knapsack Check validity of this new partial solution(using propagators) - If valid
void GScipSetRandomSeed(GScipParameters *parameters, int random_seed)
void GScipSetCatchCtrlC(const bool catch_ctrl_c, GScipParameters *const parameters)
Sets the misc/catchctrlc property.
std::string ProtoEnumToString(ProtoEnumType enum_value)
void GScipSetMaxNumThreads(int num_threads, GScipParameters *parameters)
CHECK fails if num_threads < 1.
void GScipSetTimeLimit(absl::Duration time_limit, GScipParameters *parameters)
std::function< void(GScipMessageType type, absl::string_view message)> GScipMessageHandler
inline ::absl::StatusOr< absl::Duration > DecodeGoogleApiProto(const google::protobuf::Duration &proto)
inline ::absl::StatusOr< google::protobuf::Duration > EncodeGoogleApiProto(absl::Duration d)
StatusBuilder InvalidArgumentErrorBuilder()
#define MATH_OPT_REGISTER_SOLVER(solver_type, solver_factory)
void Sort()
Sort the elements lexicographically by (constraint ID, variable ID).
Initialization arguments.