OR-Tools  8.2
deviation.cc
Go to the documentation of this file.
1// Copyright 2010-2018 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 <algorithm>
15#include <cstdlib>
16#include <memory>
17#include <string>
18#include <vector>
19
20#include "absl/strings/str_format.h"
26
27namespace operations_research {
28// Deviation Constraint, a constraint for the average absolute
29// deviation to the mean. See paper: Bound Consistent Deviation
30// Constraint, Pierre Schaus et. al., CP07
31namespace {
32class Deviation : public Constraint {
33 public:
34 Deviation(Solver* const solver, const std::vector<IntVar*>& vars,
35 IntVar* const deviation_var, int64 total_sum)
36 : Constraint(solver),
37 vars_(vars),
38 size_(vars.size()),
39 deviation_var_(deviation_var),
40 total_sum_(total_sum),
41 scaled_vars_assigned_value_(new int64[size_]),
42 scaled_vars_min_(new int64[size_]),
43 scaled_vars_max_(new int64[size_]),
44 scaled_sum_max_(0),
45 scaled_sum_min_(0),
46 maximum_(new int64[size_]),
47 overlaps_sup_(new int64[size_]),
48 active_sum_(0),
49 active_sum_rounded_down_(0),
50 active_sum_rounded_up_(0),
51 active_sum_nearest_(0) {
52 CHECK(deviation_var != nullptr);
53 }
54
55 ~Deviation() override {}
56
57 void Post() override {
58 Solver* const s = solver();
59 Demon* const demon = s->MakeConstraintInitialPropagateCallback(this);
60 for (int i = 0; i < size_; ++i) {
61 vars_[i]->WhenRange(demon);
62 }
63 deviation_var_->WhenRange(demon);
64 s->AddConstraint(s->MakeSumEquality(vars_, total_sum_));
65 }
66
67 void InitialPropagate() override {
68 const int64 delta_min = BuildMinimalDeviationAssignment();
69 deviation_var_->SetMin(delta_min);
70 PropagateBounds(delta_min);
71 }
72
73 std::string DebugString() const override {
74 return absl::StrFormat("Deviation([%s], deviation_var = %s, sum = %d)",
75 JoinDebugStringPtr(vars_, ", "),
76 deviation_var_->DebugString(), total_sum_);
77 }
78
79 void Accept(ModelVisitor* const visitor) const override {
80 visitor->BeginVisitConstraint(ModelVisitor::kDeviation, this);
81 visitor->VisitIntegerVariableArrayArgument(ModelVisitor::kVarsArgument,
82 vars_);
83 visitor->VisitIntegerExpressionArgument(ModelVisitor::kTargetArgument,
84 deviation_var_);
85 visitor->VisitIntegerArgument(ModelVisitor::kValueArgument, total_sum_);
86 visitor->EndVisitConstraint(ModelVisitor::kDeviation, this);
87 }
88
89 private:
90 // Builds an assignment with minimal deviation and assign it to
91 // scaled_vars_assigned_value_. It returns the minimal deviation:
92 // sum_i |scaled_vars_assigned_value_[i] - total_sum_|.
93 int64 BuildMinimalDeviationAssignment() {
94 RepairGreedySum(BuildGreedySum(true));
95 int64 minimal_deviation = 0;
96 for (int i = 0; i < size_; ++i) {
97 minimal_deviation +=
98 std::abs(scaled_vars_assigned_value_[i] - total_sum_);
99 }
100 return minimal_deviation;
101 }
102
103 // Propagates the upper and lower bounds of x[i]'s.
104 // It assumes the constraint is consistent:
105 // - the sum constraint is consistent
106 // - min deviation smaller than max allowed deviation
107 // min_delta is the minimum possible deviation
108 void PropagateBounds(int64 min_delta) {
109 PropagateBounds(min_delta, true); // Filter upper bounds.
110 PropagateBounds(min_delta, false); // Filter lower bounds.
111 }
112
113 // Prunes the upper/lower-bound of vars. We apply a mirroing of the
114 // domains wrt 0 to prune the lower bounds such that we can use the
115 // same algo to prune both sides of the domains. upperBounds = true
116 // to prune the upper bounds of vars, false to prune the lower
117 // bounds.
118 void PropagateBounds(int64 min_delta, bool upper_bound) {
119 // Builds greedy assignment.
120 const int64 greedy_sum = BuildGreedySum(upper_bound);
121 // Repairs assignment and store information to be used when pruning.
122 RepairSumAndComputeInfo(greedy_sum);
123 // Does the actual pruning.
124 PruneVars(min_delta, upper_bound);
125 }
126
127 // Cache min and max values of variables.
128 void ComputeData(bool upper_bound) {
129 scaled_sum_max_ = 0;
130 scaled_sum_min_ = 0;
131 for (int i = 0; i < size_; ++i) {
132 scaled_vars_max_[i] =
133 size_ * (upper_bound ? vars_[i]->Max() : -vars_[i]->Min());
134 scaled_vars_min_[i] =
135 size_ * (upper_bound ? vars_[i]->Min() : -vars_[i]->Max());
136 scaled_sum_max_ += scaled_vars_max_[i];
137 scaled_sum_min_ += scaled_vars_min_[i];
138 }
139 active_sum_ = (!upper_bound ? -total_sum_ : total_sum_);
140 // down is <= sum.
141 active_sum_rounded_down_ =
142 size_ * MathUtil::FloorOfRatio<int64>(active_sum_, size_);
143 // up is > sum, always.
144 active_sum_rounded_up_ = active_sum_rounded_down_ + size_;
145 active_sum_nearest_ = (active_sum_rounded_up_ - active_sum_ <=
146 active_sum_ - active_sum_rounded_down_)
147 ? active_sum_rounded_up_
148 : active_sum_rounded_down_;
149 }
150
151 // Builds an approximate sum in a greedy way.
152 int64 BuildGreedySum(bool upper_bound) {
153 // Update data structure.
154 ComputeData(upper_bound);
155
156 // Number of constraint should be consistent.
157 DCHECK_GE(size_ * active_sum_, scaled_sum_min_);
158 DCHECK_LE(size_ * active_sum_, scaled_sum_max_);
159
160 int64 sum = 0;
161 // Greedily assign variable to nearest value to average.
162 overlaps_.clear();
163 for (int i = 0; i < size_; ++i) {
164 if (scaled_vars_min_[i] >= active_sum_) {
165 scaled_vars_assigned_value_[i] = scaled_vars_min_[i];
166 } else if (scaled_vars_max_[i] <= active_sum_) {
167 scaled_vars_assigned_value_[i] = scaled_vars_max_[i];
168 } else {
169 // Overlapping variable scaled_vars_min_[i] < active_sum_ <
170 // scaled_vars_max_[i].
171 scaled_vars_assigned_value_[i] = active_sum_nearest_;
172 if (active_sum_ % size_ != 0) {
173 overlaps_.push_back(i);
174 }
175 }
176 sum += scaled_vars_assigned_value_[i];
177 }
178 DCHECK_EQ(0, active_sum_rounded_down_ % size_);
179 DCHECK_LE(active_sum_rounded_down_, active_sum_);
180 DCHECK_LT(active_sum_ - active_sum_rounded_down_, size_);
181
182 return sum;
183 }
184
185 bool Overlap(int var_index) const {
186 return scaled_vars_min_[var_index] < active_sum_ &&
187 scaled_vars_max_[var_index] > active_sum_;
188 }
189
190 // Repairs the greedy sum obtained above to get the correct sum.
191 void RepairGreedySum(int64 greedy_sum) {
192 // Useful constant: scaled version of the sum.
193 const int64 scaled_total_sum = size_ * active_sum_;
194 // Step used to make the repair.
195 const int64 delta = greedy_sum > scaled_total_sum ? -size_ : size_;
196
197 // Change overlapping variables as long as the sum is not
198 // satisfied and there are overlapping vars, we use that ones to
199 // repair.
200 for (int j = 0; j < overlaps_.size() && greedy_sum != scaled_total_sum;
201 j++) {
202 scaled_vars_assigned_value_[overlaps_[j]] += delta;
203 greedy_sum += delta;
204 }
205 // Change other variables if the sum is still not satisfied.
206 for (int i = 0; i < size_ && greedy_sum != scaled_total_sum; ++i) {
207 const int64 old_scaled_vars_i = scaled_vars_assigned_value_[i];
208 if (greedy_sum < scaled_total_sum) {
209 // Increase scaled_vars_assigned_value_[i] as much as
210 // possible to fix the too low sum.
211 scaled_vars_assigned_value_[i] += scaled_total_sum - greedy_sum;
212 scaled_vars_assigned_value_[i] =
213 std::min(scaled_vars_assigned_value_[i], scaled_vars_max_[i]);
214 } else {
215 // Decrease scaled_vars_assigned_value_[i] as much as
216 // possible to fix the too high sum.
217 scaled_vars_assigned_value_[i] -= (greedy_sum - scaled_total_sum);
218 scaled_vars_assigned_value_[i] =
219 std::max(scaled_vars_assigned_value_[i], scaled_vars_min_[i]);
220 }
221 // Maintain the sum.
222 greedy_sum += scaled_vars_assigned_value_[i] - old_scaled_vars_i;
223 }
224 }
225
226 // Computes the maximum values of variables in the case the repaired
227 // greedy sum is actually the active sum.
228 void ComputeMaxWhenNoRepair() {
229 int num_overlap_sum_rounded_up = 0;
230 if (active_sum_nearest_ == active_sum_rounded_up_) {
231 num_overlap_sum_rounded_up = overlaps_.size();
232 }
233 for (int i = 0; i < size_; ++i) {
234 maximum_[i] = scaled_vars_assigned_value_[i];
235 if (Overlap(i) && active_sum_nearest_ == active_sum_rounded_up_ &&
236 active_sum_ % size_ != 0) {
237 overlaps_sup_[i] = num_overlap_sum_rounded_up - 1;
238 } else {
239 overlaps_sup_[i] = num_overlap_sum_rounded_up;
240 }
241 }
242 }
243
244 // Returns the number of variables overlapping the average value,
245 // assigned to // the average value rounded up that we can/need to
246 // move.
247 int ComputeNumOverlapsVariableRoundedUp() {
248 if (active_sum_ % size_ == 0) {
249 return 0;
250 }
251 int num_overlap_sum_rounded_up = 0;
252 for (int i = 0; i < size_; ++i) {
253 if (scaled_vars_assigned_value_[i] > scaled_vars_min_[i] &&
254 scaled_vars_assigned_value_[i] == active_sum_rounded_up_) {
255 num_overlap_sum_rounded_up++;
256 }
257 }
258 return num_overlap_sum_rounded_up;
259 }
260
261 // Returns whether we can push the greedy sum across the scaled
262 // total sum in the same direction as going from the nearest rounded
263 // sum to the farthest one.
264 bool CanPushSumAcrossMean(int64 greedy_sum, int64 scaled_total_sum) const {
265 return (greedy_sum > scaled_total_sum &&
266 active_sum_nearest_ == active_sum_rounded_up_) ||
267 (greedy_sum < scaled_total_sum &&
268 active_sum_nearest_ == active_sum_rounded_down_);
269 }
270
271 // Repairs the sum and store intermediate information to be used
272 // during pruning.
273 void RepairSumAndComputeInfo(int64 greedy_sum) {
274 const int64 scaled_total_sum = size_ * active_sum_;
275 // Computation of key values for the pruning:
276 // - overlaps_sup_
277 // - maximum_[i]
278 if (greedy_sum == scaled_total_sum) { // No repair needed.
279 ComputeMaxWhenNoRepair();
280 } else { // Repair and compute maximums.
281 // Try to repair the sum greedily.
282 if (CanPushSumAcrossMean(greedy_sum, scaled_total_sum)) {
283 const int64 delta = greedy_sum > scaled_total_sum ? -size_ : size_;
284 for (int j = 0; j < overlaps_.size() && greedy_sum != scaled_total_sum;
285 ++j) {
286 scaled_vars_assigned_value_[overlaps_[j]] += delta;
287 greedy_sum += delta;
288 }
289 }
290
291 const int num_overlap_sum_rounded_up =
292 ComputeNumOverlapsVariableRoundedUp();
293
294 if (greedy_sum == scaled_total_sum) {
295 // Greedy sum is repaired.
296 for (int i = 0; i < size_; ++i) {
297 if (Overlap(i) && num_overlap_sum_rounded_up > 0) {
298 maximum_[i] = active_sum_rounded_up_;
299 overlaps_sup_[i] = num_overlap_sum_rounded_up - 1;
300 } else {
301 maximum_[i] = scaled_vars_assigned_value_[i];
302 overlaps_sup_[i] = num_overlap_sum_rounded_up;
303 }
304 }
305 } else if (greedy_sum > scaled_total_sum) {
306 // scaled_vars_assigned_value_[i] = active_sum_rounded_down_ or
307 // scaled_vars_assigned_value_[i] <= total_sum
308 // (there is no more num_overlap_sum_rounded_up).
309 for (int i = 0; i < size_; ++i) {
310 maximum_[i] = scaled_vars_assigned_value_[i];
311 overlaps_sup_[i] = 0;
312 }
313 } else { // greedy_sum < scaled_total_sum.
314 for (int i = 0; i < size_; ++i) {
315 if (Overlap(i) && num_overlap_sum_rounded_up > 0) {
316 overlaps_sup_[i] = num_overlap_sum_rounded_up - 1;
317 } else {
318 overlaps_sup_[i] = num_overlap_sum_rounded_up;
319 }
320
321 if (scaled_vars_assigned_value_[i] < scaled_vars_max_[i]) {
322 maximum_[i] =
323 scaled_vars_assigned_value_[i] + scaled_total_sum - greedy_sum;
324 } else {
325 maximum_[i] = scaled_vars_assigned_value_[i];
326 }
327 }
328 }
329 }
330 }
331
332 // Propagates onto variables with all computed data.
333 void PruneVars(int64 min_delta, bool upper_bound) {
334 // Pruning of upper bound of vars_[i] for var_index in [1..n].
335 const int64 increase_down_up = (active_sum_rounded_up_ - active_sum_) -
336 (active_sum_ - active_sum_rounded_down_);
337 for (int var_index = 0; var_index < size_; ++var_index) {
338 // Not bound, and a compatible new max.
339 if (scaled_vars_max_[var_index] != scaled_vars_min_[var_index] &&
340 maximum_[var_index] < scaled_vars_max_[var_index]) {
341 const int64 new_max =
342 ComputeNewMax(var_index, min_delta, increase_down_up);
343 PruneBound(var_index, new_max, upper_bound);
344 }
345 }
346 }
347
348 // Computes new max for a variable.
349 int64 ComputeNewMax(int var_index, int64 min_delta, int64 increase_down_up) {
350 int64 maximum_value = maximum_[var_index];
351 int64 current_min_delta = min_delta;
352
353 if (overlaps_sup_[var_index] > 0 &&
354 (current_min_delta +
355 overlaps_sup_[var_index] * (size_ - increase_down_up) >=
356 deviation_var_->Max())) {
357 const int64 delta = deviation_var_->Max() - current_min_delta;
358 maximum_value += (size_ * delta) / (size_ - increase_down_up);
359 return MathUtil::FloorOfRatio<int64>(maximum_value, size_);
360 } else {
361 if (maximum_value == active_sum_rounded_down_ &&
362 active_sum_rounded_down_ < active_sum_) {
363 DCHECK_EQ(0, overlaps_sup_[var_index]);
364 current_min_delta += size_ + increase_down_up;
365 if (current_min_delta > deviation_var_->Max()) {
366 DCHECK_EQ(0, maximum_value % size_);
367 return maximum_value / size_;
368 }
369 maximum_value += size_;
370 }
371 current_min_delta +=
372 overlaps_sup_[var_index] * (size_ - increase_down_up);
373 maximum_value += size_ * overlaps_sup_[var_index];
374 // Slope of 2 x n.
375 const int64 delta = deviation_var_->Max() - current_min_delta;
376 maximum_value += delta / 2; // n * delta / (2 * n);
377 return MathUtil::FloorOfRatio<int64>(maximum_value, size_);
378 }
379 }
380
381 // Sets maximum on var or on its opposite.
382 void PruneBound(int var_index, int64 bound, bool upper_bound) {
383 if (upper_bound) {
384 vars_[var_index]->SetMax(bound);
385 } else {
386 vars_[var_index]->SetMin(-bound);
387 }
388 }
389
390 std::vector<IntVar*> vars_;
391 const int size_;
392 IntVar* const deviation_var_;
393 const int64 total_sum_;
394 std::unique_ptr<int64[]> scaled_vars_assigned_value_;
395 std::unique_ptr<int64[]> scaled_vars_min_;
396 std::unique_ptr<int64[]> scaled_vars_max_;
397 int64 scaled_sum_max_;
398 int64 scaled_sum_min_;
399 // Stores the variables overlapping the mean value.
400 std::vector<int> overlaps_;
401 std::unique_ptr<int64[]> maximum_;
402 std::unique_ptr<int64[]> overlaps_sup_;
403 // These values are updated by ComputeData().
404 int64 active_sum_;
405 int64 active_sum_rounded_down_;
406 int64 active_sum_rounded_up_;
407 int64 active_sum_nearest_;
408};
409} // namespace
410
411Constraint* Solver::MakeDeviation(const std::vector<IntVar*>& vars,
412 IntVar* const deviation_var,
413 int64 total_sum) {
414 return RevAlloc(new Deviation(this, vars, deviation_var, total_sum));
415}
416} // namespace operations_research
int64 min
Definition: alldiff_cst.cc:138
const std::vector< IntVar * > vars_
Definition: alldiff_cst.cc:43
int64 max
Definition: alldiff_cst.cc:139
#define CHECK(condition)
Definition: base/logging.h:495
#define DCHECK_LE(val1, val2)
Definition: base/logging.h:887
#define DCHECK_GE(val1, val2)
Definition: base/logging.h:889
#define DCHECK_LT(val1, val2)
Definition: base/logging.h:888
#define DCHECK_EQ(val1, val2)
Definition: base/logging.h:885
A constraint is the main modeling object.
The class IntVar is a subset of IntExpr.
Constraint * MakeDeviation(const std::vector< IntVar * > &vars, IntVar *const deviation_var, int64 total_sum)
Deviation constraint: sum_i |n * vars[i] - total_sum| <= deviation_var and sum_i vars[i] == total_sum...
Definition: deviation.cc:411
T * RevAlloc(T *object)
Registers the given object as being reversible.
int64_t int64
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
std::string JoinDebugStringPtr(const std::vector< T > &v, const std::string &separator)
Definition: string_array.h:45
int64 delta
Definition: resource.cc:1684
int64 bound