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

Last change on this file since 1389 was 1361, checked in by bjarni, 10 years ago

Added Lou's annotations to Cbc_ampl.cpp, CbcSolver?.cpp, CbcSolver?.hpp, CbcStrategy?.cpp, CbcTreeLocal?.cpp, ClpAmplStuff?.cpp, CoinSolve?.cpp, and unitTestClp.cpp

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