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

Last change on this file since 2280 was 2280, checked in by forrest, 5 years ago

allow heuristics to see if integers are 'optional'
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.2 KB
Line 
1/* $Id: CbcHeuristicRINS.cpp 2280 2016-06-14 14:39:54Z 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    shallowDepth_ = 0;
32    lastNode_ = -999999;
33    howOften_ = 100;
34    decayFactor_ = 0.5;
35    used_ = NULL;
36    whereFrom_ = 1 + 8 + 16 + 255 * 256;
37    whereFrom_ = 1 + 8 + 255 * 256;
38}
39
40// Constructor with model - assumed before cuts
41
42CbcHeuristicRINS::CbcHeuristicRINS(CbcModel & model)
43        : CbcHeuristic(model)
44{
45    numberSolutions_ = 0;
46    numberSuccesses_ = 0;
47    numberTries_ = 0;
48    stateOfFixing_ = 0;
49    shallowDepth_ = 0;
50    lastNode_ = -999999;
51    howOften_ = 100;
52    decayFactor_ = 0.5;
53    assert(model.solver());
54    int numberColumns = model.solver()->getNumCols();
55    used_ = new char[numberColumns];
56    memset(used_, 0, numberColumns);
57    whereFrom_ = 1 + 8 + 16 + 255 * 256;
58    whereFrom_ = 1 + 8 + 255 * 256;
59}
60
61// Destructor
62CbcHeuristicRINS::~CbcHeuristicRINS ()
63{
64    delete [] used_;
65}
66
67// Clone
68CbcHeuristic *
69CbcHeuristicRINS::clone() const
70{
71    return new CbcHeuristicRINS(*this);
72}
73
74// Assignment operator
75CbcHeuristicRINS &
76CbcHeuristicRINS::operator=( const CbcHeuristicRINS & rhs)
77{
78    if (this != &rhs) {
79        CbcHeuristic::operator=(rhs);
80        numberSolutions_ = rhs.numberSolutions_;
81        howOften_ = rhs.howOften_;
82        numberSuccesses_ = rhs.numberSuccesses_;
83        numberTries_ = rhs.numberTries_;
84        stateOfFixing_ = rhs.stateOfFixing_;
85        lastNode_ = rhs.lastNode_;
86        delete [] used_;
87        if (model_ && rhs.used_) {
88            int numberColumns = model_->solver()->getNumCols();
89            used_ = new char[numberColumns];
90            memcpy(used_, rhs.used_, numberColumns);
91        } else {
92            used_ = NULL;
93        }
94    }
95    return *this;
96}
97
98// Create C++ lines to get to current state
99void
100CbcHeuristicRINS::generateCpp( FILE * fp)
101{
102    CbcHeuristicRINS other;
103    fprintf(fp, "0#include \"CbcHeuristicRINS.hpp\"\n");
104    fprintf(fp, "3  CbcHeuristicRINS heuristicRINS(*cbcModel);\n");
105    CbcHeuristic::generateCpp(fp, "heuristicRINS");
106    if (howOften_ != other.howOften_)
107        fprintf(fp, "3  heuristicRINS.setHowOften(%d);\n", howOften_);
108    else
109        fprintf(fp, "4  heuristicRINS.setHowOften(%d);\n", howOften_);
110    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicRINS);\n");
111}
112
113// Copy constructor
114CbcHeuristicRINS::CbcHeuristicRINS(const CbcHeuristicRINS & rhs)
115        :
116        CbcHeuristic(rhs),
117        numberSolutions_(rhs.numberSolutions_),
118        howOften_(rhs.howOften_),
119        numberSuccesses_(rhs.numberSuccesses_),
120        numberTries_(rhs.numberTries_),
121        stateOfFixing_(rhs.stateOfFixing_),
122        lastNode_(rhs.lastNode_)
123{
124    if (model_ && rhs.used_) {
125        int numberColumns = model_->solver()->getNumCols();
126        used_ = new char[numberColumns];
127        memcpy(used_, rhs.used_, numberColumns);
128    } else {
129        used_ = NULL;
130    }
131}
132// Resets stuff if model changes
133void
134CbcHeuristicRINS::resetModel(CbcModel * /*model*/)
135{
136    //CbcHeuristic::resetModel(model);
137    delete [] used_;
138    stateOfFixing_ = 0;
139    if (model_ && used_) {
140        int numberColumns = model_->solver()->getNumCols();
141        used_ = new char[numberColumns];
142        memset(used_, 0, numberColumns);
143    } else {
144        used_ = NULL;
145    }
146}
147/*
148  First tries setting a variable to better value.  If feasible then
149  tries setting others.  If not feasible then tries swaps
150  Returns 1 if solution, 0 if not */
151int
152CbcHeuristicRINS::solution(double & solutionValue,
153                           double * betterSolution)
154{
155    numCouldRun_++;
156    int returnCode = 0;
157    const double * bestSolution = model_->bestSolution();
158    if (!bestSolution)
159        return 0; // No solution found yet
160#ifdef HEURISTIC_INFORM
161    printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
162           heuristicName(),numRuns_,numCouldRun_,when_);
163#endif
164    if (numberSolutions_ < model_->getSolutionCount()) {
165        // new solution - add info
166        numberSolutions_ = model_->getSolutionCount();
167
168        OsiSolverInterface * solver = model_->solver();
169        int numberIntegers = model_->numberIntegers();
170        const int * integerVariable = model_->integerVariable();
171
172        int i;
173        for (i = 0; i < numberIntegers; i++) {
174            int iColumn = integerVariable[i];
175            if (!isHeuristicInteger(solver,iColumn))
176              continue;
177            const OsiObject * object = model_->object(i);
178            // get original bounds
179            double originalLower;
180            double originalUpper;
181            getIntegerInformation( object, originalLower, originalUpper);
182            double value = bestSolution[iColumn];
183            if (value < originalLower) {
184                value = originalLower;
185            } else if (value > originalUpper) {
186                value = originalUpper;
187            }
188            double nearest = floor(value + 0.5);
189            // if away from lower bound mark that fact
190            if (nearest > originalLower) {
191                used_[iColumn] = 1;
192            }
193        }
194    }
195    int numberNodes = model_->getNodeCount();
196    if (howOften_ == 100) {
197        if (numberNodes < lastNode_ + 12)
198            return 0;
199        // Do at 50 and 100
200        if ((numberNodes > 40 && numberNodes <= 50) || (numberNodes > 90 && numberNodes < 100))
201            numberNodes = howOften_;
202    }
203    // Allow for infeasible nodes - so do anyway after a bit
204    if (howOften_ >= 100 && numberNodes >= lastNode_ + 2*howOften_) {
205        numberNodes = howOften_;
206    }
207    if ((numberNodes % howOften_) == 0 && (model_->getCurrentPassNumber() <= 1 ||
208                                           model_->getCurrentPassNumber() == 999999)) {
209        lastNode_ = model_->getNodeCount();
210        OsiSolverInterface * solver = model_->solver();
211
212        int numberIntegers = model_->numberIntegers();
213        const int * integerVariable = model_->integerVariable();
214
215        const double * currentSolution = solver->getColSolution();
216        const int * used = model_->usedInSolution();
217        OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
218        int numberColumns = newSolver->getNumCols();
219        int numberContinuous = numberColumns - numberIntegers;
220
221        double primalTolerance;
222        solver->getDblParam(OsiPrimalTolerance, primalTolerance);
223
224        int i;
225        int nFix = 0;
226        for (i = 0; i < numberIntegers; i++) {
227            int iColumn = integerVariable[i];
228            if (!isHeuristicInteger(solver,iColumn))
229              continue;
230            const OsiObject * object = model_->object(i);
231            // get original bounds
232            double originalLower;
233            double originalUpper;
234            getIntegerInformation( object, originalLower, originalUpper);
235            double valueInt = bestSolution[iColumn];
236            if (valueInt < originalLower) {
237                valueInt = originalLower;
238            } else if (valueInt > originalUpper) {
239                valueInt = originalUpper;
240            }
241            if (fabs(currentSolution[iColumn] - valueInt) < 10.0*primalTolerance) {
242                double nearest = floor(valueInt + 0.5);
243                /*
244                  shallowDepth_
245                  0 - normal
246                  1 - only fix if at lb
247                  2 - only fix if not at lb
248                  3 - only fix if at lb and !used
249                */
250                bool fix=false;
251                switch (shallowDepth_) {
252                case 0:
253                  fix = true;
254                  break;
255                case 1:
256                if (nearest==originalLower) 
257                  fix = true;
258                  break;
259                case 2:
260                if (nearest!=originalLower) 
261                  fix = true;
262                  break;
263                case 3:
264                if (nearest==originalLower && !used[iColumn]) 
265                  fix = true;
266                  break;
267                }
268                if (fix) {
269                  newSolver->setColLower(iColumn, nearest);
270                  newSolver->setColUpper(iColumn, nearest);
271                  nFix++;
272                }
273            }
274        }
275        int divisor = 0;
276        if (5*nFix > numberIntegers) {
277            if (numberContinuous > 2*numberIntegers && ((nFix*10 < numberColumns &&
278                    !numRuns_ && numberTries_ > 2) ||
279                    stateOfFixing_)) {
280#define RINS_FIX_CONTINUOUS
281#ifdef RINS_FIX_CONTINUOUS
282                const double * colLower = newSolver->getColLower();
283                //const double * colUpper = newSolver->getColUpper();
284                int nAtLb = 0;
285                //double sumDj=0.0;
286                const double * dj = newSolver->getReducedCost();
287                double direction = newSolver->getObjSense();
288                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
289                    if (!isHeuristicInteger(newSolver,iColumn)) {
290                        double value = bestSolution[iColumn];
291                        if (value < colLower[iColumn] + 1.0e-8) {
292                            //double djValue = dj[iColumn]*direction;
293                            nAtLb++;
294                            //sumDj += djValue;
295                        }
296                    }
297                }
298                if (nAtLb) {
299                    // fix some continuous
300                    double * sort = new double[nAtLb];
301                    int * which = new int [nAtLb];
302                    //double threshold = CoinMax((0.01*sumDj)/static_cast<double>(nAtLb),1.0e-6);
303                    int nFix2 = 0;
304                    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
305                        if (!isHeuristicInteger(newSolver,iColumn)) {
306                            double value = bestSolution[iColumn];
307                            if (value < colLower[iColumn] + 1.0e-8) {
308                                double djValue = dj[iColumn] * direction;
309                                if (djValue > 1.0e-6) {
310                                    sort[nFix2] = -djValue;
311                                    which[nFix2++] = iColumn;
312                                }
313                            }
314                        }
315                    }
316                    CoinSort_2(sort, sort + nFix2, which);
317                    divisor = 4;
318                    if (stateOfFixing_ > 0)
319                        divisor = stateOfFixing_;
320                    else if (stateOfFixing_ < -1)
321                        divisor = (-stateOfFixing_) - 1;
322                    nFix2 = CoinMin(nFix2, (numberColumns - nFix) / divisor);
323                    for (int i = 0; i < nFix2; i++) {
324                        int iColumn = which[i];
325                        newSolver->setColUpper(iColumn, colLower[iColumn]);
326                    }
327                    delete [] sort;
328                    delete [] which;
329#ifdef CLP_INVESTIGATE2
330                    printf("%d integers have same value, and %d continuous fixed at lb\n",
331                           nFix, nFix2);
332#endif
333                }
334#endif
335            }
336            if (solutionValue==-COIN_DBL_MAX) {
337              // return fixings in betterSolution
338              const double * colLower = newSolver->getColLower();
339              const double * colUpper = newSolver->getColUpper();
340              for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
341                if (colLower[iColumn]==colUpper[iColumn])
342                  betterSolution[iColumn]=colLower[iColumn];
343                else
344                  betterSolution[iColumn]=COIN_DBL_MAX;
345              }
346              delete newSolver;
347              return 0;
348            }
349            //printf("RINS %d integers have same value\n",nFix);
350            returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
351                                             model_->getCutoff(), "CbcHeuristicRINS");
352            if (returnCode < 0) {
353                returnCode = 0; // returned on size
354                if (divisor) {
355                    stateOfFixing_ = - divisor; // say failed
356                } else if (numberContinuous > 2*numberIntegers &&
357                           !numRuns_ && numberTries_ > 2) {
358                    stateOfFixing_ = -4; //start fixing
359                }
360            } else {
361                numRuns_++;
362                if (divisor)
363                    stateOfFixing_ =  divisor; // say small enough
364            }
365            if ((returnCode&1) != 0)
366                numberSuccesses_++;
367            //printf("return code %d",returnCode);
368            if ((returnCode&2) != 0) {
369                // could add cut
370                returnCode &= ~2;
371                //printf("could add cut with %d elements (if all 0-1)\n",nFix);
372            } else {
373                //printf("\n");
374            }
375        }
376
377        numberTries_++;
378        if ((numberTries_ % 10) == 0 && numberSuccesses_*3 < numberTries_)
379            howOften_ += static_cast<int> (howOften_ * decayFactor_);
380        delete newSolver;
381    }
382    return returnCode;
383}
384// update model
385void CbcHeuristicRINS::setModel(CbcModel * model)
386{
387    model_ = model;
388    // Get a copy of original matrix
389    assert(model_->solver());
390    delete [] used_;
391    int numberColumns = model->solver()->getNumCols();
392    used_ = new char[numberColumns];
393    memset(used_, 0, numberColumns);
394}
395
396
397
Note: See TracBrowser for help on using the repository browser.