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

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

Broke out CbcHeuristicRENS, CbcHeuristicDINS and CbcHeuristicVND out of CbcHeuristicRINS.

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