source: branches/sandbox/Cbc/src/CbcHeuristicRINS.cpp @ 1286

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

Changed formatting using AStyle -A4 -p

File size: 42.7 KB
Line 
1/* $Id: CbcHeuristicRINS.cpp 1261 2009-10-30 12:45:20Z forrest $ */
2// Copyright (C) 2006, International Business Machines
3// Corporation and others.  All Rights Reserved.
4#if defined(_MSC_VER)
5// Turn off compiler warning about long names
6#  pragma warning(disable:4786)
7#endif
8#include <cassert>
9#include <cstdlib>
10#include <cmath>
11#include <cfloat>
12
13#include "OsiSolverInterface.hpp"
14#include "CbcModel.hpp"
15#include "CbcMessage.hpp"
16#include "CbcHeuristicRINS.hpp"
17#include "CbcBranchActual.hpp"
18#include "CbcStrategy.hpp"
19#include "CglPreProcess.hpp"
20
21// Default Constructor
22CbcHeuristicRINS::CbcHeuristicRINS()
23        : CbcHeuristic()
24{
25    numberSolutions_ = 0;
26    numberSuccesses_ = 0;
27    numberTries_ = 0;
28    stateOfFixing_ = 0;
29    lastNode_ = -999999;
30    howOften_ = 100;
31    decayFactor_ = 0.5;
32    used_ = NULL;
33    whereFrom_ = 1 + 8 + 255 * 256;
34}
35
36// Constructor with model - assumed before cuts
37
38CbcHeuristicRINS::CbcHeuristicRINS(CbcModel & model)
39        : CbcHeuristic(model)
40{
41    numberSolutions_ = 0;
42    numberSuccesses_ = 0;
43    numberTries_ = 0;
44    stateOfFixing_ = 0;
45    lastNode_ = -999999;
46    howOften_ = 100;
47    decayFactor_ = 0.5;
48    assert(model.solver());
49    int numberColumns = model.solver()->getNumCols();
50    used_ = new char[numberColumns];
51    memset(used_, 0, numberColumns);
52    whereFrom_ = 1 + 8 + 255 * 256;
53}
54
55// Destructor
56CbcHeuristicRINS::~CbcHeuristicRINS ()
57{
58    delete [] used_;
59}
60
61// Clone
62CbcHeuristic *
63CbcHeuristicRINS::clone() const
64{
65    return new CbcHeuristicRINS(*this);
66}
67
68// Assignment operator
69CbcHeuristicRINS &
70CbcHeuristicRINS::operator=( const CbcHeuristicRINS & 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        stateOfFixing_ = rhs.stateOfFixing_;
79        lastNode_ = rhs.lastNode_;
80        delete [] used_;
81        if (model_ && rhs.used_) {
82            int numberColumns = model_->solver()->getNumCols();
83            used_ = new char[numberColumns];
84            memcpy(used_, rhs.used_, numberColumns);
85        } else {
86            used_ = NULL;
87        }
88    }
89    return *this;
90}
91
92// Create C++ lines to get to current state
93void
94CbcHeuristicRINS::generateCpp( FILE * fp)
95{
96    CbcHeuristicRINS other;
97    fprintf(fp, "0#include \"CbcHeuristicRINS.hpp\"\n");
98    fprintf(fp, "3  CbcHeuristicRINS heuristicRINS(*cbcModel);\n");
99    CbcHeuristic::generateCpp(fp, "heuristicRINS");
100    if (howOften_ != other.howOften_)
101        fprintf(fp, "3  heuristicRINS.setHowOften(%d);\n", howOften_);
102    else
103        fprintf(fp, "4  heuristicRINS.setHowOften(%d);\n", howOften_);
104    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicRINS);\n");
105}
106
107// Copy constructor
108CbcHeuristicRINS::CbcHeuristicRINS(const CbcHeuristicRINS & rhs)
109        :
110        CbcHeuristic(rhs),
111        numberSolutions_(rhs.numberSolutions_),
112        howOften_(rhs.howOften_),
113        numberSuccesses_(rhs.numberSuccesses_),
114        numberTries_(rhs.numberTries_),
115        stateOfFixing_(rhs.stateOfFixing_),
116        lastNode_(rhs.lastNode_)
117{
118    if (model_ && rhs.used_) {
119        int numberColumns = model_->solver()->getNumCols();
120        used_ = new char[numberColumns];
121        memcpy(used_, rhs.used_, numberColumns);
122    } else {
123        used_ = NULL;
124    }
125}
126// Resets stuff if model changes
127void
128CbcHeuristicRINS::resetModel(CbcModel * /*model*/)
129{
130    //CbcHeuristic::resetModel(model);
131    delete [] used_;
132    stateOfFixing_ = 0;
133    if (model_ && used_) {
134        int numberColumns = model_->solver()->getNumCols();
135        used_ = new char[numberColumns];
136        memset(used_, 0, numberColumns);
137    } else {
138        used_ = NULL;
139    }
140}
141/*
142  First tries setting a variable to better value.  If feasible then
143  tries setting others.  If not feasible then tries swaps
144  Returns 1 if solution, 0 if not */
145int
146CbcHeuristicRINS::solution(double & solutionValue,
147                           double * betterSolution)
148{
149    numCouldRun_++;
150    int returnCode = 0;
151    const double * bestSolution = model_->bestSolution();
152    if (!bestSolution)
153        return 0; // No solution found yet
154    if (numberSolutions_ < model_->getSolutionCount()) {
155        // new solution - add info
156        numberSolutions_ = model_->getSolutionCount();
157
158        int numberIntegers = model_->numberIntegers();
159        const int * integerVariable = model_->integerVariable();
160
161        int i;
162        for (i = 0; i < numberIntegers; i++) {
163            int iColumn = integerVariable[i];
164            const OsiObject * object = model_->object(i);
165            // get original bounds
166            double originalLower;
167            double originalUpper;
168            getIntegerInformation( object, originalLower, originalUpper);
169            double value = bestSolution[iColumn];
170            if (value < originalLower) {
171                value = originalLower;
172            } else if (value > originalUpper) {
173                value = originalUpper;
174            }
175            double nearest = floor(value + 0.5);
176            // if away from lower bound mark that fact
177            if (nearest > originalLower) {
178                used_[iColumn] = 1;
179            }
180        }
181    }
182    int numberNodes = model_->getNodeCount();
183    if (howOften_ == 100) {
184        if (numberNodes < lastNode_ + 12)
185            return 0;
186        // Do at 50 and 100
187        if ((numberNodes > 40 && numberNodes <= 50) || (numberNodes > 90 && numberNodes < 100))
188            numberNodes = howOften_;
189    }
190    if ((numberNodes % howOften_) == 0 && (model_->getCurrentPassNumber() == 1 ||
191                                           model_->getCurrentPassNumber() == 999999)) {
192        lastNode_ = model_->getNodeCount();
193        OsiSolverInterface * solver = model_->solver();
194
195        int numberIntegers = model_->numberIntegers();
196        const int * integerVariable = model_->integerVariable();
197
198        const double * currentSolution = solver->getColSolution();
199        OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
200        int numberColumns = newSolver->getNumCols();
201        int numberContinuous = numberColumns - numberIntegers;
202
203        double primalTolerance;
204        solver->getDblParam(OsiPrimalTolerance, primalTolerance);
205
206        int i;
207        int nFix = 0;
208        for (i = 0; i < numberIntegers; i++) {
209            int iColumn = integerVariable[i];
210            const OsiObject * object = model_->object(i);
211            // get original bounds
212            double originalLower;
213            double originalUpper;
214            getIntegerInformation( object, originalLower, originalUpper);
215            double valueInt = bestSolution[iColumn];
216            if (valueInt < originalLower) {
217                valueInt = originalLower;
218            } else if (valueInt > originalUpper) {
219                valueInt = originalUpper;
220            }
221            if (fabs(currentSolution[iColumn] - valueInt) < 10.0*primalTolerance) {
222                double nearest = floor(valueInt + 0.5);
223                newSolver->setColLower(iColumn, nearest);
224                newSolver->setColUpper(iColumn, nearest);
225                nFix++;
226            }
227        }
228        int divisor = 0;
229        if (5*nFix > numberIntegers) {
230            if (numberContinuous > 2*numberIntegers && nFix*10 < numberColumns &&
231                    ((!numRuns_ && numberTries_ > 2) || stateOfFixing_)) {
232#define RINS_FIX_CONTINUOUS
233#ifdef RINS_FIX_CONTINUOUS
234                const double * colLower = newSolver->getColLower();
235                //const double * colUpper = newSolver->getColUpper();
236                int nAtLb = 0;
237                //double sumDj=0.0;
238                const double * dj = newSolver->getReducedCost();
239                double direction = newSolver->getObjSense();
240                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
241                    if (!newSolver->isInteger(iColumn)) {
242                        double value = bestSolution[iColumn];
243                        if (value < colLower[iColumn] + 1.0e-8) {
244                            //double djValue = dj[iColumn]*direction;
245                            nAtLb++;
246                            //sumDj += djValue;
247                        }
248                    }
249                }
250                if (nAtLb) {
251                    // fix some continuous
252                    double * sort = new double[nAtLb];
253                    int * which = new int [nAtLb];
254                    //double threshold = CoinMax((0.01*sumDj)/static_cast<double>(nAtLb),1.0e-6);
255                    int nFix2 = 0;
256                    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
257                        if (!newSolver->isInteger(iColumn)) {
258                            double value = bestSolution[iColumn];
259                            if (value < colLower[iColumn] + 1.0e-8) {
260                                double djValue = dj[iColumn] * direction;
261                                if (djValue > 1.0e-6) {
262                                    sort[nFix2] = -djValue;
263                                    which[nFix2++] = iColumn;
264                                }
265                            }
266                        }
267                    }
268                    CoinSort_2(sort, sort + nFix2, which);
269                    divisor = 4;
270                    if (stateOfFixing_ > 0)
271                        divisor = stateOfFixing_;
272                    else if (stateOfFixing_ < -1)
273                        divisor = (-stateOfFixing_) - 1;
274                    nFix2 = CoinMin(nFix2, (numberColumns - nFix) / divisor);
275                    for (int i = 0; i < nFix2; i++) {
276                        int iColumn = which[i];
277                        newSolver->setColUpper(iColumn, colLower[iColumn]);
278                    }
279                    delete [] sort;
280                    delete [] which;
281#ifdef CLP_INVESTIGATE2
282                    printf("%d integers have same value, and %d continuous fixed at lb\n",
283                           nFix, nFix2);
284#endif
285                }
286#endif
287            }
288            //printf("%d integers have same value\n",nFix);
289            returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
290                                             model_->getCutoff(), "CbcHeuristicRINS");
291            if (returnCode < 0) {
292                returnCode = 0; // returned on size
293                if (divisor)
294                    stateOfFixing_ = - divisor; // say failed
295            } else {
296                numRuns_++;
297                if (divisor)
298                    stateOfFixing_ =  divisor; // say small enough
299            }
300            if ((returnCode&1) != 0)
301                numberSuccesses_++;
302            //printf("return code %d",returnCode);
303            if ((returnCode&2) != 0) {
304                // could add cut
305                returnCode &= ~2;
306                //printf("could add cut with %d elements (if all 0-1)\n",nFix);
307            } else {
308                //printf("\n");
309            }
310        }
311
312        numberTries_++;
313        if ((numberTries_ % 10) == 0 && numberSuccesses_*3 < numberTries_)
314            howOften_ += static_cast<int> (howOften_ * decayFactor_);
315        delete newSolver;
316    }
317    return returnCode;
318}
319// update model
320void CbcHeuristicRINS::setModel(CbcModel * model)
321{
322    model_ = model;
323    // Get a copy of original matrix
324    assert(model_->solver());
325    delete [] used_;
326    int numberColumns = model->solver()->getNumCols();
327    used_ = new char[numberColumns];
328    memset(used_, 0, numberColumns);
329}
330// Default Constructor
331CbcHeuristicRENS::CbcHeuristicRENS()
332        : CbcHeuristic()
333{
334    numberTries_ = 0;
335    whereFrom_ = 256 + 1;
336}
337
338// Constructor with model - assumed before cuts
339
340CbcHeuristicRENS::CbcHeuristicRENS(CbcModel & model)
341        : CbcHeuristic(model)
342{
343    numberTries_ = 0;
344    whereFrom_ = 256 + 1;
345}
346
347// Destructor
348CbcHeuristicRENS::~CbcHeuristicRENS ()
349{
350}
351
352// Clone
353CbcHeuristic *
354CbcHeuristicRENS::clone() const
355{
356    return new CbcHeuristicRENS(*this);
357}
358
359// Assignment operator
360CbcHeuristicRENS &
361CbcHeuristicRENS::operator=( const CbcHeuristicRENS & rhs)
362{
363    if (this != &rhs) {
364        CbcHeuristic::operator=(rhs);
365        numberTries_ = rhs.numberTries_;
366    }
367    return *this;
368}
369
370// Copy constructor
371CbcHeuristicRENS::CbcHeuristicRENS(const CbcHeuristicRENS & rhs)
372        :
373        CbcHeuristic(rhs),
374        numberTries_(rhs.numberTries_)
375{
376}
377// Resets stuff if model changes
378void
379CbcHeuristicRENS::resetModel(CbcModel * )
380{
381}
382int
383CbcHeuristicRENS::solution(double & solutionValue,
384                           double * betterSolution)
385{
386    int returnCode = 0;
387    const double * bestSolution = model_->bestSolution();
388    if (numberTries_ || (when() < 2 && bestSolution))
389        return 0;
390    numberTries_++;
391    OsiSolverInterface * solver = model_->solver();
392
393    int numberIntegers = model_->numberIntegers();
394    const int * integerVariable = model_->integerVariable();
395
396    const double * currentSolution = solver->getColSolution();
397    OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
398    const double * colLower = newSolver->getColLower();
399    const double * colUpper = newSolver->getColUpper();
400
401    double primalTolerance;
402    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
403
404    int i;
405    int numberFixed = 0;
406    int numberTightened = 0;
407    int numberAtBound = 0;
408    int numberColumns = newSolver->getNumCols();
409    int numberContinuous = numberColumns - numberIntegers;
410
411    for (i = 0; i < numberIntegers; i++) {
412        int iColumn = integerVariable[i];
413        double value = currentSolution[iColumn];
414        double lower = colLower[iColumn];
415        double upper = colUpper[iColumn];
416        value = CoinMax(value, lower);
417        value = CoinMin(value, upper);
418#define RENS_FIX_ONLY_LOWER
419#ifndef RENS_FIX_ONLY_LOWER
420        if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
421            value = floor(value + 0.5);
422            if (value == lower || value == upper)
423                numberAtBound++;
424            newSolver->setColLower(iColumn, value);
425            newSolver->setColUpper(iColumn, value);
426            numberFixed++;
427        } else if (colUpper[iColumn] - colLower[iColumn] >= 2.0) {
428            numberTightened++;
429            newSolver->setColLower(iColumn, floor(value));
430            newSolver->setColUpper(iColumn, ceil(value));
431        }
432#else
433        if (fabs(value - floor(value + 0.5)) < 1.0e-8 &&
434                floor(value + 0.5) == lower) {
435            value = floor(value + 0.5);
436            numberAtBound++;
437            newSolver->setColLower(iColumn, value);
438            newSolver->setColUpper(iColumn, value);
439            numberFixed++;
440        } else if (colUpper[iColumn] - colLower[iColumn] >= 2.0) {
441            numberTightened++;
442            if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
443                value = floor(value + 0.5);
444                if (value < upper) {
445                    newSolver->setColLower(iColumn, CoinMax(value - 1.0, lower));
446                    newSolver->setColUpper(iColumn, CoinMin(value + 1.0, upper));
447                } else {
448                    newSolver->setColLower(iColumn, upper - 1.0);
449                }
450            } else {
451                newSolver->setColLower(iColumn, floor(value));
452                newSolver->setColUpper(iColumn, ceil(value));
453            }
454        }
455#endif
456    }
457    if (numberFixed > numberIntegers / 5) {
458        if (numberContinuous > numberIntegers && numberFixed < numberColumns / 5) {
459#define RENS_FIX_CONTINUOUS
460#ifdef RENS_FIX_CONTINUOUS
461            const double * colLower = newSolver->getColLower();
462            //const double * colUpper = newSolver->getColUpper();
463            int nAtLb = 0;
464            double sumDj = 0.0;
465            const double * dj = newSolver->getReducedCost();
466            double direction = newSolver->getObjSense();
467            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
468                if (!newSolver->isInteger(iColumn)) {
469                    double value = currentSolution[iColumn];
470                    if (value < colLower[iColumn] + 1.0e-8) {
471                        double djValue = dj[iColumn] * direction;
472                        nAtLb++;
473                        sumDj += djValue;
474                    }
475                }
476            }
477            if (nAtLb) {
478                // fix some continuous
479                double * sort = new double[nAtLb];
480                int * which = new int [nAtLb];
481                double threshold = CoinMax((0.01 * sumDj) / static_cast<double>(nAtLb), 1.0e-6);
482                int nFix2 = 0;
483                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
484                    if (!newSolver->isInteger(iColumn)) {
485                        double value = currentSolution[iColumn];
486                        if (value < colLower[iColumn] + 1.0e-8) {
487                            double djValue = dj[iColumn] * direction;
488                            if (djValue > threshold) {
489                                sort[nFix2] = -djValue;
490                                which[nFix2++] = iColumn;
491                            }
492                        }
493                    }
494                }
495                CoinSort_2(sort, sort + nFix2, which);
496                nFix2 = CoinMin(nFix2, (numberColumns - numberFixed) / 2);
497                for (int i = 0; i < nFix2; i++) {
498                    int iColumn = which[i];
499                    newSolver->setColUpper(iColumn, colLower[iColumn]);
500                }
501                delete [] sort;
502                delete [] which;
503#ifdef CLP_INVESTIGATE2
504                printf("%d integers fixed (%d tightened) (%d at bound), and %d continuous fixed at lb\n",
505                       numberFixed, numberTightened, numberAtBound, nFix2);
506#endif
507            }
508#endif
509        }
510#ifdef COIN_DEVELOP
511        printf("%d integers fixed and %d tightened\n", numberFixed, numberTightened);
512#endif
513        returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
514                                         model_->getCutoff(), "CbcHeuristicRENS");
515        if (returnCode < 0)
516            returnCode = 0; // returned on size
517        //printf("return code %d",returnCode);
518        if ((returnCode&2) != 0) {
519            // could add cut
520            returnCode &= ~2;
521#ifdef COIN_DEVELOP
522            if (!numberTightened && numberFixed == numberAtBound)
523                printf("could add cut with %d elements\n", numberFixed);
524#endif
525        } else {
526            //printf("\n");
527        }
528    }
529
530    delete newSolver;
531    return returnCode;
532}
533// update model
534void CbcHeuristicRENS::setModel(CbcModel * model)
535{
536    model_ = model;
537}
538
539// Default Constructor
540CbcHeuristicDINS::CbcHeuristicDINS()
541        : CbcHeuristic()
542{
543    numberSolutions_ = 0;
544    numberSuccesses_ = 0;
545    numberTries_ = 0;
546    howOften_ = 100;
547    decayFactor_ = 0.5;
548    maximumKeepSolutions_ = 5;
549    numberKeptSolutions_ = 0;
550    numberIntegers_ = -1;
551    localSpace_ = 10;
552    values_ = NULL;
553}
554
555// Constructor with model - assumed before cuts
556
557CbcHeuristicDINS::CbcHeuristicDINS(CbcModel & model)
558        : CbcHeuristic(model)
559{
560    numberSolutions_ = 0;
561    numberSuccesses_ = 0;
562    numberTries_ = 0;
563    howOften_ = 100;
564    decayFactor_ = 0.5;
565    assert(model.solver());
566    maximumKeepSolutions_ = 5;
567    numberKeptSolutions_ = 0;
568    numberIntegers_ = -1;
569    localSpace_ = 10;
570    values_ = NULL;
571}
572
573// Destructor
574CbcHeuristicDINS::~CbcHeuristicDINS ()
575{
576    for (int i = 0; i < numberKeptSolutions_; i++)
577        delete [] values_[i];
578    delete [] values_;
579}
580
581// Clone
582CbcHeuristic *
583CbcHeuristicDINS::clone() const
584{
585    return new CbcHeuristicDINS(*this);
586}
587
588// Assignment operator
589CbcHeuristicDINS &
590CbcHeuristicDINS::operator=( const CbcHeuristicDINS & rhs)
591{
592    if (this != &rhs) {
593        CbcHeuristic::operator=(rhs);
594        numberSolutions_ = rhs.numberSolutions_;
595        howOften_ = rhs.howOften_;
596        numberSuccesses_ = rhs.numberSuccesses_;
597        numberTries_ = rhs.numberTries_;
598        for (int i = 0; i < numberKeptSolutions_; i++)
599            delete [] values_[i];
600        delete [] values_;
601        maximumKeepSolutions_ = rhs.maximumKeepSolutions_;
602        numberKeptSolutions_ = rhs.numberKeptSolutions_;
603        numberIntegers_ = rhs.numberIntegers_;
604        localSpace_ = rhs.localSpace_;
605        if (model_ && rhs.values_) {
606            assert (numberIntegers_ >= 0);
607            values_ = new int * [maximumKeepSolutions_];
608            for (int i = 0; i < maximumKeepSolutions_; i++)
609                values_[i] = CoinCopyOfArray(rhs.values_[i], numberIntegers_);
610        } else {
611            values_ = NULL;
612        }
613    }
614    return *this;
615}
616
617// Create C++ lines to get to current state
618void
619CbcHeuristicDINS::generateCpp( FILE * fp)
620{
621    CbcHeuristicDINS other;
622    fprintf(fp, "0#include \"CbcHeuristicDINS.hpp\"\n");
623    fprintf(fp, "3  CbcHeuristicDINS heuristicDINS(*cbcModel);\n");
624    CbcHeuristic::generateCpp(fp, "heuristicDINS");
625    if (howOften_ != other.howOften_)
626        fprintf(fp, "3  heuristicDINS.setHowOften(%d);\n", howOften_);
627    else
628        fprintf(fp, "4  heuristicDINS.setHowOften(%d);\n", howOften_);
629    if (maximumKeepSolutions_ != other.maximumKeepSolutions_)
630        fprintf(fp, "3  heuristicDINS.setMaximumKeep(%d);\n", maximumKeepSolutions_);
631    else
632        fprintf(fp, "4  heuristicDINS.setMaximumKeep(%d);\n", maximumKeepSolutions_);
633    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicDINS);\n");
634}
635
636// Copy constructor
637CbcHeuristicDINS::CbcHeuristicDINS(const CbcHeuristicDINS & rhs)
638        :
639        CbcHeuristic(rhs),
640        numberSolutions_(rhs.numberSolutions_),
641        howOften_(rhs.howOften_),
642        numberSuccesses_(rhs.numberSuccesses_),
643        numberTries_(rhs.numberTries_),
644        maximumKeepSolutions_(rhs.maximumKeepSolutions_),
645        numberKeptSolutions_(rhs.numberKeptSolutions_),
646        numberIntegers_(rhs.numberIntegers_),
647        localSpace_(rhs.localSpace_)
648{
649    if (model_ && rhs.values_) {
650        assert (numberIntegers_ >= 0);
651        values_ = new int * [maximumKeepSolutions_];
652        for (int i = 0; i < maximumKeepSolutions_; i++)
653            values_[i] = CoinCopyOfArray(rhs.values_[i], numberIntegers_);
654    } else {
655        values_ = NULL;
656    }
657}
658// Resets stuff if model changes
659void
660CbcHeuristicDINS::resetModel(CbcModel * )
661{
662    //CbcHeuristic::resetModel(model);
663    for (int i = 0; i < numberKeptSolutions_; i++)
664        delete [] values_[i];
665    delete [] values_;
666    numberKeptSolutions_ = 0;
667    numberIntegers_ = -1;
668    numberSolutions_ = 0;
669    values_ = NULL;
670}
671/*
672  First tries setting a variable to better value.  If feasible then
673  tries setting others.  If not feasible then tries swaps
674  Returns 1 if solution, 0 if not */
675int
676CbcHeuristicDINS::solution(double & solutionValue,
677                           double * betterSolution)
678{
679    numCouldRun_++;
680    int returnCode = 0;
681    const double * bestSolution = model_->bestSolution();
682    if (!bestSolution)
683        return 0; // No solution found yet
684    if (numberSolutions_ < model_->getSolutionCount()) {
685        // new solution - add info
686        numberSolutions_ = model_->getSolutionCount();
687
688        int numberIntegers = model_->numberIntegers();
689        const int * integerVariable = model_->integerVariable();
690        if (numberIntegers_ < 0) {
691            numberIntegers_ = numberIntegers;
692            assert (!values_);
693            values_ = new int * [maximumKeepSolutions_];
694            for (int i = 0; i < maximumKeepSolutions_; i++)
695                values_[i] = NULL;
696        } else {
697            assert (numberIntegers == numberIntegers_);
698        }
699        // move solutions (0 will be most recent)
700        {
701            int * temp = values_[maximumKeepSolutions_-1];
702            for (int i = maximumKeepSolutions_ - 1; i > 0; i--)
703                values_[i] = values_[i-1];
704            if (!temp)
705                temp = new int [numberIntegers_];
706            values_[0] = temp;
707        }
708        int i;
709        for (i = 0; i < numberIntegers; i++) {
710            int iColumn = integerVariable[i];
711            double value = bestSolution[iColumn];
712            double nearest = floor(value + 0.5);
713            values_[0][i] = static_cast<int> (nearest);
714        }
715        numberKeptSolutions_ = CoinMin(numberKeptSolutions_ + 1, maximumKeepSolutions_);
716    }
717    int finalReturnCode = 0;
718    if (((model_->getNodeCount() % howOften_) == howOften_ / 2 || !model_->getNodeCount()) && (model_->getCurrentPassNumber() == 1 || model_->getCurrentPassNumber() == 999999)) {
719        OsiSolverInterface * solver = model_->solver();
720
721        int numberIntegers = model_->numberIntegers();
722        const int * integerVariable = model_->integerVariable();
723
724        const double * currentSolution = solver->getColSolution();
725        int localSpace = localSpace_;
726        // 0 means finished but no solution, 1 solution, 2 node limit
727        int status = -1;
728        double cutoff = model_->getCutoff();
729        while (status) {
730            status = 0;
731            OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
732            const double * colLower = solver->getColLower();
733            const double * colUpper = solver->getColUpper();
734
735            double primalTolerance;
736            solver->getDblParam(OsiPrimalTolerance, primalTolerance);
737            const double * continuousSolution = newSolver->getColSolution();
738            // Space for added constraint
739            double * element = new double [numberIntegers];
740            int * column = new int [numberIntegers];
741            int i;
742            int nFix = 0;
743            int nCouldFix = 0;
744            int nCouldFix2 = 0;
745            int nBound = 0;
746            int nEl = 0;
747            double bias = localSpace;
748            int okSame = numberKeptSolutions_ - 1;
749            for (i = 0; i < numberIntegers; i++) {
750                int iColumn = integerVariable[i];
751                const OsiObject * object = model_->object(i);
752                // get original bounds
753                double originalLower;
754                double originalUpper;
755                getIntegerInformation( object, originalLower, originalUpper);
756                double valueInt = bestSolution[iColumn];
757                if (valueInt < originalLower) {
758                    valueInt = originalLower;
759                } else if (valueInt > originalUpper) {
760                    valueInt = originalUpper;
761                }
762                int intValue = static_cast<int> (floor(valueInt + 0.5));
763                double currentValue = currentSolution[iColumn];
764                double currentLower = colLower[iColumn];
765                double currentUpper = colUpper[iColumn];
766                if (fabs(valueInt - currentValue) >= 0.5) {
767                    // Re-bound
768                    nBound++;
769                    if (intValue >= currentValue) {
770                        currentLower = CoinMax(currentLower, ceil(2 * currentValue - intValue));
771                        currentUpper = intValue;
772                    } else {
773                        currentLower = intValue;
774                        currentUpper = CoinMin(currentUpper, floor(2 * currentValue - intValue));
775                    }
776                    newSolver->setColLower(iColumn, currentLower);
777                    newSolver->setColUpper(iColumn, currentUpper);
778                } else {
779                    // See if can fix
780                    bool canFix = false;
781                    double continuousValue = continuousSolution[iColumn];
782                    if (fabs(currentValue - valueInt) < 10.0*primalTolerance) {
783                        if (currentUpper - currentLower > 1.0) {
784                            // General integer variable
785                            canFix = true;
786                        } else if (fabs(continuousValue - valueInt) < 10.0*primalTolerance) {
787                            int nSame = 1;
788                            //assert (intValue==values_[0][i]);
789                            for (int k = 1; k < numberKeptSolutions_; k++) {
790                                if (intValue == values_[k][i])
791                                    nSame++;
792                            }
793                            if (nSame >= okSame) {
794                                // can fix
795                                canFix = true;
796                            } else {
797                                nCouldFix++;
798                            }
799                        } else {
800                            nCouldFix2++;
801                        }
802                    }
803                    if (canFix) {
804                        newSolver->setColLower(iColumn, intValue);
805                        newSolver->setColUpper(iColumn, intValue);
806                        nFix++;
807                    } else {
808                        if (currentUpper - currentLower > 1.0) {
809                            // General integer variable
810                            currentLower = floor(currentValue);
811                            if (intValue >= currentLower && intValue <= currentLower + 1) {
812                                newSolver->setColLower(iColumn, currentLower);
813                                newSolver->setColUpper(iColumn, currentLower + 1.0);
814                            } else {
815                                // fix
816                                double value;
817                                if (intValue < currentLower)
818                                    value = currentLower;
819                                else
820                                    value = currentLower + 1;
821                                newSolver->setColLower(iColumn, value);
822                                newSolver->setColUpper(iColumn, value);
823                                nFix++;
824                            }
825                        } else {
826                            // 0-1 (ish)
827                            column[nEl] = iColumn;
828                            if (intValue == currentLower) {
829                                bias += currentLower;
830                                element[nEl++] = 1.0;
831                            } else if (intValue == currentUpper) {
832                                bias += currentUpper;
833                                element[nEl++] = -1.0;
834                            } else {
835                                printf("bad DINS logic\n");
836                                abort();
837                            }
838                        }
839                    }
840                }
841            }
842            char generalPrint[200];
843            sprintf(generalPrint,
844                    "%d fixed, %d same as cont/int, %d same as int - %d bounded %d in cut\n",
845                    nFix, nCouldFix, nCouldFix2, nBound, nEl);
846            model_->messageHandler()->message(CBC_FPUMP2, model_->messages())
847            << generalPrint
848            << CoinMessageEol;
849            if (nFix > numberIntegers / 10) {
850#if 0
851                newSolver->initialSolve();
852                printf("obj %g\n", newSolver->getObjValue());
853                for (i = 0; i < numberIntegers; i++) {
854                    int iColumn = integerVariable[i];
855                    printf("%d new bounds %g %g - solutions %g %g\n",
856                           iColumn, newSolver->getColLower()[iColumn],
857                           newSolver->getColUpper()[iColumn],
858                           bestSolution[iColumn],
859                           currentSolution[iColumn]);
860                }
861#endif
862                if (nEl > 0)
863                    newSolver->addRow(nEl, column, element, -COIN_DBL_MAX, bias);
864                //printf("%d integers have same value\n",nFix);
865                returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
866                                                 cutoff, "CbcHeuristicDINS");
867                if (returnCode < 0) {
868                    returnCode = 0; // returned on size
869                    status = 0;
870                } else {
871                    numRuns_++;
872                    if ((returnCode&2) != 0) {
873                        // could add cut as complete search
874                        returnCode &= ~2;
875                        if ((returnCode&1) != 0) {
876                            numberSuccesses_++;
877                            status = 1;
878                        } else {
879                            // no solution
880                            status = 0;
881                        }
882                    } else {
883                        if ((returnCode&1) != 0) {
884                            numberSuccesses_++;
885                            status = 1;
886                        } else {
887                            // no solution but node limit
888                            status = 2;
889                            if (nEl)
890                                localSpace -= 5;
891                            else
892                                localSpace = -1;
893                            if (localSpace < 0)
894                                status = 0;
895                        }
896                    }
897                    if ((returnCode&1) != 0) {
898                        cutoff = CoinMin(cutoff, solutionValue - model_->getCutoffIncrement());
899                        finalReturnCode = 1;
900                    }
901                }
902            }
903            delete [] element;
904            delete [] column;
905            delete newSolver;
906        }
907        numberTries_++;
908        if ((numberTries_ % 10) == 0 && numberSuccesses_*3 < numberTries_)
909            howOften_ += static_cast<int> (howOften_ * decayFactor_);
910    }
911    return finalReturnCode;
912}
913// update model
914void CbcHeuristicDINS::setModel(CbcModel * model)
915{
916    model_ = model;
917    // Get a copy of original matrix
918    assert(model_->solver());
919    for (int i = 0; i < numberKeptSolutions_; i++)
920        delete [] values_[i];
921    delete [] values_;
922    numberKeptSolutions_ = 0;
923    numberIntegers_ = -1;
924    numberSolutions_ = 0;
925    values_ = NULL;
926}
927
928// Default Constructor
929CbcHeuristicVND::CbcHeuristicVND()
930        : CbcHeuristic()
931{
932    numberSolutions_ = 0;
933    numberSuccesses_ = 0;
934    numberTries_ = 0;
935    lastNode_ = -999999;
936    howOften_ = 100;
937    decayFactor_ = 0.5;
938    baseSolution_ = NULL;
939    whereFrom_ = 1 + 8 + 255 * 256;
940    stepSize_ = 0;
941    k_ = 0;
942    kmax_ = 0;
943    nDifferent_ = 0;
944}
945
946// Constructor with model - assumed before cuts
947
948CbcHeuristicVND::CbcHeuristicVND(CbcModel & model)
949        : CbcHeuristic(model)
950{
951    numberSolutions_ = 0;
952    numberSuccesses_ = 0;
953    numberTries_ = 0;
954    lastNode_ = -999999;
955    howOften_ = 100;
956    decayFactor_ = 0.5;
957    assert(model.solver());
958    int numberColumns = model.solver()->getNumCols();
959    baseSolution_ = new double [numberColumns];
960    memset(baseSolution_, 0, numberColumns*sizeof(double));
961    whereFrom_ = 1 + 8 + 255 * 256;
962    stepSize_ = 0;
963    k_ = 0;
964    kmax_ = 0;
965    nDifferent_ = 0;
966}
967
968// Destructor
969CbcHeuristicVND::~CbcHeuristicVND ()
970{
971    delete [] baseSolution_;
972}
973
974// Clone
975CbcHeuristic *
976CbcHeuristicVND::clone() const
977{
978    return new CbcHeuristicVND(*this);
979}
980
981// Assignment operator
982CbcHeuristicVND &
983CbcHeuristicVND::operator=( const CbcHeuristicVND & rhs)
984{
985    if (this != &rhs) {
986        CbcHeuristic::operator=(rhs);
987        numberSolutions_ = rhs.numberSolutions_;
988        howOften_ = rhs.howOften_;
989        numberSuccesses_ = rhs.numberSuccesses_;
990        numberTries_ = rhs.numberTries_;
991        lastNode_ = rhs.lastNode_;
992        delete [] baseSolution_;
993        if (model_ && rhs.baseSolution_) {
994            int numberColumns = model_->solver()->getNumCols();
995            baseSolution_ = new double [numberColumns];
996            memcpy(baseSolution_, rhs.baseSolution_, numberColumns*sizeof(double));
997        } else {
998            baseSolution_ = NULL;
999        }
1000        stepSize_ = rhs.stepSize_;
1001        k_ = rhs.k_;
1002        kmax_ = rhs.kmax_;
1003        nDifferent_ = rhs.nDifferent_;
1004    }
1005    return *this;
1006}
1007
1008// Create C++ lines to get to current state
1009void
1010CbcHeuristicVND::generateCpp( FILE * fp)
1011{
1012    CbcHeuristicVND other;
1013    fprintf(fp, "0#include \"CbcHeuristicVND.hpp\"\n");
1014    fprintf(fp, "3  CbcHeuristicVND heuristicVND(*cbcModel);\n");
1015    CbcHeuristic::generateCpp(fp, "heuristicVND");
1016    if (howOften_ != other.howOften_)
1017        fprintf(fp, "3  heuristicVND.setHowOften(%d);\n", howOften_);
1018    else
1019        fprintf(fp, "4  heuristicVND.setHowOften(%d);\n", howOften_);
1020    fprintf(fp, "3  cbcModel->addHeuristic(&heuristicVND);\n");
1021}
1022
1023// Copy constructor
1024CbcHeuristicVND::CbcHeuristicVND(const CbcHeuristicVND & rhs)
1025        :
1026        CbcHeuristic(rhs),
1027        numberSolutions_(rhs.numberSolutions_),
1028        howOften_(rhs.howOften_),
1029        numberSuccesses_(rhs.numberSuccesses_),
1030        numberTries_(rhs.numberTries_),
1031        lastNode_(rhs.lastNode_)
1032{
1033    if (model_ && rhs.baseSolution_) {
1034        int numberColumns = model_->solver()->getNumCols();
1035        baseSolution_ = new double [numberColumns];
1036        memcpy(baseSolution_, rhs.baseSolution_, numberColumns*sizeof(double));
1037    } else {
1038        baseSolution_ = NULL;
1039    }
1040    stepSize_ = rhs.stepSize_;
1041    k_ = rhs.k_;
1042    kmax_ = rhs.kmax_;
1043    nDifferent_ = rhs.nDifferent_;
1044}
1045// Resets stuff if model changes
1046void
1047CbcHeuristicVND::resetModel(CbcModel * /*model*/)
1048{
1049    //CbcHeuristic::resetModel(model);
1050    delete [] baseSolution_;
1051    if (model_ && baseSolution_) {
1052        int numberColumns = model_->solver()->getNumCols();
1053        baseSolution_ = new double [numberColumns];
1054        memset(baseSolution_, 0, numberColumns*sizeof(double));
1055    } else {
1056        baseSolution_ = NULL;
1057    }
1058}
1059/*
1060  First tries setting a variable to better value.  If feasible then
1061  tries setting others.  If not feasible then tries swaps
1062  Returns 1 if solution, 0 if not */
1063int
1064CbcHeuristicVND::solution(double & solutionValue,
1065                          double * betterSolution)
1066{
1067    numCouldRun_++;
1068    int returnCode = 0;
1069    const double * bestSolution = model_->bestSolution();
1070    if (!bestSolution)
1071        return 0; // No solution found yet
1072    if (numberSolutions_ < model_->getSolutionCount()) {
1073        // new solution - add info
1074        numberSolutions_ = model_->getSolutionCount();
1075
1076        int numberIntegers = model_->numberIntegers();
1077        const int * integerVariable = model_->integerVariable();
1078
1079        int i;
1080        for (i = 0; i < numberIntegers; i++) {
1081            int iColumn = integerVariable[i];
1082            const OsiObject * object = model_->object(i);
1083            // get original bounds
1084            double originalLower;
1085            double originalUpper;
1086            getIntegerInformation( object, originalLower, originalUpper);
1087            double value = bestSolution[iColumn];
1088            if (value < originalLower) {
1089                value = originalLower;
1090            } else if (value > originalUpper) {
1091                value = originalUpper;
1092            }
1093        }
1094    }
1095    int numberNodes = model_->getNodeCount();
1096    if (howOften_ == 100) {
1097        if (numberNodes < lastNode_ + 12)
1098            return 0;
1099        // Do at 50 and 100
1100        if ((numberNodes > 40 && numberNodes <= 50) || (numberNodes > 90 && numberNodes < 100))
1101            numberNodes = howOften_;
1102    }
1103    if ((numberNodes % howOften_) == 0 && (model_->getCurrentPassNumber() == 1 ||
1104                                           model_->getCurrentPassNumber() == 999999)) {
1105        lastNode_ = model_->getNodeCount();
1106        OsiSolverInterface * solver = model_->solver();
1107
1108        int numberIntegers = model_->numberIntegers();
1109        const int * integerVariable = model_->integerVariable();
1110
1111        const double * currentSolution = solver->getColSolution();
1112        OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
1113        //const double * colLower = newSolver->getColLower();
1114        //const double * colUpper = newSolver->getColUpper();
1115
1116        double primalTolerance;
1117        solver->getDblParam(OsiPrimalTolerance, primalTolerance);
1118
1119        // Sort on distance
1120        double * distance = new double [numberIntegers];
1121        int * which = new int [numberIntegers];
1122
1123        int i;
1124        int nFix = 0;
1125        double tolerance = 10.0 * primalTolerance;
1126        for (i = 0; i < numberIntegers; i++) {
1127            int iColumn = integerVariable[i];
1128            const OsiObject * object = model_->object(i);
1129            // get original bounds
1130            double originalLower;
1131            double originalUpper;
1132            getIntegerInformation( object, originalLower, originalUpper);
1133            double valueInt = bestSolution[iColumn];
1134            if (valueInt < originalLower) {
1135                valueInt = originalLower;
1136            } else if (valueInt > originalUpper) {
1137                valueInt = originalUpper;
1138            }
1139            baseSolution_[iColumn] = currentSolution[iColumn];
1140            distance[i] = fabs(currentSolution[iColumn] - valueInt);
1141            which[i] = i;
1142            if (fabs(currentSolution[iColumn] - valueInt) < tolerance)
1143                nFix++;
1144        }
1145        CoinSort_2(distance, distance + numberIntegers, which);
1146        nDifferent_ = numberIntegers - nFix;
1147        stepSize_ = nDifferent_ / 10;
1148        k_ = stepSize_;
1149        //nFix = numberIntegers-stepSize_;
1150        for (i = 0; i < nFix; i++) {
1151            int j = which[i];
1152            int iColumn = integerVariable[j];
1153            const OsiObject * object = model_->object(i);
1154            // get original bounds
1155            double originalLower;
1156            double originalUpper;
1157            getIntegerInformation( object, originalLower, originalUpper);
1158            double valueInt = bestSolution[iColumn];
1159            if (valueInt < originalLower) {
1160                valueInt = originalLower;
1161            } else if (valueInt > originalUpper) {
1162                valueInt = originalUpper;
1163            }
1164            double nearest = floor(valueInt + 0.5);
1165            newSolver->setColLower(iColumn, nearest);
1166            newSolver->setColUpper(iColumn, nearest);
1167        }
1168        delete [] distance;
1169        delete [] which;
1170        if (nFix > numberIntegers / 5) {
1171            //printf("%d integers have samish value\n",nFix);
1172            returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
1173                                             model_->getCutoff(), "CbcHeuristicVND");
1174            if (returnCode < 0)
1175                returnCode = 0; // returned on size
1176            else
1177                numRuns_++;
1178            if ((returnCode&1) != 0)
1179                numberSuccesses_++;
1180            //printf("return code %d",returnCode);
1181            if ((returnCode&2) != 0) {
1182                // could add cut
1183                returnCode &= ~2;
1184                //printf("could add cut with %d elements (if all 0-1)\n",nFix);
1185            } else {
1186                //printf("\n");
1187            }
1188            numberTries_++;
1189            if ((numberTries_ % 10) == 0 && numberSuccesses_*3 < numberTries_)
1190                howOften_ += static_cast<int> (howOften_ * decayFactor_);
1191        }
1192
1193        delete newSolver;
1194    }
1195    return returnCode;
1196}
1197// update model
1198void CbcHeuristicVND::setModel(CbcModel * model)
1199{
1200    model_ = model;
1201    // Get a copy of original matrix
1202    assert(model_->solver());
1203    delete [] baseSolution_;
1204    int numberColumns = model->solver()->getNumCols();
1205    baseSolution_ = new double [numberColumns];
1206    memset(baseSolution_, 0, numberColumns*sizeof(double));
1207}
1208
1209
Note: See TracBrowser for help on using the repository browser.