source: branches/devel-1/PresolveImpliedFree.cpp @ 33

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

Presolve in as option

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.0 KB
RevLine 
[29]1#include <stdio.h>
2#include <math.h>
3#include <strings.h>
4
5#include "PresolveMatrix.hpp"
6#include "PresolveSubst.hpp"
7#include "PresolveIsolated.hpp"
8#include "PresolveImpliedFree.hpp"
9#include "ClpMessage.hpp"
10
11
12
13// If there is a row with a singleton column such that no matter what
14// the values of the other variables are, the constraint forces the singleton
15// column to have a feasible value, then we can drop the column and row,
16// since we just compute the value of the column from the row in postsolve.
17// This seems vaguely similar to the case of a useless constraint, but it
18// is different.  For example, if the singleton column is already free,
19// then this operation will eliminate it, but the constraint is not useless
20// (assuming the constraint is not trivial), since the variables do not imply an
21// upper or lower bound.
22//
23// If the column is not singleton, we can still do something similar if the
24// constraint is an equality constraint.  In that case, we substitute away
25// the variable in the other constraints it appears in.  This introduces
26// new coefficients, but the total number of coefficients never increases
27// if the column has only two constraints, and may not increase much even
28// if there are more.
29//
30// There is nothing to prevent us from substituting away a variable
31// in an equality from the other constraints it appears in, but since
32// that causes fill-in, it wouldn't make sense unless we could then
33// drop the equality itself.  We can't do that if the bounds on the
34// variable in equation aren't implied by the equality.
35// Another way of thinking of this is that there is nothing special
36// about an equality; just like one can't always drop an inequality constraint
37// with a column singleton, one can't always drop an equality.
38//
39// It is possible for two singleton columns to be in the same row.
40// In that case, the other one will become empty.  If it's bounds and
41// costs aren't just right, this signals an unbounded problem.
42// We don't need to check that specially here.
43//
44// invariant:  loosely packed
45const PresolveAction *implied_free_action::presolve(PresolveMatrix *prob,
46                                                     const PresolveAction *next,
47                                                    int & fill_level)
48{
49  double *colels        = prob->colels_;
50  int *hrow     = prob->hrow_;
51  const int *mcstrt     = prob->mcstrt_;
52  int *hincol   = prob->hincol_;
53  const int ncols       = prob->ncols_;
54
55  const double *clo     = prob->clo_;
56  const double *cup     = prob->cup_;
57
58  const double *rowels  = prob->rowels_;
59  const int *hcol       = prob->hcol_;
60  const int *mrstrt     = prob->mrstrt_;
61  int *hinrow   = prob->hinrow_;
62  const int nrows       = prob->nrows_;
63
64  /*const*/ double *rlo = prob->rlo_;
65  /*const*/ double *rup = prob->rup_;
66
67  double *cost          = prob->cost_;
68
69  presolvehlink *rlink = prob->rlink_;
70  presolvehlink *clink = prob->clink_;
71
72  const char *integerType = prob->integerType_;
73
74  const double tol = ZTOLDP;
75
76  int nbounds = 0;
77
78  action *actions       = new action [ncols];
79  int nactions = 0;
80
81  char *implied_free = new char[ncols];
82  bzero(implied_free, ncols*sizeof(char));
83
84  double *ilbound = new double[ncols];
85  double *iubound = new double[ncols];
86
87  double *tclo = new double[ncols];
88  double *tcup = new double[ncols];
89
90#if     PRESOLVE_TRY_TIGHTEN
91  for (int j=0; j<ncols; j++) {
92    ilbound[j] = -PRESOLVE_INF;
93    iubound[j] =  PRESOLVE_INF;
94    tclo[j] = clo[j];
95    tcup[j] = cup[j];
96  }
97
98  {
99    int ntightened;
100    do {
101      implied_bounds1(rowels, mrstrt, hrow, hinrow, clo, cup, hcol, ncols,
102                      rlo, rup, integerType, nrows,
103                      ilbound, iubound);
104
105      ntightened = 0;
106      for (int j=0; j<ncols; j++) {
107        if (tclo[j] < ilbound[j]) {
108          tclo[j] = ilbound[j];
109          ntightened++;
110        }
111        if (tcup[j] > iubound[j]) {
112          tcup[j] = iubound[j];
113          ntightened++;
114        }
115      }
116      printf("NTIGHT: %d\n", ntightened);
117    } while (ntightened);
118  }
119#endif
120
121  int numberLook = prob->numberColsToDo_;
122  int iLook;
123  int * look = prob->colsToDo_;
124  int * look2 = NULL;
125  // if gone from 2 to 3 look at all
126  if (fill_level<0) {
127    look2 = new int[ncols];
128    look=look2;
129    for (iLook=0;iLook<ncols;iLook++) 
130      look[iLook]=iLook;
131    numberLook=ncols;
132 }
133
134  for (iLook=0;iLook<numberLook;iLook++) {
135    int j=look[iLook];
136    if ((hincol[j] >= 1 && hincol[j] <= 3) &&
137        !integerType[j]) {
138      int kcs = mcstrt[j];
139      int kce = kcs + hincol[j];
140
141      for (int k=kcs; k<kce; ++k) {
142        int row = hrow[k];
143        double coeffj = colels[k];
144
145        // if its row is an equality constraint...
146        if (hinrow[row] > 1 &&  // don't bother with singleton rows
147
148            // either this is a singleton col,
149            // or this particular row is an equality
150            (hincol[j] == 1 || fabs(rlo[row] - rup[row]) < tol) &&
151
152            fabs(coeffj) > ZTOLDP) {
153
154          int krs = mrstrt[row];
155          int 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->originalModel_->messageHandler()->message(CLP_PRESOLVE_ROWINFEAS,
173                                             prob->originalModel_->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->originalModel_->messageHandler()->message(CLP_PRESOLVE_ROWINFEAS,
183                                             prob->originalModel_->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      int kcs = mcstrt[j];
217      int row = hrow[kcs];
218      double coeffj = colels[kcs];
219
220      int krs = mrstrt[row];
221      int 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 (int 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);
241
[33]242      double *save_costs = nonzero_cost ? new double[hinrow[row]] : NULL;
[29]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 (int 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
[31]294        prob->change_bias(costj * rhs / coeffj);
[29]295        // ??
296        cost[j] = 0.0;
297      }
298
299      /* remove the row from the columns in the row */
300      for (int 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  delete [] actions;
327
328  delete[]ilbound;
329  delete[]iubound;
330  delete[]tclo;
331  delete[]tcup;
332
333  if (isolated_row != -1)
334    next = isolated_constraint_action::presolve(prob, isolated_row, next);
335
336  // try more complex ones
337  if (fill_level)
338    next = subst_constraint_action::presolve(prob, implied_free, next,fill_level);
339  delete[]implied_free;
340
341  return (next);
342}
343
344
345
346const char *implied_free_action::name() const
347{
348  return ("implied_free_action");
349}
350
351void implied_free_action::postsolve(PostsolveMatrix *prob) const
352{
353  const action *const actions = actions_;
354  const int nactions = nactions_;
355
356  double *colels        = prob->colels_;
357  int *hrow             = prob->hrow_;
358  int *mcstrt           = prob->mcstrt_;
359  int *hincol           = prob->hincol_;
360  int *link             = prob->link_;
361
362  double *clo   = prob->clo_;
363  double *cup   = prob->cup_;
364
365  double *rlo   = prob->rlo_;
366  double *rup   = prob->rup_;
367
368  double *sol   = prob->sol_;
369
370  double *rcosts        = prob->rcosts_;
371  double *dcost         = prob->cost_;
372
373  double *acts  = prob->acts_;
374  double *rowduals = prob->rowduals_;
375
376  const double ztoldj   = prob->ztoldj_;
377  const double ztolzb   = prob->ztolzb_;
378
379  const double maxmin   = prob->maxmin_;
380
381  char *cdone   = prob->cdone_;
382  char *rdone   = prob->rdone_;
383  int free_list = prob->free_list_;
384
385  for (const action *f = &actions[nactions-1]; actions<=f; f--) {
386
387    int irow = f->row;
388    int icol = f->col;
389         
390    int ninrow = f->ninrow;
391    const double *rowels = f->rowels;
392    const int *rowcols = f->rowcols;
393    const double *save_costs = f->costs;
394
395    // put back coefficients in the row
396    // this includes recreating the singleton column
397    {
398      for (int k = 0; k<ninrow; k++) {
399        int jcol = rowcols[k];
400        double coeff = rowels[k];
401
402        if (save_costs)
403          dcost[jcol] = save_costs[k];
404
405        {
406          int kk = free_list;
407          free_list = link[free_list];
408
409          check_free_list(free_list);
410
411          link[kk] = mcstrt[jcol];
412          mcstrt[jcol] = kk;
413          colels[kk] = coeff;
414          hrow[kk] = irow;
415        }
416
417        if (jcol == icol) {
418          // initialize the singleton column
419          hincol[jcol] = 1;
420          clo[icol] = f->clo;
421          cup[icol] = f->cup;
422
423          cdone[icol] = IMPLIED_FREE;
424        } else {
425          hincol[jcol]++;
426        }
427      }
428      rdone[irow] = IMPLIED_FREE;
429
430      rlo[irow] = f->rlo;
431      rup[irow] = f->rup;
432    }
[33]433    delete [] save_costs;
[29]434    // coeff has now been initialized
435
436    // compute solution
437    {
438      double act = 0.0;
439      double coeff;
440
441      for (int k = 0; k<ninrow; k++)
442        if (rowcols[k] == icol)
443          coeff = rowels[k];
444        else {
445          int jcol = rowcols[k];
446          PRESOLVE_STMT(int kk = presolve_find_row2(irow, mcstrt[jcol], hincol[jcol], hrow, link));
447          PRESOLVEASSERT(colels[kk] == rowels[k]);
448          act += rowels[k] * sol[jcol];
449        }
450           
451      PRESOLVEASSERT(fabs(coeff) > ZTOLDP);
[31]452      // choose rowdual to make this col basic
453      rowduals[irow] = maxmin*dcost[icol] / coeff;
454      rcosts[icol] = 0.0;
455      if ((rlo[irow] < rup[irow] && rowduals[irow] < 0.0)
456          || rlo[irow]< -1.0e20) {
457        if (rlo[irow]<-1.0e20&&rowduals[irow]>=0.0)
458          printf("IMP %g %g %g\n",rlo[irow],rup[irow],rowduals[irow]);
[29]459        sol[icol] = (rup[irow] - act) / coeff;
460        acts[irow] = rup[irow];
461      } else {
462        sol[icol] = (rlo[irow] - act) / coeff;
463        acts[irow] = rlo[irow];
464      }
465
466
467      prob->setRowStatus(irow,PrePostsolveMatrix::atLowerBound);
468      prob->setColumnStatus(icol,PrePostsolveMatrix::basic);
469    }
470  }
471  prob->free_list_ = free_list;
472}
Note: See TracBrowser for help on using the repository browser.