source: tags/preroot/PresolveDupcol.cpp @ 778

Last change on this file since 778 was 173, checked in by forrest, 16 years ago

Improvements to presolve

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.5 KB
Line 
1#include <stdio.h>
2#include <math.h>
3
4#include "PresolveMatrix.hpp"
5#include "PresolveFixed.hpp"
6#include "PresolveDupcol.hpp"
7#include "CoinSort.hpp"
8#include "PresolveUseless.hpp"
9#include "ClpMessage.hpp"
10
11#define DSEED2 2147483647.0
12
13const char *dupcol_action::name() const
14{
15  return ("dupcol_action");
16}
17
18const char *duprow_action::name() const
19{
20  return ("duprow_action");
21}
22
23void init_random_vec(double *work, int n)
24{
25  double deseed = 12345678.0;
26
27  for (int i = 0; i < n; ++i) {
28    deseed *= 16807.;
29    int jseed = (int) (deseed /    DSEED2);
30    deseed -= (double) jseed * DSEED2;
31    double random = deseed /  DSEED2;
32
33    work[i]=random;
34  }
35}
36
37int compute_sums(const int *len, const CoinBigIndex *starts,
38                 /*const*/ int *index, /*const*/ double *elems, // index/elems are sorted
39                 const double *work,
40                 double *sums, int *sorted, int n)
41{
42  int nlook=0;
43  for (int i = 0; i < n; ++i)
44    if (len[i] > 0 /*1?*/) {
45      CoinBigIndex kcs = starts[i];
46      CoinBigIndex kce = kcs + len[i];
47
48      double value=0.0;
49
50      // sort all columns for easy comparison in next loop
51      CoinSort_2(index+kcs,index+kcs+len[i],elems+kcs);
52      //ekk_sort2(index+kcs, elems+kcs, len[i]);
53
54      for (CoinBigIndex k=kcs;k<kce;k++) {
55        int irow=index[k];
56        value += work[irow]*elems[k];
57      }
58      sums[nlook]=value;
59      sorted[nlook]=i;
60      ++nlook;
61    }
62  return (nlook);
63}
64
65dupcol_action::dupcol_action(int nactions,
66                const action *actions,
67                const PresolveAction *next) :
68    PresolveAction(next),
69    nactions_(nactions), actions_(actions)
70{
71}
72
73
74// This is just ekkredc5, adapted into the new framework.
75// the datasets scorpion.mps and allgrade.mps have duplicate columns.
76const PresolveAction *dupcol_action::presolve(PresolveMatrix *prob,
77                                               const PresolveAction *next)
78{
79  double *colels        = prob->colels_;
80  int *hrow             = prob->hrow_;
81  CoinBigIndex *mcstrt          = prob->mcstrt_;
82  int *hincol           = prob->hincol_;
83  int ncols             = prob->ncols_;
84
85  double *clo   = prob->clo_;
86  double *cup   = prob->cup_;
87  double * sol = prob->sol_;
88
89  double *rowels        = prob->rowels_;
90  int *hcol     = prob->hcol_;
91  const CoinBigIndex *mrstrt    = prob->mrstrt_;
92  int *hinrow           = prob->hinrow_;
93  int nrows             = prob->nrows_;
94
95  //  double *rlo       = prob->rlo_;
96  //  double *rup       = prob->rup_;
97
98  double *dcost = prob->cost_;
99
100  const char *integerType = prob->integerType_;
101
102  double maxmin = prob->maxmin_;
103
104  action *actions       = new action [ncols];
105  int nactions = 0;
106
107  int *fixed_down       = new int[ncols];
108  int nfixed_down       = 0;
109  int *fixed_up         = new int[ncols];
110  int nfixed_up         = 0;
111
112  double * workcol = new double[ncols];
113  double * workrow = new double[nrows];
114  int * sort = new int[ncols];
115
116  // initialize workrow to have "random" values
117  init_random_vec(workrow, nrows);
118
119  // If we sum all the coefficients of a column,
120  // then duplicate columns must have the same sum.
121  // However, other columns may also.
122  // To weed most of them out, we multiply the coefficients
123  // by the random value associated with the row.
124  int nlook = compute_sums(hincol, mcstrt, hrow, colels, workrow, workcol, sort, ncols);
125
126#if 1
127  // not tested yet
128  if (maxmin < 0.0)
129    return (next);
130#endif
131
132  CoinSort_2(workcol,workcol+nlook,sort);
133  //ekk_sortonDouble(workcol,sort,nlook);
134
135#if 0
136  // It may be the case that several columns are duplicate.
137  // If not all have the same cost, then we have to make sure
138  // that we set the most expensive one to its minimum
139  // now sort in each class by cost
140  {
141    double dval = workcol[0];
142    int first = 0;
143    for (int jj = 1; jj < nlook; jj++) {
144      while (workcol[jj]==dval)
145        jj++;
146
147      if (first + 1 < jj) {
148        double buf[jj - first];
149        for (int i=first; i<jj; ++i)
150          buf[i-first] = dcost[sort[i]]*maxmin;
151
152        CoinSort_2(buf,buf+jj-first,sort+first);
153        //ekk_sortonDouble(buf,&sort[first],jj-first);
154      }
155    }
156  }
157#endif
158
159  // it appears to be the case that this loop is finished,
160  // there may still be duplicate cols left.
161  // I haven't done anything about that yet.
162  for (int jj = 1; jj < nlook; jj++) {
163    if (workcol[jj]==workcol[jj-1]) {
164      int ithis=sort[jj];
165      int ilast=sort[jj-1];
166      CoinBigIndex kcs = mcstrt[ithis];
167      CoinBigIndex kce = kcs + hincol[ithis];
168
169      if (hincol[ithis] == hincol[ilast]&&!prob->colProhibited2(ithis)) {
170        int ishift = mcstrt[ilast] - kcs;
171        CoinBigIndex k;
172        for (k=kcs;k<kce;k++) {
173          if (hrow[k] != hrow[k+ishift] ||
174              colels[k] != colels[k+ishift]) {
175            break;
176          }
177        }
178        if (k == kce) {
179          // these really are duplicate columns
180
181          /* now check bounds and costs to see what is what */
182          double clo1=clo[ilast];
183          double cup1=cup[ilast];
184          double clo2=clo[ithis];
185          double cup2=cup[ithis];
186          double dcost1=dcost[ilast]*maxmin;
187          double dcost2=dcost[ithis]*maxmin;
188          double newSolution = sol[ilast]+sol[ithis];
189          //int itype=27;
190          if (clo1==cup1||clo2==cup2)
191            abort();
192
193          if (/*ddjs[ilast]||ddjs[ithis]*/
194              integerType[ilast] || integerType[ithis]) {
195#ifdef PRINT_DEBUG
196            printf("at least one of %d or %d integer\n",ilast-nrowsmx,
197                   ithis-nrowsmx);
198#endif
199            continue;
200          } else if (dcost1==dcost2) {
201            //
202            // SAME COSTS
203            //
204            double l_j = clo[ithis];
205            double u_j = cup[ithis];
206            double l_k = clo[ilast];
207            double u_k = cup[ilast];
208
209            if (! (l_j + u_k <= l_k + u_j)) {
210              swap(ilast, ithis);
211              swap(clo1, clo2);
212              swap(cup1, cup2);
213              //swap(dcost1, dcost2);
214            }
215
216            //PRESOLVE_STMT(printf("DUPCOL (%d,%d)\n", ithis, ilast));
217
218            /* identical cost so can add ranges */
219            clo1 += clo2;
220            cup1 += cup2;
221            if (clo1<-1.0e20) {
222              clo1=-1.0e31;
223            }
224            if (cup1>1.0e20) {
225              cup1=1.0e31;
226            }
227            /*dfix=0.0;*/
228            //itype=26;
229
230            {
231              action *s = &actions[nactions];     
232              nactions++;
233
234              s->thislo = clo[ithis];
235              s->thisup = cup[ithis];
236              s->lastlo = clo[ilast];
237              s->lastup = cup[ilast];
238              s->ithis  = ithis;
239              s->ilast  = ilast;
240
241              s->nincol = hincol[ithis];
242              s->colels = copyOfArray(&colels[mcstrt[ithis]], hincol[ithis]);
243              s->colrows= copyOfArray(&hrow[mcstrt[ithis]], hincol[ithis]);
244            }
245
246            // adjust ilast
247            clo[ilast]=clo1;
248            cup[ilast]=cup1;
249            sol[ilast] = newSolution;
250            if (prob->getColumnStatus(ilast)==PrePostsolveMatrix::basic||
251                prob->getColumnStatus(ithis)==PrePostsolveMatrix::basic)
252              prob->setColumnStatus(ilast,PrePostsolveMatrix::basic);
253
254            // delete ithis
255            dcost[ithis] = 0.0;
256            sol[ithis]=clo2;
257            {
258              CoinBigIndex kcs = mcstrt[ithis];
259              CoinBigIndex kce = kcs + hincol[ithis];
260              for (CoinBigIndex k=kcs; k<kce; ++k)
261                // delete ithis from its rows
262                presolve_delete_from_row(hrow[k], ithis, mrstrt, hinrow, hcol, rowels);
263            }
264            hincol[ithis] = 0;
265
266            /* make sure fixed one out of way */
267            sort[jj] = ilast;
268            //if (ithis==sort[jj]) {
269            //sort[jj]=sort[jj-1];
270            //}
271          } else {
272            //
273            // (dcost1!=dcost2) - DIFFERENT COSTS
274            //
275#define PRINT_DEBUG
276            /* look to see whether we can drop one ? */
277            if (clo1<-1.0e20 && clo2<-1.0e20) {
278              // both unbounded below
279              /* can't cope for now */
280#ifdef PRINT_DEBUG
281              printf("** Odd columns %d and %d\n", ilast,ithis);
282#endif
283              continue;
284            } else if (clo1<-1.0e20) {
285              printf("ILAST:  %d %d\n", ilast, ithis);
286              swap(ilast, ithis);
287              printf("ILAST1:  %d %d\n", ilast, ithis);
288              swap(clo1, clo2);
289              swap(cup1, cup2);
290              swap(dcost1, dcost2);
291            }
292            // first one (ilast) bounded below - PRESOLVEASSERT(! (clo1<-1.0e20) );
293
294            if (cup1>1.0e20) {
295              /* ilast can go to plus infinity */
296              if (clo2<-1.0e20) {
297                if (dcost1<dcost2) {
298#ifdef PRINT_DEBUG
299                  printf("** Problem seems unbounded %d and %d\n", ilast,ithis);
300#endif
301                  break;
302                } else {
303                  continue;
304                }
305              } else if (cup2>1.0e20) {
306                //PRESOLVE_STMT(printf("DUPCOL1 (%d,%d) %g\n", ithis, ilast, workcol[jj]));
307                // both have an lb, both have no ub
308                // the more expensive one would end up at its lb
309                if (dcost1>dcost2) {
310                  swap(ilast, ithis);
311                  swap(clo1, clo2);
312                  swap(cup1, cup2);
313                  swap(dcost1, dcost2);
314                }
315                sol[ilast] = newSolution;
316                if (prob->getColumnStatus(ilast)==PrePostsolveMatrix::basic||
317                    prob->getColumnStatus(ithis)==PrePostsolveMatrix::basic)
318                  prob->setColumnStatus(ilast,PrePostsolveMatrix::basic);
319               
320                sol[ithis]=clo2;
321                fixed_down[nfixed_down++] = ithis ;
322                sort[jj]                  = ilast;
323              } else {
324                /* ithis ranged - last not */
325                if (dcost2>dcost1) {
326                  //PRESOLVE_STMT(printf("DUPCOL2 (%d,%d)\n", ithis, ilast));
327                  /* set ithis to lower bound */
328                  fixed_down[nfixed_down++] = ithis;
329                  sort[jj] = ilast;
330                  sol[ilast] = newSolution;
331                  if (prob->getColumnStatus(ilast)==PrePostsolveMatrix::basic||
332                      prob->getColumnStatus(ithis)==PrePostsolveMatrix::basic)
333                    prob->setColumnStatus(ilast,PrePostsolveMatrix::basic);
334                 
335                  sol[ithis]=clo2;
336                } else {
337                  /* no good */
338                  continue;
339                }
340              }
341            } else {
342              /* ilast ranged */
343              if (clo2<-1.0e20) {
344                // ithis has no lb
345                if (dcost2>dcost1) {
346                  //PRESOLVE_STMT(printf("DUPCOL3 (%d,%d)\n", ithis, ilast));
347                  /* set ilast to upper bound */
348                  fixed_up[nfixed_up++] = ilast;
349                  sort[jj] = ithis;
350                  sol[ithis] = newSolution;
351                  if (prob->getColumnStatus(ithis)==PrePostsolveMatrix::basic||
352                      prob->getColumnStatus(ilast)==PrePostsolveMatrix::basic)
353                    prob->setColumnStatus(ithis,PrePostsolveMatrix::basic);
354                 
355                  sol[ilast]=clo1;
356                } else {
357                  /* need both */
358                  continue;
359                }
360              } else if (cup2>1.0e20) {
361                if (dcost2<dcost1) {
362                  //PRESOLVE_STMT(printf("DUPCOL4 (%d,%d)\n", ithis, ilast));
363                  /* set ilast to lower bound */
364                  //SWAP(int, ilast, ithis);
365                  fixed_down[nfixed_down++] = ilast;
366                  sort[jj] = ithis;
367                  sol[ithis] = newSolution;
368                  if (prob->getColumnStatus(ithis)==PrePostsolveMatrix::basic||
369                      prob->getColumnStatus(ilast)==PrePostsolveMatrix::basic)
370                    prob->setColumnStatus(ithis,PrePostsolveMatrix::basic);
371                 
372                  sol[ilast]=clo1;
373                } else {
374                  /* need both */
375                  continue;
376                }
377              } else {
378                /* both ranged - no hope */
379                continue;
380              }
381            }
382          }
383        }
384      }
385    }
386  }
387
388  delete[]workrow;
389  delete[]workcol;
390  delete[]sort;
391
392
393  if (nactions) {
394#if     PRESOLVE_SUMMARY
395    printf("DUPLICATE COLS:  %d\n", nactions);
396#endif
397    next = new dupcol_action(nactions, copyOfArray(actions,nactions), next);
398  }
399  deleteAction(actions,action*);
400  if (nfixed_down)
401    next = make_fixed_action::presolve(prob, fixed_down, nfixed_down,
402                                       true,
403                                       next);
404
405  if (nfixed_up)
406    next = make_fixed_action::presolve(prob, fixed_up, nfixed_up,
407                                       false,
408                                       next);
409
410  delete[]fixed_down;
411  delete[]fixed_up;
412
413  return (next);
414}
415
416void create_col(int col, int n, int *rows, double *els,
417                CoinBigIndex *mcstrt, double *colels, int *hrow, int *link,
418                CoinBigIndex *free_listp)
419{
420  CoinBigIndex free_list = *free_listp;
421  int xstart = NO_LINK;
422  for (int i=0; i<n; ++i) {
423    CoinBigIndex k = free_list;
424    free_list = link[free_list];
425
426    check_free_list(free_list);
427
428    hrow[k]   = rows[i];
429    colels[k] = els[i];
430    link[k] = xstart;
431    xstart = k;
432  }
433  mcstrt[col] = xstart;
434  *free_listp = free_list;
435}
436
437
438void dupcol_action::postsolve(PostsolveMatrix *prob) const
439{
440  const action *const actions = actions_;
441  const int nactions = nactions_;
442
443  double *clo   = prob->clo_;
444  double *cup   = prob->cup_;
445
446  double *sol   = prob->sol_;
447  double *dcost = prob->cost_;
448 
449  double *colels        = prob->colels_;
450  int *hrow             = prob->hrow_;
451  CoinBigIndex *mcstrt          = prob->mcstrt_;
452  int *hincol           = prob->hincol_;
453  int *link             = prob->link_;
454
455  double *rcosts        = prob->rcosts_;
456
457  CoinBigIndex free_list                = prob->free_list_;
458
459  for (const action *f = &actions[nactions-1]; actions<=f; f--) {
460    int icol  = f->ithis;       // was fixed
461    int icol2 = f->ilast;       // was kept
462
463    dcost[icol] = dcost[icol2];
464    clo[icol] = f->thislo;
465    cup[icol] = f->thisup;
466    clo[icol2] = f->lastlo;
467    cup[icol2] = f->lastup;
468
469    create_col(icol, f->nincol, f->colrows, f->colels,
470               mcstrt, colels, hrow, link, &free_list);
471    hincol[icol] = hincol[icol2]; // right?
472
473    int nfinite = ((fabs(f->thislo) < PRESOLVE_INF) +
474                   (fabs(f->thisup) < PRESOLVE_INF) +
475                   (fabs(f->lastlo) < PRESOLVE_INF) +
476                   (fabs(f->lastup) < PRESOLVE_INF));
477
478    if (nfinite > 1) {
479       //      double l_j = f->thislo;
480      double u_j = f->thisup;
481      double l_k = f->lastlo;
482      double u_k = f->lastup;
483      double x_k_sol = sol[icol2];
484
485      prob->setColumnStatus(icol,prob->getColumnStatus(icol2));
486      if (x_k_sol <= l_k + u_j) {
487        sol[icol2] = l_k;
488        sol[icol] = x_k_sol - l_k;
489        prob->setColumnStatus(icol2,PrePostsolveMatrix::atLowerBound);
490      } else {
491        sol[icol2] = u_k;
492        sol[icol] = x_k_sol - u_k;
493        prob->setColumnStatus(icol2,PrePostsolveMatrix::atLowerBound);
494      }
495    } else if (nfinite == 1) {
496      double x_k_sol = sol[icol2];
497
498      if (fabs(f->thislo) < PRESOLVE_INF) {
499        prob->setColumnStatus(icol,PrePostsolveMatrix::atLowerBound);
500        sol[icol] = f->thislo;
501        sol[icol2] = x_k_sol - sol[icol];
502      } else if (fabs(f->thisup) < PRESOLVE_INF) {
503        prob->setColumnStatus(icol,PrePostsolveMatrix::atUpperBound);
504        sol[icol] = f->thisup;
505        sol[icol2] = x_k_sol - sol[icol];
506      } else if (fabs(f->lastlo) < PRESOLVE_INF) {
507        prob->setColumnStatus(icol,prob->getColumnStatus(icol2));
508        prob->setColumnStatus(icol2,PrePostsolveMatrix::atLowerBound);
509        sol[icol2] = f->lastlo;
510        sol[icol] = x_k_sol - sol[icol2];
511      } else {
512        // (fabs(f->lastup) < PRESOLVE_INF)
513        prob->setColumnStatus(icol,prob->getColumnStatus(icol2));
514        prob->setColumnStatus(icol2,PrePostsolveMatrix::atUpperBound);
515        sol[icol2] = f->lastup;
516        sol[icol] = x_k_sol - sol[icol2];
517      }
518    } else {
519      // both free!  superbasic time
520      sol[icol] = 0.0;  // doesn't matter
521      prob->setColumnStatus(icol2,PrePostsolveMatrix::isFree);
522    }
523
524    // row activity doesn't change
525    // dj of both variables is the same
526    rcosts[icol] = rcosts[icol2];
527
528#ifdef DEBUG_PRESOLVE
529    const double ztolzb = prob->ztolzb_;
530    if (! (clo[icol] - ztolzb <= sol[icol] && sol[icol] <= cup[icol] + ztolzb))
531             printf("BAD DUPCOL BOUNDS:  %g %g %g\n", clo[icol], sol[icol], cup[icol]);
532    if (! (clo[icol2] - ztolzb <= sol[icol2] && sol[icol2] <= cup[icol2] + ztolzb))
533             printf("BAD DUPCOL BOUNDS:  %g %g %g\n", clo[icol2], sol[icol2], cup[icol2]);
534#endif
535  }
536  prob->free_list_ = free_list;
537}
538
539
540
541
542
543// This is just ekkredc4, adapted into the new framework.
544const PresolveAction *duprow_action::presolve(PresolveMatrix *prob,
545                                               const PresolveAction *next)
546{
547  //  CoinBigIndex *mcstrt              = prob->mcstrt_;
548  //  int *hincol               = prob->hincol_;
549  int ncols             = prob->ncols_;
550
551  double *rowels        = prob->rowels_;
552  /*const*/ int *hcol   = prob->hcol_;
553  const CoinBigIndex *mrstrt    = prob->mrstrt_;
554  int *hinrow           = prob->hinrow_;
555  int nrows             = prob->nrows_;
556
557  double *rlo   = prob->rlo_;
558  double *rup   = prob->rup_;
559
560  int nuseless_rows = 0;
561
562  double * workcol = new double[ncols];
563  double * workrow = new double[nrows];
564  int * sort = new int[nrows];
565
566  init_random_vec(workcol, ncols);
567
568  int nlook = compute_sums(hinrow, mrstrt, hcol, rowels, workcol, workrow, sort, nrows);
569  CoinSort_2(workrow,workrow+nlook,sort);
570
571  double dval = workrow[0];
572  for (int jj = 1; jj < nlook; jj++) {
573    if (workrow[jj]==dval) {
574      int ithis=sort[jj];
575      int ilast=sort[jj-1];
576      CoinBigIndex krs = mrstrt[ithis];
577      CoinBigIndex kre = krs + hinrow[ithis];
578      if (hinrow[ithis] == hinrow[ilast]) {
579        int ishift = mrstrt[ilast] - krs;
580        CoinBigIndex k;
581        for (k=krs;k<kre;k++) {
582          if (hcol[k] != hcol[k+ishift] ||
583              rowels[k] != rowels[k+ishift]) {
584            break;
585          }
586        }
587        if (k == kre) {
588          /* now check rhs to see what is what */
589          double rlo1=rlo[ilast];
590          double rup1=rup[ilast];
591          double rlo2=rlo[ithis];
592          double rup2=rup[ithis];
593
594          int idelete=-1;
595          if (rlo1<=rlo2) {
596            if (rup2<=rup1) {
597              /* this is strictly tighter than last */
598              idelete=ilast;
599            } else if (fabs(rlo1-rlo2)<1.0e-12) {
600              /* last is strictly tighter than this */
601              idelete=ithis;
602              // swap so can carry on deleting
603              sort[jj-1]=ithis;
604              sort[jj]=ilast;
605            } else {
606              /* overlapping - could merge */
607#ifdef PRINT_DEBUG
608              printf("overlapping duplicate row %g %g, %g %g\n",
609                     rlo1,rup1,rlo2,rup2);
610#endif
611              if (rup1<rlo2) {
612                // infeasible
613                prob->status_|= 1;
614                // wrong message - correct if works
615                prob->messageHandler()->message(CLP_PRESOLVE_ROWINFEAS,
616                                                prob->messages())
617                                                  <<ithis
618                                                  <<rlo[ithis]
619                                                  <<rup[ithis]
620                                                  <<CoinMessageEol;
621                break;
622              }
623            }
624          } else {
625            // rlo2<=rlo1
626            if (rup1<=rup2) {
627              /* last is strictly tighter than this */
628              idelete=ithis;
629              // swap so can carry on deleting
630              sort[jj-1]=ithis;
631              sort[jj]=ilast;
632            } else {
633              /* overlapping - could merge */
634#ifdef PRINT_DEBUG
635              printf("overlapping duplicate row %g %g, %g %g\n",
636                     rlo1,rup1,rlo2,rup2);
637#endif
638              if (rup2<rlo1) {
639                // infeasible
640                prob->status_|= 1;
641                // wrong message - correct if works
642                prob->messageHandler()->message(CLP_PRESOLVE_ROWINFEAS,
643                                                prob->messages())
644                                                  <<ithis
645                                                  <<rlo[ithis]
646                                                  <<rup[ithis]
647                                                  <<CoinMessageEol;
648                break;
649              }
650            }
651          }
652          if (idelete>=0) 
653            sort[nuseless_rows++]=idelete;
654        }
655      }
656    }
657    dval=workrow[jj];
658  }
659
660  delete[]workrow;
661  delete[]workcol;
662
663
664  if (nuseless_rows) {
665#if     PRESOLVE_SUMMARY
666    printf("DUPLICATE ROWS:  %d\n", nuseless_rows);
667#endif
668    next = useless_constraint_action::presolve(prob,
669                                               sort, nuseless_rows,
670                                               next);
671  }
672  delete[]sort;
673
674  return (next);
675}
676
677void duprow_action::postsolve(PostsolveMatrix *prob) const
678{
679  printf("STILL NO POSTSOLVE FOR DUPROW!\n");
680  abort();
681}
Note: See TracBrowser for help on using the repository browser.