source: stable/2.7/Cbc/src/CbcFollowOn.cpp

Last change on this file was 1573, checked in by lou, 9 years ago

Change to EPL license notice.

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