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/cp_model_lns.h"
15
16 #include <algorithm>
17 #include <cstdint>
18 #include <limits>
19 #include <numeric>
20 #include <vector>
21
22 #include "absl/container/flat_hash_set.h"
23 #include "absl/strings/str_join.h"
24 #include "absl/synchronization/mutex.h"
25 #include "ortools/graph/connected_components.h"
26 #include "ortools/sat/cp_model.pb.h"
27 #include "ortools/sat/cp_model_mapping.h"
28 #include "ortools/sat/cp_model_utils.h"
29 #include "ortools/sat/integer.h"
30 #include "ortools/sat/linear_programming_constraint.h"
31 #include "ortools/sat/presolve_context.h"
32 #include "ortools/sat/rins.h"
33 #include "ortools/sat/synchronization.h"
34 #include "ortools/util/saturated_arithmetic.h"
35
36 namespace operations_research {
37 namespace sat {
38
NeighborhoodGeneratorHelper(CpModelProto const * model_proto,SatParameters const * parameters,SharedResponseManager * shared_response,SharedTimeLimit * shared_time_limit,SharedBoundsManager * shared_bounds)39 NeighborhoodGeneratorHelper::NeighborhoodGeneratorHelper(
40 CpModelProto const* model_proto, SatParameters const* parameters,
41 SharedResponseManager* shared_response, SharedTimeLimit* shared_time_limit,
42 SharedBoundsManager* shared_bounds)
43 : SubSolver(""),
44 parameters_(*parameters),
45 model_proto_(*model_proto),
46 shared_time_limit_(shared_time_limit),
47 shared_bounds_(shared_bounds),
48 shared_response_(shared_response) {
49 CHECK(shared_response_ != nullptr);
50 if (shared_bounds_ != nullptr) {
51 shared_bounds_id_ = shared_bounds_->RegisterNewId();
52 }
53 *model_proto_with_only_variables_.mutable_variables() =
54 model_proto_.variables();
55 InitializeHelperData();
56 RecomputeHelperData();
57 Synchronize();
58 }
59
Synchronize()60 void NeighborhoodGeneratorHelper::Synchronize() {
61 if (shared_bounds_ != nullptr) {
62 std::vector<int> model_variables;
63 std::vector<int64_t> new_lower_bounds;
64 std::vector<int64_t> new_upper_bounds;
65 shared_bounds_->GetChangedBounds(shared_bounds_id_, &model_variables,
66 &new_lower_bounds, &new_upper_bounds);
67
68 bool new_variables_have_been_fixed = false;
69
70 {
71 absl::MutexLock domain_lock(&domain_mutex_);
72
73 for (int i = 0; i < model_variables.size(); ++i) {
74 const int var = model_variables[i];
75 const int64_t new_lb = new_lower_bounds[i];
76 const int64_t new_ub = new_upper_bounds[i];
77 if (VLOG_IS_ON(3)) {
78 const auto& domain =
79 model_proto_with_only_variables_.variables(var).domain();
80 const int64_t old_lb = domain.Get(0);
81 const int64_t old_ub = domain.Get(domain.size() - 1);
82 VLOG(3) << "Variable: " << var << " old domain: [" << old_lb << ", "
83 << old_ub << "] new domain: [" << new_lb << ", " << new_ub
84 << "]";
85 }
86 const Domain old_domain = ReadDomainFromProto(
87 model_proto_with_only_variables_.variables(var));
88 const Domain new_domain =
89 old_domain.IntersectionWith(Domain(new_lb, new_ub));
90 if (new_domain.IsEmpty()) {
91 // This can mean two things:
92 // 1/ This variable is a normal one and the problem is UNSAT or
93 // 2/ This variable is optional, and its associated literal must be
94 // set to false.
95 //
96 // Currently, we wait for any full solver to pick the crossing bounds
97 // and do the correct stuff on their own. We do not want to have empty
98 // domain in the proto as this would means INFEASIBLE. So we just
99 // ignore such bounds here.
100 //
101 // TODO(user): We could set the optional literal to false directly in
102 // the bound sharing manager. We do have to be careful that all the
103 // different solvers have the same optionality definition though.
104 continue;
105 }
106 FillDomainInProto(
107 new_domain,
108 model_proto_with_only_variables_.mutable_variables(var));
109 new_variables_have_been_fixed |= new_domain.IsFixed();
110 }
111 }
112
113 // Only trigger the computation if needed.
114 if (new_variables_have_been_fixed) {
115 RecomputeHelperData();
116 }
117 }
118 }
119
ObjectiveDomainIsConstraining() const120 bool NeighborhoodGeneratorHelper::ObjectiveDomainIsConstraining() const {
121 if (!model_proto_.has_objective()) return false;
122 if (model_proto_.objective().domain().empty()) return false;
123
124 int64_t min_activity = 0;
125 int64_t max_activity = 0;
126 const int num_terms = model_proto_.objective().vars().size();
127 for (int i = 0; i < num_terms; ++i) {
128 const int var = PositiveRef(model_proto_.objective().vars(i));
129 const int64_t coeff = model_proto_.objective().coeffs(i);
130 const auto& var_domain =
131 model_proto_with_only_variables_.variables(var).domain();
132 const int64_t v1 = coeff * var_domain[0];
133 const int64_t v2 = coeff * var_domain[var_domain.size() - 1];
134 min_activity += std::min(v1, v2);
135 max_activity += std::max(v1, v2);
136 }
137
138 const Domain obj_domain = ReadDomainFromProto(model_proto_.objective());
139 const Domain inferred_domain =
140 Domain(min_activity, max_activity)
141 .IntersectionWith(
142 Domain(std::numeric_limits<int64_t>::min(), obj_domain.Max()));
143 return !inferred_domain.IsIncludedIn(obj_domain);
144 }
145
InitializeHelperData()146 void NeighborhoodGeneratorHelper::InitializeHelperData() {
147 type_to_constraints_.clear();
148 const int num_constraints = model_proto_.constraints_size();
149 for (int c = 0; c < num_constraints; ++c) {
150 const int type = model_proto_.constraints(c).constraint_case();
151 if (type >= type_to_constraints_.size()) {
152 type_to_constraints_.resize(type + 1);
153 }
154 type_to_constraints_[type].push_back(c);
155 }
156
157 const int num_variables = model_proto_.variables().size();
158 is_in_objective_.resize(num_variables, false);
159 if (model_proto_.has_objective()) {
160 for (const int ref : model_proto_.objective().vars()) {
161 is_in_objective_[PositiveRef(ref)] = true;
162 }
163 }
164 }
165
166 // Recompute all the data when new variables have been fixed. Note that this
167 // shouldn't be called if there is no change as it is in O(problem size).
RecomputeHelperData()168 void NeighborhoodGeneratorHelper::RecomputeHelperData() {
169 absl::MutexLock graph_lock(&graph_mutex_);
170 absl::ReaderMutexLock domain_lock(&domain_mutex_);
171
172 // Do basic presolving to have a more precise graph.
173 // Here we just remove trivially true constraints.
174 //
175 // Note(user): We do that each time a new variable is fixed. It might be too
176 // much, but on the miplib and in 1200s, we do that only about 1k time on the
177 // worst case problem.
178 //
179 // TODO(user): Change API to avoid a few copy?
180 // TODO(user): We could keep the context in the class.
181 // TODO(user): We can also start from the previous simplified model instead.
182 {
183 Model local_model;
184 CpModelProto mapping_proto;
185 simplied_model_proto_.Clear();
186 *simplied_model_proto_.mutable_variables() =
187 model_proto_with_only_variables_.variables();
188 PresolveContext context(&local_model, &simplied_model_proto_,
189 &mapping_proto);
190 ModelCopy copier(&context);
191
192 // TODO(user): Not sure what to do if the model is UNSAT.
193 // This shouldn't matter as it should be dealt with elsewhere.
194 copier.ImportAndSimplifyConstraints(model_proto_, {});
195 }
196
197 // Compute the constraint <-> variable graph.
198 //
199 // TODO(user): Remove duplicate constraints?
200 const auto& constraints = simplied_model_proto_.constraints();
201 var_to_constraint_.assign(model_proto_.variables_size(), {});
202 constraint_to_var_.assign(constraints.size(), {});
203 int reduced_ct_index = 0;
204 for (int ct_index = 0; ct_index < constraints.size(); ++ct_index) {
205 // We remove the interval constraints since we should have an equivalent
206 // linear constraint somewhere else.
207 if (constraints[ct_index].constraint_case() == ConstraintProto::kInterval) {
208 continue;
209 }
210
211 for (const int var : UsedVariables(constraints[ct_index])) {
212 if (IsConstant(var)) continue;
213 constraint_to_var_[reduced_ct_index].push_back(var);
214 }
215
216 // We replace intervals by their underlying integer variables. Note that
217 // this is needed for a correct decomposition into independent part.
218 for (const int interval : UsedIntervals(constraints[ct_index])) {
219 for (const int var : UsedVariables(constraints[interval])) {
220 if (IsConstant(var)) continue;
221 constraint_to_var_[reduced_ct_index].push_back(var);
222 }
223 }
224
225 // We remove constraint of size 0 and 1 since they are not useful for LNS
226 // based on this graph.
227 if (constraint_to_var_[reduced_ct_index].size() <= 1) {
228 constraint_to_var_[reduced_ct_index].clear();
229 continue;
230 }
231
232 // Keep this constraint.
233 for (const int var : constraint_to_var_[reduced_ct_index]) {
234 var_to_constraint_[var].push_back(reduced_ct_index);
235 }
236 ++reduced_ct_index;
237 }
238 constraint_to_var_.resize(reduced_ct_index);
239
240 // We mark as active all non-constant variables.
241 // Non-active variable will never be fixed in standard LNS fragment.
242 active_variables_.clear();
243 const int num_variables = model_proto_.variables_size();
244 active_variables_set_.assign(num_variables, false);
245 for (int i = 0; i < num_variables; ++i) {
246 if (!IsConstant(i)) {
247 active_variables_.push_back(i);
248 active_variables_set_[i] = true;
249 }
250 }
251
252 // Compute connected components.
253 // Note that fixed variable are just ignored.
254 DenseConnectedComponentsFinder union_find;
255 union_find.SetNumberOfNodes(num_variables);
256 for (const std::vector<int>& var_in_constraint : constraint_to_var_) {
257 if (var_in_constraint.size() <= 1) continue;
258 for (int i = 1; i < var_in_constraint.size(); ++i) {
259 union_find.AddEdge(var_in_constraint[0], var_in_constraint[i]);
260 }
261 }
262
263 // If we have a lower bound on the objective, then this "objective constraint"
264 // might link components together.
265 if (ObjectiveDomainIsConstraining()) {
266 const auto& refs = model_proto_.objective().vars();
267 const int num_terms = refs.size();
268 for (int i = 1; i < num_terms; ++i) {
269 union_find.AddEdge(PositiveRef(refs[0]), PositiveRef(refs[i]));
270 }
271 }
272
273 // Compute all components involving non-fixed variables.
274 //
275 // TODO(user): If a component has no objective, we can fix it to any feasible
276 // solution. This will automatically be done by LNS fragment covering such
277 // component though.
278 components_.clear();
279 var_to_component_index_.assign(num_variables, -1);
280 for (int var = 0; var < num_variables; ++var) {
281 if (IsConstant(var)) continue;
282 const int root = union_find.FindRoot(var);
283 CHECK_LT(root, var_to_component_index_.size());
284 int& index = var_to_component_index_[root];
285 if (index == -1) {
286 index = components_.size();
287 components_.push_back({});
288 }
289 var_to_component_index_[var] = index;
290 components_[index].push_back(var);
291 }
292
293 // Display information about the reduced problem.
294 //
295 // TODO(user): Exploit connected component while generating fragments.
296 // TODO(user): Do not generate fragment not touching the objective.
297 std::vector<int> component_sizes;
298 for (const std::vector<int>& component : components_) {
299 component_sizes.push_back(component.size());
300 }
301 std::sort(component_sizes.begin(), component_sizes.end(),
302 std::greater<int>());
303 std::string compo_message;
304 if (component_sizes.size() > 1) {
305 if (component_sizes.size() <= 10) {
306 compo_message =
307 absl::StrCat(" compo:", absl::StrJoin(component_sizes, ","));
308 } else {
309 component_sizes.resize(10);
310 compo_message =
311 absl::StrCat(" compo:", absl::StrJoin(component_sizes, ","), ",...");
312 }
313 }
314 shared_response_->LogMessage(
315 absl::StrCat("var:", active_variables_.size(), "/", num_variables,
316 " constraints:", simplied_model_proto_.constraints().size(),
317 "/", model_proto_.constraints().size(), compo_message));
318 }
319
IsActive(int var) const320 bool NeighborhoodGeneratorHelper::IsActive(int var) const {
321 return active_variables_set_[var];
322 }
323
IsConstant(int var) const324 bool NeighborhoodGeneratorHelper::IsConstant(int var) const {
325 return model_proto_with_only_variables_.variables(var).domain_size() == 2 &&
326 model_proto_with_only_variables_.variables(var).domain(0) ==
327 model_proto_with_only_variables_.variables(var).domain(1);
328 }
329
FullNeighborhood() const330 Neighborhood NeighborhoodGeneratorHelper::FullNeighborhood() const {
331 Neighborhood neighborhood;
332 neighborhood.is_reduced = false;
333 neighborhood.is_generated = true;
334 {
335 absl::ReaderMutexLock lock(&domain_mutex_);
336 *neighborhood.delta.mutable_variables() =
337 model_proto_with_only_variables_.variables();
338 }
339 return neighborhood;
340 }
341
NoNeighborhood() const342 Neighborhood NeighborhoodGeneratorHelper::NoNeighborhood() const {
343 Neighborhood neighborhood;
344 neighborhood.is_generated = false;
345 return neighborhood;
346 }
347
GetActiveIntervals(const CpSolverResponse & initial_solution) const348 std::vector<int> NeighborhoodGeneratorHelper::GetActiveIntervals(
349 const CpSolverResponse& initial_solution) const {
350 std::vector<int> active_intervals;
351 absl::ReaderMutexLock lock(&domain_mutex_);
352 for (const int i : TypeToConstraints(ConstraintProto::kInterval)) {
353 const ConstraintProto& interval_ct = ModelProto().constraints(i);
354 // We only look at intervals that are performed in the solution. The
355 // unperformed intervals should be automatically freed during the generation
356 // phase.
357 if (interval_ct.enforcement_literal().size() == 1) {
358 const int enforcement_ref = interval_ct.enforcement_literal(0);
359 const int enforcement_var = PositiveRef(enforcement_ref);
360 const int value = initial_solution.solution(enforcement_var);
361 if (RefIsPositive(enforcement_ref) == (value == 0)) {
362 continue;
363 }
364 }
365
366 // We filter out fixed intervals. Because of presolve, if there is an
367 // enforcement literal, it cannot be fixed.
368 if (interval_ct.enforcement_literal().empty()) {
369 bool is_constant = true;
370 for (const int v : interval_ct.interval().start().vars()) {
371 if (!IsConstant(v)) {
372 is_constant = false;
373 break;
374 }
375 }
376 for (const int v : interval_ct.interval().size().vars()) {
377 if (!IsConstant(v)) {
378 is_constant = false;
379 break;
380 }
381 }
382 for (const int v : interval_ct.interval().end().vars()) {
383 if (!IsConstant(v)) {
384 is_constant = false;
385 break;
386 }
387 }
388 if (is_constant) continue;
389 }
390
391 active_intervals.push_back(i);
392 }
393 return active_intervals;
394 }
395
GetRoutingPaths(const CpSolverResponse & initial_solution) const396 std::vector<std::vector<int>> NeighborhoodGeneratorHelper::GetRoutingPaths(
397 const CpSolverResponse& initial_solution) const {
398 struct HeadAndArcLiteral {
399 int head;
400 int literal;
401 };
402
403 std::vector<std::vector<int>> result;
404 absl::flat_hash_map<int, HeadAndArcLiteral> tail_to_head_and_arc_literal;
405
406 for (const int i : TypeToConstraints(ConstraintProto::kCircuit)) {
407 const CircuitConstraintProto& ct = ModelProto().constraints(i).circuit();
408
409 // Collect arcs.
410 int min_node = std::numeric_limits<int>::max();
411 tail_to_head_and_arc_literal.clear();
412 for (int i = 0; i < ct.literals_size(); ++i) {
413 const int literal = ct.literals(i);
414 const int head = ct.heads(i);
415 const int tail = ct.tails(i);
416 const int bool_var = PositiveRef(literal);
417 const int64_t value = initial_solution.solution(bool_var);
418 // Skip unselected arcs.
419 if (RefIsPositive(literal) == (value == 0)) continue;
420 // Ignore self loops.
421 if (head == tail) continue;
422 tail_to_head_and_arc_literal[tail] = {head, bool_var};
423 min_node = std::min(tail, min_node);
424 }
425 if (tail_to_head_and_arc_literal.empty()) continue;
426
427 // Unroll the path.
428 int current_node = min_node;
429 std::vector<int> path;
430 do {
431 auto it = tail_to_head_and_arc_literal.find(current_node);
432 CHECK(it != tail_to_head_and_arc_literal.end());
433 current_node = it->second.head;
434 path.push_back(it->second.literal);
435 } while (current_node != min_node);
436 result.push_back(std::move(path));
437 }
438
439 std::vector<HeadAndArcLiteral> route_starts;
440 for (const int i : TypeToConstraints(ConstraintProto::kRoutes)) {
441 const RoutesConstraintProto& ct = ModelProto().constraints(i).routes();
442 tail_to_head_and_arc_literal.clear();
443 route_starts.clear();
444
445 // Collect route starts and arcs.
446 for (int i = 0; i < ct.literals_size(); ++i) {
447 const int literal = ct.literals(i);
448 const int head = ct.heads(i);
449 const int tail = ct.tails(i);
450 const int bool_var = PositiveRef(literal);
451 const int64_t value = initial_solution.solution(bool_var);
452 // Skip unselected arcs.
453 if (RefIsPositive(literal) == (value == 0)) continue;
454 // Ignore self loops.
455 if (head == tail) continue;
456 if (tail == 0) {
457 route_starts.push_back({head, bool_var});
458 } else {
459 tail_to_head_and_arc_literal[tail] = {head, bool_var};
460 }
461 }
462
463 // Unroll all routes.
464 for (const HeadAndArcLiteral& head_var : route_starts) {
465 std::vector<int> path;
466 int current_node = head_var.head;
467 path.push_back(head_var.literal);
468 do {
469 auto it = tail_to_head_and_arc_literal.find(current_node);
470 CHECK(it != tail_to_head_and_arc_literal.end());
471 current_node = it->second.head;
472 path.push_back(it->second.literal);
473 } while (current_node != 0);
474 result.push_back(std::move(path));
475 }
476 }
477
478 return result;
479 }
480
FixGivenVariables(const CpSolverResponse & base_solution,const absl::flat_hash_set<int> & variables_to_fix) const481 Neighborhood NeighborhoodGeneratorHelper::FixGivenVariables(
482 const CpSolverResponse& base_solution,
483 const absl::flat_hash_set<int>& variables_to_fix) const {
484 Neighborhood neighborhood;
485
486 // Fill in neighborhood.delta all variable domains.
487 {
488 absl::ReaderMutexLock domain_lock(&domain_mutex_);
489
490 const int num_variables =
491 model_proto_with_only_variables_.variables().size();
492 neighborhood.delta.mutable_variables()->Reserve(num_variables);
493 for (int i = 0; i < num_variables; ++i) {
494 const IntegerVariableProto& current_var =
495 model_proto_with_only_variables_.variables(i);
496 IntegerVariableProto* new_var = neighborhood.delta.add_variables();
497
498 // We only copy the name in debug mode.
499 if (DEBUG_MODE) new_var->set_name(current_var.name());
500
501 const Domain domain = ReadDomainFromProto(current_var);
502 const int64_t base_value = base_solution.solution(i);
503
504 // It seems better to always start from a feasible point, so if the base
505 // solution is no longer valid under the new up to date bound, we make
506 // sure to relax the domain so that it is.
507 if (!domain.Contains(base_value)) {
508 // TODO(user): this can happen when variables_to_fix.contains(i). But we
509 // should probably never consider as "active" such variable in the first
510 // place.
511 //
512 // TODO(user): This might lead to incompatibility when the base solution
513 // is not compatible with variable we fixed in a connected component! so
514 // maybe it is not great to do that. Initial investigation didn't see
515 // a big change. More work needed.
516 FillDomainInProto(domain.UnionWith(Domain(base_solution.solution(i))),
517 new_var);
518 } else if (variables_to_fix.contains(i)) {
519 new_var->add_domain(base_value);
520 new_var->add_domain(base_value);
521 } else {
522 FillDomainInProto(domain, new_var);
523 }
524 }
525 }
526
527 // Fill some statistic fields and detect if we cover a full component.
528 //
529 // TODO(user): If there is just one component, we can skip some computation.
530 {
531 absl::ReaderMutexLock graph_lock(&graph_mutex_);
532 std::vector<int> count(components_.size(), 0);
533 const int num_variables = neighborhood.delta.variables().size();
534 for (int var = 0; var < num_variables; ++var) {
535 const auto& domain = neighborhood.delta.variables(var).domain();
536 if (domain.size() != 2 || domain[0] != domain[1]) {
537 ++neighborhood.num_relaxed_variables;
538 if (is_in_objective_[var]) {
539 ++neighborhood.num_relaxed_variables_in_objective;
540 }
541 const int c = var_to_component_index_[var];
542 if (c != -1) count[c]++;
543 }
544 }
545
546 for (int i = 0; i < components_.size(); ++i) {
547 if (count[i] == components_[i].size()) {
548 neighborhood.variables_that_can_be_fixed_to_local_optimum.insert(
549 neighborhood.variables_that_can_be_fixed_to_local_optimum.end(),
550 components_[i].begin(), components_[i].end());
551 }
552 }
553 }
554
555 // If the objective domain might cut the optimal solution, we cannot exploit
556 // the connected components. We compute this outside the mutex to avoid
557 // any deadlock risk.
558 //
559 // TODO(user): We could handle some complex domain (size > 2).
560 if (model_proto_.has_objective() &&
561 (model_proto_.objective().domain().size() != 2 ||
562 shared_response_->GetInnerObjectiveLowerBound() <
563 model_proto_.objective().domain(0))) {
564 neighborhood.variables_that_can_be_fixed_to_local_optimum.clear();
565 }
566
567 AddSolutionHinting(base_solution, &neighborhood.delta);
568
569 neighborhood.is_generated = true;
570 neighborhood.is_reduced = !variables_to_fix.empty();
571 neighborhood.is_simple = true;
572
573 // TODO(user): force better objective? Note that this is already done when the
574 // hint above is successfully loaded (i.e. if it passes the presolve
575 // correctly) since the solver will try to find better solution than the
576 // current one.
577 return neighborhood;
578 }
579
AddSolutionHinting(const CpSolverResponse & initial_solution,CpModelProto * model_proto) const580 void NeighborhoodGeneratorHelper::AddSolutionHinting(
581 const CpSolverResponse& initial_solution, CpModelProto* model_proto) const {
582 // Set the current solution as a hint.
583 model_proto->clear_solution_hint();
584 const auto is_fixed = [model_proto](int var) {
585 const IntegerVariableProto& var_proto = model_proto->variables(var);
586 return var_proto.domain_size() == 2 &&
587 var_proto.domain(0) == var_proto.domain(1);
588 };
589 for (int var = 0; var < model_proto->variables_size(); ++var) {
590 if (is_fixed(var)) continue;
591
592 model_proto->mutable_solution_hint()->add_vars(var);
593 model_proto->mutable_solution_hint()->add_values(
594 initial_solution.solution(var));
595 }
596 }
597
RemoveMarkedConstraints(const std::vector<int> & constraints_to_remove) const598 Neighborhood NeighborhoodGeneratorHelper::RemoveMarkedConstraints(
599 const std::vector<int>& constraints_to_remove) const {
600 Neighborhood neighborhood = FullNeighborhood();
601
602 if (constraints_to_remove.empty()) return neighborhood;
603 neighborhood.is_reduced = false;
604 neighborhood.constraints_to_ignore = constraints_to_remove;
605 return neighborhood;
606 }
607
RelaxGivenVariables(const CpSolverResponse & initial_solution,const std::vector<int> & relaxed_variables) const608 Neighborhood NeighborhoodGeneratorHelper::RelaxGivenVariables(
609 const CpSolverResponse& initial_solution,
610 const std::vector<int>& relaxed_variables) const {
611 std::vector<bool> relaxed_variables_set(model_proto_.variables_size(), false);
612 for (const int var : relaxed_variables) relaxed_variables_set[var] = true;
613 absl::flat_hash_set<int> fixed_variables;
614 {
615 absl::ReaderMutexLock graph_lock(&graph_mutex_);
616 for (const int i : active_variables_) {
617 if (!relaxed_variables_set[i]) {
618 fixed_variables.insert(i);
619 }
620 }
621 }
622 return FixGivenVariables(initial_solution, fixed_variables);
623 }
624
FixAllVariables(const CpSolverResponse & initial_solution) const625 Neighborhood NeighborhoodGeneratorHelper::FixAllVariables(
626 const CpSolverResponse& initial_solution) const {
627 const std::vector<int>& all_variables = ActiveVariables();
628 const absl::flat_hash_set<int> fixed_variables(all_variables.begin(),
629 all_variables.end());
630 return FixGivenVariables(initial_solution, fixed_variables);
631 }
632
ReadyToGenerate() const633 bool NeighborhoodGenerator::ReadyToGenerate() const {
634 return (helper_.shared_response().SolutionsRepository().NumSolutions() > 0);
635 }
636
GetUCBScore(int64_t total_num_calls) const637 double NeighborhoodGenerator::GetUCBScore(int64_t total_num_calls) const {
638 absl::ReaderMutexLock mutex_lock(&generator_mutex_);
639 DCHECK_GE(total_num_calls, num_calls_);
640 if (num_calls_ <= 10) return std::numeric_limits<double>::infinity();
641 return current_average_ + sqrt((2 * log(total_num_calls)) / num_calls_);
642 }
643
Synchronize()644 void NeighborhoodGenerator::Synchronize() {
645 absl::MutexLock mutex_lock(&generator_mutex_);
646
647 // To make the whole update process deterministic, we currently sort the
648 // SolveData.
649 std::sort(solve_data_.begin(), solve_data_.end());
650
651 // This will be used to update the difficulty of this neighborhood.
652 int num_fully_solved_in_batch = 0;
653 int num_not_fully_solved_in_batch = 0;
654
655 for (const SolveData& data : solve_data_) {
656 AdditionalProcessingOnSynchronize(data);
657 ++num_calls_;
658
659 // INFEASIBLE or OPTIMAL means that we "fully solved" the local problem.
660 // If we didn't, then we cannot be sure that there is no improving solution
661 // in that neighborhood.
662 if (data.status == CpSolverStatus::INFEASIBLE ||
663 data.status == CpSolverStatus::OPTIMAL) {
664 ++num_fully_solved_calls_;
665 ++num_fully_solved_in_batch;
666 } else {
667 ++num_not_fully_solved_in_batch;
668 }
669
670 // It seems to make more sense to compare the new objective to the base
671 // solution objective, not the best one. However this causes issue in the
672 // logic below because on some problems the neighborhood can always lead
673 // to a better "new objective" if the base solution wasn't the best one.
674 //
675 // This might not be a final solution, but it does work ok for now.
676 const IntegerValue best_objective_improvement =
677 IsRelaxationGenerator()
678 ? IntegerValue(CapSub(data.new_objective_bound.value(),
679 data.initial_best_objective_bound.value()))
680 : IntegerValue(CapSub(data.initial_best_objective.value(),
681 data.new_objective.value()));
682 if (best_objective_improvement > 0) {
683 num_consecutive_non_improving_calls_ = 0;
684 } else {
685 ++num_consecutive_non_improving_calls_;
686 }
687
688 // TODO(user): Weight more recent data.
689 // degrade the current average to forget old learnings.
690 const double gain_per_time_unit =
691 std::max(0.0, static_cast<double>(best_objective_improvement.value())) /
692 (1.0 + data.deterministic_time);
693 if (num_calls_ <= 100) {
694 current_average_ += (gain_per_time_unit - current_average_) / num_calls_;
695 } else {
696 current_average_ = 0.9 * current_average_ + 0.1 * gain_per_time_unit;
697 }
698
699 deterministic_time_ += data.deterministic_time;
700 }
701
702 // Update the difficulty.
703 difficulty_.Update(/*num_decreases=*/num_not_fully_solved_in_batch,
704 /*num_increases=*/num_fully_solved_in_batch);
705
706 // Bump the time limit if we saw no better solution in the last few calls.
707 // This means that as the search progress, we likely spend more and more time
708 // trying to solve individual neighborhood.
709 //
710 // TODO(user): experiment with resetting the time limit if a solution is
711 // found.
712 if (num_consecutive_non_improving_calls_ > 50) {
713 num_consecutive_non_improving_calls_ = 0;
714 deterministic_limit_ *= 1.02;
715
716 // We do not want the limit to go to high. Intuitively, the goal is to try
717 // out a lot of neighborhoods, not just spend a lot of time on a few.
718 deterministic_limit_ = std::min(60.0, deterministic_limit_);
719 }
720
721 solve_data_.clear();
722 }
723
724 namespace {
725
GetRandomSubset(double relative_size,std::vector<int> * base,absl::BitGenRef random)726 void GetRandomSubset(double relative_size, std::vector<int>* base,
727 absl::BitGenRef random) {
728 if (base->empty()) return;
729
730 // TODO(user): we could generate this more efficiently than using random
731 // shuffle.
732 std::shuffle(base->begin(), base->end(), random);
733 const int target_size = std::round(relative_size * base->size());
734 base->resize(target_size);
735 }
736
737 } // namespace
738
Generate(const CpSolverResponse & initial_solution,double difficulty,absl::BitGenRef random)739 Neighborhood RelaxRandomVariablesGenerator::Generate(
740 const CpSolverResponse& initial_solution, double difficulty,
741 absl::BitGenRef random) {
742 std::vector<int> fixed_variables = helper_.ActiveVariables();
743 GetRandomSubset(1.0 - difficulty, &fixed_variables, random);
744 return helper_.FixGivenVariables(
745 initial_solution, {fixed_variables.begin(), fixed_variables.end()});
746 }
747
Generate(const CpSolverResponse & initial_solution,double difficulty,absl::BitGenRef random)748 Neighborhood RelaxRandomConstraintsGenerator::Generate(
749 const CpSolverResponse& initial_solution, double difficulty,
750 absl::BitGenRef random) {
751 if (helper_.DifficultyMeansFullNeighborhood(difficulty)) {
752 return helper_.FullNeighborhood();
753 }
754
755 std::vector<int> relaxed_variables;
756 {
757 absl::ReaderMutexLock graph_lock(&helper_.graph_mutex_);
758 const int num_active_constraints = helper_.ConstraintToVar().size();
759 std::vector<int> active_constraints(num_active_constraints);
760 for (int c = 0; c < num_active_constraints; ++c) {
761 active_constraints[c] = c;
762 }
763 std::shuffle(active_constraints.begin(), active_constraints.end(), random);
764
765 const int num_model_vars = helper_.ModelProto().variables_size();
766 std::vector<bool> visited_variables_set(num_model_vars, false);
767
768 const int num_active_vars =
769 helper_.ActiveVariablesWhileHoldingLock().size();
770 const int target_size = std::ceil(difficulty * num_active_vars);
771 CHECK_GT(target_size, 0);
772
773 for (const int constraint_index : active_constraints) {
774 for (const int var : helper_.ConstraintToVar()[constraint_index]) {
775 if (visited_variables_set[var]) continue;
776 visited_variables_set[var] = true;
777 if (helper_.IsActive(var)) {
778 relaxed_variables.push_back(var);
779 if (relaxed_variables.size() == target_size) break;
780 }
781 }
782 if (relaxed_variables.size() == target_size) break;
783 }
784 }
785
786 return helper_.RelaxGivenVariables(initial_solution, relaxed_variables);
787 }
788
Generate(const CpSolverResponse & initial_solution,double difficulty,absl::BitGenRef random)789 Neighborhood VariableGraphNeighborhoodGenerator::Generate(
790 const CpSolverResponse& initial_solution, double difficulty,
791 absl::BitGenRef random) {
792 if (helper_.DifficultyMeansFullNeighborhood(difficulty)) {
793 return helper_.FullNeighborhood();
794 }
795
796 const int num_model_vars = helper_.ModelProto().variables_size();
797 std::vector<bool> visited_variables_set(num_model_vars, false);
798 std::vector<int> relaxed_variables;
799 std::vector<int> visited_variables;
800
801 // It is important complexity wise to never scan a constraint twice!
802 const int num_model_constraints = helper_.ModelProto().constraints_size();
803 std::vector<bool> scanned_constraints(num_model_constraints, false);
804
805 std::vector<int> random_variables;
806 {
807 absl::ReaderMutexLock graph_lock(&helper_.graph_mutex_);
808
809 // The number of active variables can decrease asynchronously.
810 // We read the exact number while locked.
811 const int num_active_vars =
812 helper_.ActiveVariablesWhileHoldingLock().size();
813 const int target_size = std::ceil(difficulty * num_active_vars);
814 CHECK_GT(target_size, 0) << difficulty << " " << num_active_vars;
815
816 const int first_var =
817 helper_.ActiveVariablesWhileHoldingLock()[absl::Uniform<int>(
818 random, 0, num_active_vars)];
819
820 visited_variables_set[first_var] = true;
821 visited_variables.push_back(first_var);
822 relaxed_variables.push_back(first_var);
823
824 for (int i = 0; i < visited_variables.size(); ++i) {
825 random_variables.clear();
826 // Collect all the variables that appears in the same constraints as
827 // visited_variables[i].
828 for (const int ct : helper_.VarToConstraint()[visited_variables[i]]) {
829 if (scanned_constraints[ct]) continue;
830 scanned_constraints[ct] = true;
831 for (const int var : helper_.ConstraintToVar()[ct]) {
832 if (visited_variables_set[var]) continue;
833 visited_variables_set[var] = true;
834 random_variables.push_back(var);
835 }
836 }
837 // We always randomize to change the partial subgraph explored
838 // afterwards.
839 std::shuffle(random_variables.begin(), random_variables.end(), random);
840 for (const int var : random_variables) {
841 if (relaxed_variables.size() < target_size) {
842 visited_variables.push_back(var);
843 if (helper_.IsActive(var)) {
844 relaxed_variables.push_back(var);
845 }
846 } else {
847 break;
848 }
849 }
850 if (relaxed_variables.size() >= target_size) break;
851 }
852 }
853 return helper_.RelaxGivenVariables(initial_solution, relaxed_variables);
854 }
855
Generate(const CpSolverResponse & initial_solution,double difficulty,absl::BitGenRef random)856 Neighborhood ConstraintGraphNeighborhoodGenerator::Generate(
857 const CpSolverResponse& initial_solution, double difficulty,
858 absl::BitGenRef random) {
859 const int num_model_constraints = helper_.ModelProto().constraints_size();
860 if (num_model_constraints == 0 ||
861 helper_.DifficultyMeansFullNeighborhood(difficulty)) {
862 return helper_.FullNeighborhood();
863 }
864
865 const int num_model_vars = helper_.ModelProto().variables_size();
866 std::vector<bool> visited_variables_set(num_model_vars, false);
867 std::vector<int> relaxed_variables;
868
869 std::vector<bool> added_constraints(num_model_constraints, false);
870 std::vector<int> next_constraints;
871
872 std::vector<int> random_variables;
873 {
874 absl::ReaderMutexLock graph_lock(&helper_.graph_mutex_);
875 const int num_active_vars =
876 helper_.ActiveVariablesWhileHoldingLock().size();
877 const int target_size = std::ceil(difficulty * num_active_vars);
878 CHECK_GT(target_size, 0);
879
880 // Start by a random constraint.
881 const int num_active_constraints = helper_.ConstraintToVar().size();
882 if (num_active_constraints != 0) {
883 next_constraints.push_back(
884 absl::Uniform<int>(random, 0, num_active_constraints));
885 added_constraints[next_constraints.back()] = true;
886 }
887
888 while (relaxed_variables.size() < target_size) {
889 // Stop if we have a full connected component.
890 if (next_constraints.empty()) break;
891
892 // Pick a random unprocessed constraint.
893 const int i = absl::Uniform<int>(random, 0, next_constraints.size());
894 const int constraint_index = next_constraints[i];
895 std::swap(next_constraints[i], next_constraints.back());
896 next_constraints.pop_back();
897
898 // Add all the variable of this constraint and increase the set of next
899 // possible constraints.
900 CHECK_LT(constraint_index, num_active_constraints);
901 random_variables = helper_.ConstraintToVar()[constraint_index];
902 std::shuffle(random_variables.begin(), random_variables.end(), random);
903 for (const int var : random_variables) {
904 if (visited_variables_set[var]) continue;
905 visited_variables_set[var] = true;
906 if (helper_.IsActive(var)) {
907 relaxed_variables.push_back(var);
908 }
909 if (relaxed_variables.size() == target_size) break;
910
911 for (const int ct : helper_.VarToConstraint()[var]) {
912 if (added_constraints[ct]) continue;
913 added_constraints[ct] = true;
914 next_constraints.push_back(ct);
915 }
916 }
917 }
918 }
919 return helper_.RelaxGivenVariables(initial_solution, relaxed_variables);
920 }
921
922 namespace {
923
GetLinearExpressionValue(const LinearExpressionProto & expr,const CpSolverResponse & initial_solution)924 int64_t GetLinearExpressionValue(const LinearExpressionProto& expr,
925 const CpSolverResponse& initial_solution) {
926 int64_t result = expr.offset();
927 for (int i = 0; i < expr.vars_size(); ++i) {
928 result += expr.coeffs(i) * initial_solution.solution(expr.vars(i));
929 }
930 return result;
931 }
932
AddLinearExpressionToConstraint(const int64_t coeff,const LinearExpressionProto & expr,LinearConstraintProto * constraint,int64_t * rhs_offset)933 void AddLinearExpressionToConstraint(const int64_t coeff,
934 const LinearExpressionProto& expr,
935 LinearConstraintProto* constraint,
936 int64_t* rhs_offset) {
937 *rhs_offset -= coeff * expr.offset();
938 for (int i = 0; i < expr.vars_size(); ++i) {
939 constraint->add_vars(expr.vars(i));
940 constraint->add_coeffs(expr.coeffs(i) * coeff);
941 }
942 }
943
AddPrecedenceConstraints(const absl::Span<const int> intervals,const absl::flat_hash_set<int> & ignored_intervals,const CpSolverResponse & initial_solution,const NeighborhoodGeneratorHelper & helper,Neighborhood * neighborhood)944 void AddPrecedenceConstraints(const absl::Span<const int> intervals,
945 const absl::flat_hash_set<int>& ignored_intervals,
946 const CpSolverResponse& initial_solution,
947 const NeighborhoodGeneratorHelper& helper,
948 Neighborhood* neighborhood) {
949 // Sort all non-relaxed intervals of this constraint by current start
950 // time.
951 std::vector<std::pair<int64_t, int>> start_interval_pairs;
952 for (const int i : intervals) {
953 if (ignored_intervals.contains(i)) continue;
954 const ConstraintProto& interval_ct = helper.ModelProto().constraints(i);
955
956 // TODO(user): we ignore size zero for now.
957 const LinearExpressionProto& size_var = interval_ct.interval().size();
958 if (GetLinearExpressionValue(size_var, initial_solution) == 0) continue;
959
960 const LinearExpressionProto& start_var = interval_ct.interval().start();
961 const int64_t start_value =
962 GetLinearExpressionValue(start_var, initial_solution);
963
964 start_interval_pairs.push_back({start_value, i});
965 }
966 std::sort(start_interval_pairs.begin(), start_interval_pairs.end());
967
968 // Add precedence between the remaining intervals, forcing their order.
969 for (int i = 0; i + 1 < start_interval_pairs.size(); ++i) {
970 const LinearExpressionProto& before_start =
971 helper.ModelProto()
972 .constraints(start_interval_pairs[i].second)
973 .interval()
974 .start();
975 const LinearExpressionProto& before_end =
976 helper.ModelProto()
977 .constraints(start_interval_pairs[i].second)
978 .interval()
979 .end();
980 const LinearExpressionProto& after_start =
981 helper.ModelProto()
982 .constraints(start_interval_pairs[i + 1].second)
983 .interval()
984 .start();
985
986 // If the end was smaller we keep it that way, otherwise we just order the
987 // start variables.
988 LinearConstraintProto* linear =
989 neighborhood->delta.add_constraints()->mutable_linear();
990 linear->add_domain(std::numeric_limits<int64_t>::min());
991 int64_t rhs_offset = 0;
992 if (GetLinearExpressionValue(before_end, initial_solution) <=
993 GetLinearExpressionValue(after_start, initial_solution)) {
994 // If the end was smaller than the next start, keep it that way.
995 AddLinearExpressionToConstraint(1, before_end, linear, &rhs_offset);
996 } else {
997 // Otherwise, keep the same minimum separation. This is done in order
998 // to "simplify" the neighborhood.
999 rhs_offset = GetLinearExpressionValue(before_start, initial_solution) -
1000 GetLinearExpressionValue(after_start, initial_solution);
1001 AddLinearExpressionToConstraint(1, before_start, linear, &rhs_offset);
1002 }
1003
1004 AddLinearExpressionToConstraint(-1, after_start, linear, &rhs_offset);
1005 linear->add_domain(rhs_offset);
1006
1007 // The linear constraint should be satisfied by the current solution.
1008 if (DEBUG_MODE) {
1009 int64_t activity = 0;
1010 for (int i = 0; i < linear->vars().size(); ++i) {
1011 activity +=
1012 linear->coeffs(i) * initial_solution.solution(linear->vars(i));
1013 }
1014 CHECK_GE(activity, linear->domain(0));
1015 CHECK_LE(activity, linear->domain(1));
1016 }
1017 }
1018 }
1019 } // namespace
1020
GenerateSchedulingNeighborhoodForRelaxation(const absl::Span<const int> intervals_to_relax,const CpSolverResponse & initial_solution,const NeighborhoodGeneratorHelper & helper)1021 Neighborhood GenerateSchedulingNeighborhoodForRelaxation(
1022 const absl::Span<const int> intervals_to_relax,
1023 const CpSolverResponse& initial_solution,
1024 const NeighborhoodGeneratorHelper& helper) {
1025 Neighborhood neighborhood = helper.FullNeighborhood();
1026 neighborhood.is_reduced =
1027 (intervals_to_relax.size() <
1028 helper.TypeToConstraints(ConstraintProto::kInterval).size());
1029
1030 // We will extend the set with some interval that we cannot fix.
1031 absl::flat_hash_set<int> ignored_intervals(intervals_to_relax.begin(),
1032 intervals_to_relax.end());
1033
1034 // Fix the presence/absence of non-relaxed intervals.
1035 for (const int i : helper.TypeToConstraints(ConstraintProto::kInterval)) {
1036 DCHECK_GE(i, 0);
1037 if (ignored_intervals.contains(i)) continue;
1038
1039 const ConstraintProto& interval_ct = helper.ModelProto().constraints(i);
1040 if (interval_ct.enforcement_literal().empty()) continue;
1041
1042 CHECK_EQ(interval_ct.enforcement_literal().size(), 1);
1043 const int enforcement_ref = interval_ct.enforcement_literal(0);
1044 const int enforcement_var = PositiveRef(enforcement_ref);
1045 const int value = initial_solution.solution(enforcement_var);
1046
1047 // If the interval is not enforced, we just relax it. If it belongs to an
1048 // exactly one constraint, and the enforced interval is not relaxed, then
1049 // propagation will force this interval to stay not enforced. Otherwise,
1050 // LNS will be able to change which interval will be enforced among all
1051 // alternatives.
1052 if (RefIsPositive(enforcement_ref) == (value == 0)) {
1053 ignored_intervals.insert(i);
1054 continue;
1055 }
1056
1057 // Fix the value.
1058 neighborhood.delta.mutable_variables(enforcement_var)->clear_domain();
1059 neighborhood.delta.mutable_variables(enforcement_var)->add_domain(value);
1060 neighborhood.delta.mutable_variables(enforcement_var)->add_domain(value);
1061 }
1062
1063 for (const int c : helper.TypeToConstraints(ConstraintProto::kNoOverlap)) {
1064 AddPrecedenceConstraints(
1065 helper.ModelProto().constraints(c).no_overlap().intervals(),
1066 ignored_intervals, initial_solution, helper, &neighborhood);
1067 }
1068 for (const int c : helper.TypeToConstraints(ConstraintProto::kCumulative)) {
1069 AddPrecedenceConstraints(
1070 helper.ModelProto().constraints(c).cumulative().intervals(),
1071 ignored_intervals, initial_solution, helper, &neighborhood);
1072 }
1073 for (const int c : helper.TypeToConstraints(ConstraintProto::kNoOverlap2D)) {
1074 AddPrecedenceConstraints(
1075 helper.ModelProto().constraints(c).no_overlap_2d().x_intervals(),
1076 ignored_intervals, initial_solution, helper, &neighborhood);
1077 AddPrecedenceConstraints(
1078 helper.ModelProto().constraints(c).no_overlap_2d().y_intervals(),
1079 ignored_intervals, initial_solution, helper, &neighborhood);
1080 }
1081
1082 // Set the current solution as a hint.
1083 helper.AddSolutionHinting(initial_solution, &neighborhood.delta);
1084 neighborhood.is_generated = true;
1085
1086 return neighborhood;
1087 }
1088
Generate(const CpSolverResponse & initial_solution,double difficulty,absl::BitGenRef random)1089 Neighborhood SchedulingNeighborhoodGenerator::Generate(
1090 const CpSolverResponse& initial_solution, double difficulty,
1091 absl::BitGenRef random) {
1092 std::vector<int> intervals_to_relax =
1093 helper_.GetActiveIntervals(initial_solution);
1094 GetRandomSubset(difficulty, &intervals_to_relax, random);
1095
1096 return GenerateSchedulingNeighborhoodForRelaxation(intervals_to_relax,
1097 initial_solution, helper_);
1098 }
1099
Generate(const CpSolverResponse & initial_solution,double difficulty,absl::BitGenRef random)1100 Neighborhood SchedulingTimeWindowNeighborhoodGenerator::Generate(
1101 const CpSolverResponse& initial_solution, double difficulty,
1102 absl::BitGenRef random) {
1103 std::vector<std::pair<int64_t, int>> start_interval_pairs;
1104 const std::vector<int> active_intervals =
1105 helper_.GetActiveIntervals(initial_solution);
1106 std::vector<int> intervals_to_relax;
1107
1108 if (active_intervals.empty()) return helper_.FullNeighborhood();
1109
1110 for (const int i : active_intervals) {
1111 const ConstraintProto& interval_ct = helper_.ModelProto().constraints(i);
1112 const LinearExpressionProto& start_var = interval_ct.interval().start();
1113 const int64_t start_value =
1114 GetLinearExpressionValue(start_var, initial_solution);
1115 start_interval_pairs.push_back({start_value, i});
1116 }
1117 std::sort(start_interval_pairs.begin(), start_interval_pairs.end());
1118 const int relaxed_size = std::floor(difficulty * start_interval_pairs.size());
1119
1120 std::uniform_int_distribution<int> random_var(
1121 0, start_interval_pairs.size() - relaxed_size - 1);
1122 const int random_start_index = random_var(random);
1123
1124 // TODO(user): Consider relaxing more than one time window
1125 // intervals. This seems to help with Giza models.
1126 for (int i = random_start_index; i < relaxed_size; ++i) {
1127 intervals_to_relax.push_back(start_interval_pairs[i].second);
1128 }
1129
1130 return GenerateSchedulingNeighborhoodForRelaxation(intervals_to_relax,
1131 initial_solution, helper_);
1132 }
1133
Generate(const CpSolverResponse & initial_solution,double difficulty,absl::BitGenRef random)1134 Neighborhood RoutingRandomNeighborhoodGenerator::Generate(
1135 const CpSolverResponse& initial_solution, double difficulty,
1136 absl::BitGenRef random) {
1137 const std::vector<std::vector<int>> all_paths =
1138 helper_.GetRoutingPaths(initial_solution);
1139
1140 // Collect all unique variables.
1141 absl::flat_hash_set<int> all_path_variables;
1142 for (auto& path : all_paths) {
1143 all_path_variables.insert(path.begin(), path.end());
1144 }
1145 std::vector<int> fixed_variables(all_path_variables.begin(),
1146 all_path_variables.end());
1147 std::sort(fixed_variables.begin(), fixed_variables.end());
1148 GetRandomSubset(1.0 - difficulty, &fixed_variables, random);
1149 return helper_.FixGivenVariables(
1150 initial_solution, {fixed_variables.begin(), fixed_variables.end()});
1151 }
1152
Generate(const CpSolverResponse & initial_solution,double difficulty,absl::BitGenRef random)1153 Neighborhood RoutingPathNeighborhoodGenerator::Generate(
1154 const CpSolverResponse& initial_solution, double difficulty,
1155 absl::BitGenRef random) {
1156 std::vector<std::vector<int>> all_paths =
1157 helper_.GetRoutingPaths(initial_solution);
1158
1159 // Collect all unique variables.
1160 absl::flat_hash_set<int> all_path_variables;
1161 for (const auto& path : all_paths) {
1162 all_path_variables.insert(path.begin(), path.end());
1163 }
1164
1165 // Select variables to relax.
1166 const int num_variables_to_relax =
1167 static_cast<int>(all_path_variables.size() * difficulty);
1168 absl::flat_hash_set<int> relaxed_variables;
1169 while (relaxed_variables.size() < num_variables_to_relax) {
1170 DCHECK(!all_paths.empty());
1171 const int path_index = absl::Uniform<int>(random, 0, all_paths.size());
1172 std::vector<int>& path = all_paths[path_index];
1173 const int path_size = path.size();
1174 const int segment_length =
1175 std::min(path_size, absl::Uniform<int>(random, 4, 8));
1176 const int segment_start =
1177 absl::Uniform<int>(random, 0, path_size - segment_length);
1178 for (int i = segment_start; i < segment_start + segment_length; ++i) {
1179 relaxed_variables.insert(path[i]);
1180 }
1181
1182 // Remove segment and clean up empty paths.
1183 path.erase(path.begin() + segment_start,
1184 path.begin() + segment_start + segment_length);
1185 if (path.empty()) {
1186 std::swap(all_paths[path_index], all_paths.back());
1187 all_paths.pop_back();
1188 }
1189 }
1190
1191 // Compute the set of variables to fix.
1192 absl::flat_hash_set<int> fixed_variables;
1193 for (const int var : all_path_variables) {
1194 if (!relaxed_variables.contains(var)) fixed_variables.insert(var);
1195 }
1196 return helper_.FixGivenVariables(initial_solution, fixed_variables);
1197 }
1198
Generate(const CpSolverResponse & initial_solution,double difficulty,absl::BitGenRef random)1199 Neighborhood RoutingFullPathNeighborhoodGenerator::Generate(
1200 const CpSolverResponse& initial_solution, double difficulty,
1201 absl::BitGenRef random) {
1202 std::vector<std::vector<int>> all_paths =
1203 helper_.GetRoutingPaths(initial_solution);
1204 // Remove a corner case where all paths are empty.
1205 if (all_paths.empty()) {
1206 return helper_.NoNeighborhood();
1207 }
1208
1209 // Collect all unique variables.
1210 absl::flat_hash_set<int> all_path_variables;
1211 for (const auto& path : all_paths) {
1212 all_path_variables.insert(path.begin(), path.end());
1213 }
1214
1215 // Select variables to relax.
1216 const int num_variables_to_relax =
1217 static_cast<int>(all_path_variables.size() * difficulty);
1218 absl::flat_hash_set<int> relaxed_variables;
1219
1220 // Relax the start and end of each path to ease relocation.
1221 for (const auto& path : all_paths) {
1222 relaxed_variables.insert(path.front());
1223 relaxed_variables.insert(path.back());
1224 }
1225
1226 // Randomize paths.
1227 for (auto& path : all_paths) {
1228 std::shuffle(path.begin(), path.end(), random);
1229 }
1230
1231 // Relax all variables (if possible) in one random path.
1232 const int path_to_clean = absl::Uniform<int>(random, 0, all_paths.size());
1233 while (relaxed_variables.size() < num_variables_to_relax &&
1234 !all_paths[path_to_clean].empty()) {
1235 relaxed_variables.insert(all_paths[path_to_clean].back());
1236 all_paths[path_to_clean].pop_back();
1237 }
1238 if (all_paths[path_to_clean].empty()) {
1239 std::swap(all_paths[path_to_clean], all_paths.back());
1240 all_paths.pop_back();
1241 }
1242
1243 // Relax more variables until the target is reached.
1244 while (relaxed_variables.size() < num_variables_to_relax) {
1245 DCHECK(!all_paths.empty());
1246 const int path_index = absl::Uniform<int>(random, 0, all_paths.size());
1247 relaxed_variables.insert(all_paths[path_index].back());
1248
1249 // Remove variable and clean up empty paths.
1250 all_paths[path_index].pop_back();
1251 if (all_paths[path_index].empty()) {
1252 std::swap(all_paths[path_index], all_paths.back());
1253 all_paths.pop_back();
1254 }
1255 }
1256
1257 // Compute the set of variables to fix.
1258 absl::flat_hash_set<int> fixed_variables;
1259 for (const int var : all_path_variables) {
1260 if (!relaxed_variables.contains(var)) fixed_variables.insert(var);
1261 }
1262 return helper_.FixGivenVariables(initial_solution, fixed_variables);
1263 }
1264
ReadyToGenerate() const1265 bool RelaxationInducedNeighborhoodGenerator::ReadyToGenerate() const {
1266 if (incomplete_solutions_ != nullptr) {
1267 return incomplete_solutions_->HasNewSolution();
1268 }
1269
1270 if (response_manager_ != nullptr) {
1271 if (response_manager_->SolutionsRepository().NumSolutions() == 0) {
1272 return false;
1273 }
1274 }
1275
1276 // At least one relaxation solution should be available to generate a
1277 // neighborhood.
1278 if (lp_solutions_ != nullptr && lp_solutions_->NumSolutions() > 0) {
1279 return true;
1280 }
1281
1282 if (relaxation_solutions_ != nullptr &&
1283 relaxation_solutions_->NumSolutions() > 0) {
1284 return true;
1285 }
1286 return false;
1287 }
1288
Generate(const CpSolverResponse & initial_solution,double difficulty,absl::BitGenRef random)1289 Neighborhood RelaxationInducedNeighborhoodGenerator::Generate(
1290 const CpSolverResponse& initial_solution, double difficulty,
1291 absl::BitGenRef random) {
1292 Neighborhood neighborhood = helper_.FullNeighborhood();
1293 neighborhood.is_generated = false;
1294
1295 const bool lp_solution_available =
1296 (lp_solutions_ != nullptr && lp_solutions_->NumSolutions() > 0);
1297
1298 const bool relaxation_solution_available =
1299 (relaxation_solutions_ != nullptr &&
1300 relaxation_solutions_->NumSolutions() > 0);
1301
1302 const bool incomplete_solution_available =
1303 (incomplete_solutions_ != nullptr &&
1304 incomplete_solutions_->HasNewSolution());
1305
1306 if (!lp_solution_available && !relaxation_solution_available &&
1307 !incomplete_solution_available) {
1308 return neighborhood;
1309 }
1310
1311 RINSNeighborhood rins_neighborhood;
1312 // Randomly select the type of relaxation if both lp and relaxation solutions
1313 // are available.
1314 // TODO(user): Tune the probability value for this.
1315 std::bernoulli_distribution random_bool(0.5);
1316 const bool use_lp_relaxation =
1317 (lp_solution_available && relaxation_solution_available)
1318 ? random_bool(random)
1319 : lp_solution_available;
1320 if (use_lp_relaxation) {
1321 rins_neighborhood =
1322 GetRINSNeighborhood(response_manager_,
1323 /*relaxation_solutions=*/nullptr, lp_solutions_,
1324 incomplete_solutions_, random);
1325 neighborhood.source_info =
1326 incomplete_solution_available ? "incomplete" : "lp";
1327 } else {
1328 CHECK(relaxation_solution_available || incomplete_solution_available);
1329 rins_neighborhood = GetRINSNeighborhood(
1330 response_manager_, relaxation_solutions_,
1331 /*lp_solutions=*/nullptr, incomplete_solutions_, random);
1332 neighborhood.source_info =
1333 incomplete_solution_available ? "incomplete" : "relaxation";
1334 }
1335
1336 if (rins_neighborhood.fixed_vars.empty() &&
1337 rins_neighborhood.reduced_domain_vars.empty()) {
1338 return neighborhood;
1339 }
1340
1341 absl::ReaderMutexLock graph_lock(&helper_.graph_mutex_);
1342 // Fix the variables in the local model.
1343 for (const std::pair</*model_var*/ int, /*value*/ int64_t> fixed_var :
1344 rins_neighborhood.fixed_vars) {
1345 const int var = fixed_var.first;
1346 const int64_t value = fixed_var.second;
1347 if (var >= neighborhood.delta.variables_size()) continue;
1348 if (!helper_.IsActive(var)) continue;
1349
1350 if (!DomainInProtoContains(neighborhood.delta.variables(var), value)) {
1351 // TODO(user): Instead of aborting, pick the closest point in the domain?
1352 return neighborhood;
1353 }
1354
1355 neighborhood.delta.mutable_variables(var)->clear_domain();
1356 neighborhood.delta.mutable_variables(var)->add_domain(value);
1357 neighborhood.delta.mutable_variables(var)->add_domain(value);
1358 neighborhood.is_reduced = true;
1359 }
1360
1361 for (const std::pair</*model_var*/ int,
1362 /*domain*/ std::pair<int64_t, int64_t>>
1363 reduced_var : rins_neighborhood.reduced_domain_vars) {
1364 const int var = reduced_var.first;
1365 const int64_t lb = reduced_var.second.first;
1366 const int64_t ub = reduced_var.second.second;
1367 if (var >= neighborhood.delta.variables_size()) continue;
1368 if (!helper_.IsActive(var)) continue;
1369 Domain domain = ReadDomainFromProto(neighborhood.delta.variables(var));
1370 domain = domain.IntersectionWith(Domain(lb, ub));
1371 if (domain.IsEmpty()) {
1372 // TODO(user): Instead of aborting, pick the closest point in the domain?
1373 return neighborhood;
1374 }
1375 FillDomainInProto(domain, neighborhood.delta.mutable_variables(var));
1376 neighborhood.is_reduced = true;
1377 }
1378 neighborhood.is_generated = true;
1379 return neighborhood;
1380 }
1381
Generate(const CpSolverResponse & initial_solution,double difficulty,absl::BitGenRef random)1382 Neighborhood ConsecutiveConstraintsRelaxationNeighborhoodGenerator::Generate(
1383 const CpSolverResponse& initial_solution, double difficulty,
1384 absl::BitGenRef random) {
1385 std::vector<int> removable_constraints;
1386 const int num_constraints = helper_.ModelProto().constraints_size();
1387 removable_constraints.reserve(num_constraints);
1388 for (int c = 0; c < num_constraints; ++c) {
1389 // Removing intervals is not easy because other constraint might require
1390 // them, so for now, we don't remove them.
1391 if (helper_.ModelProto().constraints(c).constraint_case() ==
1392 ConstraintProto::kInterval) {
1393 continue;
1394 }
1395 removable_constraints.push_back(c);
1396 }
1397
1398 const int target_size =
1399 std::round((1.0 - difficulty) * removable_constraints.size());
1400
1401 const int random_start_index =
1402 absl::Uniform<int>(random, 0, removable_constraints.size());
1403 std::vector<int> removed_constraints;
1404 removed_constraints.reserve(target_size);
1405 int c = random_start_index;
1406 while (removed_constraints.size() < target_size) {
1407 removed_constraints.push_back(removable_constraints[c]);
1408 ++c;
1409 if (c == removable_constraints.size()) {
1410 c = 0;
1411 }
1412 }
1413
1414 return helper_.RemoveMarkedConstraints(removed_constraints);
1415 }
1416
1417 WeightedRandomRelaxationNeighborhoodGenerator::
WeightedRandomRelaxationNeighborhoodGenerator(NeighborhoodGeneratorHelper const * helper,const std::string & name)1418 WeightedRandomRelaxationNeighborhoodGenerator(
1419 NeighborhoodGeneratorHelper const* helper, const std::string& name)
1420 : NeighborhoodGenerator(name, helper) {
1421 std::vector<int> removable_constraints;
1422 const int num_constraints = helper_.ModelProto().constraints_size();
1423 constraint_weights_.reserve(num_constraints);
1424 // TODO(user): Experiment with different starting weights.
1425 for (int c = 0; c < num_constraints; ++c) {
1426 switch (helper_.ModelProto().constraints(c).constraint_case()) {
1427 case ConstraintProto::kCumulative:
1428 case ConstraintProto::kAllDiff:
1429 case ConstraintProto::kElement:
1430 case ConstraintProto::kRoutes:
1431 case ConstraintProto::kCircuit:
1432 constraint_weights_.push_back(3.0);
1433 num_removable_constraints_++;
1434 break;
1435 case ConstraintProto::kBoolOr:
1436 case ConstraintProto::kBoolAnd:
1437 case ConstraintProto::kBoolXor:
1438 case ConstraintProto::kIntProd:
1439 case ConstraintProto::kIntDiv:
1440 case ConstraintProto::kIntMod:
1441 case ConstraintProto::kLinMax:
1442 case ConstraintProto::kNoOverlap:
1443 case ConstraintProto::kNoOverlap2D:
1444 constraint_weights_.push_back(2.0);
1445 num_removable_constraints_++;
1446 break;
1447 case ConstraintProto::kLinear:
1448 case ConstraintProto::kTable:
1449 case ConstraintProto::kAutomaton:
1450 case ConstraintProto::kInverse:
1451 case ConstraintProto::kReservoir:
1452 case ConstraintProto::kAtMostOne:
1453 case ConstraintProto::kExactlyOne:
1454 constraint_weights_.push_back(1.0);
1455 num_removable_constraints_++;
1456 break;
1457 case ConstraintProto::CONSTRAINT_NOT_SET:
1458 case ConstraintProto::kDummyConstraint:
1459 case ConstraintProto::kInterval:
1460 // Removing intervals is not easy because other constraint might require
1461 // them, so for now, we don't remove them.
1462 constraint_weights_.push_back(0.0);
1463 break;
1464 }
1465 }
1466 }
1467
1468 void WeightedRandomRelaxationNeighborhoodGenerator::
AdditionalProcessingOnSynchronize(const SolveData & solve_data)1469 AdditionalProcessingOnSynchronize(const SolveData& solve_data) {
1470 const IntegerValue best_objective_improvement =
1471 solve_data.new_objective_bound - solve_data.initial_best_objective_bound;
1472
1473 const std::vector<int>& removed_constraints =
1474 removed_constraints_[solve_data.neighborhood_id];
1475
1476 // Heuristic: We change the weights of the removed constraints if the
1477 // neighborhood is solved (status is OPTIMAL or INFEASIBLE) or we observe an
1478 // improvement in objective bounds. Otherwise we assume that the
1479 // difficulty/time wasn't right for us to record feedbacks.
1480 //
1481 // If the objective bounds are improved, we bump up the weights. If the
1482 // objective bounds are worse and the problem status is OPTIMAL, we bump down
1483 // the weights. Otherwise if the new objective bounds are same as current
1484 // bounds (which happens a lot on some instances), we do not update the
1485 // weights as we do not have a clear signal whether the constraints removed
1486 // were good choices or not.
1487 // TODO(user): We can improve this heuristic with more experiments.
1488 if (best_objective_improvement > 0) {
1489 // Bump up the weights of all removed constraints.
1490 for (int c : removed_constraints) {
1491 if (constraint_weights_[c] <= 90.0) {
1492 constraint_weights_[c] += 10.0;
1493 } else {
1494 constraint_weights_[c] = 100.0;
1495 }
1496 }
1497 } else if (solve_data.status == CpSolverStatus::OPTIMAL &&
1498 best_objective_improvement < 0) {
1499 // Bump down the weights of all removed constraints.
1500 for (int c : removed_constraints) {
1501 if (constraint_weights_[c] > 0.5) {
1502 constraint_weights_[c] -= 0.5;
1503 }
1504 }
1505 }
1506 removed_constraints_.erase(solve_data.neighborhood_id);
1507 }
1508
Generate(const CpSolverResponse & initial_solution,double difficulty,absl::BitGenRef random)1509 Neighborhood WeightedRandomRelaxationNeighborhoodGenerator::Generate(
1510 const CpSolverResponse& initial_solution, double difficulty,
1511 absl::BitGenRef random) {
1512 const int target_size =
1513 std::round((1.0 - difficulty) * num_removable_constraints_);
1514
1515 std::vector<int> removed_constraints;
1516
1517 // Generate a random number between (0,1) = u[i] and use score[i] =
1518 // u[i]^(1/w[i]) and then select top k items with largest scores.
1519 // Reference: https://utopia.duth.gr/~pefraimi/research/data/2007EncOfAlg.pdf
1520 std::vector<std::pair<double, int>> constraint_removal_scores;
1521 std::uniform_real_distribution<double> random_var(0.0, 1.0);
1522 for (int c = 0; c < constraint_weights_.size(); ++c) {
1523 if (constraint_weights_[c] <= 0) continue;
1524 const double u = random_var(random);
1525 const double score = std::pow(u, (1 / constraint_weights_[c]));
1526 constraint_removal_scores.push_back({score, c});
1527 }
1528 std::sort(constraint_removal_scores.rbegin(),
1529 constraint_removal_scores.rend());
1530 for (int i = 0; i < target_size; ++i) {
1531 removed_constraints.push_back(constraint_removal_scores[i].second);
1532 }
1533
1534 Neighborhood result = helper_.RemoveMarkedConstraints(removed_constraints);
1535
1536 absl::MutexLock mutex_lock(&generator_mutex_);
1537 result.id = next_available_id_;
1538 next_available_id_++;
1539 removed_constraints_.insert({result.id, removed_constraints});
1540 return result;
1541 }
1542
1543 } // namespace sat
1544 } // namespace operations_research
1545