source: trunk/Cbc/src/CbcHeuristicRINS.cpp

Last change on this file was 2467, checked in by unxusr, 10 months ago

spaces after angles

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.9 KB
Line 
1/* $Id: CbcHeuristicRINS.cpp 2467 2019-01-03 21:26:29Z stefan $ */
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 CbcHeuristicRINS::generateCpp(FILE *fp)
100{
101  CbcHeuristicRINS other;
102  fprintf(fp, "0#include \"CbcHeuristicRINS.hpp\"\n");
103  fprintf(fp, "3  CbcHeuristicRINS heuristicRINS(*cbcModel);\n");
104  CbcHeuristic::generateCpp(fp, "heuristicRINS");
105  if (howOften_ != other.howOften_)
106    fprintf(fp, "3  heuristicRINS.setHowOften(%d);\n", howOften_);
107  else
108    fprintf(fp, "4  heuristicRINS.setHowOften(%d);\n", howOften_);
109  fprintf(fp, "3  cbcModel->addHeuristic(&heuristicRINS);\n");
110}
111
112// Copy constructor
113CbcHeuristicRINS::CbcHeuristicRINS(const CbcHeuristicRINS &rhs)
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 CbcHeuristicRINS::resetModel(CbcModel * /*model*/)
132{
133  //CbcHeuristic::resetModel(model);
134  delete[] used_;
135  stateOfFixing_ = 0;
136  if (model_ && used_) {
137    int numberColumns = model_->solver()->getNumCols();
138    used_ = new char[numberColumns];
139    memset(used_, 0, numberColumns);
140  } else {
141    used_ = NULL;
142  }
143}
144/*
145  First tries setting a variable to better value.  If feasible then
146  tries setting others.  If not feasible then tries swaps
147  Returns 1 if solution, 0 if not */
148int CbcHeuristicRINS::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#ifdef HEURISTIC_INFORM
157  printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
158    heuristicName(), numRuns_, numCouldRun_, when_);
159#endif
160  if (numberSolutions_ < model_->getSolutionCount()) {
161    // new solution - add info
162    numberSolutions_ = model_->getSolutionCount();
163
164    OsiSolverInterface *solver = model_->solver();
165    int numberIntegers = model_->numberIntegers();
166    const int *integerVariable = model_->integerVariable();
167
168    int i;
169    for (i = 0; i < numberIntegers; i++) {
170      int iColumn = integerVariable[i];
171      if (!isHeuristicInteger(solver, iColumn))
172        continue;
173      const OsiObject *object = model_->object(i);
174      // get original bounds
175      double originalLower;
176      double originalUpper;
177      getIntegerInformation(object, originalLower, originalUpper);
178      double value = bestSolution[iColumn];
179      if (value < originalLower) {
180        value = originalLower;
181      } else if (value > originalUpper) {
182        value = originalUpper;
183      }
184      double nearest = floor(value + 0.5);
185      // if away from lower bound mark that fact
186      if (nearest > originalLower) {
187        used_[iColumn] = 1;
188      }
189    }
190  }
191  int numberNodes = model_->getNodeCount();
192  if (howOften_ == 100) {
193    if (numberNodes < lastNode_ + 12)
194      return 0;
195    // Do at 50 and 100
196    if ((numberNodes > 40 && numberNodes <= 50) || (numberNodes > 90 && numberNodes < 100))
197      numberNodes = howOften_;
198  }
199  // Allow for infeasible nodes - so do anyway after a bit
200  if (howOften_ >= 100 && numberNodes >= lastNode_ + 2 * howOften_) {
201    numberNodes = howOften_;
202  }
203  if ((numberNodes % howOften_) == 0 && (model_->getCurrentPassNumber() <= 1 || model_->getCurrentPassNumber() == 999999)) {
204    lastNode_ = model_->getNodeCount();
205    OsiSolverInterface *solver = model_->solver();
206
207    int numberIntegers = model_->numberIntegers();
208    const int *integerVariable = model_->integerVariable();
209
210    const double *currentSolution = solver->getColSolution();
211    const int *used = model_->usedInSolution();
212    OsiSolverInterface *newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
213    int numberColumns = newSolver->getNumCols();
214    int numberContinuous = numberColumns - numberIntegers;
215
216    double primalTolerance;
217    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
218
219    int i;
220    int nFix = 0;
221    for (i = 0; i < numberIntegers; i++) {
222      int iColumn = integerVariable[i];
223      if (!isHeuristicInteger(solver, iColumn))
224        continue;
225      const OsiObject *object = model_->object(i);
226      // get original bounds
227      double originalLower;
228      double originalUpper;
229      getIntegerInformation(object, originalLower, originalUpper);
230      double valueInt = bestSolution[iColumn];
231      if (valueInt < originalLower) {
232        valueInt = originalLower;
233      } else if (valueInt > originalUpper) {
234        valueInt = originalUpper;
235      }
236      if (fabs(currentSolution[iColumn] - valueInt) < 10.0 * primalTolerance) {
237        double nearest = floor(valueInt + 0.5);
238        /*
239                  shallowDepth_
240                  0 - normal
241                  1 - only fix if at lb
242                  2 - only fix if not at lb
243                  3 - only fix if at lb and !used
244                */
245        bool fix = false;
246        switch (shallowDepth_) {
247        case 0:
248          fix = true;
249          break;
250        case 1:
251          if (nearest == originalLower)
252            fix = true;
253          break;
254        case 2:
255          if (nearest != originalLower)
256            fix = true;
257          break;
258        case 3:
259          if (nearest == originalLower && !used[iColumn])
260            fix = true;
261          break;
262        }
263        if (fix) {
264          newSolver->setColLower(iColumn, nearest);
265          newSolver->setColUpper(iColumn, nearest);
266          nFix++;
267        }
268      }
269    }
270    int divisor = 0;
271    if (5 * nFix > numberIntegers) {
272      if (numberContinuous > 2 * numberIntegers && ((nFix * 10 < numberColumns && !numRuns_ && numberTries_ > 2) || stateOfFixing_)) {
273#define RINS_FIX_CONTINUOUS
274#ifdef RINS_FIX_CONTINUOUS
275        const double *colLower = newSolver->getColLower();
276        //const double * colUpper = newSolver->getColUpper();
277        int nAtLb = 0;
278        //double sumDj=0.0;
279        const double *dj = newSolver->getReducedCost();
280        double direction = newSolver->getObjSense();
281        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
282          if (!isHeuristicInteger(newSolver, iColumn)) {
283            double value = bestSolution[iColumn];
284            if (value < colLower[iColumn] + 1.0e-8) {
285              //double djValue = dj[iColumn]*direction;
286              nAtLb++;
287              //sumDj += djValue;
288            }
289          }
290        }
291        if (nAtLb) {
292          // fix some continuous
293          double *sort = new double[nAtLb];
294          int *which = new int[nAtLb];
295          //double threshold = CoinMax((0.01*sumDj)/static_cast<double>(nAtLb),1.0e-6);
296          int nFix2 = 0;
297          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
298            if (!isHeuristicInteger(newSolver, iColumn)) {
299              double value = bestSolution[iColumn];
300              if (value < colLower[iColumn] + 1.0e-8) {
301                double djValue = dj[iColumn] * direction;
302                if (djValue > 1.0e-6) {
303                  sort[nFix2] = -djValue;
304                  which[nFix2++] = iColumn;
305                }
306              }
307            }
308          }
309          CoinSort_2(sort, sort + nFix2, which);
310          divisor = 4;
311          if (stateOfFixing_ > 0)
312            divisor = stateOfFixing_;
313          else if (stateOfFixing_ < -1)
314            divisor = (-stateOfFixing_) - 1;
315          nFix2 = CoinMin(nFix2, (numberColumns - nFix) / divisor);
316          for (int i = 0; i < nFix2; i++) {
317            int iColumn = which[i];
318            newSolver->setColUpper(iColumn, colLower[iColumn]);
319          }
320          delete[] sort;
321          delete[] which;
322#ifdef CLP_INVESTIGATE2
323          printf("%d integers have same value, and %d continuous fixed at lb\n",
324            nFix, nFix2);
325#endif
326        }
327#endif
328      }
329      if (solutionValue == -COIN_DBL_MAX) {
330        // return fixings in betterSolution
331        const double *colLower = newSolver->getColLower();
332        const double *colUpper = newSolver->getColUpper();
333        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
334          if (colLower[iColumn] == colUpper[iColumn])
335            betterSolution[iColumn] = colLower[iColumn];
336          else
337            betterSolution[iColumn] = COIN_DBL_MAX;
338        }
339        delete newSolver;
340        return 0;
341      }
342      //printf("RINS %d integers have same value\n",nFix);
343      returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
344        model_->getCutoff(), "CbcHeuristicRINS");
345      if (returnCode < 0) {
346        returnCode = 0; // returned on size
347        if (divisor) {
348          stateOfFixing_ = -divisor; // say failed
349        } else if (numberContinuous > 2 * numberIntegers && !numRuns_ && numberTries_ > 2) {
350          stateOfFixing_ = -4; //start fixing
351        }
352      } else {
353        numRuns_++;
354        if (divisor)
355          stateOfFixing_ = divisor; // say small enough
356      }
357      if ((returnCode & 1) != 0)
358        numberSuccesses_++;
359      //printf("return code %d",returnCode);
360      if ((returnCode & 2) != 0) {
361        // could add cut
362        returnCode &= ~2;
363        //printf("could add cut with %d elements (if all 0-1)\n",nFix);
364      } else {
365        //printf("\n");
366      }
367    }
368
369    numberTries_++;
370    if ((numberTries_ % 10) == 0 && numberSuccesses_ * 3 < numberTries_)
371      howOften_ += static_cast< int >(howOften_ * decayFactor_);
372    delete newSolver;
373  }
374  return returnCode;
375}
376// update model
377void CbcHeuristicRINS::setModel(CbcModel *model)
378{
379  model_ = model;
380  // Get a copy of original matrix
381  assert(model_->solver());
382  delete[] used_;
383  int numberColumns = model->solver()->getNumCols();
384  used_ = new char[numberColumns];
385  memset(used_, 0, numberColumns);
386}
387
388/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
389*/
Note: See TracBrowser for help on using the repository browser.