source: trunk/Cbc/src/CbcHeuristicRINS.cpp @ 1573

Last change on this file since 1573 was 1573, checked in by lou, 8 years ago

Change to EPL license notice.

File size: 12.0 KB
Line 
1/* $Id: CbcHeuristicRINS.cpp 1240 2009-10-02 18:41:44Z forrest $ */
2// Copyright (C) 2006, 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#if defined(_MSC_VER)
7// Turn off compiler warning about long names
8#  pragma warning(disable:4786)
9#endif
10#include <cassert>
11#include <cstdlib>
12#include <cmath>
13#include <cfloat>
14
15#include "OsiSolverInterface.hpp"
16#include "CbcModel.hpp"
17#include "CbcMessage.hpp"
18#include "CbcHeuristicRINS.hpp"
19#include "CbcBranchActual.hpp"
20#include "CbcStrategy.hpp"
21#include "CglPreProcess.hpp"
22
23// Default Constructor
24CbcHeuristicRINS::CbcHeuristicRINS()
25        : CbcHeuristic()
26{
27    numberSolutions_ = 0;
28    numberSuccesses_ = 0;
29    numberTries_ = 0;
30    stateOfFixing_ = 0;
31    lastNode_ = -999999;
32    howOften_ = 100;
33    decayFactor_ = 0.5;
34    used_ = NULL;
35    whereFrom_ = 1 + 8 + 16 + 255 * 256;
36    whereFrom_ = 1 + 8 + 255 * 256;
37}
38
39// Constructor with model - assumed before cuts
40
41CbcHeuristicRINS::CbcHeuristicRINS(CbcModel & model)
42        : CbcHeuristic(model)
43{
44    numberSolutions_ = 0;
45    numberSuccesses_ = 0;
46    numberTries_ = 0;
47    stateOfFixing_ = 0;
48    lastNode_ = -999999;
49    howOften_ = 100;
50    decayFactor_ = 0.5;
51    assert(model.solver());
52    int numberColumns = model.solver()->getNumCols();
53    used_ = new char[numberColumns];
54    memset(used_, 0, numberColumns);
55    whereFrom_ = 1 + 8 + 16 + 255 * 256;
56    whereFrom_ = 1 + 8 + 255 * 256;
57}
58
59// Destructor
60CbcHeuristicRINS::~CbcHeuristicRINS ()
61{
62    delete [] used_;
63}
64
65// Clone
66CbcHeuristic *
67CbcHeuristicRINS::clone() const
68{
69    return new CbcHeuristicRINS(*this);
70}
71
72// Assignment operator
73CbcHeuristicRINS &
74CbcHeuristicRINS::operator=( const CbcHeuristicRINS & rhs)
75{
76    if (this != &rhs) {
77        CbcHeuristic::operator=(rhs);
78        numberSolutions_ = rhs.numberSolutions_;
79        howOften_ = rhs.howOften_;
80        numberSuccesses_ = rhs.numberSuccesses_;
81        numberTries_ = rhs.numberTries_;
82        stateOfFixing_ = rhs.stateOfFixing_;
83        lastNode_ = rhs.lastNode_;
84        delete [] used_;
85        if (model_ && rhs.used_) {
86            int numberColumns = model_->solver()->getNumCols();
87            used_ = new char[numberColumns];
88            memcpy(used_, rhs.used_, numberColumns);
89        } else {
90            used_ = NULL;
91        }
92    }
93    return *this;
94}
95
96// Create C++ lines to get to current state
97void
98CbcHeuristicRINS::generateCpp( FILE * fp)
99{
100    CbcHeuristicRINS other;
101    fprintf(fp, "0#include \"CbcHeuristicRINS.hpp\"\n");
102    fprintf(fp, "3  CbcHeuristicRINS heuristicRINS(*cbcModel);\n");
103    CbcHeuristic::generateCpp(fp, "heuristicRINS");
104    if (howOften_ != other.howOften_)
105        fprintf(fp, "3  heuristicRINS.setHowOften(%d);\n", howOften_);
106    else
107        fprintf(fp, "4  heuristicRINS.setHowOften(%d);\n", howOften_);
108    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicRINS);\n");
109}
110
111// Copy constructor
112CbcHeuristicRINS::CbcHeuristicRINS(const CbcHeuristicRINS & rhs)
113        :
114        CbcHeuristic(rhs),
115        numberSolutions_(rhs.numberSolutions_),
116        howOften_(rhs.howOften_),
117        numberSuccesses_(rhs.numberSuccesses_),
118        numberTries_(rhs.numberTries_),
119        stateOfFixing_(rhs.stateOfFixing_),
120        lastNode_(rhs.lastNode_)
121{
122    if (model_ && rhs.used_) {
123        int numberColumns = model_->solver()->getNumCols();
124        used_ = new char[numberColumns];
125        memcpy(used_, rhs.used_, numberColumns);
126    } else {
127        used_ = NULL;
128    }
129}
130// Resets stuff if model changes
131void
132CbcHeuristicRINS::resetModel(CbcModel * /*model*/)
133{
134    //CbcHeuristic::resetModel(model);
135    delete [] used_;
136    stateOfFixing_ = 0;
137    if (model_ && used_) {
138        int numberColumns = model_->solver()->getNumCols();
139        used_ = new char[numberColumns];
140        memset(used_, 0, numberColumns);
141    } else {
142        used_ = NULL;
143    }
144}
145/*
146  First tries setting a variable to better value.  If feasible then
147  tries setting others.  If not feasible then tries swaps
148  Returns 1 if solution, 0 if not */
149int
150CbcHeuristicRINS::solution(double & solutionValue,
151                           double * betterSolution)
152{
153    numCouldRun_++;
154    int returnCode = 0;
155    const double * bestSolution = model_->bestSolution();
156    if (!bestSolution)
157        return 0; // No solution found yet
158    if (numberSolutions_ < model_->getSolutionCount()) {
159        // new solution - add info
160        numberSolutions_ = model_->getSolutionCount();
161
162        int numberIntegers = model_->numberIntegers();
163        const int * integerVariable = model_->integerVariable();
164
165        int i;
166        for (i = 0; i < numberIntegers; i++) {
167            int iColumn = integerVariable[i];
168            const OsiObject * object = model_->object(i);
169            // get original bounds
170            double originalLower;
171            double originalUpper;
172            getIntegerInformation( object, originalLower, originalUpper);
173            double value = bestSolution[iColumn];
174            if (value < originalLower) {
175                value = originalLower;
176            } else if (value > originalUpper) {
177                value = originalUpper;
178            }
179            double nearest = floor(value + 0.5);
180            // if away from lower bound mark that fact
181            if (nearest > originalLower) {
182                used_[iColumn] = 1;
183            }
184        }
185    }
186    int numberNodes = model_->getNodeCount();
187    if (howOften_ == 100) {
188        if (numberNodes < lastNode_ + 12)
189            return 0;
190        // Do at 50 and 100
191        if ((numberNodes > 40 && numberNodes <= 50) || (numberNodes > 90 && numberNodes < 100))
192            numberNodes = howOften_;
193    }
194    // Allow for infeasible nodes - so do anyway after a bit
195    if (howOften_ >= 100 && numberNodes >= lastNode_ + 2*howOften_) {
196        numberNodes = howOften_;
197    }
198    if ((numberNodes % howOften_) == 0 && (model_->getCurrentPassNumber() == 1 ||
199                                           model_->getCurrentPassNumber() == 999999)) {
200        lastNode_ = model_->getNodeCount();
201        OsiSolverInterface * solver = model_->solver();
202
203        int numberIntegers = model_->numberIntegers();
204        const int * integerVariable = model_->integerVariable();
205
206        const double * currentSolution = solver->getColSolution();
207        OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
208        int numberColumns = newSolver->getNumCols();
209        int numberContinuous = numberColumns - numberIntegers;
210
211        double primalTolerance;
212        solver->getDblParam(OsiPrimalTolerance, primalTolerance);
213
214        int i;
215        int nFix = 0;
216        for (i = 0; i < numberIntegers; i++) {
217            int iColumn = integerVariable[i];
218            const OsiObject * object = model_->object(i);
219            // get original bounds
220            double originalLower;
221            double originalUpper;
222            getIntegerInformation( object, originalLower, originalUpper);
223            double valueInt = bestSolution[iColumn];
224            if (valueInt < originalLower) {
225                valueInt = originalLower;
226            } else if (valueInt > originalUpper) {
227                valueInt = originalUpper;
228            }
229            if (fabs(currentSolution[iColumn] - valueInt) < 10.0*primalTolerance) {
230                double nearest = floor(valueInt + 0.5);
231                newSolver->setColLower(iColumn, nearest);
232                newSolver->setColUpper(iColumn, nearest);
233                nFix++;
234            }
235        }
236        int divisor = 0;
237        if (5*nFix > numberIntegers) {
238            if (numberContinuous > 2*numberIntegers && ((nFix*10 < numberColumns &&
239                    !numRuns_ && numberTries_ > 2) ||
240                    stateOfFixing_)) {
241#define RINS_FIX_CONTINUOUS
242#ifdef RINS_FIX_CONTINUOUS
243                const double * colLower = newSolver->getColLower();
244                //const double * colUpper = newSolver->getColUpper();
245                int nAtLb = 0;
246                //double sumDj=0.0;
247                const double * dj = newSolver->getReducedCost();
248                double direction = newSolver->getObjSense();
249                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
250                    if (!newSolver->isInteger(iColumn)) {
251                        double value = bestSolution[iColumn];
252                        if (value < colLower[iColumn] + 1.0e-8) {
253                            //double djValue = dj[iColumn]*direction;
254                            nAtLb++;
255                            //sumDj += djValue;
256                        }
257                    }
258                }
259                if (nAtLb) {
260                    // fix some continuous
261                    double * sort = new double[nAtLb];
262                    int * which = new int [nAtLb];
263                    //double threshold = CoinMax((0.01*sumDj)/static_cast<double>(nAtLb),1.0e-6);
264                    int nFix2 = 0;
265                    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
266                        if (!newSolver->isInteger(iColumn)) {
267                            double value = bestSolution[iColumn];
268                            if (value < colLower[iColumn] + 1.0e-8) {
269                                double djValue = dj[iColumn] * direction;
270                                if (djValue > 1.0e-6) {
271                                    sort[nFix2] = -djValue;
272                                    which[nFix2++] = iColumn;
273                                }
274                            }
275                        }
276                    }
277                    CoinSort_2(sort, sort + nFix2, which);
278                    divisor = 4;
279                    if (stateOfFixing_ > 0)
280                        divisor = stateOfFixing_;
281                    else if (stateOfFixing_ < -1)
282                        divisor = (-stateOfFixing_) - 1;
283                    nFix2 = CoinMin(nFix2, (numberColumns - nFix) / divisor);
284                    for (int i = 0; i < nFix2; i++) {
285                        int iColumn = which[i];
286                        newSolver->setColUpper(iColumn, colLower[iColumn]);
287                    }
288                    delete [] sort;
289                    delete [] which;
290#ifdef CLP_INVESTIGATE2
291                    printf("%d integers have same value, and %d continuous fixed at lb\n",
292                           nFix, nFix2);
293#endif
294                }
295#endif
296            }
297            //printf("%d integers have same value\n",nFix);
298            returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
299                                             model_->getCutoff(), "CbcHeuristicRINS");
300            if (returnCode < 0) {
301                returnCode = 0; // returned on size
302                if (divisor) {
303                    stateOfFixing_ = - divisor; // say failed
304                } else if (numberContinuous > 2*numberIntegers &&
305                           !numRuns_ && numberTries_ > 2) {
306                    stateOfFixing_ = -4; //start fixing
307                }
308            } else {
309                numRuns_++;
310                if (divisor)
311                    stateOfFixing_ =  divisor; // say small enough
312            }
313            if ((returnCode&1) != 0)
314                numberSuccesses_++;
315            //printf("return code %d",returnCode);
316            if ((returnCode&2) != 0) {
317                // could add cut
318                returnCode &= ~2;
319                //printf("could add cut with %d elements (if all 0-1)\n",nFix);
320            } else {
321                //printf("\n");
322            }
323        }
324
325        numberTries_++;
326        if ((numberTries_ % 10) == 0 && numberSuccesses_*3 < numberTries_)
327            howOften_ += static_cast<int> (howOften_ * decayFactor_);
328        delete newSolver;
329    }
330    return returnCode;
331}
332// update model
333void CbcHeuristicRINS::setModel(CbcModel * model)
334{
335    model_ = model;
336    // Get a copy of original matrix
337    assert(model_->solver());
338    delete [] used_;
339    int numberColumns = model->solver()->getNumCols();
340    used_ = new char[numberColumns];
341    memset(used_, 0, numberColumns);
342}
343
344
345
Note: See TracBrowser for help on using the repository browser.