Google OR-Tools v9.12
a fast and portable software suite for combinatorial optimization
Loading...
Searching...
No Matches
constraint_violation.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 <cstdint>
18#include <cstdlib>
19#include <limits>
20#include <memory>
21#include <optional>
22#include <utility>
23#include <vector>
24
25#include "absl/algorithm/container.h"
26#include "absl/container/flat_hash_map.h"
27#include "absl/container/flat_hash_set.h"
28#include "absl/log/check.h"
29#include "absl/types/span.h"
34#include "ortools/sat/cp_model.pb.h"
36#include "ortools/sat/util.h"
40
41namespace operations_research {
42namespace sat {
43
44namespace {
45
46int64_t ExprValue(const LinearExpressionProto& expr,
47 absl::Span<const int64_t> solution) {
48 int64_t result = expr.offset();
49 for (int i = 0; i < expr.vars_size(); ++i) {
50 result += solution[expr.vars(i)] * expr.coeffs(i);
51 }
52 return result;
53}
54
55LinearExpressionProto ExprDiff(const LinearExpressionProto& a,
56 const LinearExpressionProto& b) {
57 LinearExpressionProto result;
58 result.set_offset(a.offset() - b.offset());
59 result.mutable_vars()->Reserve(a.vars().size() + b.vars().size());
60 result.mutable_coeffs()->Reserve(a.vars().size() + b.vars().size());
61 for (int i = 0; i < a.vars().size(); ++i) {
62 result.add_vars(a.vars(i));
63 result.add_coeffs(a.coeffs(i));
64 }
65 for (int i = 0; i < b.vars().size(); ++i) {
66 result.add_vars(b.vars(i));
67 result.add_coeffs(-b.coeffs(i));
68 }
69 return result;
70}
71
72LinearExpressionProto LinearExprSum(LinearExpressionProto a,
73 LinearExpressionProto b) {
74 LinearExpressionProto result;
75 result.set_offset(a.offset() + b.offset());
76 result.mutable_vars()->Reserve(a.vars().size() + b.vars().size());
77 result.mutable_coeffs()->Reserve(a.vars().size() + b.vars().size());
78 for (const LinearExpressionProto& p : {a, b}) {
79 for (int i = 0; i < p.vars().size(); ++i) {
80 result.add_vars(p.vars(i));
81 result.add_coeffs(p.coeffs(i));
82 }
83 }
84 return result;
85}
86
87LinearExpressionProto NegatedLinearExpression(LinearExpressionProto a) {
88 LinearExpressionProto result = a;
89 result.set_offset(-a.offset());
90 for (int64_t& coeff : *result.mutable_coeffs()) {
91 coeff = -coeff;
92 }
93 return result;
94}
95
96int64_t ExprMin(const LinearExpressionProto& expr, const CpModelProto& model) {
97 int64_t result = expr.offset();
98 for (int i = 0; i < expr.vars_size(); ++i) {
99 const IntegerVariableProto& var_proto = model.variables(expr.vars(i));
100 if (expr.coeffs(i) > 0) {
101 result += expr.coeffs(i) * var_proto.domain(0);
102 } else {
103 result += expr.coeffs(i) * var_proto.domain(var_proto.domain_size() - 1);
104 }
105 }
106 return result;
107}
108
109int64_t ExprMax(const LinearExpressionProto& expr, const CpModelProto& model) {
110 int64_t result = expr.offset();
111 for (int i = 0; i < expr.vars_size(); ++i) {
112 const IntegerVariableProto& var_proto = model.variables(expr.vars(i));
113 if (expr.coeffs(i) > 0) {
114 result += expr.coeffs(i) * var_proto.domain(var_proto.domain_size() - 1);
115 } else {
116 result += expr.coeffs(i) * var_proto.domain(0);
117 }
118 }
119 return result;
120}
121
122bool LiteralValue(int lit, absl::Span<const int64_t> solution) {
123 if (RefIsPositive(lit)) {
124 return solution[lit] != 0;
125 } else {
126 return solution[PositiveRef(lit)] == 0;
127 }
128}
129
130} // namespace
131
132// ---- LinearIncrementalEvaluator -----
133
135 DCHECK(creation_phase_);
136 domains_.push_back(domain);
137 offsets_.push_back(0);
138 activities_.push_back(0);
139 num_false_enforcement_.push_back(0);
140 distances_.push_back(0);
141 is_violated_.push_back(false);
142 return num_constraints_++;
143}
144
146 DCHECK(creation_phase_);
147 const int var = PositiveRef(lit);
148 if (literal_entries_.size() <= var) {
149 literal_entries_.resize(var + 1);
150 }
151 literal_entries_[var].push_back(
152 {.ct_index = ct_index, .positive = RefIsPositive(lit)});
153}
154
155void LinearIncrementalEvaluator::AddLiteral(int ct_index, int lit,
156 int64_t coeff) {
157 DCHECK(creation_phase_);
158 if (RefIsPositive(lit)) {
159 AddTerm(ct_index, lit, coeff, 0);
160 } else {
161 AddTerm(ct_index, PositiveRef(lit), -coeff, coeff);
162 }
163}
164
165void LinearIncrementalEvaluator::AddTerm(int ct_index, int var, int64_t coeff,
166 int64_t offset) {
167 DCHECK(creation_phase_);
168 DCHECK_GE(var, 0);
169 if (coeff == 0) return;
170
171 if (var_entries_.size() <= var) {
172 var_entries_.resize(var + 1);
173 }
174 if (!var_entries_[var].empty() &&
175 var_entries_[var].back().ct_index == ct_index) {
176 var_entries_[var].back().coefficient += coeff;
177 if (var_entries_[var].back().coefficient == 0) {
178 var_entries_[var].pop_back();
179 }
180 } else {
181 var_entries_[var].push_back({.ct_index = ct_index, .coefficient = coeff});
182 }
183 AddOffset(ct_index, offset);
184 DCHECK(VarIsConsistent(var));
185}
186
187void LinearIncrementalEvaluator::AddOffset(int ct_index, int64_t offset) {
188 DCHECK(creation_phase_);
189 offsets_[ct_index] += offset;
190}
191
193 int ct_index, const LinearExpressionProto& expr, int64_t multiplier) {
194 DCHECK(creation_phase_);
195 AddOffset(ct_index, expr.offset() * multiplier);
196 for (int i = 0; i < expr.vars_size(); ++i) {
197 if (expr.coeffs(i) * multiplier == 0) continue;
198 AddTerm(ct_index, expr.vars(i), expr.coeffs(i) * multiplier);
199 }
200}
201
203 if (var_entries_.size() <= var) return true;
204
205 absl::flat_hash_set<int> visited;
206 for (const Entry& entry : var_entries_[var]) {
207 if (!visited.insert(entry.ct_index).second) return false;
208 }
209 return true;
210}
211
213 absl::Span<const int64_t> solution) {
214 DCHECK(!creation_phase_);
215
216 // Resets the activity as the offset and the number of false enforcement to 0.
217 activities_ = offsets_;
218 in_last_affected_variables_.resize(columns_.size(), false);
219 num_false_enforcement_.assign(num_constraints_, 0);
220
221 // Update these numbers for all columns.
222 const int num_vars = columns_.size();
223 for (int var = 0; var < num_vars; ++var) {
224 const SpanData& data = columns_[var];
225 const int64_t value = solution[var];
226
227 if (value == 0 && data.num_pos_literal > 0) {
228 const int* ct_indices = &ct_buffer_[data.start];
229 for (int k = 0; k < data.num_pos_literal; ++k) {
230 num_false_enforcement_[ct_indices[k]]++;
231 }
232 }
233
234 if (value == 1 && data.num_neg_literal > 0) {
235 const int* ct_indices = &ct_buffer_[data.start + data.num_pos_literal];
236 for (int k = 0; k < data.num_neg_literal; ++k) {
237 num_false_enforcement_[ct_indices[k]]++;
238 }
239 }
240
241 if (value != 0 && data.num_linear_entries > 0) {
242 const int* ct_indices =
243 &ct_buffer_[data.start + data.num_pos_literal + data.num_neg_literal];
244 const int64_t* coeffs = &coeff_buffer_[data.linear_start];
245 for (int k = 0; k < data.num_linear_entries; ++k) {
246 activities_[ct_indices[k]] += coeffs[k] * value;
247 }
248 }
249 }
250
251 // Cache violations (not counting enforcement).
252 for (int c = 0; c < num_constraints_; ++c) {
253 distances_[c] = domains_[c].Distance(activities_[c]);
254 is_violated_[c] = Violation(c) > 0;
255 }
256}
257
259 if (10 * last_affected_variables_.size() < columns_.size()) {
260 // Sparse.
261 in_last_affected_variables_.resize(columns_.size(), false);
262 for (const int var : last_affected_variables_) {
263 in_last_affected_variables_[var] = false;
264 }
265 } else {
266 // Dense.
267 in_last_affected_variables_.assign(columns_.size(), false);
268 }
269 last_affected_variables_.clear();
270 DCHECK(std::all_of(in_last_affected_variables_.begin(),
271 in_last_affected_variables_.end(),
272 [](bool b) { return !b; }));
273}
274
275// Tricky: Here we re-use in_last_affected_variables_ to resest
276// var_to_score_change. And in particular we need to list all variable whose
277// score changed here. Not just the one for which we have a decrease.
279 int c, absl::Span<const int64_t> jump_deltas,
280 absl::Span<double> var_to_score_change) {
281 if (c >= rows_.size()) return;
282
283 DCHECK_EQ(num_false_enforcement_[c], 0);
284 const SpanData& data = rows_[c];
285
286 // Update enforcement part. Because we only update weight of currently
287 // infeasible constraint, all change are 0 -> 1 transition and change by the
288 // same amount, which is the current distance.
289 const double enforcement_change = static_cast<double>(-distances_[c]);
290 if (enforcement_change != 0.0) {
291 int i = data.start;
292 const int end = data.num_pos_literal + data.num_neg_literal;
293 num_ops_ += end;
294 for (int k = 0; k < end; ++k, ++i) {
295 const int var = row_var_buffer_[i];
296 if (!in_last_affected_variables_[var]) {
297 var_to_score_change[var] = enforcement_change;
298 in_last_affected_variables_[var] = true;
299 last_affected_variables_.push_back(var);
300 } else {
301 var_to_score_change[var] += enforcement_change;
302 }
303 }
304 }
305
306 // Update linear part.
307 if (data.num_linear_entries > 0) {
308 const int* row_vars = &row_var_buffer_[data.start + data.num_pos_literal +
309 data.num_neg_literal];
310 const int64_t* row_coeffs = &row_coeff_buffer_[data.linear_start];
311 num_ops_ += 2 * data.num_linear_entries;
312
313 // Computing general Domain distance is slow.
314 // TODO(user): optimize even more for one sided constraints.
315 // Note(user): I tried to factor the two usage of this, but it is slower.
316 const Domain& rhs = domains_[c];
317 const int64_t rhs_min = rhs.Min();
318 const int64_t rhs_max = rhs.Max();
319 const bool is_simple = rhs.NumIntervals() == 2;
320 const auto violation = [&rhs, rhs_min, rhs_max, is_simple](int64_t v) {
321 if (v >= rhs_max) {
322 return v - rhs_max;
323 } else if (v <= rhs_min) {
324 return rhs_min - v;
325 } else {
326 return is_simple ? int64_t{0} : rhs.Distance(v);
327 }
328 };
329
330 const int64_t old_distance = distances_[c];
331 const int64_t activity = activities_[c];
332 for (int k = 0; k < data.num_linear_entries; ++k) {
333 const int var = row_vars[k];
334 const int64_t coeff = row_coeffs[k];
335 const int64_t diff =
336 violation(activity + coeff * jump_deltas[var]) - old_distance;
337 if (!in_last_affected_variables_[var]) {
338 var_to_score_change[var] = static_cast<double>(diff);
339 in_last_affected_variables_[var] = true;
340 last_affected_variables_.push_back(var);
341 } else {
342 var_to_score_change[var] += static_cast<double>(diff);
343 }
344 }
345 }
346}
347
348void LinearIncrementalEvaluator::UpdateScoreOnNewlyEnforced(
349 int c, double weight, absl::Span<const int64_t> jump_deltas,
350 absl::Span<double> jump_scores) {
351 const SpanData& data = rows_[c];
352
353 // Everyone else had a zero cost transition that now become enforced ->
354 // unenforced. So they all have better score.
355 const double weight_time_violation =
356 weight * static_cast<double>(distances_[c]);
357 if (weight_time_violation > 0.0) {
358 int i = data.start;
359 const int end = data.num_pos_literal + data.num_neg_literal;
360 num_ops_ += end;
361 for (int k = 0; k < end; ++k, ++i) {
362 const int var = row_var_buffer_[i];
363 jump_scores[var] -= weight_time_violation;
364 if (!in_last_affected_variables_[var]) {
365 in_last_affected_variables_[var] = true;
366 last_affected_variables_.push_back(var);
367 }
368 }
369 }
370
371 // Update linear part! It was zero and is now a diff.
372 {
373 int i = data.start + data.num_pos_literal + data.num_neg_literal;
374 int j = data.linear_start;
375 num_ops_ += 2 * data.num_linear_entries;
376 const int64_t old_distance = distances_[c];
377 for (int k = 0; k < data.num_linear_entries; ++k, ++i, ++j) {
378 const int var = row_var_buffer_[i];
379 const int64_t coeff = row_coeff_buffer_[j];
380 const int64_t new_distance =
381 domains_[c].Distance(activities_[c] + coeff * jump_deltas[var]);
382 jump_scores[var] +=
383 weight * static_cast<double>(new_distance - old_distance);
384 if (!in_last_affected_variables_[var]) {
385 in_last_affected_variables_[var] = true;
386 last_affected_variables_.push_back(var);
387 }
388 }
389 }
390}
391
392void LinearIncrementalEvaluator::UpdateScoreOnNewlyUnenforced(
393 int c, double weight, absl::Span<const int64_t> jump_deltas,
394 absl::Span<double> jump_scores) {
395 const SpanData& data = rows_[c];
396
397 // Everyone else had a enforced -> unenforced transition that now become zero.
398 // So they all have worst score, and we don't need to update
399 // last_affected_variables_.
400 const double weight_time_violation =
401 weight * static_cast<double>(distances_[c]);
402 if (weight_time_violation > 0.0) {
403 int i = data.start;
404 const int end = data.num_pos_literal + data.num_neg_literal;
405 num_ops_ += end;
406 for (int k = 0; k < end; ++k, ++i) {
407 const int var = row_var_buffer_[i];
408 jump_scores[var] += weight_time_violation;
409 }
410 }
411
412 // Update linear part! It had a diff and is now zero.
413 {
414 int i = data.start + data.num_pos_literal + data.num_neg_literal;
415 int j = data.linear_start;
416 num_ops_ += 2 * data.num_linear_entries;
417 const int64_t old_distance = distances_[c];
418 for (int k = 0; k < data.num_linear_entries; ++k, ++i, ++j) {
419 const int var = row_var_buffer_[i];
420 const int64_t coeff = row_coeff_buffer_[j];
421 const int64_t new_distance =
422 domains_[c].Distance(activities_[c] + coeff * jump_deltas[var]);
423 jump_scores[var] -=
424 weight * static_cast<double>(new_distance - old_distance);
425 if (!in_last_affected_variables_[var]) {
426 in_last_affected_variables_[var] = true;
427 last_affected_variables_.push_back(var);
428 }
429 }
430 }
431}
432
433// We just need to modify the old/new transition that decrease the number of
434// enforcement literal at false.
435void LinearIncrementalEvaluator::UpdateScoreOfEnforcementIncrease(
436 int c, double score_change, absl::Span<const int64_t> jump_deltas,
437 absl::Span<double> jump_scores) {
438 if (score_change == 0.0) return;
439
440 const SpanData& data = rows_[c];
441 int i = data.start;
442 num_ops_ += data.num_pos_literal;
443 for (int k = 0; k < data.num_pos_literal; ++k, ++i) {
444 const int var = row_var_buffer_[i];
445 if (jump_deltas[var] == 1) {
446 jump_scores[var] += score_change;
447 if (score_change < 0.0 && !in_last_affected_variables_[var]) {
448 in_last_affected_variables_[var] = true;
449 last_affected_variables_.push_back(var);
450 }
451 }
452 }
453 num_ops_ += data.num_neg_literal;
454 for (int k = 0; k < data.num_neg_literal; ++k, ++i) {
455 const int var = row_var_buffer_[i];
456 if (jump_deltas[var] == -1) {
457 jump_scores[var] += score_change;
458 if (score_change < 0.0 && !in_last_affected_variables_[var]) {
459 in_last_affected_variables_[var] = true;
460 last_affected_variables_.push_back(var);
461 }
462 }
463 }
464}
465
466void LinearIncrementalEvaluator::UpdateScoreOnActivityChange(
467 int c, double weight, int64_t activity_delta,
468 absl::Span<const int64_t> jump_deltas, absl::Span<double> jump_scores) {
469 if (activity_delta == 0) return;
470 const SpanData& data = rows_[c];
471
472 // In some cases, we can know that the score of all the involved variable
473 // will not change. This is the case if whatever 1 variable change the
474 // violation delta before/after is the same.
475 //
476 // TODO(user): Maintain more precise bounds.
477 // - We could easily compute on each ComputeInitialActivities() the
478 // maximum increase/decrease per variable, and take the max as each
479 // variable changes?
480 // - Know if a constraint is only <= or >= !
481 const int64_t old_activity = activities_[c];
482 const int64_t new_activity = old_activity + activity_delta;
483 int64_t min_range;
484 int64_t max_range;
485 if (new_activity > old_activity) {
486 min_range = old_activity - row_max_variations_[c];
487 max_range = new_activity + row_max_variations_[c];
488 } else {
489 min_range = new_activity - row_max_variations_[c];
490 max_range = old_activity + row_max_variations_[c];
491 }
492
493 // If the violation delta was zero and will still always be zero, we can skip.
494 if (Domain(min_range, max_range).IsIncludedIn(domains_[c])) return;
495
496 // Enforcement is always enforced -> un-enforced.
497 // So it was -weight_time_distance and is now -weight_time_new_distance.
498 const double delta =
499 -weight *
500 static_cast<double>(domains_[c].Distance(new_activity) - distances_[c]);
501 if (delta != 0.0) {
502 int i = data.start;
503 const int end = data.num_pos_literal + data.num_neg_literal;
504 num_ops_ += end;
505 for (int k = 0; k < end; ++k, ++i) {
506 const int var = row_var_buffer_[i];
507 jump_scores[var] += delta;
508 if (delta < 0.0 && !in_last_affected_variables_[var]) {
509 in_last_affected_variables_[var] = true;
510 last_affected_variables_.push_back(var);
511 }
512 }
513 }
514
515 // If we are infeasible and no move can correct it, both old_b - old_a and
516 // new_b - new_a will have the same value. We only needed to update the
517 // violation of the enforced literal.
518 if (min_range >= domains_[c].Max() || max_range <= domains_[c].Min()) return;
519
520 // Update linear part.
521 if (data.num_linear_entries > 0) {
522 const int* row_vars = &row_var_buffer_[data.start + data.num_pos_literal +
523 data.num_neg_literal];
524 const int64_t* row_coeffs = &row_coeff_buffer_[data.linear_start];
525 num_ops_ += 2 * data.num_linear_entries;
526
527 // Computing general Domain distance is slow.
528 // TODO(user): optimize even more for one sided constraints.
529 // Note(user): I tried to factor the two usage of this, but it is slower.
530 const Domain& rhs = domains_[c];
531 const int64_t rhs_min = rhs.Min();
532 const int64_t rhs_max = rhs.Max();
533 const bool is_simple = rhs.NumIntervals() == 2;
534 const auto violation = [&rhs, rhs_min, rhs_max, is_simple](int64_t v) {
535 if (v >= rhs_max) {
536 return v - rhs_max;
537 } else if (v <= rhs_min) {
538 return rhs_min - v;
539 } else {
540 return is_simple ? int64_t{0} : rhs.Distance(v);
541 }
542 };
543
544 const int64_t old_a_minus_new_a =
545 distances_[c] - domains_[c].Distance(new_activity);
546 for (int k = 0; k < data.num_linear_entries; ++k) {
547 const int var = row_vars[k];
548 const int64_t impact = row_coeffs[k] * jump_deltas[var];
549 const int64_t old_b = violation(old_activity + impact);
550 const int64_t new_b = violation(new_activity + impact);
551
552 // The old score was:
553 // weight * static_cast<double>(old_b - old_a);
554 // the new score is
555 // weight * static_cast<double>(new_b - new_a); so the diff is:
556 // weight * static_cast<double>(new_b - new_a - old_b + old_a)
557 const int64_t diff = old_a_minus_new_a + new_b - old_b;
558
559 // TODO(user): If a variable is at its lower (resp. upper) bound, then
560 // we know that the score will always move in the same direction, so we
561 // might skip the last_affected_variables_ update.
562 jump_scores[var] += weight * static_cast<double>(diff);
563 if (!in_last_affected_variables_[var]) {
564 in_last_affected_variables_[var] = true;
565 last_affected_variables_.push_back(var);
566 }
567 }
568 }
569}
570
571// Note that the code assumes that a column has no duplicates ct indices.
573 int var, int64_t delta, absl::Span<const double> weights,
574 absl::Span<const int64_t> jump_deltas, absl::Span<double> jump_scores,
575 std::vector<int>* constraints_with_changed_violation) {
576 DCHECK(!creation_phase_);
577 DCHECK_NE(delta, 0);
578 if (var >= columns_.size()) return;
579
580 const SpanData& data = columns_[var];
581 int i = data.start;
582 num_ops_ += data.num_pos_literal;
583 for (int k = 0; k < data.num_pos_literal; ++k, ++i) {
584 const int c = ct_buffer_[i];
585 const int64_t v0 = Violation(c);
586 if (delta == 1) {
587 num_false_enforcement_[c]--;
588 DCHECK_GE(num_false_enforcement_[c], 0);
589 if (num_false_enforcement_[c] == 0) {
590 UpdateScoreOnNewlyEnforced(c, weights[c], jump_deltas, jump_scores);
591 } else if (num_false_enforcement_[c] == 1) {
592 const double enforcement_change =
593 weights[c] * static_cast<double>(distances_[c]);
594 UpdateScoreOfEnforcementIncrease(c, enforcement_change, jump_deltas,
595 jump_scores);
596 }
597 } else {
598 num_false_enforcement_[c]++;
599 if (num_false_enforcement_[c] == 1) {
600 UpdateScoreOnNewlyUnenforced(c, weights[c], jump_deltas, jump_scores);
601 } else if (num_false_enforcement_[c] == 2) {
602 const double enforcement_change =
603 weights[c] * static_cast<double>(distances_[c]);
604 UpdateScoreOfEnforcementIncrease(c, -enforcement_change, jump_deltas,
605 jump_scores);
606 }
607 }
608 const int64_t v1 = Violation(c);
609 is_violated_[c] = v1 > 0;
610 if (v1 != v0) {
611 constraints_with_changed_violation->push_back(c);
612 }
613 }
614 num_ops_ += data.num_neg_literal;
615 for (int k = 0; k < data.num_neg_literal; ++k, ++i) {
616 const int c = ct_buffer_[i];
617 const int64_t v0 = Violation(c);
618 if (delta == -1) {
619 num_false_enforcement_[c]--;
620 DCHECK_GE(num_false_enforcement_[c], 0);
621 if (num_false_enforcement_[c] == 0) {
622 UpdateScoreOnNewlyEnforced(c, weights[c], jump_deltas, jump_scores);
623 } else if (num_false_enforcement_[c] == 1) {
624 const double enforcement_change =
625 weights[c] * static_cast<double>(distances_[c]);
626 UpdateScoreOfEnforcementIncrease(c, enforcement_change, jump_deltas,
627 jump_scores);
628 }
629 } else {
630 num_false_enforcement_[c]++;
631 if (num_false_enforcement_[c] == 1) {
632 UpdateScoreOnNewlyUnenforced(c, weights[c], jump_deltas, jump_scores);
633 } else if (num_false_enforcement_[c] == 2) {
634 const double enforcement_change =
635 weights[c] * static_cast<double>(distances_[c]);
636 UpdateScoreOfEnforcementIncrease(c, -enforcement_change, jump_deltas,
637 jump_scores);
638 }
639 }
640 const int64_t v1 = Violation(c);
641 is_violated_[c] = v1 > 0;
642 if (v1 != v0) {
643 constraints_with_changed_violation->push_back(c);
644 }
645 }
646 int j = data.linear_start;
647 num_ops_ += 2 * data.num_linear_entries;
648 for (int k = 0; k < data.num_linear_entries; ++k, ++i, ++j) {
649 const int c = ct_buffer_[i];
650 const int64_t v0 = Violation(c);
651 const int64_t coeff = coeff_buffer_[j];
652
653 if (num_false_enforcement_[c] == 1) {
654 // Only the 1 -> 0 are impacted.
655 // This is the same as the 1->2 transition, but the old 1->0 needs to
656 // be changed from - weight * distance to - weight * new_distance.
657 const int64_t new_distance =
658 domains_[c].Distance(activities_[c] + coeff * delta);
659 if (new_distance != distances_[c]) {
660 UpdateScoreOfEnforcementIncrease(
661 c, -weights[c] * static_cast<double>(distances_[c] - new_distance),
662 jump_deltas, jump_scores);
663 }
664 } else if (num_false_enforcement_[c] == 0) {
665 UpdateScoreOnActivityChange(c, weights[c], coeff * delta, jump_deltas,
666 jump_scores);
667 }
668
669 activities_[c] += coeff * delta;
670 distances_[c] = domains_[c].Distance(activities_[c]);
671 const int64_t v1 = Violation(c);
672 is_violated_[c] = v1 > 0;
673 if (v1 != v0) {
674 constraints_with_changed_violation->push_back(c);
675 }
676 }
677}
678
680 return activities_[c];
681}
682
684 return num_false_enforcement_[c] > 0 ? 0 : distances_[c];
685}
686
688 DCHECK_EQ(is_violated_[c], Violation(c) > 0);
689 return is_violated_[c];
690}
691
692bool LinearIncrementalEvaluator::ReduceBounds(int c, int64_t lb, int64_t ub) {
693 if (domains_[c].Min() >= lb && domains_[c].Max() <= ub) return false;
694 domains_[c] = domains_[c].IntersectionWith(Domain(lb, ub));
695 distances_[c] = domains_[c].Distance(activities_[c]);
696 return true;
697}
698
700 absl::Span<const double> weights) const {
701 double result = 0.0;
702 DCHECK_GE(weights.size(), num_constraints_);
703 for (int c = 0; c < num_constraints_; ++c) {
704 if (num_false_enforcement_[c] > 0) continue;
705 result += weights[c] * static_cast<double>(distances_[c]);
706 }
707 return result;
708}
709
710// Most of the time is spent in this function.
711//
712// TODO(user): We can safely abort early if we know that delta will be >= 0.
713// TODO(user): Maybe we can compute an absolute value instead of removing
714// old_distance.
716 absl::Span<const double> weights, int var, int64_t delta) const {
717 DCHECK_NE(delta, 0);
718 if (var >= columns_.size()) return 0.0;
719 const SpanData& data = columns_[var];
720
721 int i = data.start;
722 double result = 0.0;
723 num_ops_ += data.num_pos_literal;
724 for (int k = 0; k < data.num_pos_literal; ++k, ++i) {
725 const int c = ct_buffer_[i];
726 if (num_false_enforcement_[c] == 0) {
727 // Since delta != 0, we are sure this is an enforced -> unenforced change.
728 DCHECK_EQ(delta, -1);
729 result -= weights[c] * static_cast<double>(distances_[c]);
730 } else {
731 if (delta == 1 && num_false_enforcement_[c] == 1) {
732 result += weights[c] * static_cast<double>(distances_[c]);
733 }
734 }
735 }
736
737 num_ops_ += data.num_neg_literal;
738 for (int k = 0; k < data.num_neg_literal; ++k, ++i) {
739 const int c = ct_buffer_[i];
740 if (num_false_enforcement_[c] == 0) {
741 // Since delta != 0, we are sure this is an enforced -> unenforced change.
742 DCHECK_EQ(delta, 1);
743 result -= weights[c] * static_cast<double>(distances_[c]);
744 } else {
745 if (delta == -1 && num_false_enforcement_[c] == 1) {
746 result += weights[c] * static_cast<double>(distances_[c]);
747 }
748 }
749 }
750
751 int j = data.linear_start;
752 num_ops_ += 2 * data.num_linear_entries;
753 for (int k = 0; k < data.num_linear_entries; ++k, ++i, ++j) {
754 const int c = ct_buffer_[i];
755 if (num_false_enforcement_[c] > 0) continue;
756 const int64_t coeff = coeff_buffer_[j];
757 const int64_t old_distance = distances_[c];
758 const int64_t new_distance =
759 domains_[c].Distance(activities_[c] + coeff * delta);
760 result += weights[c] * static_cast<double>(new_distance - old_distance);
761 }
762
763 return result;
764}
765
767 if (var >= columns_.size()) return false;
768 for (const int c : VarToConstraints(var)) {
769 if (Violation(c) > 0) return true;
770 }
771 return false;
772}
773
775 int var, int64_t current_value, const Domain& var_domain) const {
776 std::vector<int64_t> result = var_domain.FlattenedIntervals();
777 if (var_domain.Size() <= 2 || var >= columns_.size()) return result;
778
779 const SpanData& data = columns_[var];
780 int i = data.start + data.num_pos_literal + data.num_neg_literal;
781 int j = data.linear_start;
782 for (int k = 0; k < data.num_linear_entries; ++k, ++i, ++j) {
783 const int c = ct_buffer_[i];
784 if (num_false_enforcement_[c] > 0) continue;
785
786 // We only consider min / max.
787 // There is a change when we cross the slack.
788 // TODO(user): Deal with holes?
789 const int64_t coeff = coeff_buffer_[j];
790 const int64_t activity = activities_[c] - current_value * coeff;
791
792 const int64_t slack_min = CapSub(domains_[c].Min(), activity);
793 const int64_t slack_max = CapSub(domains_[c].Max(), activity);
794 if (slack_min != std::numeric_limits<int64_t>::min()) {
795 const int64_t ceil_bp = MathUtil::CeilOfRatio(slack_min, coeff);
796 if (ceil_bp != result.back() && var_domain.Contains(ceil_bp)) {
797 result.push_back(ceil_bp);
798 }
799 const int64_t floor_bp = MathUtil::FloorOfRatio(slack_min, coeff);
800 if (floor_bp != result.back() && var_domain.Contains(floor_bp)) {
801 result.push_back(floor_bp);
802 }
803 }
804 if (slack_min != slack_max &&
805 slack_max != std::numeric_limits<int64_t>::min()) {
806 const int64_t ceil_bp = MathUtil::CeilOfRatio(slack_max, coeff);
807 if (ceil_bp != result.back() && var_domain.Contains(ceil_bp)) {
808 result.push_back(ceil_bp);
809 }
810 const int64_t floor_bp = MathUtil::FloorOfRatio(slack_max, coeff);
811 if (floor_bp != result.back() && var_domain.Contains(floor_bp)) {
812 result.push_back(floor_bp);
813 }
814 }
815 }
816
818 return result;
819}
820
822 absl::Span<const int64_t> var_max_variation) {
823 creation_phase_ = false;
824 if (num_constraints_ == 0) return;
825
826 // Compute the total size.
827 // Note that at this point the constraint indices are not "encoded" yet.
828 int total_size = 0;
829 int total_linear_size = 0;
830 tmp_row_sizes_.assign(num_constraints_, 0);
831 tmp_row_num_positive_literals_.assign(num_constraints_, 0);
832 tmp_row_num_negative_literals_.assign(num_constraints_, 0);
833 tmp_row_num_linear_entries_.assign(num_constraints_, 0);
834 for (const auto& column : literal_entries_) {
835 total_size += column.size();
836 for (const auto [c, is_positive] : column) {
837 tmp_row_sizes_[c]++;
838 if (is_positive) {
839 tmp_row_num_positive_literals_[c]++;
840 } else {
841 tmp_row_num_negative_literals_[c]++;
842 }
843 }
844 }
845
846 row_max_variations_.assign(num_constraints_, 0);
847 for (int var = 0; var < var_entries_.size(); ++var) {
848 const int64_t range = var_max_variation[var];
849 const auto& column = var_entries_[var];
850 total_size += column.size();
851 total_linear_size += column.size();
852 for (const auto [c, coeff] : column) {
853 tmp_row_sizes_[c]++;
854 tmp_row_num_linear_entries_[c]++;
855 row_max_variations_[c] =
856 std::max(row_max_variations_[c], range * std::abs(coeff));
857 }
858 }
859
860 // Compactify for faster WeightedViolationDelta().
861 ct_buffer_.reserve(total_size);
862 coeff_buffer_.reserve(total_linear_size);
863 columns_.resize(std::max(literal_entries_.size(), var_entries_.size()));
864 for (int var = 0; var < columns_.size(); ++var) {
865 columns_[var].start = static_cast<int>(ct_buffer_.size());
866 columns_[var].linear_start = static_cast<int>(coeff_buffer_.size());
867 if (var < literal_entries_.size()) {
868 for (const auto [c, is_positive] : literal_entries_[var]) {
869 if (is_positive) {
870 columns_[var].num_pos_literal++;
871 ct_buffer_.push_back(c);
872 }
873 }
874 for (const auto [c, is_positive] : literal_entries_[var]) {
875 if (!is_positive) {
876 columns_[var].num_neg_literal++;
877 ct_buffer_.push_back(c);
878 }
879 }
880 }
881 if (var < var_entries_.size()) {
882 for (const auto [c, coeff] : var_entries_[var]) {
883 columns_[var].num_linear_entries++;
884 ct_buffer_.push_back(c);
885 coeff_buffer_.push_back(coeff);
886 }
887 }
888 }
889
890 // We do not need var_entries_ or literal_entries_ anymore.
891 //
892 // TODO(user): We could delete them before. But at the time of this
893 // optimization, I didn't want to change the behavior of the algorithm at all.
894 gtl::STLClearObject(&var_entries_);
895 gtl::STLClearObject(&literal_entries_);
896
897 // Initialize the SpanData.
898 // Transform tmp_row_sizes_ to starts in the row_var_buffer_.
899 // Transform tmp_row_num_linear_entries_ to starts in the row_coeff_buffer_.
900 int offset = 0;
901 int linear_offset = 0;
902 rows_.resize(num_constraints_);
903 for (int c = 0; c < num_constraints_; ++c) {
904 rows_[c].num_pos_literal = tmp_row_num_positive_literals_[c];
905 rows_[c].num_neg_literal = tmp_row_num_negative_literals_[c];
906 rows_[c].num_linear_entries = tmp_row_num_linear_entries_[c];
907
908 rows_[c].start = offset;
909 offset += tmp_row_sizes_[c];
910 tmp_row_sizes_[c] = rows_[c].start;
911
912 rows_[c].linear_start = linear_offset;
913 linear_offset += tmp_row_num_linear_entries_[c];
914 tmp_row_num_linear_entries_[c] = rows_[c].linear_start;
915 }
916 DCHECK_EQ(offset, total_size);
917 DCHECK_EQ(linear_offset, total_linear_size);
918
919 // Copy data.
920 row_var_buffer_.resize(total_size);
921 row_coeff_buffer_.resize(total_linear_size);
922 for (int var = 0; var < columns_.size(); ++var) {
923 const SpanData& data = columns_[var];
924 int i = data.start;
925 for (int k = 0; k < data.num_pos_literal; ++i, ++k) {
926 const int c = ct_buffer_[i];
927 row_var_buffer_[tmp_row_sizes_[c]++] = var;
928 }
929 }
930 for (int var = 0; var < columns_.size(); ++var) {
931 const SpanData& data = columns_[var];
932 int i = data.start + data.num_pos_literal;
933 for (int k = 0; k < data.num_neg_literal; ++i, ++k) {
934 const int c = ct_buffer_[i];
935 row_var_buffer_[tmp_row_sizes_[c]++] = var;
936 }
937 }
938 for (int var = 0; var < columns_.size(); ++var) {
939 const SpanData& data = columns_[var];
940 int i = data.start + data.num_pos_literal + data.num_neg_literal;
941 int j = data.linear_start;
942 for (int k = 0; k < data.num_linear_entries; ++i, ++j, ++k) {
943 const int c = ct_buffer_[i];
944 row_var_buffer_[tmp_row_sizes_[c]++] = var;
945 row_coeff_buffer_[tmp_row_num_linear_entries_[c]++] = coeff_buffer_[j];
946 }
947 }
948
949 cached_deltas_.assign(columns_.size(), 0);
950 cached_scores_.assign(columns_.size(), 0);
951 last_affected_variables_.ClearAndReserve(columns_.size());
952}
953
955 for (const int c : VarToConstraints(var)) {
956 if (domains_[c].NumIntervals() > 2) return false;
957 }
958 return true;
959}
960
961// ----- CompiledConstraint -----
962
964 absl::Span<const int64_t> solution) {
966}
967
969 int var, int64_t old_value,
970 absl::Span<const int64_t> solution_with_new_value) {
971 violation_ += ViolationDelta(var, old_value, solution_with_new_value);
972}
973
975 absl::Span<const int64_t> solution) {
977}
978
979// ----- CompiledConstraintWithProto -----
980
982 const ConstraintProto& ct_proto)
983 : ct_proto_(ct_proto) {}
984
986 const CpModelProto& model_proto) const {
987 std::vector<int> result = sat::UsedVariables(ct_proto_);
988 for (const int i_var : UsedIntervals(ct_proto_)) {
989 const ConstraintProto& interval_proto = model_proto.constraints(i_var);
990 for (const int var : sat::UsedVariables(interval_proto)) {
991 result.push_back(var);
992 }
993 }
995 result.shrink_to_fit();
996 return result;
997}
998
999// ----- CompiledBoolXorConstraint -----
1000
1004
1006 absl::Span<const int64_t> solution) {
1007 int64_t sum_of_literals = 0;
1008 for (const int lit : ct_proto().bool_xor().literals()) {
1009 sum_of_literals += LiteralValue(lit, solution);
1010 }
1011 return 1 - (sum_of_literals % 2);
1012}
1013
1015 int /*var*/, int64_t /*old_value*/,
1016 absl::Span<const int64_t> /*solution_with_new_value*/) {
1017 return violation() == 0 ? 1 : -1;
1018}
1019
1020// ----- CompiledLinMaxConstraint -----
1021
1025
1027 absl::Span<const int64_t> solution) {
1028 const int64_t target_value =
1029 ExprValue(ct_proto().lin_max().target(), solution);
1030 int64_t max_of_expressions = std::numeric_limits<int64_t>::min();
1031 for (const LinearExpressionProto& expr : ct_proto().lin_max().exprs()) {
1032 const int64_t expr_value = ExprValue(expr, solution);
1033 max_of_expressions = std::max(max_of_expressions, expr_value);
1034 }
1035 return std::max(target_value - max_of_expressions, int64_t{0});
1036}
1037
1038// ----- CompiledIntProdConstraint -----
1039
1043
1045 absl::Span<const int64_t> solution) {
1046 const int64_t target_value =
1047 ExprValue(ct_proto().int_prod().target(), solution);
1048 int64_t prod_value = 1;
1049 for (const LinearExpressionProto& expr : ct_proto().int_prod().exprs()) {
1050 prod_value *= ExprValue(expr, solution);
1051 }
1052 return std::abs(target_value - prod_value);
1053}
1054
1055// ----- CompiledIntDivConstraint -----
1056
1060
1062 absl::Span<const int64_t> solution) {
1063 const int64_t target_value =
1064 ExprValue(ct_proto().int_div().target(), solution);
1065 DCHECK_EQ(ct_proto().int_div().exprs_size(), 2);
1066 const int64_t div_value = ExprValue(ct_proto().int_div().exprs(0), solution) /
1067 ExprValue(ct_proto().int_div().exprs(1), solution);
1068 return std::abs(target_value - div_value);
1069}
1070
1071// ----- CompiledIntModConstraint -----
1072
1076
1078 absl::Span<const int64_t> solution) {
1079 const int64_t target_value =
1080 ExprValue(ct_proto().int_mod().target(), solution);
1081 DCHECK_EQ(ct_proto().int_mod().exprs_size(), 2);
1082 // Note: The violation computation assumes the modulo is constant.
1083 const int64_t expr_value = ExprValue(ct_proto().int_mod().exprs(0), solution);
1084 const int64_t mod_value = ExprValue(ct_proto().int_mod().exprs(1), solution);
1085 const int64_t rhs = expr_value % mod_value;
1086 if ((expr_value >= 0 && target_value >= 0) ||
1087 (expr_value <= 0 && target_value <= 0)) {
1088 // Easy case.
1089 return std::min({std::abs(target_value - rhs),
1090 std::abs(target_value) + std::abs(mod_value - rhs),
1091 std::abs(rhs) + std::abs(mod_value - target_value)});
1092 } else {
1093 // Different signs.
1094 // We use the sum of the absolute value to have a better gradiant.
1095 // We could also use the min of target_move and the expr_move.
1096 return std::abs(target_value) + std::abs(expr_value);
1097 }
1098}
1099
1100// ----- CompiledAllDiffConstraint -----
1101
1105
1107 absl::Span<const int64_t> solution) {
1108 values_.clear();
1109 for (const LinearExpressionProto& expr : ct_proto().all_diff().exprs()) {
1110 values_.push_back(ExprValue(expr, solution));
1111 }
1112 std::sort(values_.begin(), values_.end());
1113
1114 int64_t value = values_[0];
1115 int counter = 1;
1116 int64_t violation = 0;
1117 for (int i = 1; i < values_.size(); ++i) {
1118 const int64_t new_value = values_[i];
1119 if (new_value == value) {
1120 counter++;
1121 } else {
1122 violation += counter * (counter - 1) / 2;
1123 counter = 1;
1124 value = new_value;
1125 }
1126 }
1127 violation += counter * (counter - 1) / 2;
1128 return violation;
1129}
1130
1131// ----- NoOverlapBetweenTwoIntervals -----
1132
1134 int interval_0, int interval_1, const CpModelProto& cp_model) {
1135 const ConstraintProto& ct0 = cp_model.constraints(interval_0);
1136 const ConstraintProto& ct1 = cp_model.constraints(interval_1);
1137
1138 // The more compact the better, hence the size + int[].
1139 num_enforcements_ =
1140 ct0.enforcement_literal().size() + ct1.enforcement_literal().size();
1141 if (num_enforcements_ > 0) {
1142 enforcements_.reset(new int[num_enforcements_]);
1143 int i = 0;
1144 for (const int lit : ct0.enforcement_literal()) enforcements_[i++] = lit;
1145 for (const int lit : ct1.enforcement_literal()) enforcements_[i++] = lit;
1146 }
1147
1148 // We prefer to use start + size instead of end so that moving "start" moves
1149 // the whole interval around (for the non-fixed duration case).
1150 end_minus_start_1_ =
1151 ExprDiff(LinearExprSum(ct0.interval().start(), ct0.interval().size()),
1152 ct1.interval().start());
1153 end_minus_start_2_ =
1154 ExprDiff(LinearExprSum(ct1.interval().start(), ct1.interval().size()),
1155 ct0.interval().start());
1156}
1157
1158// Same as NoOverlapMinRepairDistance().
1159int64_t NoOverlapBetweenTwoIntervals::ComputeViolationInternal(
1160 absl::Span<const int64_t> solution) {
1161 for (int i = 0; i < num_enforcements_; ++i) {
1162 if (!LiteralValue(enforcements_[i], solution)) return 0;
1163 }
1164 const int64_t diff1 = ExprValue(end_minus_start_1_, solution);
1165 const int64_t diff2 = ExprValue(end_minus_start_2_, solution);
1166 return std::max(std::min(diff1, diff2), int64_t{0});
1167}
1168
1170 const CpModelProto& /*model_proto*/) const {
1171 std::vector<int> result;
1172 for (int i = 0; i < num_enforcements_; ++i) {
1173 result.push_back(PositiveRef(enforcements_[i]));
1174 }
1175 for (const int var : end_minus_start_1_.vars()) {
1176 result.push_back(PositiveRef(var));
1177 }
1178 for (const int var : end_minus_start_2_.vars()) {
1179 result.push_back(PositiveRef(var));
1180 }
1182 result.shrink_to_fit();
1183 return result;
1184}
1185
1186// ----- CompiledNoOverlap2dConstraint -----
1187
1188int64_t OverlapOfTwoIntervals(const ConstraintProto& interval1,
1189 const ConstraintProto& interval2,
1190 absl::Span<const int64_t> solution) {
1191 for (const int lit : interval1.enforcement_literal()) {
1192 if (!LiteralValue(lit, solution)) return 0;
1193 }
1194 for (const int lit : interval2.enforcement_literal()) {
1195 if (!LiteralValue(lit, solution)) return 0;
1196 }
1197
1198 const int64_t start1 = ExprValue(interval1.interval().start(), solution);
1199 const int64_t end1 = ExprValue(interval1.interval().end(), solution);
1200
1201 const int64_t start2 = ExprValue(interval2.interval().start(), solution);
1202 const int64_t end2 = ExprValue(interval2.interval().end(), solution);
1203
1204 if (start1 >= end2 || start2 >= end1) return 0; // Disjoint.
1205
1206 // We force a min cost of 1 to cover the case where a interval of size 0 is in
1207 // the middle of another interval.
1208 return std::max(std::min(std::min(end2 - start2, end1 - start1),
1209 std::min(end2 - start1, end1 - start2)),
1210 int64_t{1});
1211}
1212
1213int64_t NoOverlapMinRepairDistance(const ConstraintProto& interval1,
1214 const ConstraintProto& interval2,
1215 absl::Span<const int64_t> solution) {
1216 for (const int lit : interval1.enforcement_literal()) {
1217 if (!LiteralValue(lit, solution)) return 0;
1218 }
1219 for (const int lit : interval2.enforcement_literal()) {
1220 if (!LiteralValue(lit, solution)) return 0;
1221 }
1222
1223 const int64_t start1 = ExprValue(interval1.interval().start(), solution);
1224 const int64_t end1 = ExprValue(interval1.interval().end(), solution);
1225
1226 const int64_t start2 = ExprValue(interval2.interval().start(), solution);
1227 const int64_t end2 = ExprValue(interval2.interval().end(), solution);
1228
1229 return std::max(std::min(end2 - start1, end1 - start2), int64_t{0});
1230}
1231
1233 const ConstraintProto& ct_proto, const CpModelProto& cp_model)
1234 : CompiledConstraintWithProto(ct_proto), cp_model_(cp_model) {}
1235
1237 absl::Span<const int64_t> solution) {
1238 DCHECK_GE(ct_proto().no_overlap_2d().x_intervals_size(), 2);
1239 const int size = ct_proto().no_overlap_2d().x_intervals_size();
1240
1241 int64_t violation = 0;
1242 for (int i = 0; i + 1 < size; ++i) {
1243 const ConstraintProto& x_i =
1244 cp_model_.constraints(ct_proto().no_overlap_2d().x_intervals(i));
1245 const ConstraintProto& y_i =
1246 cp_model_.constraints(ct_proto().no_overlap_2d().y_intervals(i));
1247 for (int j = i + 1; j < size; ++j) {
1248 const ConstraintProto& x_j =
1249 cp_model_.constraints(ct_proto().no_overlap_2d().x_intervals(j));
1250 const ConstraintProto& y_j =
1251 cp_model_.constraints(ct_proto().no_overlap_2d().y_intervals(j));
1252
1253 // TODO(user): Experiment with
1254 // violation +=
1255 // std::max(std::min(NoOverlapMinRepairDistance(x_i, x_j, solution),
1256 // NoOverlapMinRepairDistance(y_i, y_j, solution)),
1257 // int64_t{0});
1258 // Currently, the effect is unclear on 2d packing problems.
1259 violation +=
1260 std::max(std::min(NoOverlapMinRepairDistance(x_i, x_j, solution) *
1263 OverlapOfTwoIntervals(x_i, x_j, solution)),
1264 int64_t{0});
1265 }
1266 }
1267 return violation;
1268}
1269
1270// ----- CompiledCircuitConstraint -----
1271
1272// The violation of a circuit has three parts:
1273// 1. Flow imbalance, maintained by the linear part.
1274// 2. The number of non-skipped SCCs in the graph minus 1.
1275// 3. The number of non-skipped SCCs that cannot be reached from any other
1276// component minus 1.
1277//
1278// #3 is not necessary for correctness, but makes the function much smoother.
1279//
1280// The only difference between single and multi circuit is flow balance at the
1281// depot, so we use the same compiled constraint for both.
1283 public:
1284 explicit CompiledCircuitConstraint(const ConstraintProto& ct_proto);
1285 ~CompiledCircuitConstraint() override = default;
1286
1287 int64_t ComputeViolation(absl::Span<const int64_t> solution) override;
1288 void PerformMove(int var, int64_t old_value,
1289 absl::Span<const int64_t> new_solution) override;
1290 int64_t ViolationDelta(
1291 int var, int64_t old_value,
1292 absl::Span<const int64_t> solution_with_new_value) override;
1293
1294 private:
1295 struct SccOutput {
1296 void emplace_back(const int* start, const int* end);
1297 void reset(int num_nodes);
1298
1299 int num_components = 0;
1300 std::vector<bool> skipped;
1301 std::vector<int> root;
1302 };
1303 void InitGraph(absl::Span<const int64_t> solution);
1304 bool UpdateGraph(int var, int64_t value);
1305 int64_t ViolationForCurrentGraph();
1306
1307 absl::flat_hash_map<int, std::vector<int>> arcs_by_lit_;
1308 absl::Span<const int> literals_;
1309 absl::Span<const int> tails_;
1310 absl::Span<const int> heads_;
1311 // Stores the currently active arcs per tail node.
1312 std::vector<DenseSet<int>> graph_;
1313 SccOutput sccs_;
1314 SccOutput committed_sccs_;
1315 std::vector<bool> has_in_arc_;
1317 scc_finder_;
1318};
1319
1320void CompiledCircuitConstraint::SccOutput::emplace_back(int const* start,
1321 int const* end) {
1322 const int root_node = *start;
1323 const int size = end - start;
1324 if (size > 1) {
1325 ++num_components;
1326 }
1327 for (; start != end; ++start) {
1328 root[*start] = root_node;
1329 skipped[*start] = (size == 1);
1330 }
1331}
1332void CompiledCircuitConstraint::SccOutput::reset(int num_nodes) {
1333 num_components = 0;
1334 root.clear();
1335 root.resize(num_nodes);
1336 skipped.clear();
1337 skipped.resize(num_nodes);
1338}
1339
1341 const ConstraintProto& ct_proto)
1343 const bool routes = ct_proto.has_routes();
1344 tails_ = routes ? ct_proto.routes().tails() : ct_proto.circuit().tails();
1345 heads_ = absl::MakeConstSpan(routes ? ct_proto.routes().heads()
1346 : ct_proto.circuit().heads());
1347 literals_ = absl::MakeConstSpan(routes ? ct_proto.routes().literals()
1348 : ct_proto.circuit().literals());
1349 graph_.resize(*absl::c_max_element(tails_) + 1);
1350 for (int i = 0; i < literals_.size(); ++i) {
1351 arcs_by_lit_[literals_[i]].push_back(i);
1352 }
1353}
1354
1355void CompiledCircuitConstraint::InitGraph(absl::Span<const int64_t> solution) {
1356 for (DenseSet<int>& edges : graph_) {
1357 edges.clear();
1358 }
1359 for (int i = 0; i < tails_.size(); ++i) {
1360 if (!LiteralValue(literals_[i], solution)) continue;
1361 graph_[tails_[i]].insert(heads_[i]);
1362 }
1363}
1364
1365bool CompiledCircuitConstraint::UpdateGraph(int var, int64_t value) {
1366 bool needs_update = false;
1367 const int enabled_lit =
1368 value != 0 ? PositiveRef(var) : NegatedRef(PositiveRef(var));
1369 const int disabled_lit = NegatedRef(enabled_lit);
1370 for (const int arc : arcs_by_lit_[disabled_lit]) {
1371 const int tail = tails_[arc];
1372 const int head = heads_[arc];
1373 // Removing a self arc cannot change violation.
1374 needs_update = needs_update || tail != head;
1375 graph_[tails_[arc]].erase(heads_[arc]);
1376 }
1377 for (const int arc : arcs_by_lit_[enabled_lit]) {
1378 const int tail = tails_[arc];
1379 const int head = heads_[arc];
1380 // Adding an arc can only change violation if it connects new SCCs.
1381 needs_update = needs_update ||
1382 committed_sccs_.root[tail] != committed_sccs_.root[head];
1383 graph_[tails_[arc]].insert(heads_[arc]);
1384 }
1385 return needs_update;
1386}
1387
1389 int var, int64_t, absl::Span<const int64_t> new_solution) {
1390 UpdateGraph(var, new_solution[var]);
1391 violation_ = ViolationForCurrentGraph();
1392 std::swap(committed_sccs_, sccs_);
1393}
1394
1396 absl::Span<const int64_t> solution) {
1397 InitGraph(solution);
1398 int64_t result = ViolationForCurrentGraph();
1399 std::swap(committed_sccs_, sccs_);
1400 return result;
1401}
1402
1404 int var, int64_t old_value,
1405 absl::Span<const int64_t> solution_with_new_value) {
1406 int64_t result = 0;
1407 if (UpdateGraph(var, solution_with_new_value[var])) {
1408 result = ViolationForCurrentGraph() - violation_;
1409 }
1410 UpdateGraph(var, old_value);
1411 return result;
1412}
1413
1414int64_t CompiledCircuitConstraint::ViolationForCurrentGraph() {
1415 const int num_nodes = graph_.size();
1416 sccs_.reset(num_nodes);
1417 scc_finder_.FindStronglyConnectedComponents(num_nodes, graph_, &sccs_);
1418 // Skipping all nodes causes off-by-one errors below, so it's simpler to
1419 // handle explicitly.
1420 if (sccs_.num_components == 0) return 0;
1421 // Count the number of SCCs that have inbound cross-component arcs
1422 // as a smoother measure of progress towards strong connectivity.
1423 int num_half_connected_components = 0;
1424 has_in_arc_.clear();
1425 has_in_arc_.resize(num_nodes, false);
1426 for (int tail = 0; tail < graph_.size(); ++tail) {
1427 if (sccs_.skipped[tail]) continue;
1428 for (const int head : graph_[tail]) {
1429 const int head_root = sccs_.root[head];
1430 if (sccs_.root[tail] == head_root) continue;
1431 if (has_in_arc_[head_root]) continue;
1432 if (sccs_.skipped[head_root]) continue;
1433 has_in_arc_[head_root] = true;
1434 ++num_half_connected_components;
1435 }
1436 }
1437 const int64_t violation = sccs_.num_components - 1 + sccs_.num_components -
1438 num_half_connected_components - 1 +
1439 (ct_proto().has_routes() ? sccs_.skipped[0] : 0);
1440 VLOG(2) << "#SCCs=" << sccs_.num_components << " #nodes=" << num_nodes
1441 << " #half_connected_components=" << num_half_connected_components
1442 << " violation=" << violation;
1443 return violation;
1444}
1445
1447 const ConstraintProto& ct_proto) {
1448 const bool routes = ct_proto.has_routes();
1449 auto heads = routes ? ct_proto.routes().heads() : ct_proto.circuit().heads();
1450 auto tails = routes ? ct_proto.routes().tails() : ct_proto.circuit().tails();
1451 auto literals =
1452 routes ? ct_proto.routes().literals() : ct_proto.circuit().literals();
1453
1454 std::vector<std::vector<int>> inflow_lits;
1455 std::vector<std::vector<int>> outflow_lits;
1456 for (int i = 0; i < heads.size(); ++i) {
1457 if (heads[i] >= inflow_lits.size()) {
1458 inflow_lits.resize(heads[i] + 1);
1459 }
1460 inflow_lits[heads[i]].push_back(literals[i]);
1461 if (tails[i] >= outflow_lits.size()) {
1462 outflow_lits.resize(tails[i] + 1);
1463 }
1464 outflow_lits[tails[i]].push_back(literals[i]);
1465 }
1466 if (routes) {
1467 const int depot_net_flow = linear_evaluator.NewConstraint({0, 0});
1468 for (const int lit : inflow_lits[0]) {
1469 linear_evaluator.AddLiteral(depot_net_flow, lit, 1);
1470 }
1471 for (const int lit : outflow_lits[0]) {
1472 linear_evaluator.AddLiteral(depot_net_flow, lit, -1);
1473 }
1474 }
1475 for (int i = routes ? 1 : 0; i < inflow_lits.size(); ++i) {
1476 const int inflow_ct = linear_evaluator.NewConstraint({1, 1});
1477 for (const int lit : inflow_lits[i]) {
1478 linear_evaluator.AddLiteral(inflow_ct, lit);
1479 }
1480 }
1481 for (int i = routes ? 1 : 0; i < outflow_lits.size(); ++i) {
1482 const int outflow_ct = linear_evaluator.NewConstraint({1, 1});
1483 for (const int lit : outflow_lits[i]) {
1484 linear_evaluator.AddLiteral(outflow_ct, lit);
1485 }
1486 }
1487}
1488
1489// ----- LsEvaluator -----
1490
1491LsEvaluator::LsEvaluator(const CpModelProto& cp_model,
1492 const SatParameters& params)
1493 : cp_model_(cp_model), params_(params) {
1494 var_to_constraints_.resize(cp_model_.variables_size());
1495 jump_value_optimal_.resize(cp_model_.variables_size(), true);
1496 num_violated_constraint_per_var_ignoring_objective_.assign(
1497 cp_model_.variables_size(), 0);
1498
1499 std::vector<bool> ignored_constraints(cp_model_.constraints_size(), false);
1500 std::vector<ConstraintProto> additional_constraints;
1501 CompileConstraintsAndObjective(ignored_constraints, additional_constraints);
1502 BuildVarConstraintGraph();
1503 violated_constraints_.reserve(NumEvaluatorConstraints());
1504}
1505
1507 const CpModelProto& cp_model, const SatParameters& params,
1508 const std::vector<bool>& ignored_constraints,
1509 const std::vector<ConstraintProto>& additional_constraints)
1510 : cp_model_(cp_model), params_(params) {
1511 var_to_constraints_.resize(cp_model_.variables_size());
1512 jump_value_optimal_.resize(cp_model_.variables_size(), true);
1513 num_violated_constraint_per_var_ignoring_objective_.assign(
1514 cp_model_.variables_size(), 0);
1515 CompileConstraintsAndObjective(ignored_constraints, additional_constraints);
1516 BuildVarConstraintGraph();
1517 violated_constraints_.reserve(NumEvaluatorConstraints());
1518}
1519
1520void LsEvaluator::BuildVarConstraintGraph() {
1521 // Clear the var <-> constraint graph.
1522 for (std::vector<int>& ct_indices : var_to_constraints_) ct_indices.clear();
1523 constraint_to_vars_.resize(constraints_.size());
1524
1525 // Build the var <-> constraint graph.
1526 for (int ct_index = 0; ct_index < constraints_.size(); ++ct_index) {
1527 constraint_to_vars_[ct_index] =
1528 constraints_[ct_index]->UsedVariables(cp_model_);
1529 for (const int var : constraint_to_vars_[ct_index]) {
1530 var_to_constraints_[var].push_back(ct_index);
1531 }
1532 }
1533
1534 // Remove duplicates.
1535 for (std::vector<int>& constraints : var_to_constraints_) {
1536 gtl::STLSortAndRemoveDuplicates(&constraints);
1537 }
1538 for (std::vector<int>& vars : constraint_to_vars_) {
1540 }
1541
1542 // Scan the model to decide if a variable is linked to a convex evaluation.
1543 jump_value_optimal_.resize(cp_model_.variables_size());
1544 for (int i = 0; i < cp_model_.variables_size(); ++i) {
1545 if (!var_to_constraints_[i].empty()) {
1546 jump_value_optimal_[i] = false;
1547 continue;
1548 }
1549
1550 const IntegerVariableProto& var_proto = cp_model_.variables(i);
1551 if (var_proto.domain_size() == 2 && var_proto.domain(0) == 0 &&
1552 var_proto.domain(1) == 1) {
1553 // Boolean variables violation change is always convex.
1554 jump_value_optimal_[i] = true;
1555 continue;
1556 }
1557
1558 jump_value_optimal_[i] = linear_evaluator_.ViolationChangeIsConvex(i);
1559 }
1560}
1561
1562void LsEvaluator::CompileOneConstraint(const ConstraintProto& ct) {
1563 switch (ct.constraint_case()) {
1564 case ConstraintProto::ConstraintCase::kBoolOr: {
1565 // Encoding using enforcement literal is slightly more efficient.
1566 const int ct_index = linear_evaluator_.NewConstraint(Domain(1, 1));
1567 for (const int lit : ct.enforcement_literal()) {
1568 linear_evaluator_.AddEnforcementLiteral(ct_index, lit);
1569 }
1570 for (const int lit : ct.bool_or().literals()) {
1571 linear_evaluator_.AddEnforcementLiteral(ct_index, NegatedRef(lit));
1572 }
1573 break;
1574 }
1575 case ConstraintProto::ConstraintCase::kBoolAnd: {
1576 const int num_literals = ct.bool_and().literals_size();
1577 const int ct_index =
1578 linear_evaluator_.NewConstraint(Domain(num_literals));
1579 for (const int lit : ct.enforcement_literal()) {
1580 linear_evaluator_.AddEnforcementLiteral(ct_index, lit);
1581 }
1582 for (const int lit : ct.bool_and().literals()) {
1583 linear_evaluator_.AddLiteral(ct_index, lit);
1584 }
1585 break;
1586 }
1587 case ConstraintProto::ConstraintCase::kAtMostOne: {
1588 DCHECK(ct.enforcement_literal().empty());
1589 const int ct_index = linear_evaluator_.NewConstraint({0, 1});
1590 for (const int lit : ct.at_most_one().literals()) {
1591 linear_evaluator_.AddLiteral(ct_index, lit);
1592 }
1593 break;
1594 }
1595 case ConstraintProto::ConstraintCase::kExactlyOne: {
1596 DCHECK(ct.enforcement_literal().empty());
1597 const int ct_index = linear_evaluator_.NewConstraint({1, 1});
1598 for (const int lit : ct.exactly_one().literals()) {
1599 linear_evaluator_.AddLiteral(ct_index, lit);
1600 }
1601 break;
1602 }
1603 case ConstraintProto::ConstraintCase::kBoolXor: {
1604 constraints_.emplace_back(new CompiledBoolXorConstraint(ct));
1605 break;
1606 }
1607 case ConstraintProto::ConstraintCase::kAllDiff: {
1608 constraints_.emplace_back(new CompiledAllDiffConstraint(ct));
1609 break;
1610 }
1611 case ConstraintProto::ConstraintCase::kLinMax: {
1612 // This constraint is split into linear precedences and its max
1613 // maintenance.
1614 const LinearExpressionProto& target = ct.lin_max().target();
1615 for (const LinearExpressionProto& expr : ct.lin_max().exprs()) {
1616 const int64_t max_value =
1617 ExprMax(target, cp_model_) - ExprMin(expr, cp_model_);
1618 const int precedence_index =
1619 linear_evaluator_.NewConstraint({0, max_value});
1620 linear_evaluator_.AddLinearExpression(precedence_index, target, 1);
1621 linear_evaluator_.AddLinearExpression(precedence_index, expr, -1);
1622 }
1623
1624 // Penalty when target > all expressions.
1625 constraints_.emplace_back(new CompiledLinMaxConstraint(ct));
1626 break;
1627 }
1628 case ConstraintProto::ConstraintCase::kIntProd: {
1629 constraints_.emplace_back(new CompiledIntProdConstraint(ct));
1630 break;
1631 }
1632 case ConstraintProto::ConstraintCase::kIntDiv: {
1633 constraints_.emplace_back(new CompiledIntDivConstraint(ct));
1634 break;
1635 }
1636 case ConstraintProto::ConstraintCase::kIntMod: {
1637 DCHECK_EQ(ExprMin(ct.int_mod().exprs(1), cp_model_),
1638 ExprMax(ct.int_mod().exprs(1), cp_model_));
1639 constraints_.emplace_back(new CompiledIntModConstraint(ct));
1640 break;
1641 }
1642 case ConstraintProto::ConstraintCase::kLinear: {
1643 const Domain domain = ReadDomainFromProto(ct.linear());
1644 const int ct_index = linear_evaluator_.NewConstraint(domain);
1645 for (const int lit : ct.enforcement_literal()) {
1646 linear_evaluator_.AddEnforcementLiteral(ct_index, lit);
1647 }
1648 for (int i = 0; i < ct.linear().vars_size(); ++i) {
1649 const int var = ct.linear().vars(i);
1650 const int64_t coeff = ct.linear().coeffs(i);
1651 linear_evaluator_.AddTerm(ct_index, var, coeff);
1652 }
1653 break;
1654 }
1655 case ConstraintProto::ConstraintCase::kNoOverlap: {
1656 const int size = ct.no_overlap().intervals_size();
1657 if (size <= 1) break;
1658 if (size > params_.feasibility_jump_max_expanded_constraint_size()) {
1659 // Similar code to the kCumulative constraint.
1660 // The violation will be the area above the capacity.
1661 LinearExpressionProto one;
1662 one.set_offset(1);
1663 std::vector<std::optional<int>> is_active;
1664 std::vector<LinearExpressionProto> times;
1665 std::vector<LinearExpressionProto> demands;
1666 const int num_intervals = ct.no_overlap().intervals().size();
1667 for (int i = 0; i < num_intervals; ++i) {
1668 const ConstraintProto& interval_ct =
1669 cp_model_.constraints(ct.no_overlap().intervals(i));
1670 if (interval_ct.enforcement_literal().empty()) {
1671 is_active.push_back(std::nullopt);
1672 is_active.push_back(std::nullopt);
1673 } else {
1674 CHECK_EQ(interval_ct.enforcement_literal().size(), 1);
1675 is_active.push_back(interval_ct.enforcement_literal(0));
1676 is_active.push_back(interval_ct.enforcement_literal(0));
1677 }
1678
1679 times.push_back(interval_ct.interval().start());
1680 times.push_back(LinearExprSum(interval_ct.interval().start(),
1681 interval_ct.interval().size()));
1682 demands.push_back(one);
1683 demands.push_back(NegatedLinearExpression(one));
1684 }
1685 constraints_.emplace_back(new CompiledReservoirConstraint(
1686 std::move(one), std::move(is_active), std::move(times),
1687 std::move(demands)));
1688 } else {
1689 // We expand the no_overlap constraints into a quadratic number of
1690 // disjunctions.
1691 for (int i = 0; i + 1 < size; ++i) {
1692 const IntervalConstraintProto& interval_i =
1693 cp_model_.constraints(ct.no_overlap().intervals(i)).interval();
1694 const int64_t min_start_i = ExprMin(interval_i.start(), cp_model_);
1695 const int64_t max_end_i = ExprMax(interval_i.end(), cp_model_);
1696 for (int j = i + 1; j < size; ++j) {
1697 const IntervalConstraintProto& interval_j =
1698 cp_model_.constraints(ct.no_overlap().intervals(j)).interval();
1699 const int64_t min_start_j = ExprMin(interval_j.start(), cp_model_);
1700 const int64_t max_end_j = ExprMax(interval_j.end(), cp_model_);
1701 if (min_start_i >= max_end_j || min_start_j >= max_end_i) continue;
1702
1703 constraints_.emplace_back(new NoOverlapBetweenTwoIntervals(
1704 ct.no_overlap().intervals(i), ct.no_overlap().intervals(j),
1705 cp_model_));
1706 }
1707 }
1708 }
1709 break;
1710 }
1711 case ConstraintProto::ConstraintCase::kCumulative: {
1712 LinearExpressionProto capacity = ct.cumulative().capacity();
1713 std::vector<std::optional<int>> is_active;
1714 std::vector<LinearExpressionProto> times;
1715 std::vector<LinearExpressionProto> demands;
1716 const int num_intervals = ct.cumulative().intervals().size();
1717 for (int i = 0; i < num_intervals; ++i) {
1718 const ConstraintProto& interval_ct =
1719 cp_model_.constraints(ct.cumulative().intervals(i));
1720 if (interval_ct.enforcement_literal().empty()) {
1721 is_active.push_back(std::nullopt);
1722 is_active.push_back(std::nullopt);
1723 } else {
1724 CHECK_EQ(interval_ct.enforcement_literal().size(), 1);
1725 is_active.push_back(interval_ct.enforcement_literal(0));
1726 is_active.push_back(interval_ct.enforcement_literal(0));
1727 }
1728
1729 // Start.
1730 times.push_back(interval_ct.interval().start());
1731 demands.push_back(ct.cumulative().demands(i));
1732
1733 // End.
1734 // I tried 3 alternatives: end, max(end, start+size) and just start +
1735 // size. The most performing one was "start + size" on the multi-mode
1736 // RCPSP.
1737 //
1738 // Note that for fixed size, this do not matter. It is easy enough to
1739 // try any expression by creating a small wrapper class to use instead
1740 // of a LinearExpressionProto for time.
1741 times.push_back(LinearExprSum(interval_ct.interval().start(),
1742 interval_ct.interval().size()));
1743 demands.push_back(NegatedLinearExpression(ct.cumulative().demands(i)));
1744 }
1745
1746 constraints_.emplace_back(new CompiledReservoirConstraint(
1747 std::move(capacity), std::move(is_active), std::move(times),
1748 std::move(demands)));
1749 break;
1750 }
1751 case ConstraintProto::ConstraintCase::kNoOverlap2D: {
1752 const auto& x_intervals = ct.no_overlap_2d().x_intervals();
1753 const auto& y_intervals = ct.no_overlap_2d().y_intervals();
1754 const int size = x_intervals.size();
1755 if (size <= 1) break;
1756 if (size == 2 ||
1757 size > params_.feasibility_jump_max_expanded_constraint_size()) {
1758 CompiledNoOverlap2dConstraint* no_overlap_2d =
1759 new CompiledNoOverlap2dConstraint(ct, cp_model_);
1760 constraints_.emplace_back(no_overlap_2d);
1761 break;
1762 }
1763
1764 for (int i = 0; i + 1 < size; ++i) {
1765 const IntervalConstraintProto& x_interval_i =
1766 cp_model_.constraints(x_intervals[i]).interval();
1767 const int64_t x_min_start_i = ExprMin(x_interval_i.start(), cp_model_);
1768 const int64_t x_max_end_i = ExprMax(x_interval_i.end(), cp_model_);
1769 const IntervalConstraintProto& y_interval_i =
1770 cp_model_.constraints(y_intervals[i]).interval();
1771 const int64_t y_min_start_i = ExprMin(y_interval_i.start(), cp_model_);
1772 const int64_t y_max_end_i = ExprMax(y_interval_i.end(), cp_model_);
1773 for (int j = i + 1; j < size; ++j) {
1774 const IntervalConstraintProto& x_interval_j =
1775 cp_model_.constraints(x_intervals[j]).interval();
1776 const int64_t x_min_start_j =
1777 ExprMin(x_interval_j.start(), cp_model_);
1778 const int64_t x_max_end_j = ExprMax(x_interval_j.end(), cp_model_);
1779 const IntervalConstraintProto& y_interval_j =
1780 cp_model_.constraints(y_intervals[j]).interval();
1781 const int64_t y_min_start_j =
1782 ExprMin(y_interval_j.start(), cp_model_);
1783 const int64_t y_max_end_j = ExprMax(y_interval_j.end(), cp_model_);
1784 if (x_min_start_i >= x_max_end_j || x_min_start_j >= x_max_end_i ||
1785 y_min_start_i >= y_max_end_j || y_min_start_j >= y_max_end_i) {
1786 continue;
1787 }
1788 ConstraintProto* diffn = expanded_constraints_.add_constraints();
1789 diffn->mutable_no_overlap_2d()->add_x_intervals(x_intervals[i]);
1790 diffn->mutable_no_overlap_2d()->add_x_intervals(x_intervals[j]);
1791 diffn->mutable_no_overlap_2d()->add_y_intervals(y_intervals[i]);
1792 diffn->mutable_no_overlap_2d()->add_y_intervals(y_intervals[j]);
1793 CompiledNoOverlap2dConstraint* no_overlap_2d =
1794 new CompiledNoOverlap2dConstraint(*diffn, cp_model_);
1795 constraints_.emplace_back(no_overlap_2d);
1796 }
1797 }
1798 break;
1799 }
1800 case ConstraintProto::ConstraintCase::kCircuit:
1801 case ConstraintProto::ConstraintCase::kRoutes:
1802 constraints_.emplace_back(new CompiledCircuitConstraint(ct));
1803 AddCircuitFlowConstraints(linear_evaluator_, ct);
1804 break;
1805 default:
1806 VLOG(1) << "Not implemented: " << ct.constraint_case();
1807 break;
1808 }
1809}
1810
1811void LsEvaluator::CompileConstraintsAndObjective(
1812 const std::vector<bool>& ignored_constraints,
1813 const std::vector<ConstraintProto>& additional_constraints) {
1814 constraints_.clear();
1815
1816 // The first compiled constraint is always the objective if present.
1817 if (cp_model_.has_objective()) {
1818 const int ct_index = linear_evaluator_.NewConstraint(
1819 cp_model_.objective().domain().empty()
1821 : ReadDomainFromProto(cp_model_.objective()));
1822 DCHECK_EQ(0, ct_index);
1823 for (int i = 0; i < cp_model_.objective().vars_size(); ++i) {
1824 const int var = cp_model_.objective().vars(i);
1825 const int64_t coeff = cp_model_.objective().coeffs(i);
1826 linear_evaluator_.AddTerm(ct_index, var, coeff);
1827 }
1828 }
1829
1830 for (int c = 0; c < cp_model_.constraints_size(); ++c) {
1831 if (ignored_constraints[c]) continue;
1832 CompileOneConstraint(cp_model_.constraints(c));
1833 }
1834
1835 for (const ConstraintProto& ct : additional_constraints) {
1836 CompileOneConstraint(ct);
1837 }
1838
1839 // Make sure we have access to the data in an efficient way.
1840 std::vector<int64_t> var_max_variations(cp_model_.variables().size());
1841 for (int var = 0; var < cp_model_.variables().size(); ++var) {
1842 const auto& domain = cp_model_.variables(var).domain();
1843 var_max_variations[var] = domain[domain.size() - 1] - domain[0];
1844 }
1845 linear_evaluator_.PrecomputeCompactView(var_max_variations);
1846}
1847
1848bool LsEvaluator::ReduceObjectiveBounds(int64_t lb, int64_t ub) {
1849 if (!cp_model_.has_objective()) return false;
1850 if (linear_evaluator_.ReduceBounds(/*c=*/0, lb, ub)) {
1852 return true;
1853 }
1854 return false;
1855}
1856
1857void LsEvaluator::ComputeAllViolations(absl::Span<const int64_t> solution) {
1858 // Linear constraints.
1859 linear_evaluator_.ComputeInitialActivities(solution);
1860
1861 // Generic constraints.
1862 for (const auto& ct : constraints_) {
1863 ct->InitializeViolation(solution);
1864 }
1865
1866 RecomputeViolatedList(/*linear_only=*/false);
1867}
1868
1870 absl::Span<const int64_t> solution) {
1871 // Generic constraints.
1872 for (const auto& ct : constraints_) {
1873 ct->InitializeViolation(solution);
1874 }
1875}
1876
1878 int var, int64_t old_value, absl::Span<const int64_t> new_solution) {
1879 for (const int general_ct_index : var_to_constraints_[var]) {
1880 const int c = general_ct_index + linear_evaluator_.num_constraints();
1881 const int64_t v0 = constraints_[general_ct_index]->violation();
1882 constraints_[general_ct_index]->PerformMove(var, old_value, new_solution);
1883 const int64_t violation_delta =
1884 constraints_[general_ct_index]->violation() - v0;
1885 if (violation_delta != 0) {
1886 last_update_violation_changes_.push_back(c);
1887 }
1888 }
1889}
1890
1891void LsEvaluator::UpdateLinearScores(int var, int64_t old_value,
1892 int64_t new_value,
1893 absl::Span<const double> weights,
1894 absl::Span<const int64_t> jump_deltas,
1895 absl::Span<double> jump_scores) {
1896 DCHECK(RefIsPositive(var));
1897 if (old_value == new_value) return;
1898 last_update_violation_changes_.clear();
1899 linear_evaluator_.ClearAffectedVariables();
1900 linear_evaluator_.UpdateVariableAndScores(var, new_value - old_value, weights,
1901 jump_deltas, jump_scores,
1902 &last_update_violation_changes_);
1903}
1904
1906 // Maintain the list of violated constraints.
1907 dtime_ += 1e-8 * last_update_violation_changes_.size();
1908 for (const int c : last_update_violation_changes_) {
1910 }
1911}
1912
1914 int64_t evaluation = 0;
1915
1916 // Process the linear part.
1917 for (int i = 0; i < linear_evaluator_.num_constraints(); ++i) {
1918 evaluation += linear_evaluator_.Violation(i);
1919 DCHECK_GE(linear_evaluator_.Violation(i), 0);
1920 }
1921
1922 // Process the generic constraint part.
1923 for (const auto& ct : constraints_) {
1924 evaluation += ct->violation();
1925 DCHECK_GE(ct->violation(), 0);
1926 }
1927 return evaluation;
1928}
1929
1931 DCHECK(cp_model_.has_objective());
1932 return linear_evaluator_.Activity(/*c=*/0);
1933}
1934
1936 return linear_evaluator_.num_constraints();
1937}
1938
1940 return static_cast<int>(constraints_.size());
1941}
1942
1944 return linear_evaluator_.num_constraints() +
1945 static_cast<int>(constraints_.size());
1946}
1947
1949 int result = 0;
1950 for (int c = 0; c < linear_evaluator_.num_constraints(); ++c) {
1951 if (linear_evaluator_.Violation(c) > 0) {
1952 ++result;
1953 }
1954 }
1955 for (const auto& constraint : constraints_) {
1956 if (constraint->violation() > 0) {
1957 ++result;
1958 }
1959 }
1960 return result;
1961}
1962
1963int64_t LsEvaluator::Violation(int c) const {
1964 if (c < linear_evaluator_.num_constraints()) {
1965 return linear_evaluator_.Violation(c);
1966 } else {
1967 return constraints_[c - linear_evaluator_.num_constraints()]->violation();
1968 }
1969}
1970
1971bool LsEvaluator::IsViolated(int c) const {
1972 if (c < linear_evaluator_.num_constraints()) {
1973 return linear_evaluator_.IsViolated(c);
1974 } else {
1975 return constraints_[c - linear_evaluator_.num_constraints()]->violation() >
1976 0;
1977 }
1978}
1979
1980double LsEvaluator::WeightedViolation(absl::Span<const double> weights) const {
1981 DCHECK_EQ(weights.size(), NumEvaluatorConstraints());
1982 double result = linear_evaluator_.WeightedViolation(weights);
1983
1984 const int num_linear_constraints = linear_evaluator_.num_constraints();
1985 for (int c = 0; c < constraints_.size(); ++c) {
1986 result += static_cast<double>(constraints_[c]->violation()) *
1987 weights[num_linear_constraints + c];
1988 }
1989 return result;
1990}
1991
1993 bool linear_only, absl::Span<const double> weights, int var, int64_t delta,
1994 absl::Span<int64_t> mutable_solution) const {
1995 double result = linear_evaluator_.WeightedViolationDelta(weights, var, delta);
1996 if (linear_only) return result;
1997
1998 // We change the mutable solution here, and restore it after the evaluation.
1999 const int64_t old_value = mutable_solution[var];
2000 mutable_solution[var] += delta;
2001
2002 const int num_linear_constraints = linear_evaluator_.num_constraints();
2003 for (const int ct_index : var_to_constraints_[var]) {
2004 // We assume linear time delta computation in number of variables.
2005 // TODO(user): refine on a per constraint basis.
2006 dtime_ += 1e-8 * static_cast<double>(constraint_to_vars_[ct_index].size());
2007
2008 DCHECK_LT(ct_index, constraints_.size());
2009 const int64_t ct_delta = constraints_[ct_index]->ViolationDelta(
2010 var, old_value, mutable_solution);
2011 result += static_cast<double>(ct_delta) *
2012 weights[ct_index + num_linear_constraints];
2013 }
2014
2015 // Restore.
2016 mutable_solution[var] = old_value;
2017 return result;
2018}
2019
2021 int var) const {
2022 return jump_value_optimal_[var];
2023}
2024
2026 num_violated_constraint_per_var_ignoring_objective_.assign(
2027 cp_model_.variables_size(), 0);
2028 violated_constraints_.clear();
2029 const int num_constraints =
2031 for (int c = 0; c < num_constraints; ++c) {
2033 }
2034}
2035
2036void LsEvaluator::UpdateViolatedList(const int c) {
2037 if (Violation(c) > 0) {
2038 auto [it, inserted] = violated_constraints_.insert(c);
2039 // The constraint is violated. Add if needed.
2040 if (!inserted) return;
2041 if (IsObjectiveConstraint(c)) return;
2042 dtime_ += 1e-8 * ConstraintToVars(c).size();
2043 for (const int v : ConstraintToVars(c)) {
2044 num_violated_constraint_per_var_ignoring_objective_[v] += 1;
2045 }
2046 return;
2047 }
2048 if (violated_constraints_.erase(c) == 1) {
2049 if (IsObjectiveConstraint(c)) return;
2050 dtime_ += 1e-8 * ConstraintToVars(c).size();
2051 for (const int v : ConstraintToVars(c)) {
2052 num_violated_constraint_per_var_ignoring_objective_[v] -= 1;
2053 }
2054 }
2055}
2056
2057int64_t CompiledReservoirConstraint::BuildProfileAndReturnViolation(
2058 absl::Span<const int64_t> solution) {
2059 // Starts by filling the cache and profile_.
2060 capacity_value_ = ExprValue(capacity_, solution);
2061 const int num_events = time_values_.size();
2062 profile_.clear();
2063 for (int i = 0; i < num_events; ++i) {
2064 time_values_[i] = ExprValue(times_[i], solution);
2065 if (is_active_[i] != std::nullopt &&
2066 LiteralValue(*is_active_[i], solution) == 0) {
2067 demand_values_[i] = 0;
2068 } else {
2069 demand_values_[i] = ExprValue(demands_[i], solution);
2070 if (demand_values_[i] != 0) {
2071 profile_.push_back({time_values_[i], demand_values_[i]});
2072 }
2073 }
2074 }
2075
2076 if (profile_.empty()) return 0;
2077 absl::c_sort(profile_);
2078
2079 // Compress the profile for faster incremental evaluation.
2080 {
2081 int p = 0;
2082 for (int i = 1; i < profile_.size(); ++i) {
2083 if (profile_[i].time == profile_[p].time) {
2084 profile_[p].demand += profile_[i].demand;
2085 } else {
2086 profile_[++p] = profile_[i];
2087 }
2088 }
2089 profile_.resize(p + 1);
2090 }
2091
2092 int64_t overload = 0;
2093 int64_t current_load = 0;
2094 int64_t previous_time = std::numeric_limits<int64_t>::min();
2095 for (int i = 0; i < profile_.size(); ++i) {
2096 // At this point, current_load is the load at previous_time.
2097 const int64_t time = profile_[i].time;
2098 if (current_load > capacity_value_) {
2099 overload = CapAdd(overload, CapProd(current_load - capacity_value_,
2100 time - previous_time));
2101 }
2102
2103 current_load += profile_[i].demand;
2104 previous_time = time;
2105 }
2106 return overload;
2107}
2108
2109int64_t CompiledReservoirConstraint::IncrementalViolation(
2110 int var, absl::Span<const int64_t> solution) {
2111 const int64_t capacity = ExprValue(capacity_, solution);
2112 profile_delta_.clear();
2113 CHECK(RefIsPositive(var));
2114 for (const int i : dense_index_to_events_[var_to_dense_index_.at(var)]) {
2115 const int64_t time = ExprValue(times_[i], solution);
2116 int64_t demand = 0;
2117 if (is_active_[i] == std::nullopt ||
2118 LiteralValue(*is_active_[i], solution) == 1) {
2119 demand = ExprValue(demands_[i], solution);
2120 }
2121
2122 if (time == time_values_[i]) {
2123 if (demand != demand_values_[i]) {
2124 // Update the demand at time.
2125 profile_delta_.push_back({time, demand - demand_values_[i]});
2126 }
2127 } else {
2128 // Remove previous.
2129 if (demand_values_[i] != 0) {
2130 profile_delta_.push_back({time_values_[i], -demand_values_[i]});
2131 }
2132 // Add new.
2133 if (demand != 0) {
2134 profile_delta_.push_back({time, demand});
2135 }
2136 }
2137 }
2138
2139 // Abort early if there is no change.
2140 // This might happen because we use max(start + size, end) for the time and
2141 // even if some variable changed there, the time might not have.
2142 if (capacity == capacity_value_ && profile_delta_.empty()) {
2143 return violation_;
2144 }
2145 absl::c_sort(profile_delta_);
2146
2147 // Similar algo, but we scan the two vectors at once.
2148 int64_t overload = 0;
2149 int64_t current_load = 0;
2150 int64_t previous_time = std::numeric_limits<int64_t>::min();
2151
2152 // TODO(user): This code is the hotspot for our local search on cumulative.
2153 // It can probably be slighlty improved. We might also be able to abort early
2154 // if we know that capacity is high enough compared to the highest point of
2155 // the profile.
2156 int i = 0;
2157 int j = 0;
2158 const absl::Span<const Event> i_profile(profile_);
2159 const absl::Span<const Event> j_profile(profile_delta_);
2160 while (true) {
2161 int64_t time;
2162 if (i < i_profile.size() && j < j_profile.size()) {
2163 time = std::min(i_profile[i].time, j_profile[j].time);
2164 } else if (i < i_profile.size()) {
2165 time = i_profile[i].time;
2166 } else if (j < j_profile.size()) {
2167 time = j_profile[j].time;
2168 } else {
2169 // End of loop.
2170 break;
2171 }
2172
2173 // Update overload if needed.
2174 // At this point, current_load is the load at previous_time.
2175 if (current_load > capacity) {
2176 overload = CapAdd(overload,
2177 CapProd(current_load - capacity, time - previous_time));
2178 }
2179
2180 // Update i and current load.
2181 while (i < i_profile.size() && i_profile[i].time == time) {
2182 current_load += i_profile[i].demand;
2183 i++;
2184 }
2185
2186 // Update j and current load.
2187 while (j < j_profile.size() && j_profile[j].time == time) {
2188 current_load += j_profile[j].demand;
2189 j++;
2190 }
2191
2192 previous_time = time;
2193 }
2194 return overload;
2195}
2196
2197void CompiledReservoirConstraint::AppendVariablesForEvent(
2198 int i, std::vector<int>* result) const {
2199 if (is_active_[i] != std::nullopt) {
2200 result->push_back(PositiveRef(*is_active_[i]));
2201 }
2202 for (const int var : times_[i].vars()) {
2203 result->push_back(PositiveRef(var));
2204 }
2205 for (const int var : demands_[i].vars()) {
2206 result->push_back(PositiveRef(var));
2207 }
2208}
2209
2210void CompiledReservoirConstraint::InitializeDenseIndexToEvents() {
2211 // We scan the constraint a few times, but this is called once, so we don't
2212 // care too much.
2213 CpModelProto unused;
2214 int num_dense_indices = 0;
2215 for (const int var : UsedVariables(unused)) {
2216 var_to_dense_index_[var] = num_dense_indices++;
2217 }
2218
2219 CompactVectorVector<int, int> event_to_dense_indices;
2220 event_to_dense_indices.reserve(times_.size());
2221 const int num_events = times_.size();
2222 std::vector<int> result;
2223 for (int i = 0; i < num_events; ++i) {
2224 result.clear();
2225 AppendVariablesForEvent(i, &result);
2226
2227 // Remap and add.
2228 for (int& var : result) {
2229 var = var_to_dense_index_.at(var);
2230 }
2232 event_to_dense_indices.Add(result);
2233 }
2234
2235 // Note that because of the capacity (which might be variable) it is important
2236 // to resize this to num_dense_indices.
2237 dense_index_to_events_.ResetFromTranspose(event_to_dense_indices,
2238 num_dense_indices);
2239}
2240
2242 const CpModelProto& /*model_proto*/) const {
2243 std::vector<int> result;
2244 const int num_events = times_.size();
2245 for (int i = 0; i < num_events; ++i) {
2246 AppendVariablesForEvent(i, &result);
2247 }
2248 for (const int var : capacity_.vars()) {
2249 result.push_back(PositiveRef(var));
2250 }
2252 result.shrink_to_fit();
2253 return result;
2254}
2255
2256} // namespace sat
2257} // namespace operations_research
std::pair< iterator, bool > insert(T value)
Definition dense_set.h:56
std::vector< int64_t > FlattenedIntervals() const
bool Contains(int64_t value) const
int64_t Distance(int64_t value) const
static IntegralType CeilOfRatio(IntegralType numerator, IntegralType denominator)
Definition mathutil.h:40
static IntegralType FloorOfRatio(IntegralType numerator, IntegralType denominator)
Definition mathutil.h:54
int64_t ComputeViolation(absl::Span< const int64_t > solution) override
CompiledAllDiffConstraint(const ConstraintProto &ct_proto)
--— CompiledAllDiffConstraint --—
CompiledBoolXorConstraint(const ConstraintProto &ct_proto)
--— CompiledBoolXorConstraint --—
int64_t ViolationDelta(int, int64_t, absl::Span< const int64_t > solution_with_new_value) override
Returns the delta if var changes from old_value to solution[var].
int64_t ComputeViolation(absl::Span< const int64_t > solution) override
void PerformMove(int var, int64_t old_value, absl::Span< const int64_t > new_solution) override
Updates the violation with the new value.
int64_t ComputeViolation(absl::Span< const int64_t > solution) override
int64_t ViolationDelta(int var, int64_t old_value, absl::Span< const int64_t > solution_with_new_value) override
Returns the delta if var changes from old_value to solution[var].
CompiledConstraintWithProto(const ConstraintProto &ct_proto)
--— CompiledConstraintWithProto --—
std::vector< int > UsedVariables(const CpModelProto &model_proto) const final
This just returns the variables used by the stored ct_proto_.
int64_t violation() const
The cached violation of this constraint.
virtual void PerformMove(int var, int64_t old_value, absl::Span< const int64_t > solution_with_new_value)
Updates the violation with the new value.
virtual int64_t ComputeViolation(absl::Span< const int64_t > solution)=0
void InitializeViolation(absl::Span< const int64_t > solution)
Recomputes the violation of the constraint from scratch.
virtual int64_t ViolationDelta(int var, int64_t old_value, absl::Span< const int64_t > solution_with_new_value)
Returns the delta if var changes from old_value to solution[var].
CompiledIntDivConstraint(const ConstraintProto &ct_proto)
--— CompiledIntDivConstraint --—
int64_t ComputeViolation(absl::Span< const int64_t > solution) override
int64_t ComputeViolation(absl::Span< const int64_t > solution) override
CompiledIntModConstraint(const ConstraintProto &ct_proto)
--— CompiledIntModConstraint --—
CompiledIntProdConstraint(const ConstraintProto &ct_proto)
--— CompiledIntProdConstraint --—
int64_t ComputeViolation(absl::Span< const int64_t > solution) override
int64_t ComputeViolation(absl::Span< const int64_t > solution) override
CompiledLinMaxConstraint(const ConstraintProto &ct_proto)
--— CompiledLinMaxConstraint --—
int64_t ComputeViolation(absl::Span< const int64_t > solution) override
CompiledNoOverlap2dConstraint(const ConstraintProto &ct_proto, const CpModelProto &cp_model)
std::vector< int > UsedVariables(const CpModelProto &model_proto) const final
void UpdateScoreOnWeightUpdate(int c, absl::Span< const int64_t > jump_deltas, absl::Span< double > var_to_score_change)
Also for feasibility jump.
double WeightedViolationDelta(absl::Span< const double > weights, int var, int64_t delta) const
void ComputeInitialActivities(absl::Span< const int64_t > solution)
Compute activities.
void AddLiteral(int ct_index, int lit, int64_t coeff=1)
int NewConstraint(Domain domain)
Returns the index of the new constraint.
double WeightedViolation(absl::Span< const double > weights) const
bool VarIsConsistent(int var) const
Used to DCHECK the state of the evaluator.
void UpdateVariableAndScores(int var, int64_t delta, absl::Span< const double > weights, absl::Span< const int64_t > jump_deltas, absl::Span< double > jump_scores, std::vector< int > *constraints_with_changed_violations)
void AddLinearExpression(int ct_index, const LinearExpressionProto &expr, int64_t multiplier)
void PrecomputeCompactView(absl::Span< const int64_t > var_max_variation)
bool ViolationChangeIsConvex(int var) const
Checks if the jump value of a variable is always optimal.
std::vector< int64_t > SlopeBreakpoints(int var, int64_t current_value, const Domain &var_domain) const
void AddTerm(int ct_index, int var, int64_t coeff, int64_t offset=0)
double WeightedViolation(absl::Span< const double > weights) const
void ComputeAllNonLinearViolations(absl::Span< const int64_t > solution)
void UpdateLinearScores(int var, int64_t old_value, int64_t new_value, absl::Span< const double > weights, absl::Span< const int64_t > jump_deltas, absl::Span< double > jump_scores)
Function specific to the linear only feasibility jump.
bool VariableOnlyInLinearConstraintWithConvexViolationChange(int var) const
Indicates if the computed jump value is always the best choice.
void ComputeAllViolations(absl::Span< const int64_t > solution)
Recomputes the violations of all constraints (resp only non-linear one).
void UpdateNonLinearViolations(int var, int64_t old_value, absl::Span< const int64_t > new_solution)
Recomputes the violations of all impacted non linear constraints.
bool ReduceObjectiveBounds(int64_t lb, int64_t ub)
int64_t SumOfViolations()
Simple summation metric for the constraint and objective violations.
int64_t ObjectiveActivity() const
Returns the objective activity in the current state.
absl::Span< const int > ConstraintToVars(int c) const
LsEvaluator(const CpModelProto &cp_model, const SatParameters &params)
The cp_model must outlive this class.
double WeightedViolationDelta(bool linear_only, absl::Span< const double > weights, int var, int64_t delta, absl::Span< int64_t > mutable_solution) const
std::vector< int > UsedVariables(const CpModelProto &model_proto) const final
NoOverlapBetweenTwoIntervals(int interval_0, int interval_1, const CpModelProto &cp_model)
--— NoOverlapBetweenTwoIntervals --—
void STLSortAndRemoveDuplicates(T *v, const LessFunc &less_func)
Definition stl_util.h:58
void STLClearObject(T *obj)
Definition stl_util.h:123
void AddCircuitFlowConstraints(LinearIncrementalEvaluator &linear_evaluator, const ConstraintProto &ct_proto)
int64_t OverlapOfTwoIntervals(const ConstraintProto &interval1, const ConstraintProto &interval2, absl::Span< const int64_t > solution)
--— CompiledNoOverlap2dConstraint --—
std::vector< int > UsedVariables(const ConstraintProto &ct)
std::vector< int > UsedIntervals(const ConstraintProto &ct)
Returns the sorted list of interval used by a constraint.
int64_t NoOverlapMinRepairDistance(const ConstraintProto &interval1, const ConstraintProto &interval2, absl::Span< const int64_t > solution)
Domain ReadDomainFromProto(const ProtoWithDomain &proto)
Reads a Domain from the domain field of a proto.
int NegatedRef(int ref)
Small utility functions to deal with negative variable/literal references.
In SWIG mode, we don't want anything besides these top-level includes.
int64_t CapAdd(int64_t x, int64_t y)
Select next search node to expand Select next item_i to add this new search node to the search Generate a new search node where item_i is not in the knapsack Check validity of this new partial solution(using propagators) - If valid
int64_t CapSub(int64_t x, int64_t y)
int64_t CapProd(int64_t x, int64_t y)
int64_t Max() const
Returns the max of the domain.
Definition model.cc:346
int64_t Min() const
Returns the min of the domain.
Definition model.cc:340