tweak CP-SAT // Search; improve floating point scaling in CP-SAT

ortools/sat/cp_model_presolve.cc

@@ -33,6 +33,7 @@
 #include "ortools/base/map_util.h"
 #include "ortools/base/stl_util.h"
 #include "ortools/port/proto_utils.h"
+#include "ortools/sat/cp_model.pb.h"
 #include "ortools/sat/cp_model_checker.h"
 #include "ortools/sat/cp_model_loader.h"
 #include "ortools/sat/cp_model_objective.h"
@@ -353,6 +354,7 @@ struct PresolveContext {
   std::vector<int> tmp_literals;
   std::vector<Domain> tmp_term_domains;
   std::vector<Domain> tmp_left_domains;
+  absl::flat_hash_set<int> tmp_literal_set;
 
   // Each time a domain is modified this is set to true.
   SparseBitset<int64> modified_domains;
@@ -418,12 +420,9 @@ bool PresolveBoolOr(ConstraintProto* ct, PresolveContext* context) {
   }
 
   // Inspects the literals and deal with fixed ones.
-  //
-  // TODO(user): detect if one literal is the negation of another in which
-  // case the constraint is true. Remove duplicates too. Do the same for
-  // the PresolveBoolAnd() function.
   bool changed = false;
   context->tmp_literals.clear();
+  context->tmp_literal_set.clear();
   for (const int literal : ct->bool_or().literals()) {
     if (context->LiteralIsFalse(literal)) {
       changed = true;
@@ -441,8 +440,17 @@ bool PresolveBoolOr(ConstraintProto* ct, PresolveContext* context) {
       context->SetLiteralToTrue(literal);
       return RemoveConstraint(ct, context);
     }
-    context->tmp_literals.push_back(literal);
+    if (context->tmp_literal_set.contains(NegatedRef(literal))) {
+      context->UpdateRuleStats("bool_or: always true");
+      return RemoveConstraint(ct, context);
+    }
+
+    if (!context->tmp_literal_set.contains(literal)) {
+      context->tmp_literal_set.insert(literal);
+      context->tmp_literals.push_back(literal);
+    }
   }
+  context->tmp_literal_set.clear();
 
   if (context->tmp_literals.empty()) {
     context->UpdateRuleStats("bool_or: empty");
@@ -2181,7 +2189,7 @@ void PresolvePureSatPart(PresolveContext* context) {
   for (int i = 0; i < num_passes; ++i) {
     const int old_num_clause = postsolver.NumClauses();
     if (!presolver.Presolve(can_be_removed)) {
-      LOG(INFO) << "UNSAT during SAT presolve.";
+      VLOG(1) << "UNSAT during SAT presolve.";
       context->is_unsat = true;
       return;
     }
@@ -2579,6 +2587,206 @@ bool PresolveOneConstraint(int c, PresolveContext* context) {
   }
 }
 
+// Returns the sorted list of literals for given bool_or or at_most_one
+// constraint.
+std::vector<int> GetLiteralsFromSetPPCConstraint(ConstraintProto* ct) {
+  std::vector<int> sorted_literals;
+  if (ct->constraint_case() == ConstraintProto::ConstraintCase::kAtMostOne) {
+    for (const int literal : ct->at_most_one().literals()) {
+      sorted_literals.push_back(literal);
+    }
+  } else if (ct->constraint_case() ==
+             ConstraintProto::ConstraintCase::kBoolOr) {
+    for (const int literal : ct->bool_or().literals()) {
+      sorted_literals.push_back(literal);
+    }
+  }
+  std::sort(sorted_literals.begin(), sorted_literals.end());
+  return sorted_literals;
+}
+
+// Removes dominated constraints or fixes some variables for given pair of
+// setppc constraints. This assumes that literals in constraint c1 is subset of
+// literals in constraint c2.
+bool ProcessSetPPCSubset(int c1, int c2, const std::vector<int>& c2_minus_c1,
+                         const std::vector<int>& original_constraint_index,
+                         std::vector<bool>* marked_for_removal,
+                         PresolveContext* context) {
+  CHECK(!(*marked_for_removal)[c1]);
+  CHECK(!(*marked_for_removal)[c2]);
+  ConstraintProto* ct1 = context->working_model->mutable_constraints(
+      original_constraint_index[c1]);
+  ConstraintProto* ct2 = context->working_model->mutable_constraints(
+      original_constraint_index[c2]);
+  if (ct1->constraint_case() == ConstraintProto::ConstraintCase::kBoolOr &&
+      ct2->constraint_case() == ConstraintProto::ConstraintCase::kAtMostOne) {
+    // fix extras in c2 to 0
+    for (const int literal : c2_minus_c1) {
+      context->SetLiteralToFalse(literal);
+      context->UpdateRuleStats("setppc: fixed variables");
+    }
+    return true;
+  }
+  if (ct1->constraint_case() == ct2->constraint_case()) {
+    if (ct1->constraint_case() == ConstraintProto::ConstraintCase::kBoolOr) {
+      (*marked_for_removal)[c2] = true;
+      context->UpdateRuleStats("setppc: removed dominated constraints");
+      return false;
+    }
+    CHECK_EQ(ct1->constraint_case(),
+             ConstraintProto::ConstraintCase::kAtMostOne);
+    (*marked_for_removal)[c1] = true;
+    context->UpdateRuleStats("setppc: removed dominated constraints");
+    return false;
+  }
+  return false;
+}
+
+// SetPPC is short for set packing, partitioning and covering constraints. These
+// are sum of booleans <=, = and >= 1 respectively.
+// TODO(user,user): Process kBoolAnd too.
+bool ProcessSetPPC(PresolveContext* context, TimeLimit* time_limit) {
+  bool changed = false;
+  const int num_constraints = context->working_model->constraints_size();
+
+  // Signatures of all the constraints. In the signature the bit i is 1 if it
+  // contains a literal l such that l%64 = i.
+  std::vector<uint64> signatures;
+
+  // Graph of constraints to literals. constraint_literals[c] contains all the
+  // literals in constraint indexed by 'c' in sorted order.
+  std::vector<std::vector<int>> constraint_literals;
+
+  // Graph of literals to constraints. literals_to_constraints[l] contains the
+  // vector of constraint indices in which literal 'l' or 'neg(l)' appears.
+  std::vector<std::vector<int>> literals_to_constraints;
+
+  // vector of booleans indicating if the constraint is marked for removal. Note
+  // that we don't remove constraints while processing them but remove all the
+  // marked ones at the end.
+  std::vector<bool> marked_for_removal;
+
+  // The containers above use the local indices for setppc constraints. We store
+  // the original constraint indices corresponding to those local indices here.
+  std::vector<int> original_constraint_index;
+
+  // Fill the initial constraint <-> literal graph, compute signatures and
+  // initialize other containers defined above.
+  int num_setppc_constraints = 0;
+  for (int c = 0; c < num_constraints; ++c) {
+    ConstraintProto* ct = context->working_model->mutable_constraints(c);
+    if (ct->constraint_case() == ConstraintProto::ConstraintCase::kBoolOr ||
+        ct->constraint_case() == ConstraintProto::ConstraintCase::kAtMostOne) {
+      constraint_literals.push_back(GetLiteralsFromSetPPCConstraint(ct));
+
+      uint64 signature = 0;
+      for (const int literal : constraint_literals.back()) {
+        const int positive_literal = PositiveRef(literal);
+        signature |= (1LL << (positive_literal % 64));
+        DCHECK_GE(positive_literal, 0);
+        if (positive_literal >= literals_to_constraints.size()) {
+          literals_to_constraints.resize(positive_literal + 1);
+        }
+        literals_to_constraints[positive_literal].push_back(
+            num_setppc_constraints);
+      }
+      signatures.push_back(signature);
+      marked_for_removal.push_back(false);
+      original_constraint_index.push_back(c);
+      num_setppc_constraints++;
+    }
+  }
+  VLOG(1) << "#setppc constraints: " << num_setppc_constraints;
+
+  // Set of constraint pairs which are already compared.
+  absl::flat_hash_set<std::pair<int, int>> compared_constraints;
+  for (const std::vector<int>& literal_to_constraints :
+       literals_to_constraints) {
+    for (int index1 = 0; index1 < literal_to_constraints.size(); ++index1) {
+      if (time_limit != nullptr && time_limit->LimitReached()) {
+        return changed;
+      }
+      const int c1 = literal_to_constraints[index1];
+      if (marked_for_removal[c1]) continue;
+      const std::vector<int>& c1_literals = constraint_literals[c1];
+      ConstraintProto* ct1 = context->working_model->mutable_constraints(
+          original_constraint_index[c1]);
+      for (int index2 = index1 + 1; index2 < literal_to_constraints.size();
+           ++index2) {
+        const int c2 = literal_to_constraints[index2];
+        if (marked_for_removal[c2]) continue;
+        if (marked_for_removal[c1]) break;
+        // TODO(user): This should not happen. Investigate.
+        if (c1 == c2) continue;
+
+        CHECK_LT(c1, c2);
+        if (gtl::ContainsKey(compared_constraints,
+                             std::pair<int, int>(c1, c2))) {
+          continue;
+        }
+        compared_constraints.insert({c1, c2});
+
+        // Hard limit on number of comparisions to avoid spending too much time
+        // here.
+        if (compared_constraints.size() >= 50000) return changed;
+
+        const bool smaller = (signatures[c1] & ~signatures[c2]) == 0;
+        const bool larger = (signatures[c2] & ~signatures[c1]) == 0;
+
+        if (!(smaller || larger)) {
+          continue;
+        }
+
+        // Check if literals in c1 is subset of literals in c2 or vice versa.
+        const std::vector<int>& c2_literals = constraint_literals[c2];
+        ConstraintProto* ct2 = context->working_model->mutable_constraints(
+            original_constraint_index[c2]);
+        // TODO(user): Try avoiding computation of set differences if
+        // possible.
+        std::vector<int> c1_minus_c2;
+        gtl::STLSetDifference(c1_literals, c2_literals, &c1_minus_c2);
+        std::vector<int> c2_minus_c1;
+        gtl::STLSetDifference(c2_literals, c1_literals, &c2_minus_c1);
+
+        if (c1_minus_c2.empty() && c2_minus_c1.empty()) {
+          if (ct1->constraint_case() == ct2->constraint_case()) {
+            marked_for_removal[c2] = true;
+            context->UpdateRuleStats("setppc: removed redundent constraints");
+          }
+        } else if (c1_minus_c2.empty()) {
+          if (ProcessSetPPCSubset(c1, c2, c2_minus_c1,
+                                  original_constraint_index,
+                                  &marked_for_removal, context)) {
+            context->UpdateConstraintVariableUsage(
+                original_constraint_index[c1]);
+            context->UpdateConstraintVariableUsage(
+                original_constraint_index[c2]);
+          }
+        } else if (c2_minus_c1.empty()) {
+          if (ProcessSetPPCSubset(c2, c1, c1_minus_c2,
+                                  original_constraint_index,
+                                  &marked_for_removal, context)) {
+            context->UpdateConstraintVariableUsage(
+                original_constraint_index[c1]);
+            context->UpdateConstraintVariableUsage(
+                original_constraint_index[c2]);
+          }
+        }
+      }
+    }
+  }
+  for (int c = 0; c < num_setppc_constraints; ++c) {
+    if (marked_for_removal[c]) {
+      ConstraintProto* ct = context->working_model->mutable_constraints(
+          original_constraint_index[c]);
+      changed = RemoveConstraint(ct, context);
+      context->UpdateConstraintVariableUsage(original_constraint_index[c]);
+    }
+  }
+
+  return changed;
+}
+
 void PresolveToFixPoint(PresolveContext* context, TimeLimit* time_limit) {
   if (context->is_unsat) return;
 
@@ -2837,6 +3045,11 @@ bool PresolveCpModel(const PresolveOptions& options,
     PresolvePureSatPart(&context);
   }
 
+  // Process set packing, partitioning and covering constraint.
+  if (options.time_limit == nullptr || !options.time_limit->LimitReached()) {
+    ProcessSetPPC(&context, options.time_limit);
+  }
+
   // Extract redundant at most one constraint form the linear ones.
   //
   // TODO(user): more generally if we do some probing, the same relation will

ortools/sat/cp_model_search.cc

@@ -281,81 +281,96 @@ SatParameters DiversifySearchParameters(const SatParameters& params,
   //   - Different propatation levels for scheduling constraints
   SatParameters new_params = params;
   new_params.set_random_seed(params.random_seed() + worker_id);
+  int index = worker_id;
+
   if (cp_model.has_objective()) {
-    switch (worker_id) {
-      case 0: {  // Use default parameters and automatic search.
-        new_params.set_search_branching(SatParameters::AUTOMATIC_SEARCH);
-        new_params.set_linearization_level(1);
-        *name = "auto";
-        break;
+    if (index == 0) {  // Use default parameters and automatic search.
+      new_params.set_search_branching(SatParameters::AUTOMATIC_SEARCH);
+      new_params.set_linearization_level(1);
+      *name = "auto";
+      return new_params;
+    }
+
+    if (cp_model.search_strategy_size() > 0) {
+      if (--index == 0) {  // Use default parameters and fixed search.
+        new_params.set_search_branching(SatParameters::FIXED_SEARCH);
+        *name = "fixed";
+        return new_params;
       }
-      case 1: {  // Use default parameters and fixed search if there is one.
-        if (cp_model.search_strategy_size() > 0) {
-          new_params.set_search_branching(SatParameters::FIXED_SEARCH);
-          *name = "fixed";
-        } else {  // Use LP_SEARCH if not.
-          new_params.set_search_branching(SatParameters::LP_SEARCH);
-          *name = "lp_br";
-        }
-        break;
-      }
-      case 2: {  // Remove LP relaxation.
-        new_params.set_search_branching(SatParameters::AUTOMATIC_SEARCH);
-        new_params.set_linearization_level(0);
-        *name = "no_lp";
-        break;
-      }
-      case 3: {  // Reinforce LP relaxation.
-        new_params.set_search_branching(SatParameters::AUTOMATIC_SEARCH);
-        new_params.set_linearization_level(2);
-        *name = "max_lp";
-        break;
-      }
-      case 4: {  // Core based approach.
+    }
+
+    // TODO(user): Disable lp_br if linear part is small or empty.
+    if (--index == 0) {
+      new_params.set_search_branching(SatParameters::LP_SEARCH);
+      *name = "lp_br";
+      return new_params;
+    }
+
+    // TODO(user): Disable no_lp if linear part is small.
+    if (--index == 0) {  // Remove LP relaxation.
+      new_params.set_search_branching(SatParameters::AUTOMATIC_SEARCH);
+      new_params.set_linearization_level(0);
+      *name = "no_lp";
+      return new_params;
+    }
+
+    // TODO(user): Disable max_lp if no change in linearization against auto.
+    if (--index == 0) {  // Reinforce LP relaxation.
+      new_params.set_search_branching(SatParameters::AUTOMATIC_SEARCH);
+      new_params.set_linearization_level(2);
+      *name = "max_lp";
+      return new_params;
+    }
+
+    if (cp_model.objective().vars_size() > 1) {
+      if (--index == 0) {  // Core based approach.
         new_params.set_search_branching(SatParameters::AUTOMATIC_SEARCH);
         new_params.set_optimize_with_core(true);
         new_params.set_linearization_level(0);
         *name = "core";
-        break;
-      }
-      default: {  // LNS.
-        new_params.set_search_branching(SatParameters::AUTOMATIC_SEARCH);
-        new_params.set_use_lns(true);
-        new_params.set_lns_num_threads(1);
-        *name = absl::StrFormat("lns_%i", worker_id - 5);
+        return new_params;
       }
     }
+
+    // Use LNS for the remaining workers.
+    new_params.set_search_branching(SatParameters::AUTOMATIC_SEARCH);
+    new_params.set_use_lns(true);
+    new_params.set_lns_num_threads(1);
+    *name = absl::StrFormat("lns_%i", index);
+    return new_params;
   } else {
     // The goal here is to try fixed and free search on the first two threads.
     // Then maximize diversity on the extra threads.
-    switch (worker_id) {
-      case 0: {
+    int index = worker_id;
+
+    if (index == 0) {  // Default automatic search.
+      new_params.set_search_branching(SatParameters::AUTOMATIC_SEARCH);
+      *name = "auto";
+      return new_params;
+    }
+
+    if (cp_model.search_strategy_size() > 0) {  // Use predefined search.
+      if (--index == 0) {
         new_params.set_search_branching(SatParameters::FIXED_SEARCH);
         *name = "fixed";
-        break;
-      }
-      case 1: {
-        new_params.set_search_branching(SatParameters::AUTOMATIC_SEARCH);
-        *name = "auto";
-        break;
-      }
-      case 2: {
-        new_params.set_search_branching(SatParameters::AUTOMATIC_SEARCH);
-        new_params.set_boolean_encoding_level(0);
-        *name = "less encoding";
-        break;
-      }
-      default: {  // Randomized fixed search.
-        new_params.set_search_branching(SatParameters::FIXED_SEARCH);
-        new_params.set_randomize_search(true);
-        new_params.set_search_randomization_tolerance(worker_id - 3);
-        *name = absl::StrFormat("rnd_%i", worker_id - 3);
-        break;
+        return new_params;
       }
     }
-  }
 
-  return new_params;
+    if (--index == 0) {  // Reduce boolean encoding.
+      new_params.set_search_branching(SatParameters::AUTOMATIC_SEARCH);
+      new_params.set_boolean_encoding_level(0);
+      *name = "less encoding";
+      return new_params;
+    }
+
+    // Randomized fixed search.
+    new_params.set_search_branching(SatParameters::FIXED_SEARCH);
+    new_params.set_randomize_search(true);
+    new_params.set_search_randomization_tolerance(index);
+    *name = absl::StrFormat("rnd_%i", index);
+    return new_params;
+  }
 }
 
 // TODO(user): Better stats in the multi thread case.

ortools/sat/cp_model_solver.cc

@@ -2057,8 +2057,20 @@ CpSolverResponse SolveCpModelParallel(
   absl::Notification first_solution_found_or_search_finished;
 
   const int num_search_workers = params.num_search_workers();
+
+  // Collect per-worker parameters and names.
+  std::vector<SatParameters> worker_parameters;
+  std::vector<std::string> worker_names;
+  for (int worker_id = 0; worker_id < num_search_workers; ++worker_id) {
+    std::string worker_name;
+    const SatParameters local_params =
+        DiversifySearchParameters(params, model_proto, worker_id, &worker_name);
+    worker_names.push_back(worker_name);
+    worker_parameters.push_back(local_params);
+  }
+
   VLOG(1) << "Starting parallel search with " << num_search_workers
-          << " workers.";
+          << " workers: " << absl::StrJoin(worker_names, ", ");
 
   if (!model_proto.has_objective()) {
     {
@@ -2066,9 +2078,9 @@ CpSolverResponse SolveCpModelParallel(
       pool.StartWorkers();
 
       for (int worker_id = 0; worker_id < num_search_workers; ++worker_id) {
-        std::string worker_name;
-        const SatParameters local_params = DiversifySearchParameters(
-            params, model_proto, worker_id, &worker_name);
+        const std::string worker_name = worker_names[worker_id];
+        const SatParameters local_params = worker_parameters[worker_id];
+
         pool.Schedule([&model_proto, stopped, local_params, &best_response,
                        &mutex, worker_name, wall_timer]() {
           Model local_model;
@@ -2130,9 +2142,8 @@ CpSolverResponse SolveCpModelParallel(
     pool.StartWorkers();
 
     for (int worker_id = 0; worker_id < num_search_workers; ++worker_id) {
-      std::string worker_name;
-      const SatParameters local_params = DiversifySearchParameters(
-          params, model_proto, worker_id, &worker_name);
+      const std::string worker_name = worker_names[worker_id];
+      const SatParameters local_params = worker_parameters[worker_id];
 
       const auto solution_observer = [maximize, worker_name, &mutex,
                                       &best_response, &observer,

ortools/sat/cuts.cc

@@ -410,6 +410,7 @@ CutGenerator CreateKnapsackCoverCutGenerator(
     const std::vector<IntegerVariable>& vars, Model* model) {
   CutGenerator result;
   result.vars = vars;
+  result.type = "Knapsack";
 
   IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
   std::vector<LinearConstraint> knapsack_constraints;

ortools/sat/cuts.h

@@ -35,6 +35,7 @@ namespace sat {
 //   their negation.
 // - Only return cuts in term of the same variables or their negation.
 struct CutGenerator {
+  std::string type;
   std::vector<IntegerVariable> vars;
   std::function<std::vector<LinearConstraint>(
       const gtl::ITIVector<IntegerVariable, double>& lp_values)>

ortools/sat/linear_programming_constraint.cc

@@ -246,6 +246,12 @@ bool LinearProgrammingConstraint::IncrementalPropagate(
     const std::vector<int>& watch_indices) {
   if (!lp_solution_is_set_) return Propagate();
 
+  // At level zero, if there is still a chance to add cuts or lazy constraints,
+  // we re-run the LP.
+  if (trail_->CurrentDecisionLevel() == 0 && !lp_at_level_zero_is_final_) {
+    return Propagate();
+  }
+
   // Check whether the change breaks the current LP solution. If it does, call
   // Propagate() on the current LP.
   for (const int index : watch_indices) {
@@ -305,6 +311,10 @@ void LinearProgrammingConstraint::UpdateBoundsOfLpVariables() {
 }
 
 bool LinearProgrammingConstraint::SolveLp() {
+  if (trail_->CurrentDecisionLevel() == 0) {
+    lp_at_level_zero_is_final_ = false;
+  }
+
   const auto status = simplex_.Solve(lp_data_, time_limit_);
   if (!status.ok()) {
     LOG(WARNING) << "The LP solver encountered an error: "
@@ -641,7 +651,8 @@ bool LinearProgrammingConstraint::Propagate() {
           std::vector<LinearConstraint> cuts =
               generator.generate_cuts(expanded_lp_solution_);
           for (const LinearConstraint& cut : cuts) {
-            constraint_manager_.AddCut(cut, "TODO", expanded_lp_solution_);
+            constraint_manager_.AddCut(cut, generator.type,
+                                       expanded_lp_solution_);
           }
         }
       }
@@ -656,6 +667,10 @@ bool LinearProgrammingConstraint::Propagate() {
                   << simplex_.GetObjectiveValue()
                   << " diff: " << simplex_.GetObjectiveValue() - old_obj;
         }
+      } else {
+        if (trail_->CurrentDecisionLevel() == 0) {
+          lp_at_level_zero_is_final_ = true;
+        }
       }
     }
   }
@@ -1279,6 +1294,7 @@ CutGenerator CreateStronglyConnectedGraphCutGenerator(
     const std::vector<IntegerVariable>& vars) {
   CutGenerator result;
   result.vars = vars;
+  result.type = "StronglyConnectedGraph";
   result.generate_cuts =
       [num_nodes, tails, heads,
        vars](const gtl::ITIVector<IntegerVariable, double>& lp_values) {
@@ -1331,6 +1347,7 @@ CutGenerator CreateCVRPCutGenerator(int num_nodes,
 
   CutGenerator result;
   result.vars = vars;
+  result.type = "CVRP";
   result.generate_cuts =
       [num_nodes, tails, heads, total_demands, demands, capacity,
        vars](const gtl::ITIVector<IntegerVariable, double>& lp_values) {

ortools/sat/linear_programming_constraint.h

@@ -323,6 +323,10 @@ class LinearProgrammingConstraint : public PropagatorInterface,
   // solution is no longer optimal though.
   std::vector<double> level_zero_lp_solution_;
 
+  // True if the last time we solved the exact same LP at level zero, no cuts
+  // and no lazy constraints where added.
+  bool lp_at_level_zero_is_final_ = false;
+
   // Same as lp_solution_ but this vector is indexed differently.
   LinearProgrammingConstraintLpSolution& expanded_lp_solution_;
 

ortools/sat/lp_utils.cc

@@ -68,6 +68,7 @@ bool ConvertMPModelProtoToCpModelProto(const MPModelProto& mp_model,
   int num_truncated_bounds = 0;
   int num_small_domains = 0;
   const int64 kSmallDomainSize = 1000;
+  const double kWantedCoefficientPrecision = 1e-6;
 
   // Add the variables.
   const int num_variables = mp_model.variable_size();
@@ -164,9 +165,22 @@ bool ConvertMPModelProtoToCpModelProto(const MPModelProto& mp_model,
       lower_bounds.push_back(var_proto.domain(0));
       upper_bounds.push_back(var_proto.domain(var_proto.domain_size() - 1));
     }
-    GetBestScalingOfDoublesToInt64(coefficients, lower_bounds, upper_bounds,
-                                   kScalingTarget, &scaling_factor,
-                                   &relative_coeff_error, &scaled_sum_error);
+    scaling_factor = GetBestScalingOfDoublesToInt64(
+        coefficients, lower_bounds, upper_bounds, kScalingTarget);
+
+    // Returns the smallest factor of the form 2^i that gives us a relative
+    // coefficient precision of kWantedCoefficientPrecision and still make sure
+    // we will have no integer overflow.
+    //
+    // TODO(user): Make this faster.
+    double x = std::min(scaling_factor, 1.0);
+    for (; x <= scaling_factor; x *= 2) {
+      ComputeScalingErrors(coefficients, lower_bounds, upper_bounds, x,
+                           &relative_coeff_error, &scaled_sum_error);
+      if (relative_coeff_error < kWantedCoefficientPrecision) break;
+    }
+    scaling_factor = x;
+
     const int64 gcd = ComputeGcdOfRoundedDoubles(coefficients, scaling_factor);
     max_relative_coeff_error =
         std::max(relative_coeff_error, max_relative_coeff_error);
@@ -186,9 +200,9 @@ bool ConvertMPModelProtoToCpModelProto(const MPModelProto& mp_model,
     // Add the constraint bounds. Because we are sure the scaled constraint fit
     // on an int64, if the scaled bounds are too large, the constraint is either
     // always true or always false.
-    const Fractional lb = mp_constraint.lower_bound();
-    const Fractional scaled_lb =
-        std::round(lb * scaling_factor - scaled_sum_error);
+    Fractional lb = mp_constraint.lower_bound();
+    lb -= std::max(1.0, std::abs(lb)) * 1e-6;
+    const Fractional scaled_lb = std::ceil(lb * scaling_factor);
     if (lb == -kInfinity || scaled_lb <= kint64min) {
       arg->add_domain(kint64min);
     } else {
@@ -196,9 +210,9 @@ bool ConvertMPModelProtoToCpModelProto(const MPModelProto& mp_model,
                                 IntegerValue(gcd))
                           .value());
     }
-    const Fractional ub = mp_constraint.upper_bound();
-    const Fractional scaled_ub =
-        std::round(ub * scaling_factor + scaled_sum_error);
+    Fractional ub = mp_constraint.upper_bound();
+    ub += std::max(1.0, std::abs(ub)) * 1e-6;
+    const Fractional scaled_ub = std::floor(ub * scaling_factor);
     if (ub == kInfinity || scaled_ub >= kint64max) {
       arg->add_domain(kint64max);
     } else {
@@ -232,9 +246,22 @@ bool ConvertMPModelProtoToCpModelProto(const MPModelProto& mp_model,
     upper_bounds.push_back(var_proto.domain(var_proto.domain_size() - 1));
   }
   if (!coefficients.empty() || mp_model.objective_offset() != 0.0) {
-    GetBestScalingOfDoublesToInt64(coefficients, lower_bounds, upper_bounds,
-                                   kScalingTarget, &scaling_factor,
-                                   &relative_coeff_error, &scaled_sum_error);
+    scaling_factor = GetBestScalingOfDoublesToInt64(
+        coefficients, lower_bounds, upper_bounds, kScalingTarget);
+
+    // Returns the smallest factor of the form 2^i that gives us a relative
+    // coefficient precision of kWantedCoefficientPrecision and still make sure
+    // we will have no integer overflow.
+    //
+    // TODO(user): Make this faster.
+    double x = std::min(scaling_factor, 1.0);
+    for (; x <= scaling_factor; x *= 2) {
+      ComputeScalingErrors(coefficients, lower_bounds, upper_bounds, x,
+                           &relative_coeff_error, &scaled_sum_error);
+      if (relative_coeff_error < kWantedCoefficientPrecision) break;
+    }
+    scaling_factor = x;
+
     const int64 gcd = ComputeGcdOfRoundedDoubles(coefficients, scaling_factor);
     max_relative_coeff_error =
         std::max(relative_coeff_error, max_relative_coeff_error);
@@ -267,10 +294,10 @@ bool ConvertMPModelProtoToCpModelProto(const MPModelProto& mp_model,
   }
 
   // Test the precision of the conversion.
-  const double kRelativeTolerance = 1e-2;
+  const double kRelativeTolerance = 1e-4;
   if (max_relative_coeff_error > kRelativeTolerance) {
     LOG(WARNING) << "The relative error during double -> int64 conversion "
-                 << "is too high!";
+                 << "is too high! error:" << max_relative_coeff_error;
     return false;
   }
   return true;

ortools/util/fp_utils.cc

@@ -20,27 +20,56 @@
 namespace operations_research {
 
 namespace {
+
 void ReorderAndCapTerms(double* min, double* max) {
   if (*min > *max) std::swap(*min, *max);
   if (*min > 0.0) *min = 0.0;
   if (*max < 0.0) *max = 0.0;
 }
-}  // namespace
+
+template <bool use_bounds>
+void ComputeScalingErrors(const std::vector<double>& input,
+                          const std::vector<double>& lb,
+                          const std::vector<double>& ub, double scaling_factor,
+                          double* max_relative_coeff_error,
+                          double* max_scaled_sum_error) {
+  const double kInfinity = std::numeric_limits<double>::infinity();
+
+  double max_positive_error = 0.0;
+  double max_negative_error = 0.0;
+  *max_relative_coeff_error = 0.0;
+  const int size = input.size();
+  for (int i = 0; i < size; ++i) {
+    const double x = input[i];
+    if (x == 0.0) continue;
+    if (use_bounds && lb[i] == 0 && ub[i] == 0) continue;
+    const double scaled = std::abs(x * scaling_factor);
+    if (scaled == 0.0) {
+      *max_relative_coeff_error = kInfinity;
+      *max_scaled_sum_error = kInfinity;
+      return;
+    }
+    *max_relative_coeff_error = std::max(
+        *max_relative_coeff_error, std::abs(std::round(scaled) / scaled - 1));
+    const double error = std::round(scaled) - scaled;
+    const double error_a = error * (use_bounds ? x * lb[i] : -x);
+    const double error_b = error * (use_bounds ? x * ub[i] : x);
+    max_positive_error += std::max(0.0, std::max(error_a, error_b));
+    max_negative_error += std::max(0.0, std::max(-error_a, -error_b));
+  }
+  *max_scaled_sum_error = std::max(max_positive_error, max_negative_error);
+}
 
 template <bool use_bounds>
 void GetBestScalingOfDoublesToInt64(const std::vector<double>& input,
                                     const std::vector<double>& lb,
                                     const std::vector<double>& ub,
                                     int64 max_absolute_sum,
-                                    double* scaling_factor,
-                                    double* max_relative_coeff_error,
-                                    double* max_scaled_sum_error) {
+                                    double* scaling_factor) {
   const double kInfinity = std::numeric_limits<double>::infinity();
 
   // We start by initializing the returns value to the "error" state.
   *scaling_factor = 0;
-  *max_relative_coeff_error = kInfinity;
-  *max_scaled_sum_error = kInfinity;
 
   // Abort in the "error" state if max_absolute_sum doesn't make sense.
   if (max_absolute_sum < 0) return;
@@ -125,37 +154,27 @@ void GetBestScalingOfDoublesToInt64(const std::vector<double>& input,
     }
   }
   *scaling_factor = ldexp(1.0, factor_exponent);
-
-  // Compute the worst case errors.
-  double max_positive_error = 0.0;
-  double max_negative_error = 0.0;
-  *max_relative_coeff_error = 0.0;
-  for (int i = 0; i < size; ++i) {
-    const double x = input[i];
-    if (x == 0.0) continue;
-    if (use_bounds && lb[i] == 0 && ub[i] == 0) continue;
-    const double scaled = std::abs(ldexp(x, factor_exponent));
-    *max_relative_coeff_error = std::max(
-        *max_relative_coeff_error, std::abs(std::round(scaled) / scaled - 1));
-    const double error = std::round(scaled) - scaled;
-    const double error_a = error * (use_bounds ? x * lb[i] : -x);
-    const double error_b = error * (use_bounds ? x * ub[i] : x);
-    max_positive_error += std::max(0.0, std::max(error_a, error_b));
-    max_negative_error += std::max(0.0, std::max(-error_a, -error_b));
-  }
-  *max_scaled_sum_error = std::max(max_positive_error, max_negative_error);
 }
 
-void GetBestScalingOfDoublesToInt64(const std::vector<double>& input,
-                                    const std::vector<double>& lb,
-                                    const std::vector<double>& ub,
-                                    int64 max_absolute_sum,
-                                    double* scaling_factor,
-                                    double* max_relative_coeff_error,
-                                    double* max_scaled_sum_error) {
-  return GetBestScalingOfDoublesToInt64<true>(
-      input, lb, ub, max_absolute_sum, scaling_factor, max_relative_coeff_error,
-      max_scaled_sum_error);
+}  // namespace
+
+void ComputeScalingErrors(const std::vector<double>& input,
+                          const std::vector<double>& lb,
+                          const std::vector<double>& ub, double scaling_factor,
+                          double* max_relative_coeff_error,
+                          double* max_scaled_sum_error) {
+  ComputeScalingErrors<true>(input, lb, ub, scaling_factor,
+                             max_relative_coeff_error, max_scaled_sum_error);
+}
+
+double GetBestScalingOfDoublesToInt64(const std::vector<double>& input,
+                                      const std::vector<double>& lb,
+                                      const std::vector<double>& ub,
+                                      int64 max_absolute_sum) {
+  double scaling_factor;
+  GetBestScalingOfDoublesToInt64<true>(input, lb, ub, max_absolute_sum,
+                                       &scaling_factor);
+  return scaling_factor;
 }
 
 void GetBestScalingOfDoublesToInt64(const std::vector<double>& input,
@@ -163,9 +182,10 @@ void GetBestScalingOfDoublesToInt64(const std::vector<double>& input,
                                     double* scaling_factor,
                                     double* max_relative_coeff_error) {
   double max_scaled_sum_error;
-  return GetBestScalingOfDoublesToInt64<false>(
-      input, {}, {}, max_absolute_sum, scaling_factor, max_relative_coeff_error,
-      &max_scaled_sum_error);
+  GetBestScalingOfDoublesToInt64<false>(input, {}, {}, max_absolute_sum,
+                                        scaling_factor);
+  ComputeScalingErrors<false>(input, {}, {}, *scaling_factor,
+                              max_relative_coeff_error, &max_scaled_sum_error);
 }
 
 int64 ComputeGcdOfRoundedDoubles(const std::vector<double>& x,

ortools/util/fp_utils.h

@@ -204,22 +204,28 @@ void GetBestScalingOfDoublesToInt64(const std::vector<double>& input,
                                     double* scaling_factor,
                                     double* max_relative_coeff_error);
 
-// Same as the function above, but enforces that
+// Returns the scaling factor like above with the extra conditions:
 //  -  The sum over i of min(0, round(factor * x[i])) >= -max_sum.
 //  -  The sum over i of max(0, round(factor * x[i])) <= max_sum.
 // For any possible values of the x[i] such that x[i] is in [lb[i], ub[i]].
+double GetBestScalingOfDoublesToInt64(const std::vector<double>& input,
+                                      const std::vector<double>& lb,
+                                      const std::vector<double>& ub,
+                                      int64 max_absolute_sum);
+// This computes:
 //
-// This also computes the max_scaled_sum_error which is a bound on the maximum
-// difference between the exact scaled sum and the rounded one. One needs to
-// divide this by scaling_factor to have the maximum absolute error on the
-// original sum.
-void GetBestScalingOfDoublesToInt64(const std::vector<double>& input,
-                                    const std::vector<double>& lb,
-                                    const std::vector<double>& ub,
-                                    int64 max_absolute_sum,
-                                    double* scaling_factor,
-                                    double* max_relative_coeff_error,
-                                    double* max_scaled_sum_error);
+// The max_relative_coeff_error, which is the maximum over all coeff of
+// |round(factor * x[i]) / (factor * x[i])  - 1|.
+//
+// The max_scaled_sum_error which is a bound on the maximum difference between
+// the exact scaled sum and the rounded one. One needs to divide this by
+// scaling_factor to have the maximum absolute error on the original sum.
+void ComputeScalingErrors(const std::vector<double>& input,
+                          const std::vector<double>& lb,
+                          const std::vector<double>& ub,
+                          const double scaling_factor,
+                          double* max_relative_coeff_error,
+                          double* max_scaled_sum_error);
 
 // Returns the Greatest Common Divisor of the numbers
 // round(fabs(x[i] * scaling_factor)). The numbers 0 are ignored and if they are