1 /* $Id: CoinPresolveTighten.cpp 2083 2019-01-06 19:38:09Z unxusr $ */
2 // Copyright (C) 2002, International Business Machines
3 // Corporation and others. All Rights Reserved.
4 // This code is licensed under the terms of the Eclipse Public License (EPL).
5
6 #include <stdio.h>
7 #include <math.h>
8
9 #include "CoinPresolveMatrix.hpp"
10 #include "CoinPresolveFixed.hpp"
11 #include "CoinPresolveTighten.hpp"
12 #include "CoinPresolveUseless.hpp"
13 #include "CoinHelperFunctions.hpp"
14 #include "CoinFinite.hpp"
15
16 #if PRESOLVE_DEBUG > 0 || PRESOLVE_CONSISTENCY > 0
17 #include "CoinPresolvePsdebug.hpp"
18 #endif
19
name() const20 const char *do_tighten_action::name() const
21 {
22 return ("do_tighten_action");
23 }
24
25 /*
26 This is ekkredc2. This fairly simple transformation is not mentioned
27 in the paper. Say there is a costless variable x<t> such that all the
28 constraints it's entangled with (i.e., a<it> != 0) would be satisfied
29 as it approaches plus or minus infinity, because all its constraints
30 have only one bound, and increasing/decreasing the variable makes the
31 row activity grow away from the bound (in the right direction).
32
33 If x<j> is unbounded in that direction, it can always be made large
34 enough to satisfy the constraints, so we can just drop the variable and
35 the entangled constraints from the problem.
36
37 If x<j> *is* bounded in that direction, there is no reason not to set
38 it to that bound. This effectively weakens the constraints, which in
39 fact may be subsequently presolved away.
40
41 Note that none of the constraints may be bounded both above and below,
42 since then we don't know which way to move the variable in order to
43 satisfy the constraint.
44
45 To drop constraints, we just make them useless and let other transformations
46 take care of the rest.
47
48 Note that more than one such costless unbounded variable may be part of
49 a given constraint. In that case, the first one processed will make
50 the constraint useless, and the second will ignore it. In postsolve,
51 the first will be responsible for satisfying the constraint.
52
53 Note that if the constraints are dropped (as in the first case), then we
54 just make them useless. It will subsequently be discovered the the variable
55 does not appear in any constraints, and since it has no cost it is just set
56 to some value (either zero or a bound) and removed (by remove_empty_cols).
57
58 Oddly, pilots and baxter do *worse* when this transform is applied.
59
60 It's informative to compare this transform to the very similar transform
61 implemented in remove_dual_action. Surely they could be merged.
62 */
63
presolve(CoinPresolveMatrix * prob,const CoinPresolveAction * next)64 const CoinPresolveAction *do_tighten_action::presolve(CoinPresolveMatrix *prob,
65 const CoinPresolveAction *next)
66 {
67 double *colels = prob->colels_;
68 int *hrow = prob->hrow_;
69 CoinBigIndex *mcstrt = prob->mcstrt_;
70 int *hincol = prob->hincol_;
71 int ncols = prob->ncols_;
72
73 double *clo = prob->clo_;
74 double *cup = prob->cup_;
75
76 double *rlo = prob->rlo_;
77 double *rup = prob->rup_;
78
79 double *dcost = prob->cost_;
80
81 const unsigned char *integerType = prob->integerType_;
82
83 int *fix_cols = prob->usefulColumnInt_;
84 int nfixup_cols = 0;
85
86 int nfixdown_cols = ncols;
87
88 int *useless_rows = prob->usefulRowInt_;
89 int nuseless_rows = 0;
90
91 action *actions = new action[ncols];
92 int nactions = 0;
93
94 int numberLook = prob->numberColsToDo_;
95 int iLook;
96 int *look = prob->colsToDo_;
97 bool fixInfeasibility = ((prob->presolveOptions_ & 0x4000) != 0);
98
99 #if PRESOLVE_DEBUG > 0 || PRESOLVE_CONSISTENCY > 0
100 #if PRESOLVE_DEBUG > 0
101 std::cout
102 << "Entering do_tighten_action::presolve; considering " << numberLook
103 << " rows." << std::endl;
104 #endif
105 presolve_consistent(prob);
106 presolve_links_ok(prob);
107 presolve_check_sol(prob);
108 presolve_check_nbasic(prob);
109 #endif
110
111 #if PRESOLVE_DEBUG > 0 || COIN_PRESOLVE_TUNING > 0
112 int startEmptyRows = 0;
113 int startEmptyColumns = 0;
114 startEmptyRows = prob->countEmptyRows();
115 startEmptyColumns = prob->countEmptyCols();
116 #if COIN_PRESOLVE_TUNING > 0
117 double startTime = 0.0;
118 if (prob->tuning_)
119 startTime = CoinCpuTime();
120 #endif
121 #endif
122
123 // singleton columns are especially likely to be caught here
124 for (iLook = 0; iLook < numberLook; iLook++) {
125 int j = look[iLook];
126 // modify bounds if integer
127 if (integerType[j]) {
128 clo[j] = ceil(clo[j] - 1.0e-12);
129 cup[j] = floor(cup[j] + 1.0e-12);
130 if (clo[j] > cup[j] && !fixInfeasibility) {
131 // infeasible
132 prob->status_ |= 1;
133 prob->messageHandler()->message(COIN_PRESOLVE_COLINFEAS,
134 prob->messages())
135 << j
136 << clo[j]
137 << cup[j]
138 << CoinMessageEol;
139 }
140 }
141 if (dcost[j] == 0.0 && !prob->colProhibited2(j)) {
142 int iflag = 0; /* 1 - up is towards feasibility, -1 down is towards */
143 int nonFree = 0; // Number of non-free rows
144
145 CoinBigIndex kcs = mcstrt[j];
146 CoinBigIndex kce = kcs + hincol[j];
147
148 // check constraints
149 for (CoinBigIndex k = kcs; k < kce; ++k) {
150 int i = hrow[k];
151 double coeff = colels[k];
152 double rlb = rlo[i];
153 double rub = rup[i];
154
155 if (-1.0e28 < rlb && rub < 1.0e28) {
156 // bounded - we lose
157 iflag = 0;
158 break;
159 } else if (-1.0e28 < rlb || rub < 1.0e28) {
160 nonFree++;
161 }
162
163 PRESOLVEASSERT(fabs(coeff) > ZTOLDP);
164
165 // see what this particular row says
166 // jflag == 1 ==> up is towards feasibility
167 int jflag = (coeff > 0.0
168 ? (rub > 1.0e28 ? 1 : -1)
169 : (rlb < -1.0e28 ? 1 : -1));
170
171 if (iflag) {
172 // check that it agrees with iflag.
173 if (iflag != jflag) {
174 iflag = 0;
175 break;
176 }
177 } else {
178 // first row -- initialize iflag
179 iflag = jflag;
180 }
181 }
182 // done checking constraints
183 if (!nonFree)
184 iflag = 0; // all free anyway
185 if (iflag) {
186 if (iflag == 1 && cup[j] < 1.0e10) {
187 #if PRESOLVE_DEBUG > 1
188 printf("TIGHTEN UP: %d\n", j);
189 #endif
190 fix_cols[nfixup_cols++] = j;
191
192 } else if (iflag == -1 && clo[j] > -1.0e10) {
193 // symmetric case
194 //mpre[j] = PRESOLVE_XUP;
195
196 #if PRESOLVE_DEBUG > 1
197 printf("TIGHTEN DOWN: %d\n", j);
198 #endif
199
200 fix_cols[--nfixdown_cols] = j;
201
202 } else {
203 #if 0
204 static int limit;
205 static int which = atoi(getenv("WZ"));
206 if (which == -1)
207 ;
208 else if (limit != which) {
209 limit++;
210 continue;
211 } else
212 limit++;
213
214 printf("TIGHTEN STATS %d %g %g %d: \n", j, clo[j], cup[j], integerType[j]);
215 double *rowels = prob->rowels_;
216 int *hcol = prob->hcol_;
217 int *mrstrt = prob->mrstrt_;
218 int *hinrow = prob->hinrow_;
219 for (CoinBigIndex k=kcs; k<kce; ++k) {
220 int irow = hrow[k];
221 CoinBigIndex krs = mrstrt[irow];
222 CoinBigIndex kre = krs + hinrow[irow];
223 printf("%d %g %g %g: ",
224 irow, rlo[irow], rup[irow], colels[irow]);
225 for (CoinBigIndex kk=krs; kk<kre; ++kk)
226 printf("%d(%g) ", hcol[kk], rowels[kk]);
227 printf("\n");
228 }
229 #endif
230
231 {
232 action *s = &actions[nactions];
233 nactions++;
234 s->col = j;
235 PRESOLVE_DETAIL_PRINT(printf("pre_tighten %dC E\n", j));
236 if (integerType[j]) {
237 assert(iflag == -1 || iflag == 1);
238 iflag *= 2; // say integer
239 }
240 s->direction = iflag;
241
242 s->rows = new int[hincol[j]];
243 s->lbound = new double[hincol[j]];
244 s->ubound = new double[hincol[j]];
245 #if PRESOLVE_DEBUG > 1
246 printf("TIGHTEN FREE: %d ", j);
247 #endif
248 int nr = 0;
249 prob->addCol(j);
250 for (CoinBigIndex k = kcs; k < kce; ++k) {
251 int irow = hrow[k];
252 // ignore this if we've already made it useless
253 if (!(rlo[irow] == -PRESOLVE_INF && rup[irow] == PRESOLVE_INF)) {
254 prob->addRow(irow);
255 s->rows[nr] = irow;
256 s->lbound[nr] = rlo[irow];
257 s->ubound[nr] = rup[irow];
258 nr++;
259
260 useless_rows[nuseless_rows++] = irow;
261
262 rlo[irow] = -PRESOLVE_INF;
263 rup[irow] = PRESOLVE_INF;
264
265 #if PRESOLVE_DEBUG > 1
266 printf("%d ", irow);
267 #endif
268 }
269 }
270 s->nrows = nr;
271
272 #if PRESOLVE_DEBUG > 1
273 printf("\n");
274 #endif
275 }
276 }
277 }
278 }
279 }
280
281 #if PRESOLVE_SUMMARY > 0
282 if (nfixdown_cols < ncols || nfixup_cols || nuseless_rows) {
283 printf("NTIGHTENED: %d %d %d\n", ncols - nfixdown_cols, nfixup_cols, nuseless_rows);
284 }
285 #endif
286
287 if (nuseless_rows) {
288 next = new do_tighten_action(nactions, CoinCopyOfArray(actions, nactions), next);
289
290 next = useless_constraint_action::presolve(prob,
291 useless_rows, nuseless_rows,
292 next);
293 }
294 deleteAction(actions, action *);
295 //delete[]useless_rows;
296
297 if (nfixdown_cols < ncols) {
298 int *fixdown_cols = fix_cols + nfixdown_cols;
299 nfixdown_cols = ncols - nfixdown_cols;
300 next = make_fixed_action::presolve(prob, fixdown_cols, nfixdown_cols,
301 true,
302 next);
303 }
304 //delete[]fixdown_cols;
305
306 if (nfixup_cols) {
307 next = make_fixed_action::presolve(prob, fix_cols, nfixup_cols,
308 false,
309 next);
310 }
311 //delete[]fixup_cols;
312
313 #if COIN_PRESOLVE_TUNING > 0
314 double thisTime = 0.0;
315 if (prob->tuning_)
316 thisTime = CoinCpuTime();
317 #endif
318 #if PRESOLVE_CONSISTENCY > 0 || PRESOLVE_DEBUG > 0
319 presolve_check_sol(prob);
320 #endif
321 #if PRESOLVE_DEBUG > 0 || COIN_PRESOLVE_TUNING > 0
322 int droppedRows = prob->countEmptyRows() - startEmptyRows;
323 int droppedColumns = prob->countEmptyCols() - startEmptyColumns;
324 std::cout
325 << "Leaving do_tighten_action::presolve, " << droppedRows << " rows, "
326 << droppedColumns << " columns dropped";
327 #if COIN_PRESOLVE_TUNING > 0
328 std::cout << " in " << thisTime - startTime << "s";
329 #endif
330 std::cout << "." << std::endl;
331 #endif
332
333 return (next);
334 }
335
postsolve(CoinPostsolveMatrix * prob) const336 void do_tighten_action::postsolve(CoinPostsolveMatrix *prob) const
337 {
338 const action *const actions = actions_;
339 const int nactions = nactions_;
340
341 double *colels = prob->colels_;
342 int *hrow = prob->hrow_;
343 CoinBigIndex *mcstrt = prob->mcstrt_;
344 int *hincol = prob->hincol_;
345 CoinBigIndex *link = prob->link_;
346
347 double *clo = prob->clo_;
348 double *cup = prob->cup_;
349 double *rlo = prob->rlo_;
350 double *rup = prob->rup_;
351
352 double *sol = prob->sol_;
353 double *acts = prob->acts_;
354
355 #if PRESOLVE_DEBUG > 0 || PRESOLVE_CONSISTENCY > 0
356 char *cdone = prob->cdone_;
357 char *rdone = prob->rdone_;
358
359 presolve_check_threads(prob);
360 presolve_check_sol(prob, 2, 2, 2);
361 presolve_check_nbasic(prob);
362
363 #if PRESOLVE_DEBUG > 0
364 std::cout << "Entering do_tighten_action::postsolve." << std::endl;
365 #endif
366 #endif
367
368 for (const action *f = &actions[nactions - 1]; actions <= f; f--) {
369 int jcol = f->col;
370 int iflag = f->direction;
371 int nr = f->nrows;
372 const int *rows = f->rows;
373 const double *lbound = f->lbound;
374 const double *ubound = f->ubound;
375
376 PRESOLVEASSERT(prob->getColumnStatus(jcol) != CoinPrePostsolveMatrix::basic);
377 int i;
378 for (i = 0; i < nr; ++i) {
379 int irow = rows[i];
380
381 rlo[irow] = lbound[i];
382 rup[irow] = ubound[i];
383
384 PRESOLVEASSERT(prob->getRowStatus(irow) == CoinPrePostsolveMatrix::basic);
385 }
386
387 // We have just tightened the row bounds.
388 // That means we'll have to compute a new value
389 // for this variable that will satisfy everybody.
390 // We are supposed to be in a position where this
391 // is always possible.
392
393 // Each constraint has exactly one bound.
394 // The correction should only ever be forced to move in one direction.
395 // double orig_sol = sol[jcol];
396 double correction = 0.0;
397
398 int last_corrected = -1;
399 CoinBigIndex k = mcstrt[jcol];
400 int nk = hincol[jcol];
401 for (i = 0; i < nk; ++i) {
402 int irow = hrow[k];
403 double coeff = colels[k];
404 k = link[k];
405 double newrlo = rlo[irow];
406 double newrup = rup[irow];
407 double activity = acts[irow];
408
409 if (activity + correction * coeff < newrlo) {
410 // only one of these two should fire
411 PRESOLVEASSERT(!(activity + correction * coeff > newrup));
412
413 last_corrected = irow;
414
415 // adjust to just meet newrlo (solve for correction)
416 double new_correction = (newrlo - activity) / coeff;
417 //adjust if integer
418 if (iflag == -2 || iflag == 2) {
419 new_correction += sol[jcol];
420 if (fabs(floor(new_correction + 0.5) - new_correction) > 1.0e-4) {
421 new_correction = ceil(new_correction) - sol[jcol];
422 #ifdef COIN_DEVELOP
423 printf("integer postsolve changing correction from %g to %g - flag %d\n",
424 (newrlo - activity) / coeff, new_correction, iflag);
425 #endif
426 }
427 }
428 correction = new_correction;
429 } else if (activity + correction * coeff > newrup) {
430 last_corrected = irow;
431
432 double new_correction = (newrup - activity) / coeff;
433 //adjust if integer
434 if (iflag == -2 || iflag == 2) {
435 new_correction += sol[jcol];
436 if (fabs(floor(new_correction + 0.5) - new_correction) > 1.0e-4) {
437 new_correction = ceil(new_correction) - sol[jcol];
438 #ifdef COIN_DEVELOP
439 printf("integer postsolve changing correction from %g to %g - flag %d\n",
440 (newrup - activity) / coeff, new_correction, iflag);
441 #endif
442 }
443 }
444 correction = new_correction;
445 }
446 }
447
448 if (last_corrected >= 0) {
449 sol[jcol] += correction;
450
451 // by construction, the last row corrected (if there was one)
452 // must be at its bound, so it can be non-basic.
453 // All other rows may not be at a bound (but may if the difference
454 // is very small, causing a new correction by a tiny amount).
455
456 // now adjust the activities
457 k = mcstrt[jcol];
458 for (i = 0; i < nk; ++i) {
459 int irow = hrow[k];
460 double coeff = colels[k];
461 k = link[k];
462 // double activity = acts[irow];
463
464 acts[irow] += correction * coeff;
465 }
466 /*
467 If the col happens to get pushed to its bound, we may as well leave
468 it non-basic. Otherwise, set the status to basic.
469
470 Why do we correct the row status only when the column is made basic?
471 Need to look at preceding code. -- lh, 110528 --
472 */
473 if (fabs(sol[jcol] - clo[jcol]) > ZTOLDP && fabs(sol[jcol] - cup[jcol]) > ZTOLDP) {
474
475 prob->setColumnStatus(jcol, CoinPrePostsolveMatrix::basic);
476 if (acts[last_corrected] - rlo[last_corrected] < rup[last_corrected] - acts[last_corrected])
477 prob->setRowStatus(last_corrected,
478 CoinPrePostsolveMatrix::atUpperBound);
479 else
480 prob->setRowStatus(last_corrected,
481 CoinPrePostsolveMatrix::atLowerBound);
482 }
483 }
484 }
485
486 #if PRESOLVE_DEBUG > 0 || PRESOLVE_CONSISTENCY > 0
487 presolve_check_threads(prob);
488 presolve_check_sol(prob, 2, 2, 2);
489 presolve_check_nbasic(prob);
490 #if PRESOLVE_DEBUG > 0
491 std::cout << "Leaving do_tighten_action::postsolve." << std::endl;
492 #endif
493 #endif
494 }
495
~do_tighten_action()496 do_tighten_action::~do_tighten_action()
497 {
498 if (nactions_ > 0) {
499 for (int i = nactions_ - 1; i >= 0; --i) {
500 delete[] actions_[i].rows;
501 delete[] actions_[i].lbound;
502 delete[] actions_[i].ubound;
503 }
504 deleteAction(actions_, action *);
505 }
506 }
507
508 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
509 */
510