diff --git a/ortools/com/google/ortools/util/NestedArrayHelper.cs b/ortools/com/google/ortools/util/NestedArrayHelper.cs index 275adf4653..c401b8236e 100644 --- a/ortools/com/google/ortools/util/NestedArrayHelper.cs +++ b/ortools/com/google/ortools/util/NestedArrayHelper.cs @@ -46,4 +46,4 @@ public static class NestedArrayHelper return result; } } -} // namespace Google.OrTools \ No newline at end of file +} // namespace Google.OrTools diff --git a/ortools/graph/ebert_graph.h b/ortools/graph/ebert_graph.h index 4266fcc5ca..9dea694fd4 100644 --- a/ortools/graph/ebert_graph.h +++ b/ortools/graph/ebert_graph.h @@ -846,7 +846,7 @@ class ForwardStaticGraph // To be used in a DCHECK(). bool TailArrayComplete() const { - CHECK_NOTNULL(tail_.get()); + CHECK(tail_); for (ArcIndexType arc = kFirstArc; arc < num_arcs_; ++arc) { CHECK(CheckTailIndexValidity(arc)); CHECK(IsNodeValid((*tail_)[arc])); @@ -1684,7 +1684,7 @@ class ForwardEbertGraph // To be used in a DCHECK(). bool TailArrayComplete() const { - CHECK_NOTNULL(tail_.get()); + CHECK(tail_); for (ArcIndexType arc = kFirstArc; arc < num_arcs_; ++arc) { CHECK(CheckTailIndexValidity(arc)); CHECK(IsNodeValid((*tail_)[arc])); @@ -1778,7 +1778,7 @@ class ForwardEbertGraph // used. void SetTail(const ArcIndexType arc, const NodeIndexType tail) { DCHECK(CheckTailIndexValidity(arc)); - CHECK_NOTNULL(tail_.get()); + CHECK(tail_); representation_clean_ = false; tail_->Set(arc, tail); } diff --git a/ortools/graph/python/graph.i b/ortools/graph/python/graph.i index 15b5dc630b..c21c0ab70f 100644 --- a/ortools/graph/python/graph.i +++ b/ortools/graph/python/graph.i @@ -140,6 +140,7 @@ %include "ortools/graph/assignment.h" + %unignore operations_research::DijkstraShortestPath; %unignore operations_research::BellmanFordShortestPath; %unignore operations_research::AStarShortestPath; diff --git a/ortools/linear_solver/linear_solver.cc b/ortools/linear_solver/linear_solver.cc index c7dea9550c..b668eebb26 100644 --- a/ortools/linear_solver/linear_solver.cc +++ b/ortools/linear_solver/linear_solver.cc @@ -622,10 +622,11 @@ void MPSolver::SolveWithProto(const MPModelRequest& model_request, std::string error_message; response->set_status(solver.LoadModelFromProto(model, &error_message)); if (response->status() != MPSOLVER_MODEL_IS_VALID) { - // LOG_EVERY_N_SEC(WARNING, 1.0) - // << "Loading model from protocol buffer failed, load status = " - // << MPSolverResponseStatus_Name(response->status()) << " (" - // << response->status() << "); Error: " << error_message; + LOG(WARNING) + << "Loading model from protocol buffer failed, load status = " + << MPSolverResponseStatus_Name(response->status()) << " (" + << response->status() << "); Error: " << error_message; + return; } if (model_request.has_solver_time_limit_seconds()) { diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index 6897644d42..bfc39db38b 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -1854,8 +1854,9 @@ void PresolveCpModel(const CpModelProto& initial_model, if (context.domains[r.representative].IntersectWith(implied)) { LOG(WARNING) << "Domain of " << r.representative << " was not fully propagated using the affine relation " - << "(representative =" << r.representative << ", coeff = " - << r.coeff << ", offset = " << r.offset << ")"; + << "(representative =" << r.representative + << ", coeff = " << r.coeff << ", offset = " << r.offset + << ")"; } proto = context.mapping_model; } else { diff --git a/ortools/sat/integer.cc b/ortools/sat/integer.cc index 7bfe20fd22..acb048447f 100644 --- a/ortools/sat/integer.cc +++ b/ortools/sat/integer.cc @@ -1284,6 +1284,17 @@ std::function UnassignedVarWithLowestMinAtItsMinHeuristic( }; } +std::function SequentialSearch( + std::vector> heuristics) { + return [heuristics]() { + for (const auto& h : heuristics) { + const LiteralIndex li = h(); + if (li != kNoLiteralIndex) return li; + } + return kNoLiteralIndex; + }; +} + SatSolver::Status SolveIntegerProblemWithLazyEncoding( const std::vector& assumptions, const std::function& next_decision, Model* model) { diff --git a/ortools/sat/integer.h b/ortools/sat/integer.h index 6c2e1839af..05591fa360 100644 --- a/ortools/sat/integer.h +++ b/ortools/sat/integer.h @@ -1135,6 +1135,12 @@ std::function FirstUnassignedVarAtItsMinHeuristic( std::function UnassignedVarWithLowestMinAtItsMinHeuristic( const std::vector& vars, Model* model); +// Combines search heuristics in order: if the i-th one returns kNoLiteralIndex, +// ask the (i+1)-th. If every heuristic returned kNoLiteralIndex, +// returns kNoLiteralIndex. +std::function SequentialSearch( + std::vector> heuristics); + // Same as ExcludeCurrentSolutionAndBacktrack() but this version works for an // integer problem with optional variables. The issue is that an optional // variable that is ignored can basically take any value, and we don't really diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc index 7002dcb4d0..a82ae6ac5a 100644 --- a/ortools/sat/linear_programming_constraint.cc +++ b/ortools/sat/linear_programming_constraint.cc @@ -38,7 +38,8 @@ namespace sat { const double LinearProgrammingConstraint::kEpsilon = 1e-6; LinearProgrammingConstraint::LinearProgrammingConstraint(Model* model) - : integer_trail_(model->GetOrCreate()) { + : integer_trail_(model->GetOrCreate()), + dispatcher_(model->GetOrCreate()) { // TODO(user): Find a way to make GetOrCreate() construct it by // default. time_limit_ = model->Mutable(); @@ -79,6 +80,8 @@ glop::ColIndex LinearProgrammingConstraint::GetOrCreateMirrorVariable( integer_variables_.push_back(positive_variable); mirror_lp_variables_.push_back(lp_data_.CreateNewVariable()); lp_solution_.push_back(std::numeric_limits::infinity()); + lp_reduced_cost_.push_back(0.0); + (*dispatcher_)[positive_variable] = this; } return mirror_lp_variables_[integer_variable_to_index_[positive_variable]]; } @@ -183,6 +186,16 @@ glop::Fractional LinearProgrammingConstraint::GetVariableValueAtCpScale( return simplex_.GetVariableValue(var) / scaler_.col_scale(var); } +double LinearProgrammingConstraint::GetSolutionValue( + IntegerVariable variable) const { + return lp_solution_[FindOrDie(integer_variable_to_index_, variable)]; +} + +double LinearProgrammingConstraint::GetSolutionReducedCost( + IntegerVariable variable) const { + return lp_reduced_cost_[FindOrDie(integer_variable_to_index_, variable)]; +} + bool LinearProgrammingConstraint::Propagate() { // Copy CP state into LP state. const int num_vars = integer_variables_.size(); @@ -243,8 +256,12 @@ bool LinearProgrammingConstraint::Propagate() { } lp_data_.ScaleObjective(); } + const double objective_scale = lp_data_.objective_scaling_factor(); for (int i = 0; i < num_vars; i++) { lp_solution_[i] = GetVariableValueAtCpScale(mirror_lp_variables_[i]); + lp_reduced_cost_[i] = simplex_.GetReducedCost(mirror_lp_variables_[i]) * + scaler_.col_scale(mirror_lp_variables_[i]) * + objective_scale; } // We currently ignore the objective and return right away when we don't @@ -320,8 +337,12 @@ bool LinearProgrammingConstraint::Propagate() { // Copy current LP solution. if (simplex_.GetProblemStatus() == glop::ProblemStatus::OPTIMAL) { + const double objective_scale = lp_data_.objective_scaling_factor(); for (int i = 0; i < num_vars; i++) { lp_solution_[i] = GetVariableValueAtCpScale(mirror_lp_variables_[i]); + lp_reduced_cost_[i] = simplex_.GetReducedCost(mirror_lp_variables_[i]) * + scaler_.col_scale(mirror_lp_variables_[i]) * + objective_scale; } } return true; @@ -408,5 +429,147 @@ void LinearProgrammingConstraint::ReducedCostStrengtheningDeductions( } } +std::function HeuristicLPMostInfeasibleBinary(Model* model) { + // Gather all 0-1 variables that appear in some LP. + IntegerTrail* integer_trail = model->GetOrCreate(); + LinearProgrammingDispatcher* dispatcher = + model->GetOrCreate(); + std::vector variables; + for (const auto entry : *dispatcher) { + IntegerVariable var = entry.first; + if (integer_trail->LowerBound(var) == 0 && + integer_trail->UpperBound(var) == 1) { + variables.push_back(var); + } + } + std::sort(variables.begin(), variables.end()); + + LOG(INFO) << "HeuristicLPMostInfeasibleBinary has " << variables.size() + << " variables."; + + IntegerEncoder* integer_encoder = model->GetOrCreate(); + SatSolver* sat_solver = model->GetOrCreate(); + return [variables, dispatcher, integer_trail, integer_encoder, sat_solver]() { + const double kEpsilon = 1e-6; + // Find most fractional value. + IntegerVariable fractional_var = kNoIntegerVariable; + double fractional_distance_best = -1.0; + for (const IntegerVariable var : variables) { + // Check variable is not ignored and unfixed. + if (integer_trail->IsCurrentlyIgnored(var)) continue; + const IntegerValue lb = integer_trail->LowerBound(var); + const IntegerValue ub = integer_trail->UpperBound(var); + if (lb == ub) continue; + + // Check variable's support is fractional. + LinearProgrammingConstraint* lp = (*dispatcher)[var]; + const double lp_value = lp->GetSolutionValue(var); + const double fractional_distance = + std::min(std::ceil(lp_value - kEpsilon) - lp_value, + lp_value - std::floor(lp_value + kEpsilon)); + if (fractional_distance < kEpsilon) continue; + + // Keep variable if it is farther from integrality than the previous. + if (fractional_distance > fractional_distance_best) { + fractional_var = var; + fractional_distance_best = fractional_distance; + } + } + + if (fractional_var != kNoIntegerVariable) { + return integer_encoder + ->GetOrCreateAssociatedLiteral( + IntegerLiteral::GreaterOrEqual(fractional_var, IntegerValue(1))) + .Index(); + } + return kNoLiteralIndex; + }; +} + +std::function HeuristicLPPseudoCostBinary(Model* model) { + // Gather all 0-1 variables that appear in some LP. + IntegerTrail* integer_trail = model->GetOrCreate(); + LinearProgrammingDispatcher* dispatcher = + model->GetOrCreate(); + std::vector variables; + for (const auto entry : *dispatcher) { + IntegerVariable var = entry.first; + if (integer_trail->LowerBound(var) == 0 && + integer_trail->UpperBound(var) == 1) { + variables.push_back(var); + } + } + std::sort(variables.begin(), variables.end()); + + LOG(INFO) << "HeuristicLPPseudoCostBinary has " << variables.size() + << " variables."; + + // Store average of reduced cost from 1 to 0. The best heuristic only sets + // variables to one and cares about cost to zero, even though classic + // pseudocost will use max_var std::min(cost_to_one[var], cost_to_zero[var]). + const int num_vars = variables.size(); + std::vector cost_to_zero(num_vars, 0.0); + std::vector num_cost_to_zero(num_vars); + int num_calls = 0; + + IntegerEncoder* integer_encoder = model->GetOrCreate(); + return [=]() mutable { + const double kEpsilon = 1e-6; + + // Every 10000 calls, decay pseudocosts. + num_calls++; + if (num_calls == 10000) { + for (int i = 0; i < num_vars; i++) { + cost_to_zero[i] /= 2; + num_cost_to_zero[i] /= 2; + } + num_calls = 0; + } + + // Accumulate pseudo-costs of all unassigned variables. + for (int i = 0; i < num_vars; i++) { + const IntegerVariable var = variables[i]; + if (integer_trail->LowerBound(var) == integer_trail->UpperBound(var)) + continue; + + LinearProgrammingConstraint* lp = (*dispatcher)[var]; + const double rc = lp->GetSolutionReducedCost(var); + // Skip reduced costs that are nonzero because of numerical issues. + if (std::abs(rc) < kEpsilon) continue; + + const double value = std::round(lp->GetSolutionValue(var)); + if (value == 1.0 && rc < 0.0) { + cost_to_zero[i] -= rc; + num_cost_to_zero[i]++; + } + } + + // Select noninstantiated variable with highest pseudo-cost. + int selected_index = -1; + double best_cost = 0.0; + for (int i = 0; i < num_vars; i++) { + const IntegerVariable var = variables[i]; + if (integer_trail->LowerBound(var) == integer_trail->UpperBound(var)) { + continue; + } + + if (num_cost_to_zero[i] > 0 && + best_cost < cost_to_zero[i] / num_cost_to_zero[i]) { + best_cost = cost_to_zero[i] / num_cost_to_zero[i]; + selected_index = i; + } + } + + if (selected_index >= 0) { + const Literal decision = integer_encoder->GetOrCreateAssociatedLiteral( + IntegerLiteral::GreaterOrEqual(variables[selected_index], + IntegerValue(1))); + return decision.Index(); + } + + return kNoLiteralIndex; + }; +} + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/linear_programming_constraint.h b/ortools/sat/linear_programming_constraint.h index 87d22e3f89..20d3f4549a 100644 --- a/ortools/sat/linear_programming_constraint.h +++ b/ortools/sat/linear_programming_constraint.h @@ -62,6 +62,7 @@ namespace sat { // // TODO(user): Work with scaled version of the model, maybe by using // LPSolver instead of RevisedSimplex. +class LinearProgrammingDispatcher; class LinearProgrammingConstraint : public PropagatorInterface { public: typedef glop::RowIndex ConstraintIndex; @@ -86,6 +87,11 @@ class LinearProgrammingConstraint : public PropagatorInterface { // the arguments passed to SetObjectiveCoefficient(). void SetMainObjectiveVariable(IntegerVariable ivar) { objective_cp_ = ivar; } + // Returns the LP value and reduced cost of a variable in the current + // solution. + double GetSolutionValue(IntegerVariable variable) const; + double GetSolutionReducedCost(IntegerVariable variable) const; + // PropagatorInterface API. bool Propagate() override; @@ -156,14 +162,55 @@ class LinearProgrammingConstraint : public PropagatorInterface { // On IncrementalPropagate(), if the bound updates do not invalidate this // solution, Propagate() will not find domain reductions, no need to call it. std::vector lp_solution_; + std::vector lp_reduced_cost_; // Linear constraints cannot be created or modified after this is registered. bool lp_constraint_is_registered_ = false; // Time limit (shared with, owned by the sat solver). TimeLimit* time_limit_; + + // The dispatcher for all LP propagators of the model, allows to find which + // LinearProgrammingConstraint has a given IntegerVariable. + LinearProgrammingDispatcher* dispatcher_; }; +// A class that stores which LP propagator is associated to each variable. +// We need to give the hash_map a name so it can be used as a singleton +// in our model. +class LinearProgrammingDispatcher + : public std::unordered_map { + public: + explicit LinearProgrammingDispatcher(Model* model) {} +}; + +// Returns a LiteralIndex guided by the underlying LP constraints. +// This looks at all unassigned 0-1 variables, takes the one with +// a support value closest to 0.5, and tries to assign it to 1. +// If all 0-1 variables have an integer support, returns kNoLiteralIndex. +// Tie-breaking is done using the variable natural order. +// +// TODO(user): This fixes to 1, but for some problems fixing to 0 +// or to the std::round(support value) might work better. When this is the +// case, change behaviour automatically? +std::function HeuristicLPMostInfeasibleBinary(Model* model); + +// Returns a LiteralIndex guided by the underlying LP constraints. +// This computes the mean of reduced costs over successive calls, +// and tries to fix the variable which has the highest reduced cost. +// Tie-breaking is done using the variable natural order. +// +// TODO(user): Try to get better pseudocosts than averaging every time the +// heuristic is called. MIP solvers initialize this with strong branching, then +// keep track of the pseudocosts when doing tree search. Also, this version only +// branches on var >= 1 and keeps track of reduced costs from var = 1 to var = +// 0. This works better than the conventional MIP where the chosen variable will +// be argmax_var std::min(pseudocost_var(0->1), pseudocost_var(1->0)), probably +// because we are doing DFS search where MIP does BFS. This might depend on the +// model, more trials are necessary. We could also do exponential smoothing +// instead of decaying every N calls, i.e. pseudo = a * pseudo + (1-a) reduced. +std::function HeuristicLPPseudoCostBinary(Model* model); + } // namespace sat } // namespace operations_research diff --git a/ortools/util/rev.h b/ortools/util/rev.h index e31dd3f17f..6e5efd713e 100644 --- a/ortools/util/rev.h +++ b/ortools/util/rev.h @@ -121,7 +121,9 @@ class RevMap : ReversibleInterface { void SetLevel(int level) final; int Level() const { return first_op_index_of_next_level_.size(); } - bool ContainsKey(key_type key) const { return operations_research::ContainsKey(map_, key); } + bool ContainsKey(key_type key) const { + return operations_research::ContainsKey(map_, key); + } const mapped_type& FindOrDie(key_type key) const { return operations_research::FindOrDie(map_, key); } diff --git a/ortools/util/stats.h b/ortools/util/stats.h index 3868949b23..f4ff77e5a2 100644 --- a/ortools/util/stats.h +++ b/ortools/util/stats.h @@ -70,7 +70,10 @@ #include #include - +#ifdef HAS_PERF_SUBSYSTEM +#include "absl/ortools/base/str_replace.h" +#include "exegesis/exegesis/itineraries/perf_subsystem.h" +#endif // HAS_PERF_SUBSYSTEM #include "ortools/base/timer.h" namespace operations_research { @@ -94,7 +97,8 @@ class Stat { // Only used for display purposes. std::string Name() const { return name_; } - // Returns a human-readable formated line of the form "name: ValueAsString()". + // Returns a human-readable formatted line of the form "name: + // ValueAsString()". std::string StatString() const; // At display, stats are displayed by decreasing priority, then decreasing @@ -323,9 +327,56 @@ class DisabledScopedTimeDistributionUpdater { DISALLOW_COPY_AND_ASSIGN(DisabledScopedTimeDistributionUpdater); }; +#ifdef HAS_PERF_SUBSYSTEM +// Helper classes to count instructions during execution of a block of code and +// add print the results to logs. +// Creates new perf subsystem and start collecting 'inst_retired:any_p' event on +// creation and stops collecting on destruction. +class EnabledScopedInstructionCounter { + public: + explicit EnabledScopedInstructionCounter(const std::string& name) : name_(name) { + perf_subsystem_.CleanUp(); + perf_subsystem_.AddEvent("inst_retired:any_p:u,p"); + perf_subsystem_.AddEvent("cycles:u,p"); + perf_subsystem_.StartCollecting(); + } + EnabledScopedInstructionCounter(const EnabledScopedInstructionCounter&) = + delete; + EnabledScopedInstructionCounter& operator=( + const EnabledScopedInstructionCounter&) = delete; + ~EnabledScopedInstructionCounter() { + exegesis::PerfResult perf_result = perf_subsystem_.StopAndReadCounters(); + LOG(INFO) << name_ << ": " << perf_result.ToString(); + } + + // Used only for testing. + exegesis::PerfResult GetPerfResult() { + return perf_subsystem_.ReadCounters(); + } + + private: + exegesis::PerfSubsystem perf_subsystem_; + std::string name_; +}; +#endif // HAS_PERF_SUBSYSTEM + +class DisabledScopedInstructionCounter { + public: + explicit DisabledScopedInstructionCounter(const std::string& name) {} + DisabledScopedInstructionCounter(const DisabledScopedInstructionCounter&) = + delete; + DisabledScopedInstructionCounter& operator=( + const DisabledScopedInstructionCounter&) = delete; +}; + #ifdef OR_STATS using ScopedTimeDistributionUpdater = EnabledScopedTimeDistributionUpdater; +#ifdef HAS_PERF_SUBSYSTEM +using ScopedInstructionCounter = EnabledScopedInstructionCounter; +#else // HAS_PERF_SUBSYSTEM +using ScopedInstructionCounter = DisabledScopedInstructionCounter; +#endif // HAS_PERF_SUBSYSTEM // Simple macro to be used by a client that want to execute costly operations // only if OR_STATS is defined. @@ -342,14 +393,28 @@ using ScopedTimeDistributionUpdater = EnabledScopedTimeDistributionUpdater; operations_research::ScopedTimeDistributionUpdater scoped_time_stat( \ (stats)->LookupOrCreateTimeDistribution(__FUNCTION__)) +#ifdef HAS_PERF_SUBSYSTEM + +inline std::string RemoveOperationsResearchAndGlop(const std::string& pretty_function) { + return strings::GlobalReplaceSubstrings( + pretty_function, {{"operations_research::", ""}, {"glop::", ""}}); +} + +#define SCOPED_INSTRUCTION_COUNT \ + operations_research::ScopedInstructionCounter scoped_instruction_count( \ + RemoveOperationsResearchAndGlop(__PRETTY_FUNCTION__)) +#endif // HAS_PERF_SUBSYSTEM + #else // OR_STATS // If OR_STATS is not defined, we remove some instructions that may be time // consuming. using ScopedTimeDistributionUpdater = DisabledScopedTimeDistributionUpdater; +using ScopedInstructionCounter = DisabledScopedInstructionCounter; #define IF_STATS_ENABLED(instructions) #define SCOPED_TIME_STAT(stats) +#define SCOPED_INSTRUCTION_COUNT #endif // OR_STATS