1

2

3

4

5

6

7

8

9

10

11

12

13

14#include <algorithm>

15#include <cstdint>

16#include <cstdlib>

17#include <memory>

18#include <string>

19#include <vector>

20

21#include "absl/strings/str_format.h"

27

29

30

31

32namespace {

34 public:

35 Deviation(Solver* const solver, const std::vector<IntVar*>& vars,

36 IntVar* const deviation_var, int64_t total_sum)

37 : Constraint(solver),

38 vars_(vars),

39 size_(vars.size()),

40 deviation_var_(deviation_var),

41 total_sum_(total_sum),

42 scaled_vars_assigned_value_(new int64_t[size_]),

43 scaled_vars_min_(new int64_t[size_]),

44 scaled_vars_max_(new int64_t[size_]),

45 scaled_sum_max_(0),

46 scaled_sum_min_(0),

47 maximum_(new int64_t[size_]),

48 overlaps_sup_(new int64_t[size_]),

49 active_sum_(0),

50 active_sum_rounded_down_(0),

51 active_sum_rounded_up_(0),

52 active_sum_nearest_(0) {

53 CHECK(deviation_var != nullptr);

54 }

55

56 ~Deviation() override {}

57

58 void Post() override {

59 Solver* const s = solver();

60 Demon* const demon = s->MakeConstraintInitialPropagateCallback(this);

61 for (int i = 0; i < size_; ++i) {

62 vars_[i]->WhenRange(demon);

63 }

64 deviation_var_->WhenRange(demon);

65 s->AddConstraint(s->MakeSumEquality(vars_, total_sum_));

66 }

67

68 void InitialPropagate() override {

69 const int64_t delta_min = BuildMinimalDeviationAssignment();

70 deviation_var_->SetMin(delta_min);

71 PropagateBounds(delta_min);

72 }

73

74 std::string DebugString() const override {

75 return absl::StrFormat("Deviation([%s], deviation_var = %s, sum = %d)",

77 deviation_var_->DebugString(), total_sum_);

78 }

79

80 void Accept(ModelVisitor* const visitor) const override {

83 vars_);

85 deviation_var_);

88 }

89

90 private:

91

92

93

94 int64_t BuildMinimalDeviationAssignment() {

95 RepairGreedySum(BuildGreedySum(true));

96 int64_t minimal_deviation = 0;

97 for (int i = 0; i < size_; ++i) {

98 minimal_deviation +=

99 std::abs(scaled_vars_assigned_value_[i] - total_sum_);

100 }

101 return minimal_deviation;

102 }

103

104

105

106

107

108

109 void PropagateBounds(int64_t min_delta) {

110 PropagateBounds(min_delta, true);

111 PropagateBounds(min_delta, false);

112 }

113

114

115

116

117

118

119 void PropagateBounds(int64_t min_delta, bool upper_bound) {

120

121 const int64_t greedy_sum = BuildGreedySum(upper_bound);

122

123 RepairSumAndComputeInfo(greedy_sum);

124

125 PruneVars(min_delta, upper_bound);

126 }

127

128

129 void ComputeData(bool upper_bound) {

130 scaled_sum_max_ = 0;

131 scaled_sum_min_ = 0;

132 for (int i = 0; i < size_; ++i) {

133 scaled_vars_max_[i] =

134 size_ * (upper_bound ? vars_[i]->Max() : -vars_[i]->Min());

135 scaled_vars_min_[i] =

136 size_ * (upper_bound ? vars_[i]->Min() : -vars_[i]->Max());

137 scaled_sum_max_ += scaled_vars_max_[i];

138 scaled_sum_min_ += scaled_vars_min_[i];

139 }

140 active_sum_ = (!upper_bound ? -total_sum_ : total_sum_);

141

142 active_sum_rounded_down_ =

144

145 active_sum_rounded_up_ = active_sum_rounded_down_ + size_;

146 active_sum_nearest_ = (active_sum_rounded_up_ - active_sum_ <=

147 active_sum_ - active_sum_rounded_down_)

148 ? active_sum_rounded_up_

149 : active_sum_rounded_down_;

150 }

151

152

