diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index 9efddc3f16..bf9d59f1f8 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -688,19 +688,31 @@ bool PresolveIntProd(ConstraintProto* ct, PresolveContext* context) { if (ct->int_prod().vars_size() == 2) { int a = ct->int_prod().vars(0); int b = ct->int_prod().vars(1); - const int p = ct->int_prod().target(); + const int product = ct->int_prod().target(); if (context->IsFixed(b)) std::swap(a, b); if (context->IsFixed(a)) { - ConstraintProto* const lin = context->working_model->add_constraints(); - lin->mutable_linear()->add_vars(b); - lin->mutable_linear()->add_coeffs(context->MinOf(a)); - lin->mutable_linear()->add_vars(p); - lin->mutable_linear()->add_coeffs(-1); - lin->mutable_linear()->add_domain(0); - lin->mutable_linear()->add_domain(0); + if (b != product) { + ConstraintProto* const lin = context->working_model->add_constraints(); + lin->mutable_linear()->add_vars(b); + lin->mutable_linear()->add_coeffs(context->MinOf(a)); + lin->mutable_linear()->add_vars(product); + lin->mutable_linear()->add_coeffs(-1); + lin->mutable_linear()->add_domain(0); + lin->mutable_linear()->add_domain(0); - context->UpdateRuleStats("int_prod: linearize product by constant."); + context->UpdateRuleStats("int_prod: linearize product by constant."); + return RemoveConstraint(ct, context); + } else if (context->MinOf(a) != 1) { + context->IntersectDomainWith(product, Domain(0, 0)); + context->UpdateRuleStats("int_prod: fix variable to zero."); + } else { + context->UpdateRuleStats("int_prod: remove identity."); + return RemoveConstraint(ct, context); + } + } else if (a == b && a == product) { // x = x * x, only true for {0, 1}. + context->IntersectDomainWith(product, Domain(0, 1)); + context->UpdateRuleStats("int_prod: fix variable to zero or one."); return RemoveConstraint(ct, context); } } diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index db37fecf70..317d1d48c6 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -2079,7 +2079,7 @@ CpSolverResponse SolveCpModelParallel( } VLOG(1) << "Starting parallel search with " << num_search_workers - << " workers: " << absl::StrJoin(worker_names, ", "); + << " workers: [" << absl::StrJoin(worker_names, ", ") << "]"; if (!model_proto.has_objective()) { { diff --git a/ortools/sat/cuts.cc b/ortools/sat/cuts.cc index d0c3a33189..e742846abf 100644 --- a/ortools/sat/cuts.cc +++ b/ortools/sat/cuts.cc @@ -590,7 +590,7 @@ std::function GetSuperAdditiveRoundingFunction( const IntegerValue t = std::max( IntegerValue(1), remainder == 0 ? IntegerValue(1) - : std::min(max_scaling, CeilRatio(target, remainder))); + : std::min(max_scaling / 2, CeilRatio(target, remainder))); const IntegerValue threshold = std::max(target, t * remainder); return [t, threshold, divisor](IntegerValue coeff) { const IntegerValue ratio = FloorRatio(t * coeff, divisor); @@ -599,7 +599,33 @@ std::function GetSuperAdditiveRoundingFunction( }; } -void IntegerRoundingCut(std::vector lp_values, +// TODO(user): if 2 * rhs_remainder < divisor, multiply by a factor t before +// rounding. +std::function GetMirFunction( + IntegerValue rhs_remainder, IntegerValue divisor, + IntegerValue max_scaling) { + CHECK_GE(max_scaling, 1); + if (divisor - rhs_remainder <= max_scaling) { + return [rhs_remainder, divisor](IntegerValue coeff) { + const IntegerValue ratio = FloorRatio(coeff, divisor); + const IntegerValue remainder = coeff - ratio * divisor; + return (divisor - rhs_remainder) * ratio + + std::max(IntegerValue(0), remainder - rhs_remainder); + }; + } else { + // TODO(user): This function is not maximal, improve? + return [rhs_remainder, divisor, max_scaling](IntegerValue coeff) { + const IntegerValue ratio = FloorRatio(coeff, divisor); + const IntegerValue remainder = coeff - ratio * divisor; + return max_scaling * ratio + + std::max(FloorRatio((remainder - rhs_remainder) * max_scaling, + divisor - rhs_remainder), + IntegerValue(0)); + }; + } +} + +void IntegerRoundingCut(RoundingOptions options, std::vector lp_values, std::vector lower_bounds, std::vector upper_bounds, LinearConstraint* cut) { @@ -730,9 +756,11 @@ void IntegerRoundingCut(std::vector lp_values, // Create the super-additive function f(). const IntegerValue rhs_remainder = cut->ub - FloorRatio(cut->ub, divisor) * divisor; - const IntegerValue kMaxScaling(4); const auto f = - GetSuperAdditiveRoundingFunction(rhs_remainder, divisor, kMaxScaling); + options.use_mir + ? GetMirFunction(rhs_remainder, divisor, options.max_scaling) + : GetSuperAdditiveRoundingFunction(rhs_remainder, divisor, + options.max_scaling); // Apply f() to the cut. cut->ub = f(cut->ub); diff --git a/ortools/sat/cuts.h b/ortools/sat/cuts.h index 9da61d955c..57190899bf 100644 --- a/ortools/sat/cuts.h +++ b/ortools/sat/cuts.h @@ -65,7 +65,11 @@ struct CutGenerator { // // TODO(user): Still not clear to me if this can be improved or not. std::function GetSuperAdditiveRoundingFunction( - IntegerValue remainder, IntegerValue divisor, IntegerValue max_scaling); + IntegerValue rhs_remainder, IntegerValue divisor, IntegerValue max_scaling); + +// Same as GetSuperAdditiveRoundingFunction() but uses the classic MIR one. +std::function GetMirFunction( + IntegerValue rhs_remainder, IntegerValue divisor, IntegerValue max_scaling); // Given an upper bounded linear constraint, this function tries to transform it // to a valid cut that violate the given LP solution using integer rounding. @@ -96,7 +100,11 @@ std::function GetSuperAdditiveRoundingFunction( // more effort tunning them. In particular, one can try many heuristics and keep // the best looking cut (or more than one). This is not on the critical code // path, so we can spend more effort in finding good cuts. -void IntegerRoundingCut(std::vector lp_values, +struct RoundingOptions { + bool use_mir = false; + IntegerValue max_scaling = IntegerValue(1000); +}; +void IntegerRoundingCut(RoundingOptions options, std::vector lp_values, std::vector lower_bounds, std::vector upper_bounds, LinearConstraint* cut); diff --git a/ortools/sat/linear_constraint_manager.cc b/ortools/sat/linear_constraint_manager.cc index 4f430bd756..93e0921493 100644 --- a/ortools/sat/linear_constraint_manager.cc +++ b/ortools/sat/linear_constraint_manager.cc @@ -142,17 +142,23 @@ void LinearConstraintManager::AddCut( if (VLOG_IS_ON(1)) { IntegerValue max_magnitude(0); double activity = 0.0; + double l2_norm = 0.0; const int size = ct.vars.size(); for (int i = 0; i < size; ++i) { max_magnitude = std::max(max_magnitude, IntTypeAbs(ct.coeffs[i])); - activity += ToDouble(ct.coeffs[i]) * lp_solution[ct.vars[i]]; + + const double coeff = ToDouble(ct.coeffs[i]); + activity += coeff * lp_solution[ct.vars[i]]; + l2_norm += coeff * coeff; } + l2_norm = sqrt(l2_norm); double violation = 0.0; violation = std::max(violation, activity - ToDouble(ct.ub)); violation = std::max(violation, ToDouble(ct.lb) - activity); - LOG(INFO) << "Cut '" << type_name << "' size: " << size - << " max_magnitude: " << max_magnitude - << " violation: " << violation; + LOG(INFO) << "Cut '" << type_name << "'" + << " size=" << size << " max_magnitude=" << max_magnitude + << " norm=" << l2_norm << " violation=" << violation + << " eff=" << violation / l2_norm; } } diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc index 29e334b722..decaff3683 100644 --- a/ortools/sat/linear_programming_constraint.cc +++ b/ortools/sat/linear_programming_constraint.cc @@ -502,7 +502,9 @@ void LinearProgrammingConstraint::AddCGCuts() { } // Get the cut using some integer rounding heuristic. - IntegerRoundingCut(lp_values, var_lbs, var_ubs, &cut); + RoundingOptions options; + options.use_mir = sat_parameters_.use_mir_rounding(); + IntegerRoundingCut(options, lp_values, var_lbs, var_ubs, &cut); // Compute the activity. Warning: the cut no longer have the same size so we // cannot use lp_values. Note that the substitution below shouldn't change diff --git a/ortools/sat/lp_utils.cc b/ortools/sat/lp_utils.cc index a51de8386d..8cc3ffa257 100644 --- a/ortools/sat/lp_utils.cc +++ b/ortools/sat/lp_utils.cc @@ -201,7 +201,7 @@ bool ConvertMPModelProtoToCpModelProto(const MPModelProto& mp_model, // on an int64, if the scaled bounds are too large, the constraint is either // always true or always false. Fractional lb = mp_constraint.lower_bound(); - lb -= std::max(1.0, std::abs(lb)) * 1e-6; + lb -= std::max(1.0, std::abs(lb)) * 1e-8; const Fractional scaled_lb = std::ceil(lb * scaling_factor); if (lb == -kInfinity || scaled_lb <= kint64min) { arg->add_domain(kint64min); @@ -211,7 +211,7 @@ bool ConvertMPModelProtoToCpModelProto(const MPModelProto& mp_model, .value()); } Fractional ub = mp_constraint.upper_bound(); - ub += std::max(1.0, std::abs(ub)) * 1e-6; + ub += std::max(1.0, std::abs(ub)) * 1e-8; const Fractional scaled_ub = std::floor(ub * scaling_factor); if (ub == kInfinity || scaled_ub >= kint64max) { arg->add_domain(kint64max); diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto index f3229dddff..0ec4bc27ee 100644 --- a/ortools/sat/sat_parameters.proto +++ b/ortools/sat/sat_parameters.proto @@ -21,7 +21,7 @@ package operations_research.sat; // Contains the definitions for all the sat algorithm parameters and their // default values. // -// NEXT TAG: 118 +// NEXT TAG: 119 message SatParameters { // ========================================================================== // Branching and polarity @@ -507,6 +507,11 @@ message SatParameters { // of enum that list enabled cuts, like "cuts:[KNAPSACK,GOMORY,CIRCUIT]". optional bool add_cg_cuts = 117 [default = false]; + // Whether we use the classical MIR rounding function when computing a cut + // using integer rounding. If false we use our own variant based on the strong + // fractional rouding of Letchford and Lodi. + optional bool use_mir_rounding = 118 [default = true]; + // If true, we start by an empty LP, and only add constraints not satisfied // by the current LP solution batch by batch. A constraint that is only added // like this is known as a "lazy" constraint in the literature, except that we