source: trunk/Cbc/src/CbcHeuristicDINS.cpp

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

spaces after angles

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.3 KB
Line 
1// $Id: CbcHeuristicDINS.cpp 2467 2019-01-03 21:26:29Z 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// edwin 12/5/09 carved out of CbcHeuristicRINS
7
8#if defined(_MSC_VER)
9// Turn off compiler warning about long names
10#pragma warning(disable : 4786)
11#endif
12#include <cassert>
13#include <cstdlib>
14#include <cmath>
15#include <cfloat>
16
17#include "OsiSolverInterface.hpp"
18#include "CbcModel.hpp"
19#include "CbcMessage.hpp"
20#include "CbcHeuristicDINS.hpp"
21#include "CbcBranchActual.hpp"
22#include "CbcStrategy.hpp"
23#include "CglPreProcess.hpp"
24
25// Default Constructor
26CbcHeuristicDINS::CbcHeuristicDINS()
27  : CbcHeuristic()
28{
29  numberSolutions_ = 0;
30  numberSuccesses_ = 0;
31  numberTries_ = 0;
32  howOften_ = 100;
33  decayFactor_ = 0.5;
34  maximumKeepSolutions_ = 5;
35  numberKeptSolutions_ = 0;
36  numberIntegers_ = -1;
37  localSpace_ = 10;
38  values_ = NULL;
39}
40
41// Constructor with model - assumed before cuts
42
43CbcHeuristicDINS::CbcHeuristicDINS(CbcModel &model)
44  : CbcHeuristic(model)
45{
46  numberSolutions_ = 0;
47  numberSuccesses_ = 0;
48  numberTries_ = 0;
49  howOften_ = 100;
50  decayFactor_ = 0.5;
51  assert(model.solver());
52  maximumKeepSolutions_ = 5;
53  numberKeptSolutions_ = 0;
54  numberIntegers_ = -1;
55  localSpace_ = 10;
56  values_ = NULL;
57}
58
59// Destructor
60CbcHeuristicDINS::~CbcHeuristicDINS()
61{
62  for (int i = 0; i < numberKeptSolutions_; i++)
63    delete[] values_[i];
64  delete[] values_;
65}
66
67// Clone
68CbcHeuristic *
69CbcHeuristicDINS::clone() const
70{
71  return new CbcHeuristicDINS(*this);
72}
73
74// Assignment operator
75CbcHeuristicDINS &
76CbcHeuristicDINS::operator=(const CbcHeuristicDINS &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    for (int i = 0; i < numberKeptSolutions_; i++)
85      delete[] values_[i];
86    delete[] values_;
87    maximumKeepSolutions_ = rhs.maximumKeepSolutions_;
88    numberKeptSolutions_ = rhs.numberKeptSolutions_;
89    numberIntegers_ = rhs.numberIntegers_;
90    localSpace_ = rhs.localSpace_;
91    if (model_ && rhs.values_) {
92      assert(numberIntegers_ >= 0);
93      values_ = new int *[maximumKeepSolutions_];
94      for (int i = 0; i < maximumKeepSolutions_; i++)
95        values_[i] = CoinCopyOfArray(rhs.values_[i], numberIntegers_);
96    } else {
97      values_ = NULL;
98    }
99  }
100  return *this;
101}
102
103// Create C++ lines to get to current state
104void CbcHeuristicDINS::generateCpp(FILE *fp)
105{
106  CbcHeuristicDINS other;
107  fprintf(fp, "0#include \"CbcHeuristicDINS.hpp\"\n");
108  fprintf(fp, "3  CbcHeuristicDINS heuristicDINS(*cbcModel);\n");
109  CbcHeuristic::generateCpp(fp, "heuristicDINS");
110  if (howOften_ != other.howOften_)
111    fprintf(fp, "3  heuristicDINS.setHowOften(%d);\n", howOften_);
112  else
113    fprintf(fp, "4  heuristicDINS.setHowOften(%d);\n", howOften_);
114  if (maximumKeepSolutions_ != other.maximumKeepSolutions_)
115    fprintf(fp, "3  heuristicDINS.setMaximumKeep(%d);\n", maximumKeepSolutions_);
116  else
117    fprintf(fp, "4  heuristicDINS.setMaximumKeep(%d);\n", maximumKeepSolutions_);
118  fprintf(fp, "3  cbcModel->addHeuristic(&heuristicDINS);\n");
119}
120
121// Copy constructor
122CbcHeuristicDINS::CbcHeuristicDINS(const CbcHeuristicDINS &rhs)
123  : CbcHeuristic(rhs)
124  , numberSolutions_(rhs.numberSolutions_)
125  , howOften_(rhs.howOften_)
126  , numberSuccesses_(rhs.numberSuccesses_)
127  , numberTries_(rhs.numberTries_)
128  , maximumKeepSolutions_(rhs.maximumKeepSolutions_)
129  , numberKeptSolutions_(rhs.numberKeptSolutions_)
130  , numberIntegers_(rhs.numberIntegers_)
131  , localSpace_(rhs.localSpace_)
132{
133  if (model_ && rhs.values_) {
134    assert(numberIntegers_ >= 0);
135    values_ = new int *[maximumKeepSolutions_];
136    for (int i = 0; i < maximumKeepSolutions_; i++)
137      values_[i] = CoinCopyOfArray(rhs.values_[i], numberIntegers_);
138  } else {
139    values_ = NULL;
140  }
141}
142// Resets stuff if model changes
143void CbcHeuristicDINS::resetModel(CbcModel *)
144{
145  //CbcHeuristic::resetModel(model);
146  for (int i = 0; i < numberKeptSolutions_; i++)
147    delete[] values_[i];
148  delete[] values_;
149  numberKeptSolutions_ = 0;
150  numberIntegers_ = -1;
151  numberSolutions_ = 0;
152  values_ = NULL;
153}
154/*
155  First tries setting a variable to better value.  If feasible then
156  tries setting others.  If not feasible then tries swaps
157  Returns 1 if solution, 0 if not */
158int CbcHeuristicDINS::solution(double &solutionValue,
159  double *betterSolution)
160{
161  numCouldRun_++;
162  int returnCode = 0;
163  const double *bestSolution = model_->bestSolution();
164  if (!bestSolution)
165    return 0; // No solution found yet
166#ifdef HEURISTIC_INFORM
167  printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
168    heuristicName(), numRuns_, numCouldRun_, when_);
169#endif
170  if (numberSolutions_ < model_->getSolutionCount()) {
171    // new solution - add info
172    numberSolutions_ = model_->getSolutionCount();
173
174    int numberIntegers = model_->numberIntegers();
175    const int *integerVariable = model_->integerVariable();
176    if (numberIntegers_ < 0) {
177      numberIntegers_ = numberIntegers;
178      assert(!values_);
179      values_ = new int *[maximumKeepSolutions_];
180      for (int i = 0; i < maximumKeepSolutions_; i++)
181        values_[i] = NULL;
182    } else {
183      assert(numberIntegers >= numberIntegers_);
184    }
185    // move solutions (0 will be most recent)
186    {
187      int *temp = values_[maximumKeepSolutions_ - 1];
188      for (int i = maximumKeepSolutions_ - 1; i > 0; i--)
189        values_[i] = values_[i - 1];
190      if (!temp)
191        temp = new int[numberIntegers_];
192      values_[0] = temp;
193    }
194    int i;
195    for (i = 0; i < numberIntegers_; i++) {
196      int iColumn = integerVariable[i];
197      double value = bestSolution[iColumn];
198      double nearest = floor(value + 0.5);
199      values_[0][i] = static_cast< int >(nearest);
200    }
201    numberKeptSolutions_ = CoinMin(numberKeptSolutions_ + 1, maximumKeepSolutions_);
202  }
203  int finalReturnCode = 0;
204  if (((model_->getNodeCount() % howOften_) == howOften_ / 2 || !model_->getNodeCount()) && (model_->getCurrentPassNumber() <= 1 || model_->getCurrentPassNumber() == 999999)) {
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    int localSpace = localSpace_;
212    // 0 means finished but no solution, 1 solution, 2 node limit
213    int status = -1;
214    double cutoff = model_->getCutoff();
215    while (status) {
216      status = 0;
217      OsiSolverInterface *newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
218      const double *colLower = solver->getColLower();
219      const double *colUpper = solver->getColUpper();
220
221      double primalTolerance;
222      solver->getDblParam(OsiPrimalTolerance, primalTolerance);
223      const double *continuousSolution = newSolver->getColSolution();
224      // Space for added constraint
225      double *element = new double[numberIntegers_];
226      int *column = new int[numberIntegers_];
227      int i;
228      int nFix = 0;
229      int nCouldFix = 0;
230      int nCouldFix2 = 0;
231      int nBound = 0;
232      int nEl = 0;
233      double bias = localSpace;
234      int okSame = numberKeptSolutions_ - 1;
235      for (i = 0; i < numberIntegers_; i++) {
236        int iColumn = integerVariable[i];
237        const OsiObject *object = model_->object(i);
238        // get original bounds
239        double originalLower;
240        double originalUpper;
241        getIntegerInformation(object, originalLower, originalUpper);
242        double valueInt = bestSolution[iColumn];
243        if (valueInt < originalLower) {
244          valueInt = originalLower;
245        } else if (valueInt > originalUpper) {
246          valueInt = originalUpper;
247        }
248        int intValue = static_cast< int >(floor(valueInt + 0.5));
249        double currentValue = currentSolution[iColumn];
250        double currentLower = colLower[iColumn];
251        double currentUpper = colUpper[iColumn];
252        if (fabs(valueInt - currentValue) >= 0.5) {
253          // Re-bound
254          nBound++;
255          if (intValue >= currentValue) {
256            currentLower = CoinMax(currentLower, ceil(2 * currentValue - intValue));
257            currentUpper = intValue;
258          } else {
259            currentLower = intValue;
260            currentUpper = CoinMin(currentUpper, floor(2 * currentValue - intValue));
261          }
262          newSolver->setColLower(iColumn, currentLower);
263          newSolver->setColUpper(iColumn, currentUpper);
264        } else {
265          // See if can fix
266          bool canFix = false;
267          double continuousValue = continuousSolution[iColumn];
268          if (fabs(currentValue - valueInt) < 10.0 * primalTolerance) {
269            if (currentUpper - currentLower > 1.0) {
270              // General integer variable
271              canFix = true;
272            } else if (fabs(continuousValue - valueInt) < 10.0 * primalTolerance) {
273              int nSame = 1;
274              //assert (intValue==values_[0][i]);
275              for (int k = 1; k < numberKeptSolutions_; k++) {
276                if (intValue == values_[k][i])
277                  nSame++;
278              }
279              if (nSame >= okSame) {
280                // can fix
281                canFix = true;
282              } else {
283                nCouldFix++;
284              }
285            } else {
286              nCouldFix2++;
287            }
288          }
289          if (canFix) {
290            newSolver->setColLower(iColumn, intValue);
291            newSolver->setColUpper(iColumn, intValue);
292            nFix++;
293          } else {
294            if (currentUpper - currentLower > 1.0) {
295              // General integer variable
296              currentLower = floor(currentValue);
297              if (intValue >= currentLower && intValue <= currentLower + 1) {
298                newSolver->setColLower(iColumn, currentLower);
299                newSolver->setColUpper(iColumn, currentLower + 1.0);
300              } else {
301                // fix
302                double value;
303                if (intValue < currentLower)
304                  value = currentLower;
305                else
306                  value = currentLower + 1;
307                newSolver->setColLower(iColumn, value);
308                newSolver->setColUpper(iColumn, value);
309                nFix++;
310              }
311            } else {
312              // 0-1 (ish)
313              column[nEl] = iColumn;
314              if (intValue == currentLower) {
315                bias += currentLower;
316                element[nEl++] = 1.0;
317              } else if (intValue == currentUpper) {
318                bias += currentUpper;
319                element[nEl++] = -1.0;
320              } else {
321                printf("bad DINS logic\n");
322                abort();
323              }
324            }
325          }
326        }
327      }
328      char generalPrint[200];
329      sprintf(generalPrint,
330        "%d fixed, %d same as cont/int, %d same as int - %d bounded %d in cut\n",
331        nFix, nCouldFix, nCouldFix2, nBound, nEl);
332      model_->messageHandler()->message(CBC_FPUMP2, model_->messages())
333        << generalPrint
334        << CoinMessageEol;
335      if (nFix > numberIntegers_ / 10) {
336#ifdef JJF_ZERO
337        newSolver->initialSolve();
338        printf("obj %g\n", newSolver->getObjValue());
339        for (i = 0; i < numberIntegers_; i++) {
340          int iColumn = integerVariable[i];
341          printf("%d new bounds %g %g - solutions %g %g\n",
342            iColumn, newSolver->getColLower()[iColumn],
343            newSolver->getColUpper()[iColumn],
344            bestSolution[iColumn],
345            currentSolution[iColumn]);
346        }
347#endif
348        if (nEl > 0)
349          newSolver->addRow(nEl, column, element, -COIN_DBL_MAX, bias);
350        //printf("%d integers have same value\n",nFix);
351        returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
352          cutoff, "CbcHeuristicDINS");
353        if (returnCode < 0) {
354          returnCode = 0; // returned on size
355          status = 0;
356        } else {
357          numRuns_++;
358          if ((returnCode & 2) != 0) {
359            // could add cut as complete search
360            returnCode &= ~2;
361            if ((returnCode & 1) != 0) {
362              numberSuccesses_++;
363              status = 1;
364            } else {
365              // no solution
366              status = 0;
367            }
368          } else {
369            if ((returnCode & 1) != 0) {
370              numberSuccesses_++;
371              status = 1;
372            } else {
373              // no solution but node limit
374              status = 2;
375              if (nEl)
376                localSpace -= 5;
377              else
378                localSpace = -1;
379              if (localSpace < 0)
380                status = 0;
381            }
382          }
383          if ((returnCode & 1) != 0) {
384            cutoff = CoinMin(cutoff, solutionValue - model_->getCutoffIncrement());
385            finalReturnCode = 1;
386          }
387        }
388      }
389      delete[] element;
390      delete[] column;
391      delete newSolver;
392    }
393    numberTries_++;
394    if ((numberTries_ % 10) == 0 && numberSuccesses_ * 3 < numberTries_)
395      howOften_ += static_cast< int >(howOften_ * decayFactor_);
396  }
397  return finalReturnCode;
398}
399// update model
400void CbcHeuristicDINS::setModel(CbcModel *model)
401{
402  model_ = model;
403  // Get a copy of original matrix
404  assert(model_->solver());
405  for (int i = 0; i < numberKeptSolutions_; i++)
406    delete[] values_[i];
407  delete[] values_;
408  numberKeptSolutions_ = 0;
409  numberIntegers_ = -1;
410  numberSolutions_ = 0;
411  values_ = NULL;
412}
413
414/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
415*/
Note: See TracBrowser for help on using the repository browser.