Ignore:
Timestamp:
Dec 10, 2004 2:53:37 PM (15 years ago)
Author:
forrest
Message:

Slightly faster

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/CbcFathomDynamicProgramming.cpp

    r36 r38  
    1515#include "CoinHelperFunctions.hpp"
    1616#include "CoinPackedMatrix.hpp"
     17#include "CoinSort.hpp"
    1718// Default Constructor
    1819CbcFathomDynamicProgramming::CbcFathomDynamicProgramming()
     
    2526   lookup_(NULL),
    2627   numberActive_(0),
     28   maximumSizeAllowed_(1000000),
    2729   startBit_(NULL),
    2830   numberBits_(NULL),
    29    rhs_(NULL)
     31   rhs_(NULL),
     32   numberNonOne_(0)
    3033{
    3134
     
    4043   lookup_(NULL),
    4144   numberActive_(0),
     45   maximumSizeAllowed_(1000000),
    4246   startBit_(NULL),
    4347   numberBits_(NULL),
    44    rhs_(NULL)
     48   rhs_(NULL),
     49   numberNonOne_(0)
    4550{
    4651  type_=gutsOfCheckPossible();
     
    6772  id_ = NULL;
    6873  lookup_ = NULL;
    69   numberActive_ = 0;
    7074  startBit_ = NULL;
    7175  numberBits_ = NULL;
     
    9094  lookup_(NULL),
    9195  numberActive_(rhs.numberActive_),
     96  maximumSizeAllowed_(rhs.maximumSizeAllowed_),
    9297  startBit_(NULL),
    9398  numberBits_(NULL),
    94   rhs_(NULL)
     99  rhs_(NULL),
     100  numberNonOne_(rhs.numberNonOne_)
    95101{
    96102  if (size_) {
     
    273279  }
    274280  delete [] rhs;
     281  if (allowableSize&&size_>allowableSize) {
     282    printf("Too large - need %d entries x 16 bytes\n",size_);
     283    return -1; // too big
     284  }
    275285  if (n01==numberColumns&&!nbadcoeff)
    276286    return 0; // easiest
     
    290300{
    291301  int returnCode=0;
    292   int type=gutsOfCheckPossible(1000000);
     302  int type=gutsOfCheckPossible(maximumSizeAllowed_);
    293303  if (type>=0) {
    294304    bool gotSolution=false;
     
    309319   
    310320    int numberColumns = solver->getNumCols();
     321    int * indices = new int [numberActive_];
     322    double offset;
     323    solver->getDblParam(OsiObjOffset,offset);
     324    double fixedObj=-offset;
     325    int i;
    311326    // may be possible
    312     // for test - just if all 1
    313327    if (type==0) {
    314       int i;
    315 
    316       int * indices = new int [numberActive_];
    317       double offset;
    318       solver->getDblParam(OsiObjOffset,offset);
    319       double fixedObj=-offset;
     328      // rhs 1 and coefficients 1
    320329      for (i=0;i<numberColumns;i++) {
    321330        int j;
     
    343352        }
    344353      }
    345       int needed=0;
     354      returnCode=1;
     355    } else {
     356      // more complex
     357      int * coefficients = new int[numberActive_];
     358      // If not too many general rhs then we can be more efficient
     359      numberNonOne_=0;
     360      for (i=0;i<numberActive_;i++) {
     361        if (rhs_[i]!=1)
     362          numberNonOne_++;
     363      }
     364      if (numberNonOne_*2<numberActive_) {
     365        // put rhs >1 every second
     366        int * permute = new int[numberActive_];
     367        int * temp = new int[numberActive_];
     368        // try different ways
     369        int k=0;
     370        for (i=0;i<numberRows;i++) {
     371          int newRow = lookup_[i];
     372          if (newRow>=0&&rhs_[newRow]>1) {
     373            permute[newRow]=k;
     374            k +=2;
     375          }
     376        }
     377        // adjust so k points to last
     378        k -= 2;
     379        // and now rest
     380        int k1=1;
     381        for (i=0;i<numberRows;i++) {
     382          int newRow = lookup_[i];
     383          if (newRow>=0&&rhs_[newRow]==1) {
     384            permute[newRow]=k1;
     385            k1++;
     386            if (k1<=k)
     387              k1++;
     388          }
     389        }
     390        for (i=0;i<numberActive_;i++) {
     391          int put = permute[i];
     392          temp[put]=rhs_[i];
     393        }
     394        memcpy(rhs_,temp,numberActive_*sizeof(int));
     395        for (i=0;i<numberActive_;i++) {
     396          int put = permute[i];
     397          temp[put]=numberBits_[i];
     398        }
     399        memcpy(numberBits_,temp,numberActive_*sizeof(int));
     400        k=0;
     401        for (i=0;i<numberActive_;i++) {
     402          startBit_[i]=k;
     403          k+= numberBits_[i];
     404        }
     405        for (i=0;i<numberRows;i++) {
     406          int newRow = lookup_[i];
     407          if (newRow>=0)
     408            lookup_[i]=permute[newRow];
     409        }
     410        delete [] permute;
     411        delete [] temp;
     412        // mark new method
     413        type=2;
     414      }
     415      for (i=0;i<numberColumns;i++) {
     416        int j;
     417        double lowerValue = lower[i];
     418        assert (lowerValue==floor(lowerValue));
     419        double cost = direction * objective[i];
     420        fixedObj += lowerValue*cost;
     421        if (lowerValue==upper[i])
     422          continue;
     423        int n=0;
     424        double gap=upper[i]-lowerValue;
     425        for (j=columnStart[i];
     426             j<columnStart[i]+columnLength[i];j++) {
     427          int iRow=row[j];
     428          double value = element[j];
     429          int iValue = (int) value;
     430          int newRow = lookup_[iRow];
     431          if (newRow<0||iValue>rhs_[newRow]) {
     432            n=0;
     433            break; //can't use
     434          } else {
     435            coefficients[n]=iValue;
     436            indices[n++]=newRow;
     437            if (gap*iValue>rhs_[newRow]) {
     438              gap = rhs_[newRow]/iValue;
     439            }
     440          }
     441        }
     442        int nTotal = (int) gap;
     443        if (n) {
     444          if (type==1) {
     445            for (int k=1;k<=nTotal;k++) {
     446              addOneColumn1(i,n,indices,coefficients,cost);
     447            }
     448          } else {
     449            // may be able to take sort out with new method
     450            CoinSort_2(indices,indices+n,coefficients);
     451            for (int k=1;k<=nTotal;k++) {
     452              addOneColumn1A(i,n,indices,coefficients,cost);
     453            }
     454          }
     455        }
     456      }
     457      delete [] coefficients;
     458      returnCode=1;
     459    }
     460    int needed=0;
     461    double bestValue=COIN_DBL_MAX;
     462    int iBest=-1;
     463    if (type==0) {
    346464      int numberActive=0;
    347465      for (i=0;i<numberRows;i++) {
     
    350468          if (rowLower[i]==rowUpper[i]) {
    351469            needed += 1<<numberActive;
    352           }
    353           numberActive++;
    354         }
    355       }
    356       double bestValue=COIN_DBL_MAX;
    357       int iBest=-1;
     470            numberActive++;
     471          }
     472        }
     473      }
    358474      for (i=0;i<size_;i++) {
    359475        if ((i&needed)==needed) {
     
    365481        }
    366482      }
    367       returnCode=1;
    368       if (bestValue<COIN_DBL_MAX) {
    369         bestValue += fixedObj;
    370         printf("Can get solution of %g\n",bestValue);
    371         if (bestValue<model_->getMinimizationObjValue()) {
    372           // set up solution
    373           betterSolution = new double[numberColumns];
    374           memcpy(betterSolution,lower,numberColumns*sizeof(double));
    375           while (iBest>0) {
    376             int iColumn = id_[iBest];
    377             assert (iColumn>=0);
    378             betterSolution[iColumn]++;
    379             assert (betterSolution[iColumn]<=upper[iColumn]);
    380             iBest = back_[iBest];
    381           }
    382         } else {
    383         }
    384       }
    385       delete [] indices;
    386     }
     483    } else {
     484      int * lower = new int[numberActive_];
     485      for (i=0;i<numberRows;i++) {
     486        int newRow = lookup_[i];
     487        if (newRow>=0) {
     488          int gap=(int) (rowUpper[i]-CoinMax(0.0,rowLower[i]));
     489          lower[newRow]=rhs_[newRow]-gap;
     490          int numberBits = numberBits_[newRow];
     491          int startBit = startBit_[newRow];
     492          if (numberBits==1&&!gap) {
     493            needed |= 1<<startBit;
     494          }
     495        }
     496      }
     497      for (i=0;i<size_;i++) {
     498        if ((i&needed)==needed) {
     499        // this one may do
     500          bool good = true;
     501          for (int kk=0;kk<numberActive_;kk++) {
     502            int numberBits = numberBits_[kk];
     503            int startBit = startBit_[kk];
     504            int size = 1<<numberBits;
     505            int start = 1<<startBit;
     506            int mask = start*(size-1);
     507            int level=(i&mask)>>startBit;
     508            if (level<lower[kk]) {
     509              good=false;
     510              break;
     511            }
     512          }
     513          if (good&&cost_[i]<bestValue) {
     514            bestValue=cost_[i];
     515            iBest=i;
     516          }
     517        }
     518      }
     519      delete [] lower;
     520    }
     521    if (bestValue<COIN_DBL_MAX) {
     522      bestValue += fixedObj;
     523      printf("Can get solution of %g\n",bestValue);
     524      if (bestValue<model_->getMinimizationObjValue()) {
     525        // set up solution
     526        betterSolution = new double[numberColumns];
     527        memcpy(betterSolution,lower,numberColumns*sizeof(double));
     528        while (iBest>0) {
     529          int iColumn = id_[iBest];
     530          assert (iColumn>=0);
     531          betterSolution[iColumn]++;
     532          assert (betterSolution[iColumn]<=upper[iColumn]);
     533          iBest = back_[iBest];
     534        }
     535      }
     536    }
     537    delete [] indices;
    387538    gutsOfDelete();
    388539    if (gotSolution) {
     
    436587    mask |= 1<<iRow;
    437588  }
    438   i=0;
     589  i=size_-1-mask;
    439590  bool touched = false;
    440   while (i<size_) {
     591  while (i>=0) {
    441592    int kMask = i&mask;
    442593    if (kMask==0) {
     
    453604        }
    454605      }
    455       i++;
     606      i--;
    456607    } else {
    457608      // we can skip some
    458       int k=i;
    459       int iBit=0;
    460       k &= ~1;
    461       while ((k&kMask)!=0) {
    462         iBit++;
    463         k &= ~(1<<iBit);
    464       }
    465       // onto next
    466       k += 1<<(iBit+1);
     609      int k=(i&~mask)-1;
    467610#ifdef CBC_DEBUG
    468       for (int j=i+1;j<k;j++) {
     611      for (int j=i-1;j>k;j--) {
    469612        int jMask = j&mask;
    470613        assert (jMask!=0);
     
    476619  return touched;
    477620}
     621/* Adds one attempt of one column of type 1,
     622   returns true if was used in making any changes.
     623   At present the user has to call it once for each possible value
     624*/
     625bool
     626CbcFathomDynamicProgramming::addOneColumn1(int id,int numberElements, const int * rows,
     627                                           const int * coefficients, double cost)
     628{
     629  /* build up masks.
     630     a) mask for 1 rhs
     631     b) mask for addition
     632     c) mask so adding will overflow
     633     d) individual masks
     634  */
     635  int mask1=0;
     636  int maskAdd=0;
     637  int mask2=0;
     638  int i;
     639  int n2=0;
     640  int mask[40];
     641  int nextbit[40];
     642  int adjust[40];
     643  assert (numberElements<=40);
     644  for (i=0;i<numberElements;i++) {
     645    int iRow=rows[i];
     646    int numberBits = numberBits_[iRow];
     647    int startBit = startBit_[iRow];
     648    if (numberBits==1) {
     649      mask1 |= 1<<startBit;
     650      maskAdd |= 1<<startBit;
     651      mask2 |= 1<<startBit;
     652    } else {
     653      int value=coefficients[i];
     654      int size = 1<<numberBits;
     655      int start = 1<<startBit;
     656      assert (value<size);
     657      maskAdd |= start*value;
     658      int gap = size-rhs_[iRow]-1;
     659      assert (gap>=0);
     660      int hi2=rhs_[iRow]-value;
     661      if (hi2<size-1)
     662        hi2++;
     663      adjust[n2] = start*hi2;
     664      mask2 += start*gap;
     665      nextbit[n2]=startBit+numberBits;
     666      mask[n2++] = start*(size-1);
     667    }
     668  }
     669  i=size_-1-maskAdd;
     670  bool touched = false;
     671  while (i>=0) {
     672    int kMask = i&mask1;
     673    if (kMask==0) {
     674      bool good=true;
     675      for (int kk=n2-1;kk>=0;kk--) {
     676        int iMask = mask[kk];
     677        int jMask = iMask&mask2;
     678        int kkMask = iMask&i;
     679        kkMask += jMask;
     680        if (kkMask>iMask) {
     681          // we can skip some
     682          int k=(i&~iMask);
     683          k |= adjust[kk]; 
     684#ifdef CBC_DEBUG
     685          for (int j=i-1;j>k;j--) {
     686            int jMask = j&mask1;
     687            if (jMask==0) {
     688              bool good=true;
     689              for (int kk=n2-1;kk>=0;kk--) {
     690                int iMask = mask[kk];
     691                int jMask = iMask&mask2;
     692                int kkMask = iMask&i;
     693                kkMask += jMask;
     694                if (kkMask>iMask) {
     695                  good=false;
     696                  break;
     697                }
     698              }
     699              assert (!good);
     700            }
     701          }
     702#endif
     703          i=k;
     704          good=false;
     705          break;
     706        }
     707      }
     708      if (good) {
     709        double thisCost = cost_[i];
     710        if (thisCost!=COIN_DBL_MAX) {
     711          // possible
     712          double newCost=thisCost+cost;
     713          int next = i + maskAdd;
     714          if (cost_[next]>newCost) {
     715            cost_[next]=newCost;
     716            back_[next]=i;
     717            id_[next]=id;
     718            touched=true;
     719          }
     720        }
     721      }
     722      i--;
     723    } else {
     724      // we can skip some
     725      // we can skip some
     726      int k=(i&~mask1)-1;
     727#ifdef CBC_DEBUG
     728      for (int j=i-1;j>k;j--) {
     729        int jMask = j&mask;
     730        assert (jMask!=0);
     731      }
     732#endif
     733      i=k;
     734    }
     735  }
     736  return touched;
     737}
     738/* Adds one attempt of one column of type 1,
     739   returns true if was used in making any changes.
     740   At present the user has to call it once for each possible value
     741   This version is when there are enough 1 rhs to do faster
     742*/
     743bool
     744CbcFathomDynamicProgramming::addOneColumn1A(int id,int numberElements, const int * rows,
     745                                           const int * coefficients, double cost)
     746{
     747  /* build up masks.
     748     a) mask for 1 rhs
     749     b) mask for addition
     750     c) mask so adding will overflow
     751     d) mask for non 1 rhs
     752  */
     753  int maskA=0;
     754  int maskAdd=0;
     755  int maskC=0;
     756  int maskD=0;
     757  int i;
     758  for (i=0;i<numberElements;i++) {
     759    int iRow=rows[i];
     760    int numberBits = numberBits_[iRow];
     761    int startBit = startBit_[iRow];
     762    if (numberBits==1) {
     763      maskA |= 1<<startBit;
     764      maskAdd |= 1<<startBit;
     765    } else {
     766      int value=coefficients[i];
     767      int size = 1<<numberBits;
     768      int start = 1<<startBit;
     769      assert (value<size);
     770      maskAdd |= start*value;
     771      int gap = size-rhs_[iRow]-1;
     772      assert (gap>=0);
     773      maskC |= start*gap;
     774      maskD |= start*(size-1);
     775    }
     776  }
     777  int maskDiff = maskD-maskC;
     778  i=size_-1-maskAdd;
     779  bool touched = false;
     780  while (i>=0) {
     781    int kMask = i&maskA;
     782    if (kMask==0) {
     783      int added = i & maskD; // just bits belonging to non 1 rhs
     784      added += maskC; // will overflow mask if bad
     785      added &= (~maskD);
     786      if (added == 0) {
     787        double thisCost = cost_[i];
     788        if (thisCost!=COIN_DBL_MAX) {
     789          // possible
     790          double newCost=thisCost+cost;
     791          int next = i + maskAdd;
     792          if (cost_[next]>newCost) {
     793            cost_[next]=newCost;
     794            back_[next]=i;
     795            id_[next]=id;
     796            touched=true;
     797          }
     798        }
     799        i--;
     800      } else {
     801        // we can skip some
     802        int k = i & ~ maskD; // clear all
     803        // Put back enough - but only below where we are
     804        int kk=(numberNonOne_<<1)-2;
     805        assert (rhs_[kk]>1);
     806        int iMask=0;
     807        for(;kk>=0;kk-=2) {
     808          iMask = 1<<startBit_[kk+1];
     809          if ((added&iMask)!=0) {
     810            iMask--;
     811            break;
     812          }
     813        }
     814        assert (kk>=0);
     815        iMask &= maskDiff;
     816        k |= iMask;
     817        assert (k<i);
     818        i=k;
     819      }
     820    } else {
     821      // we can skip some
     822      // we can skip some
     823      int k=(i&~maskA)-1;
     824#ifdef CBC_DEBUG
     825      for (int j=i-1;j>k;j--) {
     826        int jMask = j&mask;
     827        assert (jMask!=0);
     828      }
     829#endif
     830      i=k;
     831    }
     832  }
     833  return touched;
     834}
    478835// update model
    479836void CbcFathomDynamicProgramming::setModel(CbcModel * model)
Note: See TracChangeset for help on using the changeset viewer.