1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the HiGHS linear optimization suite */
4 /* */
5 /* Written and engineered 2008-2021 at the University of Edinburgh */
6 /* */
7 /* Available as open-source under the MIT License */
8 /* */
9 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
10 #include "mip/HighsDomain.h"
11
12 #include <algorithm>
13 #include <cassert>
14 #include <numeric>
15 #include <queue>
16
17 #include "mip/HighsCutPool.h"
18 #include "mip/HighsMipSolverData.h"
19
activityContributionMin(double coef,const double & lb,const double & ub)20 static double activityContributionMin(double coef, const double& lb,
21 const double& ub) {
22 if (coef < 0) {
23 if (ub == HIGHS_CONST_INF) return -HIGHS_CONST_INF;
24
25 return coef * ub;
26 } else {
27 if (lb == -HIGHS_CONST_INF) return -HIGHS_CONST_INF;
28
29 return coef * lb;
30 }
31 }
32
activityContributionMax(double coef,const double & lb,const double & ub)33 static double activityContributionMax(double coef, const double& lb,
34 const double& ub) {
35 if (coef < 0) {
36 if (lb == -HIGHS_CONST_INF) return HIGHS_CONST_INF;
37
38 return coef * lb;
39 } else {
40 if (ub == HIGHS_CONST_INF) return HIGHS_CONST_INF;
41
42 return coef * ub;
43 }
44 }
45
HighsDomain(HighsMipSolver & mipsolver,HighsCutPool & cutpool)46 HighsDomain::HighsDomain(HighsMipSolver& mipsolver, HighsCutPool& cutpool)
47 : mipsolver(&mipsolver), cutpool(&cutpool), parentdomain(nullptr) {
48 colLower_ = mipsolver.model_->colLower_;
49 colUpper_ = mipsolver.model_->colUpper_;
50 changedcolsflags_.resize(mipsolver.numCol());
51 changedcols_.reserve(mipsolver.numCol());
52 }
53
computeMinActivity(int start,int end,const int * ARindex,const double * ARvalue,int & ninfmin,HighsCDouble & activitymin)54 void HighsDomain::computeMinActivity(int start, int end, const int* ARindex,
55 const double* ARvalue, int& ninfmin,
56 HighsCDouble& activitymin) {
57 activitymin = 0.0;
58 ninfmin = 0;
59 for (int j = start; j != end; ++j) {
60 int col = ARindex[j];
61 double val = ARvalue[j];
62
63 assert(col < int(colLower_.size()));
64
65 double contributionmin =
66 activityContributionMin(val, colLower_[col], colUpper_[col]);
67
68 if (contributionmin == -HIGHS_CONST_INF)
69 ++ninfmin;
70 else
71 activitymin += contributionmin;
72 }
73
74 activitymin.renormalize();
75 }
76
computeMaxActivity(int start,int end,const int * ARindex,const double * ARvalue,int & ninfmax,HighsCDouble & activitymax)77 void HighsDomain::computeMaxActivity(int start, int end, const int* ARindex,
78 const double* ARvalue, int& ninfmax,
79 HighsCDouble& activitymax) {
80 activitymax = 0.0;
81 ninfmax = 0;
82 for (int j = start; j != end; ++j) {
83 int col = ARindex[j];
84 double val = ARvalue[j];
85
86 assert(col < int(colLower_.size()));
87
88 double contributionmin =
89 activityContributionMax(val, colLower_[col], colUpper_[col]);
90
91 if (contributionmin == HIGHS_CONST_INF)
92 ++ninfmax;
93 else
94 activitymax += contributionmin;
95 }
96
97 activitymax.renormalize();
98 }
99
propagateRowUpper(const int * Rindex,const double * Rvalue,int Rlen,double Rupper,const HighsCDouble & minactivity,int ninfmin,HighsDomainChange * boundchgs)100 int HighsDomain::propagateRowUpper(const int* Rindex, const double* Rvalue,
101 int Rlen, double Rupper,
102 const HighsCDouble& minactivity, int ninfmin,
103 HighsDomainChange* boundchgs) {
104 assert(std::isfinite(double(minactivity)));
105 if (ninfmin > 1) return 0;
106 int numchgs = 0;
107 for (int i = 0; i != Rlen; ++i) {
108 HighsCDouble minresact;
109 double actcontribution = activityContributionMin(
110 Rvalue[i], colLower_[Rindex[i]], colUpper_[Rindex[i]]);
111 if (ninfmin == 1) {
112 if (actcontribution != -HIGHS_CONST_INF) continue;
113
114 minresact = minactivity;
115 } else {
116 minresact = minactivity - actcontribution;
117 }
118
119 double bound = double((Rupper - minresact) / Rvalue[i]);
120
121 if (Rvalue[i] > 0) {
122 bool accept;
123
124 if (mipsolver->variableType(Rindex[i]) != HighsVarType::CONTINUOUS) {
125 bound = std::floor(bound + mipsolver->mipdata_->feastol);
126 if (bound < colUpper_[Rindex[i]] &&
127 colUpper_[Rindex[i]] - bound >
128 1000.0 * mipsolver->mipdata_->feastol * std::abs(bound))
129 accept = true;
130 else
131 accept = false;
132 } else {
133 if (colUpper_[Rindex[i]] == HIGHS_CONST_INF)
134 accept = true;
135 else if (bound + 1000.0 * mipsolver->mipdata_->feastol <
136 colUpper_[Rindex[i]] &&
137 (colUpper_[Rindex[i]] - bound) /
138 std::max(std::abs(colUpper_[Rindex[i]]),
139 std::abs(bound)) >
140 0.3)
141 accept = true;
142 else
143 accept = false;
144 }
145
146 if (accept)
147 boundchgs[numchgs++] = {HighsBoundType::Upper, Rindex[i], bound};
148
149 } else {
150 bool accept;
151
152 if (mipsolver->variableType(Rindex[i]) != HighsVarType::CONTINUOUS) {
153 bound = std::ceil(bound - mipsolver->mipdata_->feastol);
154 if (bound > colLower_[Rindex[i]] &&
155 bound - colLower_[Rindex[i]] >
156 1000.0 * mipsolver->mipdata_->feastol * std::abs(bound))
157 accept = true;
158 else
159 accept = false;
160 } else {
161 if (colLower_[Rindex[i]] == -HIGHS_CONST_INF)
162 accept = true;
163 else if (bound - 1000.0 * mipsolver->mipdata_->feastol >
164 colLower_[Rindex[i]] &&
165 (bound - colLower_[Rindex[i]]) /
166 std::max(std::abs(colUpper_[Rindex[i]]),
167 std::abs(bound)) >
168 0.3)
169 accept = true;
170 else
171 accept = false;
172 }
173
174 if (accept)
175 boundchgs[numchgs++] = {HighsBoundType::Lower, Rindex[i], bound};
176 }
177 }
178
179 return numchgs;
180 }
181
propagateRowLower(const int * Rindex,const double * Rvalue,int Rlen,double Rlower,const HighsCDouble & maxactivity,int ninfmax,HighsDomainChange * boundchgs)182 int HighsDomain::propagateRowLower(const int* Rindex, const double* Rvalue,
183 int Rlen, double Rlower,
184 const HighsCDouble& maxactivity, int ninfmax,
185 HighsDomainChange* boundchgs) {
186 assert(std::isfinite(double(maxactivity)));
187 if (ninfmax > 1) return 0;
188 int numchgs = 0;
189 for (int i = 0; i != Rlen; ++i) {
190 HighsCDouble maxresact;
191 double actcontribution = activityContributionMax(
192 Rvalue[i], colLower_[Rindex[i]], colUpper_[Rindex[i]]);
193 if (ninfmax == 1) {
194 if (actcontribution != HIGHS_CONST_INF) continue;
195
196 maxresact = maxactivity;
197 } else {
198 maxresact = maxactivity - actcontribution;
199 }
200
201 double bound = double((Rlower - maxresact) / Rvalue[i]);
202
203 if (Rvalue[i] < 0) {
204 bool accept;
205
206 if (mipsolver->variableType(Rindex[i]) != HighsVarType::CONTINUOUS) {
207 bound = std::floor(bound + mipsolver->mipdata_->feastol);
208 if (bound < colUpper_[Rindex[i]] &&
209 colUpper_[Rindex[i]] - bound >
210 1000.0 * mipsolver->mipdata_->feastol * std::abs(bound))
211 accept = true;
212 else
213 accept = false;
214 } else {
215 if (colUpper_[Rindex[i]] == HIGHS_CONST_INF)
216 accept = true;
217 else if (bound + 1000.0 * mipsolver->mipdata_->feastol <
218 colUpper_[Rindex[i]] &&
219 (colUpper_[Rindex[i]] - bound) /
220 std::max(std::abs(colUpper_[Rindex[i]]),
221 std::abs(bound)) >
222 0.3)
223 accept = true;
224 else
225 accept = false;
226 }
227
228 if (accept)
229 boundchgs[numchgs++] = {HighsBoundType::Upper, Rindex[i], bound};
230 } else {
231 bool accept;
232
233 if (mipsolver->variableType(Rindex[i]) != HighsVarType::CONTINUOUS) {
234 bound = std::ceil(bound - mipsolver->mipdata_->feastol);
235 if (bound > colLower_[Rindex[i]] &&
236 bound - colLower_[Rindex[i]] >
237 1000.0 * mipsolver->mipdata_->feastol * std::abs(bound))
238 accept = true;
239 else
240 accept = false;
241 } else {
242 if (colLower_[Rindex[i]] == -HIGHS_CONST_INF)
243 accept = true;
244 else if (bound - 1000.0 * mipsolver->mipdata_->feastol >
245 colLower_[Rindex[i]] &&
246 (bound - colLower_[Rindex[i]]) /
247 std::max(std::abs(colUpper_[Rindex[i]]),
248 std::abs(bound)) >
249 0.3)
250 accept = true;
251 else
252 accept = false;
253 }
254
255 if (accept)
256 boundchgs[numchgs++] = {HighsBoundType::Lower, Rindex[i], bound};
257 }
258 }
259
260 return numchgs;
261 }
262
updateActivityLbChange(int col,double oldbound,double newbound)263 void HighsDomain::updateActivityLbChange(int col, double oldbound,
264 double newbound) {
265 auto mip = mipsolver->model_;
266 int start = mip->Astart_[col];
267 int end = mip->Astart_[col + 1];
268
269 for (int i = start; i != end; ++i) {
270 if (mip->Avalue_[i] > 0) {
271 double deltamin;
272 if (oldbound == -HIGHS_CONST_INF) {
273 --activitymininf_[mip->Aindex_[i]];
274 deltamin = newbound * mip->Avalue_[i];
275 } else if (newbound == -HIGHS_CONST_INF) {
276 ++activitymininf_[mip->Aindex_[i]];
277 deltamin = -oldbound * mip->Avalue_[i];
278 } else {
279 deltamin = (newbound - oldbound) * mip->Avalue_[i];
280 }
281 activitymin_[mip->Aindex_[i]] += deltamin;
282
283 if (mip->rowUpper_[mip->Aindex_[i]] != HIGHS_CONST_INF &&
284 activitymininf_[mip->Aindex_[i]] == 0 &&
285 activitymin_[mip->Aindex_[i]] - mip->rowUpper_[mip->Aindex_[i]] >
286 mipsolver->mipdata_->feastol) {
287 infeasible_ = mip->Aindex_[i] + 1;
288 }
289
290 if (deltamin > 0 && activitymininf_[mip->Aindex_[i]] <= 1 &&
291 !propagateflags_[mip->Aindex_[i]] &&
292 mip->rowUpper_[mip->Aindex_[i]] != HIGHS_CONST_INF) {
293 markPropagate(mip->Aindex_[i]);
294 // propagateflags_[mip->Aindex_[i]] = 1;
295 // propagateinds_.push_back(mip->Aindex_[i]);
296 }
297 } else {
298 double deltamax;
299 if (oldbound == -HIGHS_CONST_INF) {
300 --activitymaxinf_[mip->Aindex_[i]];
301 deltamax = newbound * mip->Avalue_[i];
302 } else if (newbound == -HIGHS_CONST_INF) {
303 ++activitymaxinf_[mip->Aindex_[i]];
304 deltamax = -oldbound * mip->Avalue_[i];
305 } else {
306 deltamax = (newbound - oldbound) * mip->Avalue_[i];
307 }
308 activitymax_[mip->Aindex_[i]] += deltamax;
309
310 if (mip->rowLower_[mip->Aindex_[i]] != -HIGHS_CONST_INF &&
311 activitymaxinf_[mip->Aindex_[i]] == 0 &&
312 mip->rowLower_[mip->Aindex_[i]] - activitymax_[mip->Aindex_[i]] >
313 mipsolver->mipdata_->feastol) {
314 infeasible_ = mip->Aindex_[i] + 1;
315 }
316
317 if (deltamax < 0 && activitymaxinf_[mip->Aindex_[i]] <= 1 &&
318 !propagateflags_[mip->Aindex_[i]] &&
319 mip->rowLower_[mip->Aindex_[i]] != -HIGHS_CONST_INF) {
320 markPropagate(mip->Aindex_[i]);
321 // propagateflags_[mip->Aindex_[i]] = 1;
322 // propagateinds_.push_back(mip->Aindex_[i]);
323 }
324 }
325 }
326
327 cutpool->getMatrix().forEachColumnEntry(col, [&](int row, double val) {
328 if (val > 0) {
329 double deltamin;
330
331 assert(row < int(activitycutversion_.size()));
332
333 if (activitycutversion_[row] != cutpool->getModificationCount(row)) {
334 int start = cutpool->getMatrix().getRowStart(row);
335 int end = cutpool->getMatrix().getRowEnd(row);
336 const int* arindex = cutpool->getMatrix().getARindex();
337 const double* arvalue = cutpool->getMatrix().getARvalue();
338
339 computeMinActivity(start, end, arindex, arvalue, activitycutsinf_[row],
340 activitycuts_[row]);
341
342 deltamin = HIGHS_CONST_INF;
343 } else {
344 if (oldbound == -HIGHS_CONST_INF) {
345 --activitycutsinf_[row];
346 deltamin = newbound * val;
347 } else if (newbound == -HIGHS_CONST_INF) {
348 ++activitycutsinf_[row];
349 deltamin = -oldbound * val;
350 } else {
351 deltamin = (newbound - oldbound) * val;
352 }
353 activitycuts_[row] += deltamin;
354 }
355
356 if (activitycutsinf_[row] == 0 &&
357 activitycuts_[row] - cutpool->getRhs()[row] >
358 mipsolver->mipdata_->feastol) {
359 infeasible_ = mip->numRow_ + row + 1;
360 }
361
362 if (deltamin > 0 && activitycutsinf_[row] <= 1 &&
363 !propagatecutflags_[row]) {
364 markPropagateCut(row);
365 // propagatecutflags_[row] = 1;
366 // propagatecutinds_.push_back(row);
367 }
368 }
369
370 return true;
371 });
372 }
373
updateActivityUbChange(int col,double oldbound,double newbound)374 void HighsDomain::updateActivityUbChange(int col, double oldbound,
375 double newbound) {
376 auto mip = mipsolver->model_;
377 int start = mip->Astart_[col];
378 int end = mip->Astart_[col + 1];
379
380 for (int i = start; i != end; ++i) {
381 if (mip->Avalue_[i] > 0) {
382 double deltamax;
383 if (oldbound == HIGHS_CONST_INF) {
384 --activitymaxinf_[mip->Aindex_[i]];
385 deltamax = newbound * mip->Avalue_[i];
386 } else if (newbound == HIGHS_CONST_INF) {
387 ++activitymaxinf_[mip->Aindex_[i]];
388 deltamax = -oldbound * mip->Avalue_[i];
389 } else {
390 deltamax = (newbound - oldbound) * mip->Avalue_[i];
391 }
392 activitymax_[mip->Aindex_[i]] += deltamax;
393
394 if (mip->rowLower_[mip->Aindex_[i]] != -HIGHS_CONST_INF &&
395 activitymaxinf_[mip->Aindex_[i]] == 0 &&
396 mip->rowLower_[mip->Aindex_[i]] - activitymax_[mip->Aindex_[i]] >
397 mipsolver->mipdata_->feastol) {
398 infeasible_ = mip->Aindex_[i] + 1;
399 }
400
401 if (deltamax < 0 && activitymaxinf_[mip->Aindex_[i]] <= 1 &&
402 !propagateflags_[mip->Aindex_[i]] &&
403 mip->rowLower_[mip->Aindex_[i]] != -HIGHS_CONST_INF) {
404 markPropagate(mip->Aindex_[i]);
405 // propagateflags_[mip->Aindex_[i]] = 1;
406 // propagateinds_.push_back(mip->Aindex_[i]);
407 }
408 } else {
409 double deltamin;
410 if (oldbound == HIGHS_CONST_INF) {
411 --activitymininf_[mip->Aindex_[i]];
412 deltamin = newbound * mip->Avalue_[i];
413 } else if (newbound == HIGHS_CONST_INF) {
414 ++activitymininf_[mip->Aindex_[i]];
415 deltamin = -oldbound * mip->Avalue_[i];
416 } else {
417 deltamin = (newbound - oldbound) * mip->Avalue_[i];
418 }
419
420 activitymin_[mip->Aindex_[i]] += deltamin;
421
422 if (mip->rowUpper_[mip->Aindex_[i]] != HIGHS_CONST_INF &&
423 activitymininf_[mip->Aindex_[i]] == 0 &&
424 activitymin_[mip->Aindex_[i]] - mip->rowUpper_[mip->Aindex_[i]] >
425 mipsolver->mipdata_->feastol) {
426 infeasible_ = mip->Aindex_[i] + 1;
427 }
428
429 if (deltamin > 0 && activitymininf_[mip->Aindex_[i]] <= 1 &&
430 !propagateflags_[mip->Aindex_[i]] &&
431 mip->rowUpper_[mip->Aindex_[i]] != HIGHS_CONST_INF) {
432 markPropagate(mip->Aindex_[i]);
433 // propagateflags_[mip->Aindex_[i]] = 1;
434 // propagateinds_.push_back(mip->Aindex_[i]);
435 }
436 }
437 }
438
439 cutpool->getMatrix().forEachColumnEntry(col, [&](int row, double val) {
440 if (val < 0) {
441 double deltamin;
442
443 assert(row < int(activitycutversion_.size()));
444
445 if (activitycutversion_[row] != cutpool->getModificationCount(row)) {
446 int start = cutpool->getMatrix().getRowStart(row);
447 int end = cutpool->getMatrix().getRowEnd(row);
448 const int* arindex = cutpool->getMatrix().getARindex();
449 const double* arvalue = cutpool->getMatrix().getARvalue();
450
451 computeMinActivity(start, end, arindex, arvalue, activitycutsinf_[row],
452 activitycuts_[row]);
453
454 activitycutversion_[row] = cutpool->getModificationCount(row);
455
456 deltamin = HIGHS_CONST_INF;
457 } else {
458 if (oldbound == HIGHS_CONST_INF) {
459 --activitycutsinf_[row];
460 deltamin = newbound * val;
461 } else if (newbound == HIGHS_CONST_INF) {
462 ++activitycutsinf_[row];
463 deltamin = -oldbound * val;
464 } else {
465 deltamin = (newbound - oldbound) * val;
466 }
467 activitycuts_[row] += deltamin;
468 }
469
470 if (activitycutsinf_[row] == 0 &&
471 activitycuts_[row] - cutpool->getRhs()[row] >
472 mipsolver->mipdata_->feastol) {
473 infeasible_ = mip->numRow_ + row + 1;
474 }
475
476 if (deltamin > 0 && activitycutsinf_[row] <= 1 &&
477 !propagatecutflags_[row]) {
478 markPropagateCut(row);
479 // propagatecutflags_[row] = 1;
480 // propagatecutinds_.push_back(row);
481 }
482 }
483
484 return true;
485 });
486 }
487
markPropagateCut(int cut)488 void HighsDomain::markPropagateCut(int cut) {
489 // todo, check if (rhs - minactivity) < amax - feastol and only mark in that
490 // case
491 if (!propagatecutflags_[cut] &&
492 (activitycutsinf_[cut] == 1 ||
493 (cutpool->getRhs()[cut] - activitycuts_[cut]) /
494 cutpool->getMaxAbsCutCoef(cut) <
495 1.0 - mipsolver->mipdata_->feastol)) {
496 propagatecutinds_.push_back(cut);
497 propagatecutflags_[cut] = 1;
498 }
499 }
500
markPropagate(int row)501 void HighsDomain::markPropagate(int row) {
502 // todo, check if std::min(maxactivity - lhs, rhs - minactivity) < amax -
503 // feastol and only mark in that case
504
505 if (!propagateflags_[row]) {
506 bool proplower = mipsolver->rowLower(row) != -HIGHS_CONST_INF &&
507 (activitymaxinf_[row] == 1 ||
508 (activitymax_[row] - mipsolver->rowLower(row)) /
509 mipsolver->mipdata_->maxAbsRowCoef[row] <
510 1.0 - mipsolver->mipdata_->feastol);
511 bool propupper = mipsolver->rowUpper(row) != HIGHS_CONST_INF &&
512 (activitymininf_[row] == 1 ||
513 (mipsolver->rowUpper(row) - activitymin_[row]) /
514 mipsolver->mipdata_->maxAbsRowCoef[row] <
515 1.0 - mipsolver->mipdata_->feastol);
516
517 if (proplower || propupper) {
518 propagateinds_.push_back(row);
519 propagateflags_[row] = 1;
520 }
521 }
522 }
523
cutAdded(int cut)524 void HighsDomain::cutAdded(int cut) {
525 int start = cutpool->getMatrix().getRowStart(cut);
526 int end = cutpool->getMatrix().getRowEnd(cut);
527 const int* arindex = cutpool->getMatrix().getARindex();
528 const double* arvalue = cutpool->getMatrix().getARvalue();
529
530 if (int(activitycuts_.size()) <= cut) {
531 activitycuts_.resize(cut + 1);
532 activitycutsinf_.resize(cut + 1);
533 propagatecutflags_.resize(cut + 1);
534 activitycutversion_.resize(cut + 1);
535 }
536
537 activitycutversion_[cut] = cutpool->getModificationCount(cut);
538 computeMinActivity(start, end, arindex, arvalue, activitycutsinf_[cut],
539 activitycuts_[cut]);
540
541 if (activitycutsinf_[cut] <= 1 && !propagatecutflags_[cut]) {
542 markPropagateCut(cut);
543 // propagatecutflags_[cut] = 1;
544 // propagatecutinds_.push_back(cut);
545 }
546
547 if (parentdomain != nullptr) parentdomain->cutAdded(cut);
548 }
549
computeRowActivities()550 void HighsDomain::computeRowActivities() {
551 activitymin_.resize(mipsolver->numRow());
552 activitymininf_.resize(mipsolver->numRow());
553 activitymax_.resize(mipsolver->numRow());
554 activitymaxinf_.resize(mipsolver->numRow());
555 propagateflags_.resize(mipsolver->numRow());
556 propagateinds_.reserve(mipsolver->numRow());
557
558 for (int i = 0; i != mipsolver->numRow(); ++i) {
559 int start = mipsolver->mipdata_->ARstart_[i];
560 int end = mipsolver->mipdata_->ARstart_[i + 1];
561
562 computeMinActivity(start, end, mipsolver->mipdata_->ARindex_.data(),
563 mipsolver->mipdata_->ARvalue_.data(), activitymininf_[i],
564 activitymin_[i]);
565 computeMaxActivity(start, end, mipsolver->mipdata_->ARindex_.data(),
566 mipsolver->mipdata_->ARvalue_.data(), activitymaxinf_[i],
567 activitymax_[i]);
568
569 if ((activitymininf_[i] <= 1 &&
570 mipsolver->rowUpper(i) != HIGHS_CONST_INF) ||
571 (activitymaxinf_[i] <= 1 &&
572 mipsolver->rowLower(i) != -HIGHS_CONST_INF)) {
573 markPropagate(i);
574 // propagateflags_[i] = 1;
575 // propagateinds_.push_back(i);
576 }
577 }
578 }
579
doChangeBound(const HighsDomainChange & boundchg)580 double HighsDomain::doChangeBound(const HighsDomainChange& boundchg) {
581 double oldbound;
582
583 if (boundchg.boundtype == HighsBoundType::Lower) {
584 oldbound = colLower_[boundchg.column];
585 colLower_[boundchg.column] = boundchg.boundval;
586 updateActivityLbChange(boundchg.column, oldbound, boundchg.boundval);
587
588 if (!changedcolsflags_[boundchg.column]) {
589 changedcolsflags_[boundchg.column] = 1;
590 changedcols_.push_back(boundchg.column);
591 }
592 } else {
593 oldbound = colUpper_[boundchg.column];
594 colUpper_[boundchg.column] = boundchg.boundval;
595 updateActivityUbChange(boundchg.column, oldbound, boundchg.boundval);
596
597 if (!changedcolsflags_[boundchg.column]) {
598 changedcolsflags_[boundchg.column] = 1;
599 changedcols_.push_back(boundchg.column);
600 }
601 }
602 assert(oldbound != boundchg.boundval);
603
604 return oldbound;
605 }
606
changeBound(HighsDomainChange boundchg,int reason)607 void HighsDomain::changeBound(HighsDomainChange boundchg, int reason) {
608 assert(boundchg.column >= 0);
609 assert(infeasible_ == 0);
610 if (boundchg.boundtype == HighsBoundType::Lower) {
611 if (boundchg.boundval <= colLower_[boundchg.column]) return;
612 if (boundchg.boundval > colUpper_[boundchg.column]) {
613 if (boundchg.boundval - colUpper_[boundchg.column] >
614 mipsolver->mipdata_->feastol) {
615 infeasible_ = reason >= 0 ? reason + 1 : reason;
616 return;
617 }
618
619 boundchg.boundval = colUpper_[boundchg.column];
620 if (boundchg.boundval == colLower_[boundchg.column]) return;
621 }
622 } else {
623 if (boundchg.boundval >= colUpper_[boundchg.column]) return;
624 if (boundchg.boundval < colLower_[boundchg.column]) {
625 if (colLower_[boundchg.column] - boundchg.boundval >
626 mipsolver->mipdata_->feastol) {
627 infeasible_ = reason >= 0 ? reason + 1 : reason;
628 return;
629 }
630
631 boundchg.boundval = colLower_[boundchg.column];
632 if (boundchg.boundval == colUpper_[boundchg.column]) return;
633 }
634 }
635
636 double oldbound = doChangeBound(boundchg);
637
638 prevboundval_.push_back(oldbound);
639 domchgstack_.push_back(boundchg);
640 domchgreason_.push_back(reason);
641 }
642
setDomainChangeStack(const std::vector<HighsDomainChange> & domchgstack)643 void HighsDomain::setDomainChangeStack(
644 const std::vector<HighsDomainChange>& domchgstack) {
645 infeasible_ = 0;
646
647 prevboundval_.clear();
648 domchgstack_.clear();
649 domchgreason_.clear();
650 int stacksize = domchgstack.size();
651 for (int k = 0; k != stacksize; ++k) {
652 if (domchgstack[k].boundtype == HighsBoundType::Lower &&
653 domchgstack[k].boundval <= colLower_[domchgstack[k].column])
654 continue;
655 if (domchgstack[k].boundtype == HighsBoundType::Upper &&
656 domchgstack[k].boundval >= colUpper_[domchgstack[k].column])
657 continue;
658
659 changeBound(domchgstack[k], -2);
660
661 if (infeasible_) break;
662 }
663 }
664
backtrack()665 HighsDomainChange HighsDomain::backtrack() {
666 int k = domchgstack_.size() - 1;
667
668 while (k >= 0) {
669 double prevbound = prevboundval_[k];
670
671 // change back to global bound
672 doChangeBound(
673 {domchgstack_[k].boundtype, domchgstack_[k].column, prevbound});
674
675 if (domchgreason_[k] == -1) break;
676
677 --k;
678 }
679
680 infeasible_ = 0;
681
682 if (k < 0) {
683 domchgstack_.clear();
684 prevboundval_.clear();
685 domchgreason_.clear();
686 return HighsDomainChange({HighsBoundType::Lower, -1, 0.0});
687 }
688
689 HighsDomainChange backtrackboundchg = domchgstack_[k];
690 domchgstack_.erase(domchgstack_.begin() + k, domchgstack_.end());
691 domchgreason_.resize(k);
692 prevboundval_.resize(k);
693
694 return backtrackboundchg;
695 }
696
propagate()697 void HighsDomain::propagate() {
698 std::vector<int> propagateinds;
699
700 #ifdef HIGHS_DEBUGSOL
701 bool debugsolactive = true;
702 HighsCDouble debugsolobj = 0;
703 for (int i = 0; i != mipsolver->numCol(); ++i) {
704 if (highsDebugSolution[i] + mipsolver->mipdata_->epsilon < colLower_[i] ||
705 highsDebugSolution[i] - mipsolver->mipdata_->epsilon > colUpper_[i]) {
706 debugsolactive = false;
707 }
708
709 debugsolobj += highsDebugSolution[i] * mipsolver->colCost(i);
710 }
711 #endif
712
713 if (propagateinds_.empty() && propagatecutinds_.empty()) return;
714
715 size_t changedboundsize = std::max(2 * mipsolver->mipdata_->ARvalue_.size(),
716 cutpool->getMatrix().nonzeroCapacity());
717 std::unique_ptr<HighsDomainChange[]> changedbounds(
718 new HighsDomainChange[changedboundsize]);
719
720 while (!propagateinds_.empty() || !propagatecutinds_.empty()) {
721 if (!propagateinds_.empty()) {
722 propagateinds.swap(propagateinds_);
723
724 int propnnz = 0;
725 int numproprows = propagateinds.size();
726 for (int i = 0; i != numproprows; ++i) {
727 int row = propagateinds[i];
728 propagateflags_[row] = 0;
729 propnnz += mipsolver->mipdata_->ARstart_[i + 1] -
730 mipsolver->mipdata_->ARstart_[i];
731 }
732
733 if (!infeasible_) {
734 propRowNumChangedBounds_.assign(numproprows, 0);
735
736 auto propagateIndex = [&](int k) {
737 // for (int k = 0; k != numproprows; ++k) {
738 int i = propagateinds[k];
739 int start = mipsolver->mipdata_->ARstart_[i];
740 int end = mipsolver->mipdata_->ARstart_[i + 1];
741 int Rlen = end - start;
742 const int* Rindex = &mipsolver->mipdata_->ARindex_[start];
743 const double* Rvalue = &mipsolver->mipdata_->ARvalue_[start];
744 int numchgs = 0;
745
746 if (mipsolver->rowUpper(i) != HIGHS_CONST_INF) {
747 // computeMinActivity(start, end, mipsolver->ARstart_.data(),
748 // mipsolver->ARvalue_.data(), activitymininf_[i],
749 // activitymin_[i]);
750 activitymin_[i].renormalize();
751 numchgs = propagateRowUpper(
752 Rindex, Rvalue, Rlen, mipsolver->rowUpper(i), activitymin_[i],
753 activitymininf_[i], &changedbounds[2 * start]);
754 }
755
756 if (mipsolver->rowLower(i) != -HIGHS_CONST_INF) {
757 // computeMaxActivity(start, end, mipsolver->ARstart_.data(),
758 // mipsolver->ARvalue_.data(), activitymaxinf_[i],
759 // activitymax_[i]);
760 activitymax_[i].renormalize();
761 numchgs += propagateRowLower(
762 Rindex, Rvalue, Rlen, mipsolver->rowLower(i), activitymax_[i],
763 activitymaxinf_[i], &changedbounds[2 * start + numchgs]);
764 }
765
766 propRowNumChangedBounds_[k] = numchgs;
767 };
768
769 // printf("numproprows (model): %d\n", numproprows);
770
771 for (int k = 0; k != numproprows; ++k) propagateIndex(k);
772
773 for (int k = 0; k != numproprows; ++k) {
774 if (propRowNumChangedBounds_[k] == 0) continue;
775 int i = propagateinds[k];
776 int start = 2 * mipsolver->mipdata_->ARstart_[i];
777 int end = start + propRowNumChangedBounds_[k];
778 for (int j = start; j != end && infeasible_ == 0; ++j)
779 changeBound(changedbounds[j], i);
780
781 if (infeasible_) break;
782 }
783 }
784
785 propagateinds.clear();
786 }
787
788 if (!propagatecutinds_.empty()) {
789 propagateinds.swap(propagatecutinds_);
790
791 int propnnz = 0;
792 int numproprows = propagateinds.size();
793
794 for (int i = 0; i != numproprows; ++i) {
795 int cut = propagateinds[i];
796 propagatecutflags_[cut] = 0;
797 propnnz += cutpool->getMatrix().getRowEnd(cut) -
798 cutpool->getMatrix().getRowStart(cut);
799 }
800
801 if (!infeasible_) {
802 propRowNumChangedBounds_.assign(numproprows, 0);
803
804 auto propagateIndex = [&](int k) {
805 int i = propagateinds[k];
806
807 int Rlen;
808 const int* Rindex;
809 const double* Rvalue;
810 cutpool->getCut(i, Rlen, Rindex, Rvalue);
811
812 if (activitycutversion_[i] != cutpool->getModificationCount(i)) {
813 activitycutversion_[i] = cutpool->getModificationCount(i);
814 int start = cutpool->getMatrix().getRowStart(i);
815 if (start == -1) {
816 activitycuts_[i] = 0;
817 return;
818 }
819 int end = cutpool->getMatrix().getRowEnd(i);
820 const int* arindex = cutpool->getMatrix().getARindex();
821 const double* arvalue = cutpool->getMatrix().getARvalue();
822
823 computeMinActivity(start, end, arindex, arvalue,
824 activitycutsinf_[i], activitycuts_[i]);
825 } else
826 activitycuts_[i].renormalize();
827
828 propRowNumChangedBounds_[k] = propagateRowUpper(
829 Rindex, Rvalue, Rlen, cutpool->getRhs()[i], activitycuts_[i],
830 activitycutsinf_[i],
831 &changedbounds[cutpool->getMatrix().getRowStart(i)]);
832 };
833
834 // printf("numproprows (cuts): %d\n", numproprows);
835
836 for (int k = 0; k != numproprows; ++k) propagateIndex(k);
837
838 for (int k = 0; k != numproprows; ++k) {
839 if (propRowNumChangedBounds_[k] == 0) continue;
840 int i = propagateinds[k];
841 cutpool->resetAge(i);
842 int start = cutpool->getMatrix().getRowStart(i);
843 int end = start + propRowNumChangedBounds_[k];
844 for (int j = start; j != end && infeasible_ == 0; ++j)
845 changeBound(changedbounds[j], i);
846
847 if (infeasible_) break;
848 }
849 }
850
851 propagateinds.clear();
852 }
853 }
854
855 #ifdef HIGHS_DEBUGSOL
856 if (debugsolactive && mipsolver->mipdata_->upper_bound >
857 debugsolobj + mipsolver->mipdata_->epsilon) {
858 assert(!infeasible_);
859 for (int i = 0; i != mipsolver->numCol(); ++i) {
860 if (highsDebugSolution[i] + mipsolver->mipdata_->epsilon < colLower_[i] ||
861 highsDebugSolution[i] - mipsolver->mipdata_->epsilon > colUpper_[i]) {
862 assert(false);
863 }
864 }
865 }
866 #endif
867 }
868
tightenCoefficients(int * inds,double * vals,int len,double & rhs) const869 void HighsDomain::tightenCoefficients(int* inds, double* vals, int len,
870 double& rhs) const {
871 HighsCDouble maxactivity = 0;
872
873 for (int i = 0; i != len; ++i) {
874 if (vals[i] > 0) {
875 if (colUpper_[inds[i]] == HIGHS_CONST_INF) return;
876
877 maxactivity += colUpper_[inds[i]] * vals[i];
878 } else {
879 if (colLower_[inds[i]] == -HIGHS_CONST_INF) return;
880
881 maxactivity += colLower_[inds[i]] * vals[i];
882 }
883 }
884
885 if (maxactivity - rhs > mipsolver->mipdata_->feastol) {
886 HighsCDouble upper = rhs;
887 HighsCDouble maxabscoef = double(maxactivity - rhs);
888 int tightened = 0;
889 for (int i = 0; i != len; ++i) {
890 if (vals[i] > maxabscoef) {
891 HighsCDouble delta = vals[i] - maxabscoef;
892 upper -= delta * colUpper_[inds[i]];
893 vals[i] = double(maxabscoef);
894 ++tightened;
895 } else if (vals[i] < -maxabscoef) {
896 HighsCDouble delta = -vals[i] - maxabscoef;
897 upper += delta * colLower_[inds[i]];
898 vals[i] = -double(maxabscoef);
899 ++tightened;
900 }
901 }
902
903 if (tightened != 0) {
904 // printf("tightened %d coefficients, rhs changed from %g to %g\n",
905 // tightened, rhs, double(upper));
906 rhs = double(upper);
907 }
908 }
909 }
910
isFixing(const HighsDomainChange & domchg) const911 bool HighsDomain::isFixing(const HighsDomainChange& domchg) const {
912 double otherbound = domchg.boundtype == HighsBoundType::Upper
913 ? colLower_[domchg.column]
914 : colUpper_[domchg.column];
915 return std::abs(domchg.boundval - otherbound) <= mipsolver->mipdata_->epsilon;
916 }