source: trunk/Cbc/src/CbcTreeLocal.cpp @ 2347

Last change on this file since 2347 was 2347, checked in by forrest, 2 years ago

get rid of divide by zero (and too much printing)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 65.3 KB
Line 
1/* $Id: CbcTreeLocal.cpp 2347 2017-11-09 17:41:54Z forrest $ */
2// Copyright (C) 2004, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include "CbcModel.hpp"
7#include "CbcNode.hpp"
8#include "CbcTreeLocal.hpp"
9#include "CoinPackedMatrix.hpp"
10#include "CoinTime.hpp"
11#include "OsiRowCutDebugger.hpp"
12#include <cassert>
13#ifdef JJF_ZERO
14// gdb doesn't always put breakpoints in this virtual function
15// just stick xxxxxx() where you want to start
16static void xxxxxx()
17{
18    printf("break\n");
19}
20#endif
21CbcTreeLocal::CbcTreeLocal()
22        : localNode_(NULL),
23        bestSolution_(NULL),
24        savedSolution_(NULL),
25        saveNumberSolutions_(0),
26        model_(NULL),
27        originalLower_(NULL),
28        originalUpper_(NULL),
29        range_(0),
30        typeCuts_(-1),
31        maxDiversification_(0),
32        diversification_(0),
33        nextStrong_(false),
34        rhs_(0.0),
35        savedGap_(0.0),
36        bestCutoff_(0.0),
37        timeLimit_(0),
38        startTime_(0),
39        nodeLimit_(0),
40        startNode_(-1),
41        searchType_(-1),
42        refine_(false)
43{
44
45}
46/* Constructor with solution.
47   range is upper bound on difference from given solution.
48   maxDiversification is maximum number of diversifications to try
49   timeLimit is seconds in subTree
50   nodeLimit is nodes in subTree
51*/
52CbcTreeLocal::CbcTreeLocal(CbcModel * model, const double * solution ,
53                           int range, int typeCuts, int maxDiversification,
54                           int timeLimit, int nodeLimit, bool refine)
55        : localNode_(NULL),
56        bestSolution_(NULL),
57        savedSolution_(NULL),
58        saveNumberSolutions_(0),
59        model_(model),
60        originalLower_(NULL),
61        originalUpper_(NULL),
62        range_(range),
63        typeCuts_(typeCuts),
64        maxDiversification_(maxDiversification),
65        diversification_(0),
66        nextStrong_(false),
67        rhs_(0.0),
68        savedGap_(0.0),
69        bestCutoff_(0.0),
70        timeLimit_(timeLimit),
71        startTime_(0),
72        nodeLimit_(nodeLimit),
73        startNode_(-1),
74        searchType_(-1),
75        refine_(refine)
76{
77
78    OsiSolverInterface * solver = model_->solver();
79    const double * lower = solver->getColLower();
80    const double * upper = solver->getColUpper();
81    //const double * solution = solver->getColSolution();
82    //const double * objective = solver->getObjCoefficients();
83    double primalTolerance;
84    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
85
86    // Get increment
87    model_->analyzeObjective();
88
89    {
90        // needed to sync cutoffs
91        double value ;
92        solver->getDblParam(OsiDualObjectiveLimit, value) ;
93        model_->setCutoff(value * solver->getObjSense());
94    }
95    bestCutoff_ = model_->getCutoff();
96    // save current gap
97    savedGap_ = model_->getDblParam(CbcModel::CbcAllowableGap);
98
99    // make sure integers found
100    model_->findIntegers(false);
101    int numberIntegers = model_->numberIntegers();
102    const int * integerVariable = model_->integerVariable();
103    int i;
104    double direction = solver->getObjSense();
105    double newSolutionValue = 1.0e50;
106    if (solution) {
107        // copy solution
108        solver->setColSolution(solution);
109        newSolutionValue = direction * solver->getObjValue();
110    }
111    originalLower_ = new double [numberIntegers];
112    originalUpper_ = new double [numberIntegers];
113    bool all01 = true;
114    int number01 = 0;
115    for (i = 0; i < numberIntegers; i++) {
116        int iColumn = integerVariable[i];
117        originalLower_[i] = lower[iColumn];
118        originalUpper_[i] = upper[iColumn];
119        if (upper[iColumn] - lower[iColumn] > 1.5)
120            all01 = false;
121        else if (upper[iColumn] - lower[iColumn] == 1.0)
122            number01++;
123    }
124    if (all01 && !typeCuts_)
125        typeCuts_ = 1; // may as well so we don't have to deal with refine
126    if (!number01 && !typeCuts_) {
127        if (model_->messageHandler()->logLevel() > 1)
128            printf("** No 0-1 variables and local search only on 0-1 - switching off\n");
129        typeCuts_ = -1;
130    } else {
131        if (model_->messageHandler()->logLevel() > 1) {
132            std::string type;
133            if (all01) {
134                printf("%d 0-1 variables normal local  cuts\n",
135                       number01);
136            } else if (typeCuts_) {
137                printf("%d 0-1 variables, %d other - general integer local cuts\n",
138                       number01, numberIntegers - number01);
139            } else {
140                printf("%d 0-1 variables, %d other - local cuts but just on 0-1 variables\n",
141                       number01, numberIntegers - number01);
142            }
143            printf("maximum diversifications %d, initial cutspace %d, max time %d seconds, max nodes %d\n",
144                   maxDiversification_, range_, timeLimit_, nodeLimit_);
145        }
146    }
147    int numberColumns = model_->getNumCols();
148    savedSolution_ = new double [numberColumns];
149    memset(savedSolution_, 0, numberColumns*sizeof(double));
150    if (solution) {
151        rhs_ = range_;
152        // Check feasible
153        int goodSolution = createCut(solution, cut_);
154        if (goodSolution >= 0) {
155            for (i = 0; i < numberIntegers; i++) {
156                int iColumn = integerVariable[i];
157                double value = floor(solution[iColumn] + 0.5);
158                // fix so setBestSolution will work
159                solver->setColLower(iColumn, value);
160                solver->setColUpper(iColumn, value);
161            }
162            model_->reserveCurrentSolution();
163            // Create cut and get total gap
164            if (newSolutionValue < bestCutoff_) {
165                model_->setBestSolution(CBC_ROUNDING, newSolutionValue, solution);
166                bestCutoff_ = model_->getCutoff();
167                // save as best solution
168                memcpy(savedSolution_, model_->bestSolution(), numberColumns*sizeof(double));
169            }
170            for (i = 0; i < numberIntegers; i++) {
171                int iColumn = integerVariable[i];
172                // restore bounds
173                solver->setColLower(iColumn, originalLower_[i]);
174                solver->setColUpper(iColumn, originalUpper_[i]);
175            }
176            // make sure can't stop on gap
177            model_->setDblParam(CbcModel::CbcAllowableGap, -1.0e50);
178        } else {
179            model_ = NULL;
180        }
181    } else {
182        // no solution
183        rhs_ = 1.0e50;
184        // make sure can't stop on gap
185        model_->setDblParam(CbcModel::CbcAllowableGap, -1.0e50);
186    }
187}
188CbcTreeLocal::~CbcTreeLocal()
189{
190    delete [] originalLower_;
191    delete [] originalUpper_;
192    delete [] bestSolution_;
193    delete [] savedSolution_;
194    delete localNode_;
195}
196// Copy constructor
197CbcTreeLocal::CbcTreeLocal ( const CbcTreeLocal & rhs)
198        : CbcTree(rhs),
199        saveNumberSolutions_(rhs.saveNumberSolutions_),
200        model_(rhs.model_),
201        range_(rhs.range_),
202        typeCuts_(rhs.typeCuts_),
203        maxDiversification_(rhs.maxDiversification_),
204        diversification_(rhs.diversification_),
205        nextStrong_(rhs.nextStrong_),
206        rhs_(rhs.rhs_),
207        savedGap_(rhs.savedGap_),
208        bestCutoff_(rhs.bestCutoff_),
209        timeLimit_(rhs.timeLimit_),
210        startTime_(rhs.startTime_),
211        nodeLimit_(rhs.nodeLimit_),
212        startNode_(rhs.startNode_),
213        searchType_(rhs.searchType_),
214        refine_(rhs.refine_)
215{
216    cut_ = rhs.cut_;
217    fixedCut_ = rhs.fixedCut_;
218    if (rhs.localNode_)
219        localNode_ = new CbcNode(*rhs.localNode_);
220    else
221        localNode_ = NULL;
222    if (rhs.originalLower_) {
223        int numberIntegers = model_->numberIntegers();
224        originalLower_ = new double [numberIntegers];
225        memcpy(originalLower_, rhs.originalLower_, numberIntegers*sizeof(double));
226        originalUpper_ = new double [numberIntegers];
227        memcpy(originalUpper_, rhs.originalUpper_, numberIntegers*sizeof(double));
228    } else {
229        originalLower_ = NULL;
230        originalUpper_ = NULL;
231    }
232    if (rhs.bestSolution_) {
233        int numberColumns = model_->getNumCols();
234        bestSolution_ = new double [numberColumns];
235        memcpy(bestSolution_, rhs.bestSolution_, numberColumns*sizeof(double));
236    } else {
237        bestSolution_ = NULL;
238    }
239    if (rhs.savedSolution_) {
240        int numberColumns = model_->getNumCols();
241        savedSolution_ = new double [numberColumns];
242        memcpy(savedSolution_, rhs.savedSolution_, numberColumns*sizeof(double));
243    } else {
244        savedSolution_ = NULL;
245    }
246}
247//----------------------------------------------------------------
248// Assignment operator
249//-------------------------------------------------------------------
250CbcTreeLocal &
251CbcTreeLocal::operator=(const CbcTreeLocal & rhs)
252{
253    if (this != &rhs) {
254        CbcTree::operator=(rhs);
255        saveNumberSolutions_ = rhs.saveNumberSolutions_;
256        cut_ = rhs.cut_;
257        fixedCut_ = rhs.fixedCut_;
258        delete localNode_;
259        if (rhs.localNode_)
260            localNode_ = new CbcNode(*rhs.localNode_);
261        else
262            localNode_ = NULL;
263        model_ = rhs.model_;
264        range_ = rhs.range_;
265        typeCuts_ = rhs.typeCuts_;
266        maxDiversification_ = rhs.maxDiversification_;
267        diversification_ = rhs.diversification_;
268        nextStrong_ = rhs.nextStrong_;
269        rhs_ = rhs.rhs_;
270        savedGap_ = rhs.savedGap_;
271        bestCutoff_ = rhs.bestCutoff_;
272        timeLimit_ = rhs.timeLimit_;
273        startTime_ = rhs.startTime_;
274        nodeLimit_ = rhs.nodeLimit_;
275        startNode_ = rhs.startNode_;
276        searchType_ = rhs.searchType_;
277        refine_ = rhs.refine_;
278        delete [] originalLower_;
279        delete [] originalUpper_;
280        if (rhs.originalLower_) {
281            int numberIntegers = model_->numberIntegers();
282            originalLower_ = new double [numberIntegers];
283            memcpy(originalLower_, rhs.originalLower_, numberIntegers*sizeof(double));
284            originalUpper_ = new double [numberIntegers];
285            memcpy(originalUpper_, rhs.originalUpper_, numberIntegers*sizeof(double));
286        } else {
287            originalLower_ = NULL;
288            originalUpper_ = NULL;
289        }
290        delete [] bestSolution_;
291        if (rhs.bestSolution_) {
292            int numberColumns = model_->getNumCols();
293            bestSolution_ = new double [numberColumns];
294            memcpy(bestSolution_, rhs.bestSolution_, numberColumns*sizeof(double));
295        } else {
296            bestSolution_ = NULL;
297        }
298        delete [] savedSolution_;
299        if (rhs.savedSolution_) {
300            int numberColumns = model_->getNumCols();
301            savedSolution_ = new double [numberColumns];
302            memcpy(savedSolution_, rhs.savedSolution_, numberColumns*sizeof(double));
303        } else {
304            savedSolution_ = NULL;
305        }
306    }
307    return *this;
308}
309// Clone
310CbcTree *
311CbcTreeLocal::clone() const
312{
313    return new CbcTreeLocal(*this);
314}
315// Pass in solution (so can be used after heuristic)
316void
317CbcTreeLocal::passInSolution(const double * solution, double solutionValue)
318{
319    int numberColumns = model_->getNumCols();
320    delete [] savedSolution_;
321    savedSolution_ = new double [numberColumns];
322    memcpy(savedSolution_, solution, numberColumns*sizeof(double));
323    rhs_ = range_;
324    // Check feasible
325    int goodSolution = createCut(solution, cut_);
326    if (goodSolution >= 0) {
327        bestCutoff_ = CoinMin(solutionValue, model_->getCutoff());
328    } else {
329        model_ = NULL;
330    }
331}
332// Return the top node of the heap
333CbcNode *
334CbcTreeLocal::top() const
335{
336#ifdef CBC_DEBUG
337    int smallest = 9999999;
338    int largest = -1;
339    double smallestD = 1.0e30;
340    double largestD = -1.0e30;
341    int n = nodes_.size();
342    for (int i = 0; i < n; i++) {
343        int nn = nodes_[i]->nodeInfo()->nodeNumber();
344        double dd = nodes_[i]->objectiveValue();
345        largest = CoinMax(largest, nn);
346        smallest = CoinMin(smallest, nn);
347        largestD = CoinMax(largestD, dd);
348        smallestD = CoinMin(smallestD, dd);
349    }
350    if (model_->messageHandler()->logLevel() > 1) {
351        printf("smallest %d, largest %d, top %d\n", smallest, largest,
352               nodes_.front()->nodeInfo()->nodeNumber());
353        printf("smallestD %g, largestD %g, top %g\n", smallestD, largestD, nodes_.front()->objectiveValue());
354    }
355#endif
356    return nodes_.front();
357}
358
359// Add a node to the heap
360void
361CbcTreeLocal::push(CbcNode * x)
362{
363    if (typeCuts_ >= 0 && !nodes_.size() && searchType_ < 0) {
364        startNode_ = model_->getNodeCount();
365        // save copy of node
366        localNode_ = new CbcNode(*x);
367
368        if (cut_.row().getNumElements()) {
369            // Add to global cuts
370            // we came in with solution
371            model_->makeGlobalCut(cut_);
372            if (model_->messageHandler()->logLevel() > 1)
373                printf("initial cut - rhs %g %g\n",
374                       cut_.lb(), cut_.ub());
375            searchType_ = 1;
376        } else {
377            // stop on first solution
378            searchType_ = 0;
379        }
380        startTime_ = static_cast<int> (CoinCpuTime());
381        saveNumberSolutions_ = model_->getSolutionCount();
382    }
383    nodes_.push_back(x);
384#ifdef CBC_DEBUG
385    if (model_->messageHandler()->logLevel() > 0)
386        printf("pushing node onto heap %d %x %x\n",
387               x->nodeInfo()->nodeNumber(), x, x->nodeInfo());
388#endif
389    std::push_heap(nodes_.begin(), nodes_.end(), comparison_);
390}
391
392// Remove the top node from the heap
393void
394CbcTreeLocal::pop()
395{
396    std::pop_heap(nodes_.begin(), nodes_.end(), comparison_);
397    nodes_.pop_back();
398}
399// Test if empty - does work if so
400bool
401CbcTreeLocal::empty()
402{
403    if (typeCuts_ < 0)
404        return !nodes_.size();
405    /* state -
406       0 iterating
407       1 subtree finished optimal solution for subtree found
408       2 subtree finished and no solution found
409       3 subtree exiting and solution found
410       4 subtree exiting and no solution found
411    */
412    int state = 0;
413    assert (searchType_ != 2);
414    if (searchType_) {
415        if (CoinCpuTime() - startTime_ > timeLimit_ || model_->getNodeCount() - startNode_ >= nodeLimit_) {
416            state = 4;
417        }
418    } else {
419        if (model_->getSolutionCount() > saveNumberSolutions_) {
420            state = 4;
421        }
422    }
423    if (!nodes_.size())
424        state = 2;
425    if (!state) {
426        return false;
427    }
428    // Finished this phase
429    int numberColumns = model_->getNumCols();
430    if (model_->getSolutionCount() > saveNumberSolutions_) {
431        if (model_->getCutoff() < bestCutoff_) {
432            // Save solution
433            if (!bestSolution_)
434                bestSolution_ = new double [numberColumns];
435            memcpy(bestSolution_, model_->bestSolution(), numberColumns*sizeof(double));
436            bestCutoff_ = model_->getCutoff();
437        }
438        state--;
439    }
440    // get rid of all nodes (safe even if already done)
441    double bestPossibleObjective;
442    cleanTree(model_, -COIN_DBL_MAX, bestPossibleObjective);
443
444    double increment = model_->getDblParam(CbcModel::CbcCutoffIncrement) ;
445    if (model_->messageHandler()->logLevel() > 1)
446        printf("local state %d after %d nodes and %d seconds, new solution %g, best solution %g, k was %g\n",
447               state,
448               model_->getNodeCount() - startNode_,
449               static_cast<int> (CoinCpuTime()) - startTime_,
450               model_->getCutoff() + increment, bestCutoff_ + increment, rhs_);
451    saveNumberSolutions_ = model_->getSolutionCount();
452    bool finished = false;
453    bool lastTry = false;
454    switch (state) {
455    case 1:
456        // solution found and subtree exhausted
457        if (rhs_ > 1.0e30) {
458            finished = true;
459        } else {
460            // find global cut and reverse
461            reverseCut(1);
462            searchType_ = 1; // first false
463            rhs_ = range_; // reset range
464            nextStrong_ = false;
465
466            // save best solution in this subtree
467            memcpy(savedSolution_, model_->bestSolution(), numberColumns*sizeof(double));
468        }
469        break;
470    case 2:
471        // solution not found and subtree exhausted
472        if (rhs_ > 1.0e30) {
473            finished = true;
474        } else {
475            // find global cut and reverse
476            reverseCut(2);
477            searchType_ = 1; // first false
478            if (diversification_ < maxDiversification_) {
479                if (nextStrong_) {
480                    diversification_++;
481                    // cut is valid so don't model_->setCutoff(1.0e50);
482                    searchType_ = 0;
483                }
484                nextStrong_ = true;
485                rhs_ += range_ / 2;
486            } else {
487                // This will be last try (may hit max time)
488                lastTry = true;
489                if (!maxDiversification_)
490                    typeCuts_ = -1; // make sure can't start again
491                model_->setCutoff(bestCutoff_);
492                if (model_->messageHandler()->logLevel() > 1)
493                    printf("Exiting local search with current set of cuts\n");
494                rhs_ = 1.0e100;
495                // Can now stop on gap
496                model_->setDblParam(CbcModel::CbcAllowableGap, savedGap_);
497            }
498        }
499        break;
500    case 3:
501        // solution found and subtree not exhausted
502        if (rhs_ < 1.0e30) {
503            if (searchType_) {
504                if (!typeCuts_ && refine_ && searchType_ == 1) {
505                    // We need to check we have best solution given these 0-1 values
506                    OsiSolverInterface * subSolver = model_->continuousSolver()->clone();
507                    CbcModel * subModel = model_->subTreeModel(subSolver);
508                    CbcTree normalTree;
509                    subModel->passInTreeHandler(normalTree);
510                    int numberIntegers = model_->numberIntegers();
511                    const int * integerVariable = model_->integerVariable();
512                    const double * solution = model_->bestSolution();
513                    int i;
514                    int numberColumns = model_->getNumCols();
515                    for (i = 0; i < numberIntegers; i++) {
516                        int iColumn = integerVariable[i];
517                        double value = floor(solution[iColumn] + 0.5);
518                        if (!typeCuts_ && originalUpper_[i] - originalLower_[i] > 1.0)
519                            continue; // skip as not 0-1
520                        if (originalLower_[i] == originalUpper_[i])
521                            continue;
522                        subSolver->setColLower(iColumn, value);
523                        subSolver->setColUpper(iColumn, value);
524                    }
525                    subSolver->initialSolve();
526                    // We can copy cutoff
527                    // But adjust
528                    subModel->setCutoff(model_->getCutoff() + model_->getDblParam(CbcModel::CbcCutoffIncrement) + 1.0e-6);
529                    subModel->setSolutionCount(0);
530                    assert (subModel->isProvenOptimal());
531                    if (!subModel->typePresolve()) {
532                        subModel->branchAndBound();
533                        if (subModel->status()) {
534                            model_->incrementSubTreeStopped();
535                        }
536                        //printf("%g %g %g %g\n",subModel->getCutoff(),model_->getCutoff(),
537                        //   subModel->getMinimizationObjValue(),model_->getMinimizationObjValue());
538                        double newCutoff = subModel->getMinimizationObjValue() -
539                                           subModel->getDblParam(CbcModel::CbcCutoffIncrement) ;
540                        if (subModel->getSolutionCount()) {
541                            if (!subModel->status())
542                                assert (subModel->isProvenOptimal());
543                            memcpy(model_->bestSolution(), subModel->bestSolution(),
544                                   numberColumns*sizeof(double));
545                            model_->setCutoff(newCutoff);
546                        }
547                    } else if (subModel->typePresolve() == 1) {
548                        CbcModel * model2 = subModel->integerPresolve(true);
549                        if (model2) {
550                            // Do complete search
551                            model2->branchAndBound();
552                            // get back solution
553                            subModel->originalModel(model2, false);
554                            if (model2->status()) {
555                                model_->incrementSubTreeStopped();
556                            }
557                            double newCutoff = model2->getMinimizationObjValue() -
558                                               model2->getDblParam(CbcModel::CbcCutoffIncrement) ;
559                            if (model2->getSolutionCount()) {
560                                if (!model2->status())
561                                    assert (model2->isProvenOptimal());
562                                memcpy(model_->bestSolution(), subModel->bestSolution(),
563                                       numberColumns*sizeof(double));
564                                model_->setCutoff(newCutoff);
565                            }
566                            delete model2;
567                        } else {
568                            // infeasible - could just be - due to cutoff
569                        }
570                    } else {
571                        // too dangerous at present
572                        assert (subModel->typePresolve() != 2);
573                    }
574                    if (model_->getCutoff() < bestCutoff_) {
575                        // Save solution
576                        if (!bestSolution_)
577                            bestSolution_ = new double [numberColumns];
578                        memcpy(bestSolution_, model_->bestSolution(), numberColumns*sizeof(double));
579                        bestCutoff_ = model_->getCutoff();
580                    }
581                    delete subModel;
582                }
583                // we have done search to make sure best general solution
584                searchType_ = 1;
585                // Reverse cut weakly
586                reverseCut(3, rhs_);
587            } else {
588                searchType_ = 1;
589                // delete last cut
590                deleteCut(cut_);
591            }
592        } else {
593            searchType_ = 1;
594        }
595        // save best solution in this subtree
596        memcpy(savedSolution_, model_->bestSolution(), numberColumns*sizeof(double));
597        nextStrong_ = false;
598        rhs_ = range_;
599        break;
600    case 4:
601        // solution not found and subtree not exhausted
602        if (maxDiversification_) {
603            if (nextStrong_) {
604                // Reverse cut weakly
605                reverseCut(4, rhs_);
606                model_->setCutoff(1.0e50);
607                diversification_++;
608                searchType_ = 0;
609            } else {
610                // delete last cut
611                deleteCut(cut_);
612                searchType_ = 1;
613            }
614            nextStrong_ = true;
615            rhs_ += range_ / 2;
616        } else {
617            // special case when using as heuristic
618            // Reverse cut weakly if lb -infinity
619            reverseCut(4, rhs_);
620            // This will be last try (may hit max time0
621            lastTry = true;
622            model_->setCutoff(bestCutoff_);
623            if (model_->messageHandler()->logLevel() > 1)
624                printf("Exiting local search with current set of cuts\n");
625            rhs_ = 1.0e100;
626            // Can now stop on gap
627            model_->setDblParam(CbcModel::CbcAllowableGap, savedGap_);
628            typeCuts_ = -1;
629        }
630        break;
631    }
632    if (rhs_ < 1.0e30 || lastTry) {
633        int goodSolution = createCut(savedSolution_, cut_);
634        if (goodSolution >= 0) {
635            // Add to global cuts
636            model_->makeGlobalCut(cut_);
637            CbcRowCuts * global = model_->globalCuts();
638            int n = global->sizeRowCuts();
639            OsiRowCut * rowCut = global->rowCutPtr(n - 1);
640            if (model_->messageHandler()->logLevel() > 1)
641                printf("inserting cut - now %d cuts, rhs %g %g, cutspace %g, diversification %d\n",
642                       n, rowCut->lb(), rowCut->ub(), rhs_, diversification_);
643            const OsiRowCutDebugger *debugger = model_->solver()->getRowCutDebuggerAlways() ;
644            if (debugger) {
645                if (debugger->invalidCut(*rowCut))
646                    printf("ZZZZTree Global cut - cuts off optimal solution!\n");
647            }
648            for (int i = 0; i < n; i++) {
649                rowCut = global->rowCutPtr(i);
650                if (model_->messageHandler()->logLevel() > 1)
651                    printf("%d - rhs %g %g\n",
652                           i, rowCut->lb(), rowCut->ub());
653            }
654        }
655        // put back node
656        startTime_ = static_cast<int> (CoinCpuTime());
657        startNode_ = model_->getNodeCount();
658        if (localNode_) {
659            // save copy of node
660            CbcNode * localNode2 = new CbcNode(*localNode_);
661            // But localNode2 now owns cuts so swap
662            //printf("pushing local node2 onto heap %d %x %x\n",localNode_->nodeNumber(),
663            //   localNode_,localNode_->nodeInfo());
664            nodes_.push_back(localNode_);
665            localNode_ = localNode2;
666            std::make_heap(nodes_.begin(), nodes_.end(), comparison_);
667        }
668    }
669    return finished;
670}
671// We may have got an intelligent tree so give it one more chance
672void
673CbcTreeLocal::endSearch()
674{
675    if (typeCuts_ >= 0) {
676        // copy best solution to model
677        int numberColumns = model_->getNumCols();
678        if (bestSolution_ && bestCutoff_ < model_->getCutoff()) {
679            memcpy(model_->bestSolution(), bestSolution_, numberColumns*sizeof(double));
680            model_->setCutoff(bestCutoff_);
681            // recompute objective value
682            const double * objCoef = model_->getObjCoefficients();
683            double objOffset = 0.0;
684            model_->continuousSolver()->getDblParam(OsiObjOffset, objOffset);
685
686            // Compute dot product of objCoef and colSol and then adjust by offset
687            double objValue = -objOffset;
688            for ( int i = 0 ; i < numberColumns ; i++ )
689                objValue += objCoef[i] * bestSolution_[i];
690            model_->setMinimizationObjValue(objValue);
691        }
692        // Can now stop on gap
693        model_->setDblParam(CbcModel::CbcAllowableGap, savedGap_);
694    }
695}
696// Create cut
697int
698CbcTreeLocal::createCut(const double * solution, OsiRowCut & rowCut)
699{
700    if (rhs_ > 1.0e20)
701        return -1;
702    OsiSolverInterface * solver = model_->solver();
703    const double * rowLower = solver->getRowLower();
704    const double * rowUpper = solver->getRowUpper();
705    //const double * solution = solver->getColSolution();
706    //const double * objective = solver->getObjCoefficients();
707    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
708    double primalTolerance;
709    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
710    // relax
711    primalTolerance *= 1000.0;
712
713    int numberRows = model_->getNumRows();
714
715    int numberIntegers = model_->numberIntegers();
716    const int * integerVariable = model_->integerVariable();
717    int i;
718
719    // Check feasible
720
721    double * rowActivity = new double[numberRows];
722    memset(rowActivity, 0, numberRows*sizeof(double));
723    solver->getMatrixByCol()->times(solution, rowActivity) ;
724    int goodSolution = 0;
725    // check was feasible
726    for (i = 0; i < numberRows; i++) {
727        if (rowActivity[i] < rowLower[i] - primalTolerance) {
728            goodSolution = -1;
729        } else if (rowActivity[i] > rowUpper[i] + primalTolerance) {
730            goodSolution = -1;
731        }
732    }
733    delete [] rowActivity;
734    for (i = 0; i < numberIntegers; i++) {
735        int iColumn = integerVariable[i];
736        double value = solution[iColumn];
737        if (fabs(floor(value + 0.5) - value) > integerTolerance) {
738            goodSolution = -1;
739        }
740    }
741    // zap cut
742    if (goodSolution == 0) {
743        // Create cut and get total gap
744        CoinPackedVector cut;
745        double rhs = rhs_;
746        double maxValue = 0.0;
747        for (i = 0; i < numberIntegers; i++) {
748            int iColumn = integerVariable[i];
749            double value = floor(solution[iColumn] + 0.5);
750            /*
751              typeCuts_ == 0 restricts to binary, 1 allows general integer. But we're
752              still restricted to being up against a bound. Consider: the notion is that
753              the cut restricts us to a k-neighbourhood. For binary variables, this
754              amounts to k variables which change value. For general integer, we could
755              end up with a single variable sucking up all of k (hence mu --- the
756              variable must swing to its other bound to look like a movement of 1).  For
757              variables in the middle of a range, we're talking about fabs(sol<j> - x<j>).
758            */
759            if (!typeCuts_ && originalUpper_[i] - originalLower_[i] > 1.0)
760                continue; // skip as not 0-1
761            if (originalLower_[i] == originalUpper_[i])
762                continue;
763            double mu = 1.0 / (originalUpper_[i] - originalLower_[i]);
764            if (value == originalLower_[i]) {
765                rhs += mu * originalLower_[i];
766                cut.insert(iColumn, 1.0);
767                maxValue += originalUpper_[i];
768            } else if (value == originalUpper_[i]) {
769                rhs -= mu * originalUpper_[i];
770                cut.insert(iColumn, -1.0);
771                maxValue += originalLower_[i];
772            }
773        }
774        if (maxValue < rhs - primalTolerance) {
775            if (model_->messageHandler()->logLevel() > 1)
776                printf("slack cut\n");
777            goodSolution = 1;
778        }
779        rowCut.setRow(cut);
780        rowCut.setLb(-COIN_DBL_MAX);
781        rowCut.setUb(rhs);
782        rowCut.setGloballyValid();
783        if (model_->messageHandler()->logLevel() > 1)
784            printf("Cut size: %i Cut rhs: %g\n", cut.getNumElements(), rhs);
785#ifdef CBC_DEBUG
786        if (model_->messageHandler()->logLevel() > 0) {
787            int k;
788            for (k = 0; k < cut.getNumElements(); k++) {
789                printf("%i    %g ", cut.getIndices()[k], cut.getElements()[k]);
790                if ((k + 1) % 5 == 0)
791                    printf("\n");
792            }
793            if (k % 5 != 0)
794                printf("\n");
795        }
796#endif
797        return goodSolution;
798    } else {
799        if (model_->messageHandler()->logLevel() > 1)
800            printf("Not a good solution\n");
801        return -1;
802    }
803}
804// Other side of last cut branch
805void
806CbcTreeLocal::reverseCut(int state, double bias)
807{
808    // find global cut
809    CbcRowCuts * global = model_->globalCuts();
810    int n = global->sizeRowCuts();
811    int i;
812    OsiRowCut * rowCut = NULL;
813    for ( i = 0; i < n; i++) {
814        rowCut = global->rowCutPtr(i);
815        if (cut_ == *rowCut) {
816            break;
817        }
818    }
819    if (!rowCut) {
820        // must have got here in odd way e.g. strong branching
821        return;
822    }
823    if (rowCut->lb() > -1.0e10)
824        return;
825    // get smallest element
826    double smallest = COIN_DBL_MAX;
827    CoinPackedVector row = cut_.row();
828    for (int k = 0; k < row.getNumElements(); k++)
829        smallest = CoinMin(smallest, fabs(row.getElements()[k]));
830    if (!typeCuts_ && !refine_) {
831        // Reverse cut very very weakly
832        if (state > 2)
833            smallest = 0.0;
834    }
835    // replace by other way
836    if (model_->messageHandler()->logLevel() > 1)
837        printf("reverseCut - changing cut %d out of %d, old rhs %g %g ",
838               i, n, rowCut->lb(), rowCut->ub());
839    rowCut->setLb(rowCut->ub() + smallest - bias);
840    rowCut->setUb(COIN_DBL_MAX);
841    if (model_->messageHandler()->logLevel() > 1)
842        printf("new rhs %g %g, bias %g smallest %g ",
843               rowCut->lb(), rowCut->ub(), bias, smallest);
844    const OsiRowCutDebugger *debugger = model_->solver()->getRowCutDebuggerAlways() ;
845    if (debugger) {
846        if (debugger->invalidCut(*rowCut))
847            printf("ZZZZTree Global cut - cuts off optimal solution!\n");
848    }
849}
850// Delete last cut branch
851void
852CbcTreeLocal::deleteCut(OsiRowCut & cut)
853{
854    // find global cut
855    CbcRowCuts * global = model_->globalCuts();
856    int n = global->sizeRowCuts();
857    int i;
858    OsiRowCut * rowCut = NULL;
859    for ( i = 0; i < n; i++) {
860        rowCut = global->rowCutPtr(i);
861        if (cut == *rowCut) {
862            break;
863        }
864    }
865    assert (i < n);
866    // delete last cut
867    if (model_->messageHandler()->logLevel() > 1)
868        printf("deleteCut - deleting cut %d out of %d, rhs %g %g\n",
869               i, n, rowCut->lb(), rowCut->ub());
870    global->eraseRowCut(i);
871}
872// Create C++ lines to get to current state
873void
874CbcTreeLocal::generateCpp( FILE * fp)
875{
876    CbcTreeLocal other;
877    fprintf(fp, "0#include \"CbcTreeLocal.hpp\"\n");
878    fprintf(fp, "5  CbcTreeLocal localTree(cbcModel,NULL);\n");
879    if (range_ != other.range_)
880        fprintf(fp, "5  localTree.setRange(%d);\n", range_);
881    if (typeCuts_ != other.typeCuts_)
882        fprintf(fp, "5  localTree.setTypeCuts(%d);\n", typeCuts_);
883    if (maxDiversification_ != other.maxDiversification_)
884        fprintf(fp, "5  localTree.setMaxDiversification(%d);\n", maxDiversification_);
885    if (timeLimit_ != other.timeLimit_)
886        fprintf(fp, "5  localTree.setTimeLimit(%d);\n", timeLimit_);
887    if (nodeLimit_ != other.nodeLimit_)
888        fprintf(fp, "5  localTree.setNodeLimit(%d);\n", nodeLimit_);
889    if (refine_ != other.refine_)
890        fprintf(fp, "5  localTree.setRefine(%s);\n", refine_ ? "true" : "false");
891    fprintf(fp, "5  cbcModel->passInTreeHandler(localTree);\n");
892}
893
894
895CbcTreeVariable::CbcTreeVariable()
896        : localNode_(NULL),
897        bestSolution_(NULL),
898        savedSolution_(NULL),
899        saveNumberSolutions_(0),
900        model_(NULL),
901        originalLower_(NULL),
902        originalUpper_(NULL),
903        range_(0),
904        typeCuts_(-1),
905        maxDiversification_(0),
906        diversification_(0),
907        nextStrong_(false),
908        rhs_(0.0),
909        savedGap_(0.0),
910        bestCutoff_(0.0),
911        timeLimit_(0),
912        startTime_(0),
913        nodeLimit_(0),
914        startNode_(-1),
915        searchType_(-1),
916        refine_(false)
917{
918
919}
920/* Constructor with solution.
921   range is upper bound on difference from given solution.
922   maxDiversification is maximum number of diversifications to try
923   timeLimit is seconds in subTree
924   nodeLimit is nodes in subTree
925*/
926CbcTreeVariable::CbcTreeVariable(CbcModel * model, const double * solution ,
927                                 int range, int typeCuts, int maxDiversification,
928                                 int timeLimit, int nodeLimit, bool refine)
929        : localNode_(NULL),
930        bestSolution_(NULL),
931        savedSolution_(NULL),
932        saveNumberSolutions_(0),
933        model_(model),
934        originalLower_(NULL),
935        originalUpper_(NULL),
936        range_(range),
937        typeCuts_(typeCuts),
938        maxDiversification_(maxDiversification),
939        diversification_(0),
940        nextStrong_(false),
941        rhs_(0.0),
942        savedGap_(0.0),
943        bestCutoff_(0.0),
944        timeLimit_(timeLimit),
945        startTime_(0),
946        nodeLimit_(nodeLimit),
947        startNode_(-1),
948        searchType_(-1),
949        refine_(refine)
950{
951
952    OsiSolverInterface * solver = model_->solver();
953    const double * lower = solver->getColLower();
954    const double * upper = solver->getColUpper();
955    //const double * solution = solver->getColSolution();
956    //const double * objective = solver->getObjCoefficients();
957    double primalTolerance;
958    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
959
960    // Get increment
961    model_->analyzeObjective();
962
963    {
964        // needed to sync cutoffs
965        double value ;
966        solver->getDblParam(OsiDualObjectiveLimit, value) ;
967        model_->setCutoff(value * solver->getObjSense());
968    }
969    bestCutoff_ = model_->getCutoff();
970    // save current gap
971    savedGap_ = model_->getDblParam(CbcModel::CbcAllowableGap);
972
973    // make sure integers found
974    model_->findIntegers(false);
975    int numberIntegers = model_->numberIntegers();
976    const int * integerVariable = model_->integerVariable();
977    int i;
978    double direction = solver->getObjSense();
979    double newSolutionValue = 1.0e50;
980    if (solution) {
981        // copy solution
982        solver->setColSolution(solution);
983        newSolutionValue = direction * solver->getObjValue();
984    }
985    originalLower_ = new double [numberIntegers];
986    originalUpper_ = new double [numberIntegers];
987    bool all01 = true;
988    int number01 = 0;
989    for (i = 0; i < numberIntegers; i++) {
990        int iColumn = integerVariable[i];
991        originalLower_[i] = lower[iColumn];
992        originalUpper_[i] = upper[iColumn];
993        if (upper[iColumn] - lower[iColumn] > 1.5)
994            all01 = false;
995        else if (upper[iColumn] - lower[iColumn] == 1.0)
996            number01++;
997    }
998    if (all01 && !typeCuts_)
999        typeCuts_ = 1; // may as well so we don't have to deal with refine
1000    if (!number01 && !typeCuts_) {
1001        if (model_->messageHandler()->logLevel() > 1)
1002            printf("** No 0-1 variables and local search only on 0-1 - switching off\n");
1003        typeCuts_ = -1;
1004    } else {
1005        if (model_->messageHandler()->logLevel() > 1) {
1006            std::string type;
1007            if (all01) {
1008                printf("%d 0-1 variables normal local  cuts\n",
1009                       number01);
1010            } else if (typeCuts_) {
1011                printf("%d 0-1 variables, %d other - general integer local cuts\n",
1012                       number01, numberIntegers - number01);
1013            } else {
1014                printf("%d 0-1 variables, %d other - local cuts but just on 0-1 variables\n",
1015                       number01, numberIntegers - number01);
1016            }
1017            printf("maximum diversifications %d, initial cutspace %d, max time %d seconds, max nodes %d\n",
1018                   maxDiversification_, range_, timeLimit_, nodeLimit_);
1019        }
1020    }
1021    int numberColumns = model_->getNumCols();
1022    savedSolution_ = new double [numberColumns];
1023    memset(savedSolution_, 0, numberColumns*sizeof(double));
1024    if (solution) {
1025        rhs_ = range_;
1026        // Check feasible
1027        int goodSolution = createCut(solution, cut_);
1028        if (goodSolution >= 0) {
1029            for (i = 0; i < numberIntegers; i++) {
1030                int iColumn = integerVariable[i];
1031                double value = floor(solution[iColumn] + 0.5);
1032                // fix so setBestSolution will work
1033                solver->setColLower(iColumn, value);
1034                solver->setColUpper(iColumn, value);
1035            }
1036            model_->reserveCurrentSolution();
1037            // Create cut and get total gap
1038            if (newSolutionValue < bestCutoff_) {
1039                model_->setBestSolution(CBC_ROUNDING, newSolutionValue, solution);
1040                bestCutoff_ = model_->getCutoff();
1041                // save as best solution
1042                memcpy(savedSolution_, model_->bestSolution(), numberColumns*sizeof(double));
1043            }
1044            for (i = 0; i < numberIntegers; i++) {
1045                int iColumn = integerVariable[i];
1046                // restore bounds
1047                solver->setColLower(iColumn, originalLower_[i]);
1048                solver->setColUpper(iColumn, originalUpper_[i]);
1049            }
1050            // make sure can't stop on gap
1051            model_->setDblParam(CbcModel::CbcAllowableGap, -1.0e50);
1052        } else {
1053            model_ = NULL;
1054        }
1055    } else {
1056        // no solution
1057        rhs_ = 1.0e50;
1058        // make sure can't stop on gap
1059        model_->setDblParam(CbcModel::CbcAllowableGap, -1.0e50);
1060    }
1061}
1062CbcTreeVariable::~CbcTreeVariable()
1063{
1064    delete [] originalLower_;
1065    delete [] originalUpper_;
1066    delete [] bestSolution_;
1067    delete [] savedSolution_;
1068    delete localNode_;
1069}
1070// Copy constructor
1071CbcTreeVariable::CbcTreeVariable ( const CbcTreeVariable & rhs)
1072        : CbcTree(rhs),
1073        saveNumberSolutions_(rhs.saveNumberSolutions_),
1074        model_(rhs.model_),
1075        range_(rhs.range_),
1076        typeCuts_(rhs.typeCuts_),
1077        maxDiversification_(rhs.maxDiversification_),
1078        diversification_(rhs.diversification_),
1079        nextStrong_(rhs.nextStrong_),
1080        rhs_(rhs.rhs_),
1081        savedGap_(rhs.savedGap_),
1082        bestCutoff_(rhs.bestCutoff_),
1083        timeLimit_(rhs.timeLimit_),
1084        startTime_(rhs.startTime_),
1085        nodeLimit_(rhs.nodeLimit_),
1086        startNode_(rhs.startNode_),
1087        searchType_(rhs.searchType_),
1088        refine_(rhs.refine_)
1089{
1090    cut_ = rhs.cut_;
1091    fixedCut_ = rhs.fixedCut_;
1092    if (rhs.localNode_)
1093        localNode_ = new CbcNode(*rhs.localNode_);
1094    else
1095        localNode_ = NULL;
1096    if (rhs.originalLower_) {
1097        int numberIntegers = model_->numberIntegers();
1098        originalLower_ = new double [numberIntegers];
1099        memcpy(originalLower_, rhs.originalLower_, numberIntegers*sizeof(double));
1100        originalUpper_ = new double [numberIntegers];
1101        memcpy(originalUpper_, rhs.originalUpper_, numberIntegers*sizeof(double));
1102    } else {
1103        originalLower_ = NULL;
1104        originalUpper_ = NULL;
1105    }
1106    if (rhs.bestSolution_) {
1107        int numberColumns = model_->getNumCols();
1108        bestSolution_ = new double [numberColumns];
1109        memcpy(bestSolution_, rhs.bestSolution_, numberColumns*sizeof(double));
1110    } else {
1111        bestSolution_ = NULL;
1112    }
1113    if (rhs.savedSolution_) {
1114        int numberColumns = model_->getNumCols();
1115        savedSolution_ = new double [numberColumns];
1116        memcpy(savedSolution_, rhs.savedSolution_, numberColumns*sizeof(double));
1117    } else {
1118        savedSolution_ = NULL;
1119    }
1120}
1121//----------------------------------------------------------------
1122// Assignment operator
1123//-------------------------------------------------------------------
1124CbcTreeVariable &
1125CbcTreeVariable::operator=(const CbcTreeVariable & rhs)
1126{
1127    if (this != &rhs) {
1128        CbcTree::operator=(rhs);
1129        saveNumberSolutions_ = rhs.saveNumberSolutions_;
1130        cut_ = rhs.cut_;
1131        fixedCut_ = rhs.fixedCut_;
1132        delete localNode_;
1133        if (rhs.localNode_)
1134            localNode_ = new CbcNode(*rhs.localNode_);
1135        else
1136            localNode_ = NULL;
1137        model_ = rhs.model_;
1138        range_ = rhs.range_;
1139        typeCuts_ = rhs.typeCuts_;
1140        maxDiversification_ = rhs.maxDiversification_;
1141        diversification_ = rhs.diversification_;
1142        nextStrong_ = rhs.nextStrong_;
1143        rhs_ = rhs.rhs_;
1144        savedGap_ = rhs.savedGap_;
1145        bestCutoff_ = rhs.bestCutoff_;
1146        timeLimit_ = rhs.timeLimit_;
1147        startTime_ = rhs.startTime_;
1148        nodeLimit_ = rhs.nodeLimit_;
1149        startNode_ = rhs.startNode_;
1150        searchType_ = rhs.searchType_;
1151        refine_ = rhs.refine_;
1152        delete [] originalLower_;
1153        delete [] originalUpper_;
1154        if (rhs.originalLower_) {
1155            int numberIntegers = model_->numberIntegers();
1156            originalLower_ = new double [numberIntegers];
1157            memcpy(originalLower_, rhs.originalLower_, numberIntegers*sizeof(double));
1158            originalUpper_ = new double [numberIntegers];
1159            memcpy(originalUpper_, rhs.originalUpper_, numberIntegers*sizeof(double));
1160        } else {
1161            originalLower_ = NULL;
1162            originalUpper_ = NULL;
1163        }
1164        delete [] bestSolution_;
1165        if (rhs.bestSolution_) {
1166            int numberColumns = model_->getNumCols();
1167            bestSolution_ = new double [numberColumns];
1168            memcpy(bestSolution_, rhs.bestSolution_, numberColumns*sizeof(double));
1169        } else {
1170            bestSolution_ = NULL;
1171        }
1172        delete [] savedSolution_;
1173        if (rhs.savedSolution_) {
1174            int numberColumns = model_->getNumCols();
1175            savedSolution_ = new double [numberColumns];
1176            memcpy(savedSolution_, rhs.savedSolution_, numberColumns*sizeof(double));
1177        } else {
1178            savedSolution_ = NULL;
1179        }
1180    }
1181    return *this;
1182}
1183// Clone
1184CbcTree *
1185CbcTreeVariable::clone() const
1186{
1187    return new CbcTreeVariable(*this);
1188}
1189// Pass in solution (so can be used after heuristic)
1190void
1191CbcTreeVariable::passInSolution(const double * solution, double solutionValue)
1192{
1193    int numberColumns = model_->getNumCols();
1194    delete [] savedSolution_;
1195    savedSolution_ = new double [numberColumns];
1196    memcpy(savedSolution_, solution, numberColumns*sizeof(double));
1197    rhs_ = range_;
1198    // Check feasible
1199    int goodSolution = createCut(solution, cut_);
1200    if (goodSolution >= 0) {
1201        bestCutoff_ = CoinMin(solutionValue, model_->getCutoff());
1202    } else {
1203        model_ = NULL;
1204    }
1205}
1206// Return the top node of the heap
1207CbcNode *
1208CbcTreeVariable::top() const
1209{
1210#ifdef CBC_DEBUG
1211    int smallest = 9999999;
1212    int largest = -1;
1213    double smallestD = 1.0e30;
1214    double largestD = -1.0e30;
1215    int n = nodes_.size();
1216    for (int i = 0; i < n; i++) {
1217        int nn = nodes_[i]->nodeInfo()->nodeNumber();
1218        double dd = nodes_[i]->objectiveValue();
1219        largest = CoinMax(largest, nn);
1220        smallest = CoinMin(smallest, nn);
1221        largestD = CoinMax(largestD, dd);
1222        smallestD = CoinMin(smallestD, dd);
1223    }
1224    if (model_->messageHandler()->logLevel() > 1) {
1225        printf("smallest %d, largest %d, top %d\n", smallest, largest,
1226               nodes_.front()->nodeInfo()->nodeNumber());
1227        printf("smallestD %g, largestD %g, top %g\n", smallestD, largestD, nodes_.front()->objectiveValue());
1228    }
1229#endif
1230    return nodes_.front();
1231}
1232
1233// Add a node to the heap
1234void
1235CbcTreeVariable::push(CbcNode * x)
1236{
1237    if (typeCuts_ >= 0 && !nodes_.size() && searchType_ < 0) {
1238        startNode_ = model_->getNodeCount();
1239        // save copy of node
1240        localNode_ = new CbcNode(*x);
1241
1242        if (cut_.row().getNumElements()) {
1243            // Add to global cuts
1244            // we came in with solution
1245            model_->makeGlobalCut(cut_);
1246            if (model_->messageHandler()->logLevel() > 1)
1247                printf("initial cut - rhs %g %g\n",
1248                       cut_.lb(), cut_.ub());
1249            searchType_ = 1;
1250        } else {
1251            // stop on first solution
1252            searchType_ = 0;
1253        }
1254        startTime_ = static_cast<int> (CoinCpuTime());
1255        saveNumberSolutions_ = model_->getSolutionCount();
1256    }
1257    nodes_.push_back(x);
1258#ifdef CBC_DEBUG
1259    if (model_->messageHandler()->logLevel() > 0)
1260        printf("pushing node onto heap %d %x %x\n",
1261               x->nodeInfo()->nodeNumber(), x, x->nodeInfo());
1262#endif
1263    std::push_heap(nodes_.begin(), nodes_.end(), comparison_);
1264}
1265
1266// Remove the top node from the heap
1267void
1268CbcTreeVariable::pop()
1269{
1270    std::pop_heap(nodes_.begin(), nodes_.end(), comparison_);
1271    nodes_.pop_back();
1272}
1273// Test if empty - does work if so
1274bool
1275CbcTreeVariable::empty()
1276{
1277    if (typeCuts_ < 0)
1278        return !nodes_.size();
1279    /* state -
1280       0 iterating
1281       1 subtree finished optimal solution for subtree found
1282       2 subtree finished and no solution found
1283       3 subtree exiting and solution found
1284       4 subtree exiting and no solution found
1285    */
1286    int state = 0;
1287    assert (searchType_ != 2);
1288    if (searchType_) {
1289        if (CoinCpuTime() - startTime_ > timeLimit_ || model_->getNodeCount() - startNode_ >= nodeLimit_) {
1290            state = 4;
1291        }
1292    } else {
1293        if (model_->getSolutionCount() > saveNumberSolutions_) {
1294            state = 4;
1295        }
1296    }
1297    if (!nodes_.size())
1298        state = 2;
1299    if (!state) {
1300        return false;
1301    }
1302    // Finished this phase
1303    int numberColumns = model_->getNumCols();
1304    if (model_->getSolutionCount() > saveNumberSolutions_) {
1305        if (model_->getCutoff() < bestCutoff_) {
1306            // Save solution
1307            if (!bestSolution_)
1308                bestSolution_ = new double [numberColumns];
1309            memcpy(bestSolution_, model_->bestSolution(), numberColumns*sizeof(double));
1310            bestCutoff_ = model_->getCutoff();
1311        }
1312        state--;
1313    }
1314    // get rid of all nodes (safe even if already done)
1315    double bestPossibleObjective;
1316    cleanTree(model_, -COIN_DBL_MAX, bestPossibleObjective);
1317
1318    double increment = model_->getDblParam(CbcModel::CbcCutoffIncrement) ;
1319    if (model_->messageHandler()->logLevel() > 1)
1320        printf("local state %d after %d nodes and %d seconds, new solution %g, best solution %g, k was %g\n",
1321               state,
1322               model_->getNodeCount() - startNode_,
1323               static_cast<int> (CoinCpuTime()) - startTime_,
1324               model_->getCutoff() + increment, bestCutoff_ + increment, rhs_);
1325    saveNumberSolutions_ = model_->getSolutionCount();
1326    bool finished = false;
1327    bool lastTry = false;
1328    switch (state) {
1329    case 1:
1330        // solution found and subtree exhausted
1331        if (rhs_ > 1.0e30) {
1332            finished = true;
1333        } else {
1334            // find global cut and reverse
1335            reverseCut(1);
1336            searchType_ = 1; // first false
1337            rhs_ = range_; // reset range
1338            nextStrong_ = false;
1339
1340            // save best solution in this subtree
1341            memcpy(savedSolution_, model_->bestSolution(), numberColumns*sizeof(double));
1342        }
1343        break;
1344    case 2:
1345        // solution not found and subtree exhausted
1346        if (rhs_ > 1.0e30) {
1347            finished = true;
1348        } else {
1349            // find global cut and reverse
1350            reverseCut(2);
1351            searchType_ = 1; // first false
1352            if (diversification_ < maxDiversification_) {
1353                if (nextStrong_) {
1354                    diversification_++;
1355                    // cut is valid so don't model_->setCutoff(1.0e50);
1356                    searchType_ = 0;
1357                }
1358                nextStrong_ = true;
1359                rhs_ += range_ / 2;
1360            } else {
1361                // This will be last try (may hit max time)
1362                lastTry = true;
1363                if (!maxDiversification_)
1364                    typeCuts_ = -1; // make sure can't start again
1365                model_->setCutoff(bestCutoff_);
1366                if (model_->messageHandler()->logLevel() > 1)
1367                    printf("Exiting local search with current set of cuts\n");
1368                rhs_ = 1.0e100;
1369                // Can now stop on gap
1370                model_->setDblParam(CbcModel::CbcAllowableGap, savedGap_);
1371            }
1372        }
1373        break;
1374    case 3:
1375        // solution found and subtree not exhausted
1376        if (rhs_ < 1.0e30) {
1377            if (searchType_) {
1378                if (!typeCuts_ && refine_ && searchType_ == 1) {
1379                    // We need to check we have best solution given these 0-1 values
1380                    OsiSolverInterface * subSolver = model_->continuousSolver()->clone();
1381                    CbcModel * subModel = model_->subTreeModel(subSolver);
1382                    CbcTree normalTree;
1383                    subModel->passInTreeHandler(normalTree);
1384                    int numberIntegers = model_->numberIntegers();
1385                    const int * integerVariable = model_->integerVariable();
1386                    const double * solution = model_->bestSolution();
1387                    int i;
1388                    int numberColumns = model_->getNumCols();
1389                    for (i = 0; i < numberIntegers; i++) {
1390                        int iColumn = integerVariable[i];
1391                        double value = floor(solution[iColumn] + 0.5);
1392                        if (!typeCuts_ && originalUpper_[i] - originalLower_[i] > 1.0)
1393                            continue; // skip as not 0-1
1394                        if (originalLower_[i] == originalUpper_[i])
1395                            continue;
1396                        subSolver->setColLower(iColumn, value);
1397                        subSolver->setColUpper(iColumn, value);
1398                    }
1399                    subSolver->initialSolve();
1400                    // We can copy cutoff
1401                    // But adjust
1402                    subModel->setCutoff(model_->getCutoff() + model_->getDblParam(CbcModel::CbcCutoffIncrement) + 1.0e-6);
1403                    subModel->setSolutionCount(0);
1404                    assert (subModel->isProvenOptimal());
1405                    if (!subModel->typePresolve()) {
1406                        subModel->branchAndBound();
1407                        if (subModel->status()) {
1408                            model_->incrementSubTreeStopped();
1409                        }
1410                        //printf("%g %g %g %g\n",subModel->getCutoff(),model_->getCutoff(),
1411                        //   subModel->getMinimizationObjValue(),model_->getMinimizationObjValue());
1412                        double newCutoff = subModel->getMinimizationObjValue() -
1413                                           subModel->getDblParam(CbcModel::CbcCutoffIncrement) ;
1414                        if (subModel->getSolutionCount()) {
1415                            if (!subModel->status())
1416                                assert (subModel->isProvenOptimal());
1417                            memcpy(model_->bestSolution(), subModel->bestSolution(),
1418                                   numberColumns*sizeof(double));
1419                            model_->setCutoff(newCutoff);
1420                        }
1421                    } else if (subModel->typePresolve() == 1) {
1422                        CbcModel * model2 = subModel->integerPresolve(true);
1423                        if (model2) {
1424                            // Do complete search
1425                            model2->branchAndBound();
1426                            // get back solution
1427                            subModel->originalModel(model2, false);
1428                            if (model2->status()) {
1429                                model_->incrementSubTreeStopped();
1430                            }
1431                            double newCutoff = model2->getMinimizationObjValue() -
1432                                               model2->getDblParam(CbcModel::CbcCutoffIncrement) ;
1433                            if (model2->getSolutionCount()) {
1434                                if (!model2->status())
1435                                    assert (model2->isProvenOptimal());
1436                                memcpy(model_->bestSolution(), subModel->bestSolution(),
1437                                       numberColumns*sizeof(double));
1438                                model_->setCutoff(newCutoff);
1439                            }
1440                            delete model2;
1441                        } else {
1442                            // infeasible - could just be - due to cutoff
1443                        }
1444                    } else {
1445                        // too dangerous at present
1446                        assert (subModel->typePresolve() != 2);
1447                    }
1448                    if (model_->getCutoff() < bestCutoff_) {
1449                        // Save solution
1450                        if (!bestSolution_)
1451                            bestSolution_ = new double [numberColumns];
1452                        memcpy(bestSolution_, model_->bestSolution(), numberColumns*sizeof(double));
1453                        bestCutoff_ = model_->getCutoff();
1454                    }
1455                    delete subModel;
1456                }
1457                // we have done search to make sure best general solution
1458                searchType_ = 1;
1459                // Reverse cut weakly
1460                reverseCut(3, rhs_);
1461            } else {
1462                searchType_ = 1;
1463                // delete last cut
1464                deleteCut(cut_);
1465            }
1466        } else {
1467            searchType_ = 1;
1468        }
1469        // save best solution in this subtree
1470        memcpy(savedSolution_, model_->bestSolution(), numberColumns*sizeof(double));
1471        nextStrong_ = false;
1472        rhs_ = range_;
1473        break;
1474    case 4:
1475        // solution not found and subtree not exhausted
1476        if (maxDiversification_) {
1477            if (nextStrong_) {
1478                // Reverse cut weakly
1479                reverseCut(4, rhs_);
1480                model_->setCutoff(1.0e50);
1481                diversification_++;
1482                searchType_ = 0;
1483            } else {
1484                // delete last cut
1485                deleteCut(cut_);
1486                searchType_ = 1;
1487            }
1488            nextStrong_ = true;
1489            rhs_ += range_ / 2;
1490        } else {
1491            // special case when using as heuristic
1492            // Reverse cut weakly if lb -infinity
1493            reverseCut(4, rhs_);
1494            // This will be last try (may hit max time0
1495            lastTry = true;
1496            model_->setCutoff(bestCutoff_);
1497            if (model_->messageHandler()->logLevel() > 1)
1498                printf("Exiting local search with current set of cuts\n");
1499            rhs_ = 1.0e100;
1500            // Can now stop on gap
1501            model_->setDblParam(CbcModel::CbcAllowableGap, savedGap_);
1502            typeCuts_ = -1;
1503        }
1504        break;
1505    }
1506    if (rhs_ < 1.0e30 || lastTry) {
1507        int goodSolution = createCut(savedSolution_, cut_);
1508        if (goodSolution >= 0) {
1509            // Add to global cuts
1510            model_->makeGlobalCut(cut_);
1511            CbcRowCuts * global = model_->globalCuts();
1512            int n = global->sizeRowCuts();
1513            OsiRowCut * rowCut = global->rowCutPtr(n - 1);
1514            if (model_->messageHandler()->logLevel() > 1)
1515                printf("inserting cut - now %d cuts, rhs %g %g, cutspace %g, diversification %d\n",
1516                       n, rowCut->lb(), rowCut->ub(), rhs_, diversification_);
1517            const OsiRowCutDebugger *debugger = model_->solver()->getRowCutDebuggerAlways() ;
1518            if (debugger) {
1519                if (debugger->invalidCut(*rowCut))
1520                    printf("ZZZZTree Global cut - cuts off optimal solution!\n");
1521            }
1522            for (int i = 0; i < n; i++) {
1523                rowCut = global->rowCutPtr(i);
1524                if (model_->messageHandler()->logLevel() > 1)
1525                    printf("%d - rhs %g %g\n",
1526                           i, rowCut->lb(), rowCut->ub());
1527            }
1528        }
1529        // put back node
1530        startTime_ = static_cast<int> (CoinCpuTime());
1531        startNode_ = model_->getNodeCount();
1532        if (localNode_) {
1533            // save copy of node
1534            CbcNode * localNode2 = new CbcNode(*localNode_);
1535            // But localNode2 now owns cuts so swap
1536            //printf("pushing local node2 onto heap %d %x %x\n",localNode_->nodeNumber(),
1537            //   localNode_,localNode_->nodeInfo());
1538            nodes_.push_back(localNode_);
1539            localNode_ = localNode2;
1540            std::make_heap(nodes_.begin(), nodes_.end(), comparison_);
1541        }
1542    }
1543    return finished;
1544}
1545// We may have got an intelligent tree so give it one more chance
1546void
1547CbcTreeVariable::endSearch()
1548{
1549    if (typeCuts_ >= 0) {
1550        // copy best solution to model
1551        int numberColumns = model_->getNumCols();
1552        if (bestSolution_ && bestCutoff_ < model_->getCutoff()) {
1553            memcpy(model_->bestSolution(), bestSolution_, numberColumns*sizeof(double));
1554            model_->setCutoff(bestCutoff_);
1555            // recompute objective value
1556            const double * objCoef = model_->getObjCoefficients();
1557            double objOffset = 0.0;
1558            model_->continuousSolver()->getDblParam(OsiObjOffset, objOffset);
1559
1560            // Compute dot product of objCoef and colSol and then adjust by offset
1561            double objValue = -objOffset;
1562            for ( int i = 0 ; i < numberColumns ; i++ )
1563                objValue += objCoef[i] * bestSolution_[i];
1564            model_->setMinimizationObjValue(objValue);
1565        }
1566        // Can now stop on gap
1567        model_->setDblParam(CbcModel::CbcAllowableGap, savedGap_);
1568    }
1569}
1570// Create cut
1571int
1572CbcTreeVariable::createCut(const double * solution, OsiRowCut & rowCut)
1573{
1574    if (rhs_ > 1.0e20)
1575        return -1;
1576    OsiSolverInterface * solver = model_->solver();
1577    const double * rowLower = solver->getRowLower();
1578    const double * rowUpper = solver->getRowUpper();
1579    //const double * solution = solver->getColSolution();
1580    //const double * objective = solver->getObjCoefficients();
1581    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
1582    double primalTolerance;
1583    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
1584    // relax
1585    primalTolerance *= 1000.0;
1586
1587    int numberRows = model_->getNumRows();
1588
1589    int numberIntegers = model_->numberIntegers();
1590    const int * integerVariable = model_->integerVariable();
1591    int i;
1592
1593    // Check feasible
1594
1595    double * rowActivity = new double[numberRows];
1596    memset(rowActivity, 0, numberRows*sizeof(double));
1597    solver->getMatrixByCol()->times(solution, rowActivity) ;
1598    int goodSolution = 0;
1599    // check was feasible
1600    for (i = 0; i < numberRows; i++) {
1601        if (rowActivity[i] < rowLower[i] - primalTolerance) {
1602            goodSolution = -1;
1603        } else if (rowActivity[i] > rowUpper[i] + primalTolerance) {
1604            goodSolution = -1;
1605        }
1606    }
1607    delete [] rowActivity;
1608    for (i = 0; i < numberIntegers; i++) {
1609        int iColumn = integerVariable[i];
1610        double value = solution[iColumn];
1611        if (fabs(floor(value + 0.5) - value) > integerTolerance) {
1612            goodSolution = -1;
1613        }
1614    }
1615    // zap cut
1616    if (goodSolution == 0) {
1617        // Create cut and get total gap
1618        CoinPackedVector cut;
1619        double rhs = rhs_;
1620        double maxValue = 0.0;
1621        for (i = 0; i < numberIntegers; i++) {
1622            int iColumn = integerVariable[i];
1623            double value = floor(solution[iColumn] + 0.5);
1624            if (!typeCuts_ && originalUpper_[i] - originalLower_[i] > 1.0)
1625                continue; // skip as not 0-1
1626            if (originalLower_[i] == originalUpper_[i])
1627                continue;
1628            double mu = 1.0 / (originalUpper_[i] - originalLower_[i]);
1629            if (value == originalLower_[i]) {
1630                rhs += mu * originalLower_[i];
1631                cut.insert(iColumn, 1.0);
1632                maxValue += originalUpper_[i];
1633            } else if (value == originalUpper_[i]) {
1634                rhs -= mu * originalUpper_[i];
1635                cut.insert(iColumn, -1.0);
1636                maxValue += originalLower_[i];
1637            }
1638        }
1639        if (maxValue < rhs - primalTolerance) {
1640            if (model_->messageHandler()->logLevel() > 1)
1641                printf("slack cut\n");
1642            goodSolution = 1;
1643        }
1644        rowCut.setRow(cut);
1645        rowCut.setLb(-COIN_DBL_MAX);
1646        rowCut.setUb(rhs);
1647        rowCut.setGloballyValid();
1648        if (model_->messageHandler()->logLevel() > 1)
1649            printf("Cut size: %i Cut rhs: %g\n", cut.getNumElements(), rhs);
1650#ifdef CBC_DEBUG
1651        if (model_->messageHandler()->logLevel() > 0) {
1652            int k;
1653            for (k = 0; k < cut.getNumElements(); k++) {
1654                printf("%i    %g ", cut.getIndices()[k], cut.getElements()[k]);
1655                if ((k + 1) % 5 == 0)
1656                    printf("\n");
1657            }
1658            if (k % 5 != 0)
1659                printf("\n");
1660        }
1661#endif
1662        return goodSolution;
1663    } else {
1664        if (model_->messageHandler()->logLevel() > 1)
1665            printf("Not a good solution\n");
1666        return -1;
1667    }
1668}
1669// Other side of last cut branch
1670void
1671CbcTreeVariable::reverseCut(int state, double bias)
1672{
1673    // find global cut
1674    CbcRowCuts * global = model_->globalCuts();
1675    int n = global->sizeRowCuts();
1676    int i;
1677    OsiRowCut * rowCut = NULL;
1678    for ( i = 0; i < n; i++) {
1679        rowCut = global->rowCutPtr(i);
1680        if (cut_ == *rowCut) {
1681            break;
1682        }
1683    }
1684    if (!rowCut) {
1685        // must have got here in odd way e.g. strong branching
1686        return;
1687    }
1688    if (rowCut->lb() > -1.0e10)
1689        return;
1690    // get smallest element
1691    double smallest = COIN_DBL_MAX;
1692    CoinPackedVector row = cut_.row();
1693    for (int k = 0; k < row.getNumElements(); k++)
1694        smallest = CoinMin(smallest, fabs(row.getElements()[k]));
1695    if (!typeCuts_ && !refine_) {
1696        // Reverse cut very very weakly
1697        if (state > 2)
1698            smallest = 0.0;
1699    }
1700    // replace by other way
1701    if (model_->messageHandler()->logLevel() > 1)
1702        printf("reverseCut - changing cut %d out of %d, old rhs %g %g ",
1703               i, n, rowCut->lb(), rowCut->ub());
1704    rowCut->setLb(rowCut->ub() + smallest - bias);
1705    rowCut->setUb(COIN_DBL_MAX);
1706    if (model_->messageHandler()->logLevel() > 1)
1707        printf("new rhs %g %g, bias %g smallest %g ",
1708               rowCut->lb(), rowCut->ub(), bias, smallest);
1709    const OsiRowCutDebugger *debugger = model_->solver()->getRowCutDebuggerAlways() ;
1710    if (debugger) {
1711        if (debugger->invalidCut(*rowCut))
1712            printf("ZZZZTree Global cut - cuts off optimal solution!\n");
1713    }
1714}
1715// Delete last cut branch
1716void
1717CbcTreeVariable::deleteCut(OsiRowCut & cut)
1718{
1719    // find global cut
1720    CbcRowCuts * global = model_->globalCuts();
1721    int n = global->sizeRowCuts();
1722    int i;
1723    OsiRowCut * rowCut = NULL;
1724    for ( i = 0; i < n; i++) {
1725        rowCut = global->rowCutPtr(i);
1726        if (cut == *rowCut) {
1727            break;
1728        }
1729    }
1730    assert (i < n);
1731    // delete last cut
1732    if (model_->messageHandler()->logLevel() > 1)
1733        printf("deleteCut - deleting cut %d out of %d, rhs %g %g\n",
1734               i, n, rowCut->lb(), rowCut->ub());
1735    global->eraseRowCut(i);
1736}
1737// Create C++ lines to get to current state
1738void
1739CbcTreeVariable::generateCpp( FILE * fp)
1740{
1741    CbcTreeVariable other;
1742    fprintf(fp, "0#include \"CbcTreeVariable.hpp\"\n");
1743    fprintf(fp, "5  CbcTreeVariable variableTree(cbcModel,NULL);\n");
1744    if (range_ != other.range_)
1745        fprintf(fp, "5  variableTree.setRange(%d);\n", range_);
1746    if (typeCuts_ != other.typeCuts_)
1747        fprintf(fp, "5  variableTree.setTypeCuts(%d);\n", typeCuts_);
1748    if (maxDiversification_ != other.maxDiversification_)
1749        fprintf(fp, "5  variableTree.setMaxDiversification(%d);\n", maxDiversification_);
1750    if (timeLimit_ != other.timeLimit_)
1751        fprintf(fp, "5  variableTree.setTimeLimit(%d);\n", timeLimit_);
1752    if (nodeLimit_ != other.nodeLimit_)
1753        fprintf(fp, "5  variableTree.setNodeLimit(%d);\n", nodeLimit_);
1754    if (refine_ != other.refine_)
1755        fprintf(fp, "5  variableTree.setRefine(%s);\n", refine_ ? "true" : "false");
1756    fprintf(fp, "5  cbcModel->passInTreeHandler(variableTree);\n");
1757}
1758
1759
1760
Note: See TracBrowser for help on using the repository browser.