source: trunk/PresolveTighten.cpp @ 198

Last change on this file since 198 was 198, checked in by forrest, 17 years ago

Fix bug

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.2 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#include <stdio.h>
4#include <math.h>
5
6#include "PresolveMatrix.hpp"
7#include "PresolveFixed.hpp"
8#include "PresolveTighten.hpp"
9#include "PresolveUseless.hpp"
10
11const char *do_tighten_action::name() const
12{
13  return ("do_tighten_action");
14}
15
16// This is ekkredc2.
17// This fairly simple transformation is not mentioned in the paper.
18// Say there is a costless variable such all its constraints
19// would be satisfied as it approaches plus or minus infinity,
20// because all its constraints have only one bound, and increasing/decreasing
21// the variable makes the row activity grow away from the bound
22// (in the right direction).
23//
24// If the variable is unbounded in that direction,
25// that means we can determine right now how large it needs
26// to get in order to satisfy the constraints, so we can
27// just drop the variable and those constraints from the problem.
28//
29// If the variable *is* bounded in that direction,
30// there is no reason not to set it to that bound.
31// This effectively weakens the constraints, and in fact
32// may be subsequently presolved away.
33//
34// Note that none of the constraints may be bounded both above and below,
35// since then we don't know which way to move the variable in order
36// to satisfy the constraint.
37//
38// To drop constraints, we just make them useless and let other
39// transformations take care of the rest.
40//
41// Note that more than one such costless unbounded variable
42// may be part of a given constraint.
43// In that case, the first one processed will make the
44// constraint useless, and the second will ignore it.
45// In postsolve, the first will be responsible for satisfying
46// the constraint.
47//
48// Note that if the constraints are dropped (as in the first case),
49// then we just make them useless.  It is subsequently discovered
50// the the variable does not appear in any constraints, and since it
51// has no cost it is just set to some value (either zero or a bound)
52// and removed (by remove_empty_cols).
53//
54// oddly, pilots and baxter do *worse* when this transform is applied.
55const PresolveAction *do_tighten_action::presolve(PresolveMatrix *prob,
56                                                   const PresolveAction *next)
57{
58  double *colels        = prob->colels_;
59  int *hrow             = prob->hrow_;
60  CoinBigIndex *mcstrt          = prob->mcstrt_;
61  int *hincol           = prob->hincol_;
62  int ncols             = prob->ncols_;
63
64  int nrows             = prob->nrows_;
65
66  double *clo   = prob->clo_;
67  double *cup   = prob->cup_;
68
69  double *rlo   = prob->rlo_;
70  double *rup   = prob->rup_;
71
72  double *dcost = prob->cost_;
73
74  // NEED TO THINK MORE ABOUT INTEGER VARS
75  //  const char *integerType = prob->integerType_;
76
77  int *fixup_cols       = new int[ncols];
78  int nfixup_cols       = 0;
79
80  int *fixdown_cols     = new int[ncols];
81  int nfixdown_cols     = 0;
82
83  int *useless_rows     = new int[nrows];
84  int nuseless_rows     = 0;
85 
86  action *actions       = new action [ncols];
87  int nactions          = 0;
88
89  int numberLook = prob->numberColsToDo_;
90  int iLook;
91  int * look = prob->colsToDo_;
92
93  // singleton columns are especially likely to be caught here
94  for (iLook=0;iLook<numberLook;iLook++) {
95    int j = look[iLook];
96    if (dcost[j]==0.0) {
97      int iflag=0; /* 1 - up is towards feasibility, -1 down is towards */
98      int nonFree=0; // Number of non-free rows
99
100      CoinBigIndex kcs = mcstrt[j];
101      CoinBigIndex kce = kcs + hincol[j];
102
103      // check constraints
104      for (CoinBigIndex k=kcs; k<kce; ++k) {
105        int i = hrow[k];
106        double coeff = colels[k];
107        double rlb = rlo[i];
108        double rub = rup[i];
109
110        if (-1.0e28 < rlb && rub < 1.0e28) {
111          // bounded - we lose
112          iflag=0;
113          break;
114        } else if (-1.0e28 < rlb || rub < 1.0e28) {
115          nonFree++;
116        }
117
118        PRESOLVEASSERT(fabs(coeff) > ZTOLDP);
119
120        // see what this particular row says
121        // jflag == 1 ==> up is towards feasibility
122        int jflag = (coeff > 0.0
123                     ? (rub >  1.0e28 ? 1 : -1)
124                     : (rlb < -1.0e28 ? 1 : -1));
125
126        if (iflag) {
127          // check that it agrees with iflag.
128          if (iflag!=jflag) {
129            iflag=0;
130            break;
131          }
132        } else {
133          // first row -- initialize iflag
134          iflag=jflag;
135        }
136      }
137      // done checking constraints
138
139      if (!nonFree)
140        iflag=0; // all free anyway
141      if (iflag) {
142        if (iflag==1 && cup[j]<1.0e10) {
143#if     DEBUG_PRESOLVE
144          printf("TIGHTEN UP:  %d\n", j);
145#endif
146          fixup_cols[nfixup_cols] = j;
147          ++nfixup_cols;
148
149        } else if (iflag==-1&&clo[j]>-1.0e10) {
150          // symmetric case
151          //mpre[j] = PRESOLVE_XUP;
152
153#if     DEBUG_PRESOLVE
154          printf("TIGHTEN DOWN:  %d\n", j);
155#endif
156
157          fixdown_cols[nfixdown_cols] = j;
158          ++nfixdown_cols;
159
160        } else {
161#if 0
162          static int limit;
163          static int which = atoi(getenv("WZ"));
164          if (which == -1)
165            ;
166          else if (limit != which) {
167            limit++;
168            continue;
169          } else
170            limit++;
171
172          printf("TIGHTEN STATS %d %g %g %d:  \n", j, clo[j], cup[j], integerType[j]);
173  double *rowels        = prob->rowels_;
174  int *hcol             = prob->hcol_;
175  int *mrstrt           = prob->mrstrt_;
176  int *hinrow           = prob->hinrow_;
177          for (CoinBigIndex k=kcs; k<kce; ++k) {
178            int irow = hrow[k];
179            CoinBigIndex krs = mrstrt[irow];
180            CoinBigIndex kre = krs + hinrow[irow];
181            printf("%d  %g %g %g:  ",
182                   irow, rlo[irow], rup[irow], colels[irow]);
183            for (CoinBigIndex kk=krs; kk<kre; ++kk)
184              printf("%d(%g) ", hcol[kk], rowels[kk]);
185            printf("\n");
186          }
187#endif
188
189          {
190            action *s = &actions[nactions];       
191            nactions++;
192            s->col = j;
193            s->direction = iflag;
194
195            s->rows =   new int[hincol[j]];
196            s->lbound = new double[hincol[j]];
197            s->ubound = new double[hincol[j]];
198#ifdef DEBUG_PRESOLVE
199            printf("TIGHTEN FREE:  %d   ", j);
200#endif
201            int nr = 0;
202            prob->addCol(j);
203            for (CoinBigIndex k=kcs; k<kce; ++k) {
204              int irow = hrow[k];
205              // ignore this if we've already made it useless
206              if (! (rlo[irow] == -PRESOLVE_INF && rup[irow] == PRESOLVE_INF)) {
207                prob->addRow(irow);
208                s->rows  [nr] = irow;
209                s->lbound[nr] = rlo[irow];
210                s->ubound[nr] = rup[irow];
211                nr++;
212
213                useless_rows[nuseless_rows++] = irow;
214
215                rlo[irow] = -PRESOLVE_INF;
216                rup[irow] = PRESOLVE_INF;
217
218#ifdef DEBUG_PRESOLVE
219                printf("%d ", irow);
220#endif
221              }
222            }
223            s->nrows = nr;
224
225#ifdef DEBUG_PRESOLVE
226            printf("\n");
227#endif
228          }
229        }
230      }
231    }
232  }
233
234
235#if     PRESOLVE_SUMMARY
236  if (nfixdown_cols || nfixup_cols || nuseless_rows) {
237    printf("NTIGHTENED:  %d %d %d\n", nfixdown_cols, nfixup_cols, nuseless_rows);
238  }
239#endif
240
241  if (nuseless_rows) {
242    next = new do_tighten_action(nactions, copyOfArray(actions,nactions), next);
243
244    next = useless_constraint_action::presolve(prob,
245                                               useless_rows, nuseless_rows,
246                                               next);
247  }
248  delete [] actions;
249  delete[]useless_rows;
250
251  if (nfixdown_cols) {
252    next = make_fixed_action::presolve(prob, fixdown_cols, nfixdown_cols,
253                                       true,
254                                       next);
255  }
256  delete[]fixdown_cols;
257
258  if (nfixup_cols) {
259    next = make_fixed_action::presolve(prob, fixup_cols, nfixup_cols,
260                                       false,
261                                       next);
262  }
263  delete[]fixup_cols;
264
265  return (next);
266}
267
268void do_tighten_action::postsolve(PostsolveMatrix *prob) const
269{
270  const action *const actions = actions_;
271  const int nactions    = nactions_;
272
273  double *colels        = prob->colels_;
274  int *hrow             = prob->hrow_;
275  CoinBigIndex *mcstrt          = prob->mcstrt_;
276  int *hincol           = prob->hincol_;
277  int *link             = prob->link_;
278  //  int ncols         = prob->ncols_;
279
280  double *clo   = prob->clo_;
281  double *cup   = prob->cup_;
282  double *rlo   = prob->rlo_;
283  double *rup   = prob->rup_;
284
285  double *sol   = prob->sol_;
286  //  double *dcost     = prob->cost_;
287  //  double *rcosts    = prob->rcosts_;
288
289  double *acts  = prob->acts_;
290  //  double *rowduals = prob->rowduals_;
291
292
293  //  const double ztolzb       = prob->ztolzb_;
294
295  //  char *cdone       = prob->cdone_;
296  //  char *rdone       = prob->rdone_;
297
298  for (const action *f = &actions[nactions-1]; actions<=f; f--) {
299    int jcol = f->col;
300    //    int iflag = f->direction;
301    int nr   = f->nrows;
302    const int *rows = f->rows;
303    const double *lbound = f->lbound;
304    const double *ubound = f->ubound;
305
306    PRESOLVEASSERT(prob->getColumnStatus(jcol)!=PrePostsolveMatrix::basic);
307    int i;
308    for (i=0;i<nr; ++i) {
309      int irow = rows[i];
310
311      rlo[irow] = lbound[i];
312      rup[irow] = ubound[i];
313
314     PRESOLVEASSERT(prob->getRowStatus(irow)==PrePostsolveMatrix::basic);
315    }
316
317    // We have just tightened the row bounds.
318    // That means we'll have to compute a new value
319    // for this variable that will satisfy everybody.
320    // We are supposed to be in a position where this
321    // is always possible.
322
323    // Each constraint has exactly one bound.
324    // The correction should only ever be forced to move in one direction.
325    //    double orig_sol = sol[jcol];
326    double correction = 0.0;
327   
328    int last_corrected = -1;
329    CoinBigIndex k = mcstrt[jcol];
330    int nk = hincol[jcol];
331    for (i=0; i<nk; ++i) {
332      int irow = hrow[k];
333      double coeff = colels[k];
334      k = link[k];
335      double newrlo = rlo[irow];
336      double newrup = rup[irow];
337      double activity = acts[irow];
338
339      if (activity + correction * coeff < newrlo) {
340        // only one of these two should fire
341        PRESOLVEASSERT( ! (activity + correction * coeff > newrup) );
342
343        last_corrected = irow;
344
345        // adjust to just meet newrlo (solve for correction)
346        double new_correction = (newrlo - activity) / coeff;
347        correction = new_correction;
348      } else if (activity + correction * coeff > newrup) {
349        last_corrected = irow;
350
351        double new_correction = (newrup - activity) / coeff;
352        correction = new_correction;
353      }
354    }
355
356    sol[jcol] += correction;
357
358    // by construction, the last row corrected (if there was one)
359    // must be at its bound, so it can be non-basic.
360    // All other rows may not be at a bound (but may if the difference
361    // is very small, causing a new correction by a tiny amount).
362
363    // now adjust the activities
364    k = mcstrt[jcol];
365    for (i=0; i<nk; ++i) {
366      int irow = hrow[k];
367      double coeff = colels[k];
368      k = link[k];
369      //      double activity = acts[irow];
370
371      acts[irow] += correction * coeff;
372
373    }
374
375    // if the col happens to get pushed to its bound,
376    // we may as well leave it non-basic.
377    if (fabs(sol[jcol] - clo[jcol]) > ZTOLDP &&
378        fabs(sol[jcol] - cup[jcol]) > ZTOLDP) {
379
380      PRESOLVEASSERT(last_corrected != -1);
381      prob->setRowStatus(last_corrected,PrePostsolveMatrix::atLowerBound);
382      prob->setColumnStatus(jcol,PrePostsolveMatrix::basic);
383    }
384  }
385}
Note: See TracBrowser for help on using the repository browser.