Ignore:
Timestamp:
Nov 13, 2009 5:50:47 PM (10 years ago)
Author:
EdwinStraver
Message:

Moved the objects CbcCutBranchingObject?, CbcBranchToFixLots?, and CbcBranchAllDifferent? to new files, out of CbcBranchCut?.

Location:
branches/sandbox/Cbc/src
Files:
6 added
2 edited

Legend:

Unmodified
Added
Removed
  • branches/sandbox/Cbc/src/CbcBranchCut.cpp

    r1286 r1305  
    125125}
    126126
    127 
    128 // Default Constructor
    129 CbcCutBranchingObject::CbcCutBranchingObject()
    130         : CbcBranchingObject()
    131 {
    132     down_ = OsiRowCut();
    133     up_ = OsiRowCut();
    134     canFix_ = false;
    135 }
    136 
    137 // Useful constructor
    138 CbcCutBranchingObject::CbcCutBranchingObject (CbcModel * model,
    139         OsiRowCut & down,
    140         OsiRowCut &up,
    141         bool canFix)
    142         : CbcBranchingObject(model, 0, -1, 0.0)
    143 {
    144     down_ = down;
    145     up_ = up;
    146     canFix_ = canFix;
    147 }
    148 
    149 // Copy constructor
    150 CbcCutBranchingObject::CbcCutBranchingObject ( const CbcCutBranchingObject & rhs) : CbcBranchingObject(rhs)
    151 {
    152     down_ = rhs.down_;
    153     up_ = rhs.up_;
    154     canFix_ = rhs.canFix_;
    155 }
    156 
    157 // Assignment operator
    158 CbcCutBranchingObject &
    159 CbcCutBranchingObject::operator=( const CbcCutBranchingObject & rhs)
    160 {
    161     if (this != &rhs) {
    162         CbcBranchingObject::operator=(rhs);
    163         down_ = rhs.down_;
    164         up_ = rhs.up_;
    165         canFix_ = rhs.canFix_;
    166     }
    167     return *this;
    168 }
    169 CbcBranchingObject *
    170 CbcCutBranchingObject::clone() const
    171 {
    172     return (new CbcCutBranchingObject(*this));
    173 }
    174 
    175 
    176 // Destructor
    177 CbcCutBranchingObject::~CbcCutBranchingObject ()
    178 {
    179 }
    180 
    181 /*
    182   Perform a branch by adjusting bounds and/or adding a cut. Note
    183   that each arm of the branch advances the object to the next arm by
    184   advancing the value of way_.
    185 
    186   Returns change in guessed objective on next branch
    187 */
    188 double
    189 CbcCutBranchingObject::branch()
    190 {
    191     decrementNumberBranchesLeft();
    192     OsiRowCut * cut;
    193     if (way_ < 0) {
    194         cut = &down_;
    195         way_ = 1;
    196     } else {
    197         cut = &up_;
    198         way_ = -1;        // Swap direction
    199     }
    200     printf("CUT %s ", (way_ == -1) ? "up" : "down");
    201     cut->print();
    202     // See if cut just fixes variables
    203     double lb = cut->lb();
    204     double ub = cut->ub();
    205     int n = cut->row().getNumElements();
    206     const int * column = cut->row().getIndices();
    207     const double * element = cut->row().getElements();
    208     OsiSolverInterface * solver = model_->solver();
    209     const double * upper = solver->getColUpper();
    210     const double * lower = solver->getColLower();
    211     double low = 0.0;
    212     double high = 0.0;
    213     for (int i = 0; i < n; i++) {
    214         int iColumn = column[i];
    215         double value = element[i];
    216         if (value > 0.0) {
    217             high += upper[iColumn] * value;
    218             low += lower[iColumn] * value;
    219         } else {
    220             high += lower[iColumn] * value;
    221             low += upper[iColumn] * value;
    222         }
    223     }
    224     // leave as cut
    225     //model_->setNextRowCut(*cut);
    226     //return 0.0;
    227     // assume cut was cunningly constructed so we need not worry too much about tolerances
    228     if (low + 1.0e-8 >= ub && canFix_) {
    229         // fix
    230         for (int i = 0; i < n; i++) {
    231             int iColumn = column[i];
    232             double value = element[i];
    233             if (value > 0.0) {
    234                 solver->setColUpper(iColumn, lower[iColumn]);
    235             } else {
    236                 solver->setColLower(iColumn, upper[iColumn]);
    237             }
    238         }
    239     } else if (high - 1.0e-8 <= lb && canFix_) {
    240         // fix
    241         for (int i = 0; i < n; i++) {
    242             int iColumn = column[i];
    243             double value = element[i];
    244             if (value > 0.0) {
    245                 solver->setColLower(iColumn, upper[iColumn]);
    246             } else {
    247                 solver->setColUpper(iColumn, lower[iColumn]);
    248             }
    249         }
    250     } else {
    251         // leave as cut
    252         model_->setNextRowCut(*cut);
    253     }
    254     return 0.0;
    255 }
    256 // Print what would happen
    257 void
    258 CbcCutBranchingObject::print()
    259 {
    260     OsiRowCut * cut;
    261     if (way_ < 0) {
    262         cut = &down_;
    263         printf("CbcCut would branch down");
    264     } else {
    265         cut = &up_;
    266         printf("CbcCut would branch up");
    267     }
    268     double lb = cut->lb();
    269     double ub = cut->ub();
    270     int n = cut->row().getNumElements();
    271     const int * column = cut->row().getIndices();
    272     const double * element = cut->row().getElements();
    273     if (n > 5) {
    274         printf(" - %d elements, lo=%g, up=%g\n", n, lb, ub);
    275     } else {
    276         printf(" - %g <=", lb);
    277         for (int i = 0; i < n; i++) {
    278             int iColumn = column[i];
    279             double value = element[i];
    280             printf(" (%d,%g)", iColumn, value);
    281         }
    282         printf(" <= %g\n", ub);
    283     }
    284 }
    285 
    286 // Return true if branch should fix variables
    287 bool
    288 CbcCutBranchingObject::boundBranch() const
    289 {
    290     return false;
    291 }
    292 
    293 /** Compare the original object of \c this with the original object of \c
    294     brObj. Assumes that there is an ordering of the original objects.
    295     This method should be invoked only if \c this and brObj are of the same
    296     type.
    297     Return negative/0/positive depending on whether \c this is
    298     smaller/same/larger than the argument.
    299 */
    300 int
    301 CbcCutBranchingObject::compareOriginalObject
    302 (const CbcBranchingObject* brObj) const
    303 {
    304     const CbcCutBranchingObject* br =
    305         dynamic_cast<const CbcCutBranchingObject*>(brObj);
    306     assert(br);
    307     const OsiRowCut& r0 = way_ == -1 ? down_ : up_;
    308     const OsiRowCut& r1 = br->way_ == -1 ? br->down_ : br->up_;
    309     return r0.row().compare(r1.row());
    310 }
    311 
    312 /** Compare the \c this with \c brObj. \c this and \c brObj must be os the
    313     same type and must have the same original object, but they may have
    314     different feasible regions.
    315     Return the appropriate CbcRangeCompare value (first argument being the
    316     sub/superset if that's the case). In case of overlap (and if \c
    317     replaceIfOverlap is true) replace the current branching object with one
    318     whose feasible region is the overlap.
    319 */
    320 
    321 CbcRangeCompare
    322 CbcCutBranchingObject::compareBranchingObject
    323 (const CbcBranchingObject* brObj, const bool replaceIfOverlap)
    324 {
    325     const CbcCutBranchingObject* br =
    326         dynamic_cast<const CbcCutBranchingObject*>(brObj);
    327     assert(br);
    328     OsiRowCut& r0 = way_ == -1 ? down_ : up_;
    329     const OsiRowCut& r1 = br->way_ == -1 ? br->down_ : br->up_;
    330     double thisBd[2];
    331     thisBd[0] = r0.lb();
    332     thisBd[1] = r0.ub();
    333     double otherBd[2];
    334     otherBd[0] = r1.lb();
    335     otherBd[1] = r1.ub();
    336     CbcRangeCompare comp = CbcCompareRanges(thisBd, otherBd, replaceIfOverlap);
    337     if (comp != CbcRangeOverlap || (comp == CbcRangeOverlap && !replaceIfOverlap)) {
    338         return comp;
    339     }
    340     r0.setLb(thisBd[0]);
    341     r0.setUb(thisBd[1]);
    342     return comp;
    343 }
    344 
    345 //##############################################################################
    346 
    347 /** Default Constructor
    348 
    349   Equivalent to an unspecified binary variable.
    350 */
    351 CbcBranchToFixLots::CbcBranchToFixLots ()
    352         : CbcBranchCut(),
    353         djTolerance_(COIN_DBL_MAX),
    354         fractionFixed_(1.0),
    355         mark_(NULL),
    356         depth_(-1),
    357         numberClean_(0),
    358         alwaysCreate_(false)
    359 {
    360 }
    361 
    362 /* Useful constructor - passed reduced cost tolerance and fraction we would like fixed.
    363    Also depth level to do at.
    364    Also passed number of 1 rows which when clean triggers fix
    365    Always does if all 1 rows cleaned up and number>0 or if fraction columns reached
    366    Also whether to create branch if can't reach fraction.
    367 */
    368 CbcBranchToFixLots::CbcBranchToFixLots (CbcModel * model, double djTolerance,
    369                                         double fractionFixed, int depth,
    370                                         int numberClean,
    371                                         const char * mark, bool alwaysCreate)
    372         : CbcBranchCut(model)
    373 {
    374     djTolerance_ = djTolerance;
    375     fractionFixed_ = fractionFixed;
    376     if (mark) {
    377         int numberColumns = model->getNumCols();
    378         mark_ = new char[numberColumns];
    379         memcpy(mark_, mark, numberColumns);
    380     } else {
    381         mark_ = NULL;
    382     }
    383     depth_ = depth;
    384     assert (model);
    385     OsiSolverInterface * solver = model_->solver();
    386     matrixByRow_ = *solver->getMatrixByRow();
    387     numberClean_ = numberClean;
    388     alwaysCreate_ = alwaysCreate;
    389 }
    390 // Copy constructor
    391 CbcBranchToFixLots::CbcBranchToFixLots ( const CbcBranchToFixLots & rhs)
    392         : CbcBranchCut(rhs)
    393 {
    394     djTolerance_ = rhs.djTolerance_;
    395     fractionFixed_ = rhs.fractionFixed_;
    396     int numberColumns = model_->getNumCols();
    397     mark_ = CoinCopyOfArray(rhs.mark_, numberColumns);
    398     matrixByRow_ = rhs.matrixByRow_;
    399     depth_ = rhs.depth_;
    400     numberClean_ = rhs.numberClean_;
    401     alwaysCreate_ = rhs.alwaysCreate_;
    402 }
    403 
    404 // Clone
    405 CbcObject *
    406 CbcBranchToFixLots::clone() const
    407 {
    408     return new CbcBranchToFixLots(*this);
    409 }
    410 
    411 // Assignment operator
    412 CbcBranchToFixLots &
    413 CbcBranchToFixLots::operator=( const CbcBranchToFixLots & rhs)
    414 {
    415     if (this != &rhs) {
    416         CbcBranchCut::operator=(rhs);
    417         djTolerance_ = rhs.djTolerance_;
    418         fractionFixed_ = rhs.fractionFixed_;
    419         int numberColumns = model_->getNumCols();
    420         delete [] mark_;
    421         mark_ = CoinCopyOfArray(rhs.mark_, numberColumns);
    422         matrixByRow_ = rhs.matrixByRow_;
    423         depth_ = rhs.depth_;
    424         numberClean_ = rhs.numberClean_;
    425         alwaysCreate_ = rhs.alwaysCreate_;
    426     }
    427     return *this;
    428 }
    429 
    430 // Destructor
    431 CbcBranchToFixLots::~CbcBranchToFixLots ()
    432 {
    433     delete [] mark_;
    434 }
    435 CbcBranchingObject *
    436 CbcBranchToFixLots::createCbcBranch(OsiSolverInterface * solver, const OsiBranchingInformation * /*info*/, int /*way*/)
    437 {
    438     // by default way must be -1
    439     //assert (way==-1);
    440     //OsiSolverInterface * solver = model_->solver();
    441     const double * solution = model_->testSolution();
    442     const double * lower = solver->getColLower();
    443     const double * upper = solver->getColUpper();
    444     const double * dj = solver->getReducedCost();
    445     int i;
    446     int numberIntegers = model_->numberIntegers();
    447     const int * integerVariable = model_->integerVariable();
    448     double integerTolerance =
    449         model_->getDblParam(CbcModel::CbcIntegerTolerance);
    450     // make smaller ?
    451     double tolerance = CoinMin(1.0e-8, integerTolerance);
    452     // How many fixed are we aiming at
    453     int wantedFixed = static_cast<int> (static_cast<double>(numberIntegers) * fractionFixed_);
    454     int nSort = 0;
    455     int numberFixed = 0;
    456     int numberColumns = solver->getNumCols();
    457     int * sort = new int[numberColumns];
    458     double * dsort = new double[numberColumns];
    459     if (djTolerance_ != -1.234567) {
    460         int type = shallWe();
    461         assert (type);
    462         // Take clean first
    463         if (type == 1) {
    464             for (i = 0; i < numberIntegers; i++) {
    465                 int iColumn = integerVariable[i];
    466                 if (upper[iColumn] > lower[iColumn]) {
    467                     if (!mark_ || !mark_[iColumn]) {
    468                         if (solution[iColumn] < lower[iColumn] + tolerance) {
    469                             if (dj[iColumn] > djTolerance_) {
    470                                 dsort[nSort] = -dj[iColumn];
    471                                 sort[nSort++] = iColumn;
    472                             }
    473                         } else if (solution[iColumn] > upper[iColumn] - tolerance) {
    474                             if (dj[iColumn] < -djTolerance_) {
    475                                 dsort[nSort] = dj[iColumn];
    476                                 sort[nSort++] = iColumn;
    477                             }
    478                         }
    479                     }
    480                 } else {
    481                     numberFixed++;
    482                 }
    483             }
    484             // sort
    485             CoinSort_2(dsort, dsort + nSort, sort);
    486             nSort = CoinMin(nSort, wantedFixed - numberFixed);
    487         } else if (type < 10) {
    488             int i;
    489             //const double * rowLower = solver->getRowLower();
    490             const double * rowUpper = solver->getRowUpper();
    491             // Row copy
    492             const double * elementByRow = matrixByRow_.getElements();
    493             const int * column = matrixByRow_.getIndices();
    494             const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
    495             const int * rowLength = matrixByRow_.getVectorLengths();
    496             const double * columnLower = solver->getColLower();
    497             const double * columnUpper = solver->getColUpper();
    498             const double * solution = solver->getColSolution();
    499             int numberColumns = solver->getNumCols();
    500             int numberRows = solver->getNumRows();
    501             for (i = 0; i < numberColumns; i++) {
    502                 sort[i] = i;
    503                 if (columnLower[i] != columnUpper[i]) {
    504                     dsort[i] = 1.0e100;
    505                 } else {
    506                     dsort[i] = 1.0e50;
    507                     numberFixed++;
    508                 }
    509             }
    510             for (i = 0; i < numberRows; i++) {
    511                 double rhsValue = rowUpper[i];
    512                 bool oneRow = true;
    513                 // check elements
    514                 int numberUnsatisfied = 0;
    515                 for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
    516                     int iColumn = column[j];
    517                     double value = elementByRow[j];
    518                     double solValue = solution[iColumn];
    519                     if (columnLower[iColumn] != columnUpper[iColumn]) {
    520                         if (solValue < 1.0 - integerTolerance && solValue > integerTolerance)
    521                             numberUnsatisfied++;
    522                         if (value != 1.0) {
    523                             oneRow = false;
    524                             break;
    525                         }
    526                     } else {
    527                         rhsValue -= value * floor(solValue + 0.5);
    528                     }
    529                 }
    530                 if (oneRow && rhsValue <= 1.0 + tolerance) {
    531                     if (!numberUnsatisfied) {
    532                         for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
    533                             int iColumn = column[j];
    534                             if (dsort[iColumn] > 1.0e50) {
    535                                 dsort[iColumn] = 0;
    536                                 nSort++;
    537                             }
    538                         }
    539                     }
    540                 }
    541             }
    542             // sort
    543             CoinSort_2(dsort, dsort + numberColumns, sort);
    544         } else {
    545             // new way
    546             for (i = 0; i < numberIntegers; i++) {
    547                 int iColumn = integerVariable[i];
    548                 if (upper[iColumn] > lower[iColumn]) {
    549                     if (!mark_ || !mark_[iColumn]) {
    550                         double distanceDown = solution[iColumn] - lower[iColumn];
    551                         double distanceUp = upper[iColumn] - solution[iColumn];
    552                         double distance = CoinMin(distanceDown, distanceUp);
    553                         if (distance > 0.001 && distance < 0.5) {
    554                             dsort[nSort] = distance;
    555                             sort[nSort++] = iColumn;
    556                         }
    557                     }
    558                 }
    559             }
    560             // sort
    561             CoinSort_2(dsort, dsort + nSort, sort);
    562             int n = 0;
    563             double sum = 0.0;
    564             for (int k = 0; k < nSort; k++) {
    565                 sum += dsort[k];
    566                 if (sum <= djTolerance_)
    567                     n = k;
    568                 else
    569                     break;
    570             }
    571             nSort = CoinMin(n, numberClean_ / 1000000);
    572         }
    573     } else {
    574 #define FIX_IF_LESS -0.1
    575         // 3 in same row and sum <FIX_IF_LESS?
    576         int numberRows = matrixByRow_.getNumRows();
    577         const double * solution = model_->testSolution();
    578         const int * column = matrixByRow_.getIndices();
    579         const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
    580         const int * rowLength = matrixByRow_.getVectorLengths();
    581         double bestSum = 1.0;
    582         int nBest = -1;
    583         int kRow = -1;
    584         OsiSolverInterface * solver = model_->solver();
    585         for (int i = 0; i < numberRows; i++) {
    586             int numberUnsatisfied = 0;
    587             double sum = 0.0;
    588             for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
    589                 int iColumn = column[j];
    590                 if (solver->isInteger(iColumn)) {
    591                     double solValue = solution[iColumn];
    592                     if (solValue > 1.0e-5 && solValue < FIX_IF_LESS) {
    593                         numberUnsatisfied++;
    594                         sum += solValue;
    595                     }
    596                 }
    597             }
    598             if (numberUnsatisfied >= 3 && sum < FIX_IF_LESS) {
    599                 // possible
    600                 if (numberUnsatisfied > nBest ||
    601                         (numberUnsatisfied == nBest && sum < bestSum)) {
    602                     nBest = numberUnsatisfied;
    603                     bestSum = sum;
    604                     kRow = i;
    605                 }
    606             }
    607         }
    608         assert (nBest > 0);
    609         for (int j = rowStart[kRow]; j < rowStart[kRow] + rowLength[kRow]; j++) {
    610             int iColumn = column[j];
    611             if (solver->isInteger(iColumn)) {
    612                 double solValue = solution[iColumn];
    613                 if (solValue > 1.0e-5 && solValue < FIX_IF_LESS) {
    614                     sort[nSort++] = iColumn;
    615                 }
    616             }
    617         }
    618     }
    619     OsiRowCut down;
    620     down.setLb(-COIN_DBL_MAX);
    621     double rhs = 0.0;
    622     for (i = 0; i < nSort; i++) {
    623         int iColumn = sort[i];
    624         double distanceDown = solution[iColumn] - lower[iColumn];
    625         double distanceUp = upper[iColumn] - solution[iColumn];
    626         if (distanceDown < distanceUp) {
    627             rhs += lower[iColumn];
    628             dsort[i] = 1.0;
    629         } else {
    630             rhs -= upper[iColumn];
    631             dsort[i] = -1.0;
    632         }
    633     }
    634     down.setUb(rhs);
    635     down.setRow(nSort, sort, dsort);
    636     down.setEffectiveness(COIN_DBL_MAX); // so will persist
    637     delete [] sort;
    638     delete [] dsort;
    639     // up is same - just with rhs changed
    640     OsiRowCut up = down;
    641     up.setLb(rhs + 1.0);
    642     up.setUb(COIN_DBL_MAX);
    643     // Say can fix one way
    644     CbcCutBranchingObject * newObject =
    645         new CbcCutBranchingObject(model_, down, up, true);
    646     if (model_->messageHandler()->logLevel() > 1)
    647         printf("creating cut in CbcBranchCut\n");
    648     return newObject;
    649 }
    650 /* Does a lot of the work,
    651    Returns 0 if no good, 1 if dj, 2 if clean, 3 if both
    652    10 if branching on ones away from bound
    653 */
    654 int
    655 CbcBranchToFixLots::shallWe() const
    656 {
    657     int returnCode = 0;
    658     OsiSolverInterface * solver = model_->solver();
    659     int numberRows = matrixByRow_.getNumRows();
    660     //if (numberRows!=solver->getNumRows())
    661     //return 0;
    662     const double * solution = model_->testSolution();
    663     const double * lower = solver->getColLower();
    664     const double * upper = solver->getColUpper();
    665     const double * dj = solver->getReducedCost();
    666     int i;
    667     int numberIntegers = model_->numberIntegers();
    668     const int * integerVariable = model_->integerVariable();
    669     if (numberClean_ > 1000000) {
    670         int wanted = numberClean_ % 1000000;
    671         int * sort = new int[numberIntegers];
    672         double * dsort = new double[numberIntegers];
    673         int nSort = 0;
    674         for (i = 0; i < numberIntegers; i++) {
    675             int iColumn = integerVariable[i];
    676             if (upper[iColumn] > lower[iColumn]) {
    677                 if (!mark_ || !mark_[iColumn]) {
    678                     double distanceDown = solution[iColumn] - lower[iColumn];
    679                     double distanceUp = upper[iColumn] - solution[iColumn];
    680                     double distance = CoinMin(distanceDown, distanceUp);
    681                     if (distance > 0.001 && distance < 0.5) {
    682                         dsort[nSort] = distance;
    683                         sort[nSort++] = iColumn;
    684                     }
    685                 }
    686             }
    687         }
    688         // sort
    689         CoinSort_2(dsort, dsort + nSort, sort);
    690         int n = 0;
    691         double sum = 0.0;
    692         for (int k = 0; k < nSort; k++) {
    693             sum += dsort[k];
    694             if (sum <= djTolerance_)
    695                 n = k;
    696             else
    697                 break;
    698         }
    699         delete [] sort;
    700         delete [] dsort;
    701         return (n >= wanted) ? 10 : 0;
    702     }
    703     double integerTolerance =
    704         model_->getDblParam(CbcModel::CbcIntegerTolerance);
    705     // make smaller ?
    706     double tolerance = CoinMin(1.0e-8, integerTolerance);
    707     // How many fixed are we aiming at
    708     int wantedFixed = static_cast<int> (static_cast<double>(numberIntegers) * fractionFixed_);
    709     if (djTolerance_ < 1.0e10) {
    710         int nSort = 0;
    711         int numberFixed = 0;
    712         for (i = 0; i < numberIntegers; i++) {
    713             int iColumn = integerVariable[i];
    714             if (upper[iColumn] > lower[iColumn]) {
    715                 if (!mark_ || !mark_[iColumn]) {
    716                     if (solution[iColumn] < lower[iColumn] + tolerance) {
    717                         if (dj[iColumn] > djTolerance_) {
    718                             nSort++;
    719                         }
    720                     } else if (solution[iColumn] > upper[iColumn] - tolerance) {
    721                         if (dj[iColumn] < -djTolerance_) {
    722                             nSort++;
    723                         }
    724                     }
    725                 }
    726             } else {
    727                 numberFixed++;
    728             }
    729         }
    730         if (numberFixed + nSort < wantedFixed && !alwaysCreate_) {
    731             returnCode = 0;
    732         } else if (numberFixed < wantedFixed) {
    733             returnCode = 1;
    734         } else {
    735             returnCode = 0;
    736         }
    737     }
    738     if (numberClean_) {
    739         // see how many rows clean
    740         int i;
    741         //const double * rowLower = solver->getRowLower();
    742         const double * rowUpper = solver->getRowUpper();
    743         // Row copy
    744         const double * elementByRow = matrixByRow_.getElements();
    745         const int * column = matrixByRow_.getIndices();
    746         const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
    747         const int * rowLength = matrixByRow_.getVectorLengths();
    748         const double * columnLower = solver->getColLower();
    749         const double * columnUpper = solver->getColUpper();
    750         const double * solution = solver->getColSolution();
    751         int numberClean = 0;
    752         bool someToDoYet = false;
    753         int numberColumns = solver->getNumCols();
    754         char * mark = new char[numberColumns];
    755         int numberFixed = 0;
    756         for (i = 0; i < numberColumns; i++) {
    757             if (columnLower[i] != columnUpper[i]) {
    758                 mark[i] = 0;
    759             } else {
    760                 mark[i] = 1;
    761                 numberFixed++;
    762             }
    763         }
    764         int numberNewFixed = 0;
    765         for (i = 0; i < numberRows; i++) {
    766             double rhsValue = rowUpper[i];
    767             bool oneRow = true;
    768             // check elements
    769             int numberUnsatisfied = 0;
    770             for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
    771                 int iColumn = column[j];
    772                 double value = elementByRow[j];
    773                 double solValue = solution[iColumn];
    774                 if (columnLower[iColumn] != columnUpper[iColumn]) {
    775                     if (solValue < 1.0 - integerTolerance && solValue > integerTolerance)
    776                         numberUnsatisfied++;
    777                     if (value != 1.0) {
    778                         oneRow = false;
    779                         break;
    780                     }
    781                 } else {
    782                     rhsValue -= value * floor(solValue + 0.5);
    783                 }
    784             }
    785             if (oneRow && rhsValue <= 1.0 + tolerance) {
    786                 if (numberUnsatisfied) {
    787                     someToDoYet = true;
    788                 } else {
    789                     numberClean++;
    790                     for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
    791                         int iColumn = column[j];
    792                         if (columnLower[iColumn] != columnUpper[iColumn] && !mark[iColumn]) {
    793                             mark[iColumn] = 1;
    794                             numberNewFixed++;
    795                         }
    796                     }
    797                 }
    798             }
    799         }
    800         delete [] mark;
    801         //printf("%d clean, %d old fixed, %d new fixed\n",
    802         //   numberClean,numberFixed,numberNewFixed);
    803         if (someToDoYet && numberClean < numberClean_
    804                 && numberNewFixed + numberFixed < wantedFixed) {
    805         } else if (numberFixed < wantedFixed) {
    806             returnCode |= 2;
    807         } else {
    808         }
    809     }
    810     return returnCode;
    811 }
    812 double
    813 CbcBranchToFixLots::infeasibility(const OsiBranchingInformation * /*info*/,
    814                                   int &preferredWay) const
    815 {
    816     preferredWay = -1;
    817     CbcNode * node = model_->currentNode();
    818     int depth;
    819     if (node)
    820         depth = CoinMax(node->depth(), 0);
    821     else
    822         return 0.0;
    823     if (depth_ < 0) {
    824         return 0.0;
    825     } else if (depth_ > 0) {
    826         if ((depth % depth_) != 0)
    827             return 0.0;
    828     }
    829     if (djTolerance_ != -1.234567) {
    830         if (!shallWe())
    831             return 0.0;
    832         else
    833             return 1.0e20;
    834     } else {
    835         // See if 3 in same row and sum <FIX_IF_LESS?
    836         int numberRows = matrixByRow_.getNumRows();
    837         const double * solution = model_->testSolution();
    838         const int * column = matrixByRow_.getIndices();
    839         const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
    840         const int * rowLength = matrixByRow_.getVectorLengths();
    841         double bestSum = 1.0;
    842         int nBest = -1;
    843         OsiSolverInterface * solver = model_->solver();
    844         for (int i = 0; i < numberRows; i++) {
    845             int numberUnsatisfied = 0;
    846             double sum = 0.0;
    847             for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
    848                 int iColumn = column[j];
    849                 if (solver->isInteger(iColumn)) {
    850                     double solValue = solution[iColumn];
    851                     if (solValue > 1.0e-5 && solValue < FIX_IF_LESS) {
    852                         numberUnsatisfied++;
    853                         sum += solValue;
    854                     }
    855                 }
    856             }
    857             if (numberUnsatisfied >= 3 && sum < FIX_IF_LESS) {
    858                 // possible
    859                 if (numberUnsatisfied > nBest ||
    860                         (numberUnsatisfied == nBest && sum < bestSum)) {
    861                     nBest = numberUnsatisfied;
    862                     bestSum = sum;
    863                 }
    864             }
    865         }
    866         if (nBest > 0)
    867             return 1.0e20;
    868         else
    869             return 0.0;
    870     }
    871 }
    872 // Redoes data when sequence numbers change
    873 void
    874 CbcBranchToFixLots::redoSequenceEtc(CbcModel * model, int numberColumns, const int * originalColumns)
    875 {
    876     model_ = model;
    877     if (mark_) {
    878         OsiSolverInterface * solver = model_->solver();
    879         int numberColumnsNow = solver->getNumCols();
    880         char * temp = new char[numberColumnsNow];
    881         memset(temp, 0, numberColumnsNow);
    882         for (int i = 0; i < numberColumns; i++) {
    883             int j = originalColumns[i];
    884             temp[i] = mark_[j];
    885         }
    886         delete [] mark_;
    887         mark_ = temp;
    888     }
    889     OsiSolverInterface * solver = model_->solver();
    890     matrixByRow_ = *solver->getMatrixByRow();
    891 }
    892 
    893 /** Default Constructor
    894 */
    895 CbcBranchAllDifferent::CbcBranchAllDifferent ()
    896         : CbcBranchCut(),
    897         numberInSet_(0),
    898         which_(NULL)
    899 {
    900 }
    901 
    902 /* Useful constructor - passed set of variables
    903 */
    904 CbcBranchAllDifferent::CbcBranchAllDifferent (CbcModel * model, int numberInSet,
    905         const int * members)
    906         : CbcBranchCut(model)
    907 {
    908     numberInSet_ = numberInSet;
    909     which_ = CoinCopyOfArray(members, numberInSet_);
    910 }
    911 // Copy constructor
    912 CbcBranchAllDifferent::CbcBranchAllDifferent ( const CbcBranchAllDifferent & rhs)
    913         : CbcBranchCut(rhs)
    914 {
    915     numberInSet_ = rhs.numberInSet_;
    916     which_ = CoinCopyOfArray(rhs.which_, numberInSet_);
    917 }
    918 
    919 // Clone
    920 CbcObject *
    921 CbcBranchAllDifferent::clone() const
    922 {
    923     return new CbcBranchAllDifferent(*this);
    924 }
    925 
    926 // Assignment operator
    927 CbcBranchAllDifferent &
    928 CbcBranchAllDifferent::operator=( const CbcBranchAllDifferent & rhs)
    929 {
    930     if (this != &rhs) {
    931         CbcBranchCut::operator=(rhs);
    932         delete [] which_;
    933         numberInSet_ = rhs.numberInSet_;
    934         which_ = CoinCopyOfArray(rhs.which_, numberInSet_);
    935     }
    936     return *this;
    937 }
    938 
    939 // Destructor
    940 CbcBranchAllDifferent::~CbcBranchAllDifferent ()
    941 {
    942     delete [] which_;
    943 }
    944 CbcBranchingObject *
    945 CbcBranchAllDifferent::createCbcBranch(OsiSolverInterface * /*solver*/
    946                                        , const OsiBranchingInformation * /*info*/,
    947                                        int /*way*/)
    948 {
    949     // by default way must be -1
    950     //assert (way==-1);
    951     const double * solution = model_->testSolution();
    952     double * values = new double[numberInSet_];
    953     int * which = new int[numberInSet_];
    954     int i;
    955     for (i = 0; i < numberInSet_; i++) {
    956         int iColumn = which_[i];
    957         values[i] = solution[iColumn];
    958         which[i] = iColumn;
    959     }
    960     CoinSort_2(values, values + numberInSet_, which);
    961     double last = -1.0;
    962     double closest = 1.0;
    963     int worst = -1;
    964     for (i = 0; i < numberInSet_; i++) {
    965         if (values[i] - last < closest) {
    966             closest = values[i] - last;
    967             worst = i - 1;
    968         }
    969         last = values[i];
    970     }
    971     assert (closest <= 0.99999);
    972     OsiRowCut down;
    973     down.setLb(-COIN_DBL_MAX);
    974     down.setUb(-1.0);
    975     int pair[2];
    976     double elements[] = {1.0, -1.0};
    977     pair[0] = which[worst];
    978     pair[1] = which[worst+1];
    979     delete [] values;
    980     delete [] which;
    981     down.setRow(2, pair, elements);
    982     // up is same - just with rhs changed
    983     OsiRowCut up = down;
    984     up.setLb(1.0);
    985     up.setUb(COIN_DBL_MAX);
    986     // Say is not a fix type branch
    987     CbcCutBranchingObject * newObject =
    988         new CbcCutBranchingObject(model_, down, up, false);
    989     if (model_->messageHandler()->logLevel() > 1)
    990         printf("creating cut in CbcBranchCut\n");
    991     return newObject;
    992 }
    993 double
    994 CbcBranchAllDifferent::infeasibility(const OsiBranchingInformation * /*info*/,
    995                                      int &preferredWay) const
    996 {
    997     preferredWay = -1;
    998     //OsiSolverInterface * solver = model_->solver();
    999     const double * solution = model_->testSolution();
    1000     //const double * lower = solver->getColLower();
    1001     //const double * upper = solver->getColUpper();
    1002     double * values = new double[numberInSet_];
    1003     int i;
    1004     for (i = 0; i < numberInSet_; i++) {
    1005         int iColumn = which_[i];
    1006         values[i] = solution[iColumn];
    1007     }
    1008     std::sort(values, values + numberInSet_);
    1009     double last = -1.0;
    1010     double closest = 1.0;
    1011     for (i = 0; i < numberInSet_; i++) {
    1012         if (values[i] - last < closest) {
    1013             closest = values[i] - last;
    1014         }
    1015         last = values[i];
    1016     }
    1017     delete [] values;
    1018     if (closest > 0.99999)
    1019         return 0.0;
    1020     else
    1021         return 0.5*(1.0 - closest);
    1022 }
  • branches/sandbox/Cbc/src/CbcBranchCut.hpp

    r1286 r1305  
    88#include "OsiRowCut.hpp"
    99#include "CoinPackedMatrix.hpp"
     10#include "CbcCutBranchingObject.hpp"
    1011
    1112/** Define a cut branching class.
     
    99100
    100101};
    101 
    102 
    103 /** Cut branching object
    104 
    105   This object can specify a two-way branch in terms of two cuts
    106 */
    107 
    108 class CbcCutBranchingObject : public CbcBranchingObject {
    109 
    110 public:
    111 
    112     /// Default constructor
    113     CbcCutBranchingObject ();
    114 
    115     /** Create a cut branching object
    116 
    117         Cut down will applied on way=-1, up on way==1
    118         Assumed down will be first so way_ set to -1
    119     */
    120     CbcCutBranchingObject (CbcModel * model, OsiRowCut & down, OsiRowCut &up, bool canFix);
    121 
    122     /// Copy constructor
    123     CbcCutBranchingObject ( const CbcCutBranchingObject &);
    124 
    125     /// Assignment operator
    126     CbcCutBranchingObject & operator= (const CbcCutBranchingObject& rhs);
    127 
    128     /// Clone
    129     virtual CbcBranchingObject * clone() const;
    130 
    131     /// Destructor
    132     virtual ~CbcCutBranchingObject ();
    133 
    134     using CbcBranchingObject::branch ;
    135     /** \brief Sets the bounds for variables or adds a cut depending on the
    136                current arm of the branch and advances the object state to the next arm.
    137            Returns change in guessed objective on next branch
    138     */
    139     virtual double branch();
    140 
    141 #if 0
    142     // No need to override. Default works fine.
    143     /** Reset every information so that the branching object appears to point to
    144         the previous child. This method does not need to modify anything in any
    145         solver. */
    146     virtual void previousBranch();
    147102#endif
    148 
    149     using CbcBranchingObject::print ;
    150     /** \brief Print something about branch - only if log level high
    151     */
    152     virtual void print();
    153 
    154     /** \brief Return true if branch should fix variables
    155     */
    156     virtual bool boundBranch() const;
    157 
    158     /** Return the type (an integer identifier) of \c this */
    159     virtual int type() const {
    160         return 200;
    161     }
    162 
    163     /** Compare the original object of \c this with the original object of \c
    164         brObj. Assumes that there is an ordering of the original objects.
    165         This method should be invoked only if \c this and brObj are of the same
    166         type.
    167         Return negative/0/positive depending on whether \c this is
    168         smaller/same/larger than the argument.
    169     */
    170     virtual int compareOriginalObject(const CbcBranchingObject* brObj) const;
    171 
    172     /** Compare the \c this with \c brObj. \c this and \c brObj must be os the
    173         same type and must have the same original object, but they may have
    174         different feasible regions.
    175         Return the appropriate CbcRangeCompare value (first argument being the
    176         sub/superset if that's the case). In case of overlap (and if \c
    177         replaceIfOverlap is true) replace the current branching object with one
    178         whose feasible region is the overlap.
    179      */
    180     virtual CbcRangeCompare compareBranchingObject
    181     (const CbcBranchingObject* brObj, const bool replaceIfOverlap = false);
    182 
    183 protected:
    184     /// Cut for the down arm (way_ = -1)
    185     OsiRowCut down_;
    186     /// Cut for the up arm (way_ = 1)
    187     OsiRowCut up_;
    188     /// True if one way can fix variables
    189     bool canFix_;
    190 };
    191 
    192 
    193 /** Define a branch class that branches so that one way variables are fixed
    194     while the other way cuts off that solution.
    195     a) On reduced cost
    196     b) When enough ==1 or <=1 rows have been satisfied (not fixed - satisfied)
    197 */
    198 
    199 
    200 class CbcBranchToFixLots : public CbcBranchCut {
    201 
    202 public:
    203 
    204     // Default Constructor
    205     CbcBranchToFixLots ();
    206 
    207     /** Useful constructor - passed reduced cost tolerance and fraction we would like fixed.
    208         Also depth level to do at.
    209         Also passed number of 1 rows which when clean triggers fix
    210         Always does if all 1 rows cleaned up and number>0 or if fraction columns reached
    211         Also whether to create branch if can't reach fraction.
    212     */
    213     CbcBranchToFixLots (CbcModel * model, double djTolerance,
    214                         double fractionFixed, int depth,
    215                         int numberClean = 0,
    216                         const char * mark = NULL,
    217                         bool alwaysCreate = false);
    218 
    219     // Copy constructor
    220     CbcBranchToFixLots ( const CbcBranchToFixLots &);
    221 
    222     /// Clone
    223     virtual CbcObject * clone() const;
    224 
    225     // Assignment operator
    226     CbcBranchToFixLots & operator=( const CbcBranchToFixLots& rhs);
    227 
    228     // Destructor
    229     ~CbcBranchToFixLots ();
    230 
    231     /** Does a lot of the work,
    232         Returns 0 if no good, 1 if dj, 2 if clean, 3 if both
    233     */
    234     int shallWe() const;
    235 
    236     /// Infeasibility - large is 0.5
    237     virtual double infeasibility(const OsiBranchingInformation * info,
    238                                  int &preferredWay) const;
    239     /** \brief Return true if object can take part in normal heuristics
    240     */
    241     virtual bool canDoHeuristics() const {
    242         return true;
    243     }
    244 
    245     /// Creates a branching object
    246     virtual CbcBranchingObject * createCbcBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) ;
    247     /// Redoes data when sequence numbers change
    248     virtual void redoSequenceEtc(CbcModel * model, int numberColumns, const int * originalColumns);
    249 
    250 
    251 protected:
    252     /// data
    253 
    254     /// Reduced cost tolerance i.e. dj has to be >= this before fixed
    255     double djTolerance_;
    256     /// We only need to make sure this fraction fixed
    257     double fractionFixed_;
    258     /// Never fix ones marked here
    259     char * mark_;
    260     /// Matrix by row
    261     CoinPackedMatrix matrixByRow_;
    262     /// Do if depth multiple of this
    263     int depth_;
    264     /// number of ==1 rows which need to be clean
    265     int numberClean_;
    266     /// If true then always create branch
    267     bool alwaysCreate_;
    268 };
    269 
    270 /** Define a branch class that branches so that it is only satsified if all
    271     members have different values
    272     So cut is x <= y-1 or x >= y+1
    273 */
    274 
    275 
    276 class CbcBranchAllDifferent : public CbcBranchCut {
    277 
    278 public:
    279 
    280     // Default Constructor
    281     CbcBranchAllDifferent ();
    282 
    283     /** Useful constructor - passed set of integer variables which must all be different
    284     */
    285     CbcBranchAllDifferent (CbcModel * model, int number, const int * which);
    286 
    287     // Copy constructor
    288     CbcBranchAllDifferent ( const CbcBranchAllDifferent &);
    289 
    290     /// Clone
    291     virtual CbcObject * clone() const;
    292 
    293     // Assignment operator
    294     CbcBranchAllDifferent & operator=( const CbcBranchAllDifferent& rhs);
    295 
    296     // Destructor
    297     ~CbcBranchAllDifferent ();
    298 
    299     /// Infeasibility - large is 0.5
    300     virtual double infeasibility(const OsiBranchingInformation * info,
    301                                  int &preferredWay) const;
    302 
    303     /// Creates a branching object
    304     virtual CbcBranchingObject * createCbcBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) ;
    305 
    306 
    307 protected:
    308     /// data
    309 
    310     /// Number of entries
    311     int numberInSet_;
    312     /// Which variables
    313     int * which_;
    314 };
    315 #endif
Note: See TracChangeset for help on using the changeset viewer.