1 // Copyright 2010-2021 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 
14 #include "ortools/sat/presolve_context.h"
15 
16 #include <algorithm>
17 #include <cstdint>
18 #include <limits>
19 #include <string>
20 
21 #include "ortools/base/map_util.h"
22 #include "ortools/base/mathutil.h"
23 #include "ortools/port/proto_utils.h"
24 #include "ortools/sat/cp_model.pb.h"
25 #include "ortools/sat/cp_model_loader.h"
26 #include "ortools/sat/lp_utils.h"
27 #include "ortools/util/saturated_arithmetic.h"
28 
29 namespace operations_research {
30 namespace sat {
31 
Get(PresolveContext * context) const32 int SavedLiteral::Get(PresolveContext* context) const {
33   return context->GetLiteralRepresentative(ref_);
34 }
35 
Get(PresolveContext * context) const36 int SavedVariable::Get(PresolveContext* context) const {
37   return context->GetVariableRepresentative(ref_);
38 }
39 
ClearStats()40 void PresolveContext::ClearStats() { stats_by_rule_name_.clear(); }
41 
NewIntVar(const Domain & domain)42 int PresolveContext::NewIntVar(const Domain& domain) {
43   IntegerVariableProto* const var = working_model->add_variables();
44   FillDomainInProto(domain, var);
45   InitializeNewDomains();
46   return working_model->variables_size() - 1;
47 }
48 
NewBoolVar()49 int PresolveContext::NewBoolVar() { return NewIntVar(Domain(0, 1)); }
50 
GetOrCreateConstantVar(int64_t cst)51 int PresolveContext::GetOrCreateConstantVar(int64_t cst) {
52   if (!constant_to_ref_.contains(cst)) {
53     constant_to_ref_[cst] = SavedVariable(working_model->variables_size());
54     IntegerVariableProto* const var_proto = working_model->add_variables();
55     var_proto->add_domain(cst);
56     var_proto->add_domain(cst);
57     InitializeNewDomains();
58   }
59   return constant_to_ref_[cst].Get(this);
60 }
61 
62 // a => b.
AddImplication(int a,int b)63 void PresolveContext::AddImplication(int a, int b) {
64   ConstraintProto* const ct = working_model->add_constraints();
65   ct->add_enforcement_literal(a);
66   ct->mutable_bool_and()->add_literals(b);
67 }
68 
69 // b => x in [lb, ub].
AddImplyInDomain(int b,int x,const Domain & domain)70 void PresolveContext::AddImplyInDomain(int b, int x, const Domain& domain) {
71   ConstraintProto* const imply = working_model->add_constraints();
72 
73   // Doing it like this seems to use slightly less memory.
74   // TODO(user): Find the best way to create such small proto.
75   imply->mutable_enforcement_literal()->Resize(1, b);
76   LinearConstraintProto* mutable_linear = imply->mutable_linear();
77   mutable_linear->mutable_vars()->Resize(1, x);
78   mutable_linear->mutable_coeffs()->Resize(1, 1);
79   FillDomainInProto(domain, mutable_linear);
80 }
81 
DomainIsEmpty(int ref) const82 bool PresolveContext::DomainIsEmpty(int ref) const {
83   return domains[PositiveRef(ref)].IsEmpty();
84 }
85 
IsFixed(int ref) const86 bool PresolveContext::IsFixed(int ref) const {
87   DCHECK_LT(PositiveRef(ref), domains.size());
88   DCHECK(!DomainIsEmpty(ref));
89   return domains[PositiveRef(ref)].IsFixed();
90 }
91 
CanBeUsedAsLiteral(int ref) const92 bool PresolveContext::CanBeUsedAsLiteral(int ref) const {
93   const int var = PositiveRef(ref);
94   return domains[var].Min() >= 0 && domains[var].Max() <= 1;
95 }
96 
LiteralIsTrue(int lit) const97 bool PresolveContext::LiteralIsTrue(int lit) const {
98   DCHECK(CanBeUsedAsLiteral(lit));
99   if (RefIsPositive(lit)) {
100     return domains[lit].Min() == 1;
101   } else {
102     return domains[PositiveRef(lit)].Max() == 0;
103   }
104 }
105 
LiteralIsFalse(int lit) const106 bool PresolveContext::LiteralIsFalse(int lit) const {
107   DCHECK(CanBeUsedAsLiteral(lit));
108   if (RefIsPositive(lit)) {
109     return domains[lit].Max() == 0;
110   } else {
111     return domains[PositiveRef(lit)].Min() == 1;
112   }
113 }
114 
MinOf(int ref) const115 int64_t PresolveContext::MinOf(int ref) const {
116   DCHECK(!DomainIsEmpty(ref));
117   return RefIsPositive(ref) ? domains[PositiveRef(ref)].Min()
118                             : -domains[PositiveRef(ref)].Max();
119 }
120 
MaxOf(int ref) const121 int64_t PresolveContext::MaxOf(int ref) const {
122   DCHECK(!DomainIsEmpty(ref));
123   return RefIsPositive(ref) ? domains[PositiveRef(ref)].Max()
124                             : -domains[PositiveRef(ref)].Min();
125 }
126 
FixedValue(int ref) const127 int64_t PresolveContext::FixedValue(int ref) const {
128   DCHECK(!DomainIsEmpty(ref));
129   CHECK(IsFixed(ref));
130   return RefIsPositive(ref) ? domains[PositiveRef(ref)].Min()
131                             : -domains[PositiveRef(ref)].Min();
132 }
133 
MinOf(const LinearExpressionProto & expr) const134 int64_t PresolveContext::MinOf(const LinearExpressionProto& expr) const {
135   int64_t result = expr.offset();
136   for (int i = 0; i < expr.vars_size(); ++i) {
137     const int64_t coeff = expr.coeffs(i);
138     if (coeff > 0) {
139       result += coeff * MinOf(expr.vars(i));
140     } else {
141       result += coeff * MaxOf(expr.vars(i));
142     }
143   }
144   return result;
145 }
146 
MaxOf(const LinearExpressionProto & expr) const147 int64_t PresolveContext::MaxOf(const LinearExpressionProto& expr) const {
148   int64_t result = expr.offset();
149   for (int i = 0; i < expr.vars_size(); ++i) {
150     const int64_t coeff = expr.coeffs(i);
151     if (coeff > 0) {
152       result += coeff * MaxOf(expr.vars(i));
153     } else {
154       result += coeff * MinOf(expr.vars(i));
155     }
156   }
157   return result;
158 }
159 
IsFixed(const LinearExpressionProto & expr) const160 bool PresolveContext::IsFixed(const LinearExpressionProto& expr) const {
161   for (int i = 0; i < expr.vars_size(); ++i) {
162     if (expr.coeffs(i) != 0 && !IsFixed(expr.vars(i))) return false;
163   }
164   return true;
165 }
166 
FixedValue(const LinearExpressionProto & expr) const167 int64_t PresolveContext::FixedValue(const LinearExpressionProto& expr) const {
168   int64_t result = expr.offset();
169   for (int i = 0; i < expr.vars_size(); ++i) {
170     if (expr.coeffs(i) == 0) continue;
171     result += expr.coeffs(i) * FixedValue(expr.vars(i));
172   }
173   return result;
174 }
175 
DomainSuperSetOf(const LinearExpressionProto & expr) const176 Domain PresolveContext::DomainSuperSetOf(
177     const LinearExpressionProto& expr) const {
178   Domain result(expr.offset());
179   for (int i = 0; i < expr.vars_size(); ++i) {
180     result = result.AdditionWith(
181         DomainOf(expr.vars(i)).MultiplicationBy(expr.coeffs(i)));
182   }
183   return result;
184 }
185 
ExpressionIsAffineBoolean(const LinearExpressionProto & expr) const186 bool PresolveContext::ExpressionIsAffineBoolean(
187     const LinearExpressionProto& expr) const {
188   if (expr.vars().size() != 1) return false;
189   return CanBeUsedAsLiteral(expr.vars(0));
190 }
191 
LiteralForExpressionMax(const LinearExpressionProto & expr) const192 int PresolveContext::LiteralForExpressionMax(
193     const LinearExpressionProto& expr) const {
194   const int ref = expr.vars(0);
195   return RefIsPositive(ref) == (expr.coeffs(0) > 0) ? ref : NegatedRef(ref);
196 }
197 
ExpressionIsSingleVariable(const LinearExpressionProto & expr) const198 bool PresolveContext::ExpressionIsSingleVariable(
199     const LinearExpressionProto& expr) const {
200   return expr.offset() == 0 && expr.vars_size() == 1 && expr.coeffs(0) == 1;
201 }
202 
ExpressionIsALiteral(const LinearExpressionProto & expr,int * literal) const203 bool PresolveContext::ExpressionIsALiteral(const LinearExpressionProto& expr,
204                                            int* literal) const {
205   if (expr.vars_size() != 1) return false;
206   const int ref = expr.vars(0);
207   const int var = PositiveRef(ref);
208   if (MinOf(var) < 0 || MaxOf(var) > 1) return false;
209 
210   if (expr.offset() == 0 && expr.coeffs(0) == 1 && RefIsPositive(ref)) {
211     if (literal != nullptr) *literal = ref;
212     return true;
213   }
214   if (expr.offset() == 1 && expr.coeffs(0) == -1 && RefIsPositive(ref)) {
215     if (literal != nullptr) *literal = NegatedRef(ref);
216     return true;
217   }
218   if (expr.offset() == 1 && expr.coeffs(0) == 1 && !RefIsPositive(ref)) {
219     if (literal != nullptr) *literal = ref;
220     return true;
221   }
222   return false;
223 }
224 
225 // Note that we only support converted intervals.
IntervalIsConstant(int ct_ref) const226 bool PresolveContext::IntervalIsConstant(int ct_ref) const {
227   const ConstraintProto& proto = working_model->constraints(ct_ref);
228   if (!proto.enforcement_literal().empty()) return false;
229   if (!proto.interval().has_start()) return false;
230   for (const int var : proto.interval().start().vars()) {
231     if (!IsFixed(var)) return false;
232   }
233   for (const int var : proto.interval().size().vars()) {
234     if (!IsFixed(var)) return false;
235   }
236   for (const int var : proto.interval().end().vars()) {
237     if (!IsFixed(var)) return false;
238   }
239   return true;
240 }
241 
IntervalDebugString(int ct_ref) const242 std::string PresolveContext::IntervalDebugString(int ct_ref) const {
243   if (IntervalIsConstant(ct_ref)) {
244     return absl::StrCat("interval_", ct_ref, "(", StartMin(ct_ref), "..",
245                         EndMax(ct_ref), ")");
246   } else if (ConstraintIsOptional(ct_ref)) {
247     const int literal =
248         working_model->constraints(ct_ref).enforcement_literal(0);
249     if (SizeMin(ct_ref) == SizeMax(ct_ref)) {
250       return absl::StrCat("interval_", ct_ref, "(lit=", literal, ", ",
251                           StartMin(ct_ref), " --(", SizeMin(ct_ref), ")--> ",
252                           EndMax(ct_ref), ")");
253     } else {
254       return absl::StrCat("interval_", ct_ref, "(lit=", literal, ", ",
255                           StartMin(ct_ref), " --(", SizeMin(ct_ref), "..",
256                           SizeMax(ct_ref), ")--> ", EndMax(ct_ref), ")");
257     }
258   } else if (SizeMin(ct_ref) == SizeMax(ct_ref)) {
259     return absl::StrCat("interval_", ct_ref, "(", StartMin(ct_ref), " --(",
260                         SizeMin(ct_ref), ")--> ", EndMax(ct_ref), ")");
261   } else {
262     return absl::StrCat("interval_", ct_ref, "(", StartMin(ct_ref), " --(",
263                         SizeMin(ct_ref), "..", SizeMax(ct_ref), ")--> ",
264                         EndMax(ct_ref), ")");
265   }
266 }
267 
StartMin(int ct_ref) const268 int64_t PresolveContext::StartMin(int ct_ref) const {
269   const IntervalConstraintProto& interval =
270       working_model->constraints(ct_ref).interval();
271   return MinOf(interval.start());
272 }
273 
StartMax(int ct_ref) const274 int64_t PresolveContext::StartMax(int ct_ref) const {
275   const IntervalConstraintProto& interval =
276       working_model->constraints(ct_ref).interval();
277   return MaxOf(interval.start());
278 }
279 
EndMin(int ct_ref) const280 int64_t PresolveContext::EndMin(int ct_ref) const {
281   const IntervalConstraintProto& interval =
282       working_model->constraints(ct_ref).interval();
283   return MinOf(interval.end());
284 }
285 
EndMax(int ct_ref) const286 int64_t PresolveContext::EndMax(int ct_ref) const {
287   const IntervalConstraintProto& interval =
288       working_model->constraints(ct_ref).interval();
289   return MaxOf(interval.end());
290 }
291 
SizeMin(int ct_ref) const292 int64_t PresolveContext::SizeMin(int ct_ref) const {
293   const IntervalConstraintProto& interval =
294       working_model->constraints(ct_ref).interval();
295   return MinOf(interval.size());
296 }
297 
SizeMax(int ct_ref) const298 int64_t PresolveContext::SizeMax(int ct_ref) const {
299   const IntervalConstraintProto& interval =
300       working_model->constraints(ct_ref).interval();
301   return MaxOf(interval.size());
302 }
303 
304 // Important: To be sure a variable can be removed, we need it to not be a
305 // representative of both affine and equivalence relation.
VariableIsNotRepresentativeOfEquivalenceClass(int var) const306 bool PresolveContext::VariableIsNotRepresentativeOfEquivalenceClass(
307     int var) const {
308   DCHECK(RefIsPositive(var));
309   if (affine_relations_.ClassSize(var) > 1 &&
310       affine_relations_.Get(var).representative == var) {
311     return false;
312   }
313   if (var_equiv_relations_.ClassSize(var) > 1 &&
314       var_equiv_relations_.Get(var).representative == var) {
315     return false;
316   }
317   return true;
318 }
319 
VariableIsRemovable(int ref) const320 bool PresolveContext::VariableIsRemovable(int ref) const {
321   const int var = PositiveRef(ref);
322   return VariableIsNotRepresentativeOfEquivalenceClass(var) &&
323          !keep_all_feasible_solutions;
324 }
325 
326 // Tricky: If this variable is equivalent to another one (but not the
327 // representative) and appear in just one constraint, then this constraint must
328 // be the affine defining one. And in this case the code using this function
329 // should do the proper stuff.
VariableIsUniqueAndRemovable(int ref) const330 bool PresolveContext::VariableIsUniqueAndRemovable(int ref) const {
331   if (!ConstraintVariableGraphIsUpToDate()) return false;
332   const int var = PositiveRef(ref);
333   return var_to_constraints_[var].size() == 1 && VariableIsRemovable(var);
334 }
335 
VariableWithCostIsUnique(int ref) const336 bool PresolveContext::VariableWithCostIsUnique(int ref) const {
337   if (!ConstraintVariableGraphIsUpToDate()) return false;
338   const int var = PositiveRef(ref);
339   return VariableIsNotRepresentativeOfEquivalenceClass(var) &&
340          var_to_constraints_[var].contains(kObjectiveConstraint) &&
341          var_to_constraints_[var].size() == 2;
342 }
343 
344 // Tricky: Same remark as for VariableIsUniqueAndRemovable().
345 //
346 // Also if the objective domain is constraining, we can't have a preferred
347 // direction, so we cannot easily remove such variable.
VariableWithCostIsUniqueAndRemovable(int ref) const348 bool PresolveContext::VariableWithCostIsUniqueAndRemovable(int ref) const {
349   if (!ConstraintVariableGraphIsUpToDate()) return false;
350   const int var = PositiveRef(ref);
351   return VariableIsRemovable(var) && !objective_domain_is_constraining_ &&
352          VariableWithCostIsUnique(var);
353 }
354 
355 // Here, even if the variable is equivalent to others, if its affine defining
356 // constraints where removed, then it is not needed anymore.
VariableIsNotUsedAnymore(int ref) const357 bool PresolveContext::VariableIsNotUsedAnymore(int ref) const {
358   if (!ConstraintVariableGraphIsUpToDate()) return false;
359   return var_to_constraints_[PositiveRef(ref)].empty();
360 }
361 
MarkVariableAsRemoved(int ref)362 void PresolveContext::MarkVariableAsRemoved(int ref) {
363   removed_variables_.insert(PositiveRef(ref));
364 }
365 
366 // Note(user): I added an indirection and a function for this to be able to
367 // display debug information when this return false. This should actually never
368 // return false in the cases where it is used.
VariableWasRemoved(int ref) const369 bool PresolveContext::VariableWasRemoved(int ref) const {
370   // It is okay to reuse removed fixed variable.
371   if (IsFixed(ref)) return false;
372   if (!removed_variables_.contains(PositiveRef(ref))) return false;
373   if (!var_to_constraints_[PositiveRef(ref)].empty()) {
374     SOLVER_LOG(logger_, "Variable ", PositiveRef(ref),
375                " was removed, yet it appears in some constraints!");
376     SOLVER_LOG(logger_, "affine relation: ",
377                AffineRelationDebugString(PositiveRef(ref)));
378     for (const int c : var_to_constraints_[PositiveRef(ref)]) {
379       SOLVER_LOG(
380           logger_, "constraint #", c, " : ",
381           c >= 0 ? working_model->constraints(c).ShortDebugString() : "");
382     }
383   }
384   return true;
385 }
386 
VariableIsOnlyUsedInEncodingAndMaybeInObjective(int ref) const387 bool PresolveContext::VariableIsOnlyUsedInEncodingAndMaybeInObjective(
388     int ref) const {
389   if (!ConstraintVariableGraphIsUpToDate()) return false;
390   const int var = PositiveRef(ref);
391   return var_to_num_linear1_[var] == var_to_constraints_[var].size() ||
392          (var_to_constraints_[var].contains(kObjectiveConstraint) &&
393           var_to_num_linear1_[var] + 1 == var_to_constraints_[var].size());
394 }
395 
DomainOf(int ref) const396 Domain PresolveContext::DomainOf(int ref) const {
397   Domain result;
398   if (RefIsPositive(ref)) {
399     result = domains[ref];
400   } else {
401     result = domains[PositiveRef(ref)].Negation();
402   }
403   return result;
404 }
405 
DomainContains(int ref,int64_t value) const406 bool PresolveContext::DomainContains(int ref, int64_t value) const {
407   if (!RefIsPositive(ref)) {
408     return domains[PositiveRef(ref)].Contains(-value);
409   }
410   return domains[ref].Contains(value);
411 }
412 
DomainContains(const LinearExpressionProto & expr,int64_t value) const413 bool PresolveContext::DomainContains(const LinearExpressionProto& expr,
414                                      int64_t value) const {
415   CHECK_LE(expr.vars_size(), 1);
416   if (IsFixed(expr)) {
417     return FixedValue(expr) == value;
418   }
419   if ((value - expr.offset()) % expr.coeffs(0) != 0) return false;
420   return DomainContains(expr.vars(0), (value - expr.offset()) / expr.coeffs(0));
421 }
422 
IntersectDomainWith(int ref,const Domain & domain,bool * domain_modified)423 ABSL_MUST_USE_RESULT bool PresolveContext::IntersectDomainWith(
424     int ref, const Domain& domain, bool* domain_modified) {
425   DCHECK(!DomainIsEmpty(ref));
426   const int var = PositiveRef(ref);
427 
428   if (RefIsPositive(ref)) {
429     if (domains[var].IsIncludedIn(domain)) {
430       return true;
431     }
432     domains[var] = domains[var].IntersectionWith(domain);
433   } else {
434     const Domain temp = domain.Negation();
435     if (domains[var].IsIncludedIn(temp)) {
436       return true;
437     }
438     domains[var] = domains[var].IntersectionWith(temp);
439   }
440 
441   if (domain_modified != nullptr) {
442     *domain_modified = true;
443   }
444   modified_domains.Set(var);
445   if (domains[var].IsEmpty()) {
446     is_unsat_ = true;
447     return false;
448   }
449 
450   // Propagate the domain of the representative right away.
451   // Note that the recursive call should only by one level deep.
452   const AffineRelation::Relation r = GetAffineRelation(var);
453   if (r.representative == var) return true;
454   return IntersectDomainWith(r.representative,
455                              DomainOf(var)
456                                  .AdditionWith(Domain(-r.offset))
457                                  .InverseMultiplicationBy(r.coeff));
458 }
459 
IntersectDomainWith(const LinearExpressionProto & expr,const Domain & domain,bool * domain_modified)460 ABSL_MUST_USE_RESULT bool PresolveContext::IntersectDomainWith(
461     const LinearExpressionProto& expr, const Domain& domain,
462     bool* domain_modified) {
463   if (expr.vars().empty()) {
464     if (domain.Contains(expr.offset())) {
465       return true;
466     } else {
467       is_unsat_ = true;
468       return false;
469     }
470   }
471   if (expr.vars().size() == 1) {  // Affine
472     return IntersectDomainWith(expr.vars(0),
473                                domain.AdditionWith(Domain(-expr.offset()))
474                                    .InverseMultiplicationBy(expr.coeffs(0)),
475                                domain_modified);
476   }
477 
478   // We don't do anything for longer expression for now.
479   return true;
480 }
481 
SetLiteralToFalse(int lit)482 ABSL_MUST_USE_RESULT bool PresolveContext::SetLiteralToFalse(int lit) {
483   const int var = PositiveRef(lit);
484   const int64_t value = RefIsPositive(lit) ? 0 : 1;
485   return IntersectDomainWith(var, Domain(value));
486 }
487 
SetLiteralToTrue(int lit)488 ABSL_MUST_USE_RESULT bool PresolveContext::SetLiteralToTrue(int lit) {
489   return SetLiteralToFalse(NegatedRef(lit));
490 }
491 
ConstraintIsInactive(int index) const492 bool PresolveContext::ConstraintIsInactive(int index) const {
493   const ConstraintProto& ct = working_model->constraints(index);
494   if (ct.constraint_case() ==
495       ConstraintProto::ConstraintCase::CONSTRAINT_NOT_SET) {
496     return true;
497   }
498   for (const int literal : ct.enforcement_literal()) {
499     if (LiteralIsFalse(literal)) return true;
500   }
501   return false;
502 }
503 
ConstraintIsOptional(int ct_ref) const504 bool PresolveContext::ConstraintIsOptional(int ct_ref) const {
505   const ConstraintProto& ct = working_model->constraints(ct_ref);
506   bool contains_one_free_literal = false;
507   for (const int literal : ct.enforcement_literal()) {
508     if (LiteralIsFalse(literal)) return false;
509     if (!LiteralIsTrue(literal)) contains_one_free_literal = true;
510   }
511   return contains_one_free_literal;
512 }
513 
UpdateRuleStats(const std::string & name,int num_times)514 void PresolveContext::UpdateRuleStats(const std::string& name, int num_times) {
515   // We only count if we are going to display it.
516   if (logger_->LoggingIsEnabled()) {
517     VLOG(2) << num_presolve_operations << " : " << name;
518     stats_by_rule_name_[name] += num_times;
519   }
520   num_presolve_operations += num_times;
521 }
522 
UpdateLinear1Usage(const ConstraintProto & ct,int c)523 void PresolveContext::UpdateLinear1Usage(const ConstraintProto& ct, int c) {
524   const int old_var = constraint_to_linear1_var_[c];
525   if (old_var >= 0) {
526     var_to_num_linear1_[old_var]--;
527   }
528   if (ct.constraint_case() == ConstraintProto::ConstraintCase::kLinear &&
529       ct.linear().vars().size() == 1) {
530     const int var = PositiveRef(ct.linear().vars(0));
531     constraint_to_linear1_var_[c] = var;
532     var_to_num_linear1_[var]++;
533   }
534 }
535 
AddVariableUsage(int c)536 void PresolveContext::AddVariableUsage(int c) {
537   const ConstraintProto& ct = working_model->constraints(c);
538   constraint_to_vars_[c] = UsedVariables(ct);
539   constraint_to_intervals_[c] = UsedIntervals(ct);
540   for (const int v : constraint_to_vars_[c]) {
541     DCHECK(!VariableWasRemoved(v));
542     var_to_constraints_[v].insert(c);
543   }
544   for (const int i : constraint_to_intervals_[c]) interval_usage_[i]++;
545   UpdateLinear1Usage(ct, c);
546 }
547 
UpdateConstraintVariableUsage(int c)548 void PresolveContext::UpdateConstraintVariableUsage(int c) {
549   if (is_unsat_) return;
550   DCHECK_EQ(constraint_to_vars_.size(), working_model->constraints_size());
551   const ConstraintProto& ct = working_model->constraints(c);
552 
553   // We don't optimize the interval usage as this is not super frequent.
554   for (const int i : constraint_to_intervals_[c]) interval_usage_[i]--;
555   constraint_to_intervals_[c] = UsedIntervals(ct);
556   for (const int i : constraint_to_intervals_[c]) interval_usage_[i]++;
557 
558   // For the variables, we avoid an erase() followed by an insert() for the
559   // variables that didn't change.
560   tmp_new_usage_ = UsedVariables(ct);
561   const std::vector<int>& old_usage = constraint_to_vars_[c];
562   const int old_size = old_usage.size();
563   int i = 0;
564   for (const int var : tmp_new_usage_) {
565     DCHECK(!VariableWasRemoved(var));
566     while (i < old_size && old_usage[i] < var) {
567       var_to_constraints_[old_usage[i]].erase(c);
568       ++i;
569     }
570     if (i < old_size && old_usage[i] == var) {
571       ++i;
572     } else {
573       var_to_constraints_[var].insert(c);
574     }
575   }
576   for (; i < old_size; ++i) var_to_constraints_[old_usage[i]].erase(c);
577   constraint_to_vars_[c] = tmp_new_usage_;
578 
579   UpdateLinear1Usage(ct, c);
580 }
581 
ConstraintVariableGraphIsUpToDate() const582 bool PresolveContext::ConstraintVariableGraphIsUpToDate() const {
583   return constraint_to_vars_.size() == working_model->constraints_size();
584 }
585 
UpdateNewConstraintsVariableUsage()586 void PresolveContext::UpdateNewConstraintsVariableUsage() {
587   if (is_unsat_) return;
588   const int old_size = constraint_to_vars_.size();
589   const int new_size = working_model->constraints_size();
590   CHECK_LE(old_size, new_size);
591   constraint_to_vars_.resize(new_size);
592   constraint_to_linear1_var_.resize(new_size, -1);
593   constraint_to_intervals_.resize(new_size);
594   interval_usage_.resize(new_size);
595   for (int c = old_size; c < new_size; ++c) {
596     AddVariableUsage(c);
597   }
598 }
599 
600 // TODO(user): Also test var_to_constraints_ !!
ConstraintVariableUsageIsConsistent()601 bool PresolveContext::ConstraintVariableUsageIsConsistent() {
602   if (is_unsat_) return true;  // We do not care in this case.
603   if (constraint_to_vars_.size() != working_model->constraints_size()) {
604     LOG(INFO) << "Wrong constraint_to_vars size!";
605     return false;
606   }
607   for (int c = 0; c < constraint_to_vars_.size(); ++c) {
608     if (constraint_to_vars_[c] !=
609         UsedVariables(working_model->constraints(c))) {
610       LOG(INFO) << "Wrong variables usage for constraint: \n"
611                 << ProtobufDebugString(working_model->constraints(c))
612                 << "old_size: " << constraint_to_vars_[c].size();
613       return false;
614     }
615   }
616   int num_in_objective = 0;
617   for (int v = 0; v < var_to_constraints_.size(); ++v) {
618     if (var_to_constraints_[v].contains(kObjectiveConstraint)) {
619       ++num_in_objective;
620       if (!objective_map_.contains(v)) {
621         LOG(INFO) << "Variable " << v
622                   << " is marked as part of the objective but isn't.";
623         return false;
624       }
625     }
626   }
627   if (num_in_objective != objective_map_.size()) {
628     LOG(INFO) << "Not all variables are marked as part of the objective";
629     return false;
630   }
631 
632   return true;
633 }
634 
635 // If a Boolean variable (one with domain [0, 1]) appear in this affine
636 // equivalence class, then we want its representative to be Boolean. Note that
637 // this is always possible because a Boolean variable can never be equal to a
638 // multiple of another if std::abs(coeff) is greater than 1 and if it is not
639 // fixed to zero. This is important because it allows to simply use the same
640 // representative for any referenced literals.
641 //
642 // Note(user): When both domain contains [0,1] and later the wrong variable
643 // become usable as boolean, then we have a bug. Because of that, the code
644 // for GetLiteralRepresentative() is not as simple as it should be.
AddRelation(int x,int y,int64_t c,int64_t o,AffineRelation * repo)645 bool PresolveContext::AddRelation(int x, int y, int64_t c, int64_t o,
646                                   AffineRelation* repo) {
647   // When the coefficient is larger than one, then if later one variable becomes
648   // Boolean, it must be the representative.
649   if (std::abs(c) != 1) return repo->TryAdd(x, y, c, o);
650 
651   CHECK(!VariableWasRemoved(x));
652   CHECK(!VariableWasRemoved(y));
653 
654   // To avoid integer overflow, we always want to use the representative with
655   // the smallest domain magnitude. Otherwise we might express a variable in say
656   // [0, 3] as ([x, x + 3] - x) for an arbitrary large x, and substituting
657   // something like this in a linear expression could break our overflow
658   // precondition.
659   //
660   // Note that if either rep_x or rep_y can be used as a literal, then it will
661   // also be the variable with the smallest domain magnitude (1 or 0 if fixed).
662   const int rep_x = repo->Get(x).representative;
663   const int rep_y = repo->Get(y).representative;
664   const int64_t m_x = std::max(std::abs(MinOf(rep_x)), std::abs(MaxOf(rep_x)));
665   const int64_t m_y = std::max(std::abs(MinOf(rep_y)), std::abs(MaxOf(rep_y)));
666   bool allow_rep_x = m_x < m_y;
667   bool allow_rep_y = m_y < m_x;
668   if (m_x == m_y) {
669     // If both magnitude are the same, we prefer a positive domain.
670     // This is important so we don't use [-1, 0] as a representative for [0, 1].
671     allow_rep_x = MinOf(rep_x) >= MinOf(rep_y);
672     allow_rep_y = MinOf(rep_y) >= MinOf(rep_x);
673   }
674   if (allow_rep_x && allow_rep_y) {
675     // If both representative are okay, we force the choice to the variable
676     // with lower index. This is needed because we have two "equivalence"
677     // relations, and we want the same representative in both.
678     if (rep_x < rep_y) {
679       allow_rep_y = false;
680     } else {
681       allow_rep_x = false;
682     }
683   }
684   return repo->TryAdd(x, y, c, o, allow_rep_x, allow_rep_y);
685 }
686 
687 // Note that we just add the relation to the var_equiv_relations_, not to the
688 // affine one. This is enough, and should prevent overflow in the affine
689 // relation class: if we keep chaining variable fixed to zero, the coefficient
690 // in the relation can overflow. For instance if x = 200 y and z = 200 t,
691 // nothing prevent us if all end up being zero, to say y = z, which will result
692 // in x = 200^2 t. If we make a few bad choices like this, then we can have an
693 // overflow.
ExploitFixedDomain(int var)694 void PresolveContext::ExploitFixedDomain(int var) {
695   DCHECK(RefIsPositive(var));
696   DCHECK(IsFixed(var));
697   const int64_t min = MinOf(var);
698   if (constant_to_ref_.contains(min)) {
699     const int rep = constant_to_ref_[min].Get(this);
700     if (RefIsPositive(rep)) {
701       if (rep != var) {
702         AddRelation(var, rep, 1, 0, &var_equiv_relations_);
703       }
704     } else {
705       if (PositiveRef(rep) == var) {
706         CHECK_EQ(min, 0);
707       } else {
708         AddRelation(var, PositiveRef(rep), -1, 0, &var_equiv_relations_);
709       }
710     }
711   } else {
712     constant_to_ref_[min] = SavedVariable(var);
713   }
714 }
715 
PropagateAffineRelation(int ref)716 bool PresolveContext::PropagateAffineRelation(int ref) {
717   const int var = PositiveRef(ref);
718   const AffineRelation::Relation r = GetAffineRelation(var);
719   if (r.representative == var) return true;
720 
721   // Propagate domains both ways.
722   // var = coeff * rep + offset
723   if (!IntersectDomainWith(r.representative,
724                            DomainOf(var)
725                                .AdditionWith(Domain(-r.offset))
726                                .InverseMultiplicationBy(r.coeff))) {
727     return false;
728   }
729   if (!IntersectDomainWith(var, DomainOf(r.representative)
730                                     .MultiplicationBy(r.coeff)
731                                     .AdditionWith(Domain(r.offset)))) {
732     return false;
733   }
734 
735   return true;
736 }
737 
RemoveAllVariablesFromAffineRelationConstraint()738 void PresolveContext::RemoveAllVariablesFromAffineRelationConstraint() {
739   for (auto& ref_map : var_to_constraints_) {
740     ref_map.erase(kAffineRelationConstraint);
741   }
742 }
743 
744 // We only call that for a non representative variable that is only used in
745 // the kAffineRelationConstraint. Such variable can be ignored and should never
746 // be seen again in the presolve.
RemoveVariableFromAffineRelation(int var)747 void PresolveContext::RemoveVariableFromAffineRelation(int var) {
748   const int rep = GetAffineRelation(var).representative;
749 
750   CHECK(RefIsPositive(var));
751   CHECK_NE(var, rep);
752   CHECK_EQ(var_to_constraints_[var].size(), 1);
753   CHECK(var_to_constraints_[var].contains(kAffineRelationConstraint));
754   CHECK(var_to_constraints_[rep].contains(kAffineRelationConstraint));
755 
756   // We shouldn't reuse this variable again!
757   MarkVariableAsRemoved(var);
758 
759   var_to_constraints_[var].erase(kAffineRelationConstraint);
760   affine_relations_.IgnoreFromClassSize(var);
761   var_equiv_relations_.IgnoreFromClassSize(var);
762 
763   // If the representative is left alone, we can remove it from the special
764   // affine relation constraint too.
765   if (affine_relations_.ClassSize(rep) == 1 &&
766       var_equiv_relations_.ClassSize(rep) == 1) {
767     var_to_constraints_[rep].erase(kAffineRelationConstraint);
768   }
769 
770   if (VLOG_IS_ON(2)) {
771     LOG(INFO) << "Removing affine relation: " << AffineRelationDebugString(var);
772   }
773 }
774 
CanonicalizeVariable(int ref)775 void PresolveContext::CanonicalizeVariable(int ref) {
776   const int var = GetAffineRelation(ref).representative;
777   const int64_t min = MinOf(var);
778   if (min == 0 || IsFixed(var)) return;  // Nothing to do.
779 
780   const int new_var = NewIntVar(DomainOf(var).AdditionWith(Domain(-min)));
781   CHECK(StoreAffineRelation(var, new_var, 1, min, /*debug_no_recursion=*/true));
782   UpdateRuleStats("variables: canonicalize domain");
783   UpdateNewConstraintsVariableUsage();
784 }
785 
ScaleFloatingPointObjective()786 bool PresolveContext::ScaleFloatingPointObjective() {
787   DCHECK(working_model->has_floating_point_objective());
788   DCHECK(!working_model->has_objective());
789   const auto& objective = working_model->floating_point_objective();
790   std::vector<std::pair<int, double>> terms;
791   for (int i = 0; i < objective.vars_size(); ++i) {
792     DCHECK(RefIsPositive(objective.vars(i)));
793     terms.push_back({objective.vars(i), objective.coeffs(i)});
794   }
795   const double offset = objective.offset();
796   const bool maximize = objective.maximize();
797   working_model->clear_floating_point_objective();
798 
799   // We need the domains up to date before scaling.
800   WriteVariableDomainsToProto();
801   return ScaleAndSetObjective(params_, terms, offset, maximize, working_model,
802                               logger_);
803 }
804 
CanonicalizeAffineVariable(int ref,int64_t coeff,int64_t mod,int64_t rhs)805 bool PresolveContext::CanonicalizeAffineVariable(int ref, int64_t coeff,
806                                                  int64_t mod, int64_t rhs) {
807   CHECK_NE(mod, 0);
808   CHECK_NE(coeff, 0);
809 
810   const int64_t gcd = std::gcd(coeff, mod);
811   if (gcd != 1) {
812     if (rhs % gcd != 0) {
813       return NotifyThatModelIsUnsat(
814           absl::StrCat("Infeasible ", coeff, " * X = ", rhs, " % ", mod));
815     }
816     coeff /= gcd;
817     mod /= gcd;
818     rhs /= gcd;
819   }
820 
821   // We just abort in this case as there is no point introducing a new variable.
822   if (std::abs(mod) == 1) return true;
823 
824   int var = ref;
825   if (!RefIsPositive(var)) {
826     var = NegatedRef(ref);
827     coeff = -coeff;
828     rhs = -rhs;
829   }
830 
831   // From var * coeff % mod = rhs
832   // We have var = mod * X + offset.
833   const int64_t offset = ProductWithModularInverse(coeff, mod, rhs);
834 
835   // Lets create a new integer variable and add the affine relation.
836   const Domain new_domain =
837       DomainOf(var).AdditionWith(Domain(-offset)).InverseMultiplicationBy(mod);
838   if (new_domain.IsEmpty()) {
839     return NotifyThatModelIsUnsat(
840         "Empty domain in CanonicalizeAffineVariable()");
841   }
842   if (new_domain.IsFixed()) {
843     UpdateRuleStats("variables: fixed value due to affine relation");
844     return IntersectDomainWith(
845         var, new_domain.ContinuousMultiplicationBy(mod).AdditionWith(
846                  Domain(offset)));
847   }
848 
849   // We make sure the new variable has a domain starting at zero to minimize
850   // future overflow issues. If it end up Boolean, it is also nice to be able to
851   // use it as such.
852   //
853   // A potential problem with this is that it messes up the natural variable
854   // order chosen by the modeler. We try to correct that when mapping variables
855   // at the end of the presolve.
856   const int64_t min_value = new_domain.Min();
857   const int new_var = NewIntVar(new_domain.AdditionWith(Domain(-min_value)));
858   CHECK(StoreAffineRelation(var, new_var, mod, offset + mod * min_value,
859                             /*debug_no_recursion=*/true));
860   UpdateRuleStats("variables: canonicalize affine domain");
861   UpdateNewConstraintsVariableUsage();
862   return true;
863 }
864 
StoreAffineRelation(int ref_x,int ref_y,int64_t coeff,int64_t offset,bool debug_no_recursion)865 bool PresolveContext::StoreAffineRelation(int ref_x, int ref_y, int64_t coeff,
866                                           int64_t offset,
867                                           bool debug_no_recursion) {
868   CHECK_NE(coeff, 0);
869   if (is_unsat_) return false;
870 
871   // TODO(user): I am not 100% sure why, but sometimes the representative is
872   // fixed but that is not propagated to ref_x or ref_y and this causes issues.
873   if (!PropagateAffineRelation(ref_x)) return false;
874   if (!PropagateAffineRelation(ref_y)) return false;
875 
876   if (IsFixed(ref_x)) {
877     const int64_t lhs = DomainOf(ref_x).FixedValue() - offset;
878     if (lhs % std::abs(coeff) != 0) {
879       return NotifyThatModelIsUnsat();
880     }
881     UpdateRuleStats("affine: fixed");
882     return IntersectDomainWith(ref_y, Domain(lhs / coeff));
883   }
884 
885   if (IsFixed(ref_y)) {
886     const int64_t value_x = DomainOf(ref_y).FixedValue() * coeff + offset;
887     UpdateRuleStats("affine: fixed");
888     return IntersectDomainWith(ref_x, Domain(value_x));
889   }
890 
891   // If both are already in the same class, we need to make sure the relations
892   // are compatible.
893   const AffineRelation::Relation rx = GetAffineRelation(ref_x);
894   const AffineRelation::Relation ry = GetAffineRelation(ref_y);
895   if (rx.representative == ry.representative) {
896     // x = rx.coeff * rep + rx.offset;
897     // y = ry.coeff * rep + ry.offset;
898     // And x == coeff * ry.coeff * rep + (coeff * ry.offset + offset).
899     //
900     // So we get the relation a * rep == b with a and b defined here:
901     const int64_t a = coeff * ry.coeff - rx.coeff;
902     const int64_t b = coeff * ry.offset + offset - rx.offset;
903     if (a == 0) {
904       if (b != 0) return NotifyThatModelIsUnsat();
905       return true;
906     }
907     if (b % a != 0) {
908       return NotifyThatModelIsUnsat();
909     }
910     UpdateRuleStats("affine: unique solution");
911     const int64_t unique_value = -b / a;
912     if (!IntersectDomainWith(rx.representative, Domain(unique_value))) {
913       return false;
914     }
915     if (!IntersectDomainWith(ref_x,
916                              Domain(unique_value * rx.coeff + rx.offset))) {
917       return false;
918     }
919     if (!IntersectDomainWith(ref_y,
920                              Domain(unique_value * ry.coeff + ry.offset))) {
921       return false;
922     }
923     return true;
924   }
925 
926   // ref_x = coeff * ref_y + offset;
927   // rx.coeff * rep_x + rx.offset =
928   //    coeff * (ry.coeff * rep_y + ry.offset) + offset
929   //
930   // We have a * rep_x + b * rep_y == o
931   int64_t a = rx.coeff;
932   int64_t b = coeff * ry.coeff;
933   int64_t o = coeff * ry.offset + offset - rx.offset;
934   CHECK_NE(a, 0);
935   CHECK_NE(b, 0);
936   {
937     const int64_t gcd = MathUtil::GCD64(std::abs(a), std::abs(b));
938     if (gcd != 1) {
939       a /= gcd;
940       b /= gcd;
941       if (o % gcd != 0) return NotifyThatModelIsUnsat();
942       o /= gcd;
943     }
944   }
945 
946   // In this (rare) case, we need to canonicalize one of the variable that will
947   // become the representative for both.
948   if (std::abs(a) > 1 && std::abs(b) > 1) {
949     UpdateRuleStats("affine: created common representative");
950     if (!CanonicalizeAffineVariable(rx.representative, a, std::abs(b),
951                                     offset)) {
952       return false;
953     }
954 
955     // Re-add the relation now that a will resolve to a multiple of b.
956     return StoreAffineRelation(ref_x, ref_y, coeff, offset,
957                                /*debug_no_recursion=*/true);
958   }
959 
960   // Canonicalize to x = c * y + o
961   int x, y;
962   int64_t c;
963   bool negate = false;
964   if (std::abs(a) == 1) {
965     x = rx.representative;
966     y = ry.representative;
967     c = b;
968     negate = a < 0;
969   } else {
970     CHECK_EQ(std::abs(b), 1);
971     x = ry.representative;
972     y = rx.representative;
973     c = a;
974     negate = b < 0;
975   }
976   if (negate) {
977     c = -c;
978     o = -o;
979   }
980   CHECK(RefIsPositive(x));
981   CHECK(RefIsPositive(y));
982 
983   // Lets propagate domains first.
984   if (!IntersectDomainWith(
985           y, DomainOf(x).AdditionWith(Domain(-o)).InverseMultiplicationBy(c))) {
986     return false;
987   }
988   if (!IntersectDomainWith(
989           x,
990           DomainOf(y).ContinuousMultiplicationBy(c).AdditionWith(Domain(o)))) {
991     return false;
992   }
993 
994   // To avoid corner cases where replacing x by y in a linear expression
995   // can cause overflow, we might want to canonicalize y first to avoid
996   // cases like x = c * [large_value, ...] - large_value.
997   //
998   // TODO(user): we can do better for overflow by not always choosing the
999   // min at zero, do the best things if it becomes needed.
1000   if (std::abs(o) > std::max(std::abs(MinOf(x)), std::abs(MaxOf(x)))) {
1001     // Both these function recursively call StoreAffineRelation() but shouldn't
1002     // be able to cascade (CHECKED).
1003     CHECK(!debug_no_recursion);
1004     CanonicalizeVariable(y);
1005     return StoreAffineRelation(x, y, c, o, /*debug_no_recursion=*/true);
1006   }
1007 
1008   // TODO(user): can we force the rep and remove GetAffineRelation()?
1009   CHECK(AddRelation(x, y, c, o, &affine_relations_));
1010   if ((c == 1 || c == -1) && o == 0) {
1011     CHECK(AddRelation(x, y, c, o, &var_equiv_relations_));
1012   }
1013 
1014   UpdateRuleStats("affine: new relation");
1015 
1016   // Lets propagate again the new relation. We might as well do it as early
1017   // as possible and not all call site do it.
1018   //
1019   // TODO(user): I am not sure this is needed given the propagation above.
1020   if (!PropagateAffineRelation(ref_x)) return false;
1021   if (!PropagateAffineRelation(ref_y)) return false;
1022 
1023   // These maps should only contains representative, so only need to remap
1024   // either x or y.
1025   const int rep = GetAffineRelation(x).representative;
1026 
1027   // The domain didn't change, but this notification allows to re-process any
1028   // constraint containing these variables. Note that we do not need to
1029   // retrigger a propagation of the constraint containing a variable whose
1030   // representative didn't change.
1031   if (x != rep) modified_domains.Set(x);
1032   if (y != rep) modified_domains.Set(y);
1033 
1034   var_to_constraints_[x].insert(kAffineRelationConstraint);
1035   var_to_constraints_[y].insert(kAffineRelationConstraint);
1036   return true;
1037 }
1038 
StoreBooleanEqualityRelation(int ref_a,int ref_b)1039 bool PresolveContext::StoreBooleanEqualityRelation(int ref_a, int ref_b) {
1040   if (is_unsat_) return false;
1041 
1042   CHECK(!VariableWasRemoved(ref_a));
1043   CHECK(!VariableWasRemoved(ref_b));
1044   CHECK(!DomainOf(ref_a).IsEmpty());
1045   CHECK(!DomainOf(ref_b).IsEmpty());
1046   CHECK(CanBeUsedAsLiteral(ref_a));
1047   CHECK(CanBeUsedAsLiteral(ref_b));
1048 
1049   if (ref_a == ref_b) return true;
1050   if (ref_a == NegatedRef(ref_b)) return IntersectDomainWith(ref_a, Domain(0));
1051 
1052   const int var_a = PositiveRef(ref_a);
1053   const int var_b = PositiveRef(ref_b);
1054   if (RefIsPositive(ref_a) == RefIsPositive(ref_b)) {
1055     // a = b
1056     return StoreAffineRelation(var_a, var_b, /*coeff=*/1, /*offset=*/0);
1057   }
1058   // a = 1 - b
1059   return StoreAffineRelation(var_a, var_b, /*coeff=*/-1, /*offset=*/1);
1060 }
1061 
StoreAbsRelation(int target_ref,int ref)1062 bool PresolveContext::StoreAbsRelation(int target_ref, int ref) {
1063   const auto insert_status = abs_relations_.insert(
1064       std::make_pair(target_ref, SavedVariable(PositiveRef(ref))));
1065   if (!insert_status.second) {
1066     // Tricky: overwrite if the old value refer to a now unused variable.
1067     const int candidate = insert_status.first->second.Get(this);
1068     if (removed_variables_.contains(candidate)) {
1069       insert_status.first->second = SavedVariable(PositiveRef(ref));
1070       return true;
1071     }
1072     return false;
1073   }
1074   return true;
1075 }
1076 
GetAbsRelation(int target_ref,int * ref)1077 bool PresolveContext::GetAbsRelation(int target_ref, int* ref) {
1078   auto it = abs_relations_.find(target_ref);
1079   if (it == abs_relations_.end()) return false;
1080 
1081   // Tricky: In some rare case the stored relation can refer to a deleted
1082   // variable, so we need to ignore it.
1083   //
1084   // TODO(user): Incorporate this as part of SavedVariable/SavedLiteral so we
1085   // make sure we never forget about this.
1086   const int candidate = it->second.Get(this);
1087   if (removed_variables_.contains(candidate)) {
1088     abs_relations_.erase(it);
1089     return false;
1090   }
1091   *ref = candidate;
1092   CHECK(!VariableWasRemoved(*ref));
1093   return true;
1094 }
1095 
GetLiteralRepresentative(int ref) const1096 int PresolveContext::GetLiteralRepresentative(int ref) const {
1097   const AffineRelation::Relation r = GetAffineRelation(PositiveRef(ref));
1098 
1099   CHECK(CanBeUsedAsLiteral(ref));
1100   if (!CanBeUsedAsLiteral(r.representative)) {
1101     // Note(user): This can happen is some corner cases where the affine
1102     // relation where added before the variable became usable as Boolean. When
1103     // this is the case, the domain will be of the form [x, x + 1] and should be
1104     // later remapped to a Boolean variable.
1105     return ref;
1106   }
1107 
1108   // We made sure that the affine representative can always be used as a
1109   // literal. However, if some variable are fixed, we might not have only
1110   // (coeff=1 offset=0) or (coeff=-1 offset=1) and we might have something like
1111   // (coeff=8 offset=0) which is only valid for both variable at zero...
1112   //
1113   // What is sure is that depending on the value, only one mapping can be valid
1114   // because r.coeff can never be zero.
1115   const bool positive_possible = (r.offset == 0 || r.coeff + r.offset == 1);
1116   const bool negative_possible = (r.offset == 1 || r.coeff + r.offset == 0);
1117   DCHECK_NE(positive_possible, negative_possible);
1118   if (RefIsPositive(ref)) {
1119     return positive_possible ? r.representative : NegatedRef(r.representative);
1120   } else {
1121     return positive_possible ? NegatedRef(r.representative) : r.representative;
1122   }
1123 }
1124 
GetVariableRepresentative(int ref) const1125 int PresolveContext::GetVariableRepresentative(int ref) const {
1126   const AffineRelation::Relation r = var_equiv_relations_.Get(PositiveRef(ref));
1127   CHECK_EQ(std::abs(r.coeff), 1);
1128   CHECK_EQ(r.offset, 0);
1129   return RefIsPositive(ref) == (r.coeff == 1) ? r.representative
1130                                               : NegatedRef(r.representative);
1131 }
1132 
1133 // This makes sure that the affine relation only uses one of the
1134 // representative from the var_equiv_relations_.
GetAffineRelation(int ref) const1135 AffineRelation::Relation PresolveContext::GetAffineRelation(int ref) const {
1136   AffineRelation::Relation r = affine_relations_.Get(PositiveRef(ref));
1137   AffineRelation::Relation o = var_equiv_relations_.Get(r.representative);
1138   r.representative = o.representative;
1139   if (o.coeff == -1) r.coeff = -r.coeff;
1140   if (!RefIsPositive(ref)) {
1141     r.coeff *= -1;
1142     r.offset *= -1;
1143   }
1144   return r;
1145 }
1146 
RefDebugString(int ref) const1147 std::string PresolveContext::RefDebugString(int ref) const {
1148   return absl::StrCat(RefIsPositive(ref) ? "X" : "-X", PositiveRef(ref),
1149                       DomainOf(ref).ToString());
1150 }
1151 
AffineRelationDebugString(int ref) const1152 std::string PresolveContext::AffineRelationDebugString(int ref) const {
1153   const AffineRelation::Relation r = GetAffineRelation(ref);
1154   return absl::StrCat(RefDebugString(ref), " = ", r.coeff, " * ",
1155                       RefDebugString(r.representative), " + ", r.offset);
1156 }
1157 
1158 // Create the internal structure for any new variables in working_model.
InitializeNewDomains()1159 void PresolveContext::InitializeNewDomains() {
1160   for (int i = domains.size(); i < working_model->variables_size(); ++i) {
1161     domains.emplace_back(ReadDomainFromProto(working_model->variables(i)));
1162     if (domains.back().IsEmpty()) {
1163       is_unsat_ = true;
1164       return;
1165     }
1166     if (IsFixed(i)) ExploitFixedDomain(i);
1167   }
1168   modified_domains.Resize(domains.size());
1169   var_to_constraints_.resize(domains.size());
1170   var_to_num_linear1_.resize(domains.size());
1171   var_to_ub_only_constraints.resize(domains.size());
1172   var_to_lb_only_constraints.resize(domains.size());
1173 }
1174 
CanonicalizeDomainOfSizeTwo(int var)1175 void PresolveContext::CanonicalizeDomainOfSizeTwo(int var) {
1176   CHECK(RefIsPositive(var));
1177   CHECK_EQ(DomainOf(var).Size(), 2);
1178   const int64_t var_min = MinOf(var);
1179   const int64_t var_max = MaxOf(var);
1180 
1181   if (is_unsat_) return;
1182 
1183   absl::flat_hash_map<int64_t, SavedLiteral>& var_map = encoding_[var];
1184 
1185   // Find encoding for min if present.
1186   auto min_it = var_map.find(var_min);
1187   if (min_it != var_map.end()) {
1188     const int old_var = PositiveRef(min_it->second.Get(this));
1189     if (removed_variables_.contains(old_var)) {
1190       var_map.erase(min_it);
1191       min_it = var_map.end();
1192     }
1193   }
1194 
1195   // Find encoding for max if present.
1196   auto max_it = var_map.find(var_max);
1197   if (max_it != var_map.end()) {
1198     const int old_var = PositiveRef(max_it->second.Get(this));
1199     if (removed_variables_.contains(old_var)) {
1200       var_map.erase(max_it);
1201       max_it = var_map.end();
1202     }
1203   }
1204 
1205   // Insert missing encoding.
1206   int min_literal;
1207   int max_literal;
1208   if (min_it != var_map.end() && max_it != var_map.end()) {
1209     min_literal = min_it->second.Get(this);
1210     max_literal = max_it->second.Get(this);
1211     if (min_literal != NegatedRef(max_literal)) {
1212       UpdateRuleStats("variables with 2 values: merge encoding literals");
1213       StoreBooleanEqualityRelation(min_literal, NegatedRef(max_literal));
1214       if (is_unsat_) return;
1215     }
1216     min_literal = GetLiteralRepresentative(min_literal);
1217     max_literal = GetLiteralRepresentative(max_literal);
1218     if (!IsFixed(min_literal)) CHECK_EQ(min_literal, NegatedRef(max_literal));
1219   } else if (min_it != var_map.end() && max_it == var_map.end()) {
1220     UpdateRuleStats("variables with 2 values: register other encoding");
1221     min_literal = min_it->second.Get(this);
1222     max_literal = NegatedRef(min_literal);
1223     var_map[var_max] = SavedLiteral(max_literal);
1224   } else if (min_it == var_map.end() && max_it != var_map.end()) {
1225     UpdateRuleStats("variables with 2 values: register other encoding");
1226     max_literal = max_it->second.Get(this);
1227     min_literal = NegatedRef(max_literal);
1228     var_map[var_min] = SavedLiteral(min_literal);
1229   } else {
1230     UpdateRuleStats("variables with 2 values: create encoding literal");
1231     max_literal = NewBoolVar();
1232     min_literal = NegatedRef(max_literal);
1233     var_map[var_min] = SavedLiteral(min_literal);
1234     var_map[var_max] = SavedLiteral(max_literal);
1235   }
1236 
1237   if (IsFixed(min_literal) || IsFixed(max_literal)) {
1238     CHECK(IsFixed(min_literal));
1239     CHECK(IsFixed(max_literal));
1240     UpdateRuleStats("variables with 2 values: fixed encoding");
1241     if (LiteralIsTrue(min_literal)) {
1242       return static_cast<void>(IntersectDomainWith(var, Domain(var_min)));
1243     } else {
1244       return static_cast<void>(IntersectDomainWith(var, Domain(var_max)));
1245     }
1246   }
1247 
1248   // Add affine relation.
1249   if (GetAffineRelation(var).representative != PositiveRef(min_literal)) {
1250     UpdateRuleStats("variables with 2 values: new affine relation");
1251     if (RefIsPositive(max_literal)) {
1252       (void)StoreAffineRelation(var, PositiveRef(max_literal),
1253                                 var_max - var_min, var_min);
1254     } else {
1255       (void)StoreAffineRelation(var, PositiveRef(max_literal),
1256                                 var_min - var_max, var_max);
1257     }
1258   }
1259 }
1260 
InsertVarValueEncodingInternal(int literal,int var,int64_t value,bool add_constraints)1261 void PresolveContext::InsertVarValueEncodingInternal(int literal, int var,
1262                                                      int64_t value,
1263                                                      bool add_constraints) {
1264   CHECK(RefIsPositive(var));
1265   CHECK(!VariableWasRemoved(literal));
1266   CHECK(!VariableWasRemoved(var));
1267   absl::flat_hash_map<int64_t, SavedLiteral>& var_map = encoding_[var];
1268 
1269   // The code below is not 100% correct if this is not the case.
1270   DCHECK(DomainOf(var).Contains(value));
1271 
1272   // If an encoding already exist, make the two Boolean equals.
1273   const auto [it, inserted] =
1274       var_map.insert(std::make_pair(value, SavedLiteral(literal)));
1275   if (!inserted) {
1276     const int previous_literal = it->second.Get(this);
1277 
1278     // Ticky and rare: I have only observed this on the LNS of
1279     // radiation_m18_12_05_sat.fzn. The value was encoded, but maybe we never
1280     // used the involved variables / constraints, so it was removed (with the
1281     // encoding constraints) from the model already! We have to be careful.
1282     if (VariableWasRemoved(previous_literal)) {
1283       it->second = SavedLiteral(literal);
1284     } else {
1285       if (literal != previous_literal) {
1286         UpdateRuleStats(
1287             "variables: merge equivalent var value encoding literals");
1288         StoreBooleanEqualityRelation(literal, previous_literal);
1289       }
1290     }
1291     return;
1292   }
1293 
1294   if (DomainOf(var).Size() == 2) {
1295     // TODO(user): There is a bug here if the var == value was not in the
1296     // domain, it will just be ignored.
1297     CanonicalizeDomainOfSizeTwo(var);
1298   } else {
1299     VLOG(2) << "Insert lit(" << literal << ") <=> var(" << var
1300             << ") == " << value;
1301     eq_half_encoding_[var][value].insert(literal);
1302     neq_half_encoding_[var][value].insert(NegatedRef(literal));
1303     if (add_constraints) {
1304       UpdateRuleStats("variables: add encoding constraint");
1305       AddImplyInDomain(literal, var, Domain(value));
1306       AddImplyInDomain(NegatedRef(literal), var, Domain(value).Complement());
1307     }
1308   }
1309 }
1310 
InsertHalfVarValueEncoding(int literal,int var,int64_t value,bool imply_eq)1311 bool PresolveContext::InsertHalfVarValueEncoding(int literal, int var,
1312                                                  int64_t value, bool imply_eq) {
1313   if (is_unsat_) return false;
1314   CHECK(RefIsPositive(var));
1315 
1316   // Creates the linking sets on demand.
1317   // Insert the enforcement literal in the half encoding map.
1318   auto& direct_set =
1319       imply_eq ? eq_half_encoding_[var][value] : neq_half_encoding_[var][value];
1320   if (!direct_set.insert(literal).second) return false;  // Already there.
1321 
1322   VLOG(2) << "Collect lit(" << literal << ") implies var(" << var
1323           << (imply_eq ? ") == " : ") != ") << value;
1324   UpdateRuleStats("variables: detect half reified value encoding");
1325 
1326   // Note(user): We don't expect a lot of literals in these sets, so doing
1327   // a scan should be okay.
1328   auto& other_set =
1329       imply_eq ? neq_half_encoding_[var][value] : eq_half_encoding_[var][value];
1330   for (const int other : other_set) {
1331     if (GetLiteralRepresentative(other) != NegatedRef(literal)) continue;
1332 
1333     UpdateRuleStats("variables: detect fully reified value encoding");
1334     const int imply_eq_literal = imply_eq ? literal : NegatedRef(literal);
1335     InsertVarValueEncodingInternal(imply_eq_literal, var, value,
1336                                    /*add_constraints=*/false);
1337     break;
1338   }
1339 
1340   return true;
1341 }
1342 
CanonicalizeEncoding(int * ref,int64_t * value)1343 bool PresolveContext::CanonicalizeEncoding(int* ref, int64_t* value) {
1344   const AffineRelation::Relation r = GetAffineRelation(*ref);
1345   if ((*value - r.offset) % r.coeff != 0) return false;
1346   *ref = r.representative;
1347   *value = (*value - r.offset) / r.coeff;
1348   return true;
1349 }
1350 
InsertVarValueEncoding(int literal,int ref,int64_t value)1351 bool PresolveContext::InsertVarValueEncoding(int literal, int ref,
1352                                              int64_t value) {
1353   if (!CanonicalizeEncoding(&ref, &value)) {
1354     return SetLiteralToFalse(literal);
1355   }
1356   literal = GetLiteralRepresentative(literal);
1357   InsertVarValueEncodingInternal(literal, ref, value, /*add_constraints=*/true);
1358   return true;
1359 }
1360 
StoreLiteralImpliesVarEqValue(int literal,int var,int64_t value)1361 bool PresolveContext::StoreLiteralImpliesVarEqValue(int literal, int var,
1362                                                     int64_t value) {
1363   if (!CanonicalizeEncoding(&var, &value)) return false;
1364   literal = GetLiteralRepresentative(literal);
1365   return InsertHalfVarValueEncoding(literal, var, value, /*imply_eq=*/true);
1366 }
1367 
StoreLiteralImpliesVarNEqValue(int literal,int var,int64_t value)1368 bool PresolveContext::StoreLiteralImpliesVarNEqValue(int literal, int var,
1369                                                      int64_t value) {
1370   if (!CanonicalizeEncoding(&var, &value)) return false;
1371   literal = GetLiteralRepresentative(literal);
1372   return InsertHalfVarValueEncoding(literal, var, value, /*imply_eq=*/false);
1373 }
1374 
HasVarValueEncoding(int ref,int64_t value,int * literal)1375 bool PresolveContext::HasVarValueEncoding(int ref, int64_t value,
1376                                           int* literal) {
1377   CHECK(!VariableWasRemoved(ref));
1378   if (!CanonicalizeEncoding(&ref, &value)) return false;
1379   const absl::flat_hash_map<int64_t, SavedLiteral>& var_map = encoding_[ref];
1380   const auto it = var_map.find(value);
1381   if (it != var_map.end()) {
1382     if (literal != nullptr) {
1383       *literal = it->second.Get(this);
1384     }
1385     return true;
1386   }
1387   return false;
1388 }
1389 
IsFullyEncoded(int ref) const1390 bool PresolveContext::IsFullyEncoded(int ref) const {
1391   const int var = PositiveRef(ref);
1392   const int64_t size = domains[var].Size();
1393   if (size <= 2) return true;
1394   const auto& it = encoding_.find(var);
1395   return it == encoding_.end() ? false : size <= it->second.size();
1396 }
1397 
IsFullyEncoded(const LinearExpressionProto & expr) const1398 bool PresolveContext::IsFullyEncoded(const LinearExpressionProto& expr) const {
1399   CHECK_LE(expr.vars_size(), 1);
1400   if (IsFixed(expr)) return true;
1401   return IsFullyEncoded(expr.vars(0));
1402 }
1403 
GetOrCreateVarValueEncoding(int ref,int64_t value)1404 int PresolveContext::GetOrCreateVarValueEncoding(int ref, int64_t value) {
1405   CHECK(!VariableWasRemoved(ref));
1406   if (!CanonicalizeEncoding(&ref, &value)) return GetOrCreateConstantVar(0);
1407 
1408   // Positive after CanonicalizeEncoding().
1409   const int var = ref;
1410 
1411   // Returns the false literal if the value is not in the domain.
1412   if (!domains[var].Contains(value)) {
1413     return GetOrCreateConstantVar(0);
1414   }
1415 
1416   // Returns the associated literal if already present.
1417   absl::flat_hash_map<int64_t, SavedLiteral>& var_map = encoding_[var];
1418   auto it = var_map.find(value);
1419   if (it != var_map.end()) {
1420     const int lit = it->second.Get(this);
1421     if (VariableWasRemoved(lit)) {
1422       // If the variable was already removed, for now we create a new one.
1423       // This should be rare hopefully.
1424       var_map.erase(value);
1425     } else {
1426       return lit;
1427     }
1428   }
1429 
1430   // Special case for fixed domains.
1431   if (domains[var].Size() == 1) {
1432     const int true_literal = GetOrCreateConstantVar(1);
1433     var_map[value] = SavedLiteral(true_literal);
1434     return true_literal;
1435   }
1436 
1437   // Special case for domains of size 2.
1438   const int64_t var_min = MinOf(var);
1439   const int64_t var_max = MaxOf(var);
1440   if (domains[var].Size() == 2) {
1441     // Checks if the other value is already encoded.
1442     const int64_t other_value = value == var_min ? var_max : var_min;
1443     auto other_it = var_map.find(other_value);
1444     if (other_it != var_map.end()) {
1445       const int literal = NegatedRef(other_it->second.Get(this));
1446       if (VariableWasRemoved(literal)) {
1447         // If the variable was already removed, for now we create a new one.
1448         // This should be rare hopefully.
1449         var_map.erase(other_value);
1450       } else {
1451         // Update the encoding map. The domain could have been reduced to size
1452         // two after the creation of the first literal.
1453         var_map[value] = SavedLiteral(literal);
1454         return literal;
1455       }
1456     }
1457 
1458     if (var_min == 0 && var_max == 1) {
1459       const int representative = GetLiteralRepresentative(var);
1460       var_map[1] = SavedLiteral(representative);
1461       var_map[0] = SavedLiteral(NegatedRef(representative));
1462       return value == 1 ? representative : NegatedRef(representative);
1463     } else {
1464       const int literal = NewBoolVar();
1465       InsertVarValueEncoding(literal, var, var_max);
1466       const int representative = GetLiteralRepresentative(literal);
1467       return value == var_max ? representative : NegatedRef(representative);
1468     }
1469   }
1470 
1471   const int literal = NewBoolVar();
1472   InsertVarValueEncoding(literal, var, value);
1473   return GetLiteralRepresentative(literal);
1474 }
1475 
GetOrCreateAffineValueEncoding(const LinearExpressionProto & expr,int64_t value)1476 int PresolveContext::GetOrCreateAffineValueEncoding(
1477     const LinearExpressionProto& expr, int64_t value) {
1478   DCHECK_LE(expr.vars_size(), 1);
1479   if (IsFixed(expr)) {
1480     if (FixedValue(expr) == value) {
1481       return GetOrCreateConstantVar(1);
1482     } else {
1483       return GetOrCreateConstantVar(0);
1484     }
1485   }
1486 
1487   if ((value - expr.offset()) % expr.coeffs(0) != 0) {
1488     return GetOrCreateConstantVar(0);
1489   }
1490 
1491   return GetOrCreateVarValueEncoding(expr.vars(0),
1492                                      (value - expr.offset()) / expr.coeffs(0));
1493 }
1494 
ReadObjectiveFromProto()1495 void PresolveContext::ReadObjectiveFromProto() {
1496   const CpObjectiveProto& obj = working_model->objective();
1497 
1498   objective_offset_ = obj.offset();
1499   objective_scaling_factor_ = obj.scaling_factor();
1500   if (objective_scaling_factor_ == 0.0) {
1501     objective_scaling_factor_ = 1.0;
1502   }
1503 
1504   objective_integer_offset_ = obj.integer_offset();
1505   objective_integer_scaling_factor_ = obj.integer_scaling_factor();
1506   if (objective_integer_scaling_factor_ == 0) {
1507     objective_integer_scaling_factor_ = 1;
1508   }
1509 
1510   if (!obj.domain().empty()) {
1511     // We might relax this in CanonicalizeObjective() when we will compute
1512     // the possible objective domain from the domains of the variables.
1513     objective_domain_is_constraining_ = true;
1514     objective_domain_ = ReadDomainFromProto(obj);
1515   } else {
1516     objective_domain_is_constraining_ = false;
1517     objective_domain_ = Domain::AllValues();
1518   }
1519 
1520   // This is an upper bound of the higher magnitude that can be reach by
1521   // summing an objective partial sum. Because of the model validation, this
1522   // shouldn't overflow, and we make sure it stays this way.
1523   objective_overflow_detection_ = 0;
1524 
1525   objective_map_.clear();
1526   for (int i = 0; i < obj.vars_size(); ++i) {
1527     const int ref = obj.vars(i);
1528     const int64_t var_max_magnitude =
1529         std::max(std::abs(MinOf(ref)), std::abs(MaxOf(ref)));
1530 
1531     // Skipping var fixed to zero allow to avoid some overflow in situation
1532     // were we can deal with it.
1533     if (var_max_magnitude == 0) continue;
1534 
1535     const int64_t coeff = obj.coeffs(i);
1536     objective_overflow_detection_ += var_max_magnitude * std::abs(coeff);
1537 
1538     const int var = PositiveRef(ref);
1539     objective_map_[var] += RefIsPositive(ref) ? coeff : -coeff;
1540     if (objective_map_[var] == 0) {
1541       objective_map_.erase(var);
1542       var_to_constraints_[var].erase(kObjectiveConstraint);
1543     } else {
1544       var_to_constraints_[var].insert(kObjectiveConstraint);
1545     }
1546   }
1547 }
1548 
CanonicalizeObjective(bool simplify_domain)1549 bool PresolveContext::CanonicalizeObjective(bool simplify_domain) {
1550   int64_t offset_change = 0;
1551 
1552   // We replace each entry by its affine representative.
1553   // Note that the non-deterministic loop is fine, but because we iterate
1554   // one the map while modifying it, it is safer to do a copy rather than to
1555   // try to handle that in one pass.
1556   tmp_entries_.clear();
1557   for (const auto& entry : objective_map_) {
1558     tmp_entries_.push_back(entry);
1559   }
1560 
1561   // TODO(user): This is a bit duplicated with the presolve linear code.
1562   // We also do not propagate back any domain restriction from the objective to
1563   // the variables if any.
1564   for (const auto& entry : tmp_entries_) {
1565     const int var = entry.first;
1566     const auto it = objective_map_.find(var);
1567     if (it == objective_map_.end()) continue;
1568     const int64_t coeff = it->second;
1569 
1570     // If a variable only appear in objective, we can fix it!
1571     // Note that we don't care if it was in affine relation, because if none
1572     // of the relations are left, then we can still fix it.
1573     if (!keep_all_feasible_solutions && !objective_domain_is_constraining_ &&
1574         ConstraintVariableGraphIsUpToDate() &&
1575         var_to_constraints_[var].size() == 1 &&
1576         var_to_constraints_[var].contains(kObjectiveConstraint)) {
1577       UpdateRuleStats("objective: variable not used elsewhere");
1578       if (coeff > 0) {
1579         if (!IntersectDomainWith(var, Domain(MinOf(var)))) {
1580           return false;
1581         }
1582       } else {
1583         if (!IntersectDomainWith(var, Domain(MaxOf(var)))) {
1584           return false;
1585         }
1586       }
1587     }
1588 
1589     if (IsFixed(var)) {
1590       offset_change += coeff * MinOf(var);
1591       var_to_constraints_[var].erase(kObjectiveConstraint);
1592       objective_map_.erase(var);
1593       continue;
1594     }
1595 
1596     const AffineRelation::Relation r = GetAffineRelation(var);
1597     if (r.representative == var) continue;
1598 
1599     objective_map_.erase(var);
1600     var_to_constraints_[var].erase(kObjectiveConstraint);
1601 
1602     // Do the substitution.
1603     offset_change += coeff * r.offset;
1604     const int64_t new_coeff = objective_map_[r.representative] +=
1605         coeff * r.coeff;
1606 
1607     // Process new term.
1608     if (new_coeff == 0) {
1609       objective_map_.erase(r.representative);
1610       var_to_constraints_[r.representative].erase(kObjectiveConstraint);
1611     } else {
1612       var_to_constraints_[r.representative].insert(kObjectiveConstraint);
1613       if (IsFixed(r.representative)) {
1614         offset_change += new_coeff * MinOf(r.representative);
1615         var_to_constraints_[r.representative].erase(kObjectiveConstraint);
1616         objective_map_.erase(r.representative);
1617       }
1618     }
1619   }
1620 
1621   Domain implied_domain(0);
1622   int64_t gcd(0);
1623 
1624   // We need to sort the entries to be deterministic.
1625   tmp_entries_.clear();
1626   for (const auto& entry : objective_map_) {
1627     tmp_entries_.push_back(entry);
1628   }
1629   std::sort(tmp_entries_.begin(), tmp_entries_.end());
1630   for (const auto& entry : tmp_entries_) {
1631     const int var = entry.first;
1632     const int64_t coeff = entry.second;
1633     gcd = MathUtil::GCD64(gcd, std::abs(coeff));
1634     implied_domain =
1635         implied_domain.AdditionWith(DomainOf(var).MultiplicationBy(coeff))
1636             .RelaxIfTooComplex();
1637   }
1638 
1639   // This is the new domain.
1640   // Note that the domain never include the offset.
1641   objective_domain_ = objective_domain_.AdditionWith(Domain(-offset_change))
1642                           .IntersectionWith(implied_domain);
1643 
1644   // Depending on the use case, we cannot do that.
1645   if (simplify_domain) {
1646     objective_domain_ =
1647         objective_domain_.SimplifyUsingImpliedDomain(implied_domain);
1648   }
1649 
1650   // Update the offset.
1651   objective_offset_ += offset_change;
1652   objective_integer_offset_ +=
1653       offset_change * objective_integer_scaling_factor_;
1654 
1655   // Maybe divide by GCD.
1656   if (gcd > 1) {
1657     for (auto& entry : objective_map_) {
1658       entry.second /= gcd;
1659     }
1660     objective_domain_ = objective_domain_.InverseMultiplicationBy(gcd);
1661     objective_offset_ /= static_cast<double>(gcd);
1662     objective_scaling_factor_ *= static_cast<double>(gcd);
1663     objective_integer_scaling_factor_ *= gcd;
1664   }
1665 
1666   if (objective_domain_.IsEmpty()) return false;
1667 
1668   // Detect if the objective domain do not limit the "optimal" objective value.
1669   // If this is true, then we can apply any reduction that reduce the objective
1670   // value without any issues.
1671   objective_domain_is_constraining_ =
1672       !implied_domain
1673            .IntersectionWith(Domain(std::numeric_limits<int64_t>::min(),
1674                                     objective_domain_.Max()))
1675            .IsIncludedIn(objective_domain_);
1676   return true;
1677 }
1678 
RemoveVariableFromObjective(int var)1679 void PresolveContext::RemoveVariableFromObjective(int var) {
1680   objective_map_.erase(var);
1681   var_to_constraints_[var].erase(kObjectiveConstraint);
1682 }
1683 
AddToObjective(int var,int64_t value)1684 void PresolveContext::AddToObjective(int var, int64_t value) {
1685   int64_t& map_ref = objective_map_[var];
1686   map_ref += value;
1687   if (map_ref == 0) {
1688     objective_map_.erase(var);
1689     var_to_constraints_[var].erase(kObjectiveConstraint);
1690   } else {
1691     var_to_constraints_[var].insert(kObjectiveConstraint);
1692   }
1693 }
1694 
AddToObjectiveOffset(int64_t value)1695 void PresolveContext::AddToObjectiveOffset(int64_t value) {
1696   // Tricky: The objective domain is without the offset, so we need to shift it.
1697   objective_offset_ += static_cast<double>(value);
1698   objective_integer_offset_ += value * objective_integer_scaling_factor_;
1699   objective_domain_ = objective_domain_.AdditionWith(Domain(-value));
1700 }
1701 
SubstituteVariableInObjective(int var_in_equality,int64_t coeff_in_equality,const ConstraintProto & equality,std::vector<int> * new_vars_in_objective)1702 bool PresolveContext::SubstituteVariableInObjective(
1703     int var_in_equality, int64_t coeff_in_equality,
1704     const ConstraintProto& equality, std::vector<int>* new_vars_in_objective) {
1705   CHECK(equality.enforcement_literal().empty());
1706   CHECK(RefIsPositive(var_in_equality));
1707 
1708   if (new_vars_in_objective != nullptr) new_vars_in_objective->clear();
1709 
1710   // We can only "easily" substitute if the objective coefficient is a multiple
1711   // of the one in the constraint.
1712   const int64_t coeff_in_objective =
1713       gtl::FindOrDie(objective_map_, var_in_equality);
1714   CHECK_NE(coeff_in_equality, 0);
1715   CHECK_EQ(coeff_in_objective % coeff_in_equality, 0);
1716 
1717   const int64_t multiplier = coeff_in_objective / coeff_in_equality;
1718 
1719   // Abort if the new objective seems to violate our overflow preconditions.
1720   int64_t change = 0;
1721   for (int i = 0; i < equality.linear().vars().size(); ++i) {
1722     int var = equality.linear().vars(i);
1723     if (PositiveRef(var) == var_in_equality) continue;
1724     int64_t coeff = equality.linear().coeffs(i);
1725     change +=
1726         std::abs(coeff) * std::max(std::abs(MinOf(var)), std::abs(MaxOf(var)));
1727   }
1728   const int64_t new_value =
1729       CapAdd(CapProd(std::abs(multiplier), change),
1730              objective_overflow_detection_ -
1731                  std::abs(coeff_in_equality) *
1732                      std::max(std::abs(MinOf(var_in_equality)),
1733                               std::abs(MaxOf(var_in_equality))));
1734   if (new_value == std::numeric_limits<int64_t>::max()) return false;
1735   objective_overflow_detection_ = new_value;
1736 
1737   // Compute the objective offset change.
1738   Domain offset = ReadDomainFromProto(equality.linear());
1739   DCHECK_EQ(offset.Min(), offset.Max());
1740   bool exact = true;
1741   offset = offset.MultiplicationBy(multiplier, &exact);
1742   CHECK(exact);
1743   CHECK(!offset.IsEmpty());
1744 
1745   // We also need to make sure the integer_offset will not overflow.
1746   {
1747     int64_t temp = CapProd(offset.Min(), objective_integer_scaling_factor_);
1748     if (temp == std::numeric_limits<int64_t>::max()) return false;
1749     if (temp == std::numeric_limits<int64_t>::min()) return false;
1750     temp = CapAdd(temp, objective_integer_offset_);
1751     if (temp == std::numeric_limits<int64_t>::max()) return false;
1752     if (temp == std::numeric_limits<int64_t>::min()) return false;
1753   }
1754 
1755   // Perform the substitution.
1756   for (int i = 0; i < equality.linear().vars().size(); ++i) {
1757     int var = equality.linear().vars(i);
1758     int64_t coeff = equality.linear().coeffs(i);
1759     if (!RefIsPositive(var)) {
1760       var = NegatedRef(var);
1761       coeff = -coeff;
1762     }
1763     if (var == var_in_equality) continue;
1764 
1765     int64_t& map_ref = objective_map_[var];
1766     if (map_ref == 0 && new_vars_in_objective != nullptr) {
1767       new_vars_in_objective->push_back(var);
1768     }
1769     map_ref -= coeff * multiplier;
1770 
1771     if (map_ref == 0) {
1772       objective_map_.erase(var);
1773       var_to_constraints_[var].erase(kObjectiveConstraint);
1774     } else {
1775       var_to_constraints_[var].insert(kObjectiveConstraint);
1776     }
1777   }
1778 
1779   objective_map_.erase(var_in_equality);
1780   var_to_constraints_[var_in_equality].erase(kObjectiveConstraint);
1781 
1782   // Tricky: The objective domain is without the offset, so we need to shift it.
1783   objective_offset_ += static_cast<double>(offset.Min());
1784   objective_integer_offset_ += offset.Min() * objective_integer_scaling_factor_;
1785   objective_domain_ = objective_domain_.AdditionWith(Domain(-offset.Min()));
1786 
1787   // Because we can assume that the constraint we used was constraining
1788   // (otherwise it would have been removed), the objective domain should be now
1789   // constraining.
1790   objective_domain_is_constraining_ = true;
1791 
1792   if (objective_domain_.IsEmpty()) {
1793     return NotifyThatModelIsUnsat();
1794   }
1795   return true;
1796 }
1797 
ExploitExactlyOneInObjective(absl::Span<const int> exactly_one)1798 bool PresolveContext::ExploitExactlyOneInObjective(
1799     absl::Span<const int> exactly_one) {
1800   if (objective_map_.empty()) return false;
1801 
1802   int64_t min_coeff = std::numeric_limits<int64_t>::max();
1803   for (const int ref : exactly_one) {
1804     const auto it = objective_map_.find(PositiveRef(ref));
1805     if (it == objective_map_.end()) return false;
1806 
1807     const int64_t coeff = it->second;
1808     if (RefIsPositive(ref)) {
1809       min_coeff = std::min(min_coeff, coeff);
1810     } else {
1811       // Objective = coeff * var = coeff * (1 - ref);
1812       min_coeff = std::min(min_coeff, -coeff);
1813     }
1814   }
1815 
1816   int64_t offset = min_coeff;
1817   for (const int ref : exactly_one) {
1818     const int var = PositiveRef(ref);
1819     int64_t& map_ref = objective_map_.at(var);
1820     if (RefIsPositive(ref)) {
1821       map_ref -= min_coeff;
1822       if (map_ref == 0) {
1823         objective_map_.erase(var);
1824         var_to_constraints_[var].erase(kObjectiveConstraint);
1825       }
1826     } else {
1827       // Term = coeff * (1 - X) = coeff  - coeff * X;
1828       // So -coeff -> -coeff  -min_coeff
1829       // And Term = coeff + min_coeff - min_coeff - (coeff + min_coeff) * X
1830       //          = (coeff + min_coeff) * (1 - X) - min_coeff;
1831       map_ref += min_coeff;
1832       if (map_ref == 0) {
1833         objective_map_.erase(var);
1834         var_to_constraints_[var].erase(kObjectiveConstraint);
1835       }
1836       offset -= min_coeff;
1837     }
1838   }
1839 
1840   // Note that the domain never include the offset, so we need to update it.
1841   if (offset != 0) {
1842     objective_offset_ += offset;
1843     objective_integer_offset_ += offset * objective_integer_scaling_factor_;
1844     objective_domain_ = objective_domain_.AdditionWith(Domain(-offset));
1845   }
1846 
1847   return true;
1848 }
1849 
WriteObjectiveToProto() const1850 void PresolveContext::WriteObjectiveToProto() const {
1851   // We need to sort the entries to be deterministic.
1852   std::vector<std::pair<int, int64_t>> entries;
1853   for (const auto& entry : objective_map_) {
1854     entries.push_back(entry);
1855   }
1856   std::sort(entries.begin(), entries.end());
1857 
1858   CpObjectiveProto* mutable_obj = working_model->mutable_objective();
1859   mutable_obj->set_offset(objective_offset_);
1860   mutable_obj->set_scaling_factor(objective_scaling_factor_);
1861   mutable_obj->set_integer_offset(objective_integer_offset_);
1862   if (objective_integer_scaling_factor_ == 1) {
1863     mutable_obj->set_integer_scaling_factor(0);  // Default.
1864   } else {
1865     mutable_obj->set_integer_scaling_factor(objective_integer_scaling_factor_);
1866   }
1867   FillDomainInProto(objective_domain_, mutable_obj);
1868   mutable_obj->clear_vars();
1869   mutable_obj->clear_coeffs();
1870   for (const auto& entry : entries) {
1871     mutable_obj->add_vars(entry.first);
1872     mutable_obj->add_coeffs(entry.second);
1873   }
1874 }
1875 
WriteVariableDomainsToProto() const1876 void PresolveContext::WriteVariableDomainsToProto() const {
1877   for (int i = 0; i < working_model->variables_size(); ++i) {
1878     FillDomainInProto(DomainOf(i), working_model->mutable_variables(i));
1879   }
1880 }
1881 
GetOrCreateReifiedPrecedenceLiteral(const LinearExpressionProto & time_i,const LinearExpressionProto & time_j,int active_i,int active_j)1882 int PresolveContext::GetOrCreateReifiedPrecedenceLiteral(
1883     const LinearExpressionProto& time_i, const LinearExpressionProto& time_j,
1884     int active_i, int active_j) {
1885   CHECK(!LiteralIsFalse(active_i));
1886   CHECK(!LiteralIsFalse(active_j));
1887   DCHECK(ExpressionIsAffine(time_i));
1888   DCHECK(ExpressionIsAffine(time_j));
1889 
1890   const std::tuple<int, int64_t, int, int64_t, int64_t, int, int> key =
1891       GetReifiedPrecedenceKey(time_i, time_j, active_i, active_j);
1892   const auto& it = reified_precedences_cache_.find(key);
1893   if (it != reified_precedences_cache_.end()) return it->second;
1894 
1895   const int result = NewBoolVar();
1896   reified_precedences_cache_[key] = result;
1897 
1898   // result => (time_i <= time_j) && active_i && active_j.
1899   ConstraintProto* const lesseq = working_model->add_constraints();
1900   lesseq->add_enforcement_literal(result);
1901   if (!IsFixed(time_i)) {
1902     lesseq->mutable_linear()->add_vars(time_i.vars(0));
1903     lesseq->mutable_linear()->add_coeffs(-time_i.coeffs(0));
1904   }
1905   if (!IsFixed(time_j)) {
1906     lesseq->mutable_linear()->add_vars(time_j.vars(0));
1907     lesseq->mutable_linear()->add_coeffs(time_j.coeffs(0));
1908   }
1909 
1910   const int64_t offset =
1911       (IsFixed(time_i) ? FixedValue(time_i) : time_i.offset()) -
1912       (IsFixed(time_j) ? FixedValue(time_j) : time_j.offset());
1913   lesseq->mutable_linear()->add_domain(offset);
1914   lesseq->mutable_linear()->add_domain(std::numeric_limits<int64_t>::max());
1915   if (!LiteralIsTrue(active_i)) {
1916     AddImplication(result, active_i);
1917   }
1918   if (!LiteralIsTrue(active_j)) {
1919     AddImplication(result, active_j);
1920   }
1921 
1922   // Not(result) && active_i && active_j => (time_i > time_j)
1923   ConstraintProto* const greater = working_model->add_constraints();
1924   if (!IsFixed(time_i)) {
1925     greater->mutable_linear()->add_vars(time_i.vars(0));
1926     greater->mutable_linear()->add_coeffs(-time_i.coeffs(0));
1927   }
1928   if (!IsFixed(time_j)) {
1929     greater->mutable_linear()->add_vars(time_j.vars(0));
1930     greater->mutable_linear()->add_coeffs(time_j.coeffs(0));
1931   }
1932   greater->mutable_linear()->add_domain(std::numeric_limits<int64_t>::min());
1933   greater->mutable_linear()->add_domain(offset - 1);
1934 
1935   // Manages enforcement literal.
1936   greater->add_enforcement_literal(NegatedRef(result));
1937   if (!LiteralIsTrue(active_i)) {
1938     greater->add_enforcement_literal(active_i);
1939   }
1940   if (!LiteralIsTrue(active_j)) {
1941     greater->add_enforcement_literal(active_j);
1942   }
1943 
1944   // This is redundant but should improves performance.
1945   //
1946   // If GetOrCreateReifiedPrecedenceLiteral(time_j, time_i, active_j, active_i)
1947   // (the reverse precedence) has been called too, then we can link the two
1948   // precedence literals, and the two active literals together.
1949   const auto& rev_it = reified_precedences_cache_.find(
1950       GetReifiedPrecedenceKey(time_j, time_i, active_j, active_i));
1951   if (rev_it != reified_precedences_cache_.end()) {
1952     auto* const bool_or = working_model->add_constraints()->mutable_bool_or();
1953     bool_or->add_literals(result);
1954     bool_or->add_literals(rev_it->second);
1955     bool_or->add_literals(NegatedRef(active_i));
1956     bool_or->add_literals(NegatedRef(active_j));
1957   }
1958 
1959   return result;
1960 }
1961 
1962 std::tuple<int, int64_t, int, int64_t, int64_t, int, int>
GetReifiedPrecedenceKey(const LinearExpressionProto & time_i,const LinearExpressionProto & time_j,int active_i,int active_j)1963 PresolveContext::GetReifiedPrecedenceKey(const LinearExpressionProto& time_i,
1964                                          const LinearExpressionProto& time_j,
1965                                          int active_i, int active_j) {
1966   const int var_i =
1967       IsFixed(time_i) ? std::numeric_limits<int>::min() : time_i.vars(0);
1968   const int64_t coeff_i = IsFixed(time_i) ? 0 : time_i.coeffs(0);
1969   const int var_j =
1970       IsFixed(time_j) ? std::numeric_limits<int>::min() : time_j.vars(0);
1971   const int64_t coeff_j = IsFixed(time_j) ? 0 : time_j.coeffs(0);
1972   const int64_t offset =
1973       (IsFixed(time_i) ? FixedValue(time_i) : time_i.offset()) -
1974       (IsFixed(time_j) ? FixedValue(time_j) : time_j.offset());
1975   // In all formulas, active_i and active_j are symmetrical, we can sort the
1976   // active literals.
1977   if (active_j < active_i) std::swap(active_i, active_j);
1978   return std::make_tuple(var_i, coeff_i, var_j, coeff_j, offset, active_i,
1979                          active_j);
1980 }
1981 
ClearPrecedenceCache()1982 void PresolveContext::ClearPrecedenceCache() {
1983   reified_precedences_cache_.clear();
1984 }
1985 
LogInfo()1986 void PresolveContext::LogInfo() {
1987   SOLVER_LOG(logger_, "");
1988   SOLVER_LOG(logger_, "Presolve summary:");
1989   SOLVER_LOG(logger_, "  - ", NumAffineRelations(),
1990              " affine relations were detected.");
1991   SOLVER_LOG(logger_, "  - ", NumEquivRelations(),
1992              " variable equivalence relations were detected.");
1993   std::map<std::string, int> sorted_rules(stats_by_rule_name_.begin(),
1994                                           stats_by_rule_name_.end());
1995   for (const auto& entry : sorted_rules) {
1996     if (entry.second == 1) {
1997       SOLVER_LOG(logger_, "  - rule '", entry.first, "' was applied 1 time.");
1998     } else {
1999       SOLVER_LOG(logger_, "  - rule '", entry.first, "' was applied ",
2000                  entry.second, " times.");
2001     }
2002   }
2003 }
2004 
LoadModelForProbing(PresolveContext * context,Model * local_model)2005 bool LoadModelForProbing(PresolveContext* context, Model* local_model) {
2006   if (context->ModelIsUnsat()) return false;
2007 
2008   // Update the domain in the current CpModelProto.
2009   context->WriteVariableDomainsToProto();
2010   const CpModelProto& model_proto = *(context->working_model);
2011 
2012   // Load the constraints in a local model.
2013   //
2014   // TODO(user): The model we load does not contain affine relations! But
2015   // ideally we should be able to remove all of them once we allow more complex
2016   // constraints to contains linear expression.
2017   //
2018   // TODO(user): remove code duplication with cp_model_solver. Here we also do
2019   // not run the heuristic to decide which variable to fully encode.
2020   //
2021   // TODO(user): Maybe do not load slow to propagate constraints? for instance
2022   // we do not use any linear relaxation here.
2023   Model model;
2024   local_model->Register<SolverLogger>(context->logger());
2025 
2026   // Adapt some of the parameters during this probing phase.
2027   auto* local_param = local_model->GetOrCreate<SatParameters>();
2028   *local_param = context->params();
2029   local_param->set_use_implied_bounds(false);
2030 
2031   local_model->GetOrCreate<TimeLimit>()->MergeWithGlobalTimeLimit(
2032       context->time_limit());
2033   local_model->Register<ModelRandomGenerator>(context->random());
2034   auto* encoder = local_model->GetOrCreate<IntegerEncoder>();
2035   encoder->DisableImplicationBetweenLiteral();
2036   auto* mapping = local_model->GetOrCreate<CpModelMapping>();
2037 
2038   // Important: Because the model_proto do not contains affine relation or the
2039   // objective, we cannot call DetectOptionalVariables() ! This might wrongly
2040   // detect optionality and derive bad conclusion.
2041   LoadVariables(model_proto, /*view_all_booleans_as_integers=*/false,
2042                 local_model);
2043   ExtractEncoding(model_proto, local_model);
2044   auto* sat_solver = local_model->GetOrCreate<SatSolver>();
2045   for (const ConstraintProto& ct : model_proto.constraints()) {
2046     if (mapping->ConstraintIsAlreadyLoaded(&ct)) continue;
2047     CHECK(LoadConstraint(ct, local_model));
2048     if (sat_solver->IsModelUnsat()) {
2049       return context->NotifyThatModelIsUnsat(absl::StrCat(
2050           "after loading constraint during probing ", ct.ShortDebugString()));
2051     }
2052   }
2053   encoder->AddAllImplicationsBetweenAssociatedLiterals();
2054   if (!sat_solver->Propagate()) {
2055     return context->NotifyThatModelIsUnsat(
2056         "during probing initial propagation");
2057   }
2058 
2059   return true;
2060 }
2061 
2062 }  // namespace sat
2063 }  // namespace operations_research
2064