source: stable/2.6/Cbc/src/CbcHeuristicRENS.cpp @ 2122

Last change on this file since 2122 was 1499, checked in by forrest, 9 years ago

a few more rens heuristics

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