source: trunk/PresolveTighten.cpp @ 56

Last change on this file since 56 was 56, checked in by forrest, 18 years ago

Idiot and modifications for Visual C++

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