From 6bb8272b2ad5459543dcd6d31f88617ba4e42c09 Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Mon, 16 Dec 2019 15:33:52 +0100 Subject: [PATCH] [CP-SAT] improve rounding methods in linear cuts --- ortools/sat/cuts.cc | 70 ++++++++++++++------ ortools/sat/cuts.h | 7 +- ortools/sat/linear_programming_constraint.cc | 1 - ortools/sat/python/cp_model.py | 2 + ortools/sat/sat_parameters.proto | 8 +-- 5 files changed, 57 insertions(+), 31 deletions(-) diff --git a/ortools/sat/cuts.cc b/ortools/sat/cuts.cc index a4cb6c2bc0..a6bd4b51c7 100644 --- a/ortools/sat/cuts.cc +++ b/ortools/sat/cuts.cc @@ -592,28 +592,35 @@ CutGenerator CreateKnapsackCoverCutGenerator( } std::function GetSuperAdditiveRoundingFunction( - bool use_letchford_lodi_version, IntegerValue rhs_remainder, - IntegerValue divisor, IntegerValue max_scaling) { - // Compute the larger t <= max_scaling such that + IntegerValue rhs_remainder, IntegerValue divisor, IntegerValue max_t, + IntegerValue max_scaling) { + CHECK_GE(max_t, 1); + CHECK_GE(max_scaling, 1); + + // Compute the larger t <= max_t such that // t * rhs_remainder >= divisor / 2. const IntegerValue t = rhs_remainder == 0 - ? max_scaling - : std::min(max_scaling, CeilRatio(divisor / 2, rhs_remainder)); + ? max_t + : std::min(max_t, CeilRatio(divisor / 2, rhs_remainder)); // Adjust after the multiplication by t. + // + // Note(user): the modulo is only needed when t is large which is currently + // not possible, but I left it here to not forget to do that in experiments. rhs_remainder *= t; - max_scaling /= t; + rhs_remainder %= divisor; - // This is the only difference compared to a discretized MIR function. - if (use_letchford_lodi_version && max_scaling > 2) max_scaling = 2; + // Make sure we don't have an integer overflow below. Note that we assume that + // divisor and the maximum coeff magnitude are not too different (maybe a + // factor 1000 at most) so that the final result will never overflow. + max_scaling = std::min(max_scaling, kint64max / divisor); - CHECK_GE(max_scaling, 1); const IntegerValue size = divisor - rhs_remainder; - if (max_scaling == 1) { + if (max_scaling == 1 || size == 1) { // TODO(user): Use everywhere a two step computation to avoid overflow? - // First divide by divisor, then multiply by t. For now, we limit the - // max_scaling so that we never have an overflow instead. + // First divide by divisor, then multiply by t. For now, we limit t so that + // we never have an overflow instead. return [t, divisor](IntegerValue coeff) { return FloorRatio(t * coeff, divisor); }; @@ -624,6 +631,22 @@ std::function GetSuperAdditiveRoundingFunction( const IntegerValue diff = remainder - rhs_remainder; return size * ratio + std::max(IntegerValue(0), diff); }; + } else if (max_scaling.value() * rhs_remainder.value() < divisor) { + // Because of our max_t limitation, the rhs_remainder might stay small. + // + // If it is "too small" we cannot use the code below because it will not be + // valid. So we just divide divisor into max_scaling bucket. The + // rhs_remainder will be in the bucket 0. + // + // Note(user): This seems the same as just increasing t, modulo integer + // overflows. Maybe we should just always do the computation like this so + // that we can use larger t even if coeff is close to kint64max. + return [t, divisor, max_scaling](IntegerValue coeff) { + const IntegerValue ratio = FloorRatio(t * coeff, divisor); + const IntegerValue remainder = t * coeff - ratio * divisor; + const IntegerValue bucket = FloorRatio(remainder * max_scaling, divisor); + return max_scaling * ratio + bucket; + }; } else { // We divide (size = divisor - rhs_remainder) into (max_scaling - 1) buckets // and increase the function by 1 / max_scaling for each of them. @@ -632,7 +655,7 @@ std::function GetSuperAdditiveRoundingFunction( // functions that do not dominate each others. So potentially, a max scaling // as low as 2 could lead to the better cut (this is exactly the Letchford & // Lodi function). - /// + // // Another intersting fact, is that if we want to compute the maximum alpha // for a constraint with 2 terms like: // divisor * Y + (ratio * divisor + remainder) * X @@ -777,7 +800,7 @@ void IntegerRoundingCut(RoundingOptions options, // we need to include it too since we will apply the rounding function on it. max_magnitude = std::max(max_magnitude, IntTypeAbs(cut->ub)); - // Make sure that when we multiply the rhs or the coefficient by max_scaling, + // Make sure that when we multiply the rhs or the coefficient by a factor t, // we do not have an integer overflow. Actually, we need a bit more room // because we might round down a value to the next multiple of // max_magnitude. @@ -787,10 +810,9 @@ void IntegerRoundingCut(RoundingOptions options, *cut = LinearConstraint(IntegerValue(0), IntegerValue(0)); return; } - const IntegerValue max_scaling = - std::min(options.max_scaling, threshold / max_magnitude); + const IntegerValue max_t = threshold / max_magnitude; - // There is no point trying twice the same divisor or s divisor that is too + // There is no point trying twice the same divisor or a divisor that is too // small. Note that we use a higher threshold than the remainder_threshold // because we can boost the remainder thanks to our adjusting heuristic below // and also because this allows to have cuts with a small range of @@ -879,8 +901,8 @@ void IntegerRoundingCut(RoundingOptions options, temp_ub - FloorRatio(temp_ub, divisor) * divisor; if (rhs_remainder == 0) continue; - const auto f = GetSuperAdditiveRoundingFunction( - !options.use_mir, rhs_remainder, divisor, max_scaling); + const auto f = GetSuperAdditiveRoundingFunction(rhs_remainder, divisor, + max_t, options.max_scaling); // As we round coefficients, we will compute the loss compared to the // current scaled constraint activity. As soon as this loss crosses the @@ -967,10 +989,16 @@ void IntegerRoundingCut(RoundingOptions options, } // Create the super-additive function f(). + // + // TODO(user): Try out different rounding function and keep the best. We can + // change max_t and max_scaling. It might not be easy to choose which cut is + // the best, but we can at least know for sure if one dominate the other + // completely. That is, if for all coeff f(coeff)/f(divisor) is greater than + // or equal to the same value for another function f. const IntegerValue rhs_remainder = cut->ub - FloorRatio(cut->ub, best_divisor) * best_divisor; - const auto f = GetSuperAdditiveRoundingFunction( - !options.use_mir, rhs_remainder, best_divisor, max_scaling); + const auto f = GetSuperAdditiveRoundingFunction(rhs_remainder, best_divisor, + max_t, options.max_scaling); // Apply f() to the cut. // diff --git a/ortools/sat/cuts.h b/ortools/sat/cuts.h index 7b3089e1e0..272552c7e0 100644 --- a/ortools/sat/cuts.h +++ b/ortools/sat/cuts.h @@ -96,7 +96,7 @@ class ImpliedBoundsProcessor { // // Algorithm: // - We first scale by a factor t so that rhs_remainder >= divisor / 2. -// - Then, if use_letchford_lodi_version is true, we use the function described +// - Then, if max_scaling == 2, we use the function described // in "Strenghtening Chvatal-Gomory cuts and Gomory fractional cuts", Adam N. // Letchfrod, Andrea Lodi. // - Otherwise, we use a generalization of this which is a discretized version @@ -109,8 +109,8 @@ class ImpliedBoundsProcessor { // it could be nice to try to generate a cut using different values of // max_scaling. std::function GetSuperAdditiveRoundingFunction( - bool use_letchford_lodi_version, IntegerValue rhs_remainder, - IntegerValue divisor, IntegerValue max_scaling); + IntegerValue rhs_remainder, IntegerValue divisor, IntegerValue max_t, + 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. @@ -142,7 +142,6 @@ std::function GetSuperAdditiveRoundingFunction( // 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. struct RoundingOptions { - bool use_mir = false; IntegerValue max_scaling = IntegerValue(60); }; void IntegerRoundingCut(RoundingOptions options, diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc index e33479cb2c..a1ebaf28e6 100644 --- a/ortools/sat/linear_programming_constraint.cc +++ b/ortools/sat/linear_programming_constraint.cc @@ -670,7 +670,6 @@ void LinearProgrammingConstraint::AddCutFromConstraints( // Get the cut using some integer rounding heuristic. RoundingOptions options; - options.use_mir = sat_parameters_.use_mir_rounding(); options.max_scaling = sat_parameters_.max_integer_rounding_scaling(); IntegerRoundingCut(options, lp_values, var_lbs, var_ubs, &cut); diff --git a/ortools/sat/python/cp_model.py b/ortools/sat/python/cp_model.py index ba7b800725..f7c0ac9982 100644 --- a/ortools/sat/python/cp_model.py +++ b/ortools/sat/python/cp_model.py @@ -1529,6 +1529,7 @@ def EvaluateLinearExpr(expression, solution): if not isinstance(expression, LinearExpr): raise TypeError('Cannot interpret %s as a linear expression.' % expression) + value = 0 to_process = [(expression, 1)] while to_process: @@ -1754,6 +1755,7 @@ class CpSolverSolutionCallback(pywrapsat.SolutionCallback): if not isinstance(expression, LinearExpr): raise TypeError('Cannot interpret %s as a linear expression.' % expression) + value = 0 to_process = [(expression, 1)] while to_process: diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto index 86dfc3660d..d706e3abfe 100644 --- a/ortools/sat/sat_parameters.proto +++ b/ortools/sat/sat_parameters.proto @@ -565,16 +565,14 @@ message SatParameters { // have a cut generator. optional int32 max_all_diff_cut_size = 148 [default = 7]; - // 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]; - // In the integer rounding procedure used for MIR and Gomory cut, the maximum // "scaling" we use (must be positive). The lower this is, the lower the // integer coefficients of the cut will be. Note that cut generated by lower // values are not necessarily worse than cut generated by larger value. There // is no strict dominance relationship. + // + // Setting this to 2 result in the "strong fractional rouding" of Letchford + // and Lodi. optional int32 max_integer_rounding_scaling = 119 [default = 600]; // If true, we start by an empty LP, and only add constraints not satisfied