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

Last change on this file since 1899 was 1899, checked in by stefan, 6 years ago

fixup svn properties

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