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

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

Changed formatting using AStyle -A4 -p

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 64.6 KB
Line 
1/* $Id: CbcTreeLocal.cpp 1286 2009-11-09 23:33:07Z EdwinStraver $ */
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            if (!typeCuts_ && originalUpper_[i] - originalLower_[i] > 1.0)
750                continue; // skip as not 0-1
751            if (originalLower_[i] == originalUpper_[i])
752                continue;
753            double mu = 1.0 / (originalUpper_[i] - originalLower_[i]);
754            if (value == originalLower_[i]) {
755                rhs += mu * originalLower_[i];
756                cut.insert(iColumn, 1.0);
757                maxValue += originalUpper_[i];
758            } else if (value == originalUpper_[i]) {
759                rhs -= mu * originalUpper_[i];
760                cut.insert(iColumn, -1.0);
761                maxValue += originalLower_[i];
762            }
763        }
764        if (maxValue < rhs - primalTolerance) {
765            if (model_->messageHandler()->logLevel() > 0)
766                printf("slack cut\n");
767            goodSolution = 1;
768        }
769        rowCut.setRow(cut);
770        rowCut.setLb(-COIN_DBL_MAX);
771        rowCut.setUb(rhs);
772        rowCut.setGloballyValid();
773        if (model_->messageHandler()->logLevel() > 0)
774            printf("Cut size: %i Cut rhs: %g\n", cut.getNumElements(), rhs);
775#ifdef CBC_DEBUG
776        if (model_->messageHandler()->logLevel() > 0) {
777            int k;
778            for (k = 0; k < cut.getNumElements(); k++) {
779                printf("%i    %g ", cut.getIndices()[k], cut.getElements()[k]);
780                if ((k + 1) % 5 == 0)
781                    printf("\n");
782            }
783            if (k % 5 != 0)
784                printf("\n");
785        }
786#endif
787        return goodSolution;
788    } else {
789        if (model_->messageHandler()->logLevel() > 0)
790            printf("Not a good solution\n");
791        return -1;
792    }
793}
794// Other side of last cut branch
795void
796CbcTreeLocal::reverseCut(int state, double bias)
797{
798    // find global cut
799    OsiCuts * global = model_->globalCuts();
800    int n = global->sizeRowCuts();
801    int i;
802    OsiRowCut * rowCut = NULL;
803    for ( i = 0; i < n; i++) {
804        rowCut = global->rowCutPtr(i);
805        if (cut_ == *rowCut) {
806            break;
807        }
808    }
809    if (!rowCut) {
810        // must have got here in odd way e.g. strong branching
811        return;
812    }
813    if (rowCut->lb() > -1.0e10)
814        return;
815    // get smallest element
816    double smallest = COIN_DBL_MAX;
817    CoinPackedVector row = cut_.row();
818    for (int k = 0; k < row.getNumElements(); k++)
819        smallest = CoinMin(smallest, fabs(row.getElements()[k]));
820    if (!typeCuts_ && !refine_) {
821        // Reverse cut very very weakly
822        if (state > 2)
823            smallest = 0.0;
824    }
825    // replace by other way
826    if (model_->messageHandler()->logLevel() > 0)
827        printf("reverseCut - changing cut %d out of %d, old rhs %g %g ",
828               i, n, rowCut->lb(), rowCut->ub());
829    rowCut->setLb(rowCut->ub() + smallest - bias);
830    rowCut->setUb(COIN_DBL_MAX);
831    if (model_->messageHandler()->logLevel() > 0)
832        printf("new rhs %g %g, bias %g smallest %g ",
833               rowCut->lb(), rowCut->ub(), bias, smallest);
834    const OsiRowCutDebugger *debugger = model_->solver()->getRowCutDebuggerAlways() ;
835    if (debugger) {
836        if (debugger->invalidCut(*rowCut))
837            printf("ZZZZTree Global cut - cuts off optimal solution!\n");
838    }
839}
840// Delete last cut branch
841void
842CbcTreeLocal::deleteCut(OsiRowCut & cut)
843{
844    // find global cut
845    OsiCuts * global = model_->globalCuts();
846    int n = global->sizeRowCuts();
847    int i;
848    OsiRowCut * rowCut = NULL;
849    for ( i = 0; i < n; i++) {
850        rowCut = global->rowCutPtr(i);
851        if (cut == *rowCut) {
852            break;
853        }
854    }
855    assert (i < n);
856    // delete last cut
857    if (model_->messageHandler()->logLevel() > 0)
858        printf("deleteCut - deleting cut %d out of %d, rhs %g %g\n",
859               i, n, rowCut->lb(), rowCut->ub());
860    global->eraseRowCut(i);
861}
862// Create C++ lines to get to current state
863void
864CbcTreeLocal::generateCpp( FILE * fp)
865{
866    CbcTreeLocal other;
867    fprintf(fp, "0#include \"CbcTreeLocal.hpp\"\n");
868    fprintf(fp, "5  CbcTreeLocal localTree(cbcModel,NULL);\n");
869    if (range_ != other.range_)
870        fprintf(fp, "5  localTree.setRange(%d);\n", range_);
871    if (typeCuts_ != other.typeCuts_)
872        fprintf(fp, "5  localTree.setTypeCuts(%d);\n", typeCuts_);
873    if (maxDiversification_ != other.maxDiversification_)
874        fprintf(fp, "5  localTree.setMaxDiversification(%d);\n", maxDiversification_);
875    if (timeLimit_ != other.timeLimit_)
876        fprintf(fp, "5  localTree.setTimeLimit(%d);\n", timeLimit_);
877    if (nodeLimit_ != other.nodeLimit_)
878        fprintf(fp, "5  localTree.setNodeLimit(%d);\n", nodeLimit_);
879    if (refine_ != other.refine_)
880        fprintf(fp, "5  localTree.setRefine(%s);\n", refine_ ? "true" : "false");
881    fprintf(fp, "5  cbcModel->passInTreeHandler(localTree);\n");
882}
883
884
885CbcTreeVariable::CbcTreeVariable()
886        : localNode_(NULL),
887        bestSolution_(NULL),
888        savedSolution_(NULL),
889        saveNumberSolutions_(0),
890        model_(NULL),
891        originalLower_(NULL),
892        originalUpper_(NULL),
893        range_(0),
894        typeCuts_(-1),
895        maxDiversification_(0),
896        diversification_(0),
897        nextStrong_(false),
898        rhs_(0.0),
899        savedGap_(0.0),
900        bestCutoff_(0.0),
901        timeLimit_(0),
902        startTime_(0),
903        nodeLimit_(0),
904        startNode_(-1),
905        searchType_(-1),
906        refine_(false)
907{
908
909}
910/* Constructor with solution.
911   range is upper bound on difference from given solution.
912   maxDiversification is maximum number of diversifications to try
913   timeLimit is seconds in subTree
914   nodeLimit is nodes in subTree
915*/
916CbcTreeVariable::CbcTreeVariable(CbcModel * model, const double * solution ,
917                                 int range, int typeCuts, int maxDiversification,
918                                 int timeLimit, int nodeLimit, bool refine)
919        : localNode_(NULL),
920        bestSolution_(NULL),
921        savedSolution_(NULL),
922        saveNumberSolutions_(0),
923        model_(model),
924        originalLower_(NULL),
925        originalUpper_(NULL),
926        range_(range),
927        typeCuts_(typeCuts),
928        maxDiversification_(maxDiversification),
929        diversification_(0),
930        nextStrong_(false),
931        rhs_(0.0),
932        savedGap_(0.0),
933        bestCutoff_(0.0),
934        timeLimit_(timeLimit),
935        startTime_(0),
936        nodeLimit_(nodeLimit),
937        startNode_(-1),
938        searchType_(-1),
939        refine_(refine)
940{
941
942    OsiSolverInterface * solver = model_->solver();
943    const double * lower = solver->getColLower();
944    const double * upper = solver->getColUpper();
945    //const double * solution = solver->getColSolution();
946    //const double * objective = solver->getObjCoefficients();
947    double primalTolerance;
948    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
949
950    // Get increment
951    model_->analyzeObjective();
952
953    {
954        // needed to sync cutoffs
955        double value ;
956        solver->getDblParam(OsiDualObjectiveLimit, value) ;
957        model_->setCutoff(value * solver->getObjSense());
958    }
959    bestCutoff_ = model_->getCutoff();
960    // save current gap
961    savedGap_ = model_->getDblParam(CbcModel::CbcAllowableGap);
962
963    // make sure integers found
964    model_->findIntegers(false);
965    int numberIntegers = model_->numberIntegers();
966    const int * integerVariable = model_->integerVariable();
967    int i;
968    double direction = solver->getObjSense();
969    double newSolutionValue = 1.0e50;
970    if (solution) {
971        // copy solution
972        solver->setColSolution(solution);
973        newSolutionValue = direction * solver->getObjValue();
974    }
975    originalLower_ = new double [numberIntegers];
976    originalUpper_ = new double [numberIntegers];
977    bool all01 = true;
978    int number01 = 0;
979    for (i = 0; i < numberIntegers; i++) {
980        int iColumn = integerVariable[i];
981        originalLower_[i] = lower[iColumn];
982        originalUpper_[i] = upper[iColumn];
983        if (upper[iColumn] - lower[iColumn] > 1.5)
984            all01 = false;
985        else if (upper[iColumn] - lower[iColumn] == 1.0)
986            number01++;
987    }
988    if (all01 && !typeCuts_)
989        typeCuts_ = 1; // may as well so we don't have to deal with refine
990    if (!number01 && !typeCuts_) {
991        if (model_->messageHandler()->logLevel() > 0)
992            printf("** No 0-1 variables and local search only on 0-1 - switching off\n");
993        typeCuts_ = -1;
994    } else {
995        if (model_->messageHandler()->logLevel() > 0) {
996            std::string type;
997            if (all01) {
998                printf("%d 0-1 variables normal local  cuts\n",
999                       number01);
1000            } else if (typeCuts_) {
1001                printf("%d 0-1 variables, %d other - general integer local cuts\n",
1002                       number01, numberIntegers - number01);
1003            } else {
1004                printf("%d 0-1 variables, %d other - local cuts but just on 0-1 variables\n",
1005                       number01, numberIntegers - number01);
1006            }
1007            printf("maximum diversifications %d, initial cutspace %d, max time %d seconds, max nodes %d\n",
1008                   maxDiversification_, range_, timeLimit_, nodeLimit_);
1009        }
1010    }
1011    int numberColumns = model_->getNumCols();
1012    savedSolution_ = new double [numberColumns];
1013    memset(savedSolution_, 0, numberColumns*sizeof(double));
1014    if (solution) {
1015        rhs_ = range_;
1016        // Check feasible
1017        int goodSolution = createCut(solution, cut_);
1018        if (goodSolution >= 0) {
1019            for (i = 0; i < numberIntegers; i++) {
1020                int iColumn = integerVariable[i];
1021                double value = floor(solution[iColumn] + 0.5);
1022                // fix so setBestSolution will work
1023                solver->setColLower(iColumn, value);
1024                solver->setColUpper(iColumn, value);
1025            }
1026            model_->reserveCurrentSolution();
1027            // Create cut and get total gap
1028            if (newSolutionValue < bestCutoff_) {
1029                model_->setBestSolution(CBC_ROUNDING, newSolutionValue, solution);
1030                bestCutoff_ = model_->getCutoff();
1031                // save as best solution
1032                memcpy(savedSolution_, model_->bestSolution(), numberColumns*sizeof(double));
1033            }
1034            for (i = 0; i < numberIntegers; i++) {
1035                int iColumn = integerVariable[i];
1036                // restore bounds
1037                solver->setColLower(iColumn, originalLower_[i]);
1038                solver->setColUpper(iColumn, originalUpper_[i]);
1039            }
1040            // make sure can't stop on gap
1041            model_->setDblParam(CbcModel::CbcAllowableGap, -1.0e50);
1042        } else {
1043            model_ = NULL;
1044        }
1045    } else {
1046        // no solution
1047        rhs_ = 1.0e50;
1048        // make sure can't stop on gap
1049        model_->setDblParam(CbcModel::CbcAllowableGap, -1.0e50);
1050    }
1051}
1052CbcTreeVariable::~CbcTreeVariable()
1053{
1054    delete [] originalLower_;
1055    delete [] originalUpper_;
1056    delete [] bestSolution_;
1057    delete [] savedSolution_;
1058    delete localNode_;
1059}
1060// Copy constructor
1061CbcTreeVariable::CbcTreeVariable ( const CbcTreeVariable & rhs)
1062        : CbcTree(rhs),
1063        saveNumberSolutions_(rhs.saveNumberSolutions_),
1064        model_(rhs.model_),
1065        range_(rhs.range_),
1066        typeCuts_(rhs.typeCuts_),
1067        maxDiversification_(rhs.maxDiversification_),
1068        diversification_(rhs.diversification_),
1069        nextStrong_(rhs.nextStrong_),
1070        rhs_(rhs.rhs_),
1071        savedGap_(rhs.savedGap_),
1072        bestCutoff_(rhs.bestCutoff_),
1073        timeLimit_(rhs.timeLimit_),
1074        startTime_(rhs.startTime_),
1075        nodeLimit_(rhs.nodeLimit_),
1076        startNode_(rhs.startNode_),
1077        searchType_(rhs.searchType_),
1078        refine_(rhs.refine_)
1079{
1080    cut_ = rhs.cut_;
1081    fixedCut_ = rhs.fixedCut_;
1082    if (rhs.localNode_)
1083        localNode_ = new CbcNode(*rhs.localNode_);
1084    else
1085        localNode_ = NULL;
1086    if (rhs.originalLower_) {
1087        int numberIntegers = model_->numberIntegers();
1088        originalLower_ = new double [numberIntegers];
1089        memcpy(originalLower_, rhs.originalLower_, numberIntegers*sizeof(double));
1090        originalUpper_ = new double [numberIntegers];
1091        memcpy(originalUpper_, rhs.originalUpper_, numberIntegers*sizeof(double));
1092    } else {
1093        originalLower_ = NULL;
1094        originalUpper_ = NULL;
1095    }
1096    if (rhs.bestSolution_) {
1097        int numberColumns = model_->getNumCols();
1098        bestSolution_ = new double [numberColumns];
1099        memcpy(bestSolution_, rhs.bestSolution_, numberColumns*sizeof(double));
1100    } else {
1101        bestSolution_ = NULL;
1102    }
1103    if (rhs.savedSolution_) {
1104        int numberColumns = model_->getNumCols();
1105        savedSolution_ = new double [numberColumns];
1106        memcpy(savedSolution_, rhs.savedSolution_, numberColumns*sizeof(double));
1107    } else {
1108        savedSolution_ = NULL;
1109    }
1110}
1111//----------------------------------------------------------------
1112// Assignment operator
1113//-------------------------------------------------------------------
1114CbcTreeVariable &
1115CbcTreeVariable::operator=(const CbcTreeVariable & rhs)
1116{
1117    if (this != &rhs) {
1118        CbcTree::operator=(rhs);
1119        saveNumberSolutions_ = rhs.saveNumberSolutions_;
1120        cut_ = rhs.cut_;
1121        fixedCut_ = rhs.fixedCut_;
1122        delete localNode_;
1123        if (rhs.localNode_)
1124            localNode_ = new CbcNode(*rhs.localNode_);
1125        else
1126            localNode_ = NULL;
1127        model_ = rhs.model_;
1128        range_ = rhs.range_;
1129        typeCuts_ = rhs.typeCuts_;
1130        maxDiversification_ = rhs.maxDiversification_;
1131        diversification_ = rhs.diversification_;
1132        nextStrong_ = rhs.nextStrong_;
1133        rhs_ = rhs.rhs_;
1134        savedGap_ = rhs.savedGap_;
1135        bestCutoff_ = rhs.bestCutoff_;
1136        timeLimit_ = rhs.timeLimit_;
1137        startTime_ = rhs.startTime_;
1138        nodeLimit_ = rhs.nodeLimit_;
1139        startNode_ = rhs.startNode_;
1140        searchType_ = rhs.searchType_;
1141        refine_ = rhs.refine_;
1142        delete [] originalLower_;
1143        delete [] originalUpper_;
1144        if (rhs.originalLower_) {
1145            int numberIntegers = model_->numberIntegers();
1146            originalLower_ = new double [numberIntegers];
1147            memcpy(originalLower_, rhs.originalLower_, numberIntegers*sizeof(double));
1148            originalUpper_ = new double [numberIntegers];
1149            memcpy(originalUpper_, rhs.originalUpper_, numberIntegers*sizeof(double));
1150        } else {
1151            originalLower_ = NULL;
1152            originalUpper_ = NULL;
1153        }
1154        delete [] bestSolution_;
1155        if (rhs.bestSolution_) {
1156            int numberColumns = model_->getNumCols();
1157            bestSolution_ = new double [numberColumns];
1158            memcpy(bestSolution_, rhs.bestSolution_, numberColumns*sizeof(double));
1159        } else {
1160            bestSolution_ = NULL;
1161        }
1162        delete [] savedSolution_;
1163        if (rhs.savedSolution_) {
1164            int numberColumns = model_->getNumCols();
1165            savedSolution_ = new double [numberColumns];
1166            memcpy(savedSolution_, rhs.savedSolution_, numberColumns*sizeof(double));
1167        } else {
1168            savedSolution_ = NULL;
1169        }
1170    }
1171    return *this;
1172}
1173// Clone
1174CbcTree *
1175CbcTreeVariable::clone() const
1176{
1177    return new CbcTreeVariable(*this);
1178}
1179// Pass in solution (so can be used after heuristic)
1180void
1181CbcTreeVariable::passInSolution(const double * solution, double solutionValue)
1182{
1183    int numberColumns = model_->getNumCols();
1184    delete [] savedSolution_;
1185    savedSolution_ = new double [numberColumns];
1186    memcpy(savedSolution_, solution, numberColumns*sizeof(double));
1187    rhs_ = range_;
1188    // Check feasible
1189    int goodSolution = createCut(solution, cut_);
1190    if (goodSolution >= 0) {
1191        bestCutoff_ = CoinMin(solutionValue, model_->getCutoff());
1192    } else {
1193        model_ = NULL;
1194    }
1195}
1196// Return the top node of the heap
1197CbcNode *
1198CbcTreeVariable::top() const
1199{
1200#ifdef CBC_DEBUG
1201    int smallest = 9999999;
1202    int largest = -1;
1203    double smallestD = 1.0e30;
1204    double largestD = -1.0e30;
1205    int n = nodes_.size();
1206    for (int i = 0; i < n; i++) {
1207        int nn = nodes_[i]->nodeInfo()->nodeNumber();
1208        double dd = nodes_[i]->objectiveValue();
1209        largest = CoinMax(largest, nn);
1210        smallest = CoinMin(smallest, nn);
1211        largestD = CoinMax(largestD, dd);
1212        smallestD = CoinMin(smallestD, dd);
1213    }
1214    if (model_->messageHandler()->logLevel() > 0) {
1215        printf("smallest %d, largest %d, top %d\n", smallest, largest,
1216               nodes_.front()->nodeInfo()->nodeNumber());
1217        printf("smallestD %g, largestD %g, top %g\n", smallestD, largestD, nodes_.front()->objectiveValue());
1218    }
1219#endif
1220    return nodes_.front();
1221}
1222
1223// Add a node to the heap
1224void
1225CbcTreeVariable::push(CbcNode * x)
1226{
1227    if (typeCuts_ >= 0 && !nodes_.size() && searchType_ < 0) {
1228        startNode_ = model_->getNodeCount();
1229        // save copy of node
1230        localNode_ = new CbcNode(*x);
1231
1232        if (cut_.row().getNumElements()) {
1233            // Add to global cuts
1234            // we came in with solution
1235            model_->globalCuts()->insert(cut_);
1236            if (model_->messageHandler()->logLevel() > 0)
1237                printf("initial cut - rhs %g %g\n",
1238                       cut_.lb(), cut_.ub());
1239            searchType_ = 1;
1240        } else {
1241            // stop on first solution
1242            searchType_ = 0;
1243        }
1244        startTime_ = static_cast<int> (CoinCpuTime());
1245        saveNumberSolutions_ = model_->getSolutionCount();
1246    }
1247    nodes_.push_back(x);
1248#ifdef CBC_DEBUG
1249    if (model_->messageHandler()->logLevel() > 0)
1250        printf("pushing node onto heap %d %x %x\n",
1251               x->nodeInfo()->nodeNumber(), x, x->nodeInfo());
1252#endif
1253    std::push_heap(nodes_.begin(), nodes_.end(), comparison_);
1254}
1255
1256// Remove the top node from the heap
1257void
1258CbcTreeVariable::pop()
1259{
1260    std::pop_heap(nodes_.begin(), nodes_.end(), comparison_);
1261    nodes_.pop_back();
1262}
1263// Test if empty - does work if so
1264bool
1265CbcTreeVariable::empty()
1266{
1267    if (typeCuts_ < 0)
1268        return !nodes_.size();
1269    /* state -
1270       0 iterating
1271       1 subtree finished optimal solution for subtree found
1272       2 subtree finished and no solution found
1273       3 subtree exiting and solution found
1274       4 subtree exiting and no solution found
1275    */
1276    int state = 0;
1277    assert (searchType_ != 2);
1278    if (searchType_) {
1279        if (CoinCpuTime() - startTime_ > timeLimit_ || model_->getNodeCount() - startNode_ >= nodeLimit_) {
1280            state = 4;
1281        }
1282    } else {
1283        if (model_->getSolutionCount() > saveNumberSolutions_) {
1284            state = 4;
1285        }
1286    }
1287    if (!nodes_.size())
1288        state = 2;
1289    if (!state) {
1290        return false;
1291    }
1292    // Finished this phase
1293    int numberColumns = model_->getNumCols();
1294    if (model_->getSolutionCount() > saveNumberSolutions_) {
1295        if (model_->getCutoff() < bestCutoff_) {
1296            // Save solution
1297            if (!bestSolution_)
1298                bestSolution_ = new double [numberColumns];
1299            memcpy(bestSolution_, model_->bestSolution(), numberColumns*sizeof(double));
1300            bestCutoff_ = model_->getCutoff();
1301        }
1302        state--;
1303    }
1304    // get rid of all nodes (safe even if already done)
1305    double bestPossibleObjective;
1306    cleanTree(model_, -COIN_DBL_MAX, bestPossibleObjective);
1307
1308    double increment = model_->getDblParam(CbcModel::CbcCutoffIncrement) ;
1309    if (model_->messageHandler()->logLevel() > 0)
1310        printf("local state %d after %d nodes and %d seconds, new solution %g, best solution %g, k was %g\n",
1311               state,
1312               model_->getNodeCount() - startNode_,
1313               static_cast<int> (CoinCpuTime()) - startTime_,
1314               model_->getCutoff() + increment, bestCutoff_ + increment, rhs_);
1315    saveNumberSolutions_ = model_->getSolutionCount();
1316    bool finished = false;
1317    bool lastTry = false;
1318    switch (state) {
1319    case 1:
1320        // solution found and subtree exhausted
1321        if (rhs_ > 1.0e30) {
1322            finished = true;
1323        } else {
1324            // find global cut and reverse
1325            reverseCut(1);
1326            searchType_ = 1; // first false
1327            rhs_ = range_; // reset range
1328            nextStrong_ = false;
1329
1330            // save best solution in this subtree
1331            memcpy(savedSolution_, model_->bestSolution(), numberColumns*sizeof(double));
1332        }
1333        break;
1334    case 2:
1335        // solution not found and subtree exhausted
1336        if (rhs_ > 1.0e30) {
1337            finished = true;
1338        } else {
1339            // find global cut and reverse
1340            reverseCut(2);
1341            searchType_ = 1; // first false
1342            if (diversification_ < maxDiversification_) {
1343                if (nextStrong_) {
1344                    diversification_++;
1345                    // cut is valid so don't model_->setCutoff(1.0e50);
1346                    searchType_ = 0;
1347                }
1348                nextStrong_ = true;
1349                rhs_ += range_ / 2;
1350            } else {
1351                // This will be last try (may hit max time)
1352                lastTry = true;
1353                if (!maxDiversification_)
1354                    typeCuts_ = -1; // make sure can't start again
1355                model_->setCutoff(bestCutoff_);
1356                if (model_->messageHandler()->logLevel() > 0)
1357                    printf("Exiting local search with current set of cuts\n");
1358                rhs_ = 1.0e100;
1359                // Can now stop on gap
1360                model_->setDblParam(CbcModel::CbcAllowableGap, savedGap_);
1361            }
1362        }
1363        break;
1364    case 3:
1365        // solution found and subtree not exhausted
1366        if (rhs_ < 1.0e30) {
1367            if (searchType_) {
1368                if (!typeCuts_ && refine_ && searchType_ == 1) {
1369                    // We need to check we have best solution given these 0-1 values
1370                    OsiSolverInterface * subSolver = model_->continuousSolver()->clone();
1371                    CbcModel * subModel = model_->subTreeModel(subSolver);
1372                    CbcTree normalTree;
1373                    subModel->passInTreeHandler(normalTree);
1374                    int numberIntegers = model_->numberIntegers();
1375                    const int * integerVariable = model_->integerVariable();
1376                    const double * solution = model_->bestSolution();
1377                    int i;
1378                    int numberColumns = model_->getNumCols();
1379                    for (i = 0; i < numberIntegers; i++) {
1380                        int iColumn = integerVariable[i];
1381                        double value = floor(solution[iColumn] + 0.5);
1382                        if (!typeCuts_ && originalUpper_[i] - originalLower_[i] > 1.0)
1383                            continue; // skip as not 0-1
1384                        if (originalLower_[i] == originalUpper_[i])
1385                            continue;
1386                        subSolver->setColLower(iColumn, value);
1387                        subSolver->setColUpper(iColumn, value);
1388                    }
1389                    subSolver->initialSolve();
1390                    // We can copy cutoff
1391                    // But adjust
1392                    subModel->setCutoff(model_->getCutoff() + model_->getDblParam(CbcModel::CbcCutoffIncrement) + 1.0e-6);
1393                    subModel->setSolutionCount(0);
1394                    assert (subModel->isProvenOptimal());
1395                    if (!subModel->typePresolve()) {
1396                        subModel->branchAndBound();
1397                        if (subModel->status()) {
1398                            model_->incrementSubTreeStopped();
1399                        }
1400                        //printf("%g %g %g %g\n",subModel->getCutoff(),model_->getCutoff(),
1401                        //   subModel->getMinimizationObjValue(),model_->getMinimizationObjValue());
1402                        double newCutoff = subModel->getMinimizationObjValue() -
1403                                           subModel->getDblParam(CbcModel::CbcCutoffIncrement) ;
1404                        if (subModel->getSolutionCount()) {
1405                            if (!subModel->status())
1406                                assert (subModel->isProvenOptimal());
1407                            memcpy(model_->bestSolution(), subModel->bestSolution(),
1408                                   numberColumns*sizeof(double));
1409                            model_->setCutoff(newCutoff);
1410                        }
1411                    } else if (subModel->typePresolve() == 1) {
1412                        CbcModel * model2 = subModel->integerPresolve(true);
1413                        if (model2) {
1414                            // Do complete search
1415                            model2->branchAndBound();
1416                            // get back solution
1417                            subModel->originalModel(model2, false);
1418                            if (model2->status()) {
1419                                model_->incrementSubTreeStopped();
1420                            }
1421                            double newCutoff = model2->getMinimizationObjValue() -
1422                                               model2->getDblParam(CbcModel::CbcCutoffIncrement) ;
1423                            if (model2->getSolutionCount()) {
1424                                if (!model2->status())
1425                                    assert (model2->isProvenOptimal());
1426                                memcpy(model_->bestSolution(), subModel->bestSolution(),
1427                                       numberColumns*sizeof(double));
1428                                model_->setCutoff(newCutoff);
1429                            }
1430                            delete model2;
1431                        } else {
1432                            // infeasible - could just be - due to cutoff
1433                        }
1434                    } else {
1435                        // too dangerous at present
1436                        assert (subModel->typePresolve() != 2);
1437                    }
1438                    if (model_->getCutoff() < bestCutoff_) {
1439                        // Save solution
1440                        if (!bestSolution_)
1441                            bestSolution_ = new double [numberColumns];
1442                        memcpy(bestSolution_, model_->bestSolution(), numberColumns*sizeof(double));
1443                        bestCutoff_ = model_->getCutoff();
1444                    }
1445                    delete subModel;
1446                }
1447                // we have done search to make sure best general solution
1448                searchType_ = 1;
1449                // Reverse cut weakly
1450                reverseCut(3, rhs_);
1451            } else {
1452                searchType_ = 1;
1453                // delete last cut
1454                deleteCut(cut_);
1455            }
1456        } else {
1457            searchType_ = 1;
1458        }
1459        // save best solution in this subtree
1460        memcpy(savedSolution_, model_->bestSolution(), numberColumns*sizeof(double));
1461        nextStrong_ = false;
1462        rhs_ = range_;
1463        break;
1464    case 4:
1465        // solution not found and subtree not exhausted
1466        if (maxDiversification_) {
1467            if (nextStrong_) {
1468                // Reverse cut weakly
1469                reverseCut(4, rhs_);
1470                model_->setCutoff(1.0e50);
1471                diversification_++;
1472                searchType_ = 0;
1473            } else {
1474                // delete last cut
1475                deleteCut(cut_);
1476                searchType_ = 1;
1477            }
1478            nextStrong_ = true;
1479            rhs_ += range_ / 2;
1480        } else {
1481            // special case when using as heuristic
1482            // Reverse cut weakly if lb -infinity
1483            reverseCut(4, rhs_);
1484            // This will be last try (may hit max time0
1485            lastTry = true;
1486            model_->setCutoff(bestCutoff_);
1487            if (model_->messageHandler()->logLevel() > 0)
1488                printf("Exiting local search with current set of cuts\n");
1489            rhs_ = 1.0e100;
1490            // Can now stop on gap
1491            model_->setDblParam(CbcModel::CbcAllowableGap, savedGap_);
1492            typeCuts_ = -1;
1493        }
1494        break;
1495    }
1496    if (rhs_ < 1.0e30 || lastTry) {
1497        int goodSolution = createCut(savedSolution_, cut_);
1498        if (goodSolution >= 0) {
1499            // Add to global cuts
1500            model_->globalCuts()->insert(cut_);
1501            OsiCuts * global = model_->globalCuts();
1502            int n = global->sizeRowCuts();
1503            OsiRowCut * rowCut = global->rowCutPtr(n - 1);
1504            if (model_->messageHandler()->logLevel() > 0)
1505                printf("inserting cut - now %d cuts, rhs %g %g, cutspace %g, diversification %d\n",
1506                       n, rowCut->lb(), rowCut->ub(), rhs_, diversification_);
1507            const OsiRowCutDebugger *debugger = model_->solver()->getRowCutDebuggerAlways() ;
1508            if (debugger) {
1509                if (debugger->invalidCut(*rowCut))
1510                    printf("ZZZZTree Global cut - cuts off optimal solution!\n");
1511            }
1512            for (int i = 0; i < n; i++) {
1513                rowCut = global->rowCutPtr(i);
1514                if (model_->messageHandler()->logLevel() > 0)
1515                    printf("%d - rhs %g %g\n",
1516                           i, rowCut->lb(), rowCut->ub());
1517            }
1518        }
1519        // put back node
1520        startTime_ = static_cast<int> (CoinCpuTime());
1521        startNode_ = model_->getNodeCount();
1522        if (localNode_) {
1523            // save copy of node
1524            CbcNode * localNode2 = new CbcNode(*localNode_);
1525            // But localNode2 now owns cuts so swap
1526            //printf("pushing local node2 onto heap %d %x %x\n",localNode_->nodeNumber(),
1527            //   localNode_,localNode_->nodeInfo());
1528            nodes_.push_back(localNode_);
1529            localNode_ = localNode2;
1530            std::make_heap(nodes_.begin(), nodes_.end(), comparison_);
1531        }
1532    }
1533    return finished;
1534}
1535// We may have got an intelligent tree so give it one more chance
1536void
1537CbcTreeVariable::endSearch()
1538{
1539    if (typeCuts_ >= 0) {
1540        // copy best solution to model
1541        int numberColumns = model_->getNumCols();
1542        if (bestSolution_ && bestCutoff_ < model_->getCutoff()) {
1543            memcpy(model_->bestSolution(), bestSolution_, numberColumns*sizeof(double));
1544            model_->setCutoff(bestCutoff_);
1545            // recompute objective value
1546            const double * objCoef = model_->getObjCoefficients();
1547            double objOffset = 0.0;
1548            model_->continuousSolver()->getDblParam(OsiObjOffset, objOffset);
1549
1550            // Compute dot product of objCoef and colSol and then adjust by offset
1551            double objValue = -objOffset;
1552            for ( int i = 0 ; i < numberColumns ; i++ )
1553                objValue += objCoef[i] * bestSolution_[i];
1554            model_->setMinimizationObjValue(objValue);
1555        }
1556        // Can now stop on gap
1557        model_->setDblParam(CbcModel::CbcAllowableGap, savedGap_);
1558    }
1559}
1560// Create cut
1561int
1562CbcTreeVariable::createCut(const double * solution, OsiRowCut & rowCut)
1563{
1564    if (rhs_ > 1.0e20)
1565        return -1;
1566    OsiSolverInterface * solver = model_->solver();
1567    const double * rowLower = solver->getRowLower();
1568    const double * rowUpper = solver->getRowUpper();
1569    //const double * solution = solver->getColSolution();
1570    //const double * objective = solver->getObjCoefficients();
1571    double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
1572    double primalTolerance;
1573    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
1574    // relax
1575    primalTolerance *= 1000.0;
1576
1577    int numberRows = model_->getNumRows();
1578
1579    int numberIntegers = model_->numberIntegers();
1580    const int * integerVariable = model_->integerVariable();
1581    int i;
1582
1583    // Check feasible
1584
1585    double * rowActivity = new double[numberRows];
1586    memset(rowActivity, 0, numberRows*sizeof(double));
1587    solver->getMatrixByCol()->times(solution, rowActivity) ;
1588    int goodSolution = 0;
1589    // check was feasible
1590    for (i = 0; i < numberRows; i++) {
1591        if (rowActivity[i] < rowLower[i] - primalTolerance) {
1592            goodSolution = -1;
1593        } else if (rowActivity[i] > rowUpper[i] + primalTolerance) {
1594            goodSolution = -1;
1595        }
1596    }
1597    delete [] rowActivity;
1598    for (i = 0; i < numberIntegers; i++) {
1599        int iColumn = integerVariable[i];
1600        double value = solution[iColumn];
1601        if (fabs(floor(value + 0.5) - value) > integerTolerance) {
1602            goodSolution = -1;
1603        }
1604    }
1605    // zap cut
1606    if (goodSolution == 0) {
1607        // Create cut and get total gap
1608        CoinPackedVector cut;
1609        double rhs = rhs_;
1610        double maxValue = 0.0;
1611        for (i = 0; i < numberIntegers; i++) {
1612            int iColumn = integerVariable[i];
1613            double value = floor(solution[iColumn] + 0.5);
1614            if (!typeCuts_ && originalUpper_[i] - originalLower_[i] > 1.0)
1615                continue; // skip as not 0-1
1616            if (originalLower_[i] == originalUpper_[i])
1617                continue;
1618            double mu = 1.0 / (originalUpper_[i] - originalLower_[i]);
1619            if (value == originalLower_[i]) {
1620                rhs += mu * originalLower_[i];
1621                cut.insert(iColumn, 1.0);
1622                maxValue += originalUpper_[i];
1623            } else if (value == originalUpper_[i]) {
1624                rhs -= mu * originalUpper_[i];
1625                cut.insert(iColumn, -1.0);
1626                maxValue += originalLower_[i];
1627            }
1628        }
1629        if (maxValue < rhs - primalTolerance) {
1630            if (model_->messageHandler()->logLevel() > 0)
1631                printf("slack cut\n");
1632            goodSolution = 1;
1633        }
1634        rowCut.setRow(cut);
1635        rowCut.setLb(-COIN_DBL_MAX);
1636        rowCut.setUb(rhs);
1637        rowCut.setGloballyValid();
1638        if (model_->messageHandler()->logLevel() > 0)
1639            printf("Cut size: %i Cut rhs: %g\n", cut.getNumElements(), rhs);
1640#ifdef CBC_DEBUG
1641        if (model_->messageHandler()->logLevel() > 0) {
1642            int k;
1643            for (k = 0; k < cut.getNumElements(); k++) {
1644                printf("%i    %g ", cut.getIndices()[k], cut.getElements()[k]);
1645                if ((k + 1) % 5 == 0)
1646                    printf("\n");
1647            }
1648            if (k % 5 != 0)
1649                printf("\n");
1650        }
1651#endif
1652        return goodSolution;
1653    } else {
1654        if (model_->messageHandler()->logLevel() > 0)
1655            printf("Not a good solution\n");
1656        return -1;
1657    }
1658}
1659// Other side of last cut branch
1660void
1661CbcTreeVariable::reverseCut(int state, double bias)
1662{
1663    // find global cut
1664    OsiCuts * global = model_->globalCuts();
1665    int n = global->sizeRowCuts();
1666    int i;
1667    OsiRowCut * rowCut = NULL;
1668    for ( i = 0; i < n; i++) {
1669        rowCut = global->rowCutPtr(i);
1670        if (cut_ == *rowCut) {
1671            break;
1672        }
1673    }
1674    if (!rowCut) {
1675        // must have got here in odd way e.g. strong branching
1676        return;
1677    }
1678    if (rowCut->lb() > -1.0e10)
1679        return;
1680    // get smallest element
1681    double smallest = COIN_DBL_MAX;
1682    CoinPackedVector row = cut_.row();
1683    for (int k = 0; k < row.getNumElements(); k++)
1684        smallest = CoinMin(smallest, fabs(row.getElements()[k]));
1685    if (!typeCuts_ && !refine_) {
1686        // Reverse cut very very weakly
1687        if (state > 2)
1688            smallest = 0.0;
1689    }
1690    // replace by other way
1691    if (model_->messageHandler()->logLevel() > 0)
1692        printf("reverseCut - changing cut %d out of %d, old rhs %g %g ",
1693               i, n, rowCut->lb(), rowCut->ub());
1694    rowCut->setLb(rowCut->ub() + smallest - bias);
1695    rowCut->setUb(COIN_DBL_MAX);
1696    if (model_->messageHandler()->logLevel() > 0)
1697        printf("new rhs %g %g, bias %g smallest %g ",
1698               rowCut->lb(), rowCut->ub(), bias, smallest);
1699    const OsiRowCutDebugger *debugger = model_->solver()->getRowCutDebuggerAlways() ;
1700    if (debugger) {
1701        if (debugger->invalidCut(*rowCut))
1702            printf("ZZZZTree Global cut - cuts off optimal solution!\n");
1703    }
1704}
1705// Delete last cut branch
1706void
1707CbcTreeVariable::deleteCut(OsiRowCut & cut)
1708{
1709    // find global cut
1710    OsiCuts * global = model_->globalCuts();
1711    int n = global->sizeRowCuts();
1712    int i;
1713    OsiRowCut * rowCut = NULL;
1714    for ( i = 0; i < n; i++) {
1715        rowCut = global->rowCutPtr(i);
1716        if (cut == *rowCut) {
1717            break;
1718        }
1719    }
1720    assert (i < n);
1721    // delete last cut
1722    if (model_->messageHandler()->logLevel() > 0)
1723        printf("deleteCut - deleting cut %d out of %d, rhs %g %g\n",
1724               i, n, rowCut->lb(), rowCut->ub());
1725    global->eraseRowCut(i);
1726}
1727// Create C++ lines to get to current state
1728void
1729CbcTreeVariable::generateCpp( FILE * fp)
1730{
1731    CbcTreeVariable other;
1732    fprintf(fp, "0#include \"CbcTreeVariable.hpp\"\n");
1733    fprintf(fp, "5  CbcTreeVariable variableTree(cbcModel,NULL);\n");
1734    if (range_ != other.range_)
1735        fprintf(fp, "5  variableTree.setRange(%d);\n", range_);
1736    if (typeCuts_ != other.typeCuts_)
1737        fprintf(fp, "5  variableTree.setTypeCuts(%d);\n", typeCuts_);
1738    if (maxDiversification_ != other.maxDiversification_)
1739        fprintf(fp, "5  variableTree.setMaxDiversification(%d);\n", maxDiversification_);
1740    if (timeLimit_ != other.timeLimit_)
1741        fprintf(fp, "5  variableTree.setTimeLimit(%d);\n", timeLimit_);
1742    if (nodeLimit_ != other.nodeLimit_)
1743        fprintf(fp, "5  variableTree.setNodeLimit(%d);\n", nodeLimit_);
1744    if (refine_ != other.refine_)
1745        fprintf(fp, "5  variableTree.setRefine(%s);\n", refine_ ? "true" : "false");
1746    fprintf(fp, "5  cbcModel->passInTreeHandler(variableTree);\n");
1747}
1748
1749
1750
Note: See TracBrowser for help on using the repository browser.