153 int64_t BuildGreedySum(bool upper_bound) {

154

155 ComputeData(upper_bound);

156

157

158 DCHECK_GE(size_ * active_sum_, scaled_sum_min_);

159 DCHECK_LE(size_ * active_sum_, scaled_sum_max_);

160

161 int64_t sum = 0;

162

163 overlaps_.clear();

164 for (int i = 0; i < size_; ++i) {

165 if (scaled_vars_min_[i] >= active_sum_) {

166 scaled_vars_assigned_value_[i] = scaled_vars_min_[i];

167 } else if (scaled_vars_max_[i] <= active_sum_) {

168 scaled_vars_assigned_value_[i] = scaled_vars_max_[i];

169 } else {

170

171

172 scaled_vars_assigned_value_[i] = active_sum_nearest_;

173 if (active_sum_ % size_ != 0) {

174 overlaps_.push_back(i);

175 }

176 }

177 sum += scaled_vars_assigned_value_[i];

178 }

179 DCHECK_EQ(0, active_sum_rounded_down_ % size_);

180 DCHECK_LE(active_sum_rounded_down_, active_sum_);

181 DCHECK_LT(active_sum_ - active_sum_rounded_down_, size_);

182

183 return sum;

184 }

185

186 bool Overlap(int var_index) const {

187 return scaled_vars_min_[var_index] < active_sum_ &&

188 scaled_vars_max_[var_index] > active_sum_;

189 }

190

191

192 void RepairGreedySum(int64_t greedy_sum) {

193

194 const int64_t scaled_total_sum = size_ * active_sum_;

195

196 const int64_t delta = greedy_sum > scaled_total_sum ? -size_ : size_;

197

198

199

200

201 for (int j = 0; j < overlaps_.size() && greedy_sum != scaled_total_sum;

202 j++) {

203 scaled_vars_assigned_value_[overlaps_[j]] += delta;

204 greedy_sum += delta;

205 }

206

207 for (int i = 0; i < size_ && greedy_sum != scaled_total_sum; ++i) {

208 const int64_t old_scaled_vars_i = scaled_vars_assigned_value_[i];

209 if (greedy_sum < scaled_total_sum) {

210

211

212 scaled_vars_assigned_value_[i] += scaled_total_sum - greedy_sum;

213 scaled_vars_assigned_value_[i] =

214 std::min(scaled_vars_assigned_value_[i], scaled_vars_max_[i]);

215 } else {

216

217

218 scaled_vars_assigned_value_[i] -= (greedy_sum - scaled_total_sum);

219 scaled_vars_assigned_value_[i] =

220 std::max(scaled_vars_assigned_value_[i], scaled_vars_min_[i]);

221 }

222

223 greedy_sum += scaled_vars_assigned_value_[i] - old_scaled_vars_i;

224 }

225 }

226

227

228

229 void ComputeMaxWhenNoRepair() {

230 int num_overlap_sum_rounded_up = 0;

231 if (active_sum_nearest_ == active_sum_rounded_up_) {

232 num_overlap_sum_rounded_up = overlaps_.size();

233 }

234 for (int i = 0; i < size_; ++i) {

235 maximum_[i] = scaled_vars_assigned_value_[i];

236 if (Overlap(i) && active_sum_nearest_ == active_sum_rounded_up_ &&

237 active_sum_ % size_ != 0) {

238 overlaps_sup_[i] = num_overlap_sum_rounded_up - 1;

239 } else {

240 overlaps_sup_[i] = num_overlap_sum_rounded_up;

241 }

242 }

243 }

244

245

246

247

248 int ComputeNumOverlapsVariableRoundedUp() {

249 if (active_sum_ % size_ == 0) {

250 return 0;

251 }

252 int num_overlap_sum_rounded_up = 0;

253 for (int i = 0; i < size_; ++i) {

254 if (scaled_vars_assigned_value_[i] > scaled_vars_min_[i] &&

255 scaled_vars_assigned_value_[i] == active_sum_rounded_up_) {

256 num_overlap_sum_rounded_up++;

257 }

258 }

259 return num_overlap_sum_rounded_up;

260 }

261

262

263

264

265 bool CanPushSumAcrossMean(int64_t greedy_sum,

266 int64_t scaled_total_sum) const {

267 return (greedy_sum > scaled_total_sum &&

268 active_sum_nearest_ == active_sum_rounded_up_) ||

269 (greedy_sum < scaled_total_sum &&

270 active_sum_nearest_ == active_sum_rounded_down_);

271 }

272

273

274

