Ignore:
Timestamp:
Aug 30, 2012 11:43:19 AM (7 years ago)
Author:
forrest
Message:

minor changes to implement Aboca

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/ClpPresolve.cpp

    r1834 r1878  
    1313
    1414#include "CoinHelperFunctions.hpp"
     15#if CLP_HAS_ABC
     16#include "CoinAbcCommon.hpp"
     17#endif
    1518
    1619#include "CoinPackedMatrix.hpp"
     
    5558     nrows_(0),
    5659     nelems_(0),
     60#ifdef ABC_INHERIT
     61     numberPasses_(20),
     62#else
    5763     numberPasses_(5),
     64#endif
    5865     substitution_(3),
    5966#ifndef CLP_NO_STD
     
    274281#ifndef CLP_NO_STD
    275282     }
    276 #endif
     283#endif 
    277284     // put back duals
    278285     CoinMemcpyN(prob.rowduals_,        nrows_, originalModel_->dualRowSolution());
     
    444451}
    445452#endif
     453static int tightenDoubletons2(CoinPresolveMatrix * prob)
     454{
     455  // column-major representation
     456  const int ncols = prob->ncols_ ;
     457  const CoinBigIndex *const mcstrt = prob->mcstrt_ ;
     458  const int *const hincol = prob->hincol_ ;
     459  const int *const hrow = prob->hrow_ ;
     460  double * colels = prob->colels_ ;
     461  double * cost = prob->cost_ ;
     462
     463  // column type, bounds, solution, and status
     464  const unsigned char *const integerType = prob->integerType_ ;
     465  double * clo = prob->clo_ ;
     466  double * cup = prob->cup_ ;
     467  // row-major representation
     468  //const int nrows = prob->nrows_ ;
     469  const CoinBigIndex *const mrstrt = prob->mrstrt_ ;
     470  const int *const hinrow = prob->hinrow_ ;
     471  const int *const hcol = prob->hcol_ ;
     472  double * rowels = prob->rowels_ ;
     473
     474  // row bounds
     475  double *const rlo = prob->rlo_ ;
     476  double *const rup = prob->rup_ ;
     477
     478  // tolerances
     479  //const double ekkinf2 = PRESOLVE_SMALL_INF ;
     480  //const double ekkinf = ekkinf2*1.0e8 ;
     481  //const double ztolcbarj = prob->ztoldj_ ;
     482  //const CoinRelFltEq relEq(prob->ztolzb_) ;
     483  int numberChanged=0;
     484  double bound[2];
     485  double alpha[2]={0.0,0.0};
     486  double offset=0.0;
     487
     488  for (int icol=0;icol<ncols;icol++) {
     489    if (hincol[icol]==2) {
     490      CoinBigIndex start=mcstrt[icol];
     491      int row0 = hrow[start];
     492      if (hinrow[row0]!=2)
     493        continue;
     494      int row1 = hrow[start+1];
     495      if (hinrow[row1]!=2)
     496        continue;
     497      double element0 = colels[start];
     498      double rowUpper0=rup[row0];
     499      bool swapSigns0=false;
     500      if (rlo[row0]>-1.0e30) {
     501        if (rup[row0]>1.0e30) {
     502          swapSigns0=true;
     503          rowUpper0=-rlo[row0];
     504          element0=-element0;
     505        } else {
     506          // range or equality
     507          continue;
     508        }
     509      } else if (rup[row0]>1.0e30) {
     510        // free
     511        continue;
     512      }
     513#if 0
     514      // skip here for speed
     515      // skip if no cost (should be able to get rid of)
     516      if (!cost[icol]) {
     517        printf("should be able to get rid of %d with no cost\n",icol);
     518        continue;
     519      }
     520      // skip if negative cost for now
     521      if (cost[icol]<0.0) {
     522        printf("code for negative cost\n");
     523        continue;
     524      }
     525#endif
     526      double element1 = colels[start+1];
     527      double rowUpper1=rup[row1];
     528      bool swapSigns1=false;
     529      if (rlo[row1]>-1.0e30) {
     530        if (rup[row1]>1.0e30) {
     531          swapSigns1=true;
     532          rowUpper1=-rlo[row1];
     533          element1=-element1;
     534        } else {
     535          // range or equality
     536          continue;
     537        }
     538      } else if (rup[row1]>1.0e30) {
     539        // free
     540        continue;
     541      }
     542      double lowerX=clo[icol];
     543      double upperX=cup[icol];
     544      int otherCol=-1;
     545      CoinBigIndex startRow=mrstrt[row0];
     546      for (CoinBigIndex j=startRow;j<startRow+2;j++) {
     547        int jcol=hcol[j];
     548        if (jcol!=icol) {
     549          alpha[0]=swapSigns0 ? -rowels[j] :rowels[j];
     550          otherCol=jcol;
     551        }
     552      }
     553      startRow=mrstrt[row1];
     554      bool possible=true;
     555      for (CoinBigIndex j=startRow;j<startRow+2;j++) {
     556        int jcol=hcol[j];
     557        if (jcol!=icol) {
     558          if (jcol==otherCol) {
     559            alpha[1]=swapSigns1 ? -rowels[j] :rowels[j];
     560          } else {
     561            possible=false;
     562          }
     563        }
     564      }
     565      if (possible) {
     566        // skip if no cost (should be able to get rid of)
     567        if (!cost[icol]) {
     568          printf("should be able to get rid of %d with no cost\n",icol);
     569          continue;
     570        }
     571        // skip if negative cost for now
     572        if (cost[icol]<0.0) {
     573          printf("code for negative cost\n");
     574          continue;
     575        }
     576        bound[0]=clo[otherCol];
     577        bound[1]=cup[otherCol];
     578        double lowestLowest=COIN_DBL_MAX;
     579        double highestLowest=-COIN_DBL_MAX;
     580        double lowestHighest=COIN_DBL_MAX;
     581        double highestHighest=-COIN_DBL_MAX;
     582        int binding0=0;
     583        int binding1=0;
     584        for (int k=0;k<2;k++) {
     585          bool infLow0=false;
     586          bool infLow1=false;
     587          double sum0=0.0;
     588          double sum1=0.0;
     589          double value=bound[k];
     590          if (fabs(value)<1.0e30) {
     591            sum0+=alpha[0]*value;
     592            sum1+=alpha[1]*value;
     593          } else {
     594            if (alpha[0]>0.0) {
     595              if (value<0.0)
     596                infLow0 =true;
     597            } else if (alpha[0]<0.0) {
     598              if (value>0.0)
     599                infLow0 =true;
     600            }
     601            if (alpha[1]>0.0) {
     602              if (value<0.0)
     603                infLow1 =true;
     604            } else if (alpha[1]<0.0) {
     605              if (value>0.0)
     606                infLow1 =true;
     607            }
     608          }
     609          /* Got sums
     610           */
     611          double thisLowest0=-COIN_DBL_MAX;
     612          double thisHighest0=COIN_DBL_MAX;
     613          if (element0>0.0) {
     614            // upper bound unless inf&2 !=0
     615            if (!infLow0)
     616              thisHighest0 = (rowUpper0-sum0)/element0;
     617          } else {
     618            // lower bound unless inf&2 !=0
     619            if (!infLow0)
     620              thisLowest0 = (rowUpper0-sum0)/element0;
     621          }
     622          double thisLowest1=-COIN_DBL_MAX;
     623          double thisHighest1=COIN_DBL_MAX;
     624          if (element1>0.0) {
     625            // upper bound unless inf&2 !=0
     626            if (!infLow1)
     627              thisHighest1 = (rowUpper1-sum1)/element1;
     628          } else {
     629            // lower bound unless inf&2 !=0
     630            if (!infLow1)
     631              thisLowest1 = (rowUpper1-sum1)/element1;
     632          }
     633          if (thisLowest0>thisLowest1+1.0e-12) {
     634            if (thisLowest0>lowerX+1.0e-12)
     635              binding0|= 1<<k;
     636          } else if (thisLowest1>thisLowest0+1.0e-12) {
     637            if (thisLowest1>lowerX+1.0e-12)
     638              binding1|= 1<<k;
     639            thisLowest0=thisLowest1;
     640          }
     641          if (thisHighest0<thisHighest1-1.0e-12) {
     642            if (thisHighest0<upperX-1.0e-12)
     643              binding0|= 1<<k;
     644          } else if (thisHighest1<thisHighest0-1.0e-12) {
     645            if (thisHighest1<upperX-1.0e-12)
     646              binding1|= 1<<k;
     647            thisHighest0=thisHighest1;
     648          }
     649          lowestLowest=CoinMin(lowestLowest,thisLowest0);
     650          highestHighest=CoinMax(highestHighest,thisHighest0);
     651          lowestHighest=CoinMin(lowestHighest,thisHighest0);
     652          highestLowest=CoinMax(highestLowest,thisLowest0);
     653        }
     654        // see if any good
     655        //#define PRINT_VALUES
     656        if (!binding0||!binding1) {
     657          printf("Row redundant for column %d\n",icol);
     658        } else {
     659#ifdef PRINT_VALUES
     660          printf("Column %d bounds %g,%g lowest %g,%g highest %g,%g\n",
     661                 icol,lowerX,upperX,lowestLowest,lowestHighest,
     662                 highestLowest,highestHighest);
     663#endif
     664          // if integer adjust
     665          if (integerType[icol]) {
     666            lowestLowest=ceil(lowestLowest-1.0e-5);
     667            highestLowest=ceil(highestLowest-1.0e-5);
     668            lowestHighest=floor(lowestHighest+1.0e-5);
     669            highestHighest=floor(highestHighest+1.0e-5);
     670          }
     671          // if costed may be able to adjust
     672          if (cost[icol]>=0.0) {
     673            if (highestLowest<upperX&&highestLowest>=lowerX&&highestHighest<1.0e30) {
     674              highestHighest=CoinMin(highestHighest,highestLowest);
     675            }
     676          }
     677          if (cost[icol]<=0.0) {
     678            if (lowestHighest>lowerX&&lowestHighest<=upperX&&lowestHighest>-1.0e30) {
     679              lowestLowest=CoinMax(lowestLowest,lowestHighest);
     680            }
     681          }
     682#if 1
     683          if (lowestLowest>lowerX+1.0e-8) {
     684#ifdef PRINT_VALUES
     685            printf("Can increase lower bound on %d from %g to %g\n",
     686                   icol,lowerX,lowestLowest);
     687#endif
     688            lowerX=lowestLowest;
     689          }
     690          if (highestHighest<upperX-1.0e-8) {
     691#ifdef PRINT_VALUES
     692            printf("Can decrease upper bound on %d from %g to %g\n",
     693                   icol,upperX,highestHighest);
     694#endif
     695            upperX=highestHighest;
     696           
     697          }
     698#endif
     699          // see if we can move costs
     700          double xValue;
     701          double yValue0;
     702          double yValue1;
     703          double newLower=COIN_DBL_MAX;
     704          double newUpper=-COIN_DBL_MAX;
     705          double ranges0[2];
     706          double ranges1[2];
     707          double costEqual;
     708          double slope[2];
     709          assert (binding0+binding1==3);
     710          // get where equal
     711          xValue=(rowUpper0*element1-rowUpper1*element0)/(alpha[0]*element1-alpha[1]*element0);
     712          yValue0=(rowUpper0-xValue*alpha[0])/element0;
     713          yValue1=(rowUpper1-xValue*alpha[1])/element1;
     714          newLower=CoinMin(newLower,CoinMax(yValue0,yValue1));
     715          newUpper=CoinMax(newUpper,CoinMax(yValue0,yValue1));
     716          double xValueEqual=xValue;
     717          double yValueEqual=yValue0;
     718          costEqual = xValue*cost[otherCol]+yValueEqual*cost[icol];
     719          if (binding0==1) {
     720            ranges0[0]=bound[0];
     721            ranges0[1]=yValue0;
     722            ranges1[0]=yValue0;
     723            ranges1[1]=bound[1];
     724            // take x 1.0 down
     725            double x=xValue-1.0;
     726            double y=(rowUpper0-x*alpha[0])/element0;
     727            double costTotal = x*cost[otherCol]+y*cost[icol];
     728            slope[0] = costEqual-costTotal;
     729            // take x 1.0 up
     730            x=xValue+1.0;
     731            y=(rowUpper1-x*alpha[1])/element0;
     732            costTotal = x*cost[otherCol]+y*cost[icol];
     733            slope[1] = costTotal-costEqual;
     734          } else {
     735            ranges1[0]=bound[0];
     736            ranges1[1]=yValue0;
     737            ranges0[0]=yValue0;
     738            ranges0[1]=bound[1];
     739            // take x 1.0 down
     740            double x=xValue-1.0;
     741            double y=(rowUpper1-x*alpha[1])/element0;
     742            double costTotal = x*cost[otherCol]+y*cost[icol];
     743            slope[1] = costEqual-costTotal;
     744            // take x 1.0 up
     745            x=xValue+1.0;
     746            y=(rowUpper0-x*alpha[0])/element0;
     747            costTotal = x*cost[otherCol]+y*cost[icol];
     748            slope[0] = costTotal-costEqual;
     749          }
     750#ifdef PRINT_VALUES
     751          printf("equal value of %d is %g, value of %d is max(%g,%g) - %g\n",
     752                 otherCol,xValue,icol,yValue0,yValue1,CoinMax(yValue0,yValue1));
     753          printf("Cost at equality %g for constraint 0 ranges %g -> %g slope %g for constraint 1 ranges %g -> %g slope %g\n",
     754                 costEqual,ranges0[0],ranges0[1],slope[0],ranges1[0],ranges1[1],slope[1]);
     755#endif
     756          xValue=bound[0];
     757          yValue0=(rowUpper0-xValue*alpha[0])/element0;
     758          yValue1=(rowUpper1-xValue*alpha[1])/element1;
     759#ifdef PRINT_VALUES
     760          printf("value of %d is %g, value of %d is max(%g,%g) - %g\n",
     761                 otherCol,xValue,icol,yValue0,yValue1,CoinMax(yValue0,yValue1));
     762#endif
     763          newLower=CoinMin(newLower,CoinMax(yValue0,yValue1));
     764          // cost>0 so will be at lower
     765          //double yValueAtBound0=newLower;
     766          newUpper=CoinMax(newUpper,CoinMax(yValue0,yValue1));
     767          xValue=bound[1];
     768          yValue0=(rowUpper0-xValue*alpha[0])/element0;
     769          yValue1=(rowUpper1-xValue*alpha[1])/element1;
     770#ifdef PRINT_VALUES
     771          printf("value of %d is %g, value of %d is max(%g,%g) - %g\n",
     772                 otherCol,xValue,icol,yValue0,yValue1,CoinMax(yValue0,yValue1));
     773#endif
     774          newLower=CoinMin(newLower,CoinMax(yValue0,yValue1));
     775          // cost>0 so will be at lower
     776          //double yValueAtBound1=newLower;
     777          newUpper=CoinMax(newUpper,CoinMax(yValue0,yValue1));
     778          lowerX=CoinMax(lowerX,newLower-1.0e-12*fabs(newLower));
     779          upperX=CoinMin(upperX,newUpper+1.0e-12*fabs(newUpper));
     780          // Now make duplicate row
     781          // keep row 0 so need to adjust costs so same
     782#ifdef PRINT_VALUES
     783          printf("Costs for x %g,%g,%g are %g,%g,%g\n",
     784                 xValueEqual-1.0,xValueEqual,xValueEqual+1.0,
     785                 costEqual-slope[0],costEqual,costEqual+slope[1]);
     786#endif
     787          double costOther=cost[otherCol]+slope[1];
     788          double costThis=cost[icol]+slope[1]*(element0/alpha[0]);
     789          xValue=xValueEqual;
     790          yValue0=CoinMax((rowUpper0-xValue*alpha[0])/element0,lowerX);
     791          double thisOffset=costEqual-(costOther*xValue+costThis*yValue0);
     792          offset += thisOffset;
     793#ifdef PRINT_VALUES
     794          printf("new cost at equal %g\n",costOther*xValue+costThis*yValue0+thisOffset);
     795#endif
     796          xValue=xValueEqual-1.0;
     797          yValue0=CoinMax((rowUpper0-xValue*alpha[0])/element0,lowerX);
     798#ifdef PRINT_VALUES
     799          printf("new cost at -1 %g\n",costOther*xValue+costThis*yValue0+thisOffset);
     800#endif
     801          assert(fabs((costOther*xValue+costThis*yValue0+thisOffset)-(costEqual-slope[0]))<1.0e-5);
     802          xValue=xValueEqual+1.0;
     803          yValue0=CoinMax((rowUpper0-xValue*alpha[0])/element0,lowerX);
     804#ifdef PRINT_VALUES
     805          printf("new cost at +1 %g\n",costOther*xValue+costThis*yValue0+thisOffset);
     806#endif
     807          assert(fabs((costOther*xValue+costThis*yValue0+thisOffset)-(costEqual+slope[1]))<1.0e-5);
     808          numberChanged++;
     809          //      continue;
     810          cost[otherCol] = costOther;
     811          cost[icol] = costThis;
     812          clo[icol]=lowerX;
     813          cup[icol]=upperX;
     814          int startCol[2];
     815          int endCol[2];
     816          startCol[0]=mcstrt[icol];
     817          endCol[0]=startCol[0]+2;
     818          startCol[1]=mcstrt[otherCol];
     819          endCol[1]=startCol[1]+hincol[otherCol];
     820          double values[2]={0.0,0.0};
     821          for (int k=0;k<2;k++) {
     822            for (CoinBigIndex i=startCol[k];i<endCol[k];i++) {
     823              if (hrow[i]==row0)
     824                values[k]=colels[i];
     825            }
     826            for (CoinBigIndex i=startCol[k];i<endCol[k];i++) {
     827              if (hrow[i]==row1)
     828                colels[i]=values[k];
     829            }
     830          }
     831          for (CoinBigIndex i=mrstrt[row1];i<mrstrt[row1]+2;i++) {
     832            if (hcol[i]==icol)
     833              rowels[i]=values[0];
     834            else
     835              rowels[i]=values[1];
     836          }
     837        }
     838      }
     839    }
     840  }
     841#if ABC_NORMAL_DEBUG>0
     842  if (offset)
     843    printf("Cost offset %g\n",offset);
     844#endif
     845  return numberChanged;
     846}
    446847//#define COIN_PRESOLVE_BUG
    447848#ifdef COIN_PRESOLVE_BUG
     
    557958          }
    558959
     960            // transfer costs (may want to do it in OsiPresolve)
     961            // need a transfer back at end of postsolve transferCosts(prob);
    559962
    560963          int iLoop;
     
    569972               paction_ = dupcol_action::presolve(prob, paction_);
    570973          }
     974#ifdef ABC_INHERIT
     975          if (doTwoxTwo()) {
     976            possibleSkip;
     977            paction_ = twoxtwo_action::presolve(prob, paction_);
     978          }
     979#endif
    571980          if (duprow) {
    572981            possibleSkip;
    573                paction_ = duprow_action::presolve(prob, paction_);
     982            int nTightened=tightenDoubletons2(prob);
     983            if (nTightened)
     984              printf("%d doubletons tightened\n",nTightened);
     985            paction_ = duprow_action::presolve(prob, paction_);
    574986          }
    575987          if (doGubrow()) {
     
    583995          int lastDropped = 0;
    584996          prob->pass_ = 0;
     997#ifdef ABC_INHERIT
     998          int numberRowsStart=nrows_-prob->countEmptyRows();
     999          int numberColumnsStart=ncols_-prob->countEmptyCols();
     1000          int numberRowsLeft=numberRowsStart;
     1001          int numberColumnsLeft=numberColumnsStart;
     1002          bool lastPassWasGood=true;
     1003#if ABC_NORMAL_DEBUG
     1004          printf("Original rows,columns %d,%d starting first pass with %d,%d\n",
     1005                 nrows_,ncols_,numberRowsLeft,numberColumnsLeft);
     1006#endif
     1007#endif
    5851008          if (numberPasses_<=5)
    5861009              prob->presolveOptions_ |= 0x10000; // say more lightweight
     
    8671290               if (paction_ == paction0 || stopLoop)
    8681291                    break;
    869 
     1292#ifdef ABC_INHERIT
     1293               // see whether to stop anyway
     1294               int numberRowsNow=nrows_-prob->countEmptyRows();
     1295               int numberColumnsNow=ncols_-prob->countEmptyCols();
     1296#if ABC_NORMAL_DEBUG
     1297               printf("Original rows,columns %d,%d - last %d,%d end of pass %d has %d,%d\n",
     1298                      nrows_,ncols_,numberRowsLeft,numberColumnsLeft,iLoop+1,numberRowsNow,
     1299                      numberColumnsNow);
     1300#endif
     1301               int rowsDeleted=numberRowsLeft-numberRowsNow;
     1302               int columnsDeleted=numberColumnsLeft-numberColumnsNow;
     1303               if (iLoop>15) {
     1304                 if (rowsDeleted*100<numberRowsStart&&
     1305                     columnsDeleted*100<numberColumnsStart)
     1306                   break;
     1307                 lastPassWasGood=true;
     1308               } else if (rowsDeleted*100<numberRowsStart&&rowsDeleted<500&&
     1309                          columnsDeleted*100<numberColumnsStart&&columnsDeleted<500) {
     1310                 if (!lastPassWasGood)
     1311                   break;
     1312                 else
     1313                   lastPassWasGood=false;
     1314               } else {
     1315                 lastPassWasGood=true;
     1316               }
     1317               numberRowsLeft=numberRowsNow;
     1318               numberColumnsLeft=numberColumnsNow;
     1319#endif
    8701320          }
    8711321     }
     
    9431393               }
    9441394          }
     1395     }
     1396     if (prob.maxmin_<0) {
     1397       //for (int i=0;i<presolvedModel_->numberRows();i++)
     1398       //prob.rowduals_[i]=-prob.rowduals_[i];
     1399       for (int i=0;i<ncols_;i++) {
     1400         prob.cost_[i]=-prob.cost_[i];
     1401         //prob.rcosts_[i]=-prob.rcosts_[i];
     1402       }
     1403       prob.maxmin_=1.0;
    9451404     }
    9461405     const CoinPresolveAction *paction = paction_;
     
    12451704     mcstrt_[0] = 0;
    12461705     ClpDisjointCopyN(m->getVectorLengths(), ncols_,  hincol_);
     1706     if (si->getObjSense() < 0.0) {
     1707       for (int i=0;i<ncols_;i++)
     1708         cost_[i]=-cost_[i];
     1709       maxmin_=1.0;
     1710     }
    12471711     for (icol = 0; icol < ncols_; icol++) {
    12481712          CoinBigIndex j;
     
    14161880#endif
    14171881{
     1882     if (si->getObjSense() < 0.0) {
     1883       for (int i=0;i<ncols_;i++)
     1884         cost_[i]=-cost_[i];
     1885       dobias_=-dobias_;
     1886     }
    14181887     si->loadProblem(ncols_, nrows_, mcstrt_, hrow_, colels_, hincol_,
    14191888                     clo_, cup_, cost_, rlo_, rup_);
     
    14361905#endif
    14371906     si->setDblParam(ClpObjOffset, originalOffset_ - dobias_);
     1907     if (si->getObjSense() < 0.0) {
     1908       // put back
     1909       for (int i=0;i<ncols_;i++)
     1910         cost_[i]=-cost_[i];
     1911       dobias_=-dobias_;
     1912       maxmin_=-1.0;
     1913     }
    14381914
    14391915}
     
    19042380               if (rowObjective_) {
    19052381                    int iRow;
     2382#ifndef NDEBUG
    19062383                    int k = -1;
     2384#endif
    19072385                    int nObj = 0;
    19082386                    for (iRow = 0; iRow < nrowsNow; iRow++) {
    19092387                         int kRow = originalRow_[iRow];
     2388#ifndef NDEBUG
    19102389                         assert (kRow > k);
    19112390                         k = kRow;
     2391#endif
    19122392                         rowObjective_[iRow] = rowObjective_[kRow];
    19132393                         if (rowObjective_[iRow])
Note: See TracChangeset for help on using the changeset viewer.