source: trunk/Cbc/src/CbcHeuristicDINS.cpp @ 2094

Last change on this file since 2094 was 2094, checked in by forrest, 4 years ago

for memory leaks and heuristics and some experimental stuff

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