source: branches/sandbox/Cbc/src/CbcHeuristicDINS.cpp @ 1389

Last change on this file since 1389 was 1368, checked in by EdwinStraver, 10 years ago

Broke out CbcHeuristicRENS, CbcHeuristicDINS and CbcHeuristicVND out of CbcHeuristicRINS.

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