source: trunk/Cbc/src/CbcFollowOn.cpp @ 1422

Last change on this file since 1422 was 1393, checked in by lou, 10 years ago

Mark #if 0 with JJF_ZERO and #if 1 with JJF_ONE as a historical reference
point.

File size: 16.7 KB
Line 
1// Edwin 11/10/2009-- carved out of CbcBranchActual
2#if defined(_MSC_VER)
3// Turn off compiler warning about long names
4#  pragma warning(disable:4786)
5#endif
6#include <cassert>
7#include <cstdlib>
8#include <cmath>
9#include <cfloat>
10//#define CBC_DEBUG
11
12#include "CoinTypes.hpp"
13#include "OsiSolverInterface.hpp"
14#include "OsiSolverBranch.hpp"
15#include "CbcModel.hpp"
16#include "CbcMessage.hpp"
17#include "CbcFollowOn.hpp"
18#include "CbcBranchActual.hpp"
19#include "CoinSort.hpp"
20#include "CoinError.hpp"
21
22//##############################################################################
23
24// Default Constructor
25CbcFollowOn::CbcFollowOn ()
26        : CbcObject(),
27        rhs_(NULL)
28{
29}
30
31// Useful constructor
32CbcFollowOn::CbcFollowOn (CbcModel * model)
33        : CbcObject(model)
34{
35    assert (model);
36    OsiSolverInterface * solver = model_->solver();
37    matrix_ = *solver->getMatrixByCol();
38    matrix_.removeGaps();
39    matrix_.setExtraGap(0.0);
40    matrixByRow_ = *solver->getMatrixByRow();
41    int numberRows = matrix_.getNumRows();
42
43    rhs_ = new int[numberRows];
44    int i;
45    const double * rowLower = solver->getRowLower();
46    const double * rowUpper = solver->getRowUpper();
47    // Row copy
48    const double * elementByRow = matrixByRow_.getElements();
49    const int * column = matrixByRow_.getIndices();
50    const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
51    const int * rowLength = matrixByRow_.getVectorLengths();
52    for (i = 0; i < numberRows; i++) {
53        rhs_[i] = 0;
54        double value = rowLower[i];
55        if (value == rowUpper[i]) {
56            if (floor(value) == value && value >= 1.0 && value < 10.0) {
57                // check elements
58                bool good = true;
59                for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
60                    int iColumn = column[j];
61                    if (!solver->isBinary(iColumn))
62                        good = false;
63                    double elValue = elementByRow[j];
64                    if (floor(elValue) != elValue || value < 1.0)
65                        good = false;
66                }
67                if (good)
68                    rhs_[i] = static_cast<int> (value);
69            }
70        }
71    }
72}
73
74// Copy constructor
75CbcFollowOn::CbcFollowOn ( const CbcFollowOn & rhs)
76        : CbcObject(rhs),
77        matrix_(rhs.matrix_),
78        matrixByRow_(rhs.matrixByRow_)
79{
80    int numberRows = matrix_.getNumRows();
81    rhs_ = CoinCopyOfArray(rhs.rhs_, numberRows);
82}
83
84// Clone
85CbcObject *
86CbcFollowOn::clone() const
87{
88    return new CbcFollowOn(*this);
89}
90
91// Assignment operator
92CbcFollowOn &
93CbcFollowOn::operator=( const CbcFollowOn & rhs)
94{
95    if (this != &rhs) {
96        CbcObject::operator=(rhs);
97        delete [] rhs_;
98        matrix_ = rhs.matrix_;
99        matrixByRow_ = rhs.matrixByRow_;
100        int numberRows = matrix_.getNumRows();
101        rhs_ = CoinCopyOfArray(rhs.rhs_, numberRows);
102    }
103    return *this;
104}
105
106// Destructor
107CbcFollowOn::~CbcFollowOn ()
108{
109    delete [] rhs_;
110}
111// As some computation is needed in more than one place - returns row
112int
113CbcFollowOn::gutsOfFollowOn(int & otherRow, int & preferredWay) const
114{
115    int whichRow = -1;
116    otherRow = -1;
117    int numberRows = matrix_.getNumRows();
118
119    int i;
120    // For sorting
121    int * sort = new int [numberRows];
122    int * isort = new int [numberRows];
123    // Column copy
124    //const double * element = matrix_.getElements();
125    const int * row = matrix_.getIndices();
126    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
127    const int * columnLength = matrix_.getVectorLengths();
128    // Row copy
129    const double * elementByRow = matrixByRow_.getElements();
130    const int * column = matrixByRow_.getIndices();
131    const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
132    const int * rowLength = matrixByRow_.getVectorLengths();
133    OsiSolverInterface * solver = model_->solver();
134    const double * columnLower = solver->getColLower();
135    const double * columnUpper = solver->getColUpper();
136    const double * solution = solver->getColSolution();
137    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
138    int nSort = 0;
139    for (i = 0; i < numberRows; i++) {
140        if (rhs_[i]) {
141            // check elements
142            double smallest = 1.0e10;
143            double largest = 0.0;
144            int rhsValue = rhs_[i];
145            int number1 = 0;
146            int numberUnsatisfied = 0;
147            for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
148                int iColumn = column[j];
149                double value = elementByRow[j];
150                double solValue = solution[iColumn];
151                if (columnLower[iColumn] != columnUpper[iColumn]) {
152                    smallest = CoinMin(smallest, value);
153                    largest = CoinMax(largest, value);
154                    if (value == 1.0)
155                        number1++;
156                    if (solValue < 1.0 - integerTolerance && solValue > integerTolerance)
157                        numberUnsatisfied++;
158                } else {
159                    rhsValue -= static_cast<int>(value * floor(solValue + 0.5));
160                }
161            }
162            if (numberUnsatisfied > 1) {
163                if (smallest < largest) {
164                    // probably no good but check a few things
165                    assert (largest <= rhsValue);
166                    if (number1 == 1 && largest == rhsValue)
167                        printf("could fix\n");
168                } else if (largest == rhsValue) {
169                    sort[nSort] = i;
170                    isort[nSort++] = -numberUnsatisfied;
171                }
172            }
173        }
174    }
175    if (nSort > 1) {
176        CoinSort_2(isort, isort + nSort, sort);
177        CoinZeroN(isort, numberRows);
178        double * other = new double[numberRows];
179        CoinZeroN(other, numberRows);
180        int * which = new int[numberRows];
181        //#define COUNT
182#ifndef COUNT
183        bool beforeSolution = model_->getSolutionCount() == 0;
184#endif
185        for (int k = 0; k < nSort - 1; k++) {
186            i = sort[k];
187            int numberUnsatisfied = 0;
188            int n = 0;
189            int j;
190            for (j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
191                int iColumn = column[j];
192                if (columnLower[iColumn] != columnUpper[iColumn]) {
193                    double solValue = solution[iColumn] - columnLower[iColumn];
194                    if (solValue < 1.0 - integerTolerance && solValue > integerTolerance) {
195                        numberUnsatisfied++;
196                        for (int jj = columnStart[iColumn]; jj < columnStart[iColumn] + columnLength[iColumn]; jj++) {
197                            int iRow = row[jj];
198                            if (rhs_[iRow]) {
199                                other[iRow] += solValue;
200                                if (isort[iRow]) {
201                                    isort[iRow]++;
202                                } else {
203                                    isort[iRow] = 1;
204                                    which[n++] = iRow;
205                                }
206                            }
207                        }
208                    }
209                }
210            }
211            double total = 0.0;
212            // Take out row
213            double sumThis = other[i];
214            other[i] = 0.0;
215            assert (numberUnsatisfied == isort[i]);
216            // find one nearest half if solution, one if before solution
217            int iBest = -1;
218            double dtarget = 0.5 * total;
219#ifdef COUNT
220            int target = (numberUnsatisfied + 1) >> 1;
221            int best = numberUnsatisfied;
222#else
223            double best;
224            if (beforeSolution)
225                best = dtarget;
226            else
227                best = 1.0e30;
228#endif
229            for (j = 0; j < n; j++) {
230                int iRow = which[j];
231                double dvalue = other[iRow];
232                other[iRow] = 0.0;
233#ifdef COUNT
234                int value = isort[iRow];
235#endif
236                isort[iRow] = 0;
237                if (fabs(dvalue) < 1.0e-8 || fabs(sumThis - dvalue) < 1.0e-8)
238                    continue;
239                if (dvalue < integerTolerance || dvalue > 1.0 - integerTolerance)
240                    continue;
241#ifdef COUNT
242                if (abs(value - target) < best && value != numberUnsatisfied) {
243                    best = abs(value - target);
244                    iBest = iRow;
245                    if (dvalue < dtarget)
246                        preferredWay = 1;
247                    else
248                        preferredWay = -1;
249                }
250#else
251                if (beforeSolution) {
252                    if (fabs(dvalue - dtarget) > best) {
253                        best = fabs(dvalue - dtarget);
254                        iBest = iRow;
255                        if (dvalue < dtarget)
256                            preferredWay = 1;
257                        else
258                            preferredWay = -1;
259                    }
260                } else {
261                    if (fabs(dvalue - dtarget) < best) {
262                        best = fabs(dvalue - dtarget);
263                        iBest = iRow;
264                        if (dvalue < dtarget)
265                            preferredWay = 1;
266                        else
267                            preferredWay = -1;
268                    }
269                }
270#endif
271            }
272            if (iBest >= 0) {
273                whichRow = i;
274                otherRow = iBest;
275                break;
276            }
277        }
278        delete [] which;
279        delete [] other;
280    }
281    delete [] sort;
282    delete [] isort;
283    return whichRow;
284}
285double
286CbcFollowOn::infeasibility(const OsiBranchingInformation * /*info*/,
287                           int &preferredWay) const
288{
289    int otherRow = 0;
290    int whichRow = gutsOfFollowOn(otherRow, preferredWay);
291    if (whichRow < 0)
292        return 0.0;
293    else
294        return 2.0* model_->getDblParam(CbcModel::CbcIntegerTolerance);
295}
296
297// This looks at solution and sets bounds to contain solution
298void
299CbcFollowOn::feasibleRegion()
300{
301}
302
303CbcBranchingObject *
304CbcFollowOn::createCbcBranch(OsiSolverInterface * solver, const OsiBranchingInformation * /*info*/, int way)
305{
306    int otherRow = 0;
307    int preferredWay;
308    int whichRow = gutsOfFollowOn(otherRow, preferredWay);
309    assert(way == preferredWay);
310    assert (whichRow >= 0);
311    int numberColumns = matrix_.getNumCols();
312
313    // Column copy
314    //const double * element = matrix_.getElements();
315    const int * row = matrix_.getIndices();
316    const CoinBigIndex * columnStart = matrix_.getVectorStarts();
317    const int * columnLength = matrix_.getVectorLengths();
318    // Row copy
319    //const double * elementByRow = matrixByRow_.getElements();
320    const int * column = matrixByRow_.getIndices();
321    const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
322    const int * rowLength = matrixByRow_.getVectorLengths();
323    //OsiSolverInterface * solver = model_->solver();
324    const double * columnLower = solver->getColLower();
325    const double * columnUpper = solver->getColUpper();
326    //const double * solution = solver->getColSolution();
327    int nUp = 0;
328    int nDown = 0;
329    int * upList = new int[numberColumns];
330    int * downList = new int[numberColumns];
331    int j;
332    for (j = rowStart[whichRow]; j < rowStart[whichRow] + rowLength[whichRow]; j++) {
333        int iColumn = column[j];
334        if (columnLower[iColumn] != columnUpper[iColumn]) {
335            bool up = true;
336            for (int jj = columnStart[iColumn]; jj < columnStart[iColumn] + columnLength[iColumn]; jj++) {
337                int iRow = row[jj];
338                if (iRow == otherRow) {
339                    up = false;
340                    break;
341                }
342            }
343            if (up)
344                upList[nUp++] = iColumn;
345            else
346                downList[nDown++] = iColumn;
347        }
348    }
349    //printf("way %d\n",way);
350    // create object
351    //printf("would fix %d down and %d up\n",nDown,nUp);
352    CbcBranchingObject * branch
353    = new CbcFixingBranchingObject(model_, way,
354                                   nDown, downList, nUp, upList);
355    delete [] upList;
356    delete [] downList;
357    return branch;
358}
359
360//##############################################################################
361
362// Default Constructor
363CbcFixingBranchingObject::CbcFixingBranchingObject()
364        : CbcBranchingObject()
365{
366    numberDown_ = 0;
367    numberUp_ = 0;
368    downList_ = NULL;
369    upList_ = NULL;
370}
371
372// Useful constructor
373CbcFixingBranchingObject::CbcFixingBranchingObject (CbcModel * model,
374        int way ,
375        int numberOnDownSide, const int * down,
376        int numberOnUpSide, const int * up)
377        : CbcBranchingObject(model, 0, way, 0.5)
378{
379    numberDown_ = numberOnDownSide;
380    numberUp_ = numberOnUpSide;
381    downList_ = CoinCopyOfArray(down, numberDown_);
382    upList_ = CoinCopyOfArray(up, numberUp_);
383}
384
385// Copy constructor
386CbcFixingBranchingObject::CbcFixingBranchingObject ( const CbcFixingBranchingObject & rhs) : CbcBranchingObject(rhs)
387{
388    numberDown_ = rhs.numberDown_;
389    numberUp_ = rhs.numberUp_;
390    downList_ = CoinCopyOfArray(rhs.downList_, numberDown_);
391    upList_ = CoinCopyOfArray(rhs.upList_, numberUp_);
392}
393
394// Assignment operator
395CbcFixingBranchingObject &
396CbcFixingBranchingObject::operator=( const CbcFixingBranchingObject & rhs)
397{
398    if (this != &rhs) {
399        CbcBranchingObject::operator=(rhs);
400        delete [] downList_;
401        delete [] upList_;
402        numberDown_ = rhs.numberDown_;
403        numberUp_ = rhs.numberUp_;
404        downList_ = CoinCopyOfArray(rhs.downList_, numberDown_);
405        upList_ = CoinCopyOfArray(rhs.upList_, numberUp_);
406    }
407    return *this;
408}
409CbcBranchingObject *
410CbcFixingBranchingObject::clone() const
411{
412    return (new CbcFixingBranchingObject(*this));
413}
414
415
416// Destructor
417CbcFixingBranchingObject::~CbcFixingBranchingObject ()
418{
419    delete [] downList_;
420    delete [] upList_;
421}
422double
423CbcFixingBranchingObject::branch()
424{
425    decrementNumberBranchesLeft();
426    OsiSolverInterface * solver = model_->solver();
427    const double * columnLower = solver->getColLower();
428    int i;
429    // *** for way - up means fix all those in up section
430    if (way_ < 0) {
431#ifdef FULL_PRINT
432        printf("Down Fix ");
433#endif
434        //printf("Down Fix %d\n",numberDown_);
435        for (i = 0; i < numberDown_; i++) {
436            int iColumn = downList_[i];
437            model_->solver()->setColUpper(iColumn, columnLower[iColumn]);
438#ifdef FULL_PRINT
439            printf("Setting bound on %d to lower bound\n", iColumn);
440#endif
441        }
442        way_ = 1;         // Swap direction
443    } else {
444#ifdef FULL_PRINT
445        printf("Up Fix ");
446#endif
447        //printf("Up Fix %d\n",numberUp_);
448        for (i = 0; i < numberUp_; i++) {
449            int iColumn = upList_[i];
450            model_->solver()->setColUpper(iColumn, columnLower[iColumn]);
451#ifdef FULL_PRINT
452            printf("Setting bound on %d to lower bound\n", iColumn);
453#endif
454        }
455        way_ = -1;        // Swap direction
456    }
457#ifdef FULL_PRINT
458    printf("\n");
459#endif
460    return 0.0;
461}
462void
463CbcFixingBranchingObject::print()
464{
465    int i;
466    // *** for way - up means fix all those in up section
467    if (way_ < 0) {
468        printf("Down Fix ");
469        for (i = 0; i < numberDown_; i++) {
470            int iColumn = downList_[i];
471            printf("%d ", iColumn);
472        }
473    } else {
474        printf("Up Fix ");
475        for (i = 0; i < numberUp_; i++) {
476            int iColumn = upList_[i];
477            printf("%d ", iColumn);
478        }
479    }
480    printf("\n");
481}
482
483/** Compare the original object of \c this with the original object of \c
484    brObj. Assumes that there is an ordering of the original objects.
485    This method should be invoked only if \c this and brObj are of the same
486    type.
487    Return negative/0/positive depending on whether \c this is
488    smaller/same/larger than the argument.
489*/
490int
491CbcFixingBranchingObject::compareOriginalObject
492(const CbcBranchingObject* /*brObj*/) const
493{
494    throw("must implement");
495}
496
497/** Compare the \c this with \c brObj. \c this and \c brObj must be os the
498    same type and must have the same original object, but they may have
499    different feasible regions.
500    Return the appropriate CbcRangeCompare value (first argument being the
501    sub/superset if that's the case). In case of overlap (and if \c
502    replaceIfOverlap is true) replace the current branching object with one
503    whose feasible region is the overlap.
504   */
505CbcRangeCompare
506CbcFixingBranchingObject::compareBranchingObject
507(const CbcBranchingObject* /*brObj*/, const bool /*replaceIfOverlap*/)
508{
509#ifdef JJF_ZERO //ndef NDEBUG
510    const CbcFixingBranchingObject* br =
511        dynamic_cast<const CbcFixingBranchingObject*>(brObj);
512    assert(br);
513#endif
514    // If two FixingBranchingObject's have the same base object then it's pretty
515    // much guaranteed
516    throw("must implement");
517}
518
519//##############################################################################
Note: See TracBrowser for help on using the repository browser.