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

Last change on this file since 1315 was 1315, checked in by forrest, 10 years ago

final changes before cleaning

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