source: trunk/PresolveImpliedFree.cpp @ 58

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

Too agressive with costed slacks

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.2 KB
Line 
1#include <stdio.h>
2#include <math.h>
3
4#include "PresolveMatrix.hpp"
5#include "PresolveSubst.hpp"
6#include "PresolveIsolated.hpp"
7#include "PresolveImpliedFree.hpp"
8#include "ClpMessage.hpp"
9
10
11
12// If there is a row with a singleton column such that no matter what
13// the values of the other variables are, the constraint forces the singleton
14// column to have a feasible value, then we can drop the column and row,
15// since we just compute the value of the column from the row in postsolve.
16// This seems vaguely similar to the case of a useless constraint, but it
17// is different.  For example, if the singleton column is already free,
18// then this operation will eliminate it, but the constraint is not useless
19// (assuming the constraint is not trivial), since the variables do not imply an
20// upper or lower bound.
21//
22// If the column is not singleton, we can still do something similar if the
23// constraint is an equality constraint.  In that case, we substitute away
24// the variable in the other constraints it appears in.  This introduces
25// new coefficients, but the total number of coefficients never increases
26// if the column has only two constraints, and may not increase much even
27// if there are more.
28//
29// There is nothing to prevent us from substituting away a variable
30// in an equality from the other constraints it appears in, but since
31// that causes fill-in, it wouldn't make sense unless we could then
32// drop the equality itself.  We can't do that if the bounds on the
33// variable in equation aren't implied by the equality.
34// Another way of thinking of this is that there is nothing special
35// about an equality; just like one can't always drop an inequality constraint
36// with a column singleton, one can't always drop an equality.
37//
38// It is possible for two singleton columns to be in the same row.
39// In that case, the other one will become empty.  If it's bounds and
40// costs aren't just right, this signals an unbounded problem.
41// We don't need to check that specially here.
42//
43// invariant:  loosely packed
44const PresolveAction *implied_free_action::presolve(PresolveMatrix *prob,
45                                                     const PresolveAction *next,
46                                                    int & fill_level)
47{
48  double *colels        = prob->colels_;
49  int *hrow     = prob->hrow_;
50  const CoinBigIndex *mcstrt    = prob->mcstrt_;
51  int *hincol   = prob->hincol_;
52  const int ncols       = prob->ncols_;
53
54  const double *clo     = prob->clo_;
55  const double *cup     = prob->cup_;
56
57  const double *rowels  = prob->rowels_;
58  const int *hcol       = prob->hcol_;
59  const CoinBigIndex *mrstrt    = prob->mrstrt_;
60  int *hinrow   = prob->hinrow_;
61  //  const int nrows   = prob->nrows_;
62
63  /*const*/ double *rlo = prob->rlo_;
64  /*const*/ double *rup = prob->rup_;
65
66  double *cost          = prob->cost_;
67
68  presolvehlink *rlink = prob->rlink_;
69  presolvehlink *clink = prob->clink_;
70
71  const char *integerType = prob->integerType_;
72
73  const double tol = prob->feasibilityTolerance_;
74
75  //  int nbounds = 0;
76
77  action *actions       = new action [ncols];
78  int nactions = 0;
79
80  char *implied_free = new char[ncols];
81  bzero(implied_free, ncols*sizeof(char));
82
83  double *ilbound = new double[ncols];
84  double *iubound = new double[ncols];
85
86  double *tclo = new double[ncols];
87  double *tcup = new double[ncols];
88
89#if     PRESOLVE_TRY_TIGHTEN
90  for (int j=0; j<ncols; j++) {
91    ilbound[j] = -PRESOLVE_INF;
92    iubound[j] =  PRESOLVE_INF;
93    tclo[j] = clo[j];
94    tcup[j] = cup[j];
95  }
96
97  {
98    int ntightened;
99    do {
100      implied_bounds1(rowels, mrstrt, hrow, hinrow, clo, cup, hcol, ncols,
101                      rlo, rup, integerType, nrows,
102                      ilbound, iubound);
103
104      ntightened = 0;
105      for (int j=0; j<ncols; j++) {
106        if (tclo[j] < ilbound[j]) {
107          tclo[j] = ilbound[j];
108          ntightened++;
109        }
110        if (tcup[j] > iubound[j]) {
111          tcup[j] = iubound[j];
112          ntightened++;
113        }
114      }
115      printf("NTIGHT: %d\n", ntightened);
116    } while (ntightened);
117  }
118#endif
119
120  int numberLook = prob->numberColsToDo_;
121  int iLook;
122  int * look = prob->colsToDo_;
123  int * look2 = NULL;
124  // if gone from 2 to 3 look at all
125  if (fill_level<0) {
126    look2 = new int[ncols];
127    look=look2;
128    for (iLook=0;iLook<ncols;iLook++) 
129      look[iLook]=iLook;
130    numberLook=ncols;
131 }
132
133  for (iLook=0;iLook<numberLook;iLook++) {
134    int j=look[iLook];
135    if ((hincol[j] >= 1 && hincol[j] <= 3) &&
136        !integerType[j]) {
137      CoinBigIndex kcs = mcstrt[j];
138      CoinBigIndex kce = kcs + hincol[j];
139
140      for (CoinBigIndex k=kcs; k<kce; ++k) {
141        int row = hrow[k];
142        double coeffj = colels[k];
143
144        // if its row is an equality constraint...
145        if (hinrow[row] > 1 &&  // don't bother with singleton rows
146
147            // either this is a singleton col,
148            // (for momement - only if no cost)
149            // or this particular row is an equality
150            ((hincol[j] == 1 && !cost[j]) || fabs(rlo[row] - rup[row]) < tol) &&
151
152            fabs(coeffj) > ZTOLDP) {
153
154          CoinBigIndex krs = mrstrt[row];
155          CoinBigIndex kre = krs + hinrow[row];
156
157          double maxup, maxdown, ilow, iup;
158          implied_bounds(rowels, clo, cup, hcol,
159                         krs, kre,
160                         &maxup, &maxdown,
161                         j, rlo[row], rup[row], &ilow, &iup);
162
163#if     PRESOLVE_TRY_TIGHTEN
164          if ((clo[j] <= ilbound[j] && iubound[j] <= cup[j]) &&
165              ! (clo[j] <= ilow && iup <= cup[j]))
166            printf("TIGHTER:  %6d %6d\n", row, j);
167#endif
168
169          if (maxup < PRESOLVE_INF && maxup + tol < rlo[row]) {
170            /* there is an upper bound and it can't be reached */
171            prob->status_|= 1;
172            prob->messageHandler()->message(CLP_PRESOLVE_ROWINFEAS,
173                                             prob->messages())
174                                               <<row
175                                               <<rlo[row]
176                                               <<rup[row]
177                                               <<CoinMessageEol;
178            break;
179          } else if (-PRESOLVE_INF < maxdown && rup[row] < maxdown - tol) {
180            /* there is a lower bound and it can't be reached */
181            prob->status_|= 1;
182            prob->messageHandler()->message(CLP_PRESOLVE_ROWINFEAS,
183                                             prob->messages())
184                                               <<row
185                                               <<rlo[row]
186                                               <<rup[row]
187                                               <<CoinMessageEol;
188            break;
189          } else if (clo[j] <= ilow && iup <= cup[j]) {
190
191            // both column bounds implied by the constraints of the problem
192            implied_free[j] = hincol[j];
193            break;
194          }
195        }
196      }
197    }
198  }
199  // implied_free[j] == hincol[j] && hincol[j] > 0 ==> j is implied free
200
201#if 0
202  // DEBUG
203  static int nfree = 0;
204  static int maxfree = atoi(getenv("MAXFREE"));
205#endif
206
207  int isolated_row = -1;
208
209  // first pick off the easy ones
210  // note that this will only deal with columns that were originally
211  // singleton; it will not deal with doubleton columns that become
212  // singletons as a result of dropping rows.
213  for (iLook=0;iLook<numberLook;iLook++) {
214    int j=look[iLook];
215    if (hincol[j] == 1 && implied_free[j] == 1) {
216      CoinBigIndex kcs = mcstrt[j];
217      int row = hrow[kcs];
218      double coeffj = colels[kcs];
219
220      CoinBigIndex krs = mrstrt[row];
221      CoinBigIndex kre = krs + hinrow[row];
222
223#if 0
224      if (nfree >= maxfree)
225        continue;
226      nfree++;
227#endif
228
229      // isolated rows are weird
230      {
231        int n = 0;
232        for (CoinBigIndex k=krs; k<kre; ++k)
233          n += hincol[hcol[k]];
234        if (n==hinrow[row]) {
235          isolated_row = row;
236          break;
237        }
238      }
239
240      const bool nonzero_cost = (cost[j] != 0.0&&fabs(rup[row]-rlo[row])<=tol);
241
242      double *save_costs = nonzero_cost ? new double[hinrow[row]] : NULL;
243
244      {
245        action *s = &actions[nactions++];
246             
247        s->row = row;
248        s->col = j;
249
250        s->clo = clo[j];
251        s->cup = cup[j];
252        s->rlo = rlo[row];
253        s->rup = rup[row];
254
255        s->ninrow = hinrow[row];
256        s->rowels = presolve_duparray(&rowels[krs], hinrow[row]);
257        s->rowcols = presolve_duparray(&hcol[krs], hinrow[row]);
258        s->costs = save_costs;
259      }
260
261      if (nonzero_cost) {
262        double rhs = rlo[row];
263        double costj = cost[j];
264
265#if     DEBUG_PRESOLVE
266        printf("FREE COSTS:  %g  ", costj);
267#endif
268        for (CoinBigIndex k=krs; k<kre; k++) {
269          int jcol = hcol[k];
270          save_costs[k-krs] = cost[jcol];
271
272          if (jcol != j) {
273            double coeff = rowels[k];
274
275#if     DEBUG_PRESOLVE
276            printf("%g %g   ", cost[jcol], coeff/coeffj);
277#endif
278            /*
279             * Similar to eliminating doubleton:
280             *   cost1 x = cost1 (c - b y) / a = (c cost1)/a - (b cost1)/a
281             *   cost[icoly] += cost[icolx] * (-coeff2 / coeff1);
282             */
283            cost[jcol] += costj * (-coeff / coeffj);
284          }
285        }
286#if     DEBUG_PRESOLVE
287        printf("\n");
288
289        /* similar to doubleton */
290        printf("BIAS??????? %g %g %g %g\n",
291               costj * rhs / coeffj,
292               costj, rhs, coeffj);
293#endif
294        prob->change_bias(costj * rhs / coeffj);
295        // ??
296        cost[j] = 0.0;
297      }
298
299      /* remove the row from the columns in the row */
300      for (CoinBigIndex k=krs; k<kre; k++) {
301        int jcol=hcol[k];
302        prob->addCol(jcol);
303        presolve_delete_from_row(jcol, row, mcstrt, hincol, hrow, colels);
304      }
305      PRESOLVE_REMOVE_LINK(rlink, row);
306      hinrow[row] = 0;
307
308      // just to make things squeeky
309      rlo[row] = 0.0;
310      rup[row] = 0.0;
311
312      PRESOLVE_REMOVE_LINK(clink, j);
313      hincol[j] = 0;
314
315      implied_free[j] = 0;      // probably unnecessary
316    }
317  }
318
319  delete [] look2;
320  if (nactions) {
321#if     PRESOLVE_SUMMARY
322    printf("NIMPLIED FREE:  %d\n", nactions);
323#endif
324    next = new implied_free_action(nactions, copyOfArray(actions,nactions), next);
325  }
326  deleteAction(actions);
327
328  delete[]ilbound;
329  delete[]iubound;
330  delete[]tclo;
331  delete[]tcup;
332
333  if (isolated_row != -1) {
334    const PresolveAction *nextX = isolated_constraint_action::presolve(prob, 
335                                                isolated_row, next);
336    if (nextX)
337      next = nextX; // may fail
338  }
339
340  // try more complex ones
341  if (fill_level)
342    next = subst_constraint_action::presolve(prob, implied_free, next,fill_level);
343  delete[]implied_free;
344
345  return (next);
346}
347
348
349
350const char *implied_free_action::name() const
351{
352  return ("implied_free_action");
353}
354
355void implied_free_action::postsolve(PostsolveMatrix *prob) const
356{
357  const action *const actions = actions_;
358  const int nactions = nactions_;
359
360  double *colels        = prob->colels_;
361  int *hrow             = prob->hrow_;
362  CoinBigIndex *mcstrt          = prob->mcstrt_;
363  int *hincol           = prob->hincol_;
364  int *link             = prob->link_;
365
366  double *clo   = prob->clo_;
367  double *cup   = prob->cup_;
368
369  double *rlo   = prob->rlo_;
370  double *rup   = prob->rup_;
371
372  double *sol   = prob->sol_;
373
374  double *rcosts        = prob->rcosts_;
375  double *dcost         = prob->cost_;
376
377  double *acts  = prob->acts_;
378  double *rowduals = prob->rowduals_;
379
380  //  const double ztoldj       = prob->ztoldj_;
381  //  const double ztolzb       = prob->ztolzb_;
382
383  const double maxmin   = prob->maxmin_;
384
385  char *cdone   = prob->cdone_;
386  char *rdone   = prob->rdone_;
387  CoinBigIndex free_list = prob->free_list_;
388
389  for (const action *f = &actions[nactions-1]; actions<=f; f--) {
390
391    int irow = f->row;
392    int icol = f->col;
393         
394    int ninrow = f->ninrow;
395    const double *rowels = f->rowels;
396    const int *rowcols = f->rowcols;
397    const double *save_costs = f->costs;
398
399    // put back coefficients in the row
400    // this includes recreating the singleton column
401    {
402      for (int k = 0; k<ninrow; k++) {
403        int jcol = rowcols[k];
404        double coeff = rowels[k];
405
406        if (save_costs)
407          dcost[jcol] = save_costs[k];
408
409        {
410          CoinBigIndex kk = free_list;
411          free_list = link[free_list];
412
413          check_free_list(free_list);
414
415          link[kk] = mcstrt[jcol];
416          mcstrt[jcol] = kk;
417          colels[kk] = coeff;
418          hrow[kk] = irow;
419        }
420
421        if (jcol == icol) {
422          // initialize the singleton column
423          hincol[jcol] = 1;
424          clo[icol] = f->clo;
425          cup[icol] = f->cup;
426
427          cdone[icol] = IMPLIED_FREE;
428        } else {
429          hincol[jcol]++;
430        }
431      }
432      rdone[irow] = IMPLIED_FREE;
433
434      rlo[irow] = f->rlo;
435      rup[irow] = f->rup;
436    }
437    delete [] save_costs;
438    // coeff has now been initialized
439
440    // compute solution
441    {
442      double act = 0.0;
443      double coeff = 0.0;
444
445      for (int k = 0; k<ninrow; k++)
446        if (rowcols[k] == icol)
447          coeff = rowels[k];
448        else {
449          int jcol = rowcols[k];
450          PRESOLVE_STMT(CoinBigIndex kk = presolve_find_row2(irow, mcstrt[jcol], hincol[jcol], hrow, link));
451          act += rowels[k] * sol[jcol];
452        }
453           
454      PRESOLVEASSERT(fabs(coeff) > ZTOLDP);
455      // choose rowdual to make this col basic
456      rowduals[irow] = maxmin*dcost[icol] / coeff;
457      rcosts[icol] = 0.0;
458      if ((rlo[irow] < rup[irow] && rowduals[irow] < 0.0)
459          || rlo[irow]< -1.0e20) {
460        if (rlo[irow]<-1.0e20&&rowduals[irow]>=0.0)
461          printf("IMP %g %g %g\n",rlo[irow],rup[irow],rowduals[irow]);
462        sol[icol] = (rup[irow] - act) / coeff;
463        acts[irow] = rup[irow];
464      } else {
465        sol[icol] = (rlo[irow] - act) / coeff;
466        acts[irow] = rlo[irow];
467      }
468
469
470      prob->setRowStatus(irow,PrePostsolveMatrix::atLowerBound);
471      prob->setColumnStatus(icol,PrePostsolveMatrix::basic);
472    }
473  }
474  prob->free_list_ = free_list;
475}
Note: See TracBrowser for help on using the repository browser.