source: stable/2.8/Cbc/src/CbcHeuristicDINS.cpp @ 1856

Last change on this file since 1856 was 1854, checked in by stefan, 7 years ago

fix svn keywords property

  • Property svn:keywords set to Author Date Id Revision
File size: 15.6 KB
Line 
1// $Id: CbcHeuristicDINS.cpp 1854 2013-01-28 00:02:55Z 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// 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    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
Note: See TracBrowser for help on using the repository browser.