always-cautious/ortools-clone
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
73e4f754456f512edc450286dbe6929f4abeb5d8
ortools-clone/ortools/sat/cuts.cc

2389 lines
92 KiB
C++
Raw Blame History

// Copyright 2010-2022 Google LLC
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include "ortools/sat/cuts.h"
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <functional>
#include <limits>
#include <memory>
#include <string>
#include <utility>
#include <vector>
#include "absl/container/btree_set.h"
#include "absl/container/flat_hash_map.h"
#include "absl/container/flat_hash_set.h"
#include "ortools/base/logging.h"
#include "ortools/base/stl_util.h"
#include "ortools/base/strong_vector.h"
#include "ortools/sat/clause.h"
#include "ortools/sat/implied_bounds.h"
#include "ortools/sat/integer.h"
#include "ortools/sat/linear_constraint.h"
#include "ortools/sat/linear_constraint_manager.h"
#include "ortools/sat/model.h"
#include "ortools/sat/sat_base.h"
#include "ortools/util/saturated_arithmetic.h"
#include "ortools/util/sorted_interval_list.h"
#include "ortools/util/strong_integers.h"
#include "ortools/util/time_limit.h"
namespace operations_research {
namespace sat {
std::string CutTerm::DebugString() const {
return absl::StrCat("coeff=", coeff.value(), " lp=", lp_value,
" range=", bound_diff.value());
}
bool CutTerm::Complement(IntegerValue* rhs) {
// We replace coeff * X by coeff * (X - bound_diff + bound_diff)
// which gives -coeff * complement(X) + coeff * bound_diff;
if (!AddProductTo(-coeff, bound_diff, rhs)) return false;
// We keep the same expression variable.
for (int i = 0; i < 2; ++i) {
expr_coeffs[i] = -expr_coeffs[i];
}
expr_offset = bound_diff - expr_offset;
// Note that this is not involutive because of floating point error. Fix?
lp_value = ToDouble(bound_diff) - lp_value;
coeff = -coeff;
return true;
}
// To try to minimize the risk of overflow, we switch to the bound closer
// to the lp_value. Since most of our base constraint for cut are tight,
// hopefully this is not too bad.
bool CutData::AppendOneTerm(IntegerVariable var, IntegerValue coeff,
double lp_value, IntegerValue lb, IntegerValue ub) {
if (coeff == 0) return true;
const IntegerValue bound_diff = ub - lb;
// Complement the variable so that it is always closer to its lb.
bool complement = false;
const double lb_dist = std::abs(lp_value - ToDouble(lb));
const double ub_dist = std::abs(lp_value - ToDouble(ub));
if (ub_dist < lb_dist) {
complement = true;
}
// See formula below, the constant term is either coeff * lb or coeff * ub.
if (!AddProductTo(-coeff, complement ? ub : lb, &rhs)) {
return false;
}
// Deal with fixed variable, no need to shift back in this case, we can
// just remove the term.
if (bound_diff == 0) return true;
CutTerm entry;
entry.expr_vars[0] = var;
entry.expr_coeffs[1] = 0;
entry.bound_diff = bound_diff;
if (complement) {
// X = -(UB - X) + UB
entry.expr_coeffs[0] = -IntegerValue(1);
entry.expr_offset = ub;
entry.coeff = -coeff;
entry.lp_value = ToDouble(ub) - lp_value;
} else {
// C = (X - LB) + LB
entry.expr_coeffs[0] = IntegerValue(1);
entry.expr_offset = -lb;
entry.coeff = coeff;
entry.lp_value = lp_value - ToDouble(lb);
}
terms.push_back(entry);
return true;
}
bool CutData::FillFromLinearConstraint(
const LinearConstraint& base_ct,
const absl::StrongVector<IntegerVariable, double>& lp_values,
IntegerTrail* integer_trail) {
rhs = base_ct.ub;
terms.clear();
const int num_terms = base_ct.vars.size();
for (int i = 0; i < num_terms; ++i) {
const IntegerVariable var = base_ct.vars[i];
if (!AppendOneTerm(var, base_ct.coeffs[i], lp_values[base_ct.vars[i]],
integer_trail->LevelZeroLowerBound(var),
integer_trail->LevelZeroUpperBound(var))) {
return false;
}
}
return true;
}
bool CutData::FillFromParallelVectors(
const LinearConstraint& base_ct, const std::vector<double>& lp_values,
const std::vector<IntegerValue>& lower_bounds,
const std::vector<IntegerValue>& upper_bounds) {
rhs = base_ct.ub;
terms.clear();
const int size = lp_values.size();
if (size == 0) return true;
CHECK_EQ(lower_bounds.size(), size);
CHECK_EQ(upper_bounds.size(), size);
CHECK_EQ(base_ct.vars.size(), size);
CHECK_EQ(base_ct.coeffs.size(), size);
CHECK_EQ(base_ct.lb, kMinIntegerValue);
for (int i = 0; i < size; ++i) {
if (!AppendOneTerm(base_ct.vars[i], base_ct.coeffs[i], lp_values[i],
lower_bounds[i], upper_bounds[i])) {
return false;
}
}
return true;
}
void CutData::Canonicalize() {
num_relevant_entries = 0;
max_magnitude = IntTypeAbs(rhs);
for (int i = 0; i < terms.size(); ++i) {
CutTerm& entry = terms[i];
max_magnitude = std::max(max_magnitude, IntTypeAbs(entry.coeff));
if (entry.HasRelevantLpValue()) {
std::swap(terms[num_relevant_entries], entry);
++num_relevant_entries;
}
}
// Sort by larger lp_value first.
std::sort(terms.begin(), terms.begin() + num_relevant_entries,
[](const CutTerm& a, const CutTerm& b) {
return a.lp_value > b.lp_value;
});
}
void CutDataBuilder::ClearIndices() {
num_merges_ = 0;
direct_index_.clear();
complemented_index_.clear();
}
void CutDataBuilder::RegisterAllBooleansTerms(const CutData& cut) {
ClearIndices();
const int size = cut.terms.size();
for (int i = 0; i < size; ++i) {
const CutTerm& term = cut.terms[i];
if (term.bound_diff != 1) continue;
if (!term.IsSimple()) continue;
if (term.expr_coeffs[0] > 0) {
direct_index_[term.expr_vars[0]] = i;
} else {
complemented_index_[term.expr_vars[0]] = i;
}
}
}
void CutDataBuilder::AddOrMergeTerm(const CutTerm& term, IntegerValue t,
CutData* cut) {
DCHECK(term.IsSimple());
const IntegerVariable var = term.expr_vars[0];
const int new_index = cut->terms.size();
const auto [it, inserted] =
term.expr_coeffs[0] > 0 ? direct_index_.insert({var, new_index})
: complemented_index_.insert({var, new_index});
const int entry_index = it->second;
if (inserted) {
cut->terms.push_back(term);
} else {
// We can only merge the term if term.coeff + old_coeff do not overflow and
// if t * new_coeff do not overflow.
//
// If we cannot merge the term, we will keep them separate. The produced cut
// will be less strong, but can still be used.
const int64_t new_coeff =
CapAdd(cut->terms[entry_index].coeff.value(), term.coeff.value());
const int64_t overflow_check = CapProd(t.value(), new_coeff);
if (AtMinOrMaxInt64(new_coeff) || AtMinOrMaxInt64(overflow_check)) {
// If we cannot merge the term, we keep them separate.
cut->terms.push_back(term);
} else {
++num_merges_;
cut->terms[entry_index].coeff = IntegerValue(new_coeff);
}
}
}
bool CutDataBuilder::ConvertToLinearConstraint(const CutData& cut,
LinearConstraint* output) {
tmp_map_.clear();
IntegerValue new_rhs = cut.rhs;
for (const CutTerm& term : cut.terms) {
for (int i = 0; i < 2; ++i) {
if (term.expr_coeffs[i] == 0) continue;
if (!AddProductTo(term.coeff, term.expr_coeffs[i],
&tmp_map_[term.expr_vars[i]])) {
return false;
}
}
if (!AddProductTo(-term.coeff, term.expr_offset, &new_rhs)) {
return false;
}
}
output->ClearTerms();
output->lb = kMinIntegerValue;
output->ub = new_rhs;
for (const auto [var, coeff] : tmp_map_) {
if (coeff == 0) continue;
output->vars.push_back(var);
output->coeffs.push_back(coeff);
}
DivideByGCD(output);
return true;
}
namespace {
// Minimum amount of violation of the cut constraint by the solution. This
// is needed to avoid numerical issues and adding cuts with minor effect.
const double kMinCutViolation = 1e-4;
IntegerValue CapProdI(IntegerValue a, IntegerValue b) {
return IntegerValue(CapProd(a.value(), b.value()));
}
IntegerValue CapSubI(IntegerValue a, IntegerValue b) {
return IntegerValue(CapSub(a.value(), b.value()));
}
IntegerValue CapAddI(IntegerValue a, IntegerValue b) {
return IntegerValue(CapAdd(a.value(), b.value()));
}
bool ProdOverflow(IntegerValue t, IntegerValue value) {
return AtMinOrMaxInt64(CapProd(t.value(), value.value()));
}
} // namespace
// Compute the larger t <= max_t such that t * rhs_remainder >= divisor / 2.
//
// This is just a separate function as it is slightly faster to compute the
// result only once.
IntegerValue GetFactorT(IntegerValue rhs_remainder, IntegerValue divisor,
IntegerValue max_magnitude) {
// Make sure that when we multiply the rhs or the coefficient by a factor t,
// we do not have an integer overflow. Note that the rhs should be counted
// in max_magnitude since we will apply f() on it.
IntegerValue max_t(std::numeric_limits<int64_t>::max());
if (max_magnitude != 0) {
max_t = max_t / max_magnitude;
}
return rhs_remainder == 0
? max_t
: std::min(max_t, CeilRatio(divisor / 2, rhs_remainder));
}
std::function<IntegerValue(IntegerValue)> GetSuperAdditiveRoundingFunction(
IntegerValue rhs_remainder, IntegerValue divisor, IntegerValue t,
IntegerValue max_scaling) {
DCHECK_GE(max_scaling, 1);
DCHECK_GE(t, 1);
// Adjust after the multiplication by t.
rhs_remainder *= t;
DCHECK_LT(rhs_remainder, divisor);
// 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, std::numeric_limits<int64_t>::max() / divisor);
const IntegerValue size = divisor - rhs_remainder;
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 t so that
// we never have an overflow instead.
return [t, divisor](IntegerValue coeff) {
return FloorRatio(t * coeff, divisor);
};
} else if (size <= max_scaling) {
return [size, rhs_remainder, t, divisor](IntegerValue coeff) {
const IntegerValue t_coeff = t * coeff;
const IntegerValue ratio = FloorRatio(t_coeff, divisor);
const IntegerValue remainder = PositiveRemainder(t_coeff, divisor);
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 t_coeff = t * coeff;
const IntegerValue ratio = FloorRatio(t_coeff, divisor);
const IntegerValue remainder = PositiveRemainder(t_coeff, 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.
//
// Note that for different values of max_scaling, we get a family of
// 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 interesting fact, is that if we want to compute the maximum alpha
// for a constraint with 2 terms like:
// divisor * Y + (ratio * divisor + remainder) * X
// <= rhs_ratio * divisor + rhs_remainder
// so that we have the cut:
// Y + (ratio + alpha) * X <= rhs_ratio
// This is the same as computing the maximum alpha such that for all integer
// X > 0 we have CeilRatio(alpha * divisor * X, divisor)
// <= CeilRatio(remainder * X - rhs_remainder, divisor).
// We can prove that this alpha is of the form (n - 1) / n, and it will
// be reached by such function for a max_scaling of n.
//
// TODO(user): This function is not always maximal when
// size % (max_scaling - 1) == 0. Improve?
return [size, rhs_remainder, t, divisor, max_scaling](IntegerValue coeff) {
const IntegerValue t_coeff = t * coeff;
const IntegerValue ratio = FloorRatio(t_coeff, divisor);
const IntegerValue remainder = PositiveRemainder(t_coeff, divisor);
const IntegerValue diff = remainder - rhs_remainder;
const IntegerValue bucket =
diff > 0 ? CeilRatio(diff * (max_scaling - 1), size)
: IntegerValue(0);
return max_scaling * ratio + bucket;
};
}
}
IntegerRoundingCutHelper::~IntegerRoundingCutHelper() {
if (!VLOG_IS_ON(1)) return;
if (shared_stats_ == nullptr) return;
std::vector<std::pair<std::string, int64_t>> stats;
stats.push_back({"rounding_cut/num_initial_ibs_", total_num_initial_ibs_});
stats.push_back(
{"rounding_cut/num_initial_merges_", total_num_initial_merges_});
stats.push_back({"rounding_cut/num_pos_lifts", total_num_pos_lifts_});
stats.push_back({"rounding_cut/num_neg_lifts", total_num_neg_lifts_});
stats.push_back(
{"rounding_cut/num_post_complements", total_num_post_complements_});
stats.push_back({"rounding_cut/num_overflows", total_num_overflow_abort_});
stats.push_back({"rounding_cut/num_adjusts", total_num_coeff_adjust_});
stats.push_back({"rounding_cut/num_merges", total_num_merges_});
stats.push_back({"rounding_cut/num_bumps", total_num_bumps_});
stats.push_back(
{"rounding_cut/num_final_complements", total_num_final_complements_});
stats.push_back({"rounding_cut/num_dominating_f", total_num_dominating_f_});
shared_stats_->AddStats(stats);
}
double IntegerRoundingCutHelper::GetScaledViolation(
IntegerValue divisor, IntegerValue max_scaling,
IntegerValue remainder_threshold, const CutData& cut) {
IntegerValue rhs = cut.rhs;
IntegerValue max_magnitude = cut.max_magnitude;
const IntegerValue initial_rhs_remainder = PositiveRemainder(rhs, divisor);
if (initial_rhs_remainder < remainder_threshold) return 0.0;
// We will adjust coefficient that are just under an exact multiple of
// divisor to an exact multiple. This is meant to get rid of small errors
// that appears due to rounding error in our exact computation of the
// initial constraint given to this class.
//
// Each adjustement will cause the initial_rhs_remainder to increase, and we
// do not want to increase it above divisor. Our threshold below guarantees
// this. Note that the higher the rhs_remainder becomes, the more the
// function f() has a chance to reduce the violation, so it is not always a
// good idea to use all the slack we have between initial_rhs_remainder and
// divisor.
//
// TODO(user): We could see if for a fixed function f, the increase is
// interesting?
// before: f(rhs) - f(coeff) * lp_value
// after: f(rhs + increase * bound_diff) - f(coeff + increase) * lp_value.
adjusted_coeffs_.clear();
const IntegerValue adjust_threshold =
(divisor - initial_rhs_remainder - 1) /
IntegerValue(std::max(1000, cut.num_relevant_entries));
if (adjust_threshold > 0) {
// Even before we finish the adjust, we can have a lower bound on the
// activily loss using this divisor, and so we can abort early. This is
// similar to what is done below.
double max_violation = ToDouble(initial_rhs_remainder);
for (int i = 0; i < cut.num_relevant_entries; ++i) {
const CutTerm& entry = cut.terms[i];
const IntegerValue remainder = PositiveRemainder(entry.coeff, divisor);
if (remainder == 0) continue;
if (remainder <= initial_rhs_remainder) {
// We do not know exactly f() yet, but it will always round to the
// floor of the division by divisor in this case.
max_violation -= ToDouble(remainder) * entry.lp_value;
if (max_violation <= 1e-3) return 0.0;
continue;
}
// Adjust coeff of the form k * divisor - epsilon.
const IntegerValue adjust = divisor - remainder;
const IntegerValue prod = CapProdI(adjust, entry.bound_diff);
if (prod <= adjust_threshold) {
rhs += prod;
const IntegerValue new_coeff = entry.coeff + adjust;
adjusted_coeffs_.push_back({i, new_coeff});
max_magnitude = std::max(max_magnitude, IntTypeAbs(new_coeff));
}
}
}
max_magnitude = std::max(max_magnitude, IntTypeAbs(rhs));
const IntegerValue rhs_remainder = PositiveRemainder(rhs, divisor);
const IntegerValue t = GetFactorT(rhs_remainder, divisor, max_magnitude);
const auto f =
GetSuperAdditiveRoundingFunction(rhs_remainder, divisor, t, 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
// slack, then we known that there is no violation and we can abort early.
//
// TODO(user): modulo the scaling, we could compute the exact threshold
// using our current best cut. Note that we also have to account the change
// in slack due to the adjust code above.
const double scaling = ToDouble(f(divisor)) / ToDouble(divisor);
double max_violation = scaling * ToDouble(rhs_remainder);
// Apply f() to the cut and compute the cut violation. Note that it is
// okay to just look at the relevant indices since the other have a lp
// value which is almost zero. Doing it like this is faster, and even if
// the max_magnitude might be off it should still be relevant enough.
double violation = -ToDouble(f(rhs));
double l2_norm = 0.0;
int adjusted_coeffs_index = 0;
for (int i = 0; i < cut.num_relevant_entries; ++i) {
const CutTerm& entry = cut.terms[i];
// Adjust coeff according to our previous computation if needed.
IntegerValue coeff = entry.coeff;
if (adjusted_coeffs_index < adjusted_coeffs_.size() &&
adjusted_coeffs_[adjusted_coeffs_index].first == i) {
coeff = adjusted_coeffs_[adjusted_coeffs_index].second;
adjusted_coeffs_index++;
}
if (coeff == 0) continue;
const IntegerValue new_coeff = f(coeff);
const double new_coeff_double = ToDouble(new_coeff);
const double lp_value = entry.lp_value;
// TODO(user): Shall we compute the norm after slack are substituted back?
// it might be widely different. Another reason why this might not be
// the best measure.
l2_norm += new_coeff_double * new_coeff_double;
violation += new_coeff_double * lp_value;
max_violation -= (scaling * ToDouble(coeff) - new_coeff_double) * lp_value;
if (max_violation <= 1e-3) return 0.0;
}
if (l2_norm == 0.0) return 0.0;
// Here we scale by the L2 norm over the "relevant" positions. This seems
// to work slighly better in practice.
//
// Note(user): The non-relevant position have an LP value of zero. If their
// coefficient is positive, it seems good not to take it into account in the
// norm since the larger this coeff is, the stronger the cut. If the coeff
// is negative though, a large coeff means a small increase from zero of the
// lp value will make the cut satisfied, so we might want to look at them.
return violation / sqrt(l2_norm);
}
bool IntegerRoundingCutHelper::HasComplementedImpliedBound(
const CutTerm& entry, ImpliedBoundsProcessor* ib_processor) {
if (ib_processor == nullptr) return false;
if (!entry.IsSimple()) return false;
if (entry.bound_diff == 1) return false;
const ImpliedBoundsProcessor::BestImpliedBoundInfo info =
ib_processor->GetCachedImpliedBoundInfo(
entry.expr_coeffs[0] > 0 ? NegationOf(entry.expr_vars[0])
: entry.expr_vars[0]);
return info.bool_var != kNoIntegerVariable;
}
// TODO(user): This is slow, 50% of run time on a2c1s1.pb.gz. Optimize!
bool IntegerRoundingCutHelper::ComputeCut(
RoundingOptions options, const CutData& base_ct,
ImpliedBoundsProcessor* ib_processor) {
// Try IB before heuristic?
// This should be better except it can mess up the norm and the divisors.
best_cut_ = base_ct;
if (options.use_ib_before_heuristic && ib_processor != nullptr) {
cut_builder_.RegisterAllBooleansTerms(best_cut_);
const int old_size = static_cast<int>(best_cut_.terms.size());
bool abort = true;
for (int i = 0; i < old_size; ++i) {
if (best_cut_.terms[i].bound_diff <= 1) continue;
if (!best_cut_.terms[i].HasRelevantLpValue()) continue;
if (options.prefer_positive_ib && best_cut_.terms[i].coeff < 0) {
// We complement the term before trying the implied bound.
if (best_cut_.terms[i].Complement(&best_cut_.rhs)) {
if (ib_processor->TryToExpandWithLowerImpliedbound(
IntegerValue(1), i,
/*complement=*/true, &best_cut_, &cut_builder_)) {
++total_num_initial_ibs_;
abort = false;
continue;
}
best_cut_.terms[i].Complement(&best_cut_.rhs);
}
}
if (ib_processor->TryToExpandWithLowerImpliedbound(
IntegerValue(1), i,
/*complement=*/true, &best_cut_, &cut_builder_)) {
abort = false;
++total_num_initial_ibs_;
}
}
total_num_initial_merges_ += cut_builder_.NumMergesSinceLastClear();
// TODO(user): We assume that this is called with and without the option
// use_ib_before_heuristic, so that we can abort if no IB has been applied
// since then we will redo the computation. This is not really clean.
if (abort) return false;
}
// Our heuristic will try to generate a few different cuts, and we will keep
// the most violated one scaled by the l2 norm of the relevant position.
//
// TODO(user): Experiment for the best value of this initial violation
// threshold. Note also that we use the l2 norm on the restricted position
// here. Maybe we should change that? On that note, the L2 norm usage seems
// a bit weird to me since it grows with the number of term in the cut. And
// often, we already have a good cut, and we make it stronger by adding
// extra terms that do not change its activity.
//
// The discussion above only concern the best_scaled_violation initial
// value. The remainder_threshold allows to not consider cuts for which the
// final efficacity is clearly lower than 1e-3 (it is a bound, so we could
// generate cuts with a lower efficacity than this).
//
// TODO(user): If the rhs is small and close to zero, we might want to
// consider different way of complementing the variables.
best_cut_.Canonicalize();
const IntegerValue remainder_threshold(
std::max(IntegerValue(1), best_cut_.max_magnitude / 1000));
if (best_cut_.rhs >= 0 && best_cut_.rhs < remainder_threshold) {
return false;
}
// 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
// coefficients.
divisors_.clear();
for (const CutTerm& entry : best_cut_.terms) {
// Note that because of the slacks, initial coeff are here too.
const IntegerValue magnitude = IntTypeAbs(entry.coeff);
if (magnitude <= remainder_threshold) continue;
divisors_.push_back(magnitude);
}
if (divisors_.empty()) return false;
gtl::STLSortAndRemoveDuplicates(&divisors_, std::greater<IntegerValue>());
// Note that most of the time is spend here since we call this function on
// many linear equation, and just a few of them have a good enough scaled
// violation. We can spend more time afterwards to tune the cut.
//
// TODO(user): Avoid quadratic algorithm? Note that we are quadratic in
// relevant positions not the full cut size, but this is still too much on
// some problems.
IntegerValue best_divisor(0);
double best_scaled_violation = 1e-3;
for (const IntegerValue divisor : divisors_) {
// Note that the function will abort right away if PositiveRemainder() is
// not good enough, so it is quick for bad divisor.
const double violation = GetScaledViolation(divisor, options.max_scaling,
remainder_threshold, best_cut_);
if (violation > best_scaled_violation) {
best_scaled_violation = violation;
best_adjusted_coeffs_ = adjusted_coeffs_;
best_divisor = divisor;
}
}
if (best_divisor == 0) return false;
// Try best_divisor divided by small number.
for (int div = 2; div < 9; ++div) {
const IntegerValue divisor = best_divisor / IntegerValue(div);
if (divisor <= 1) continue;
const double violation = GetScaledViolation(divisor, options.max_scaling,
remainder_threshold, best_cut_);
if (violation > best_scaled_violation) {
best_scaled_violation = violation;
best_adjusted_coeffs_ = adjusted_coeffs_;
best_divisor = divisor;
}
}
// Re try complementation on the transformed cut.
for (CutTerm& entry : best_cut_.terms) {
if (!entry.HasRelevantLpValue()) break;
if (entry.coeff % best_divisor == 0) continue;
// Temporary try to complement this variable.
if (!entry.Complement(&best_cut_.rhs)) continue;
const double violation = GetScaledViolation(
best_divisor, options.max_scaling, remainder_threshold, best_cut_);
if (violation > best_scaled_violation) {
// keep the change.
++total_num_post_complements_;
best_scaled_violation = violation;
best_adjusted_coeffs_ = adjusted_coeffs_;
} else {
// Restore.
entry.Complement(&best_cut_.rhs);
}
}
// Adjust coefficients as computed by the best GetScaledViolation().
for (const auto [index, new_coeff] : best_adjusted_coeffs_) {
++total_num_coeff_adjust_;
CutTerm& entry = best_cut_.terms[index];
const IntegerValue remainder = new_coeff - entry.coeff;
CHECK_GT(remainder, 0);
entry.coeff = new_coeff;
best_cut_.rhs += remainder * entry.bound_diff;
best_cut_.max_magnitude =
std::max(best_cut_.max_magnitude, IntTypeAbs(new_coeff));
}
best_cut_.max_magnitude =
std::max(best_cut_.max_magnitude, IntTypeAbs(best_cut_.rhs));
// Create the base super-additive function f().
const IntegerValue rhs_remainder =
PositiveRemainder(best_cut_.rhs, best_divisor);
IntegerValue factor_t =
GetFactorT(rhs_remainder, best_divisor, best_cut_.max_magnitude);
auto f = GetSuperAdditiveRoundingFunction(rhs_remainder, best_divisor,
factor_t, options.max_scaling);
// Look amongst all our possible function f() for one that dominate greedily
// our current best one. Note that we prefer lower scaling factor since that
// result in a cut with lower coefficients.
//
// We only look at relevant position and ignore the other. Not sure this is
// the best approach.
remainders_.clear();
for (const CutTerm& entry : best_cut_.terms) {
if (!entry.HasRelevantLpValue()) break;
const IntegerValue coeff = entry.coeff;
const IntegerValue r = PositiveRemainder(coeff, best_divisor);
if (r > rhs_remainder) remainders_.push_back(r);
}
gtl::STLSortAndRemoveDuplicates(&remainders_);
if (remainders_.size() <= 100) {
best_rs_.clear();
for (const IntegerValue r : remainders_) {
best_rs_.push_back(f(r));
}
IntegerValue best_d = f(best_divisor);
// Note that the complexity seems high 100 * 2 * options.max_scaling, but
// this only run on cuts that are already efficient and the inner loop tend
// to abort quickly. I didn't see this code in the cpu profile so far.
for (const IntegerValue t :
{IntegerValue(1),
GetFactorT(rhs_remainder, best_divisor, best_cut_.max_magnitude)}) {
for (IntegerValue s(2); s <= options.max_scaling; ++s) {
const auto g =
GetSuperAdditiveRoundingFunction(rhs_remainder, best_divisor, t, s);
int num_strictly_better = 0;
rs_.clear();
const IntegerValue d = g(best_divisor);
for (int i = 0; i < best_rs_.size(); ++i) {
const IntegerValue temp = g(remainders_[i]);
if (temp * best_d < best_rs_[i] * d) break;
if (temp * best_d > best_rs_[i] * d) num_strictly_better++;
rs_.push_back(temp);
}
if (rs_.size() == best_rs_.size() && num_strictly_better > 0) {
++total_num_dominating_f_;
f = g;
factor_t = t;
best_rs_ = rs_;
best_d = d;
}
}
}
}
// Use implied bounds to "lift" Booleans into the cut.
// This should lead to stronger cuts even if the norms migth be worse.
num_ib_used_ = 0;
if (ib_processor != nullptr) {
cut_builder_.RegisterAllBooleansTerms(best_cut_);
const int old_size = best_cut_.terms.size();
for (int i = 0; i < old_size; ++i) {
CutTerm& term = best_cut_.terms[i];
// We only want to expand non-Boolean and non-slack term!
if (term.bound_diff <= 1) continue;
if (!term.IsSimple()) continue;
if (ib_processor->TryToExpandWithLowerImpliedbound(
factor_t, i, /*complement=*/false, &best_cut_, &cut_builder_)) {
++num_ib_used_;
++total_num_pos_lifts_;
continue;
}
// Use the implied bound on (-X) if it is beneficial to do so.
// Like complementing, this is not always good.
//
// We have comp(X) = diff - X = diff * B + S
// X = diff * (1 - B) - S.
// So if we applies f, we will get:
// f(coeff * diff) * (1 - B) + f(-coeff) * S
// and substituing S = diff * (1 - B) - X, we get:
// -f(-coeff) * X + [f(coeff * diff) + f(-coeff) * diff] (1 - B).
//
// TODO(user): Note that while the violation might be higher, if the slack
// becomes large this will result in a less powerfull cut. Shall we do
// that? It is a bit the same problematic with complementing.
//
// TODO(user): If the slack is close to zero, then this transformation
// will always increase the violation. So we could potentially do it in
// Before our divisor selection heuristic. But the norm of the final cut
// will increase too.
if (!HasComplementedImpliedBound(term, ib_processor)) continue;
const ImpliedBoundsProcessor::BestImpliedBoundInfo info =
ib_processor->GetCachedImpliedBoundInfo(
term.expr_coeffs[0] > 0 ? NegationOf(term.expr_vars[0])
: term.expr_vars[0]);
const IntegerValue lb = -term.expr_offset;
const IntegerValue bound_diff = info.implied_bound - lb;
// We do not want overflow when computing f().
if (ProdOverflow(factor_t, CapProdI(term.coeff, bound_diff))) {
continue;
}
// We only consider IB that span the full range here.
if (bound_diff != term.bound_diff) continue;
// Note that -f(-coeff) >= f(coeff) but coeff_b <= 0.
const IntegerValue coeff_b =
f(term.coeff * bound_diff) + f(-term.coeff) * bound_diff;
CHECK_LE(coeff_b, 0);
const double lp1 = ToDouble(f(term.coeff)) * term.lp_value;
const double lp2 = -ToDouble(f(-term.coeff)) * term.lp_value +
ToDouble(coeff_b) * (1 - info.bool_lp_value);
if (lp2 > lp1 + 1e-2) {
// Create the Boolean term for (1 - B) in X = diff * (1 - B) - S
// We reverse the is_positive meaning here since we have (1 - B).
CutTerm bool_term;
bool_term.coeff = bound_diff * term.coeff;
bool_term.expr_vars[0] = info.bool_var;
bool_term.expr_coeffs[1] = 0;
bool_term.bound_diff = IntegerValue(1);
bool_term.lp_value = 1.0 - info.bool_lp_value;
if (!info.is_positive) {
bool_term.expr_coeffs[0] = IntegerValue(1);
bool_term.expr_offset = IntegerValue(0);
} else {
bool_term.expr_coeffs[0] = IntegerValue(-1);
bool_term.expr_offset = IntegerValue(1);
}
// Create the slack term in X = diff * (1 - B) - S
CutTerm slack_term;
slack_term.coeff = -term.coeff;
slack_term.expr_vars[0] = term.expr_vars[0];
slack_term.expr_coeffs[0] = -term.expr_coeffs[0];
slack_term.expr_vars[1] = bool_term.expr_vars[0];
slack_term.expr_coeffs[1] = bound_diff * bool_term.expr_coeffs[0];
slack_term.expr_offset =
bound_diff * bool_term.expr_offset - term.expr_offset;
slack_term.lp_value = info.SlackLpValue(lb);
slack_term.bound_diff = term.bound_diff;
// Commit the change.
++num_ib_used_;
++total_num_neg_lifts_;
term = slack_term;
cut_builder_.AddOrMergeTerm(bool_term, factor_t, &best_cut_);
}
}
total_num_merges_ += cut_builder_.NumMergesSinceLastClear();
}
// More complementation, but for the same f.
// If we can do that, it probably means our heuristics above are not great.
for (int i = 0; i < 3; ++i) {
const int64_t saved = total_num_final_complements_;
for (CutTerm& entry : best_cut_.terms) {
// Complementing an entry gives:
// [a * X <= b] -> [-a * (diff - X) <= b - a * diff]
//
// We will compare what happen when we apply f:
// [f(b) - f(a) * lp(X)] -> [f(b - a * diff) - f(-a) * (diff - lp(X))].
//
// If lp(X) is zero, then the transformation is always worse.
// Because f(b - a * diff) >= f(b) + f(-a) * diff by super-additivity.
//
// However the larger is X, the better it gets since at diff, we have
// f(b) >= f(b - a * diff) + f(a * diff) >= f(b - a * diff) + f(a) * diff.
//
// TODO(user): It is still unclear if we have a * X + b * (1 - X) <= rhs
// for a Boolean X, what is the best way to apply f and if we should merge
// the terms. If there is no other terms, best is probably
// f(rhs - a) * X + f(rhs - b) * (1 - X).
if (entry.coeff % best_divisor == 0) continue;
if (!entry.HasRelevantLpValue()) continue;
// Avoid potential overflow here.
const IntegerValue prod(CapProdI(entry.bound_diff, entry.coeff));
if (ProdOverflow(factor_t, prod)) continue;
if (ProdOverflow(factor_t, CapSubI(best_cut_.rhs, prod))) continue;
const double lp1 = ToDouble(f(best_cut_.rhs)) -
ToDouble(f(entry.coeff)) * entry.lp_value;
const double lp2 = ToDouble(f(best_cut_.rhs - prod)) -
ToDouble(f(-entry.coeff)) *
(ToDouble(entry.bound_diff) - entry.lp_value);
if (lp2 + 1e-2 < lp1) {
if (!entry.Complement(&best_cut_.rhs)) continue;
++total_num_final_complements_;
}
}
if (total_num_final_complements_ == saved) break;
}
// Apply f() to the best_cut_ with a potential improvement for one Boolean:
//
// If we have a Boolean X, and a cut: terms + a * X <= b;
// By setting X to true or false, we have two inequalities:
// terms <= b if X == 0
// terms <= b - a if X == 1
// We can apply f to both inequalities and recombine:
// f(terms) <= f(b) * (1 - X) + f(b - a) * X
// Which change the final coeff of X from f(a) to [f(b) - f(b - a)].
// This can only improve the cut since f(b) >= f(b - a) + f(a)
//
// Note that we re-Canonicalize after our possible complementation so that the
// "improvement" is applied to larger lp_value first.
best_cut_.Canonicalize();
bool improved = false;
const IntegerValue rhs = best_cut_.rhs;
const IntegerValue f_rhs = f(best_cut_.rhs);
best_cut_.rhs = f_rhs;
for (CutTerm& entry : best_cut_.terms) {
const IntegerValue f_coeff = f(entry.coeff);
if (!improved && entry.bound_diff == 1 &&
!ProdOverflow(factor_t, CapSubI(rhs, entry.coeff))) {
const IntegerValue alternative = f_rhs - f(rhs - entry.coeff);
DCHECK_GE(alternative, f_coeff);
if (alternative > f_coeff) {
++total_num_bumps_;
improved = true;
entry.coeff = alternative;
continue;
}
}
entry.coeff = f_coeff;
}
if (!cut_builder_.ConvertToLinearConstraint(best_cut_, &cut_)) {
++total_num_overflow_abort_;
return false;
}
return true;
}
CoverCutHelper::~CoverCutHelper() {
if (!VLOG_IS_ON(1)) return;
if (shared_stats_ == nullptr) return;
std::vector<std::pair<std::string, int64_t>> stats;
stats.push_back({"cover_cut/num_overflows", total_num_overflow_abort_});
stats.push_back({"cover_cut/num_lifting", total_num_lifting_});
stats.push_back({"cover_cut/num_implied_bounds", total_num_ibs_});
shared_stats_->AddStats(stats);
}
// Try a simple cover heuristic.
// Look for violated CUT of the form: sum (UB - X) or (X - LB) >= 1.
int CoverCutHelper::GetCoverSize(int relevant_size, IntegerValue* rhs) {
relevant_size =
std::partition(
base_ct_.terms.begin(), base_ct_.terms.begin() + relevant_size,
[](const CutTerm& t) { return t.LpDistToMaxValue() < 0.9999; }) -
base_ct_.terms.begin();
std::sort(base_ct_.terms.begin(), base_ct_.terms.begin() + relevant_size,
[](const CutTerm& a, const CutTerm& b) {
const double dist_a = a.LpDistToMaxValue();
const double dist_b = b.LpDistToMaxValue();
if (dist_a == dist_b) {
// Prefer low coefficients if the distance is the same.
return a.coeff < b.coeff;
}
return dist_a < dist_b;
});
double activity = 0.0;
int cover_size = relevant_size;
*rhs = base_ct_.rhs;
for (int i = 0; i < relevant_size; ++i) {
const CutTerm& term = base_ct_.terms[i];
activity += term.LpDistToMaxValue();
// As an heuristic we select all the term so that the sum of distance
// to the upper bound is <= 1.0. If the corresponding rhs is negative, then
// we will have a cut of violation at least 0.0. Note that this violation
// can be improved by the lifting.
//
// TODO(user): experiment with different threshold (even greater than one).
// Or come up with an algo that incorporate the lifting into the heuristic.
if (activity > 1.0) {
cover_size = i; // before this entry.
break;
}
if (!AddProductTo(-term.coeff, term.bound_diff, rhs)) {
// Abort early if we run into overflow.
// In that case, rhs must be negative, and we can try this cover still.
cover_size = i;
CHECK_LT(*rhs, 0);
break;
}
}
// If the rhs is now negative, we have a cut.
//
// Note(user): past this point, now that a given "base" cover has been chosen,
// we basically compute the cut (of the form sum X <= bound) with the maximum
// possible violation. Note also that we lift as much as possible, so we don't
// necessarily optimize for the cut efficacity though. But we do get a
// stronger cut.
if (*rhs >= 0) return 0;
if (cover_size == 0) return 0;
// Transform to a minimal cover. We want to greedily remove the largest coeff
// first, so we have more chance for the "lifting" below which can increase
// the cut violation. If the coeff are the same, we prefer to remove high
// distance from upper bound first.
std::sort(base_ct_.terms.begin(), base_ct_.terms.begin() + cover_size,
[](const CutTerm& a, const CutTerm& b) {
if (a.coeff == b.coeff) {
return a.LpDistToMaxValue() > b.LpDistToMaxValue();
}
return a.coeff > b.coeff;
});
for (int i = 0; i < cover_size; ++i) {
const CutTerm& t = base_ct_.terms[i];
if (t.bound_diff * t.coeff + *rhs >= 0) continue;
*rhs += t.bound_diff * t.coeff;
std::swap(base_ct_.terms[i], base_ct_.terms[--cover_size]);
}
CHECK_GT(cover_size, 0);
return cover_size;
}
bool CoverCutHelper::MakeAllTermsPositive() {
// Make sure each coeff is positive.
//
// TODO(user): maybe we should do it all at once to avoid some overflow
// condition.
for (CutTerm& term : base_ct_.terms) {
if (term.coeff >= 0) continue;
if (!term.Complement(&base_ct_.rhs)) {
++total_num_overflow_abort_;
return false;
}
}
// We should have aborted early if the base constraint was already infeasible.
CHECK_GE(base_ct_.rhs, 0);
return true;
}
bool CoverCutHelper::TrySimpleKnapsack(const CutData& input,
ImpliedBoundsProcessor* ib_processor) {
cut_.Clear();
base_ct_ = input;
if (!MakeAllTermsPositive()) return false;
if (ib_processor != nullptr) {
cut_builder_.RegisterAllBooleansTerms(base_ct_);
const int old_size = static_cast<int>(base_ct_.terms.size());
for (int i = 0; i < old_size; ++i) {
// We only look at non-Boolean with an lp value not close to the upper
// bound.
const CutTerm& term = base_ct_.terms[i];
if (term.bound_diff <= 1) continue;
if (term.lp_value + 1e-4 > static_cast<double>(term.bound_diff.value())) {
continue;
}
if (ib_processor->TryToExpandWithLowerImpliedbound(
IntegerValue(1), i,
/*complement=*/false, &base_ct_, &cut_builder_)) {
++total_num_ibs_;
}
}
}
IntegerValue rhs;
const int base_size = static_cast<int>(base_ct_.terms.size());
const int cover_size = GetCoverSize(base_size, &rhs);
if (cover_size == 0) return false;
// The cut is just that the sum of variable cannot be at their max value.
base_ct_.rhs = IntegerValue(-1);
IntegerValue max_coeff(0);
for (int i = 0; i < cover_size; ++i) {
max_coeff = std::max(max_coeff, base_ct_.terms[i].coeff);
base_ct_.terms[i].coeff = IntegerValue(1);
base_ct_.rhs += base_ct_.terms[i].bound_diff;
}
CHECK_GT(max_coeff, 0);
// In case the max_coeff variable is not binary, it might be possible to
// tighten the cut a bit more.
//
// Note(user): I never observed this on the miplib so far.
if (max_coeff < -rhs) {
const IntegerValue m = FloorRatio(-rhs - 1, max_coeff);
rhs += max_coeff * m;
base_ct_.rhs -= m;
}
CHECK_LT(rhs, 0);
IntegerValue max_base_magnitude = max_coeff;
max_base_magnitude = std::max(max_base_magnitude, IntTypeAbs(base_ct_.rhs));
for (int i = cover_size; i < base_size; ++i) {
max_base_magnitude = std::max(max_base_magnitude, base_ct_.terms[i].coeff);
}
const IntegerValue max_scaling(std::min(
IntegerValue(60), FloorRatio(kMaxIntegerValue, max_base_magnitude)));
// Lift all at once the variables not used in the cover.
//
// We have a cut of the form sum_i X_i <= b that we will lift into
// sum_i scaling X_i + sum f(base_coeff_j) X_j <= b * scaling.
//
// Using the super additivity of f() and how we construct it, for all N >= 0
// we know that: sum_j base_coeff_j X_j <= N * max_coeff + (max_coeff - slack)
// implies that: sum_j f(base_coeff_j) X_j <= N * scaling.
// So by inverting the implication we have:
// 1/ lift > N * scaling => lift_sum > N * max_coeff + (max_coeff - slack)
// We also have:
// 2/ cut > b -(N+1) => original sum + (N+1) * max_coeff >= rhs + slack
//
// Now if we assume cut * scaling + lift > b * scaling,
// by taking the largest N >=0 such that lift > N * scaling, we have
// lift <= (N + 1) * scaling, so cut * scaling > (b - (N+1)) * scaling
//
// And adding scaling * 2/ + 1/ we prove what we want:
// cut * scaling + lift > b * scaling => original_sum + lift_sum > rhs.
const IntegerValue slack = -rhs;
const IntegerValue remainder = max_coeff - slack;
const auto f = GetSuperAdditiveRoundingFunction(remainder, max_coeff,
IntegerValue(1), max_scaling);
const IntegerValue scaling = f(max_coeff);
if (scaling > 1) {
for (int i = 0; i < cover_size; ++i) {
base_ct_.terms[i].coeff *= scaling;
}
base_ct_.rhs *= scaling;
}
num_lifting_ = 0;
for (int i = cover_size; i < base_size; ++i) {
const IntegerValue positive_coeff = base_ct_.terms[i].coeff;
const IntegerValue new_coeff = f(positive_coeff);
base_ct_.terms[i].coeff = new_coeff;
if (new_coeff != 0) ++num_lifting_;
}
total_num_lifting_ += num_lifting_;
if (!cut_builder_.ConvertToLinearConstraint(base_ct_, &cut_)) {
cut_.Clear();
++total_num_overflow_abort_;
return false;
}
if (scaling > 1) DivideByGCD(&cut_);
return true;
}
bool CoverCutHelper::TryWithLetchfordSouliLifting(
const CutData& input, ImpliedBoundsProcessor* ib_processor) {
cut_.Clear();
base_ct_ = input;
if (!MakeAllTermsPositive()) return false;
// Perform IB expansion with no restriction, all coeff should still be
// positive.
//
// TODO(user): Merge Boolean terms that are complement of each other.
if (ib_processor != nullptr) {
cut_builder_.RegisterAllBooleansTerms(base_ct_);
const int old_size = static_cast<int>(base_ct_.terms.size());
for (int i = 0; i < old_size; ++i) {
if (base_ct_.terms[i].bound_diff <= 1) continue;
if (ib_processor->TryToExpandWithLowerImpliedbound(
IntegerValue(1), i,
/*complement=*/false, &base_ct_, &cut_builder_)) {
++total_num_ibs_;
}
}
}
// TODO(user): we currently only deal with Boolean in the cover. Fix.
const int num_bools =
std::partition(base_ct_.terms.begin(), base_ct_.terms.end(),
[](const CutTerm& t) { return t.bound_diff == 1; }) -
base_ct_.terms.begin();
if (num_bools == 0) return false;
IntegerValue rhs;
const int cover_size = GetCoverSize(num_bools, &rhs);
if (cover_size == 0) return false;
// Collect the weight in the cover.
IntegerValue sum(0);
std::vector<IntegerValue> cover_weights;
for (int i = 0; i < cover_size; ++i) {
CHECK_EQ(base_ct_.terms[i].bound_diff, 1);
CHECK_GT(base_ct_.terms[i].coeff, 0);
cover_weights.push_back(base_ct_.terms[i].coeff);
sum = CapAddI(sum, base_ct_.terms[i].coeff);
}
if (AtMinOrMaxInt64(sum.value())) {
++total_num_overflow_abort_;
return false;
}
CHECK_GT(sum, base_ct_.rhs);
// Compute the correct threshold so that if we round down larger weights to
// p/q. We have sum of the weight in cover == base_rhs.
IntegerValue p(0);
IntegerValue q(0);
IntegerValue previous_sum(0);
std::sort(cover_weights.begin(), cover_weights.end());
for (int i = 0; i < cover_size; ++i) {
q = IntegerValue(cover_weights.size() - i);
if (previous_sum + cover_weights[i] * q > base_ct_.rhs) {
p = base_ct_.rhs - previous_sum;
break;
}
previous_sum += cover_weights[i];
}
CHECK_GE(q, 1);
// Compute thresholds.
// For the first q values, thresholds[i] is the smallest integer such that
// q * threshold[i] > p * (i + 1).
std::vector<IntegerValue> thresholds;
for (int i = 0; i < q; ++i) {
// TODO(user): compute this in an overflow-safe way.
if (CapProd(p.value(), i + 1) >= std::numeric_limits<int64_t>::max() - 1) {
++total_num_overflow_abort_;
return false;
}
thresholds.push_back(CeilRatio(p * (i + 1) + 1, q));
}
// For the other values, we just add the weights.
std::reverse(cover_weights.begin(), cover_weights.end());
for (int i = q.value(); i < cover_size; ++i) {
thresholds.push_back(thresholds.back() + cover_weights[i]);
}
CHECK_EQ(thresholds.back(), base_ct_.rhs + 1);
// Generate the cut.
//
// Our algo is quadratic in worst case, but large coefficients should be
// rare, and in practice we don't really see this.
//
// Note that this work for non-Boolean since we can just "theorically" split
// them as a sum of Booleans :) Probably a cleaner proof exist by just using
// the super-additivity of the lifting function on [0, rhs].
temp_cut_.rhs = IntegerValue(cover_size - 1);
temp_cut_.terms.clear();
num_lifting_ = 0;
const int base_size = static_cast<int>(base_ct_.terms.size());
for (int i = 0; i < base_size; ++i) {
const CutTerm& term = base_ct_.terms[i];
const IntegerValue coeff = term.coeff;
IntegerValue cut_coeff(1);
if (coeff < thresholds[0]) {
if (i >= cover_size) continue;
} else {
// Find the largest index <= coeff.
//
// TODO(user): For exact multiple of p/q we can increase the coeff by 1/2.
// See section in the paper on getting maximal super additive function.
for (int i = 1; i < cover_size; ++i) {
if (coeff < thresholds[i]) break;
cut_coeff = IntegerValue(i + 1);
}
if (cut_coeff != 0 && i >= cover_size) ++num_lifting_;
if (cut_coeff > 1 && i < cover_size) ++num_lifting_; // happen?
}
temp_cut_.terms.push_back(term);
temp_cut_.terms.back().coeff = cut_coeff;
}
if (!cut_builder_.ConvertToLinearConstraint(temp_cut_, &cut_)) {
cut_.Clear();
++total_num_overflow_abort_;
return false;
}
return true;
}
CutGenerator CreatePositiveMultiplicationCutGenerator(AffineExpression z,
AffineExpression x,
AffineExpression y,
int linearization_level,
Model* model) {
CutGenerator result;
if (z.var != kNoIntegerVariable) result.vars.push_back(z.var);
if (x.var != kNoIntegerVariable) result.vars.push_back(x.var);
if (y.var != kNoIntegerVariable) result.vars.push_back(y.var);
IntegerTrail* const integer_trail = model->GetOrCreate<IntegerTrail>();
Trail* trail = model->GetOrCreate<Trail>();
result.generate_cuts =
[z, x, y, linearization_level, model, trail, integer_trail](
const absl::StrongVector<IntegerVariable, double>& lp_values,
LinearConstraintManager* manager) {
if (trail->CurrentDecisionLevel() > 0 && linearization_level == 1) {
return true;
}
const int64_t x_lb = integer_trail->LevelZeroLowerBound(x).value();
const int64_t x_ub = integer_trail->LevelZeroUpperBound(x).value();
const int64_t y_lb = integer_trail->LevelZeroLowerBound(y).value();
const int64_t y_ub = integer_trail->LevelZeroUpperBound(y).value();
// if x or y are fixed, the McCormick equations are exact.
if (x_lb == x_ub || y_lb == y_ub) return true;
// Check for overflow with the product of expression bounds and the
// product of one expression bound times the constant part of the other
// expression.
const int64_t x_max_amp = std::max(std::abs(x_lb), std::abs(x_ub));
const int64_t y_max_amp = std::max(std::abs(y_lb), std::abs(y_ub));
constexpr int64_t kMaxSafeInteger = (int64_t{1} << 53) - 1;
if (CapProd(y_max_amp, x_max_amp) > kMaxSafeInteger) return true;
if (CapProd(y_max_amp, std::abs(x.constant.value())) >
kMaxSafeInteger) {
return true;
}
if (CapProd(x_max_amp, std::abs(y.constant.value())) >
kMaxSafeInteger) {
return true;
}
const double x_lp_value = x.LpValue(lp_values);
const double y_lp_value = y.LpValue(lp_values);
const double z_lp_value = z.LpValue(lp_values);
// TODO(user): As the bounds change monotonically, these cuts
// dominate any previous one. try to keep a reference to the cut and
// replace it. Alternatively, add an API for a level-zero bound change
// callback.
// Cut -z + x_coeff * x + y_coeff* y <= rhs
auto try_add_above_cut = [&](int64_t x_coeff, int64_t y_coeff,
int64_t rhs) {
if (-z_lp_value + x_lp_value * x_coeff + y_lp_value * y_coeff >=
rhs + kMinCutViolation) {
LinearConstraintBuilder cut(model, /*lb=*/kMinIntegerValue,
/*ub=*/IntegerValue(rhs));
cut.AddTerm(z, IntegerValue(-1));
if (x_coeff != 0) cut.AddTerm(x, IntegerValue(x_coeff));
if (y_coeff != 0) cut.AddTerm(y, IntegerValue(y_coeff));
manager->AddCut(cut.Build(), "PositiveProduct", lp_values);
}
};
// Cut -z + x_coeff * x + y_coeff* y >= rhs
auto try_add_below_cut = [&](int64_t x_coeff, int64_t y_coeff,
int64_t rhs) {
if (-z_lp_value + x_lp_value * x_coeff + y_lp_value * y_coeff <=
rhs - kMinCutViolation) {
LinearConstraintBuilder cut(model, /*lb=*/IntegerValue(rhs),
/*ub=*/kMaxIntegerValue);
cut.AddTerm(z, IntegerValue(-1));
if (x_coeff != 0) cut.AddTerm(x, IntegerValue(x_coeff));
if (y_coeff != 0) cut.AddTerm(y, IntegerValue(y_coeff));
manager->AddCut(cut.Build(), "PositiveProduct", lp_values);
}
};
// McCormick relaxation of bilinear constraints. These 4 cuts are the
// exact facets of the x * y polyhedron for a bounded x and y.
//
// Each cut correspond to plane that contains two of the line
// (x=x_lb), (x=x_ub), (y=y_lb), (y=y_ub). The easiest to
// understand them is to draw the x*y curves and see the 4
// planes that correspond to the convex hull of the graph.
try_add_above_cut(y_lb, x_lb, x_lb * y_lb);
try_add_above_cut(y_ub, x_ub, x_ub * y_ub);
try_add_below_cut(y_ub, x_lb, x_lb * y_ub);
try_add_below_cut(y_lb, x_ub, x_ub * y_lb);
return true;
};
return result;
}
LinearConstraint ComputeHyperplanAboveSquare(AffineExpression x,
AffineExpression square,
IntegerValue x_lb,
IntegerValue x_ub, Model* model) {
const IntegerValue above_slope = x_ub + x_lb;
LinearConstraintBuilder above_hyperplan(model, kMinIntegerValue,
-x_lb * x_ub);
above_hyperplan.AddTerm(square, 1);
above_hyperplan.AddTerm(x, IntegerValue(-above_slope));
return above_hyperplan.Build();
}
LinearConstraint ComputeHyperplanBelowSquare(AffineExpression x,
AffineExpression square,
IntegerValue x_value,
Model* model) {
const IntegerValue below_slope = 2 * x_value + 1;
LinearConstraintBuilder below_hyperplan(model, -x_value - x_value * x_value,
kMaxIntegerValue);
below_hyperplan.AddTerm(square, 1);
below_hyperplan.AddTerm(x, -below_slope);
return below_hyperplan.Build();
}
CutGenerator CreateSquareCutGenerator(AffineExpression y, AffineExpression x,
int linearization_level, Model* model) {
CutGenerator result;
if (x.var != kNoIntegerVariable) result.vars.push_back(x.var);
if (y.var != kNoIntegerVariable) result.vars.push_back(y.var);
Trail* trail = model->GetOrCreate<Trail>();
IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
result.generate_cuts =
[y, x, linearization_level, trail, integer_trail, model](
const absl::StrongVector<IntegerVariable, double>& lp_values,
LinearConstraintManager* manager) {
if (trail->CurrentDecisionLevel() > 0 && linearization_level == 1) {
return true;
}
const IntegerValue x_ub = integer_trail->LevelZeroUpperBound(x);
const IntegerValue x_lb = integer_trail->LevelZeroLowerBound(x);
DCHECK_GE(x_lb, 0);
if (x_lb == x_ub) return true;
// Check for potential overflows.
if (x_ub > (int64_t{1} << 31)) return true;
DCHECK_GE(x_lb, 0);
manager->AddCut(ComputeHyperplanAboveSquare(x, y, x_lb, x_ub, model),
"SquareUpper", lp_values);
const IntegerValue x_floor =
static_cast<int64_t>(std::floor(x.LpValue(lp_values)));
manager->AddCut(ComputeHyperplanBelowSquare(x, y, x_floor, model),
"SquareLower", lp_values);
return true;
};
return result;
}
ImpliedBoundsProcessor::BestImpliedBoundInfo
ImpliedBoundsProcessor::GetCachedImpliedBoundInfo(IntegerVariable var) const {
auto it = cache_.find(var);
if (it != cache_.end()) {
BestImpliedBoundInfo result = it->second;
if (result.bool_var == kNoIntegerVariable) return BestImpliedBoundInfo();
if (integer_trail_->IsFixed(result.bool_var)) return BestImpliedBoundInfo();
return result;
}
return BestImpliedBoundInfo();
}
ImpliedBoundsProcessor::BestImpliedBoundInfo
ImpliedBoundsProcessor::ComputeBestImpliedBound(
IntegerVariable var,
const absl::StrongVector<IntegerVariable, double>& lp_values) {
auto it = cache_.find(var);
if (it != cache_.end()) return it->second;
BestImpliedBoundInfo result;
double result_slack_lp_value = std::numeric_limits<double>::infinity();
const IntegerValue lb = integer_trail_->LevelZeroLowerBound(var);
for (const ImpliedBoundEntry& entry :
implied_bounds_->GetImpliedBounds(var)) {
// Only process entries with a Boolean variable currently part of the LP
// we are considering for this cut.
//
// TODO(user): the more we use cuts, the less it make sense to have a
// lot of small independent LPs.
if (!lp_vars_.contains(PositiveVariable(entry.literal_view))) {
continue;
}
// The equation is X = lb + diff * Bool + Slack where Bool is in [0, 1]
// and slack in [0, ub - lb].
const IntegerValue diff = entry.lower_bound - lb;
CHECK_GE(diff, 0);
const double bool_lp_value = entry.is_positive
? lp_values[entry.literal_view]
: 1.0 - lp_values[entry.literal_view];
const double slack_lp_value =
lp_values[var] - ToDouble(lb) - bool_lp_value * ToDouble(diff);
// If the implied bound equation is not respected, we just add it
// to implied_bound_cuts, and skip the entry for now.
if (slack_lp_value < -1e-4) {
LinearConstraint ib_cut;
ib_cut.lb = kMinIntegerValue;
std::vector<std::pair<IntegerVariable, IntegerValue>> terms;
if (entry.is_positive) {
// X >= Indicator * (bound - lb) + lb
terms.push_back({entry.literal_view, diff});
terms.push_back({var, IntegerValue(-1)});
ib_cut.ub = -lb;
} else {
// X >= -Indicator * (bound - lb) + bound
terms.push_back({entry.literal_view, -diff});
terms.push_back({var, IntegerValue(-1)});
ib_cut.ub = -entry.lower_bound;
}
CleanTermsAndFillConstraint(&terms, &ib_cut);
ib_cut_pool_.AddCut(std::move(ib_cut), "IB", lp_values);
continue;
}
// We look for tight implied bounds, and amongst the tightest one, we
// prefer larger coefficient in front of the Boolean.
if (slack_lp_value + 1e-4 < result_slack_lp_value ||
(slack_lp_value < result_slack_lp_value + 1e-4 &&
entry.lower_bound > result.implied_bound)) {
result_slack_lp_value = slack_lp_value;
result.var_lp_value = lp_values[var];
result.bool_lp_value = bool_lp_value;
result.implied_bound = entry.lower_bound;
result.is_positive = entry.is_positive;
result.bool_var = entry.literal_view;
}
}
cache_[var] = result;
return result;
}
void ImpliedBoundsProcessor::RecomputeCacheAndSeparateSomeImpliedBoundCuts(
const absl::StrongVector<IntegerVariable, double>& lp_values) {
cache_.clear();
for (const IntegerVariable var :
implied_bounds_->VariablesWithImpliedBounds()) {
if (!lp_vars_.contains(PositiveVariable(var))) continue;
ComputeBestImpliedBound(var, lp_values);
}
}
// Important: The cut_builder_ must have been reset.
bool ImpliedBoundsProcessor::TryToExpandWithLowerImpliedbound(
IntegerValue factor_t, int i, bool complement, CutData* cut,
CutDataBuilder* builder) {
CutTerm& term = cut->terms[i];
// We only want to expand non-Boolean and non-slack term!
if (term.bound_diff <= 1) return false;
if (!term.IsSimple()) return false;
CHECK_EQ(IntTypeAbs(term.expr_coeffs[0]), 1);
// Try lower bounded direction for implied bound.
// This kind should always be beneficial if it exists:
//
// Because X = bound_diff * B + S
// We can replace coeff * X by the expression before applying f:
// = f(coeff * bound_diff) * B + f(coeff) * [X - bound_diff * B]
// = f(coeff) * X + (f(coeff * bound_diff) - f(coeff) * bound_diff] * B
// So we can "lift" B into the cut with a non-negative coefficient.
//
// Note that this lifting is really the same as if we used that implied
// bound before since in f(coeff * bound_diff) * B + f(coeff) * S, if we
// replace S by its value [X - bound_diff * B] we get the same result.
//
// TODO(user): Ignore if bound_diff == 1 ? But we can still merge B with
// another entry if it exists, so it can still be good in this case.
//
// TODO(user): Only do it if coeff_b > 0 ? But again we could still merge
// B with an existing Boolean for a better cut even if coeff_b == 0.
const IntegerVariable ib_var = term.expr_coeffs[0] > 0
? term.expr_vars[0]
: NegationOf(term.expr_vars[0]);
const ImpliedBoundsProcessor::BestImpliedBoundInfo info =
GetCachedImpliedBoundInfo(ib_var);
const IntegerValue lb = -term.expr_offset;
const IntegerValue bound_diff = info.implied_bound - lb;
if (bound_diff <= 0) return false;
if (info.bool_var == kNoIntegerVariable) return false;
if (ProdOverflow(factor_t, CapProdI(term.coeff, bound_diff))) return false;
// We have X = info.diff * Boolean + slack.
CutTerm bool_term;
bool_term.coeff = term.coeff * bound_diff;
bool_term.expr_vars[0] = info.bool_var;
bool_term.expr_coeffs[1] = 0;
bool_term.bound_diff = IntegerValue(1);
bool_term.lp_value = info.bool_lp_value;
if (info.is_positive) {
bool_term.expr_coeffs[0] = IntegerValue(1);
bool_term.expr_offset = IntegerValue(0);
} else {
bool_term.expr_coeffs[0] = IntegerValue(-1);
bool_term.expr_offset = IntegerValue(1);
}
// Create slack.
// The expression is term.exp - bound_diff * bool_term
// The variable shouldn't be the same.
DCHECK_NE(term.expr_vars[0], bool_term.expr_vars[0]);
CutTerm slack_term;
slack_term.expr_vars[0] = term.expr_vars[0];
slack_term.expr_coeffs[0] = term.expr_coeffs[0];
slack_term.expr_vars[1] = bool_term.expr_vars[0];
slack_term.expr_coeffs[1] = -bound_diff * bool_term.expr_coeffs[0];
slack_term.expr_offset =
term.expr_offset - bound_diff * bool_term.expr_offset;
slack_term.lp_value = info.SlackLpValue(lb);
slack_term.coeff = term.coeff;
slack_term.bound_diff = term.bound_diff;
// It should be good to use IB, but sometime we have things like
// 7.3 = 2 * bool@1 + 5.3 and the expanded Boolean is at its upper bound.
// It is always good to complement such variable.
//
// Note that here we do more and just complement anything closer to UB.
//
// TODO(user): Because of merges, we might have entry with a coefficient of
// zero than are not useful. Remove them.
if (complement) {
if (bool_term.lp_value > 0.5) {
bool_term.Complement(&cut->rhs);
}
if (slack_term.lp_value >
0.5 * static_cast<double>(slack_term.bound_diff.value())) {
slack_term.Complement(&cut->rhs);
}
}
term = slack_term;
builder->AddOrMergeTerm(bool_term, factor_t, cut);
return true;
}
std::string SingleNodeFlow::DebugString() const {
return absl::StrCat("#in:", in_flow.size(), " #out:", out_flow.size(),
" demand:", demand.value(), " #bool:", num_bool,
" #lb:", num_to_lb, " #ub:", num_to_ub);
}
bool FlowCoverCutHelper::TryXminusLB(IntegerVariable var, double lp_value,
IntegerValue lb, IntegerValue ub,
IntegerValue coeff,
ImpliedBoundsProcessor* ib_helper,
SingleNodeFlow* result) const {
const ImpliedBoundsProcessor::BestImpliedBoundInfo ib =
ib_helper->GetCachedImpliedBoundInfo(NegationOf(var));
if (ib.bool_var == kNoIntegerVariable) return false;
if (ib.implied_bound != -lb) return false;
// We have -var >= (ub - lb) bool - ub;
// so (var - lb) <= -(ub - lb) * bool + ub - lb;
// and (var - lb) <= bound_diff * (1 - bool).
FlowInfo info;
if (ib.is_positive) {
info.bool_lp_value = 1 - ib.bool_lp_value;
info.bool_expr.vars.push_back(ib.bool_var);
info.bool_expr.coeffs.push_back(-1);
info.bool_expr.offset = 1;
} else {
info.bool_lp_value = ib.bool_lp_value;
info.bool_expr.vars.push_back(ib.bool_var);
info.bool_expr.coeffs.push_back(1);
}
info.capacity = IntTypeAbs(coeff) * (ub - lb);
info.flow_lp_value = ToDouble(IntTypeAbs(coeff)) * (lp_value - ToDouble(lb));
info.flow_expr.vars.push_back(var);
info.flow_expr.coeffs.push_back(IntTypeAbs(coeff));
info.flow_expr.offset = -lb * IntTypeAbs(coeff);
// We use (var - lb) so sign is preserved
result->demand -= coeff * lb;
if (coeff > 0) {
result->in_flow.push_back(info);
} else {
result->out_flow.push_back(info);
}
return true;
}
bool FlowCoverCutHelper::TryUBminusX(IntegerVariable var, double lp_value,
IntegerValue lb, IntegerValue ub,
IntegerValue coeff,
ImpliedBoundsProcessor* ib_helper,
SingleNodeFlow* result) const {
const ImpliedBoundsProcessor::BestImpliedBoundInfo ib =
ib_helper->GetCachedImpliedBoundInfo(var);
if (ib.bool_var == kNoIntegerVariable) return false;
if (ib.implied_bound != ub) return false;
// We have var >= (ub - lb) bool + lb.
// so ub - var <= ub - (ub - lb) * bool - lb.
// and (ub - var) <= bound_diff * (1 - bool).
FlowInfo info;
if (ib.is_positive) {
info.bool_lp_value = 1 - ib.bool_lp_value;
info.bool_expr.vars.push_back(ib.bool_var);
info.bool_expr.coeffs.push_back(-1);
info.bool_expr.offset = 1;
} else {
info.bool_lp_value = ib.bool_lp_value;
info.bool_expr.vars.push_back(ib.bool_var);
info.bool_expr.coeffs.push_back(1);
}
info.capacity = IntTypeAbs(coeff) * (ub - lb);
info.flow_lp_value = ToDouble(IntTypeAbs(coeff)) * (ToDouble(ub) - lp_value);
info.flow_expr.vars.push_back(var);
info.flow_expr.coeffs.push_back(-IntTypeAbs(coeff));
info.flow_expr.offset = ub * IntTypeAbs(coeff);
// We reverse the sign because we use (ub - var) here.
// So coeff * var = -coeff * (ub - var) + coeff * ub;
result->demand -= coeff * ub;
if (coeff > 0) {
result->out_flow.push_back(info);
} else {
result->in_flow.push_back(info);
}
return true;
}
SingleNodeFlow FlowCoverCutHelper::ComputeFlowCoverRelaxation(
const LinearConstraint& base_ct,
const absl::StrongVector<IntegerVariable, double>& lp_values,
IntegerTrail* integer_trail, ImpliedBoundsProcessor* ib_helper) {
SingleNodeFlow result;
result.demand = base_ct.ub;
// Stats.
int num_bool = 0;
int num_to_ub = 0;
int num_to_lb = 0;
const int size = base_ct.vars.size();
for (int i = 0; i < size; ++i) {
// We can either use (X - LB) or (UB - X) for a variable in [0, capacity].
const IntegerVariable var = base_ct.vars[i];
const IntegerValue coeff = base_ct.coeffs[i];
// Hack: abort if coefficient in the base constraint are too large.
if (IntTypeAbs(coeff) > 1'000'000) {
result.clear();
return result;
}
const IntegerValue lb = integer_trail->LevelZeroLowerBound(var);
const IntegerValue ub = integer_trail->LevelZeroUpperBound(var);
const IntegerValue capacity(
CapProd(IntTypeAbs(coeff).value(), (ub - lb).value()));
if (capacity >= kMaxIntegerValue) {
// Abort.
result.clear();
return result;
}
if (lb == ub) {
// Fixed variable shouldn't really appear here.
result.demand -= coeff * lb;
continue;
}
// We have a Boolean, this is an easy case.
if (ub - lb == 1) {
++num_bool;
FlowInfo info;
info.bool_lp_value = (lp_values[var] - ToDouble(lb));
info.capacity = capacity;
info.bool_expr.vars.push_back(var);
info.bool_expr.coeffs.push_back(1);
info.bool_expr.offset = -lb;
info.flow_lp_value = ToDouble(capacity) * info.bool_lp_value;
info.flow_expr.vars.push_back(var);
info.flow_expr.coeffs.push_back(info.capacity);
info.flow_expr.offset = -lb * info.capacity;
// coeff * x = coeff * (x - lb) + coeff * lb;
result.demand -= coeff * lb;
if (coeff > 0) {
result.in_flow.push_back(info);
} else {
result.out_flow.push_back(info);
}
continue;
}
// TODO(user): Improve our logic to decide what implied bounds to use. We
// rely on the best implied bounds, not necessarily one implying var at its
// level zero bound like we need here.
const double lp = lp_values[var];
const bool prefer_lb = (lp - ToDouble(lb)) > (ToDouble(ub) - lp);
if (prefer_lb) {
if (TryXminusLB(var, lp, lb, ub, coeff, ib_helper, &result)) {
++num_to_lb;
continue;
}
if (TryUBminusX(var, lp, lb, ub, coeff, ib_helper, &result)) {
++num_to_ub;
continue;
}
} else {
if (TryUBminusX(var, lp, lb, ub, coeff, ib_helper, &result)) {
++num_to_ub;
continue;
}
if (TryXminusLB(var, lp, lb, ub, coeff, ib_helper, &result)) {
++num_to_lb;
continue;
}
}
// Abort.
// TODO(user): Technically we can always use a arc usage Boolean fixed to 1.
result.clear();
return result;
}
result.num_bool = num_bool;
result.num_to_ub = num_to_ub;
result.num_to_lb = num_to_lb;
return result;
}
// Reference: "Lifted flow cover inequalities for mixed 0-1 integer programs".
// Zonghao Gu, George L. Nemhauser, Martin W.P. Savelsbergh. 1999.
bool FlowCoverCutHelper::GenerateCut(const SingleNodeFlow& data) {
if (data.empty()) return false;
const double tolerance = 1e-2;
// We are looking for two subsets CI (in-flow subset) and CO (out-flow subset)
// so that sum_CI capa - sum_CO capa = demand + slack, slack > 0.
//
// Moreover we want to maximize sum_CI bool_lp_value + sum_CO bool_lp_value.
std::vector<bool> in_cover(data.in_flow.size(), false);
std::vector<bool> out_cover(data.out_flow.size(), false);
// Start by selecting all the possible in_flow (except low bool value) and
// all the out_flow with a bool value close to one.
IntegerValue slack;
{
IntegerValue sum_in(0);
IntegerValue sum_out(0);
for (int i = 0; i < data.in_flow.size(); ++i) {
const FlowInfo& info = data.in_flow[i];
if (info.bool_lp_value > tolerance) {
in_cover[i] = true;
sum_in += info.capacity;
}
}
for (int i = 0; i < data.out_flow.size(); ++i) {
const FlowInfo& info = data.out_flow[i];
if (info.bool_lp_value > 1 - tolerance) {
out_cover[i] = true;
sum_out += info.capacity;
}
}
// This is the best slack we can hope for.
slack = sum_in - sum_out - data.demand;
}
if (slack <= 0) return false;
// Now greedily remove item from the in_cover and add_item to the out_cover
// as long as we have remaining slack. We prefer item with a high score an
// low slack variation.
//
// Note that this is just the classic greedy heuristic of a knapsack problem.
if (slack > 1) {
struct Item {
bool correspond_to_in_flow;
int index;
double score;
};
std::vector<Item> actions;
for (int i = 0; i < data.in_flow.size(); ++i) {
if (!in_cover[i]) continue;
const FlowInfo& info = data.in_flow[i];
if (info.bool_lp_value > 1 - tolerance) continue; // Do not remove these.
actions.push_back(
{true, i, (1 - info.bool_lp_value) / ToDouble(info.capacity)});
}
for (int i = 0; i < data.out_flow.size(); ++i) {
if (out_cover[i]) continue;
const FlowInfo& info = data.out_flow[i];
if (info.bool_lp_value < tolerance) continue; // Do not add these.
actions.push_back(
{false, i, info.bool_lp_value / ToDouble(info.capacity)});
}
// Sort by decreasing score.
std::sort(actions.begin(), actions.end(),
[](const Item& a, const Item& b) { return a.score > b.score; });
// Greedily remove/add item as long as we have slack.
for (const Item& item : actions) {
if (item.correspond_to_in_flow) {
const IntegerValue delta = data.in_flow[item.index].capacity;
if (delta >= slack) continue;
slack -= delta;
in_cover[item.index] = false;
} else {
const IntegerValue delta = data.out_flow[item.index].capacity;
if (delta >= slack) continue;
slack -= delta;
out_cover[item.index] = true;
}
}
}
// The non-lifted simple generalized flow cover inequality (SGFCI) cut will be
// demand - sum_CI flow_i - sum_CI++ (capa_i - slack)(1 - bool_i)
// + sum_CO capa_i + sum_L- slack * bool_i + sum_L-- flow_i >=0
//
// Where CI++ are the arc with capa > slack in CI.
// And L is O \ CO. L- arc with capa > slack and L-- the other.
//
// TODO(user): Also try to generate the extended generalized flow cover
// inequality (EGFCI).
CHECK_GT(slack, 0);
// For display only.
slack_ = slack;
num_in_ignored_ = 0;
num_in_flow_ = 0;
num_in_bin_ = 0;
num_out_capa_ = 0;
num_out_flow_ = 0;
num_out_bin_ = 0;
cut_builder_.Clear();
for (int i = 0; i < data.in_flow.size(); ++i) {
const FlowInfo& info = data.in_flow[i];
if (!in_cover[i]) {
num_in_ignored_++;
continue;
}
num_in_flow_++;
cut_builder_.AddLinearExpression(info.flow_expr, -1);
if (info.capacity > slack) {
num_in_bin_++;
const IntegerValue coeff = info.capacity - slack;
cut_builder_.AddConstant(-coeff);
cut_builder_.AddLinearExpression(info.bool_expr, coeff);
}
}
for (int i = 0; i < data.out_flow.size(); ++i) {
const FlowInfo& info = data.out_flow[i];
if (out_cover[i]) {
num_out_capa_++;
cut_builder_.AddConstant(info.capacity);
} else if (info.capacity > slack) {
num_out_bin_++;
cut_builder_.AddLinearExpression(info.bool_expr, slack);
} else {
num_out_flow_++;
cut_builder_.AddLinearExpression(info.flow_expr);
}
}
// TODO(user): Lift the cut.
cut_ = cut_builder_.BuildConstraint(-data.demand, kMaxIntegerValue);
return true;
}
namespace {
int64_t SumOfKMinValues(const absl::btree_set<int64_t>& values, int k) {
int count = 0;
int64_t sum = 0;
for (const int64_t value : values) {
sum += value;
if (++count >= k) return sum;
}
return sum;
}
// Copies the data before sorting.
int64_t SumOfDiffMins(std::vector<int64_t> values) {
std::sort(values.begin(), values.end());
int64_t tmp_value = std::numeric_limits<int64_t>::min();
int64_t result = 0;
for (const int64_t value : values) {
// Make sure values are different.
tmp_value = std::max(tmp_value + 1, value);
result += tmp_value;
}
return result;
}
void TryToGenerateAllDiffCut(
const std::vector<std::pair<double, AffineExpression>>& sorted_exprs_lp,
const IntegerTrail& integer_trail,
const absl::StrongVector<IntegerVariable, double>& lp_values,
LinearConstraintManager* manager, Model* model) {
const int num_exprs = sorted_exprs_lp.size();
std::vector<AffineExpression> current_set_exprs;
absl::btree_set<int64_t> min_values;
absl::btree_set<int64_t> negated_max_values;
std::vector<int64_t> expr_mins;
std::vector<int64_t> negated_expr_maxes;
double sum = 0.0;
TopNCuts top_n_cuts(5);
for (const auto& [expr_lp, expr] : sorted_exprs_lp) {
sum += expr_lp;
expr_mins.push_back(integer_trail.LevelZeroLowerBound(expr).value());
negated_expr_maxes.push_back(
-integer_trail.LevelZeroUpperBound(expr).value());
if (integer_trail.IsFixed(expr)) {
const int64_t value = integer_trail.FixedValue(expr).value();
min_values.insert(value);
negated_max_values.insert(-value);
} else {
int count = 0;
const int64_t coeff = expr.coeff.value();
const int64_t constant = expr.constant.value();
for (const int64_t value :
integer_trail.InitialVariableDomain(expr.var).Values()) {
if (coeff > 0) {
min_values.insert(value * coeff + constant);
} else {
negated_max_values.insert(-(value * coeff + constant));
}
if (++count >= num_exprs) break;
}
count = 0;
for (const int64_t value :
integer_trail.InitialVariableDomain(expr.var).Negation().Values()) {
if (coeff > 0) {
negated_max_values.insert(value * coeff - constant);
} else {
min_values.insert(-value * coeff + constant);
}
if (++count >= num_exprs) break;
}
}
current_set_exprs.push_back(expr);
CHECK_EQ(current_set_exprs.size(), expr_mins.size());
CHECK_EQ(current_set_exprs.size(), negated_expr_maxes.size());
const int64_t required_min_sum_smallest =
SumOfKMinValues(min_values, current_set_exprs.size());
const int64_t required_max_sum_smallest =
-SumOfKMinValues(negated_max_values, current_set_exprs.size());
const int64_t required_min_sum_alldiff = SumOfDiffMins(expr_mins);
const int64_t required_max_sum_alldiff = -SumOfDiffMins(negated_expr_maxes);
const int64_t required_min_sum =
std::max(required_min_sum_smallest, required_min_sum_alldiff);
const int64_t required_max_sum =
std::min(required_max_sum_smallest, required_max_sum_alldiff);
const std::string min_suffix =
required_min_sum_alldiff > required_min_sum_smallest ? "a"
: (required_min_sum_alldiff == required_min_sum_smallest) ? "e"
: "d";
const std::string max_suffix =
required_max_sum_alldiff < required_max_sum_smallest ? "a"
: (required_max_sum_alldiff == required_max_sum_smallest) ? "e"
: "d";
if (sum < static_cast<double>(required_min_sum) - kMinCutViolation ||
sum > static_cast<double>(required_max_sum) + kMinCutViolation) {
LinearConstraintBuilder cut(model, required_min_sum, required_max_sum);
for (AffineExpression expr : current_set_exprs) {
cut.AddTerm(expr, IntegerValue(1));
}
top_n_cuts.AddCut(cut.Build(),
absl::StrCat("AllDiff_", min_suffix, max_suffix),
lp_values);
// NOTE: We can extend the current set but it is more helpful to generate
// the cut on a different set of variables so we reset the counters.
sum = 0.0;
current_set_exprs.clear();
min_values.clear();
negated_max_values.clear();
expr_mins.clear();
negated_expr_maxes.clear();
}
}
top_n_cuts.TransferToManager(lp_values, manager);
}
} // namespace
CutGenerator CreateAllDifferentCutGenerator(
const std::vector<AffineExpression>& exprs, Model* model) {
CutGenerator result;
IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
for (const AffineExpression& expr : exprs) {
if (!integer_trail->IsFixed(expr)) {
result.vars.push_back(expr.var);
}
}
gtl::STLSortAndRemoveDuplicates(&result.vars);
Trail* trail = model->GetOrCreate<Trail>();
result.generate_cuts =
[exprs, integer_trail, trail, model](
const absl::StrongVector<IntegerVariable, double>& lp_values,
LinearConstraintManager* manager) {
// These cuts work at all levels but the generator adds too many cuts on
// some instances and degrade the performance so we only use it at level
// 0.
if (trail->CurrentDecisionLevel() > 0) return true;
std::vector<std::pair<double, AffineExpression>> sorted_exprs;
for (const AffineExpression expr : exprs) {
if (integer_trail->LevelZeroLowerBound(expr) ==
integer_trail->LevelZeroUpperBound(expr)) {
continue;
}
sorted_exprs.push_back(std::make_pair(expr.LpValue(lp_values), expr));
}
std::sort(sorted_exprs.begin(), sorted_exprs.end(),
[](std::pair<double, AffineExpression>& a,
const std::pair<double, AffineExpression>& b) {
return a.first < b.first;
});
TryToGenerateAllDiffCut(sorted_exprs, *integer_trail, lp_values,
manager, model);
// Other direction.
std::reverse(sorted_exprs.begin(), sorted_exprs.end());
TryToGenerateAllDiffCut(sorted_exprs, *integer_trail, lp_values,
manager, model);
return true;
};
VLOG(2) << "Created all_diff cut generator of size: " << exprs.size();
return result;
}
namespace {
// Returns max((w2i - w1i)*Li, (w2i - w1i)*Ui).
IntegerValue MaxCornerDifference(const IntegerVariable var,
const IntegerValue w1_i,
const IntegerValue w2_i,
const IntegerTrail& integer_trail) {
const IntegerValue lb = integer_trail.LevelZeroLowerBound(var);
const IntegerValue ub = integer_trail.LevelZeroUpperBound(var);
return std::max((w2_i - w1_i) * lb, (w2_i - w1_i) * ub);
}
// This is the coefficient of zk in the cut, where k = max_index.
// MPlusCoefficient_ki = max((wki - wI(i)i) * Li,
// (wki - wI(i)i) * Ui)
// = max corner difference for variable i,
// target expr I(i), max expr k.
// The coefficient of zk is Sum(i=1..n)(MPlusCoefficient_ki) + bk
IntegerValue MPlusCoefficient(
const std::vector<IntegerVariable>& x_vars,
const std::vector<LinearExpression>& exprs,
const absl::StrongVector<IntegerVariable, int>& variable_partition,
const int max_index, const IntegerTrail& integer_trail) {
IntegerValue coeff = exprs[max_index].offset;
// TODO(user): This algo is quadratic since GetCoefficientOfPositiveVar()
// is linear. This can be optimized (better complexity) if needed.
for (const IntegerVariable var : x_vars) {
const int target_index = variable_partition[var];
if (max_index != target_index) {
coeff += MaxCornerDifference(
var, GetCoefficientOfPositiveVar(var, exprs[target_index]),
GetCoefficientOfPositiveVar(var, exprs[max_index]), integer_trail);
}
}
return coeff;
}
// Compute the value of
// rhs = wI(i)i * xi + Sum(k=1..d)(MPlusCoefficient_ki * zk)
// for variable xi for given target index I(i).
double ComputeContribution(
const IntegerVariable xi_var, const std::vector<IntegerVariable>& z_vars,
const std::vector<LinearExpression>& exprs,
const absl::StrongVector<IntegerVariable, double>& lp_values,
const IntegerTrail& integer_trail, const int target_index) {
CHECK_GE(target_index, 0);
CHECK_LT(target_index, exprs.size());
const LinearExpression& target_expr = exprs[target_index];
const double xi_value = lp_values[xi_var];
const IntegerValue wt_i = GetCoefficientOfPositiveVar(xi_var, target_expr);
double contrib = ToDouble(wt_i) * xi_value;
for (int expr_index = 0; expr_index < exprs.size(); ++expr_index) {
if (expr_index == target_index) continue;
const LinearExpression& max_expr = exprs[expr_index];
const double z_max_value = lp_values[z_vars[expr_index]];
const IntegerValue corner_value = MaxCornerDifference(
xi_var, wt_i, GetCoefficientOfPositiveVar(xi_var, max_expr),
integer_trail);
contrib += ToDouble(corner_value) * z_max_value;
}
return contrib;
}
} // namespace
CutGenerator CreateLinMaxCutGenerator(
const IntegerVariable target, const std::vector<LinearExpression>& exprs,
const std::vector<IntegerVariable>& z_vars, Model* model) {
CutGenerator result;
std::vector<IntegerVariable> x_vars;
result.vars = {target};
const int num_exprs = exprs.size();
for (int i = 0; i < num_exprs; ++i) {
result.vars.push_back(z_vars[i]);
x_vars.insert(x_vars.end(), exprs[i].vars.begin(), exprs[i].vars.end());
}
gtl::STLSortAndRemoveDuplicates(&x_vars);
// All expressions should only contain positive variables.
DCHECK(std::all_of(x_vars.begin(), x_vars.end(), [](IntegerVariable var) {
return VariableIsPositive(var);
}));
result.vars.insert(result.vars.end(), x_vars.begin(), x_vars.end());
IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
result.generate_cuts =
[x_vars, z_vars, target, num_exprs, exprs, integer_trail, model](
const absl::StrongVector<IntegerVariable, double>& lp_values,
LinearConstraintManager* manager) {
absl::StrongVector<IntegerVariable, int> variable_partition(
lp_values.size(), -1);
absl::StrongVector<IntegerVariable, double> variable_partition_contrib(
lp_values.size(), std::numeric_limits<double>::infinity());
for (int expr_index = 0; expr_index < num_exprs; ++expr_index) {
for (const IntegerVariable var : x_vars) {
const double contribution = ComputeContribution(
var, z_vars, exprs, lp_values, *integer_trail, expr_index);
const double prev_contribution = variable_partition_contrib[var];
if (contribution < prev_contribution) {
variable_partition[var] = expr_index;
variable_partition_contrib[var] = contribution;
}
}
}
LinearConstraintBuilder cut(model, /*lb=*/IntegerValue(0),
/*ub=*/kMaxIntegerValue);
double violation = lp_values[target];
cut.AddTerm(target, IntegerValue(-1));
for (const IntegerVariable xi_var : x_vars) {
const int input_index = variable_partition[xi_var];
const LinearExpression& expr = exprs[input_index];
const IntegerValue coeff = GetCoefficientOfPositiveVar(xi_var, expr);
if (coeff != IntegerValue(0)) {
cut.AddTerm(xi_var, coeff);
}
violation -= ToDouble(coeff) * lp_values[xi_var];
}
for (int expr_index = 0; expr_index < num_exprs; ++expr_index) {
const IntegerVariable z_var = z_vars[expr_index];
const IntegerValue z_coeff = MPlusCoefficient(
x_vars, exprs, variable_partition, expr_index, *integer_trail);
if (z_coeff != IntegerValue(0)) {
cut.AddTerm(z_var, z_coeff);
}
violation -= ToDouble(z_coeff) * lp_values[z_var];
}
if (violation > 1e-2) {
manager->AddCut(cut.Build(), "LinMax", lp_values);
}
return true;
};
return result;
}
namespace {
IntegerValue EvaluateMaxAffine(
const std::vector<std::pair<IntegerValue, IntegerValue>>& affines,
IntegerValue x) {
IntegerValue y = kMinIntegerValue;
for (const auto& p : affines) {
y = std::max(y, x * p.first + p.second);
}
return y;
}
} // namespace
bool BuildMaxAffineUpConstraint(
const LinearExpression& target, IntegerVariable var,
const std::vector<std::pair<IntegerValue, IntegerValue>>& affines,
Model* model, LinearConstraintBuilder* builder) {
auto* integer_trail = model->GetOrCreate<IntegerTrail>();
const IntegerValue x_min = integer_trail->LevelZeroLowerBound(var);
const IntegerValue x_max = integer_trail->LevelZeroUpperBound(var);
const IntegerValue y_at_min = EvaluateMaxAffine(affines, x_min);
const IntegerValue y_at_max = EvaluateMaxAffine(affines, x_max);
const IntegerValue delta_x = x_max - x_min;
const IntegerValue delta_y = y_at_max - y_at_min;
// target <= y_at_min + (delta_y / delta_x) * (var - x_min)
// delta_x * target <= delta_x * y_at_min + delta_y * (var - x_min)
// -delta_y * var + delta_x * target <= delta_x * y_at_min - delta_y * x_min
//
// Checks the rhs for overflows.
if (AtMinOrMaxInt64(CapProd(delta_x.value(), y_at_min.value())) ||
AtMinOrMaxInt64(CapProd(delta_y.value(), x_min.value()))) {
return false;
}
builder->ResetBounds(kMinIntegerValue, delta_x * y_at_min - delta_y * x_min);
builder->AddLinearExpression(target, delta_x);
builder->AddTerm(var, -delta_y);
// Prevent to create constraints that can overflow.
if (!ValidateLinearConstraintForOverflow(builder->Build(), *integer_trail)) {
VLOG(2) << "Linear constraint can cause overflow: " << builder->Build();
return false;
}
return true;
}
CutGenerator CreateMaxAffineCutGenerator(
LinearExpression target, IntegerVariable var,
std::vector<std::pair<IntegerValue, IntegerValue>> affines,
const std::string cut_name, Model* model) {
CutGenerator result;
result.vars = target.vars;
result.vars.push_back(var);
gtl::STLSortAndRemoveDuplicates(&result.vars);
IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
result.generate_cuts =
[target, var, affines, cut_name, integer_trail, model](
const absl::StrongVector<IntegerVariable, double>& lp_values,
LinearConstraintManager* manager) {
if (integer_trail->IsFixed(var)) return true;
LinearConstraintBuilder builder(model);
if (BuildMaxAffineUpConstraint(target, var, affines, model, &builder)) {
manager->AddCut(builder.Build(), cut_name, lp_values);
}
return true;
};
return result;
}
CutGenerator CreateCliqueCutGenerator(
const std::vector<IntegerVariable>& base_variables, Model* model) {
// Filter base_variables to only keep the one with a literal view, and
// do the conversion.
std::vector<IntegerVariable> variables;
std::vector<Literal> literals;
absl::flat_hash_map<LiteralIndex, IntegerVariable> positive_map;
absl::flat_hash_map<LiteralIndex, IntegerVariable> negative_map;
auto* integer_trail = model->GetOrCreate<IntegerTrail>();
auto* encoder = model->GetOrCreate<IntegerEncoder>();
for (const IntegerVariable var : base_variables) {
if (integer_trail->LowerBound(var) != IntegerValue(0)) continue;
if (integer_trail->UpperBound(var) != IntegerValue(1)) continue;
const LiteralIndex literal_index = encoder->GetAssociatedLiteral(
IntegerLiteral::GreaterOrEqual(var, IntegerValue(1)));
if (literal_index != kNoLiteralIndex) {
variables.push_back(var);
literals.push_back(Literal(literal_index));
positive_map[literal_index] = var;
negative_map[Literal(literal_index).NegatedIndex()] = var;
}
}
CutGenerator result;
result.vars = variables;
auto* implication_graph = model->GetOrCreate<BinaryImplicationGraph>();
result.generate_cuts =
[variables, literals, implication_graph, positive_map, negative_map,
model](const absl::StrongVector<IntegerVariable, double>& lp_values,
LinearConstraintManager* manager) {
std::vector<double> packed_values;
for (int i = 0; i < literals.size(); ++i) {
packed_values.push_back(lp_values[variables[i]]);
}
const std::vector<std::vector<Literal>> at_most_ones =
implication_graph->GenerateAtMostOnesWithLargeWeight(literals,
packed_values);
for (const std::vector<Literal>& at_most_one : at_most_ones) {
// We need to express such "at most one" in term of the initial
// variables, so we do not use the
// LinearConstraintBuilder::AddLiteralTerm() here.
LinearConstraintBuilder builder(
model, IntegerValue(std::numeric_limits<int64_t>::min()),
IntegerValue(1));
for (const Literal l : at_most_one) {
if (positive_map.contains(l.Index())) {
builder.AddTerm(positive_map.at(l.Index()), IntegerValue(1));
} else {
// Add 1 - X to the linear constraint.
builder.AddTerm(negative_map.at(l.Index()), IntegerValue(-1));
builder.AddConstant(IntegerValue(1));
}
}
manager->AddCut(builder.Build(), "Clique", lp_values);
}
return true;
};
return result;
}
} // namespace sat
} // namespace operations_research
Reference in New Issue View Git Blame Copy Permalink