Google OR-Tools v9.12
a fast and portable software suite for combinatorial optimization
Loading...
Searching...
No Matches
linear_constraint.cc
Go to the documentation of this file.
1// Copyright 2010-2025 Google LLC
2// Licensed under the Apache License, Version 2.0 (the "License");
3// you may not use this file except in compliance with the License.
4// You may obtain a copy of the License at
5//
6// http://www.apache.org/licenses/LICENSE-2.0
7//
8// Unless required by applicable law or agreed to in writing, software
9// distributed under the License is distributed on an "AS IS" BASIS,
10// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11// See the License for the specific language governing permissions and
12// limitations under the License.
13
15
16#include <algorithm>
17#include <cmath>
18#include <cstddef>
19#include <cstdint>
20#include <limits>
21#include <string>
22#include <utility>
23#include <vector>
24
25#include "absl/base/attributes.h"
26#include "absl/container/flat_hash_set.h"
27#include "absl/log/check.h"
28#include "absl/strings/str_cat.h"
29#include "absl/types/span.h"
32#include "ortools/sat/integer.h"
37
38namespace operations_research {
39namespace sat {
40
41void LinearConstraintBuilder::AddTerm(IntegerVariable var, IntegerValue coeff) {
42 if (coeff == 0) return;
43 // We can either add var or NegationOf(var), and we always choose the
44 // positive one.
45 if (VariableIsPositive(var)) {
46 terms_.push_back({var, coeff});
47 } else {
48 terms_.push_back({NegationOf(var), -coeff});
49 }
50}
51
53 IntegerValue coeff) {
54 if (coeff == 0) return;
55 // We can either add var or NegationOf(var), and we always choose the
56 // positive one.
57 if (expr.var != kNoIntegerVariable) {
58 if (VariableIsPositive(expr.var)) {
59 terms_.push_back({expr.var, coeff * expr.coeff});
60 } else {
61 terms_.push_back({NegationOf(expr.var), -coeff * expr.coeff});
62 }
63 }
64 DCHECK(!ProdOverflow(coeff, expr.constant));
65 offset_ += coeff * expr.constant;
66}
67
69 const LinearExpression& expr) {
70 AddLinearExpression(expr, IntegerValue(1));
71}
72
74 IntegerValue coeff) {
75 for (int i = 0; i < expr.vars.size(); ++i) {
76 // We must use positive variables.
77 if (VariableIsPositive(expr.vars[i])) {
78 terms_.push_back({expr.vars[i], expr.coeffs[i] * coeff});
79 } else {
80 terms_.push_back({NegationOf(expr.vars[i]), -expr.coeffs[i] * coeff});
81 }
82 }
83 offset_ += expr.offset * coeff;
84}
85
87 absl::Span<const LiteralValueValue> product) {
88 if (product.empty()) return true;
89
90 IntegerValue product_min = kMaxIntegerValue;
91 // TODO(user): Checks the value of literals.
92 for (const LiteralValueValue& term : product) {
93 product_min = std::min(product_min, term.left_value * term.right_value);
94 }
95
96 for (const LiteralValueValue& term : product) {
97 IntegerValue coeff = term.left_value * term.right_value - product_min;
98 if (coeff == 0) continue;
99 if (!AddLiteralTerm(term.literal, coeff)) {
100 return false;
101 }
102 }
103 AddConstant(product_min);
104 return true;
105}
106
108 AffineExpression left, AffineExpression right, IntegerTrail* integer_trail,
109 bool* is_quadratic) {
110 if (integer_trail->IsFixed(left)) {
111 AddTerm(right, integer_trail->FixedValue(left));
112 } else if (integer_trail->IsFixed(right)) {
113 AddTerm(left, integer_trail->FixedValue(right));
114 } else {
115 const IntegerValue left_min = integer_trail->LowerBound(left);
116 const IntegerValue right_min = integer_trail->LowerBound(right);
117 AddTerm(left, right_min);
118 AddTerm(right, left_min);
119 // Substract the energy counted twice.
120 AddConstant(-left_min * right_min);
121 if (is_quadratic != nullptr) *is_quadratic = true;
122 }
123}
124
125void LinearConstraintBuilder::AddConstant(IntegerValue value) {
126 offset_ += value;
127}
128
130 Literal lit, IntegerValue coeff) {
131 DCHECK(encoder_ != nullptr);
132 IntegerVariable var = kNoIntegerVariable;
133 bool view_is_direct = true;
134 if (!encoder_->LiteralOrNegationHasView(lit, &var, &view_is_direct)) {
135 return false;
136 }
137
138 if (view_is_direct) {
139 AddTerm(var, coeff);
140 } else {
141 AddTerm(var, -coeff);
142 offset_ += coeff;
143 }
144 return true;
145}
146
150
152 IntegerValue ub) {
153 LinearConstraint result;
154 result.lb = lb > kMinIntegerValue ? lb - offset_ : lb;
155 result.ub = ub < kMaxIntegerValue ? ub - offset_ : ub;
156 CleanTermsAndFillConstraint(&terms_, &result);
157 return result;
158}
159
161 IntegerValue lb, IntegerValue ub, LinearConstraint* ct) {
162 ct->lb = lb > kMinIntegerValue ? lb - offset_ : lb;
163 ct->ub = ub < kMaxIntegerValue ? ub - offset_ : ub;
165}
166
168 LinearExpression result;
169 CleanTermsAndFillConstraint(&terms_, &result);
170 result.offset = offset_;
171 return result;
172}
173
175 const LinearConstraint& constraint,
177 int i = 0;
178 const int size = constraint.num_terms;
179 const int shifted_size = size - 3;
180 double a0 = 0.0;
181 double a1 = 0.0;
182 double a2 = 0.0;
183 double a3 = 0.0;
184 const double* view = values.data();
185 for (; i < shifted_size; i += 4) {
186 a0 += static_cast<double>(constraint.coeffs[i].value()) *
187 view[constraint.vars[i].value()];
188 a1 += static_cast<double>(constraint.coeffs[i + 1].value()) *
189 view[constraint.vars[i + 1].value()];
190 a2 += static_cast<double>(constraint.coeffs[i + 2].value()) *
191 view[constraint.vars[i + 2].value()];
192 a3 += static_cast<double>(constraint.coeffs[i + 3].value()) *
193 view[constraint.vars[i + 3].value()];
194 }
195 double activity = a0 + a1 + a2 + a3;
196 if (i < size) {
197 activity += static_cast<double>(constraint.coeffs[i].value()) *
198 view[constraint.vars[i].value()];
199 if (i + 1 < size) {
200 activity += static_cast<double>(constraint.coeffs[i + 1].value()) *
201 view[constraint.vars[i + 1].value()];
202 if (i + 2 < size) {
203 activity += static_cast<double>(constraint.coeffs[i + 2].value()) *
204 view[constraint.vars[i + 2].value()];
205 }
206 }
207 }
208 return activity;
209}
210
212 double sum = 0.0;
213 for (int i = 0; i < ct.num_terms; ++i) {
214 sum += ToDouble(ct.coeffs[i]) * ToDouble(ct.coeffs[i]);
215 }
216 return std::sqrt(sum);
217}
218
219IntegerValue ComputeInfinityNorm(const LinearConstraint& ct) {
220 IntegerValue result(0);
221 for (int i = 0; i < ct.num_terms; ++i) {
222 result = std::max(result, IntTypeAbs(ct.coeffs[i]));
223 }
224 return result;
225}
226
227double ScalarProduct(const LinearConstraint& ct1, const LinearConstraint& ct2) {
228 if (ct1.num_terms == 0 || ct2.num_terms == 0) return 0.0;
229 DCHECK(std::is_sorted(ct1.vars.get(), ct1.vars.get() + ct1.num_terms));
230 DCHECK(std::is_sorted(ct2.vars.get(), ct2.vars.get() + ct2.num_terms));
231 double scalar_product = 0.0;
232 int index_1 = 0;
233 int index_2 = 0;
234 IntegerVariable var1 = ct1.vars[index_1];
235 IntegerVariable var2 = ct2.vars[index_2];
236 while (true) {
237 if (var1 == var2) {
238 scalar_product += static_cast<double>(ct1.coeffs[index_1].value()) *
239 static_cast<double>(ct2.coeffs[index_2].value());
240 if (++index_1 == ct1.num_terms) break;
241 if (++index_2 == ct2.num_terms) break;
242 var1 = ct1.vars[index_1];
243 var2 = ct2.vars[index_2];
244 } else if (var1 > var2) {
245 if (++index_2 == ct2.num_terms) break;
246 var2 = ct2.vars[index_2];
247 } else {
248 if (++index_1 == ct1.num_terms) break;
249 var1 = ct1.vars[index_1];
250 }
251 }
252 return scalar_product;
253}
254
255namespace {
256
257// TODO(user): Template for any integer type and expose this?
258IntegerValue ComputeGcd(absl::Span<const IntegerValue> values) {
259 if (values.empty()) return IntegerValue(1);
260 int64_t gcd = 0;
261 for (const IntegerValue value : values) {
262 gcd = MathUtil::GCD64(gcd, std::abs(value.value()));
263 if (gcd == 1) break;
264 }
265 if (gcd < 0) return IntegerValue(1); // Can happen with kint64min.
266 return IntegerValue(gcd);
267}
268
269} // namespace
270
271void DivideByGCD(LinearConstraint* constraint) {
272 if (constraint->num_terms == 0) return;
273 const IntegerValue gcd = ComputeGcd(
274 {constraint->coeffs.get(), static_cast<size_t>(constraint->num_terms)});
275 if (gcd == 1) return;
276
277 if (constraint->lb > kMinIntegerValue) {
278 constraint->lb = CeilRatio(constraint->lb, gcd);
279 }
280 if (constraint->ub < kMaxIntegerValue) {
281 constraint->ub = FloorRatio(constraint->ub, gcd);
282 }
283 for (int i = 0; i < constraint->num_terms; ++i) {
284 constraint->coeffs[i] /= gcd;
285 }
286}
287
289 int new_size = 0;
290 const int size = constraint->num_terms;
291 for (int i = 0; i < size; ++i) {
292 if (constraint->coeffs[i] == 0) continue;
293 constraint->vars[new_size] = constraint->vars[i];
294 constraint->coeffs[new_size] = constraint->coeffs[i];
295 ++new_size;
296 }
297 constraint->resize(new_size);
298}
299
301 const int size = constraint->num_terms;
302 for (int i = 0; i < size; ++i) {
303 const IntegerValue coeff = constraint->coeffs[i];
304 if (coeff < 0) {
305 constraint->coeffs[i] = -coeff;
306 constraint->vars[i] = NegationOf(constraint->vars[i]);
307 }
308 }
309}
310
312 const int size = constraint->num_terms;
313 for (int i = 0; i < size; ++i) {
314 const IntegerVariable var = constraint->vars[i];
315 if (!VariableIsPositive(var)) {
316 constraint->coeffs[i] = -constraint->coeffs[i];
317 constraint->vars[i] = NegationOf(var);
318 }
319 }
320}
321
324 double result = ToDouble(offset);
325 for (int i = 0; i < vars.size(); ++i) {
326 result += ToDouble(coeffs[i]) * lp_values[vars[i]];
327 }
328 return result;
329}
330
331IntegerValue LinearExpression::LevelZeroMin(IntegerTrail* integer_trail) const {
332 IntegerValue result = offset;
333 for (int i = 0; i < vars.size(); ++i) {
334 DCHECK_GE(coeffs[i], 0);
335 result += coeffs[i] * integer_trail->LevelZeroLowerBound(vars[i]);
336 }
337 return result;
338}
339
340IntegerValue LinearExpression::Min(const IntegerTrail& integer_trail) const {
341 IntegerValue result = offset;
342 for (int i = 0; i < vars.size(); ++i) {
343 if (coeffs[i] > 0) {
344 result += coeffs[i] * integer_trail.LowerBound(vars[i]);
345 } else {
346 result += coeffs[i] * integer_trail.UpperBound(vars[i]);
347 }
348 }
349 return result;
350}
351
352IntegerValue LinearExpression::Max(const IntegerTrail& integer_trail) const {
353 IntegerValue result = offset;
354 for (int i = 0; i < vars.size(); ++i) {
355 if (coeffs[i] > 0) {
356 result += coeffs[i] * integer_trail.UpperBound(vars[i]);
357 } else {
358 result += coeffs[i] * integer_trail.LowerBound(vars[i]);
359 }
360 }
361 return result;
362}
363
364std::string LinearExpression::DebugString() const {
365 if (vars.empty()) return absl::StrCat(offset.value());
366 std::string result;
367 for (int i = 0; i < vars.size(); ++i) {
368 absl::StrAppend(&result, i > 0 ? " " : "",
370 }
371 if (offset != 0) {
372 absl::StrAppend(&result, " + ", offset.value());
373 }
374 return result;
375}
376
378 absl::flat_hash_set<IntegerVariable> seen_variables;
379 const int size = ct.num_terms;
380 for (int i = 0; i < size; ++i) {
381 if (VariableIsPositive(ct.vars[i])) {
382 if (!seen_variables.insert(ct.vars[i]).second) return false;
383 } else {
384 if (!seen_variables.insert(NegationOf(ct.vars[i])).second) return false;
385 }
386 }
387 return true;
388}
389
391 LinearExpression canonical_expr;
392 canonical_expr.offset = expr.offset;
393 for (int i = 0; i < expr.vars.size(); ++i) {
394 if (expr.coeffs[i] < 0) {
395 canonical_expr.vars.push_back(NegationOf(expr.vars[i]));
396 canonical_expr.coeffs.push_back(-expr.coeffs[i]);
397 } else {
398 canonical_expr.vars.push_back(expr.vars[i]);
399 canonical_expr.coeffs.push_back(expr.coeffs[i]);
400 }
401 }
402 return canonical_expr;
403}
404
405// TODO(user): Avoid duplication with PossibleIntegerOverflow() in the checker?
406// At least make sure the code is the same.
408 const IntegerTrail& integer_trail) {
409 int64_t positive_sum(0);
410 int64_t negative_sum(0);
411 for (int i = 0; i < constraint.num_terms; ++i) {
412 const IntegerVariable var = constraint.vars[i];
413 const IntegerValue coeff = constraint.coeffs[i];
414 const IntegerValue lb = integer_trail.LevelZeroLowerBound(var);
415 const IntegerValue ub = integer_trail.LevelZeroUpperBound(var);
416
417 int64_t min_prod = CapProd(coeff.value(), lb.value());
418 int64_t max_prod = CapProd(coeff.value(), ub.value());
419 if (min_prod > max_prod) std::swap(min_prod, max_prod);
420
421 positive_sum = CapAdd(positive_sum, std::max(int64_t{0}, max_prod));
422 negative_sum = CapAdd(negative_sum, std::min(int64_t{0}, min_prod));
423 }
424
425 const int64_t limit = std::numeric_limits<int64_t>::max();
426 if (positive_sum >= limit) return false;
427 if (negative_sum <= -limit) return false;
428 if (CapSub(positive_sum, negative_sum) >= limit) return false;
429
430 return true;
431}
432
434 LinearExpression result;
435 result.vars = NegationOf(expr.vars);
436 result.coeffs = expr.coeffs;
437 result.offset = -expr.offset;
438 return result;
439}
440
442 LinearExpression result;
443 result.offset = expr.offset;
444 for (int i = 0; i < expr.vars.size(); ++i) {
445 if (VariableIsPositive(expr.vars[i])) {
446 result.vars.push_back(expr.vars[i]);
447 result.coeffs.push_back(expr.coeffs[i]);
448 } else {
449 result.vars.push_back(NegationOf(expr.vars[i]));
450 result.coeffs.push_back(-expr.coeffs[i]);
451 }
452 }
453 return result;
454}
455
456IntegerValue GetCoefficient(const IntegerVariable var,
457 const LinearExpression& expr) {
458 for (int i = 0; i < expr.vars.size(); ++i) {
459 if (expr.vars[i] == var) {
460 return expr.coeffs[i];
461 } else if (expr.vars[i] == NegationOf(var)) {
462 return -expr.coeffs[i];
463 }
464 }
465 return IntegerValue(0);
466}
467
468IntegerValue GetCoefficientOfPositiveVar(const IntegerVariable var,
469 const LinearExpression& expr) {
470 CHECK(VariableIsPositive(var));
471 for (int i = 0; i < expr.vars.size(); ++i) {
472 if (expr.vars[i] == var) {
473 return expr.coeffs[i];
474 }
475 }
476 return IntegerValue(0);
477}
478
479bool PossibleOverflow(const IntegerTrail& integer_trail,
480 const LinearConstraint& constraint) {
481 IntegerValue min_activity(0);
482 IntegerValue max_activity(0);
483 const int size = constraint.num_terms;
484 for (int i = 0; i < size; ++i) {
485 const IntegerVariable var = constraint.vars[i];
486 const IntegerValue coeff = constraint.coeffs[i];
487 CHECK_NE(coeff, 0);
488 const IntegerValue lb = integer_trail.LevelZeroLowerBound(var);
489 const IntegerValue ub = integer_trail.LevelZeroUpperBound(var);
490 if (coeff > 0) {
491 if (!AddProductTo(lb, coeff, &min_activity)) return true;
492 if (!AddProductTo(ub, coeff, &max_activity)) return true;
493 } else {
494 if (!AddProductTo(ub, coeff, &min_activity)) return true;
495 if (!AddProductTo(lb, coeff, &max_activity)) return true;
496 }
497 }
498 return AtMinOrMaxInt64(CapSub(max_activity.value(), min_activity.value()));
499}
500
501} // namespace sat
502} // namespace operations_research
static int64_t GCD64(int64_t x, int64_t y)
Definition mathutil.h:108
IntegerValue LowerBound(IntegerVariable i) const
Returns the current lower/upper bound of the given integer variable.
Definition integer.h:1301
IntegerValue UpperBound(IntegerVariable i) const
Definition integer.h:1305
IntegerValue FixedValue(IntegerVariable i) const
Checks that the variable is fixed and returns its value.
Definition integer.h:1313
bool IsFixed(IntegerVariable i) const
Checks if the variable is fixed.
Definition integer.h:1309
IntegerValue LevelZeroUpperBound(IntegerVariable var) const
Definition integer.h:1401
IntegerValue LevelZeroLowerBound(IntegerVariable var) const
Returns globally valid lower/upper bound on the given integer variable.
Definition integer.h:1396
void AddQuadraticLowerBound(AffineExpression left, AffineExpression right, IntegerTrail *integer_trail, bool *is_quadratic=nullptr)
ABSL_MUST_USE_RESULT bool AddLiteralTerm(Literal lit, IntegerValue coeff=IntegerValue(1))
void AddTerm(IntegerVariable var, IntegerValue coeff)
void AddLinearExpression(const LinearExpression &expr)
void AddConstant(IntegerValue value)
Adds the corresponding term to the current linear expression.
LinearConstraint BuildConstraint(IntegerValue lb, IntegerValue ub)
bool BuildIntoConstraintAndCheckOverflow(IntegerValue lb, IntegerValue ub, LinearConstraint *ct)
ABSL_MUST_USE_RESULT bool AddDecomposedProduct(absl::Span< const LiteralValueValue > product)
value_type * data()
– Pass-through methods to STL vector ----------------------------------—
double ComputeActivity(const LinearConstraint &constraint, const util_intops::StrongVector< IntegerVariable, double > &values)
IntegerValue FloorRatio(IntegerValue dividend, IntegerValue positive_divisor)
bool AddProductTo(IntegerValue a, IntegerValue b, IntegerValue *result)
Computes result += a * b, and return false iff there is an overflow.
void DivideByGCD(LinearConstraint *constraint)
constexpr IntegerValue kMaxIntegerValue(std::numeric_limits< IntegerValue::ValueType >::max() - 1)
IntType IntTypeAbs(IntType t)
double ComputeL2Norm(const LinearConstraint &ct)
Returns sqrt(sum square(coeff)).
IntegerValue CeilRatio(IntegerValue dividend, IntegerValue positive_divisor)
bool ProdOverflow(IntegerValue t, IntegerValue value)
void RemoveZeroTerms(LinearConstraint *constraint)
Removes the entries with a coefficient of zero.
IntegerValue ComputeInfinityNorm(const LinearConstraint &ct)
Returns the maximum absolute value of the coefficients.
std::string IntegerTermDebugString(IntegerVariable var, IntegerValue coeff)
std::vector< IntegerVariable > NegationOf(absl::Span< const IntegerVariable > vars)
Returns the vector of the negated variables.
Definition integer.cc:52
constexpr IntegerValue kMinIntegerValue(-kMaxIntegerValue.value())
const IntegerVariable kNoIntegerVariable(-1)
bool PossibleOverflow(const IntegerTrail &integer_trail, const LinearConstraint &constraint)
void CleanTermsAndFillConstraint(std::vector< std::pair< IntegerVariable, IntegerValue > > *terms, LinearExpression *output)
void MakeAllCoefficientsPositive(LinearConstraint *constraint)
Makes all coefficients positive by transforming a variable to its negation.
LinearExpression CanonicalizeExpr(const LinearExpression &expr)
bool MergePositiveVariableTermsAndCheckForOverflow(std::vector< std::pair< IntegerVariable, IntegerValue > > *terms, LinearConstraint *output)
LinearExpression PositiveVarExpr(const LinearExpression &expr)
Returns the same expression with positive variables.
IntegerValue GetCoefficient(const IntegerVariable var, const LinearExpression &expr)
void MakeAllVariablesPositive(LinearConstraint *constraint)
Makes all variables "positive" by transforming a variable to its negation.
IntegerValue GetCoefficientOfPositiveVar(const IntegerVariable var, const LinearExpression &expr)
bool ValidateLinearConstraintForOverflow(const LinearConstraint &constraint, const IntegerTrail &integer_trail)
bool NoDuplicateVariable(const LinearConstraint &ct)
Returns false if duplicate variables are found in ct.
bool VariableIsPositive(IntegerVariable i)
double ToDouble(IntegerValue value)
double ScalarProduct(const LinearConstraint &ct1, const LinearConstraint &ct2)
In SWIG mode, we don't want anything besides these top-level includes.
bool AtMinOrMaxInt64(int64_t x)
Checks if x is equal to the min or the max value of an int64_t.
int64_t CapAdd(int64_t x, int64_t y)
int64_t CapSub(int64_t x, int64_t y)
int64_t CapProd(int64_t x, int64_t y)
std::unique_ptr< IntegerValue[]> coeffs
std::unique_ptr< IntegerVariable[]> vars
IntegerValue LevelZeroMin(IntegerTrail *integer_trail) const
IntegerValue Min(const IntegerTrail &integer_trail) const
double LpValue(const util_intops::StrongVector< IntegerVariable, double > &lp_values) const
Return[s] the evaluation of the linear expression.
IntegerValue Max(const IntegerTrail &integer_trail) const