source: trunk/Cbc/src/CbcBranchToFixLots.cpp @ 1514

Last change on this file since 1514 was 1432, checked in by bjarni, 10 years ago

Added extra return at end of each source file where needed, to remove possible linefeed conflicts (NightlyBuild? errors)

File size: 20.4 KB
Line 
1// Edwin 11/13/2009-- carved out of CbcBranchCut
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 "OsiSolverInterface.hpp"
13#include "CbcModel.hpp"
14#include "CbcMessage.hpp"
15#include "CbcBranchCut.hpp"
16#include "CoinSort.hpp"
17#include "CoinError.hpp"
18#include "CbcBranchToFixLots.hpp"
19
20/** Default Constructor
21
22  Equivalent to an unspecified binary variable.
23*/
24CbcBranchToFixLots::CbcBranchToFixLots ()
25        : CbcBranchCut(),
26        djTolerance_(COIN_DBL_MAX),
27        fractionFixed_(1.0),
28        mark_(NULL),
29        depth_(-1),
30        numberClean_(0),
31        alwaysCreate_(false)
32{
33}
34
35/* Useful constructor - passed reduced cost tolerance and fraction we would like fixed.
36   Also depth level to do at.
37   Also passed number of 1 rows which when clean triggers fix
38   Always does if all 1 rows cleaned up and number>0 or if fraction columns reached
39   Also whether to create branch if can't reach fraction.
40*/
41CbcBranchToFixLots::CbcBranchToFixLots (CbcModel * model, double djTolerance,
42                                        double fractionFixed, int depth,
43                                        int numberClean,
44                                        const char * mark, bool alwaysCreate)
45        : CbcBranchCut(model)
46{
47    djTolerance_ = djTolerance;
48    fractionFixed_ = fractionFixed;
49    if (mark) {
50        int numberColumns = model->getNumCols();
51        mark_ = new char[numberColumns];
52        memcpy(mark_, mark, numberColumns);
53    } else {
54        mark_ = NULL;
55    }
56    depth_ = depth;
57    assert (model);
58    OsiSolverInterface * solver = model_->solver();
59    matrixByRow_ = *solver->getMatrixByRow();
60    numberClean_ = numberClean;
61    alwaysCreate_ = alwaysCreate;
62}
63// Copy constructor
64CbcBranchToFixLots::CbcBranchToFixLots ( const CbcBranchToFixLots & rhs)
65        : CbcBranchCut(rhs)
66{
67    djTolerance_ = rhs.djTolerance_;
68    fractionFixed_ = rhs.fractionFixed_;
69    int numberColumns = model_->getNumCols();
70    mark_ = CoinCopyOfArray(rhs.mark_, numberColumns);
71    matrixByRow_ = rhs.matrixByRow_;
72    depth_ = rhs.depth_;
73    numberClean_ = rhs.numberClean_;
74    alwaysCreate_ = rhs.alwaysCreate_;
75}
76
77// Clone
78CbcObject *
79CbcBranchToFixLots::clone() const
80{
81    return new CbcBranchToFixLots(*this);
82}
83
84// Assignment operator
85CbcBranchToFixLots &
86CbcBranchToFixLots::operator=( const CbcBranchToFixLots & rhs)
87{
88    if (this != &rhs) {
89        CbcBranchCut::operator=(rhs);
90        djTolerance_ = rhs.djTolerance_;
91        fractionFixed_ = rhs.fractionFixed_;
92        int numberColumns = model_->getNumCols();
93        delete [] mark_;
94        mark_ = CoinCopyOfArray(rhs.mark_, numberColumns);
95        matrixByRow_ = rhs.matrixByRow_;
96        depth_ = rhs.depth_;
97        numberClean_ = rhs.numberClean_;
98        alwaysCreate_ = rhs.alwaysCreate_;
99    }
100    return *this;
101}
102
103// Destructor
104CbcBranchToFixLots::~CbcBranchToFixLots ()
105{
106    delete [] mark_;
107}
108CbcBranchingObject *
109CbcBranchToFixLots::createCbcBranch(OsiSolverInterface * solver, const OsiBranchingInformation * /*info*/, int /*way*/)
110{
111    // by default way must be -1
112    //assert (way==-1);
113    //OsiSolverInterface * solver = model_->solver();
114    const double * solution = model_->testSolution();
115    const double * lower = solver->getColLower();
116    const double * upper = solver->getColUpper();
117    const double * dj = solver->getReducedCost();
118    int i;
119    int numberIntegers = model_->numberIntegers();
120    const int * integerVariable = model_->integerVariable();
121    double integerTolerance =
122        model_->getDblParam(CbcModel::CbcIntegerTolerance);
123    // make smaller ?
124    double tolerance = CoinMin(1.0e-8, integerTolerance);
125    // How many fixed are we aiming at
126    int wantedFixed = static_cast<int> (static_cast<double>(numberIntegers) * fractionFixed_);
127    int nSort = 0;
128    int numberFixed = 0;
129    int numberColumns = solver->getNumCols();
130    int * sort = new int[numberColumns];
131    double * dsort = new double[numberColumns];
132    if (djTolerance_ != -1.234567) {
133        int type = shallWe();
134        assert (type);
135        // Take clean first
136        if (type == 1) {
137            for (i = 0; i < numberIntegers; i++) {
138                int iColumn = integerVariable[i];
139                if (upper[iColumn] > lower[iColumn]) {
140                    if (!mark_ || !mark_[iColumn]) {
141                        if (solution[iColumn] < lower[iColumn] + tolerance) {
142                            if (dj[iColumn] > djTolerance_) {
143                                dsort[nSort] = -dj[iColumn];
144                                sort[nSort++] = iColumn;
145                            }
146                        } else if (solution[iColumn] > upper[iColumn] - tolerance) {
147                            if (dj[iColumn] < -djTolerance_) {
148                                dsort[nSort] = dj[iColumn];
149                                sort[nSort++] = iColumn;
150                            }
151                        }
152                    }
153                } else {
154                    numberFixed++;
155                }
156            }
157            // sort
158            CoinSort_2(dsort, dsort + nSort, sort);
159            nSort = CoinMin(nSort, wantedFixed - numberFixed);
160        } else if (type < 10) {
161            int i;
162            //const double * rowLower = solver->getRowLower();
163            const double * rowUpper = solver->getRowUpper();
164            // Row copy
165            const double * elementByRow = matrixByRow_.getElements();
166            const int * column = matrixByRow_.getIndices();
167            const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
168            const int * rowLength = matrixByRow_.getVectorLengths();
169            const double * columnLower = solver->getColLower();
170            const double * columnUpper = solver->getColUpper();
171            const double * solution = solver->getColSolution();
172            int numberColumns = solver->getNumCols();
173            int numberRows = solver->getNumRows();
174            for (i = 0; i < numberColumns; i++) {
175                sort[i] = i;
176                if (columnLower[i] != columnUpper[i]) {
177                    dsort[i] = 1.0e100;
178                } else {
179                    dsort[i] = 1.0e50;
180                    numberFixed++;
181                }
182            }
183            for (i = 0; i < numberRows; i++) {
184                double rhsValue = rowUpper[i];
185                bool oneRow = true;
186                // check elements
187                int numberUnsatisfied = 0;
188                for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
189                    int iColumn = column[j];
190                    double value = elementByRow[j];
191                    double solValue = solution[iColumn];
192                    if (columnLower[iColumn] != columnUpper[iColumn]) {
193                        if (solValue < 1.0 - integerTolerance && solValue > integerTolerance)
194                            numberUnsatisfied++;
195                        if (value != 1.0) {
196                            oneRow = false;
197                            break;
198                        }
199                    } else {
200                        rhsValue -= value * floor(solValue + 0.5);
201                    }
202                }
203                if (oneRow && rhsValue <= 1.0 + tolerance) {
204                    if (!numberUnsatisfied) {
205                        for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
206                            int iColumn = column[j];
207                            if (dsort[iColumn] > 1.0e50) {
208                                dsort[iColumn] = 0;
209                                nSort++;
210                            }
211                        }
212                    }
213                }
214            }
215            // sort
216            CoinSort_2(dsort, dsort + numberColumns, sort);
217        } else {
218            // new way
219            for (i = 0; i < numberIntegers; i++) {
220                int iColumn = integerVariable[i];
221                if (upper[iColumn] > lower[iColumn]) {
222                    if (!mark_ || !mark_[iColumn]) {
223                        double distanceDown = solution[iColumn] - lower[iColumn];
224                        double distanceUp = upper[iColumn] - solution[iColumn];
225                        double distance = CoinMin(distanceDown, distanceUp);
226                        if (distance > 0.001 && distance < 0.5) {
227                            dsort[nSort] = distance;
228                            sort[nSort++] = iColumn;
229                        }
230                    }
231                }
232            }
233            // sort
234            CoinSort_2(dsort, dsort + nSort, sort);
235            int n = 0;
236            double sum = 0.0;
237            for (int k = 0; k < nSort; k++) {
238                sum += dsort[k];
239                if (sum <= djTolerance_)
240                    n = k;
241                else
242                    break;
243            }
244            nSort = CoinMin(n, numberClean_ / 1000000);
245        }
246    } else {
247#define FIX_IF_LESS -0.1
248        // 3 in same row and sum <FIX_IF_LESS?
249        int numberRows = matrixByRow_.getNumRows();
250        const double * solution = model_->testSolution();
251        const int * column = matrixByRow_.getIndices();
252        const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
253        const int * rowLength = matrixByRow_.getVectorLengths();
254        double bestSum = 1.0;
255        int nBest = -1;
256        int kRow = -1;
257        OsiSolverInterface * solver = model_->solver();
258        for (int i = 0; i < numberRows; i++) {
259            int numberUnsatisfied = 0;
260            double sum = 0.0;
261            for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
262                int iColumn = column[j];
263                if (solver->isInteger(iColumn)) {
264                    double solValue = solution[iColumn];
265                    if (solValue > 1.0e-5 && solValue < FIX_IF_LESS) {
266                        numberUnsatisfied++;
267                        sum += solValue;
268                    }
269                }
270            }
271            if (numberUnsatisfied >= 3 && sum < FIX_IF_LESS) {
272                // possible
273                if (numberUnsatisfied > nBest ||
274                        (numberUnsatisfied == nBest && sum < bestSum)) {
275                    nBest = numberUnsatisfied;
276                    bestSum = sum;
277                    kRow = i;
278                }
279            }
280        }
281        assert (nBest > 0);
282        for (int j = rowStart[kRow]; j < rowStart[kRow] + rowLength[kRow]; j++) {
283            int iColumn = column[j];
284            if (solver->isInteger(iColumn)) {
285                double solValue = solution[iColumn];
286                if (solValue > 1.0e-5 && solValue < FIX_IF_LESS) {
287                    sort[nSort++] = iColumn;
288                }
289            }
290        }
291    }
292    OsiRowCut down;
293    down.setLb(-COIN_DBL_MAX);
294    double rhs = 0.0;
295    for (i = 0; i < nSort; i++) {
296        int iColumn = sort[i];
297        double distanceDown = solution[iColumn] - lower[iColumn];
298        double distanceUp = upper[iColumn] - solution[iColumn];
299        if (distanceDown < distanceUp) {
300            rhs += lower[iColumn];
301            dsort[i] = 1.0;
302        } else {
303            rhs -= upper[iColumn];
304            dsort[i] = -1.0;
305        }
306    }
307    down.setUb(rhs);
308    down.setRow(nSort, sort, dsort);
309    down.setEffectiveness(COIN_DBL_MAX); // so will persist
310    delete [] sort;
311    delete [] dsort;
312    // up is same - just with rhs changed
313    OsiRowCut up = down;
314    up.setLb(rhs + 1.0);
315    up.setUb(COIN_DBL_MAX);
316    // Say can fix one way
317    CbcCutBranchingObject * newObject =
318        new CbcCutBranchingObject(model_, down, up, true);
319    if (model_->messageHandler()->logLevel() > 1)
320        printf("creating cut in CbcBranchCut\n");
321    return newObject;
322}
323/* Does a lot of the work,
324   Returns 0 if no good, 1 if dj, 2 if clean, 3 if both
325   10 if branching on ones away from bound
326*/
327int
328CbcBranchToFixLots::shallWe() const
329{
330    int returnCode = 0;
331    OsiSolverInterface * solver = model_->solver();
332    int numberRows = matrixByRow_.getNumRows();
333    //if (numberRows!=solver->getNumRows())
334    //return 0;
335    const double * solution = model_->testSolution();
336    const double * lower = solver->getColLower();
337    const double * upper = solver->getColUpper();
338    const double * dj = solver->getReducedCost();
339    int i;
340    int numberIntegers = model_->numberIntegers();
341    const int * integerVariable = model_->integerVariable();
342    if (numberClean_ > 1000000) {
343        int wanted = numberClean_ % 1000000;
344        int * sort = new int[numberIntegers];
345        double * dsort = new double[numberIntegers];
346        int nSort = 0;
347        for (i = 0; i < numberIntegers; i++) {
348            int iColumn = integerVariable[i];
349            if (upper[iColumn] > lower[iColumn]) {
350                if (!mark_ || !mark_[iColumn]) {
351                    double distanceDown = solution[iColumn] - lower[iColumn];
352                    double distanceUp = upper[iColumn] - solution[iColumn];
353                    double distance = CoinMin(distanceDown, distanceUp);
354                    if (distance > 0.001 && distance < 0.5) {
355                        dsort[nSort] = distance;
356                        sort[nSort++] = iColumn;
357                    }
358                }
359            }
360        }
361        // sort
362        CoinSort_2(dsort, dsort + nSort, sort);
363        int n = 0;
364        double sum = 0.0;
365        for (int k = 0; k < nSort; k++) {
366            sum += dsort[k];
367            if (sum <= djTolerance_)
368                n = k;
369            else
370                break;
371        }
372        delete [] sort;
373        delete [] dsort;
374        return (n >= wanted) ? 10 : 0;
375    }
376    double integerTolerance =
377        model_->getDblParam(CbcModel::CbcIntegerTolerance);
378    // make smaller ?
379    double tolerance = CoinMin(1.0e-8, integerTolerance);
380    // How many fixed are we aiming at
381    int wantedFixed = static_cast<int> (static_cast<double>(numberIntegers) * fractionFixed_);
382    if (djTolerance_ < 1.0e10) {
383        int nSort = 0;
384        int numberFixed = 0;
385        for (i = 0; i < numberIntegers; i++) {
386            int iColumn = integerVariable[i];
387            if (upper[iColumn] > lower[iColumn]) {
388                if (!mark_ || !mark_[iColumn]) {
389                    if (solution[iColumn] < lower[iColumn] + tolerance) {
390                        if (dj[iColumn] > djTolerance_) {
391                            nSort++;
392                        }
393                    } else if (solution[iColumn] > upper[iColumn] - tolerance) {
394                        if (dj[iColumn] < -djTolerance_) {
395                            nSort++;
396                        }
397                    }
398                }
399            } else {
400                numberFixed++;
401            }
402        }
403        if (numberFixed + nSort < wantedFixed && !alwaysCreate_) {
404            returnCode = 0;
405        } else if (numberFixed < wantedFixed) {
406            returnCode = 1;
407        } else {
408            returnCode = 0;
409        }
410    }
411    if (numberClean_) {
412        // see how many rows clean
413        int i;
414        //const double * rowLower = solver->getRowLower();
415        const double * rowUpper = solver->getRowUpper();
416        // Row copy
417        const double * elementByRow = matrixByRow_.getElements();
418        const int * column = matrixByRow_.getIndices();
419        const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
420        const int * rowLength = matrixByRow_.getVectorLengths();
421        const double * columnLower = solver->getColLower();
422        const double * columnUpper = solver->getColUpper();
423        const double * solution = solver->getColSolution();
424        int numberClean = 0;
425        bool someToDoYet = false;
426        int numberColumns = solver->getNumCols();
427        char * mark = new char[numberColumns];
428        int numberFixed = 0;
429        for (i = 0; i < numberColumns; i++) {
430            if (columnLower[i] != columnUpper[i]) {
431                mark[i] = 0;
432            } else {
433                mark[i] = 1;
434                numberFixed++;
435            }
436        }
437        int numberNewFixed = 0;
438        for (i = 0; i < numberRows; i++) {
439            double rhsValue = rowUpper[i];
440            bool oneRow = true;
441            // check elements
442            int numberUnsatisfied = 0;
443            for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
444                int iColumn = column[j];
445                double value = elementByRow[j];
446                double solValue = solution[iColumn];
447                if (columnLower[iColumn] != columnUpper[iColumn]) {
448                    if (solValue < 1.0 - integerTolerance && solValue > integerTolerance)
449                        numberUnsatisfied++;
450                    if (value != 1.0) {
451                        oneRow = false;
452                        break;
453                    }
454                } else {
455                    rhsValue -= value * floor(solValue + 0.5);
456                }
457            }
458            if (oneRow && rhsValue <= 1.0 + tolerance) {
459                if (numberUnsatisfied) {
460                    someToDoYet = true;
461                } else {
462                    numberClean++;
463                    for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
464                        int iColumn = column[j];
465                        if (columnLower[iColumn] != columnUpper[iColumn] && !mark[iColumn]) {
466                            mark[iColumn] = 1;
467                            numberNewFixed++;
468                        }
469                    }
470                }
471            }
472        }
473        delete [] mark;
474        //printf("%d clean, %d old fixed, %d new fixed\n",
475        //   numberClean,numberFixed,numberNewFixed);
476        if (someToDoYet && numberClean < numberClean_
477                && numberNewFixed + numberFixed < wantedFixed) {
478        } else if (numberFixed < wantedFixed) {
479            returnCode |= 2;
480        } else {
481        }
482    }
483    return returnCode;
484}
485double
486CbcBranchToFixLots::infeasibility(const OsiBranchingInformation * /*info*/,
487                                  int &preferredWay) const
488{
489    preferredWay = -1;
490    CbcNode * node = model_->currentNode();
491    int depth;
492    if (node)
493        depth = CoinMax(node->depth(), 0);
494    else
495        return 0.0;
496    if (depth_ < 0) {
497        return 0.0;
498    } else if (depth_ > 0) {
499        if ((depth % depth_) != 0)
500            return 0.0;
501    }
502    if (djTolerance_ != -1.234567) {
503        if (!shallWe())
504            return 0.0;
505        else
506            return 1.0e20;
507    } else {
508        // See if 3 in same row and sum <FIX_IF_LESS?
509        int numberRows = matrixByRow_.getNumRows();
510        const double * solution = model_->testSolution();
511        const int * column = matrixByRow_.getIndices();
512        const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
513        const int * rowLength = matrixByRow_.getVectorLengths();
514        double bestSum = 1.0;
515        int nBest = -1;
516        OsiSolverInterface * solver = model_->solver();
517        for (int i = 0; i < numberRows; i++) {
518            int numberUnsatisfied = 0;
519            double sum = 0.0;
520            for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
521                int iColumn = column[j];
522                if (solver->isInteger(iColumn)) {
523                    double solValue = solution[iColumn];
524                    if (solValue > 1.0e-5 && solValue < FIX_IF_LESS) {
525                        numberUnsatisfied++;
526                        sum += solValue;
527                    }
528                }
529            }
530            if (numberUnsatisfied >= 3 && sum < FIX_IF_LESS) {
531                // possible
532                if (numberUnsatisfied > nBest ||
533                        (numberUnsatisfied == nBest && sum < bestSum)) {
534                    nBest = numberUnsatisfied;
535                    bestSum = sum;
536                }
537            }
538        }
539        if (nBest > 0)
540            return 1.0e20;
541        else
542            return 0.0;
543    }
544}
545// Redoes data when sequence numbers change
546void
547CbcBranchToFixLots::redoSequenceEtc(CbcModel * model, int numberColumns, const int * originalColumns)
548{
549    model_ = model;
550    if (mark_) {
551        OsiSolverInterface * solver = model_->solver();
552        int numberColumnsNow = solver->getNumCols();
553        char * temp = new char[numberColumnsNow];
554        memset(temp, 0, numberColumnsNow);
555        for (int i = 0; i < numberColumns; i++) {
556            int j = originalColumns[i];
557            temp[i] = mark_[j];
558        }
559        delete [] mark_;
560        mark_ = temp;
561    }
562    OsiSolverInterface * solver = model_->solver();
563    matrixByRow_ = *solver->getMatrixByRow();
564}
565
Note: See TracBrowser for help on using the repository browser.