cleanups and minor syncs; sat: add LP guided heuristics

ortools/com/google/ortools/util/NestedArrayHelper.cs

@@ -46,4 +46,4 @@ public static class NestedArrayHelper
     return result;
   }
 }
-}  // namespace Google.OrTools
+}  // namespace Google.OrTools

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);
   }

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;

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()) {

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 {

ortools/sat/integer.cc

@@ -1284,6 +1284,17 @@ std::function<LiteralIndex()> UnassignedVarWithLowestMinAtItsMinHeuristic(
   };
 }
 
+std::function<LiteralIndex()> SequentialSearch(
+    std::vector<std::function<LiteralIndex()>> 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<Literal>& assumptions,
     const std::function<LiteralIndex()>& next_decision, Model* model) {

ortools/sat/integer.h

@@ -1135,6 +1135,12 @@ std::function<LiteralIndex()> FirstUnassignedVarAtItsMinHeuristic(
 std::function<LiteralIndex()> UnassignedVarWithLowestMinAtItsMinHeuristic(
     const std::vector<IntegerVariable>& 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<LiteralIndex()> SequentialSearch(
+    std::vector<std::function<LiteralIndex()>> 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

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<IntegerTrail>()) {
+    : integer_trail_(model->GetOrCreate<IntegerTrail>()),
+      dispatcher_(model->GetOrCreate<LinearProgrammingDispatcher>()) {
   // TODO(user): Find a way to make GetOrCreate<TimeLimit>() construct it by
   // default.
   time_limit_ = model->Mutable<TimeLimit>();
@@ -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<double>::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<LiteralIndex()> HeuristicLPMostInfeasibleBinary(Model* model) {
+  // Gather all 0-1 variables that appear in some LP.
+  IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
+  LinearProgrammingDispatcher* dispatcher =
+      model->GetOrCreate<LinearProgrammingDispatcher>();
+  std::vector<IntegerVariable> 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<IntegerEncoder>();
+  SatSolver* sat_solver = model->GetOrCreate<SatSolver>();
+  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<LiteralIndex()> HeuristicLPPseudoCostBinary(Model* model) {
+  // Gather all 0-1 variables that appear in some LP.
+  IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
+  LinearProgrammingDispatcher* dispatcher =
+      model->GetOrCreate<LinearProgrammingDispatcher>();
+  std::vector<IntegerVariable> 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<double> cost_to_zero(num_vars, 0.0);
+  std::vector<int> num_cost_to_zero(num_vars);
+  int num_calls = 0;
+
+  IntegerEncoder* integer_encoder = model->GetOrCreate<IntegerEncoder>();
+  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

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<double> lp_solution_;
+  std::vector<double> 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<IntegerVariable, LinearProgrammingConstraint*> {
+ 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<LiteralIndex()> 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<LiteralIndex()> HeuristicLPPseudoCostBinary(Model* model);
+
 }  // namespace sat
 }  // namespace operations_research
 

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);
   }

ortools/util/stats.h

@@ -70,7 +70,10 @@
 
 #include <map>
 #include <string>
-
+#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