Changeset 39


Ignore:
Timestamp:
Dec 13, 2004 11:53:08 AM (15 years ago)
Author:
forrest
Message:

fixing bug

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/CbcFathomDynamicProgramming.cpp

    r38 r39  
    2323   cost_(NULL),
    2424   back_(NULL),
    25    id_(NULL),
    2625   lookup_(NULL),
     26   indices_(NULL),
    2727   numberActive_(0),
    2828   maximumSizeAllowed_(1000000),
     
    3030   numberBits_(NULL),
    3131   rhs_(NULL),
    32    numberNonOne_(0)
     32   coefficients_(NULL),
     33   target_(0),
     34   numberNonOne_(0),
     35   bitPattern_(0),
     36   algorithm_(-1)
    3337{
    3438
     
    4044   cost_(NULL),
    4145   back_(NULL),
    42    id_(NULL),
    4346   lookup_(NULL),
     47   indices_(NULL),
    4448   numberActive_(0),
    4549   maximumSizeAllowed_(1000000),
     
    4751   numberBits_(NULL),
    4852   rhs_(NULL),
    49    numberNonOne_(0)
    50 {
    51   type_=gutsOfCheckPossible();
     53   coefficients_(NULL),
     54   target_(0),
     55   numberNonOne_(0),
     56   bitPattern_(0),
     57   algorithm_(-1)
     58{
     59  type_=checkPossible();
    5260}
    5361
     
    6371  delete [] cost_;
    6472  delete [] back_;
    65   delete [] id_;
    6673  delete [] lookup_;
     74  delete [] indices_;
    6775  delete [] startBit_;
    6876  delete [] numberBits_;
    6977  delete [] rhs_;
     78  delete [] coefficients_;
    7079  cost_ = NULL;
    7180  back_ = NULL;
    72   id_ = NULL;
    7381  lookup_ = NULL;
     82  indices_ =NULL;
    7483  startBit_ = NULL;
    7584  numberBits_ = NULL;
    7685  rhs_ = NULL;
     86  coefficients_ = NULL;
    7787}
    7888// Clone
     
    91101  cost_(NULL),
    92102  back_(NULL),
    93   id_(NULL),
    94103  lookup_(NULL),
     104  indices_(NULL),
    95105  numberActive_(rhs.numberActive_),
    96106  maximumSizeAllowed_(rhs.maximumSizeAllowed_),
     
    98108  numberBits_(NULL),
    99109  rhs_(NULL),
    100   numberNonOne_(rhs.numberNonOne_)
     110  coefficients_(NULL),
     111  target_(rhs.target_),
     112  numberNonOne_(rhs.numberNonOne_),
     113  bitPattern_(rhs.bitPattern_),
     114  algorithm_(rhs.algorithm_)
    101115{
    102116  if (size_) {
    103117    cost_=CoinCopyOfArray(rhs.cost_,size_);
    104118    back_=CoinCopyOfArray(rhs.back_,size_);
    105     id_=CoinCopyOfArray(rhs.id_,size_);
    106119    int numberRows=model_->getNumRows();
    107120    lookup_=CoinCopyOfArray(rhs.lookup_,numberRows);
    108121    startBit_=CoinCopyOfArray(rhs.startBit_,numberActive_);
     122    indices_=CoinCopyOfArray(rhs.indices_,numberActive_);
    109123    numberBits_=CoinCopyOfArray(rhs.numberBits_,numberActive_);
    110124    rhs_=CoinCopyOfArray(rhs.rhs_,numberActive_);
     125    coefficients_=CoinCopyOfArray(rhs.coefficients_,numberActive_);
    111126  }
    112127}
    113128// Returns type
    114129int
    115 CbcFathomDynamicProgramming::gutsOfCheckPossible(int allowableSize)
     130CbcFathomDynamicProgramming::checkPossible(int allowableSize)
    116131{
    117132  assert(model_->solver());
     
    227242      n01++;
    228243  }
     244  if (n01==numberColumns&&!nbadcoeff)
     245    algorithm_ = 0; // easiest
     246  else
     247    algorithm_ = 1;
    229248  if (allowableSize&&size_<=allowableSize) {
    230249    numberActive_=numberActive;
    231     cost_ = new double [size_];
    232     CoinFillN(cost_,size_,COIN_DBL_MAX);
     250    indices_ = new int [numberActive_];
     251    cost_ = new float [size_];
     252    CoinFillN(cost_,size_,FLT_MAX);
    233253    // but do nothing is okay
    234254    cost_[0]=0.0;
    235255    back_ = new int[size_];
    236256    CoinFillN(back_,size_,-1);
    237     id_ = new int[size_];
    238     CoinFillN(id_,size_,-1);
    239257    startBit_=new int[numberActive_];
    240258    numberBits_=new int[numberActive_];
     
    277295      }
    278296    }
     297    const double * rowLower = solver->getRowLower();
     298    if (algorithm_==0) {
     299      // rhs 1 and coefficients 1
     300      // Get first possible solution for printing
     301      target_=-1;
     302      int needed=0;
     303      int numberActive=0;
     304      for (i=0;i<numberRows;i++) {
     305        int newRow = lookup_[i];
     306        if (newRow>=0) {
     307          if (rowLower[i]==rowUpper[i]) {
     308            needed += 1<<numberActive;
     309            numberActive++;
     310          }
     311        }
     312      }
     313      for (i=0;i<size_;i++) {
     314        if ((i&needed)==needed) {
     315          break;
     316        }
     317      }
     318      target_=i;
     319    } else {
     320      coefficients_ = new int[numberActive_];
     321      // If not too many general rhs then we can be more efficient
     322      numberNonOne_=0;
     323      for (i=0;i<numberActive_;i++) {
     324        if (rhs_[i]!=1)
     325          numberNonOne_++;
     326      }
     327      if (numberNonOne_*2<numberActive_) {
     328        // put rhs >1 every second
     329        int * permute = new int[numberActive_];
     330        int * temp = new int[numberActive_];
     331        // try different ways
     332        int k=0;
     333        for (i=0;i<numberRows;i++) {
     334          int newRow = lookup_[i];
     335          if (newRow>=0&&rhs_[newRow]>1) {
     336            permute[newRow]=k;
     337            k +=2;
     338          }
     339        }
     340        // adjust so k points to last
     341        k -= 2;
     342        // and now rest
     343        int k1=1;
     344        for (i=0;i<numberRows;i++) {
     345          int newRow = lookup_[i];
     346          if (newRow>=0&&rhs_[newRow]==1) {
     347            permute[newRow]=k1;
     348            k1++;
     349            if (k1<=k)
     350              k1++;
     351          }
     352        }
     353        for (i=0;i<numberActive_;i++) {
     354          int put = permute[i];
     355          temp[put]=rhs_[i];
     356        }
     357        memcpy(rhs_,temp,numberActive_*sizeof(int));
     358        for (i=0;i<numberActive_;i++) {
     359          int put = permute[i];
     360          temp[put]=numberBits_[i];
     361        }
     362        memcpy(numberBits_,temp,numberActive_*sizeof(int));
     363        k=0;
     364        for (i=0;i<numberActive_;i++) {
     365          startBit_[i]=k;
     366          k+= numberBits_[i];
     367        }
     368        for (i=0;i<numberRows;i++) {
     369          int newRow = lookup_[i];
     370          if (newRow>=0)
     371            lookup_[i]=permute[newRow];
     372        }
     373        delete [] permute;
     374        delete [] temp;
     375        // mark new method
     376        algorithm_=2;
     377      }
     378      // Get first possible solution for printing
     379      target_=-1;
     380      int needed=0;
     381      int * lower2 = new int[numberActive_];
     382      for (i=0;i<numberRows;i++) {
     383        int newRow = lookup_[i];
     384        if (newRow>=0) {
     385          int gap=(int) (rowUpper[i]-CoinMax(0.0,rowLower[i]));
     386          lower2[newRow]=rhs_[newRow]-gap;
     387          int numberBits = numberBits_[newRow];
     388          int startBit = startBit_[newRow];
     389          if (numberBits==1&&!gap) {
     390            needed |= 1<<startBit;
     391          }
     392        }
     393      }
     394      for (i=0;i<size_;i++) {
     395        if ((i&needed)==needed) {
     396          // this one may do
     397          bool good = true;
     398          for (int kk=0;kk<numberActive_;kk++) {
     399            int numberBits = numberBits_[kk];
     400            int startBit = startBit_[kk];
     401            int size = 1<<numberBits;
     402            int start = 1<<startBit;
     403            int mask = start*(size-1);
     404            int level=(i&mask)>>startBit;
     405            if (level<lower2[kk]) {
     406              good=false;
     407              break;
     408            }
     409          }
     410          if (good) {
     411            break;
     412          }
     413        }
     414      }
     415      delete [] lower2;
     416      target_=i;
     417    }
    279418  }
    280419  delete [] rhs;
    281420  if (allowableSize&&size_>allowableSize) {
    282     printf("Too large - need %d entries x 16 bytes\n",size_);
     421    printf("Too large - need %d entries x 8 bytes\n",size_);
    283422    return -1; // too big
    284   }
    285   if (n01==numberColumns&&!nbadcoeff)
    286     return 0; // easiest
    287   else
    288     return 1;
     423  } else {
     424    return algorithm_;
     425  }
    289426}
    290427
     
    294431{
    295432  model_=model;
    296   type_=gutsOfCheckPossible();
     433  type_=checkPossible();
    297434}
    298435int
     
    300437{
    301438  int returnCode=0;
    302   int type=gutsOfCheckPossible(maximumSizeAllowed_);
    303   if (type>=0) {
    304     bool gotSolution=false;
     439  checkPossible(maximumSizeAllowed_);
     440  if (algorithm_>=0) {
    305441    OsiSolverInterface * solver = model_->solver();
    306442    const double * lower = solver->getColLower();
     
    319455   
    320456    int numberColumns = solver->getNumCols();
    321     int * indices = new int [numberActive_];
    322457    double offset;
    323458    solver->getDblParam(OsiObjOffset,offset);
     
    325460    int i;
    326461    // may be possible
    327     if (type==0) {
    328       // rhs 1 and coefficients 1
    329       for (i=0;i<numberColumns;i++) {
    330         int j;
    331         double lowerValue = lower[i];
    332         assert (lowerValue==floor(lowerValue));
    333         double cost = direction * objective[i];
    334         fixedObj += lowerValue*cost;
    335         if (lowerValue==upper[i])
    336           continue;
    337         int n=0;
    338         for (j=columnStart[i];
    339              j<columnStart[i]+columnLength[i];j++) {
    340           int iRow=row[j];
    341           double value = element[j];
    342           int newRow = lookup_[iRow];
    343           if (newRow<0||value>rhs_[newRow]) {
    344             n=0;
    345             break; //can't use
    346           } else {
    347             indices[n++]=newRow;
    348           }
    349         }
    350         if (n) {
    351           addOneColumn0(i,n,indices,cost);
    352         }
    353       }
    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     }
     462    double bestAtTarget = FLT_MAX;
     463    for (i=0;i<numberColumns;i++) {
     464      if (size_>10000000&&(i%100)==0)
     465        printf("column %d\n",i);
     466      double lowerValue = lower[i];
     467      assert (lowerValue==floor(lowerValue));
     468      float cost = direction * objective[i];
     469      fixedObj += lowerValue*cost;
     470      int gap = (int) (upper[i]-lowerValue);
     471      CoinBigIndex start = columnStart[i];
     472      tryColumn(columnLength[i],row+start,element+start,cost,gap);
     473      if (cost_[target_]<bestAtTarget) {
     474        printf("At column %d new best objective of %g\n",i,cost_[target_]);
     475        bestAtTarget = cost_[target_];
     476      }
     477    }
     478    returnCode=1;
    460479    int needed=0;
    461     double bestValue=COIN_DBL_MAX;
     480    float bestValue=FLT_MAX;
    462481    int iBest=-1;
    463     if (type==0) {
     482    if (algorithm_==0) {
    464483      int numberActive=0;
    465484      for (i=0;i<numberRows;i++) {
     
    519538      delete [] lower;
    520539    }
    521     if (bestValue<COIN_DBL_MAX) {
     540    if (bestValue<FLT_MAX) {
    522541      bestValue += fixedObj;
    523542      printf("Can get solution of %g\n",bestValue);
     
    527546        memcpy(betterSolution,lower,numberColumns*sizeof(double));
    528547        while (iBest>0) {
    529           int iColumn = id_[iBest];
     548          int n=decodeBitPattern(iBest-back_[iBest],indices_,numberRows);
     549          // Search for cheapest
     550          float bestCost=FLT_MAX;
     551          int iColumn=-1;
     552          for (i=0;i<numberColumns;i++) {
     553            if (n==columnLength[i]) {
     554              bool good=true;
     555              for (int j=columnStart[i];
     556                   j<columnStart[i]+columnLength[i];j++) {
     557                int iRow=row[j];
     558                double value = element[j];
     559                int iValue = (int) value;
     560                if(iValue!=indices_[iRow]) {
     561                  good=false;
     562                  break;
     563                }
     564              }
     565              if (good && objective[i]<bestCost) {
     566                bestCost=objective[i];
     567                iColumn=i;
     568              }
     569            }
     570          }
    530571          assert (iColumn>=0);
    531572          betterSolution[iColumn]++;
     
    534575        }
    535576      }
    536     }
    537     delete [] indices;
    538     gutsOfDelete();
    539     if (gotSolution) {
    540       int i;
    541577      // paranoid check
    542578      double * rowActivity = new double [numberRows];
     
    570606      delete [] rowActivity;
    571607    }
     608    gutsOfDelete();
    572609  }
    573610  return returnCode;
     611}
     612/* Tries a column
     613   returns true if was used in making any changes.
     614*/
     615bool
     616CbcFathomDynamicProgramming::tryColumn(int numberElements, const int * rows,
     617                                          const double * coefficients, float cost,
     618                                          int upper)
     619{
     620  bool touched=false;
     621  int n=0;
     622  if (algorithm_==0) {
     623    for (int j=0;j<numberElements;j++) {
     624      int iRow=rows[j];
     625      double value = coefficients[j];
     626      int newRow = lookup_[iRow];
     627      if (newRow<0||value>rhs_[newRow]) {
     628        n=0;
     629        break; //can't use
     630      } else {
     631        indices_[n++]=newRow;
     632      }
     633    }
     634    if (n&&upper) {
     635      touched=addOneColumn0(n,indices_,cost);
     636    }
     637  } else {
     638    for (int j=0;j<numberElements;j++) {
     639      int iRow=rows[j];
     640      double value = coefficients[j];
     641      int iValue = (int) value;
     642      int newRow = lookup_[iRow];
     643      if (newRow<0||iValue>rhs_[newRow]) {
     644        n=0;
     645        break; //can't use
     646      } else {
     647        coefficients_[n]=iValue;
     648        indices_[n++]=newRow;
     649        if (upper*iValue>rhs_[newRow]) {
     650          upper = rhs_[newRow]/iValue;
     651        }
     652      }
     653    }
     654    if (n) {
     655      if (algorithm_==1) {
     656        for (int k=1;k<=upper;k++) {
     657          bool t = addOneColumn1(n,indices_,coefficients_,cost);
     658          if (t)
     659            touched=true;
     660        }
     661      } else {
     662        CoinSort_2(indices_,indices_+n,coefficients_);
     663        for (int k=1;k<=upper;k++) {
     664          bool t = addOneColumn1A(n,indices_,coefficients_,cost);
     665          if (t)
     666            touched=true;
     667        }
     668      }
     669    }
     670  }
     671  return touched;
    574672}
    575673/* Adds one column if type 0,
     
    577675*/
    578676bool
    579 CbcFathomDynamicProgramming::addOneColumn0(int id,int numberElements, const int * rows,
    580                      double cost)
     677CbcFathomDynamicProgramming::addOneColumn0(int numberElements, const int * rows,
     678                     float cost)
    581679{
    582680  // build up mask
     
    587685    mask |= 1<<iRow;
    588686  }
     687  bitPattern_=mask;
    589688  i=size_-1-mask;
    590689  bool touched = false;
     
    592691    int kMask = i&mask;
    593692    if (kMask==0) {
    594       double thisCost = cost_[i];
    595       if (thisCost!=COIN_DBL_MAX) {
     693      float thisCost = cost_[i];
     694      if (thisCost!=FLT_MAX) {
    596695        // possible
    597         double newCost=thisCost+cost;
     696        float newCost=thisCost+cost;
    598697        int next = i + mask;
    599698        if (cost_[next]>newCost) {
    600699          cost_[next]=newCost;
    601700          back_[next]=i;
    602           id_[next]=id;
    603701          touched=true;
    604702        }
     
    607705    } else {
    608706      // we can skip some
    609       int k=(i&~mask)-1;
     707      int k=(i&~mask);
    610708#ifdef CBC_DEBUG
    611709      for (int j=i-1;j>k;j--) {
     
    624722*/
    625723bool
    626 CbcFathomDynamicProgramming::addOneColumn1(int id,int numberElements, const int * rows,
    627                                            const int * coefficients, double cost)
     724CbcFathomDynamicProgramming::addOneColumn1(int numberElements, const int * rows,
     725                                           const int * coefficients, float cost)
    628726{
    629727  /* build up masks.
     
    667765    }
    668766  }
     767  bitPattern_=maskAdd;
    669768  i=size_-1-maskAdd;
    670769  bool touched = false;
     
    708807      if (good) {
    709808        double thisCost = cost_[i];
    710         if (thisCost!=COIN_DBL_MAX) {
     809        if (thisCost!=FLT_MAX) {
    711810          // possible
    712811          double newCost=thisCost+cost;
     
    715814            cost_[next]=newCost;
    716815            back_[next]=i;
    717             id_[next]=id;
    718816            touched=true;
    719817          }
     
    724822      // we can skip some
    725823      // we can skip some
    726       int k=(i&~mask1)-1;
     824      int k=(i&~mask1);
    727825#ifdef CBC_DEBUG
    728826      for (int j=i-1;j>k;j--) {
     
    742840*/
    743841bool
    744 CbcFathomDynamicProgramming::addOneColumn1A(int id,int numberElements, const int * rows,
    745                                            const int * coefficients, double cost)
     842CbcFathomDynamicProgramming::addOneColumn1A(int numberElements, const int * rows,
     843                                           const int * coefficients, float cost)
    746844{
    747845  /* build up masks.
     
    775873    }
    776874  }
     875  bitPattern_=maskAdd;
    777876  int maskDiff = maskD-maskC;
    778877  i=size_-1-maskAdd;
    779878  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) {
     879  if (!maskD) {
     880    // Just ones
     881    while (i>=0) {
     882      int kMask = i&maskA;
     883      if (kMask==0) {
     884        float thisCost = cost_[i];
     885        if (thisCost!=FLT_MAX) {
    789886          // possible
    790           double newCost=thisCost+cost;
     887          float newCost=thisCost+cost;
    791888          int next = i + maskAdd;
    792889          if (cost_[next]>newCost) {
    793890            cost_[next]=newCost;
    794891            back_[next]=i;
    795             id_[next]=id;
    796892            touched=true;
    797893          }
     
    799895        i--;
    800896      } 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);
     897        // we can skip some
     898        int k=(i&~maskA);
    818899        i=k;
    819900      }
    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;
     901    }
     902  } else {
     903    // More general
     904    while (i>=0) {
     905      int kMask = i&maskA;
     906      if (kMask==0) {
     907        int added = i & maskD; // just bits belonging to non 1 rhs
     908        added += maskC; // will overflow mask if bad
     909        added &= (~maskD);
     910        if (added == 0) {
     911          float thisCost = cost_[i];
     912          if (thisCost!=FLT_MAX) {
     913            // possible
     914            float newCost=thisCost+cost;
     915            int next = i + maskAdd;
     916            if (cost_[next]>newCost) {
     917              cost_[next]=newCost;
     918              back_[next]=i;
     919              touched=true;
     920            }
     921          }
     922          i--;
     923        } else {
     924          // we can skip some
     925          int k = i & ~ maskD; // clear all
     926          // Put back enough - but only below where we are
     927          int kk=(numberNonOne_<<1)-2;
     928          assert (rhs_[kk]>1);
     929          int iMask=0;
     930          for(;kk>=0;kk-=2) {
     931            iMask = 1<<startBit_[kk+1];
     932            if ((added&iMask)!=0) {
     933              iMask--;
     934              break;
     935            }
     936          }
     937          assert (kk>=0);
     938          iMask &= maskDiff;
     939          k |= iMask;
     940          assert (k<i);
     941          i=k;
     942        }
     943      } else {
     944        // we can skip some
     945        int k=(i&~maskA);
     946        i=k;
     947      }
    831948    }
    832949  }
     
    837954{
    838955  model_ = model;
    839   type_=gutsOfCheckPossible();
     956  type_=checkPossible();
     957}
     958// Gets bit pattern from original column
     959int CbcFathomDynamicProgramming::bitPattern(int numberElements, const int * rows,
     960                                            const int * coefficients)
     961{
     962  int i;
     963  int mask=0;
     964  switch (algorithm_) {
     965    // just ones
     966  case 0:
     967    for (i=0;i<numberElements;i++) {
     968      int iRow=rows[i];
     969      iRow = lookup_[iRow];
     970      if (iRow>=0)
     971        mask |= 1<<iRow;
     972    }
     973    break;
     974    //
     975  case 1:
     976  case 2:
     977    for (i=0;i<numberElements;i++) {
     978      int iRow=rows[i];
     979      iRow = lookup_[iRow];
     980      if (iRow>=0) {
     981        int startBit = startBit_[iRow];
     982        int value=coefficients[i];
     983        int start = 1<<startBit;
     984        mask |= start*value;
     985      }
     986    }
     987    break;
     988  }
     989  return mask;
     990}
     991// Fills in original column (dense) from bit pattern
     992int CbcFathomDynamicProgramming::decodeBitPattern(int bitPattern,
     993                                                   int * values,
     994                                                   int numberRows)
     995{
     996  int i;
     997  int n=0;
     998  switch (algorithm_) {
     999    // just ones
     1000  case 0:
     1001    for (i=0;i<numberRows;i++) {
     1002      values[i]=0;
     1003      int iRow = lookup_[i];
     1004      if (iRow>=0) {
     1005        if ((bitPattern&(1<<iRow))!=0) {
     1006          values[i]=1;
     1007          n++;
     1008        }
     1009      }
     1010    }
     1011    break;
     1012    //
     1013  case 1:
     1014  case 2:
     1015    for (i=0;i<numberRows;i++) {
     1016      values[i]=0;
     1017      int iRow = lookup_[i];
     1018      if (iRow>=0) {
     1019        int startBit = startBit_[iRow];
     1020        int numberBits = numberBits_[iRow];
     1021        int iValue = bitPattern>>startBit;
     1022        iValue &= ((1<<numberBits)-1);
     1023        if (iValue) {
     1024          values[i]=iValue;
     1025          n++;
     1026        }
     1027      }
     1028    }
     1029    break;
     1030  }
     1031  return n;
    8401032}
    8411033
  • trunk/include/CbcFathomDynamicProgramming.hpp

    r38 r39  
    77
    88//#############################################################################
    9 /** FathomDynamicProgramming base class.
     9/** FathomDynamicProgramming class.
    1010
    1111    The idea is that after some branching the problem will be effectively smaller than
     
    1313    fathom this branch quickly.
    1414
    15     One method is to presolve the problem to give a much smaller new problem and then do branch
    16     and cut on that.  Another might be dynamic programming.
    17 
     15    This is a dynamic programming implementation which is very fast for some
     16    specialized problems.  It expects small integral rhs, an all integer problem
     17    and positive integral coefficients. At present it can not do general set covering
     18    problems just set partitioning.  It can find multiple optima for various rhs
     19    combinations.
     20   
     21    The main limiting factor is size of state space.  Each 1 rhs doubles the size of the problem.
     22    2 or 3 rhs quadruples, 4,5,6,7 by 8 etc.
    1823 */
    1924
     
    5459  inline void setMaximumSize(int value)
    5560  { maximumSizeAllowed_=value;};
     61  /// Returns type of algorithm and sets up arrays
     62  int checkPossible(int allowableSize=0);
     63  // set algorithm
     64  inline void setAlgorithm(int value)
     65  { algorithm_=value;};
     66  /** Tries a column
     67      returns true if was used in making any changes.
     68  */
     69  bool tryColumn(int numberElements, const int * rows,
     70                    const double * coefficients, float cost,
     71                    int upper=INT_MAX);
     72  /// Returns cost array
     73  inline const float * cost() const
     74  { return cost_;};
     75  /// Returns back array
     76  inline const int * back() const
     77  { return back_;};
     78  /// Gets bit pattern for target result
     79  inline int target() const
     80  { return target_;};
     81  /// Sets bit pattern for target result
     82  inline void setTarget(int value)
     83  { target_=value;};
    5684private:
    57   /// Returns type
    58   int gutsOfCheckPossible(int allowableSize=0);
    5985  /// Does deleteions
    6086  void gutsOfDelete();
     
    6389      returns true if was used in making any changes
    6490  */
    65   bool addOneColumn0(int id,int numberElements, const int * rows,
    66                      double cost);
     91  bool addOneColumn0(int numberElements, const int * rows,
     92                     float cost);
    6793  /** Adds one attempt of one column of type 1,
    6894      returns true if was used in making any changes.
    6995      At present the user has to call it once for each possible value
    7096  */
    71   bool addOneColumn1(int id,int numberElements, const int * rows,
    72                      const int * coefficients, double cost);
     97  bool addOneColumn1(int numberElements, const int * rows,
     98                     const int * coefficients, float cost);
    7399  /** Adds one attempt of one column of type 1,
    74100      returns true if was used in making any changes.
     
    76102      This version is when there are enough 1 rhs to do faster
    77103  */
    78   bool addOneColumn1A(int id,int numberElements, const int * rows,
    79                      const int * coefficients, double cost);
     104  bool addOneColumn1A(int numberElements, const int * rows,
     105                     const int * coefficients, float cost);
     106  /// Gets bit pattern from original column
     107  int bitPattern(int numberElements, const int * rows,
     108                     const int * coefficients);
     109  /// Gets bit pattern from original column
     110  int bitPattern(int numberElements, const int * rows,
     111                     const double * coefficients);
     112  /// Fills in original column (dense) from bit pattern - returning number nonzero
     113  int decodeBitPattern(int bitPattern, int * values, int numberRows);
    80114
    81115protected:
     
    87121  */
    88122  int type_;
    89   /// Space for states (? could be float)
    90   double * cost_;
     123  /// Space for states (float as who cares)
     124  float * cost_;
    91125  /// Which state produced this cheapest one
    92126  int * back_;
    93   /// Some id as to which variable was used
    94   int * id_;
    95127  /// Some rows may be satisified so we need a lookup
    96128  int * lookup_;
     129  /// Space for sorted indices
     130  int * indices_;
    97131  /// Number of active rows
    98132  int numberActive_;
     
    105139  /// Effective rhs
    106140  int * rhs_;
     141  /// Space for sorted coefficients
     142  int * coefficients_;
     143  /// Target pattern
     144  int target_;
    107145  /// Number of Non 1 rhs
    108146  int numberNonOne_;
     147  /// Current bit pattern
     148  int bitPattern_;
     149  /// Current algorithm
     150  int algorithm_;
    109151private:
    110152 
Note: See TracChangeset for help on using the changeset viewer.