source: branches/devel-1/PresolveDupcol.cpp @ 29

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

Presolve (no changes to Makefile)

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