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

Last change on this file since 1422 was 1357, checked in by coin, 9 years ago

run 'astyle -A4 -p' and dos2unix

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}
Note: See TracBrowser for help on using the repository browser.