diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index 1452cf41d1..9efddc3f16 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/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 tmp_literals; std::vector tmp_term_domains; std::vector tmp_left_domains; + absl::flat_hash_set tmp_literal_set; // Each time a domain is modified this is set to true. SparseBitset 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 GetLiteralsFromSetPPCConstraint(ConstraintProto* ct) { + std::vector 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& c2_minus_c1, + const std::vector& original_constraint_index, + std::vector* 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 signatures; + + // Graph of constraints to literals. constraint_literals[c] contains all the + // literals in constraint indexed by 'c' in sorted order. + std::vector> 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> 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 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 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> compared_constraints; + for (const std::vector& 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& 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(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& 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 c1_minus_c2; + gtl::STLSetDifference(c1_literals, c2_literals, &c1_minus_c2); + std::vector 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 diff --git a/ortools/sat/cp_model_search.cc b/ortools/sat/cp_model_search.cc index a65bb81727..f49cfc8999 100644 --- a/ortools/sat/cp_model_search.cc +++ b/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. diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index 823367b3ec..215021e9ca 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/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 worker_parameters; + std::vector 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, diff --git a/ortools/sat/cuts.cc b/ortools/sat/cuts.cc index 2e7b503c2c..d0c3a33189 100644 --- a/ortools/sat/cuts.cc +++ b/ortools/sat/cuts.cc @@ -410,6 +410,7 @@ CutGenerator CreateKnapsackCoverCutGenerator( const std::vector& vars, Model* model) { CutGenerator result; result.vars = vars; + result.type = "Knapsack"; IntegerTrail* integer_trail = model->GetOrCreate(); std::vector knapsack_constraints; diff --git a/ortools/sat/cuts.h b/ortools/sat/cuts.h index 80fd473ecb..9da61d955c 100644 --- a/ortools/sat/cuts.h +++ b/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 vars; std::function( const gtl::ITIVector& lp_values)> diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc index deaa263129..29e334b722 100644 --- a/ortools/sat/linear_programming_constraint.cc +++ b/ortools/sat/linear_programming_constraint.cc @@ -246,6 +246,12 @@ bool LinearProgrammingConstraint::IncrementalPropagate( const std::vector& 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 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& vars) { CutGenerator result; result.vars = vars; + result.type = "StronglyConnectedGraph"; result.generate_cuts = [num_nodes, tails, heads, vars](const gtl::ITIVector& 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& lp_values) { diff --git a/ortools/sat/linear_programming_constraint.h b/ortools/sat/linear_programming_constraint.h index 2be08907ae..267ec09329 100644 --- a/ortools/sat/linear_programming_constraint.h +++ b/ortools/sat/linear_programming_constraint.h @@ -323,6 +323,10 @@ class LinearProgrammingConstraint : public PropagatorInterface, // solution is no longer optimal though. std::vector 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_; diff --git a/ortools/sat/lp_utils.cc b/ortools/sat/lp_utils.cc index 64fb903838..a51de8386d 100644 --- a/ortools/sat/lp_utils.cc +++ b/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; diff --git a/ortools/util/fp_utils.cc b/ortools/util/fp_utils.cc index a3484d62e0..f7867a5621 100644 --- a/ortools/util/fp_utils.cc +++ b/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 +void ComputeScalingErrors(const std::vector& input, + const std::vector& lb, + const std::vector& ub, double scaling_factor, + double* max_relative_coeff_error, + double* max_scaled_sum_error) { + const double kInfinity = std::numeric_limits::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 void GetBestScalingOfDoublesToInt64(const std::vector& input, const std::vector& lb, const std::vector& 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::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& 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& input, - const std::vector& lb, - const std::vector& ub, - int64 max_absolute_sum, - double* scaling_factor, - double* max_relative_coeff_error, - double* max_scaled_sum_error) { - return GetBestScalingOfDoublesToInt64( - input, lb, ub, max_absolute_sum, scaling_factor, max_relative_coeff_error, - max_scaled_sum_error); +} // namespace + +void ComputeScalingErrors(const std::vector& input, + const std::vector& lb, + const std::vector& ub, double scaling_factor, + double* max_relative_coeff_error, + double* max_scaled_sum_error) { + ComputeScalingErrors(input, lb, ub, scaling_factor, + max_relative_coeff_error, max_scaled_sum_error); +} + +double GetBestScalingOfDoublesToInt64(const std::vector& input, + const std::vector& lb, + const std::vector& ub, + int64 max_absolute_sum) { + double scaling_factor; + GetBestScalingOfDoublesToInt64(input, lb, ub, max_absolute_sum, + &scaling_factor); + return scaling_factor; } void GetBestScalingOfDoublesToInt64(const std::vector& input, @@ -163,9 +182,10 @@ void GetBestScalingOfDoublesToInt64(const std::vector& input, double* scaling_factor, double* max_relative_coeff_error) { double max_scaled_sum_error; - return GetBestScalingOfDoublesToInt64( - input, {}, {}, max_absolute_sum, scaling_factor, max_relative_coeff_error, - &max_scaled_sum_error); + GetBestScalingOfDoublesToInt64(input, {}, {}, max_absolute_sum, + scaling_factor); + ComputeScalingErrors(input, {}, {}, *scaling_factor, + max_relative_coeff_error, &max_scaled_sum_error); } int64 ComputeGcdOfRoundedDoubles(const std::vector& x, diff --git a/ortools/util/fp_utils.h b/ortools/util/fp_utils.h index cf239c99fa..569eeddf47 100644 --- a/ortools/util/fp_utils.h +++ b/ortools/util/fp_utils.h @@ -204,22 +204,28 @@ void GetBestScalingOfDoublesToInt64(const std::vector& 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& input, + const std::vector& lb, + const std::vector& 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& input, - const std::vector& lb, - const std::vector& 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& input, + const std::vector& lb, + const std::vector& 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