source: branches/sandbox/Cbc/src/CbcHeuristicRENS.cpp @ 1389

Last change on this file since 1389 was 1368, checked in by EdwinStraver, 10 years ago

Broke out CbcHeuristicRENS, CbcHeuristicDINS and CbcHeuristicVND out of CbcHeuristicRINS.

File size: 10.2 KB
Line 
1// edwin 12/5/09 carved out of CbcHeuristicRINS
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
11#include "OsiSolverInterface.hpp"
12#include "CbcModel.hpp"
13#include "CbcMessage.hpp"
14#include "CbcHeuristicRENS.hpp"
15#include "CbcBranchActual.hpp"
16#include "CbcStrategy.hpp"
17#include "CglPreProcess.hpp"
18
19// Default Constructor
20CbcHeuristicRENS::CbcHeuristicRENS()
21        : CbcHeuristic()
22{
23    numberTries_ = 0;
24    whereFrom_ = 256 + 1;
25}
26
27// Constructor with model - assumed before cuts
28
29CbcHeuristicRENS::CbcHeuristicRENS(CbcModel & model)
30        : CbcHeuristic(model)
31{
32    numberTries_ = 0;
33    whereFrom_ = 256 + 1;
34}
35
36// Destructor
37CbcHeuristicRENS::~CbcHeuristicRENS ()
38{
39}
40
41// Clone
42CbcHeuristic *
43CbcHeuristicRENS::clone() const
44{
45    return new CbcHeuristicRENS(*this);
46}
47
48// Assignment operator
49CbcHeuristicRENS &
50CbcHeuristicRENS::operator=( const CbcHeuristicRENS & rhs)
51{
52    if (this != &rhs) {
53        CbcHeuristic::operator=(rhs);
54        numberTries_ = rhs.numberTries_;
55    }
56    return *this;
57}
58
59// Copy constructor
60CbcHeuristicRENS::CbcHeuristicRENS(const CbcHeuristicRENS & rhs)
61        :
62        CbcHeuristic(rhs),
63        numberTries_(rhs.numberTries_)
64{
65}
66// Resets stuff if model changes
67void
68CbcHeuristicRENS::resetModel(CbcModel * )
69{
70}
71int
72CbcHeuristicRENS::solution(double & solutionValue,
73                           double * betterSolution)
74{
75    int returnCode = 0;
76    const double * bestSolution = model_->bestSolution();
77    if (numberTries_ || (when() < 2 && bestSolution))
78        return 0;
79    numberTries_++;
80    OsiSolverInterface * solver = model_->solver();
81
82    int numberIntegers = model_->numberIntegers();
83    const int * integerVariable = model_->integerVariable();
84
85    const double * currentSolution = solver->getColSolution();
86    OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
87    const double * colLower = newSolver->getColLower();
88    const double * colUpper = newSolver->getColUpper();
89
90    double primalTolerance;
91    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
92
93    int i;
94    int numberFixed = 0;
95    int numberTightened = 0;
96    int numberAtBound = 0;
97    int numberColumns = newSolver->getNumCols();
98    int numberContinuous = numberColumns - numberIntegers;
99
100    for (i = 0; i < numberIntegers; i++) {
101        int iColumn = integerVariable[i];
102        double value = currentSolution[iColumn];
103        double lower = colLower[iColumn];
104        double upper = colUpper[iColumn];
105        value = CoinMax(value, lower);
106        value = CoinMin(value, upper);
107#define RENS_FIX_ONLY_LOWER
108#ifndef RENS_FIX_ONLY_LOWER
109        if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
110            value = floor(value + 0.5);
111            if (value == lower || value == upper)
112                numberAtBound++;
113            newSolver->setColLower(iColumn, value);
114            newSolver->setColUpper(iColumn, value);
115            numberFixed++;
116        } else if (colUpper[iColumn] - colLower[iColumn] >= 2.0) {
117            numberTightened++;
118            newSolver->setColLower(iColumn, floor(value));
119            newSolver->setColUpper(iColumn, ceil(value));
120        }
121#else
122        if (fabs(value - floor(value + 0.5)) < 1.0e-8 &&
123                floor(value + 0.5) == lower) {
124            value = floor(value + 0.5);
125            numberAtBound++;
126            newSolver->setColLower(iColumn, value);
127            newSolver->setColUpper(iColumn, value);
128            numberFixed++;
129        } else if (colUpper[iColumn] - colLower[iColumn] >= 2.0) {
130            numberTightened++;
131            if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
132                value = floor(value + 0.5);
133                if (value < upper) {
134                    newSolver->setColLower(iColumn, CoinMax(value - 1.0, lower));
135                    newSolver->setColUpper(iColumn, CoinMin(value + 1.0, upper));
136                } else {
137                    newSolver->setColLower(iColumn, upper - 1.0);
138                }
139            } else {
140                newSolver->setColLower(iColumn, floor(value));
141                newSolver->setColUpper(iColumn, ceil(value));
142            }
143        }
144#endif
145    }
146    if (numberFixed > numberIntegers / 5) {
147        if (numberContinuous > numberIntegers && numberFixed < numberColumns / 5) {
148#define RENS_FIX_CONTINUOUS
149#ifdef RENS_FIX_CONTINUOUS
150            const double * colLower = newSolver->getColLower();
151            //const double * colUpper = newSolver->getColUpper();
152            int nAtLb = 0;
153            double sumDj = 0.0;
154            const double * dj = newSolver->getReducedCost();
155            double direction = newSolver->getObjSense();
156            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
157                if (!newSolver->isInteger(iColumn)) {
158                    double value = currentSolution[iColumn];
159                    if (value < colLower[iColumn] + 1.0e-8) {
160                        double djValue = dj[iColumn] * direction;
161                        nAtLb++;
162                        sumDj += djValue;
163                    }
164                }
165            }
166            if (nAtLb) {
167                // fix some continuous
168                double * sort = new double[nAtLb];
169                int * which = new int [nAtLb];
170                double threshold = CoinMax((0.01 * sumDj) / static_cast<double>(nAtLb), 1.0e-6);
171                int nFix2 = 0;
172                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
173                    if (!newSolver->isInteger(iColumn)) {
174                        double value = currentSolution[iColumn];
175                        if (value < colLower[iColumn] + 1.0e-8) {
176                            double djValue = dj[iColumn] * direction;
177                            if (djValue > threshold) {
178                                sort[nFix2] = -djValue;
179                                which[nFix2++] = iColumn;
180                            }
181                        }
182                    }
183                }
184                CoinSort_2(sort, sort + nFix2, which);
185                nFix2 = CoinMin(nFix2, (numberColumns - numberFixed) / 2);
186                for (int i = 0; i < nFix2; i++) {
187                    int iColumn = which[i];
188                    newSolver->setColUpper(iColumn, colLower[iColumn]);
189                }
190                delete [] sort;
191                delete [] which;
192#ifdef CLP_INVESTIGATE2
193                printf("%d integers fixed (%d tightened) (%d at bound), and %d continuous fixed at lb\n",
194                       numberFixed, numberTightened, numberAtBound, nFix2);
195#endif
196            }
197#endif
198        }
199#ifdef COIN_DEVELOP
200        printf("%d integers fixed and %d tightened\n", numberFixed, numberTightened);
201#endif
202        returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
203                                         model_->getCutoff(), "CbcHeuristicRENS");
204        if (returnCode < 0) {
205            returnCode = 0; // returned on size
206#ifdef RENS_FIX_CONTINUOUS
207            if (numberContinuous > numberIntegers && numberFixed >= numberColumns / 5) {
208                const double * colLower = newSolver->getColLower();
209                //const double * colUpper = newSolver->getColUpper();
210                int nAtLb = 0;
211                double sumDj = 0.0;
212                const double * dj = newSolver->getReducedCost();
213                double direction = newSolver->getObjSense();
214                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
215                    if (!newSolver->isInteger(iColumn)) {
216                        double value = currentSolution[iColumn];
217                        if (value < colLower[iColumn] + 1.0e-8) {
218                            double djValue = dj[iColumn] * direction;
219                            nAtLb++;
220                            sumDj += djValue;
221                        }
222                    }
223                }
224                if (nAtLb) {
225                    // fix some continuous
226                    double * sort = new double[nAtLb];
227                    int * which = new int [nAtLb];
228                    double threshold = CoinMax((0.01 * sumDj) / static_cast<double>(nAtLb), 1.0e-6);
229                    int nFix2 = 0;
230                    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
231                        if (!newSolver->isInteger(iColumn)) {
232                            double value = currentSolution[iColumn];
233                            if (value < colLower[iColumn] + 1.0e-8) {
234                                double djValue = dj[iColumn] * direction;
235                                if (djValue > threshold) {
236                                    sort[nFix2] = -djValue;
237                                    which[nFix2++] = iColumn;
238                                }
239                            }
240                        }
241                    }
242                    CoinSort_2(sort, sort + nFix2, which);
243                    nFix2 = CoinMin(nFix2, (numberColumns - numberFixed) / 2);
244                    for (int i = 0; i < nFix2; i++) {
245                        int iColumn = which[i];
246                        newSolver->setColUpper(iColumn, colLower[iColumn]);
247                    }
248                    delete [] sort;
249                    delete [] which;
250#ifdef CLP_INVESTIGATE2
251                    printf("%d integers fixed (%d tightened) (%d at bound), and %d continuous fixed at lb\n",
252                           numberFixed, numberTightened, numberAtBound, nFix2);
253#endif
254                }
255                returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
256                                                 model_->getCutoff(), "CbcHeuristicRENS");
257#endif
258            }
259        }
260        //printf("return code %d",returnCode);
261        if ((returnCode&2) != 0) {
262            // could add cut
263            returnCode &= ~2;
264#ifdef COIN_DEVELOP
265            if (!numberTightened && numberFixed == numberAtBound)
266                printf("could add cut with %d elements\n", numberFixed);
267#endif
268        } else {
269            //printf("\n");
270        }
271    }
272
273    delete newSolver;
274    return returnCode;
275}
276// update model
277void CbcHeuristicRENS::setModel(CbcModel * model)
278{
279    model_ = model;
280}
Note: See TracBrowser for help on using the repository browser.