source: branches/sandbox/Cbc/src/CbcBranchCut.cpp @ 1286

Last change on this file since 1286 was 1286, checked in by EdwinStraver, 10 years ago

Changed formatting using AStyle -A4 -p

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 32.6 KB
RevLine 
[1271]1/* $Id: CbcBranchCut.cpp 1286 2009-11-09 23:33:07Z EdwinStraver $ */
[5]2// Copyright (C) 2004, International Business Machines
3// Corporation and others.  All Rights Reserved.
4#if defined(_MSC_VER)
5// Turn off compiler warning about long names
6#  pragma warning(disable:4786)
7#endif
8#include <cassert>
[904]9#include <cstdlib>
[5]10#include <cmath>
11#include <cfloat>
12//#define CBC_DEBUG
13
14#include "OsiSolverInterface.hpp"
15#include "CbcModel.hpp"
16#include "CbcMessage.hpp"
17#include "CbcBranchCut.hpp"
18#include "CoinSort.hpp"
19#include "CoinError.hpp"
20
21
22/** Default Constructor
23
24*/
25CbcBranchCut::CbcBranchCut ()
[1286]26        : CbcObject()
[5]27{
28}
29
30/* Constructor so model can be passed up
[1286]31*/
[5]32CbcBranchCut::CbcBranchCut (CbcModel * model)
[1286]33        : CbcObject(model)
[5]34{
35}
[1286]36// Copy constructor
[5]37CbcBranchCut::CbcBranchCut ( const CbcBranchCut & rhs)
[1286]38        : CbcObject(rhs)
[5]39
40{
41}
42
43// Clone
44CbcObject *
45CbcBranchCut::clone() const
46{
[1286]47    return new CbcBranchCut(*this);
[5]48}
49
[1286]50// Assignment operator
51CbcBranchCut &
[1271]52CbcBranchCut::operator=( const CbcBranchCut& /*rhs*/)
[5]53{
[1286]54    return *this;
[5]55}
56
[1286]57// Destructor
[5]58CbcBranchCut::~CbcBranchCut ()
59{
60}
[1286]61double
[1271]62CbcBranchCut::infeasibility(const OsiBranchingInformation * /*info*/,
[1286]63                            int &preferredWay) const
[5]64{
[1286]65    throw CoinError("Use of base class", "infeasibility", "CbcBranchCut");
66    preferredWay = -1;
67    return 0.0;
[5]68}
69
70// This looks at solution and sets bounds to contain solution
71/** More precisely: it first forces the variable within the existing
72    bounds, and then tightens the bounds to fix the variable at the
73    nearest integer value.
74*/
[1286]75void
[5]76CbcBranchCut::feasibleRegion()
77{
78}
[6]79/* Return true if branch created by object should fix variables
80 */
[1286]81bool
82CbcBranchCut::boundBranch() const
[5]83{
[1286]84    return false;
[5]85}
[1286]86CbcBranchingObject *
87CbcBranchCut::createCbcBranch(OsiSolverInterface * /*solver*/, const OsiBranchingInformation * /*info*/, int /*way*/)
88{
89    throw CoinError("Use of base class", "createCbcBranch", "CbcBranchCut");
90    return new CbcCutBranchingObject();
91}
[5]92
93
94/* Given valid solution (i.e. satisfied) and reduced costs etc
95   returns a branching object which would give a new feasible
96   point in direction reduced cost says would be cheaper.
97   If no feasible point returns null
98*/
[1286]99CbcBranchingObject *
[5]100CbcBranchCut::preferredNewFeasible() const
101{
[1286]102    throw CoinError("Use of base class", "preferredNewFeasible", "CbcBranchCut");
103    return new CbcCutBranchingObject();
[5]104}
[1286]105
[5]106/* Given valid solution (i.e. satisfied) and reduced costs etc
107   returns a branching object which would give a new feasible
108   point in direction opposite to one reduced cost says would be cheaper.
109   If no feasible point returns null
110*/
[1286]111CbcBranchingObject *
112CbcBranchCut::notPreferredNewFeasible() const
[5]113{
[1286]114    throw CoinError("Use of base class", "notPreferredNewFeasible", "CbcBranchCut");
115    return new CbcCutBranchingObject();
[5]116}
[1286]117
[5]118/*
119  Bounds may be tightened, so it may be good to be able to refresh the local
120  copy of the original bounds.
121 */
[1286]122void
[5]123CbcBranchCut::resetBounds()
124{
125}
126
127
[1286]128// Default Constructor
[5]129CbcCutBranchingObject::CbcCutBranchingObject()
[1286]130        : CbcBranchingObject()
[5]131{
[1286]132    down_ = OsiRowCut();
133    up_ = OsiRowCut();
134    canFix_ = false;
[5]135}
136
137// Useful constructor
[1286]138CbcCutBranchingObject::CbcCutBranchingObject (CbcModel * model,
139        OsiRowCut & down,
140        OsiRowCut &up,
141        bool canFix)
142        : CbcBranchingObject(model, 0, -1, 0.0)
[5]143{
[1286]144    down_ = down;
145    up_ = up;
146    canFix_ = canFix;
[5]147}
148
[1286]149// Copy constructor
150CbcCutBranchingObject::CbcCutBranchingObject ( const CbcCutBranchingObject & rhs) : CbcBranchingObject(rhs)
[5]151{
152    down_ = rhs.down_;
153    up_ = rhs.up_;
[216]154    canFix_ = rhs.canFix_;
[5]155}
[1286]156
157// Assignment operator
158CbcCutBranchingObject &
159CbcCutBranchingObject::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}
169CbcBranchingObject *
[5]170CbcCutBranchingObject::clone() const
[1286]171{
172    return (new CbcCutBranchingObject(*this));
[5]173}
174
175
[1286]176// Destructor
[5]177CbcCutBranchingObject::~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*/
188double
[640]189CbcCutBranchingObject::branch()
[5]190{
[1286]191    decrementNumberBranchesLeft();
192    OsiRowCut * cut;
193    if (way_ < 0) {
194        cut = &down_;
195        way_ = 1;
[5]196    } else {
[1286]197        cut = &up_;
198        way_ = -1;        // Swap direction
[5]199    }
[1286]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        }
[5]223    }
[1286]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);
[5]253    }
[1286]254    return 0.0;
[5]255}
[1286]256// Print what would happen
[5]257void
[640]258CbcCutBranchingObject::print()
[5]259{
[1286]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");
[5]267    }
[1286]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    }
[5]284}
285
286// Return true if branch should fix variables
[1286]287bool
[5]288CbcCutBranchingObject::boundBranch() const
289{
[1286]290    return false;
[5]291}
292
[912]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
[1286]296    type.
[912]297    Return negative/0/positive depending on whether \c this is
298    smaller/same/larger than the argument.
299*/
300int
301CbcCutBranchingObject::compareOriginalObject
302(const CbcBranchingObject* brObj) const
303{
[1286]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());
[912]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
321CbcRangeCompare
322CbcCutBranchingObject::compareBranchingObject
323(const CbcBranchingObject* brObj, const bool replaceIfOverlap)
324{
[1286]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]);
[912]342    return comp;
343}
344
345//##############################################################################
346
[5]347/** Default Constructor
348
349  Equivalent to an unspecified binary variable.
350*/
[11]351CbcBranchToFixLots::CbcBranchToFixLots ()
[1286]352        : CbcBranchCut(),
353        djTolerance_(COIN_DBL_MAX),
354        fractionFixed_(1.0),
355        mark_(NULL),
356        depth_(-1),
357        numberClean_(0),
358        alwaysCreate_(false)
[5]359{
360}
361
362/* Useful constructor - passed reduced cost tolerance and fraction we would like fixed.
363   Also depth level to do at.
[11]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
[5]366   Also whether to create branch if can't reach fraction.
[1286]367*/
[11]368CbcBranchToFixLots::CbcBranchToFixLots (CbcModel * model, double djTolerance,
[1286]369                                        double fractionFixed, int depth,
370                                        int numberClean,
371                                        const char * mark, bool alwaysCreate)
372        : CbcBranchCut(model)
[5]373{
[1286]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;
[5]389}
[1286]390// Copy constructor
[11]391CbcBranchToFixLots::CbcBranchToFixLots ( const CbcBranchToFixLots & rhs)
[1286]392        : CbcBranchCut(rhs)
[5]393{
[1286]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_;
[5]402}
403
404// Clone
405CbcObject *
[11]406CbcBranchToFixLots::clone() const
[5]407{
[1286]408    return new CbcBranchToFixLots(*this);
[5]409}
410
[1286]411// Assignment operator
412CbcBranchToFixLots &
413CbcBranchToFixLots::operator=( const CbcBranchToFixLots & rhs)
[5]414{
[1286]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;
[5]428}
429
[1286]430// Destructor
[11]431CbcBranchToFixLots::~CbcBranchToFixLots ()
[5]432{
[1286]433    delete [] mark_;
[5]434}
[1286]435CbcBranchingObject *
436CbcBranchToFixLots::createCbcBranch(OsiSolverInterface * solver, const OsiBranchingInformation * /*info*/, int /*way*/)
[5]437{
[1286]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        }
[1271]573    } else {
574#define FIX_IF_LESS -0.1
[1286]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        }
[5]618    }
[1286]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        }
[1121]633    }
[1286]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;
[5]649}
[11]650/* Does a lot of the work,
651   Returns 0 if no good, 1 if dj, 2 if clean, 3 if both
[1121]652   10 if branching on ones away from bound
[11]653*/
[1286]654int
[11]655CbcBranchToFixLots::shallWe() const
[5]656{
[1286]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();
[11]666    int i;
[1286]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;
[11]702    }
[1286]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        }
[11]737    }
[1286]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        }
[5]809    }
[1286]810    return returnCode;
[11]811}
[1286]812double
[1271]813CbcBranchToFixLots::infeasibility(const OsiBranchingInformation * /*info*/,
[1286]814                                  int &preferredWay) const
[11]815{
[1286]816    preferredWay = -1;
817    CbcNode * node = model_->currentNode();
818    int depth;
819    if (node)
820        depth = CoinMax(node->depth(), 0);
[1271]821    else
[1286]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;
[1271]828    }
[1286]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    }
[5]871}
[1271]872// Redoes data when sequence numbers change
[1286]873void
[1271]874CbcBranchToFixLots::redoSequenceEtc(CbcModel * model, int numberColumns, const int * originalColumns)
875{
[1286]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    }
[1271]889    OsiSolverInterface * solver = model_->solver();
[1286]890    matrixByRow_ = *solver->getMatrixByRow();
[1271]891}
[216]892
893/** Default Constructor
894*/
895CbcBranchAllDifferent::CbcBranchAllDifferent ()
[1286]896        : CbcBranchCut(),
897        numberInSet_(0),
898        which_(NULL)
[216]899{
900}
901
902/* Useful constructor - passed set of variables
[1286]903*/
[216]904CbcBranchAllDifferent::CbcBranchAllDifferent (CbcModel * model, int numberInSet,
[1286]905        const int * members)
906        : CbcBranchCut(model)
[216]907{
[1286]908    numberInSet_ = numberInSet;
909    which_ = CoinCopyOfArray(members, numberInSet_);
[216]910}
[1286]911// Copy constructor
[216]912CbcBranchAllDifferent::CbcBranchAllDifferent ( const CbcBranchAllDifferent & rhs)
[1286]913        : CbcBranchCut(rhs)
[216]914{
[1286]915    numberInSet_ = rhs.numberInSet_;
916    which_ = CoinCopyOfArray(rhs.which_, numberInSet_);
[216]917}
918
919// Clone
920CbcObject *
921CbcBranchAllDifferent::clone() const
922{
[1286]923    return new CbcBranchAllDifferent(*this);
[216]924}
925
[1286]926// Assignment operator
927CbcBranchAllDifferent &
928CbcBranchAllDifferent::operator=( const CbcBranchAllDifferent & rhs)
[216]929{
[1286]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;
[216]937}
938
[1286]939// Destructor
[216]940CbcBranchAllDifferent::~CbcBranchAllDifferent ()
941{
[1286]942    delete [] which_;
[216]943}
[1286]944CbcBranchingObject *
[1271]945CbcBranchAllDifferent::createCbcBranch(OsiSolverInterface * /*solver*/
[1286]946                                       , const OsiBranchingInformation * /*info*/,
947                                       int /*way*/)
[216]948{
[1286]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;
[216]959    }
[1286]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;
[216]992}
[1286]993double
[1271]994CbcBranchAllDifferent::infeasibility(const OsiBranchingInformation * /*info*/,
[1286]995                                     int &preferredWay) const
[216]996{
[1286]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];
[216]1007    }
[1286]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);
[216]1022}
Note: See TracBrowser for help on using the repository browser.