275 void RepairSumAndComputeInfo(int64_t greedy_sum) {

276 const int64_t scaled_total_sum = size_ * active_sum_;

277

278

279

280 if (greedy_sum == scaled_total_sum) {

281 ComputeMaxWhenNoRepair();

282 } else {

283

284 if (CanPushSumAcrossMean(greedy_sum, scaled_total_sum)) {

285 const int64_t delta = greedy_sum > scaled_total_sum ? -size_ : size_;

286 for (int j = 0; j < overlaps_.size() && greedy_sum != scaled_total_sum;

287 ++j) {

288 scaled_vars_assigned_value_[overlaps_[j]] += delta;

289 greedy_sum += delta;

290 }

291 }

292

293 const int num_overlap_sum_rounded_up =

294 ComputeNumOverlapsVariableRoundedUp();

295

296 if (greedy_sum == scaled_total_sum) {

297

298 for (int i = 0; i < size_; ++i) {

299 if (Overlap(i) && num_overlap_sum_rounded_up > 0) {

300 maximum_[i] = active_sum_rounded_up_;

301 overlaps_sup_[i] = num_overlap_sum_rounded_up - 1;

302 } else {

303 maximum_[i] = scaled_vars_assigned_value_[i];

304 overlaps_sup_[i] = num_overlap_sum_rounded_up;

305 }

306 }

307 } else if (greedy_sum > scaled_total_sum) {

308

309

310

311 for (int i = 0; i < size_; ++i) {

312 maximum_[i] = scaled_vars_assigned_value_[i];

313 overlaps_sup_[i] = 0;

314 }

315 } else {

316 for (int i = 0; i < size_; ++i) {

317 if (Overlap(i) && num_overlap_sum_rounded_up > 0) {

318 overlaps_sup_[i] = num_overlap_sum_rounded_up - 1;

319 } else {

320 overlaps_sup_[i] = num_overlap_sum_rounded_up;

321 }

322

323 if (scaled_vars_assigned_value_[i] < scaled_vars_max_[i]) {

324 maximum_[i] =

325 scaled_vars_assigned_value_[i] + scaled_total_sum - greedy_sum;

326 } else {

327 maximum_[i] = scaled_vars_assigned_value_[i];

328 }

329 }

330 }

331 }

332 }

333

334

335 void PruneVars(int64_t min_delta, bool upper_bound) {

336

337 const int64_t increase_down_up = (active_sum_rounded_up_ - active_sum_) -

338 (active_sum_ - active_sum_rounded_down_);

339 for (int var_index = 0; var_index < size_; ++var_index) {

340

341 if (scaled_vars_max_[var_index] != scaled_vars_min_[var_index] &&

342 maximum_[var_index] < scaled_vars_max_[var_index]) {

343 const int64_t new_max =

344 ComputeNewMax(var_index, min_delta, increase_down_up);

345 PruneBound(var_index, new_max, upper_bound);

346 }

347 }

348 }

349

350

351 int64_t ComputeNewMax(int var_index, int64_t min_delta,

352 int64_t increase_down_up) {

353 int64_t maximum_value = maximum_[var_index];

354 int64_t current_min_delta = min_delta;

355

356 if (overlaps_sup_[var_index] > 0 &&

357 (current_min_delta +

358 overlaps_sup_[var_index] * (size_ - increase_down_up) >=

359 deviation_var_->Max())) {

360 const int64_t delta = deviation_var_->Max() - current_min_delta;

361 maximum_value += (size_ * delta) / (size_ - increase_down_up);

363 } else {

364 if (maximum_value == active_sum_rounded_down_ &&

365 active_sum_rounded_down_ < active_sum_) {

366 DCHECK_EQ(0, overlaps_sup_[var_index]);

367 current_min_delta += size_ + increase_down_up;

368 if (current_min_delta > deviation_var_->Max()) {

369 DCHECK_EQ(0, maximum_value % size_);

370 return maximum_value / size_;

371 }

372 maximum_value += size_;

373 }

374 current_min_delta +=

375 overlaps_sup_[var_index] * (size_ - increase_down_up);

376 maximum_value += size_ * overlaps_sup_[var_index];

377

378 const int64_t delta = deviation_var_->Max() - current_min_delta;

379 maximum_value += delta / 2;

381 }

382 }

383

384

385 void PruneBound(int var_index, int64_t bound, bool upper_bound) {

386 if (upper_bound) {

387 vars_[var_index]->SetMax(bound);

388 } else {

389 vars_[var_index]->SetMin(-bound);

390 }

391 }

392

393 std::vector<IntVar*> vars_;

394 const int size_;

395 IntVar* const deviation_var_;

396 const int64_t total_sum_;

397 std::unique_ptr<int64_t[]> scaled_vars_assigned_value_;

398 std::unique_ptr<int64_t[]> scaled_vars_min_;

399 std::unique_ptr<int64_t[]> scaled_vars_max_;

400 int64_t scaled_sum_max_;

401 int64_t scaled_sum_min_;

402

403 std::vector<int> overlaps_;

404 std::unique_ptr<int64_t[]> maximum_;

405 std::unique_ptr<int64_t[]> overlaps_sup_;

406

407 int64_t active_sum_;

408 int64_t active_sum_rounded_down_;

409 int64_t active_sum_rounded_up_;

410 int64_t active_sum_nearest_;

411};

412}

413

419}