source: tags/release_2002-11-05/PresolveTighten.cpp @ 778

Last change on this file since 778 was 50, checked in by ladanyi, 17 years ago

devel-1 merged into HEAD

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