source: branches/sandbox/Cbc/src/CbcFollowOn.cpp @ 1293

Last change on this file since 1293 was 1293, checked in by EdwinStraver, 10 years ago
File size: 12.3 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}
Note: See TracBrowser for help on using the repository browser.