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