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

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

Changed formatting using AStyle -A4 -p

File size: 221.7 KB
Line 
1/* $Id: CbcNode.cpp 1271 2009-11-05 15:57:25Z forrest $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4#if defined(_MSC_VER)
5// Turn off compiler warning about long names
6#  pragma warning(disable:4786)
7#endif
8
9#include "CbcConfig.h"
10
11#include <string>
12//#define CBC_DEBUG 1
13//#define CHECK_CUT_COUNTS
14//#define CHECK_NODE
15//#define CBC_CHECK_BASIS
16#include <cassert>
17#include <cfloat>
18#define CUTS
19#include "OsiSolverInterface.hpp"
20#include "OsiChooseVariable.hpp"
21#include "OsiAuxInfo.hpp"
22#include "OsiSolverBranch.hpp"
23#include "CoinWarmStartBasis.hpp"
24#include "CoinTime.hpp"
25#include "CbcModel.hpp"
26#include "CbcNode.hpp"
27#include "CbcStatistics.hpp"
28#include "CbcStrategy.hpp"
29#include "CbcBranchActual.hpp"
30#include "CbcBranchDynamic.hpp"
31#include "OsiRowCut.hpp"
32#include "OsiRowCutDebugger.hpp"
33#include "OsiCuts.hpp"
34#include "CbcCountRowCut.hpp"
35#include "CbcFeasibilityBase.hpp"
36#include "CbcMessage.hpp"
37#ifdef COIN_HAS_CLP
38#include "OsiClpSolverInterface.hpp"
39#include "ClpSimplexOther.hpp"
40#endif
41using namespace std;
42#include "CglCutGenerator.hpp"
43// Default Constructor
44CbcNodeInfo::CbcNodeInfo ()
45        :
46        numberPointingToThis_(0),
47        parent_(NULL),
48        parentBranch_(NULL),
49        owner_(NULL),
50        numberCuts_(0),
51        nodeNumber_(0),
52        cuts_(NULL),
53        numberRows_(0),
54        numberBranchesLeft_(0),
55        active_(7)
56{
57#ifdef CHECK_NODE
58    printf("CbcNodeInfo %x Constructor\n", this);
59#endif
60}
61
62void
63CbcNodeInfo::setParentBasedData()
64{
65    if (parent_) {
66        numberRows_ = parent_->numberRows_ + parent_->numberCuts_;
67        //parent_->increment();
68        if (parent_->owner()) {
69            const OsiBranchingObject* br = parent_->owner()->branchingObject();
70            assert(br);
71            parentBranch_ = br->clone();
72        }
73    }
74}
75
76void
77CbcNodeInfo::unsetParentBasedData()
78{
79    if (parent_) {
80        numberRows_ = 0;
81        if (parent_->owner()) {
82            delete parentBranch_;
83            parentBranch_ = NULL;
84        }
85    }
86}
87
88#if 0
89// Constructor given parent
90CbcNodeInfo::CbcNodeInfo (CbcNodeInfo * parent)
91        :
92        numberPointingToThis_(2),
93        parent_(parent),
94        parentBranch_(NULL),
95        owner_(NULL),
96        numberCuts_(0),
97        nodeNumber_(0),
98        cuts_(NULL),
99        numberRows_(0),
100        numberBranchesLeft_(2),
101        active_(7)
102{
103#ifdef CHECK_NODE
104    printf("CbcNodeInfo %x Constructor from parent %x\n", this, parent_);
105#endif
106    //setParentBasedData();
107}
108#endif
109
110// Copy Constructor
111CbcNodeInfo::CbcNodeInfo (const CbcNodeInfo & rhs)
112        :
113        numberPointingToThis_(rhs.numberPointingToThis_),
114        parent_(rhs.parent_),
115        parentBranch_(NULL),
116        owner_(rhs.owner_),
117        numberCuts_(rhs.numberCuts_),
118        nodeNumber_(rhs.nodeNumber_),
119        cuts_(NULL),
120        numberRows_(rhs.numberRows_),
121        numberBranchesLeft_(rhs.numberBranchesLeft_),
122        active_(rhs.active_)
123{
124#ifdef CHECK_NODE
125    printf("CbcNodeInfo %x Copy constructor\n", this);
126#endif
127    if (numberCuts_) {
128        cuts_ = new CbcCountRowCut * [numberCuts_];
129        int n = 0;
130        for (int i = 0; i < numberCuts_; i++) {
131            CbcCountRowCut * thisCut = rhs.cuts_[i];
132            if (thisCut) {
133                // I think this is correct - new one should take priority
134                thisCut->setInfo(this, n);
135                thisCut->increment(numberBranchesLeft_);
136                cuts_[n++] = thisCut;
137            }
138        }
139        numberCuts_ = n;
140    }
141    if (rhs.parentBranch_) {
142        parentBranch_ = rhs.parentBranch_->clone();
143    }
144}
145// Constructor given parent and owner
146CbcNodeInfo::CbcNodeInfo (CbcNodeInfo * parent, CbcNode * owner)
147        :
148        numberPointingToThis_(2),
149        parent_(parent),
150        parentBranch_(NULL),
151        owner_(owner),
152        numberCuts_(0),
153        nodeNumber_(0),
154        cuts_(NULL),
155        numberRows_(0),
156        numberBranchesLeft_(2),
157        active_(7)
158{
159#ifdef CHECK_NODE
160    printf("CbcNodeInfo %x Constructor from parent %x\n", this, parent_);
161#endif
162    //setParentBasedData();
163}
164
165/**
166   Take care to detach from the owning CbcNode and decrement the reference
167   count in the parent.  If this is the last nodeInfo object pointing to the
168   parent, make a recursive call to delete the parent.
169*/
170CbcNodeInfo::~CbcNodeInfo()
171{
172#ifdef CHECK_NODE
173    printf("CbcNodeInfo %x Destructor parent %x\n", this, parent_);
174#endif
175
176    assert(!numberPointingToThis_);
177    // But there may be some left (max nodes?)
178    for (int i = 0; i < numberCuts_; i++) {
179        if (cuts_[i]) {
180#ifndef GLOBAL_CUTS_JUST_POINTERS
181            delete cuts_[i];
182#else
183            if (cuts_[i]->globallyValidAsInteger() != 2)
184                delete cuts_[i];
185#endif
186        }
187    }
188    delete [] cuts_;
189    if (owner_)
190        owner_->nullNodeInfo();
191    if (parent_) {
192        int numberLinks = parent_->decrement();
193        if (!numberLinks) delete parent_;
194    }
195    delete parentBranch_;
196}
197
198
199//#define ALLCUTS
200void
201CbcNodeInfo::decrementCuts(int change)
202{
203    int i;
204    // get rid of all remaining if negative
205    int changeThis;
206    if (change < 0)
207        changeThis = numberBranchesLeft_;
208    else
209        changeThis = change;
210    // decrement cut counts
211    for (i = 0; i < numberCuts_; i++) {
212        if (cuts_[i]) {
213            int number = cuts_[i]->decrement(changeThis);
214            if (!number) {
215                //printf("info %x del cut %d %x\n",this,i,cuts_[i]);
216#ifndef GLOBAL_CUTS_JUST_POINTERS
217                delete cuts_[i];
218#else
219                if (cuts_[i]->globallyValidAsInteger() != 2)
220                    delete cuts_[i];
221#endif
222                cuts_[i] = NULL;
223            }
224        }
225    }
226}
227void
228CbcNodeInfo::incrementCuts(int change)
229{
230    int i;
231    assert (change > 0);
232    // increment cut counts
233    for (i = 0; i < numberCuts_; i++) {
234        if (cuts_[i])
235            cuts_[i]->increment(change);
236    }
237}
238void
239CbcNodeInfo::decrementParentCuts(CbcModel * model, int change)
240{
241    if (parent_) {
242        // get rid of all remaining if negative
243        int changeThis;
244        if (change < 0)
245            changeThis = numberBranchesLeft_;
246        else
247            changeThis = change;
248        int i;
249        // Get over-estimate of space needed for basis
250        CoinWarmStartBasis & dummy = model->workingBasis();
251        dummy.setSize(0, numberRows_ + numberCuts_);
252        buildRowBasis(dummy);
253        /* everything is zero (i.e. free) so we can use to see
254           if latest basis */
255        CbcNodeInfo * thisInfo = parent_;
256        while (thisInfo)
257            thisInfo = thisInfo->buildRowBasis(dummy);
258        // decrement cut counts
259        thisInfo = parent_;
260        int numberRows = numberRows_;
261        while (thisInfo) {
262            for (i = thisInfo->numberCuts_ - 1; i >= 0; i--) {
263                CoinWarmStartBasis::Status status = dummy.getArtifStatus(--numberRows);
264#ifdef ALLCUTS
265                status = CoinWarmStartBasis::isFree;
266#endif
267                if (thisInfo->cuts_[i]) {
268                    int number = 1;
269                    if (status != CoinWarmStartBasis::basic) {
270                        // tight - drop 1 or 2
271                        if (change < 0)
272                            number = thisInfo->cuts_[i]->decrement(changeThis);
273                        else
274                            number = thisInfo->cuts_[i]->decrement(change);
275                    }
276                    if (!number) {
277#ifndef GLOBAL_CUTS_JUST_POINTERS
278                        delete thisInfo->cuts_[i];
279#else
280                        if (thisInfo->cuts_[i]->globallyValidAsInteger() != 2)
281                            delete thisInfo->cuts_[i];
282#endif
283                        thisInfo->cuts_[i] = NULL;
284                    }
285                }
286            }
287            thisInfo = thisInfo->parent_;
288        }
289    }
290}
291#if 0
292void
293CbcNodeInfo::incrementParentCuts(CbcModel * model, int change)
294{
295    if (parent_) {
296        int i;
297        // Get over-estimate of space needed for basis
298        CoinWarmStartBasis & dummy = model->workingBasis();
299        dummy.setSize(0, numberRows_ + numberCuts_);
300        /* everything is zero (i.e. free) so we can use to see
301           if latest basis */
302        buildRowBasis(dummy);
303        CbcNodeInfo * thisInfo = parent_;
304        while (thisInfo)
305            thisInfo = thisInfo->buildRowBasis(dummy);
306        // increment cut counts
307        thisInfo = parent_;
308        int numberRows = numberRows_;
309        while (thisInfo) {
310            for (i = thisInfo->numberCuts_ - 1; i >= 0; i--) {
311                CoinWarmStartBasis::Status status = dummy.getArtifStatus(--numberRows);
312#ifdef ALLCUTS
313                status = CoinWarmStartBasis::isFree;
314#endif
315                if (thisInfo->cuts_[i] && status != CoinWarmStartBasis::basic) {
316                    thisInfo->cuts_[i]->increment(change);
317                }
318            }
319            thisInfo = thisInfo->parent_;
320        }
321    }
322}
323#endif
324/*
325  Append cuts to the cuts_ array in a nodeInfo. The initial reference count
326  is set to numberToBranchOn, which will normally be the number of arms
327  defined for the CbcBranchingObject attached to the CbcNode that owns this
328  CbcNodeInfo.
329*/
330void
331CbcNodeInfo::addCuts (OsiCuts & cuts, int numberToBranchOn,
332                      /*int * whichGenerator,*/int numberPointingToThis)
333{
334    int numberCuts = cuts.sizeRowCuts();
335    if (numberCuts) {
336        int i;
337        if (!numberCuts_) {
338            cuts_ = new CbcCountRowCut * [numberCuts];
339        } else {
340            CbcCountRowCut ** temp = new CbcCountRowCut * [numberCuts+numberCuts_];
341            memcpy(temp, cuts_, numberCuts_*sizeof(CbcCountRowCut *));
342            delete [] cuts_;
343            cuts_ = temp;
344        }
345        for (i = 0; i < numberCuts; i++) {
346            CbcCountRowCut * thisCut = new CbcCountRowCut(*cuts.rowCutPtr(i),
347                    this, numberCuts_,
348                    -1, numberPointingToThis);
349            thisCut->increment(numberToBranchOn);
350            cuts_[numberCuts_++] = thisCut;
351#ifdef CBC_DEBUG
352#if CBC_DEBUG>1
353            int n = thisCut->row().getNumElements();
354            printf("Cut %d has %d entries, rhs %g %g =>", i, n, thisCut->lb(),
355                   thisCut->ub());
356            int j;
357            const int * index = thisCut->row().getIndices();
358            const double * element = thisCut->row().getElements();
359            for (j = 0; j < n; j++) {
360                printf(" (%d,%g)", index[j], element[j]);
361                assert(fabs(element[j]) > 1.00e-12);
362            }
363            printf("\n");
364#else
365            int n = thisCut->row().getNumElements();
366            int j;
367            const double * element = thisCut->row().getElements();
368            for (j = 0; j < n; j++) {
369                assert(fabs(element[j]) > 1.00e-12);
370            }
371#endif
372#endif
373        }
374    }
375}
376
377void
378CbcNodeInfo::addCuts(int numberCuts, CbcCountRowCut ** cut,
379                     int numberToBranchOn)
380{
381    if (numberCuts) {
382        int i;
383        if (!numberCuts_) {
384            cuts_ = new CbcCountRowCut * [numberCuts];
385        } else {
386            CbcCountRowCut ** temp = new CbcCountRowCut * [numberCuts+numberCuts_];
387            memcpy(temp, cuts_, numberCuts_*sizeof(CbcCountRowCut *));
388            delete [] cuts_;
389            cuts_ = temp;
390        }
391        for (i = 0; i < numberCuts; i++) {
392            CbcCountRowCut * thisCut = cut[i];
393            thisCut->setInfo(this, numberCuts_);
394            //printf("info %x cut %d %x\n",this,i,thisCut);
395            thisCut->increment(numberToBranchOn);
396            cuts_[numberCuts_++] = thisCut;
397#ifdef CBC_DEBUG
398            int n = thisCut->row().getNumElements();
399#if CBC_DEBUG>1
400            printf("Cut %d has %d entries, rhs %g %g =>", i, n, thisCut->lb(),
401                   thisCut->ub());
402#endif
403            int j;
404#if CBC_DEBUG>1
405            const int * index = thisCut->row().getIndices();
406#endif
407            const double * element = thisCut->row().getElements();
408            for (j = 0; j < n; j++) {
409#if CBC_DEBUG>1
410                printf(" (%d,%g)", index[j], element[j]);
411#endif
412                assert(fabs(element[j]) > 1.00e-12);
413            }
414            printf("\n");
415#endif
416        }
417    }
418}
419
420// delete cuts
421void
422CbcNodeInfo::deleteCuts(int numberToDelete, CbcCountRowCut ** cuts)
423{
424    int i;
425    int j;
426    int last = -1;
427    for (i = 0; i < numberToDelete; i++) {
428        CbcCountRowCut * next = cuts[i];
429        for (j = last + 1; j < numberCuts_; j++) {
430            if (next == cuts_[j])
431                break;
432        }
433        if (j == numberCuts_) {
434            // start from beginning
435            for (j = 0; j < last; j++) {
436                if (next == cuts_[j])
437                    break;
438            }
439            assert(j < last);
440        }
441        last = j;
442        int number = cuts_[j]->decrement();
443        if (!number) {
444#ifndef GLOBAL_CUTS_JUST_POINTERS
445            delete cuts_[j];
446#else
447            if (cuts_[j]->globallyValidAsInteger() != 2)
448                delete cuts_[j];
449#endif
450        }
451        cuts_[j] = NULL;
452    }
453    j = 0;
454    for (i = 0; i < numberCuts_; i++) {
455        if (cuts_[i])
456            cuts_[j++] = cuts_[i];
457    }
458    numberCuts_ = j;
459}
460
461// delete cuts
462void
463CbcNodeInfo::deleteCuts(int numberToDelete, int * which)
464{
465    int i;
466    for (i = 0; i < numberToDelete; i++) {
467        int iCut = which[i];
468        int number = cuts_[iCut]->decrement();
469        if (!number) {
470#ifndef GLOBAL_CUTS_JUST_POINTERS
471            delete cuts_[iCut];
472#else
473            if (cuts_[iCut]->globallyValidAsInteger() != 2)
474                delete cuts_[iCut];
475#endif
476        }
477        cuts_[iCut] = NULL;
478    }
479    int j = 0;
480    for (i = 0; i < numberCuts_; i++) {
481        if (cuts_[i])
482            cuts_[j++] = cuts_[i];
483    }
484    numberCuts_ = j;
485}
486
487// Really delete a cut
488void
489CbcNodeInfo::deleteCut(int whichOne)
490{
491    assert(whichOne < numberCuts_);
492    cuts_[whichOne] = NULL;
493}
494/* Deactivate node information.
495   1 - bounds
496   2 - cuts
497   4 - basis!
498*/
499void
500CbcNodeInfo::deactivate(int mode)
501{
502    active_ &= (~mode);
503}
504
505CbcFullNodeInfo::CbcFullNodeInfo() :
506        CbcNodeInfo(),
507        basis_(),
508        numberIntegers_(0),
509        lower_(NULL),
510        upper_(NULL)
511{
512}
513CbcFullNodeInfo::CbcFullNodeInfo(CbcModel * model,
514                                 int numberRowsAtContinuous) :
515        CbcNodeInfo(NULL, model->currentNode())
516{
517    OsiSolverInterface * solver = model->solver();
518    numberRows_ = numberRowsAtContinuous;
519    numberIntegers_ = model->numberIntegers();
520    int numberColumns = model->getNumCols();
521    lower_ = new double [numberColumns];
522    upper_ = new double [numberColumns];
523    const double * lower = solver->getColLower();
524    const double * upper = solver->getColUpper();
525    int i;
526
527    for (i = 0; i < numberColumns; i++) {
528        lower_[i] = lower[i];
529        upper_[i] = upper[i];
530    }
531
532    basis_ =  dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
533}
534
535CbcFullNodeInfo::CbcFullNodeInfo(const CbcFullNodeInfo & rhs) :
536        CbcNodeInfo(rhs)
537{
538    basis_ = dynamic_cast<CoinWarmStartBasis *>(rhs.basis_->clone()) ;
539    numberIntegers_ = rhs.numberIntegers_;
540    lower_ = NULL;
541    upper_ = NULL;
542    if (rhs.lower_ != NULL) {
543        int numberColumns = basis_->getNumStructural();
544        lower_ = new double [numberColumns];
545        upper_ = new double [numberColumns];
546        assert (upper_ != NULL);
547        memcpy(lower_, rhs.lower_, numberColumns*sizeof(double));
548        memcpy(upper_, rhs.upper_, numberColumns*sizeof(double));
549    }
550}
551
552CbcNodeInfo *
553CbcFullNodeInfo::clone() const
554{
555    return (new CbcFullNodeInfo(*this));
556}
557
558CbcFullNodeInfo::~CbcFullNodeInfo ()
559{
560    delete basis_ ;
561    delete [] lower_;
562    delete [] upper_;
563}
564
565/*
566  The basis supplied as a parameter is deleted and replaced with a new basis
567  appropriate for the node, and lower and upper bounds on variables are
568  reset according to the stored bounds arrays. Any cuts associated with this
569  node are added to the list in addCuts, but not actually added to the
570  constraint system in the model.
571
572  Why pass in a basis at all? The short answer is ``We need the parameter to
573  pass out a basis, so might as well use it to pass in the size.''
574
575  A longer answer is that in practice we take a memory allocation hit up in
576  addCuts1 (the only place applyToModel is called) when we setSize() the
577  basis that's passed in. It's immediately tossed here in favour of a clone
578  of the basis attached to this nodeInfo. This can probably be fixed, given
579  a bit of thought.
580*/
581
582void CbcFullNodeInfo::applyToModel (CbcModel *model,
583                                    CoinWarmStartBasis *&basis,
584                                    CbcCountRowCut **addCuts,
585                                    int &currentNumberCuts) const
586
587{
588    OsiSolverInterface *solver = model->solver() ;
589
590    // branch - do bounds
591    assert (active_ == 7 || active_ == 15);
592    int i;
593    solver->setColLower(lower_);
594    solver->setColUpper(upper_);
595    int numberColumns = model->getNumCols();
596    // move basis - but make sure size stays
597    // for bon-min - should not be needed int numberRows = model->getNumRows();
598    int numberRows = basis->getNumArtificial();
599    delete basis ;
600    if (basis_) {
601        basis = dynamic_cast<CoinWarmStartBasis *>(basis_->clone()) ;
602        basis->resize(numberRows, numberColumns);
603#ifdef CBC_CHECK_BASIS
604        std::cout << "Basis (after applying root " << this << ") " << std::endl ;
605        basis->print() ;
606#endif
607    } else {
608        // We have a solver without a basis
609        basis = NULL;
610    }
611    for (i = 0; i < numberCuts_; i++)
612        addCuts[currentNumberCuts+i] = cuts_[i];
613    currentNumberCuts += numberCuts_;
614    assert(!parent_);
615    return ;
616}
617// Just apply bounds to one variable (1=>infeasible)
618int
619CbcFullNodeInfo::applyBounds(int iColumn, double & lower, double & upper, int force)
620{
621    if ((force && 1) == 0) {
622        if (lower > lower_[iColumn])
623            printf("%d odd lower going from %g to %g\n", iColumn, lower, lower_[iColumn]);
624        lower = lower_[iColumn];
625    } else {
626        lower_[iColumn] = lower;
627    }
628    if ((force && 2) == 0) {
629        if (upper < upper_[iColumn])
630            printf("%d odd upper going from %g to %g\n", iColumn, upper, upper_[iColumn]);
631        upper = upper_[iColumn];
632    } else {
633        upper_[iColumn] = upper;
634    }
635    return (upper_[iColumn] >= lower_[iColumn]) ? 0 : 1;
636}
637
638/* Builds up row basis backwards (until original model).
639   Returns NULL or previous one to apply .
640   Depends on Free being 0 and impossible for cuts
641*/
642CbcNodeInfo *
643CbcFullNodeInfo::buildRowBasis(CoinWarmStartBasis & basis ) const
644{
645    const unsigned int * saved =
646        reinterpret_cast<const unsigned int *> (basis_->getArtificialStatus());
647    unsigned int * now =
648        reinterpret_cast<unsigned int *> (basis.getArtificialStatus());
649    int number = basis_->getNumArtificial() >> 4;;
650    int i;
651    for (i = 0; i < number; i++) {
652        if (!now[i])
653            now[i] = saved[i];
654    }
655    return NULL;
656}
657
658
659// Default constructor
660CbcPartialNodeInfo::CbcPartialNodeInfo()
661
662        : CbcNodeInfo(),
663        basisDiff_(NULL),
664        variables_(NULL),
665        newBounds_(NULL),
666        numberChangedBounds_(0)
667
668{ /* this space intentionally left blank */ }
669
670// Constructor from current state
671CbcPartialNodeInfo::CbcPartialNodeInfo (CbcNodeInfo *parent, CbcNode *owner,
672                                        int numberChangedBounds,
673                                        const int *variables,
674                                        const double *boundChanges,
675                                        const CoinWarmStartDiff *basisDiff)
676        : CbcNodeInfo(parent, owner)
677{
678    basisDiff_ = basisDiff->clone() ;
679#ifdef CBC_CHECK_BASIS
680    std::cout << "Constructor (" << this << ") " << std::endl ;
681#endif
682
683    numberChangedBounds_ = numberChangedBounds;
684    int size = numberChangedBounds_ * (sizeof(double) + sizeof(int));
685    char * temp = new char [size];
686    newBounds_ = reinterpret_cast<double *> (temp);
687    variables_ = reinterpret_cast<int *> (newBounds_ + numberChangedBounds_);
688
689    int i ;
690    for (i = 0; i < numberChangedBounds_; i++) {
691        variables_[i] = variables[i];
692        newBounds_[i] = boundChanges[i];
693    }
694}
695
696CbcPartialNodeInfo::CbcPartialNodeInfo (const CbcPartialNodeInfo & rhs)
697
698        : CbcNodeInfo(rhs)
699
700{
701    basisDiff_ = rhs.basisDiff_->clone() ;
702
703#ifdef CBC_CHECK_BASIS
704    std::cout << "Copy constructor (" << this << ") from " << this << std::endl ;
705#endif
706    numberChangedBounds_ = rhs.numberChangedBounds_;
707    int size = numberChangedBounds_ * (sizeof(double) + sizeof(int));
708    char * temp = new char [size];
709    newBounds_ = reinterpret_cast<double *> (temp);
710    variables_ = reinterpret_cast<int *> (newBounds_ + numberChangedBounds_);
711
712    int i ;
713    for (i = 0; i < numberChangedBounds_; i++) {
714        variables_[i] = rhs.variables_[i];
715        newBounds_[i] = rhs.newBounds_[i];
716    }
717}
718
719CbcNodeInfo *
720CbcPartialNodeInfo::clone() const
721{
722    return (new CbcPartialNodeInfo(*this));
723}
724
725
726CbcPartialNodeInfo::~CbcPartialNodeInfo ()
727{
728    delete basisDiff_ ;
729    delete [] newBounds_;
730}
731
732
733/**
734   The basis supplied as a parameter is incrementally modified, and lower and
735   upper bounds on variables in the model are incrementally modified. Any
736   cuts associated with this node are added to the list in addCuts.
737*/
738
739void CbcPartialNodeInfo::applyToModel (CbcModel *model,
740                                       CoinWarmStartBasis *&basis,
741                                       CbcCountRowCut **addCuts,
742                                       int &currentNumberCuts) const
743
744{
745    OsiSolverInterface *solver = model->solver();
746    if ((active_&4) != 0) {
747        basis->applyDiff(basisDiff_) ;
748#ifdef CBC_CHECK_BASIS
749        std::cout << "Basis (after applying " << this << ") " << std::endl ;
750        basis->print() ;
751#endif
752    }
753
754    // branch - do bounds
755    int i;
756    if ((active_&1) != 0) {
757        for (i = 0; i < numberChangedBounds_; i++) {
758            int variable = variables_[i];
759            int k = variable & 0x3fffffff;
760            if ((variable&0x80000000) == 0) {
761                // lower bound changing
762                //#define CBC_PRINT2
763#ifdef CBC_PRINT2
764                if (solver->getColLower()[k] != newBounds_[i])
765                    printf("lower change for column %d - from %g to %g\n", k, solver->getColLower()[k], newBounds_[i]);
766#endif
767#ifndef NDEBUG
768                if ((variable&0x40000000) == 0 && false) {
769                    double oldValue = solver->getColLower()[k];
770                    assert (newBounds_[i] > oldValue - 1.0e-8);
771                    if (newBounds_[i] < oldValue + 1.0e-8)
772                        printf("bad null lower change for column %d - bound %g\n", k, oldValue);
773                }
774#endif
775                solver->setColLower(k, newBounds_[i]);
776            } else {
777                // upper bound changing
778#ifdef CBC_PRINT2
779                if (solver->getColUpper()[k] != newBounds_[i])
780                    printf("upper change for column %d - from %g to %g\n", k, solver->getColUpper()[k], newBounds_[i]);
781#endif
782#ifndef NDEBUG
783                if ((variable&0x40000000) == 0 && false) {
784                    double oldValue = solver->getColUpper()[k];
785                    assert (newBounds_[i] < oldValue + 1.0e-8);
786                    if (newBounds_[i] > oldValue - 1.0e-8)
787                        printf("bad null upper change for column %d - bound %g\n", k, oldValue);
788                }
789#endif
790                solver->setColUpper(k, newBounds_[i]);
791            }
792        }
793    }
794    if ((active_&2) != 0) {
795        for (i = 0; i < numberCuts_; i++) {
796            addCuts[currentNumberCuts+i] = cuts_[i];
797            if (cuts_[i] && model->messageHandler()->logLevel() > 4) {
798                cuts_[i]->print();
799            }
800        }
801
802        currentNumberCuts += numberCuts_;
803    }
804    return ;
805}
806// Just apply bounds to one variable (1=>infeasible)
807int
808CbcPartialNodeInfo::applyBounds(int iColumn, double & lower, double & upper, int force)
809{
810    // branch - do bounds
811    int i;
812    int found = 0;
813    double newLower = -COIN_DBL_MAX;
814    double newUpper = COIN_DBL_MAX;
815    for (i = 0; i < numberChangedBounds_; i++) {
816        int variable = variables_[i];
817        int k = variable & 0x3fffffff;
818        if (k == iColumn) {
819            if ((variable&0x80000000) == 0) {
820                // lower bound changing
821                found |= 1;
822                newLower = CoinMax(newLower, newBounds_[i]);
823                if ((force&1) == 0) {
824                    if (lower > newBounds_[i])
825                        printf("%d odd lower going from %g to %g\n", iColumn, lower, newBounds_[i]);
826                    lower = newBounds_[i];
827                } else {
828                    newBounds_[i] = lower;
829                    variables_[i] |= 0x40000000; // say can go odd way
830                }
831            } else {
832                // upper bound changing
833                found |= 2;
834                newUpper = CoinMin(newUpper, newBounds_[i]);
835                if ((force&2) == 0) {
836                    if (upper < newBounds_[i])
837                        printf("%d odd upper going from %g to %g\n", iColumn, upper, newBounds_[i]);
838                    upper = newBounds_[i];
839                } else {
840                    newBounds_[i] = upper;
841                    variables_[i] |= 0x40000000; // say can go odd way
842                }
843            }
844        }
845    }
846    newLower = CoinMax(newLower, lower);
847    newUpper = CoinMin(newUpper, upper);
848    int nAdd = 0;
849    if ((force&2) != 0 && (found&2) == 0) {
850        // need to add new upper
851        nAdd++;
852    }
853    if ((force&1) != 0 && (found&1) == 0) {
854        // need to add new lower
855        nAdd++;
856    }
857    if (nAdd) {
858        int size = (numberChangedBounds_ + nAdd) * (sizeof(double) + sizeof(int));
859        char * temp = new char [size];
860        double * newBounds = reinterpret_cast<double *> (temp);
861        int * variables = reinterpret_cast<int *> (newBounds + numberChangedBounds_ + nAdd);
862
863        int i ;
864        for (i = 0; i < numberChangedBounds_; i++) {
865            variables[i] = variables_[i];
866            newBounds[i] = newBounds_[i];
867        }
868        delete [] newBounds_;
869        newBounds_ = newBounds;
870        variables_ = variables;
871        if ((force&2) != 0 && (found&2) == 0) {
872            // need to add new upper
873            int variable = iColumn | 0x80000000;
874            variables_[numberChangedBounds_] = variable;
875            newBounds_[numberChangedBounds_++] = newUpper;
876        }
877        if ((force&1) != 0 && (found&1) == 0) {
878            // need to add new lower
879            int variable = iColumn;
880            variables_[numberChangedBounds_] = variable;
881            newBounds_[numberChangedBounds_++] = newLower;
882        }
883    }
884
885    return (newUpper >= newLower) ? 0 : 1;
886}
887
888/* Builds up row basis backwards (until original model).
889   Returns NULL or previous one to apply .
890   Depends on Free being 0 and impossible for cuts
891*/
892
893CbcNodeInfo *
894CbcPartialNodeInfo::buildRowBasis(CoinWarmStartBasis & basis ) const
895
896{
897    basis.applyDiff(basisDiff_) ;
898
899    return parent_ ;
900}
901CbcNode::CbcNode() :
902        nodeInfo_(NULL),
903        objectiveValue_(1.0e100),
904        guessedObjectiveValue_(1.0e100),
905        sumInfeasibilities_(0.0),
906        branch_(NULL),
907        depth_(-1),
908        numberUnsatisfied_(0),
909        nodeNumber_(-1),
910        state_(0)
911{
912#ifdef CHECK_NODE
913    printf("CbcNode %x Constructor\n", this);
914#endif
915}
916// Print
917void
918CbcNode::print() const
919{
920    printf("number %d obj %g depth %d sumun %g nunsat %d state %d\n",
921           nodeNumber_, objectiveValue_, depth_, sumInfeasibilities_, numberUnsatisfied_, state_);
922}
923CbcNode::CbcNode(CbcModel * model,
924                 CbcNode * lastNode) :
925        nodeInfo_(NULL),
926        objectiveValue_(1.0e100),
927        guessedObjectiveValue_(1.0e100),
928        sumInfeasibilities_(0.0),
929        branch_(NULL),
930        depth_(-1),
931        numberUnsatisfied_(0),
932        nodeNumber_(-1),
933        state_(0)
934{
935#ifdef CHECK_NODE
936    printf("CbcNode %x Constructor from model\n", this);
937#endif
938    model->setObjectiveValue(this, lastNode);
939
940    if (lastNode) {
941        if (lastNode->nodeInfo_) {
942            lastNode->nodeInfo_->increment();
943        }
944    }
945    nodeNumber_ = model->getNodeCount();
946}
947
948#define CBC_NEW_CREATEINFO
949#ifdef CBC_NEW_CREATEINFO
950
951/*
952  New createInfo, with basis manipulation hidden inside mergeBasis. Allows
953  solvers to override and carry over all information from one basis to
954  another.
955*/
956
957void
958CbcNode::createInfo (CbcModel *model,
959                     CbcNode *lastNode,
960                     const CoinWarmStartBasis *lastws,
961                     const double *lastLower, const double *lastUpper,
962                     int numberOldActiveCuts, int numberNewCuts)
963
964{
965    OsiSolverInterface *solver = model->solver();
966    CbcStrategy *strategy = model->strategy();
967    /*
968      The root --- no parent. Create full basis and bounds information.
969    */
970    if (!lastNode) {
971        if (!strategy)
972            nodeInfo_ = new CbcFullNodeInfo(model, solver->getNumRows());
973        else
974            nodeInfo_ = strategy->fullNodeInfo(model, solver->getNumRows());
975    } else {
976        /*
977          Not the root. Create an edit from the parent's basis & bound information.
978          This is not quite as straightforward as it seems. We need to reintroduce
979          cuts we may have dropped out of the basis, in the correct position, because
980          this whole process is strictly positional. Start by grabbing the current
981          basis.
982        */
983        bool mustDeleteBasis;
984        const CoinWarmStartBasis *ws =
985            dynamic_cast<const CoinWarmStartBasis*>(solver->getPointerToWarmStart(mustDeleteBasis));
986        assert(ws != NULL); // make sure not volume
987        //int numberArtificials = lastws->getNumArtificial();
988        int numberColumns = solver->getNumCols();
989        int numberRowsAtContinuous = model->numberRowsAtContinuous();
990        int currentNumberCuts = model->currentNumberCuts();
991#   ifdef CBC_CHECK_BASIS
992        std::cout
993            << "Before expansion: orig " << numberRowsAtContinuous
994            << ", old " << numberOldActiveCuts
995            << ", new " << numberNewCuts
996            << ", current " << currentNumberCuts << "." << std::endl ;
997        ws->print();
998#   endif
999        /*
1000          Clone the basis and resize it to hold the structural constraints, plus
1001          all the cuts: old cuts, both active and inactive (currentNumberCuts),
1002          and new cuts (numberNewCuts). This will become the expanded basis.
1003        */
1004        CoinWarmStartBasis *expanded =
1005            dynamic_cast<CoinWarmStartBasis *>(ws->clone()) ;
1006        int iCompact = numberRowsAtContinuous + numberOldActiveCuts + numberNewCuts ;
1007        // int nPartial = numberRowsAtContinuous+currentNumberCuts;
1008        int iFull = numberRowsAtContinuous + currentNumberCuts + numberNewCuts;
1009        // int maxBasisLength = ((iFull+15)>>4)+((numberColumns+15)>>4);
1010        // printf("l %d full %d\n",maxBasisLength,iFull);
1011        expanded->resize(iFull, numberColumns);
1012#   ifdef CBC_CHECK_BASIS
1013        std::cout
1014            << "\tFull basis " << iFull << " rows, "
1015            << numberColumns << " columns; compact "
1016            << iCompact << " rows." << std::endl ;
1017#   endif
1018        /*
1019          Now flesh out the expanded basis. The clone already has the
1020          correct status information for the variables and for the structural
1021          (numberRowsAtContinuous) constraints. Any indices beyond nPartial must be
1022          cuts created while processing this node --- they can be copied en bloc
1023          into the correct position in the expanded basis. The space reserved for
1024          xferRows is a gross overestimate.
1025        */
1026        CoinWarmStartBasis::XferVec xferRows ;
1027        xferRows.reserve(iFull - numberRowsAtContinuous + 1) ;
1028        if (numberNewCuts) {
1029            xferRows.push_back(
1030                CoinWarmStartBasis::XferEntry(iCompact - numberNewCuts,
1031                                              iFull - numberNewCuts, numberNewCuts)) ;
1032        }
1033        /*
1034          From nPartial down, record the entries we want to copy from the current
1035          basis (the entries for the active cuts; non-zero in the list returned
1036          by addedCuts). Fill the expanded basis with entries showing a status of
1037          basic for the deactivated (loose) cuts.
1038        */
1039        CbcCountRowCut **cut = model->addedCuts();
1040        iFull -= (numberNewCuts + 1) ;
1041        iCompact -= (numberNewCuts + 1) ;
1042        int runLen = 0 ;
1043        CoinWarmStartBasis::XferEntry entry(-1, -1, -1) ;
1044        while (iFull >= numberRowsAtContinuous) {
1045            for ( ; iFull >= numberRowsAtContinuous &&
1046                    cut[iFull-numberRowsAtContinuous] ; iFull--)
1047                runLen++ ;
1048            if (runLen) {
1049                iCompact -= runLen ;
1050                entry.first = iCompact + 1 ;
1051                entry.second = iFull + 1 ;
1052                entry.third = runLen ;
1053                runLen = 0 ;
1054                xferRows.push_back(entry) ;
1055            }
1056            for ( ; iFull >= numberRowsAtContinuous &&
1057                    !cut[iFull-numberRowsAtContinuous] ; iFull--)
1058                expanded->setArtifStatus(iFull, CoinWarmStartBasis::basic);
1059        }
1060        /*
1061          Finally, call mergeBasis to copy over entries from the current basis to
1062          the expanded basis. Since we cloned the expanded basis from the active basis
1063          and haven't changed the number of variables, only row status entries need
1064          to be copied.
1065        */
1066        expanded->mergeBasis(ws, &xferRows, 0) ;
1067
1068#ifdef CBC_CHECK_BASIS
1069        std::cout << "Expanded basis:" << std::endl ;
1070        expanded->print() ;
1071        std::cout << "Diffing against:" << std::endl ;
1072        lastws->print() ;
1073#endif
1074        assert (expanded->getNumArtificial() >= lastws->getNumArtificial());
1075#ifdef CLP_INVESTIGATE
1076        if (!expanded->fullBasis()) {
1077            int iFull = numberRowsAtContinuous + currentNumberCuts + numberNewCuts;
1078            printf("cont %d old %d new %d current %d full inc %d full %d\n",
1079                   numberRowsAtContinuous, numberOldActiveCuts, numberNewCuts,
1080                   currentNumberCuts, iFull, iFull - numberNewCuts);
1081        }
1082#endif
1083
1084        /*
1085          Now that we have two bases in proper positional correspondence, creating
1086          the actual diff is dead easy.
1087
1088          Note that we're going to compare the expanded basis here to the stripped
1089          basis (lastws) produced by addCuts. It doesn't affect the correctness (the
1090          diff process has no knowledge of the meaning of an entry) but it does
1091          mean that we'll always generate a whack of diff entries because the expanded
1092          basis is considerably larger than the stripped basis.
1093        */
1094        CoinWarmStartDiff *basisDiff = expanded->generateDiff(lastws) ;
1095
1096        /*
1097          Diff the bound vectors. It's assumed the number of structural variables
1098          is not changing. For branching objects that change bounds on integer
1099          variables, we should see at least one bound change as a consequence
1100          of applying the branch that generated this subproblem from its parent.
1101          This need not hold for other types of branching objects (hyperplane
1102          branches, for example).
1103        */
1104        const double * lower = solver->getColLower();
1105        const double * upper = solver->getColUpper();
1106
1107        double *boundChanges = new double [2*numberColumns] ;
1108        int *variables = new int [2*numberColumns] ;
1109        int numberChangedBounds = 0;
1110
1111        int i;
1112        for (i = 0; i < numberColumns; i++) {
1113            if (lower[i] != lastLower[i]) {
1114                variables[numberChangedBounds] = i;
1115                boundChanges[numberChangedBounds++] = lower[i];
1116            }
1117            if (upper[i] != lastUpper[i]) {
1118                variables[numberChangedBounds] = i | 0x80000000;
1119                boundChanges[numberChangedBounds++] = upper[i];
1120            }
1121#ifdef CBC_DEBUG
1122            if (lower[i] != lastLower[i]) {
1123                std::cout
1124                    << "lower on " << i << " changed from "
1125                    << lastLower[i] << " to " << lower[i] << std::endl ;
1126            }
1127            if (upper[i] != lastUpper[i]) {
1128                std::cout
1129                    << "upper on " << i << " changed from "
1130                    << lastUpper[i] << " to " << upper[i] << std::endl ;
1131            }
1132#endif
1133        }
1134#ifdef CBC_DEBUG
1135        std::cout << numberChangedBounds << " changed bounds." << std::endl ;
1136#endif
1137        //if (lastNode->branchingObject()->boundBranch())
1138        //assert (numberChangedBounds);
1139        /*
1140          Hand the lot over to the CbcPartialNodeInfo constructor, then clean up and
1141          return.
1142        */
1143        if (!strategy)
1144            nodeInfo_ =
1145                new CbcPartialNodeInfo(lastNode->nodeInfo_, this, numberChangedBounds,
1146                                       variables, boundChanges, basisDiff) ;
1147        else
1148            nodeInfo_ =
1149                strategy->partialNodeInfo(model, lastNode->nodeInfo_, this,
1150                                          numberChangedBounds, variables, boundChanges,
1151                                          basisDiff) ;
1152        delete basisDiff ;
1153        delete [] boundChanges;
1154        delete [] variables;
1155        delete expanded ;
1156        if  (mustDeleteBasis)
1157            delete ws;
1158    }
1159    // Set node number
1160    nodeInfo_->setNodeNumber(model->getNodeCount2());
1161    state_ |= 2; // say active
1162}
1163
1164#else   // CBC_NEW_CREATEINFO
1165
1166/*
1167  Original createInfo, with bare manipulation of basis vectors. Fails if solver
1168  maintains additional information in basis.
1169*/
1170
1171void
1172CbcNode::createInfo (CbcModel *model,
1173                     CbcNode *lastNode,
1174                     const CoinWarmStartBasis *lastws,
1175                     const double *lastLower, const double *lastUpper,
1176                     int numberOldActiveCuts, int numberNewCuts)
1177{
1178    OsiSolverInterface * solver = model->solver();
1179    CbcStrategy * strategy = model->strategy();
1180    /*
1181      The root --- no parent. Create full basis and bounds information.
1182    */
1183    if (!lastNode) {
1184        if (!strategy)
1185            nodeInfo_ = new CbcFullNodeInfo(model, solver->getNumRows());
1186        else
1187            nodeInfo_ = strategy->fullNodeInfo(model, solver->getNumRows());
1188    }
1189    /*
1190      Not the root. Create an edit from the parent's basis & bound information.
1191      This is not quite as straightforward as it seems. We need to reintroduce
1192      cuts we may have dropped out of the basis, in the correct position, because
1193      this whole process is strictly positional. Start by grabbing the current
1194      basis.
1195    */
1196    else {
1197        bool mustDeleteBasis;
1198        const CoinWarmStartBasis* ws =
1199            dynamic_cast<const CoinWarmStartBasis*>(solver->getPointerToWarmStart(mustDeleteBasis));
1200        assert(ws != NULL); // make sure not volume
1201        //int numberArtificials = lastws->getNumArtificial();
1202        int numberColumns = solver->getNumCols();
1203
1204        const double * lower = solver->getColLower();
1205        const double * upper = solver->getColUpper();
1206
1207        int i;
1208        /*
1209        Create a clone and resize it to hold all the structural constraints, plus
1210        all the cuts: old cuts, both active and inactive (currentNumberCuts), and
1211        new cuts (numberNewCuts).
1212
1213        TODO: You'd think that the set of constraints (logicals) in the expanded
1214        basis should match the set represented in lastws. At least, that's
1215        what I thought. But at the point I first looked hard at this bit of
1216        code, it turned out that lastws was the stripped basis produced at
1217        the end of addCuts(), rather than the raw basis handed back by
1218        addCuts1(). The expanded basis here is equivalent to the raw basis of
1219        addCuts1(). I said ``whoa, that's not good, I must have introduced a
1220        bug'' and went back to John's code to see where I'd gone wrong.
1221        And discovered the same `error' in his code.
1222
1223        After a bit of thought, my conclusion is that correctness is not
1224        affected by whether lastws is the stripped or raw basis. The diffs
1225        have no semantics --- just a set of changes that need to be made
1226        to convert lastws into expanded. I think the only effect is that we
1227        store a lot more diffs (everything in expanded that's not covered by
1228        the stripped basis). But I need to give this more thought. There
1229        may well be some subtle error cases.
1230
1231        In the mean time, I've twiddled addCuts() to set lastws to the raw
1232        basis. Makes me (Lou) less nervous to compare apples to apples.
1233        */
1234        CoinWarmStartBasis *expanded =
1235            dynamic_cast<CoinWarmStartBasis *>(ws->clone()) ;
1236        int numberRowsAtContinuous = model->numberRowsAtContinuous();
1237        int iFull = numberRowsAtContinuous + model->currentNumberCuts() +
1238                    numberNewCuts;
1239        //int numberArtificialsNow = iFull;
1240        //int maxBasisLength = ((iFull+15)>>4)+((numberColumns+15)>>4);
1241        //printf("l %d full %d\n",maxBasisLength,iFull);
1242        if (expanded)
1243            expanded->resize(iFull, numberColumns);
1244#ifdef CBC_CHECK_BASIS
1245        printf("Before expansion: orig %d, old %d, new %d, current %d\n",
1246               numberRowsAtContinuous, numberOldActiveCuts, numberNewCuts,
1247               model->currentNumberCuts()) ;
1248        ws->print();
1249#endif
1250        /*
1251        Now fill in the expanded basis. Any indices beyond nPartial must
1252        be cuts created while processing this node --- they can be copied directly
1253        into the expanded basis. From nPartial down, pull the status of active cuts
1254        from ws, interleaving with a B entry for the deactivated (loose) cuts.
1255        */
1256        int numberDropped = model->currentNumberCuts() - numberOldActiveCuts;
1257        int iCompact = iFull - numberDropped;
1258        CbcCountRowCut ** cut = model->addedCuts();
1259        int nPartial = model->currentNumberCuts() + numberRowsAtContinuous;
1260        iFull--;
1261        for (; iFull >= nPartial; iFull--) {
1262            CoinWarmStartBasis::Status status = ws->getArtifStatus(--iCompact);
1263            //assert (status != CoinWarmStartBasis::basic); // may be permanent cut
1264            expanded->setArtifStatus(iFull, status);
1265        }
1266        for (; iFull >= numberRowsAtContinuous; iFull--) {
1267            if (cut[iFull-numberRowsAtContinuous]) {
1268                CoinWarmStartBasis::Status status = ws->getArtifStatus(--iCompact);
1269                // If no cut generator being used then we may have basic variables
1270                //if (model->getMaximumCutPasses()&&
1271                //  status == CoinWarmStartBasis::basic)
1272                //printf("cut basic\n");
1273                expanded->setArtifStatus(iFull, status);
1274            } else {
1275                expanded->setArtifStatus(iFull, CoinWarmStartBasis::basic);
1276            }
1277        }
1278#ifdef CBC_CHECK_BASIS
1279        printf("Expanded basis\n");
1280        expanded->print() ;
1281        printf("Diffing against\n") ;
1282        lastws->print() ;
1283#endif
1284        /*
1285        Now that we have two bases in proper positional correspondence, creating
1286        the actual diff is dead easy.
1287        */
1288
1289        CoinWarmStartDiff *basisDiff = expanded->generateDiff(lastws) ;
1290        /*
1291        Diff the bound vectors. It's assumed the number of structural variables is
1292        not changing. Assuming that branching objects all involve integer variables,
1293        we should see at least one bound change as a consequence of processing this
1294        subproblem. Different types of branching objects could break this assertion.
1295        Not true at all - we have not applied current branch - JJF.
1296        */
1297        double *boundChanges = new double [2*numberColumns] ;
1298        int *variables = new int [2*numberColumns] ;
1299        int numberChangedBounds = 0;
1300        for (i = 0; i < numberColumns; i++) {
1301            if (lower[i] != lastLower[i]) {
1302                variables[numberChangedBounds] = i;
1303                boundChanges[numberChangedBounds++] = lower[i];
1304            }
1305            if (upper[i] != lastUpper[i]) {
1306                variables[numberChangedBounds] = i | 0x80000000;
1307                boundChanges[numberChangedBounds++] = upper[i];
1308            }
1309#ifdef CBC_DEBUG
1310            if (lower[i] != lastLower[i])
1311                printf("lower on %d changed from %g to %g\n",
1312                       i, lastLower[i], lower[i]);
1313            if (upper[i] != lastUpper[i])
1314                printf("upper on %d changed from %g to %g\n",
1315                       i, lastUpper[i], upper[i]);
1316#endif
1317        }
1318#ifdef CBC_DEBUG
1319        printf("%d changed bounds\n", numberChangedBounds) ;
1320#endif
1321        //if (lastNode->branchingObject()->boundBranch())
1322        //assert (numberChangedBounds);
1323        /*
1324        Hand the lot over to the CbcPartialNodeInfo constructor, then clean up and
1325        return.
1326        */
1327        if (!strategy)
1328            nodeInfo_ =
1329                new CbcPartialNodeInfo(lastNode->nodeInfo_, this, numberChangedBounds,
1330                                       variables, boundChanges, basisDiff) ;
1331        else
1332            nodeInfo_ = strategy->partialNodeInfo(model, lastNode->nodeInfo_, this, numberChangedBounds,
1333                                                  variables, boundChanges, basisDiff) ;
1334        delete basisDiff ;
1335        delete [] boundChanges;
1336        delete [] variables;
1337        delete expanded ;
1338        if  (mustDeleteBasis)
1339            delete ws;
1340    }
1341    // Set node number
1342    nodeInfo_->setNodeNumber(model->getNodeCount2());
1343    state_ |= 2; // say active
1344}
1345
1346#endif  // CBC_NEW_CREATEINFO
1347/*
1348  The routine scans through the object list of the model looking for objects
1349  that indicate infeasibility. It tests each object using strong branching
1350  and selects the one with the least objective degradation.  A corresponding
1351  branching object is left attached to lastNode.
1352
1353  If strong branching is disabled, a candidate object is chosen essentially
1354  at random (whatever object ends up in pos'n 0 of the candidate array).
1355
1356  If a branching candidate is found to be monotone, bounds are set to fix the
1357  variable and the routine immediately returns (the caller is expected to
1358  reoptimize).
1359
1360  If a branching candidate is found to result in infeasibility in both
1361  directions, the routine immediately returns an indication of infeasibility.
1362
1363  Returns:  0   both branch directions are feasible
1364  -1    branching variable is monotone
1365  -2    infeasible
1366
1367  Original comments:
1368  Here could go cuts etc etc
1369  For now just fix on objective from strong branching.
1370*/
1371
1372int CbcNode::chooseBranch (CbcModel *model, CbcNode *lastNode, int numberPassesLeft)
1373
1374{
1375    if (lastNode)
1376        depth_ = lastNode->depth_ + 1;
1377    else
1378        depth_ = 0;
1379    delete branch_;
1380    branch_ = NULL;
1381    OsiSolverInterface * solver = model->solver();
1382# ifdef COIN_HAS_CLP
1383    OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
1384    int saveClpOptions = 0;
1385    if (osiclp) {
1386        // for faster hot start
1387        saveClpOptions = osiclp->specialOptions();
1388        osiclp->setSpecialOptions(saveClpOptions | 8192);
1389    }
1390# else
1391    OsiSolverInterface *osiclp = NULL ;
1392# endif
1393    double saveObjectiveValue = solver->getObjValue();
1394    double objectiveValue = CoinMax(solver->getObjSense() * saveObjectiveValue, objectiveValue_);
1395    const double * lower = solver->getColLower();
1396    const double * upper = solver->getColUpper();
1397    // See what user thinks
1398    int anyAction = model->problemFeasibility()->feasible(model, 0);
1399    if (anyAction) {
1400        // will return -2 if infeasible , 0 if treat as integer
1401        return anyAction - 1;
1402    }
1403    double integerTolerance =
1404        model->getDblParam(CbcModel::CbcIntegerTolerance);
1405    // point to useful information
1406    OsiBranchingInformation usefulInfo = model->usefulInformation();
1407    // and modify
1408    usefulInfo.depth_ = depth_;
1409    int i;
1410    bool beforeSolution = model->getSolutionCount() == 0;
1411    int numberStrong = model->numberStrong();
1412    // switch off strong if hotstart
1413    const double * hotstartSolution = model->hotstartSolution();
1414    const int * hotstartPriorities = model->hotstartPriorities();
1415    int numberObjects = model->numberObjects();
1416    int numberColumns = model->getNumCols();
1417    double * saveUpper = new double[numberColumns];
1418    double * saveLower = new double[numberColumns];
1419    for (i = 0; i < numberColumns; i++) {
1420        saveLower[i] = lower[i];
1421        saveUpper[i] = upper[i];
1422    }
1423
1424    // Save solution in case heuristics need good solution later
1425
1426    double * saveSolution = new double[numberColumns];
1427    memcpy(saveSolution, solver->getColSolution(), numberColumns*sizeof(double));
1428    model->reserveCurrentSolution(saveSolution);
1429    if (hotstartSolution) {
1430        numberStrong = 0;
1431        if ((model->moreSpecialOptions()&1024) != 0) {
1432            int nBad = 0;
1433            int nUnsat = 0;
1434            int nDiff = 0;
1435            for (int i = 0; i < numberObjects; i++) {
1436                OsiObject * object = model->modifiableObject(i);
1437                const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
1438                if (thisOne) {
1439                    int iColumn = thisOne->columnNumber();
1440                    double targetValue = hotstartSolution[iColumn];
1441                    double value = saveSolution[iColumn];
1442                    if (fabs(value - floor(value + 0.5)) > 1.0e-6) {
1443                        nUnsat++;
1444#ifdef CLP_INVESTIGATE
1445                        printf("H %d is %g target %g\n", iColumn, value, targetValue);
1446#endif
1447                    } else if (fabs(targetValue - value) > 1.0e-6) {
1448                        nDiff++;
1449                    }
1450                    if (targetValue < saveLower[iColumn] ||
1451                            targetValue > saveUpper[iColumn]) {
1452#ifdef CLP_INVESTIGATE
1453                        printf("%d has target %g and current bounds %g and %g\n",
1454                               iColumn, targetValue, saveLower[iColumn], saveUpper[iColumn]);
1455#endif
1456                        nBad++;
1457                    }
1458                }
1459            }
1460#ifdef CLP_INVESTIGATE
1461            printf("Hot %d unsatisfied, %d outside limits, %d different\n",
1462                   nUnsat, nBad, nDiff);
1463#endif
1464            if (nBad) {
1465                // switch off as not possible
1466                hotstartSolution = NULL;
1467                model->setHotstartSolution(NULL, NULL);
1468                usefulInfo.hotstartSolution_ = NULL;
1469            }
1470        }
1471    }
1472    int numberStrongDone = 0;
1473    int numberUnfinished = 0;
1474    int numberStrongInfeasible = 0;
1475    int numberStrongIterations = 0;
1476    int saveNumberStrong = numberStrong;
1477    bool checkFeasibility = numberObjects > model->numberIntegers();
1478    int maximumStrong = CoinMax(CoinMin(numberStrong, numberObjects), 1);
1479    /*
1480      Get a branching decision object. Use the default decision criteria unless
1481      the user has loaded a decision method into the model.
1482    */
1483    CbcBranchDecision *decision = model->branchingMethod();
1484    CbcDynamicPseudoCostBranchingObject * dynamicBranchingObject =
1485        dynamic_cast<CbcDynamicPseudoCostBranchingObject *>(decision);
1486    if (!decision || dynamicBranchingObject)
1487        decision = new CbcBranchDefaultDecision();
1488    decision->initialize(model);
1489    CbcStrongInfo * choice = new CbcStrongInfo[maximumStrong];
1490    // May go round twice if strong branching fixes all local candidates
1491    bool finished = false;
1492    double estimatedDegradation = 0.0;
1493    while (!finished) {
1494        finished = true;
1495        // Some objects may compute an estimate of best solution from here
1496        estimatedDegradation = 0.0;
1497        //int numberIntegerInfeasibilities=0; // without odd ones
1498        numberStrongDone = 0;
1499        numberUnfinished = 0;
1500        numberStrongInfeasible = 0;
1501        numberStrongIterations = 0;
1502
1503        // We may go round this loop twice (only if we think we have solution)
1504        for (int iPass = 0; iPass < 2; iPass++) {
1505
1506            // compute current state
1507            //int numberObjectInfeasibilities; // just odd ones
1508            //model->feasibleSolution(
1509            //                      numberIntegerInfeasibilities,
1510            //                      numberObjectInfeasibilities);
1511            // Some objects may compute an estimate of best solution from here
1512            estimatedDegradation = 0.0;
1513            numberUnsatisfied_ = 0;
1514            // initialize sum of "infeasibilities"
1515            sumInfeasibilities_ = 0.0;
1516            int bestPriority = COIN_INT_MAX;
1517            /*
1518              Scan for branching objects that indicate infeasibility. Choose the best
1519              maximumStrong candidates, using priority as the first criteria, then
1520              integer infeasibility.
1521
1522              The algorithm is to fill the choice array with a set of good candidates (by
1523              infeasibility) with priority bestPriority.  Finding a candidate with
1524              priority better (less) than bestPriority flushes the choice array. (This
1525              serves as initialization when the first candidate is found.)
1526
1527              A new candidate is added to choices only if its infeasibility exceeds the
1528              current max infeasibility (mostAway). When a candidate is added, it
1529              replaces the candidate with the smallest infeasibility (tracked by
1530              iSmallest).
1531            */
1532            int iSmallest = 0;
1533            double mostAway = 1.0e-100;
1534            for (i = 0 ; i < maximumStrong ; i++)
1535                choice[i].possibleBranch = NULL ;
1536            numberStrong = 0;
1537            bool canDoOneHot = false;
1538            for (i = 0; i < numberObjects; i++) {
1539                OsiObject * object = model->modifiableObject(i);
1540                int preferredWay;
1541                double infeasibility = object->infeasibility(&usefulInfo, preferredWay);
1542                int priorityLevel = object->priority();
1543                if (hotstartSolution) {
1544                    // we are doing hot start
1545                    const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
1546                    if (thisOne) {
1547                        int iColumn = thisOne->columnNumber();
1548                        bool canDoThisHot = true;
1549                        double targetValue = hotstartSolution[iColumn];
1550                        if (saveUpper[iColumn] > saveLower[iColumn]) {
1551                            double value = saveSolution[iColumn];
1552                            if (hotstartPriorities)
1553                                priorityLevel = hotstartPriorities[iColumn];
1554                            //double originalLower = thisOne->originalLower();
1555                            //double originalUpper = thisOne->originalUpper();
1556                            // switch off if not possible
1557                            if (targetValue >= saveLower[iColumn] && targetValue <= saveUpper[iColumn]) {
1558                                /* priority outranks rest always if negative
1559                                   otherwise can be downgraded if at correct level.
1560                                   Infeasibility may be increased to choose 1.0 values first.
1561                                   choose one near wanted value
1562                                */
1563                                if (fabs(value - targetValue) > integerTolerance) {
1564                                    //if (infeasibility>0.01)
1565                                    //infeasibility = fabs(1.0e6-fabs(value-targetValue));
1566                                    //else
1567                                    infeasibility = fabs(value - targetValue);
1568                                    //if (targetValue==1.0)
1569                                    //infeasibility += 1.0;
1570                                    if (value > targetValue) {
1571                                        preferredWay = -1;
1572                                    } else {
1573                                        preferredWay = 1;
1574                                    }
1575                                    priorityLevel = CoinAbs(priorityLevel);
1576                                } else if (priorityLevel < 0) {
1577                                    priorityLevel = CoinAbs(priorityLevel);
1578                                    if (targetValue == saveLower[iColumn]) {
1579                                        infeasibility = integerTolerance + 1.0e-12;
1580                                        preferredWay = -1;
1581                                    } else if (targetValue == saveUpper[iColumn]) {
1582                                        infeasibility = integerTolerance + 1.0e-12;
1583                                        preferredWay = 1;
1584                                    } else {
1585                                        // can't
1586                                        priorityLevel += 10000000;
1587                                        canDoThisHot = false;
1588                                    }
1589                                } else {
1590                                    priorityLevel += 10000000;
1591                                    canDoThisHot = false;
1592                                }
1593                            } else {
1594                                // switch off if not possible
1595                                canDoThisHot = false;
1596                            }
1597                            if (canDoThisHot)
1598                                canDoOneHot = true;
1599                        } else if (targetValue < saveLower[iColumn] || targetValue > saveUpper[iColumn]) {
1600                        }
1601                    } else {
1602                        priorityLevel += 10000000;
1603                    }
1604                }
1605                if (infeasibility) {
1606                    // Increase estimated degradation to solution
1607                    estimatedDegradation += CoinMin(object->upEstimate(), object->downEstimate());
1608                    numberUnsatisfied_++;
1609                    sumInfeasibilities_ += infeasibility;
1610                    // Better priority? Flush choices.
1611                    if (priorityLevel < bestPriority) {
1612                        int j;
1613                        iSmallest = 0;
1614                        for (j = 0; j < maximumStrong; j++) {
1615                            choice[j].upMovement = 0.0;
1616                            delete choice[j].possibleBranch;
1617                            choice[j].possibleBranch = NULL;
1618                        }
1619                        bestPriority = priorityLevel;
1620                        mostAway = 1.0e-100;
1621                        numberStrong = 0;
1622                    } else if (priorityLevel > bestPriority) {
1623                        continue;
1624                    }
1625                    // Check for suitability based on infeasibility.
1626                    if (infeasibility > mostAway) {
1627                        //add to list
1628                        choice[iSmallest].upMovement = infeasibility;
1629                        delete choice[iSmallest].possibleBranch;
1630                        CbcObject * obj =
1631                            dynamic_cast <CbcObject *>(object) ;
1632                        assert (obj);
1633                        choice[iSmallest].possibleBranch = obj->createCbcBranch(solver, &usefulInfo, preferredWay);
1634                        numberStrong = CoinMax(numberStrong, iSmallest + 1);
1635                        // Save which object it was
1636                        choice[iSmallest].objectNumber = i;
1637                        int j;
1638                        iSmallest = -1;
1639                        mostAway = 1.0e50;
1640                        for (j = 0; j < maximumStrong; j++) {
1641                            if (choice[j].upMovement < mostAway) {
1642                                mostAway = choice[j].upMovement;
1643                                iSmallest = j;
1644                            }
1645                        }
1646                    }
1647                }
1648            }
1649            if (!canDoOneHot && hotstartSolution) {
1650                // switch off as not possible
1651                hotstartSolution = NULL;
1652                model->setHotstartSolution(NULL, NULL);
1653                usefulInfo.hotstartSolution_ = NULL;
1654            }
1655            if (numberUnsatisfied_) {
1656                // some infeasibilities - go to next steps
1657#ifdef CLP_INVESTIGATE
1658                if (hotstartSolution) {
1659                    int k = choice[0].objectNumber;
1660                    OsiObject * object = model->modifiableObject(k);
1661                    const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
1662                    assert (thisOne);
1663                    int iColumn = thisOne->columnNumber();
1664                    double targetValue = hotstartSolution[iColumn];
1665                    double value = saveSolution[iColumn];
1666                    printf("Branch on %d has target %g (value %g) and current bounds %g and %g\n",
1667                           iColumn, targetValue, value, saveLower[iColumn], saveUpper[iColumn]);
1668                }
1669#endif
1670                break;
1671            } else if (!iPass) {
1672                // looks like a solution - get paranoid
1673                bool roundAgain = false;
1674                // get basis
1675                CoinWarmStartBasis * ws = dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
1676                if (!ws)
1677                    break;
1678                for (i = 0; i < numberColumns; i++) {
1679                    double value = saveSolution[i];
1680                    if (value < lower[i]) {
1681                        saveSolution[i] = lower[i];
1682                        roundAgain = true;
1683                        ws->setStructStatus(i, CoinWarmStartBasis::atLowerBound);
1684                    } else if (value > upper[i]) {
1685                        saveSolution[i] = upper[i];
1686                        roundAgain = true;
1687                        ws->setStructStatus(i, CoinWarmStartBasis::atUpperBound);
1688                    }
1689                }
1690                if (roundAgain && saveNumberStrong) {
1691                    // restore basis
1692                    solver->setWarmStart(ws);
1693                    delete ws;
1694                    solver->resolve();
1695                    memcpy(saveSolution, solver->getColSolution(), numberColumns*sizeof(double));
1696                    model->reserveCurrentSolution(saveSolution);
1697                    if (!solver->isProvenOptimal()) {
1698                        // infeasible
1699                        anyAction = -2;
1700                        break;
1701                    }
1702                } else {
1703                    delete ws;
1704                    break;
1705                }
1706            }
1707        }
1708        /* Some solvers can do the strong branching calculations faster if
1709           they do them all at once.  At present only Clp does for ordinary
1710           integers but I think this coding would be easy to modify
1711        */
1712        bool allNormal = true; // to say if we can do fast strong branching
1713        // Say which one will be best
1714        int bestChoice = 0;
1715        double worstInfeasibility = 0.0;
1716        for (i = 0; i < numberStrong; i++) {
1717            choice[i].numIntInfeasUp = numberUnsatisfied_;
1718            choice[i].numIntInfeasDown = numberUnsatisfied_;
1719            choice[i].fix = 0; // say not fixed
1720            if (!dynamic_cast <const CbcSimpleInteger *> (model->object(choice[i].objectNumber)))
1721                allNormal = false; // Something odd so lets skip clever fast branching
1722            if ( !model->object(choice[i].objectNumber)->boundBranch())
1723                numberStrong = 0; // switch off
1724            if ( choice[i].possibleBranch->numberBranches() > 2)
1725                numberStrong = 0; // switch off
1726            // Do best choice in case switched off
1727            if (choice[i].upMovement > worstInfeasibility) {
1728                worstInfeasibility = choice[i].upMovement;
1729                bestChoice = i;
1730            }
1731        }
1732        // If we have hit max time don't do strong branching
1733        bool hitMaxTime = ( CoinCpuTime() - model->getDblParam(CbcModel::CbcStartSeconds) >
1734                            model->getDblParam(CbcModel::CbcMaximumSeconds));
1735        // also give up if we are looping round too much
1736        if (hitMaxTime || numberPassesLeft <= 0)
1737            numberStrong = 0;
1738        /*
1739          Is strong branching enabled? If so, set up and do it. Otherwise, we'll
1740          fall through to simple branching.
1741
1742          Setup for strong branching involves saving the current basis (for restoration
1743          afterwards) and setting up for hot starts.
1744        */
1745        if (numberStrong && saveNumberStrong) {
1746
1747            bool solveAll = false; // set true to say look at all even if some fixed (experiment)
1748            solveAll = true;
1749            // worth trying if too many times
1750            // Save basis
1751            CoinWarmStart * ws = solver->getWarmStart();
1752            // save limit
1753            int saveLimit;
1754            solver->getIntParam(OsiMaxNumIterationHotStart, saveLimit);
1755            if (beforeSolution && saveLimit < 100)
1756                solver->setIntParam(OsiMaxNumIterationHotStart, 100); // go to end
1757#     ifdef COIN_HAS_CLP
1758            /* If we are doing all strong branching in one go then we create new arrays
1759               to store information.  If clp NULL then doing old way.
1760               Going down -
1761               outputSolution[2*i] is final solution.
1762               outputStuff[2*i] is status (0 - finished, 1 infeas, other unknown
1763               outputStuff[2*i+numberStrong] is number iterations
1764               On entry newUpper[i] is new upper bound, on exit obj change
1765               Going up -
1766               outputSolution[2*i+1] is final solution.
1767               outputStuff[2*i+1] is status (0 - finished, 1 infeas, other unknown
1768               outputStuff[2*i+1+numberStrong] is number iterations
1769            On entry newLower[i] is new lower bound, on exit obj change
1770            */
1771            ClpSimplex * clp = NULL;
1772            double * newLower = NULL;
1773            double * newUpper = NULL;
1774            double ** outputSolution = NULL;
1775            int * outputStuff = NULL;
1776            // Go back to normal way if user wants it
1777            if (osiclp && (osiclp->specialOptions()&16) != 0 && osiclp->specialOptions() > 0)
1778                allNormal = false;
1779            if (osiclp && !allNormal) {
1780                // say do fast
1781                int easy = 1;
1782                osiclp->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, &easy) ;
1783            }
1784            if (osiclp && allNormal) {
1785                clp = osiclp->getModelPtr();
1786                // Clp - do a different way
1787                newLower = new double[numberStrong];
1788                newUpper = new double[numberStrong];
1789                outputSolution = new double * [2*numberStrong];
1790                outputStuff = new int [4*numberStrong];
1791                int * which = new int[numberStrong];
1792                int startFinishOptions;
1793                int specialOptions = osiclp->specialOptions();
1794                int clpOptions = clp->specialOptions();
1795                int returnCode = 0;
1796#define CRUNCH
1797#ifdef CRUNCH
1798                // Crunch down problem
1799                int numberRows = clp->numberRows();
1800                // Use dual region
1801                double * rhs = clp->dualRowSolution();
1802                int * whichRow = new int[3*numberRows];
1803                int * whichColumn = new int[2*numberColumns];
1804                int nBound;
1805                ClpSimplex * small = static_cast<ClpSimplexOther *> (clp)->crunch(rhs, whichRow, whichColumn, nBound, true);
1806                if (!small) {
1807                    anyAction = -2;
1808                    //printf("XXXX Inf by inspection\n");
1809                    delete [] whichColumn;
1810                    whichColumn = NULL;
1811                    delete [] whichRow;
1812                    whichRow = NULL;
1813                    break;
1814                } else {
1815                    clp = small;
1816                }
1817#else
1818                int saveLogLevel = clp->logLevel();
1819                int saveMaxIts = clp->maximumIterations();
1820#endif
1821                clp->setLogLevel(0);
1822                if ((specialOptions&1) == 0) {
1823                    startFinishOptions = 0;
1824                    clp->setSpecialOptions(clpOptions | (64 | 1024));
1825                } else {
1826                    startFinishOptions = 1 + 2 + 4;
1827                    //startFinishOptions=1+4; // for moment re-factorize
1828                    if ((specialOptions&4) == 0)
1829                        clp->setSpecialOptions(clpOptions | (64 | 128 | 512 | 1024 | 4096));
1830                    else
1831                        clp->setSpecialOptions(clpOptions | (64 | 128 | 512 | 1024 | 2048 | 4096));
1832                }
1833                // User may want to clean up before strong branching
1834                if ((clp->specialOptions()&32) != 0) {
1835                    clp->primal(1);
1836                    if (clp->numberIterations())
1837                        model->messageHandler()->message(CBC_ITERATE_STRONG, *model->messagesPointer())
1838                        << clp->numberIterations()
1839                        << CoinMessageEol;
1840                }
1841                clp->setMaximumIterations(saveLimit);
1842#ifdef CRUNCH
1843                int * backColumn = whichColumn + numberColumns;
1844#endif
1845                for (i = 0; i < numberStrong; i++) {
1846                    int iObject = choice[i].objectNumber;
1847                    const OsiObject * object = model->object(iObject);
1848                    const CbcSimpleInteger * simple = static_cast <const CbcSimpleInteger *> (object);
1849                    int iSequence = simple->columnNumber();
1850                    newLower[i] = ceil(saveSolution[iSequence]);
1851                    newUpper[i] = floor(saveSolution[iSequence]);
1852#ifdef CRUNCH
1853                    iSequence = backColumn[iSequence];
1854                    assert (iSequence >= 0);
1855#endif
1856                    which[i] = iSequence;
1857                    outputSolution[2*i] = new double [numberColumns];
1858                    outputSolution[2*i+1] = new double [numberColumns];
1859                }
1860                //clp->writeMps("bad");
1861                returnCode = clp->strongBranching(numberStrong, which,
1862                                                  newLower, newUpper, outputSolution,
1863                                                  outputStuff, outputStuff + 2 * numberStrong, !solveAll, false,
1864                                                  startFinishOptions);
1865#ifndef CRUNCH
1866                clp->setSpecialOptions(clpOptions); // restore
1867                clp->setMaximumIterations(saveMaxIts);
1868                clp->setLogLevel(saveLogLevel);
1869#endif
1870                if (returnCode == -2) {
1871                    // bad factorization!!!
1872                    // Doing normal way
1873                    // Mark hot start
1874                    solver->markHotStart();
1875                    clp = NULL;
1876                } else {
1877#ifdef CRUNCH
1878                    // extract solution
1879                    //bool checkSol=true;
1880                    for (i = 0; i < numberStrong; i++) {
1881                        int iObject = choice[i].objectNumber;
1882                        const OsiObject * object = model->object(iObject);
1883                        const CbcSimpleInteger * simple = static_cast <const CbcSimpleInteger *> (object);
1884                        int iSequence = simple->columnNumber();
1885                        which[i] = iSequence;
1886                        double * sol = outputSolution[2*i];
1887                        double * sol2 = outputSolution[2*i+1];
1888                        //bool x=true;
1889                        //bool x2=true;
1890                        for (int iColumn = numberColumns - 1; iColumn >= 0; iColumn--) {
1891                            int jColumn = backColumn[iColumn];
1892                            if (jColumn >= 0) {
1893                                sol[iColumn] = sol[jColumn];
1894                                sol2[iColumn] = sol2[jColumn];
1895                            } else {
1896                                sol[iColumn] = saveSolution[iColumn];
1897                                sol2[iColumn] = saveSolution[iColumn];
1898                            }
1899                        }
1900                    }
1901#endif
1902                }
1903#ifdef CRUNCH
1904                delete [] whichColumn;
1905                delete [] whichRow;
1906                delete small;
1907#endif
1908                delete [] which;
1909            } else {
1910                // Doing normal way
1911                // Mark hot start
1912                solver->markHotStart();
1913            }
1914#     else      /* COIN_HAS_CLP */
1915
1916            OsiSolverInterface *clp = NULL ;
1917            double **outputSolution = NULL ;
1918            int *outputStuff = NULL ;
1919            double * newLower = NULL ;
1920            double * newUpper = NULL ;
1921
1922            solver->markHotStart();
1923
1924#     endif     /* COIN_HAS_CLP */
1925            /*
1926              Open a loop to do the strong branching LPs. For each candidate variable,
1927              solve an LP with the variable forced down, then up. If a direction turns
1928              out to be infeasible or monotonic (i.e., over the dual objective cutoff),
1929              force the objective change to be big (1.0e100). If we determine the problem
1930              is infeasible, or find a monotone variable, escape the loop.
1931
1932              TODO: The `restore bounds' part might be better encapsulated as an
1933            unbranch() method. Branching objects more exotic than simple integers
1934            or cliques might not restrict themselves to variable bounds.
1935
1936              TODO: Virtuous solvers invalidate the current solution (or give bogus
1937            results :-) when the bounds are changed out from under them. So we
1938            need to do all the work associated with finding a new solution before
1939            restoring the bounds.
1940            */
1941            for (i = 0 ; i < numberStrong ; i++) {
1942                double objectiveChange ;
1943                double newObjectiveValue = 1.0e100;
1944                // status is 0 finished, 1 infeasible and other
1945                int iStatus;
1946                /*
1947                  Try the down direction first. (Specify the initial branching alternative as
1948                  down with a call to way(-1). Each subsequent call to branch() performs the
1949                  specified branch and advances the branch object state to the next branch
1950                  alternative.)
1951                */
1952                if (!clp) {
1953                    choice[i].possibleBranch->way(-1) ;
1954                    choice[i].possibleBranch->branch() ;
1955                    bool feasible = true;
1956                    if (checkFeasibility) {
1957                        // check branching did not make infeasible
1958                        int iColumn;
1959                        int numberColumns = solver->getNumCols();
1960                        const double * columnLower = solver->getColLower();
1961                        const double * columnUpper = solver->getColUpper();
1962                        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1963                            if (columnLower[iColumn] > columnUpper[iColumn] + 1.0e-5)
1964                                feasible = false;
1965                        }
1966                    }
1967                    if (feasible) {
1968                        solver->solveFromHotStart() ;
1969                        numberStrongDone++;
1970                        numberStrongIterations += solver->getIterationCount();
1971                        /*
1972                        We now have an estimate of objective degradation that we can use for strong
1973                        branching. If we're over the cutoff, the variable is monotone up.
1974                        If we actually made it to optimality, check for a solution, and if we have
1975                        a good one, call setBestSolution to process it. Note that this may reduce the
1976                        cutoff, so we check again to see if we can declare this variable monotone.
1977                        */
1978                        if (solver->isProvenOptimal())
1979                            iStatus = 0; // optimal
1980                        else if (solver->isIterationLimitReached()
1981                                 && !solver->isDualObjectiveLimitReached())
1982                            iStatus = 2; // unknown
1983                        else
1984                            iStatus = 1; // infeasible
1985                        newObjectiveValue = solver->getObjSense() * solver->getObjValue();
1986                        choice[i].numItersDown = solver->getIterationCount();
1987                    } else {
1988                        iStatus = 1; // infeasible
1989                        newObjectiveValue = 1.0e100;
1990                        choice[i].numItersDown = 0;
1991                    }
1992                } else {
1993                    iStatus = outputStuff[2*i];
1994                    choice[i].numItersDown = outputStuff[2*numberStrong+2*i];
1995                    numberStrongDone++;
1996                    numberStrongIterations += choice[i].numItersDown;
1997                    newObjectiveValue = objectiveValue + newUpper[i];
1998                    solver->setColSolution(outputSolution[2*i]);
1999                }
2000                objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_, 0.0);
2001                if (!iStatus) {
2002                    choice[i].finishedDown = true ;
2003                    if (newObjectiveValue >= model->getCutoff()) {
2004                        objectiveChange = 1.0e100; // say infeasible
2005                        numberStrongInfeasible++;
2006                    } else {
2007                        // See if integer solution
2008                        if (model->feasibleSolution(choice[i].numIntInfeasDown,
2009                                                    choice[i].numObjInfeasDown)
2010                                && model->problemFeasibility()->feasible(model, -1) >= 0) {
2011                            model->setBestSolution(CBC_STRONGSOL,
2012                                                   newObjectiveValue,
2013                                                   solver->getColSolution()) ;
2014                            // only needed for odd solvers
2015                            newObjectiveValue = solver->getObjSense() * solver->getObjValue();
2016                            objectiveChange = CoinMax(newObjectiveValue - objectiveValue_, 0.0) ;
2017                            model->setLastHeuristic(NULL);
2018                            model->incrementUsed(solver->getColSolution());
2019                            if (newObjectiveValue >= model->getCutoff()) {      //  *new* cutoff
2020                                objectiveChange = 1.0e100 ;
2021                                numberStrongInfeasible++;
2022                            }
2023                        }
2024                    }
2025                } else if (iStatus == 1) {
2026                    objectiveChange = 1.0e100 ;
2027                    numberStrongInfeasible++;
2028                } else {
2029                    // Can't say much as we did not finish
2030                    choice[i].finishedDown = false ;
2031                    numberUnfinished++;
2032                }
2033                choice[i].downMovement = objectiveChange ;
2034
2035                // restore bounds
2036                if (!clp) {
2037                    for (int j = 0; j < numberColumns; j++) {
2038                        if (saveLower[j] != lower[j])
2039                            solver->setColLower(j, saveLower[j]);
2040                        if (saveUpper[j] != upper[j])
2041                            solver->setColUpper(j, saveUpper[j]);
2042                    }
2043                }
2044                //printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
2045                //     choice[i].objectNumber,iStatus,newObjectiveValue,choice[i].numItersDown,
2046                //     choice[i].downMovement,choice[i].finishedDown,choice[i].numIntInfeasDown,
2047                //     choice[i].numObjInfeasDown);
2048
2049                // repeat the whole exercise, forcing the variable up
2050                if (!clp) {
2051                    bool feasible = true;
2052                    // If odd branching then maybe just one possibility
2053                    if (choice[i].possibleBranch->numberBranchesLeft() > 0) {
2054                        choice[i].possibleBranch->branch();
2055                        if (checkFeasibility) {
2056                            // check branching did not make infeasible
2057                            int iColumn;
2058                            int numberColumns = solver->getNumCols();
2059                            const double * columnLower = solver->getColLower();
2060                            const double * columnUpper = solver->getColUpper();
2061                            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2062                                if (columnLower[iColumn] > columnUpper[iColumn] + 1.0e-5)
2063                                    feasible = false;
2064                            }
2065                        }
2066                    } else {
2067                        // second branch infeasible
2068                        feasible = false;
2069                    }
2070                    if (feasible) {
2071                        solver->solveFromHotStart() ;
2072                        numberStrongDone++;
2073                        numberStrongIterations += solver->getIterationCount();
2074                        /*
2075                        We now have an estimate of objective degradation that we can use for strong
2076                        branching. If we're over the cutoff, the variable is monotone up.
2077                        If we actually made it to optimality, check for a solution, and if we have
2078                        a good one, call setBestSolution to process it. Note that this may reduce the
2079                        cutoff, so we check again to see if we can declare this variable monotone.
2080                        */
2081                        if (solver->isProvenOptimal())
2082                            iStatus = 0; // optimal
2083                        else if (solver->isIterationLimitReached()
2084                                 && !solver->isDualObjectiveLimitReached())
2085                            iStatus = 2; // unknown
2086                        else
2087                            iStatus = 1; // infeasible
2088                        newObjectiveValue = solver->getObjSense() * solver->getObjValue();
2089                        choice[i].numItersUp = solver->getIterationCount();
2090                    } else {
2091                        iStatus = 1; // infeasible
2092                        newObjectiveValue = 1.0e100;
2093                        choice[i].numItersDown = 0;
2094                    }
2095                } else {
2096                    iStatus = outputStuff[2*i+1];
2097                    choice[i].numItersUp = outputStuff[2*numberStrong+2*i+1];
2098                    numberStrongDone++;
2099                    numberStrongIterations += choice[i].numItersUp;
2100                    newObjectiveValue = objectiveValue + newLower[i];
2101                    solver->setColSolution(outputSolution[2*i+1]);
2102                }
2103                objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_, 0.0);
2104                if (!iStatus) {
2105                    choice[i].finishedUp = true ;
2106                    if (newObjectiveValue >= model->getCutoff()) {
2107                        objectiveChange = 1.0e100; // say infeasible
2108                        numberStrongInfeasible++;
2109                    } else {
2110                        // See if integer solution
2111                        if (model->feasibleSolution(choice[i].numIntInfeasUp,
2112                                                    choice[i].numObjInfeasUp)
2113                                && model->problemFeasibility()->feasible(model, -1) >= 0) {
2114                            model->setBestSolution(CBC_STRONGSOL,
2115                                                   newObjectiveValue,
2116                                                   solver->getColSolution()) ;
2117                            // only needed for odd solvers
2118                            newObjectiveValue = solver->getObjSense() * solver->getObjValue();
2119                            objectiveChange = CoinMax(newObjectiveValue - objectiveValue_, 0.0) ;
2120                            model->setLastHeuristic(NULL);
2121                            model->incrementUsed(solver->getColSolution());
2122                            if (newObjectiveValue >= model->getCutoff()) {      //  *new* cutoff
2123                                objectiveChange = 1.0e100 ;
2124                                numberStrongInfeasible++;
2125                            }
2126                        }
2127                    }
2128                } else if (iStatus == 1) {
2129                    objectiveChange = 1.0e100 ;
2130                    numberStrongInfeasible++;
2131                } else {
2132                    // Can't say much as we did not finish
2133                    choice[i].finishedUp = false ;
2134                    numberUnfinished++;
2135                }
2136                choice[i].upMovement = objectiveChange ;
2137
2138                // restore bounds
2139                if (!clp) {
2140                    for (int j = 0; j < numberColumns; j++) {
2141                        if (saveLower[j] != lower[j])
2142                            solver->setColLower(j, saveLower[j]);
2143                        if (saveUpper[j] != upper[j])
2144                            solver->setColUpper(j, saveUpper[j]);
2145                    }
2146                }
2147
2148                //printf("Up on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
2149                //     choice[i].objectNumber,iStatus,newObjectiveValue,choice[i].numItersUp,
2150                //     choice[i].upMovement,choice[i].finishedUp,choice[i].numIntInfeasUp,
2151                //     choice[i].numObjInfeasUp);
2152
2153                /*
2154                  End of evaluation for this candidate variable. Possibilities are:
2155                  * Both sides below cutoff; this variable is a candidate for branching.
2156                  * Both sides infeasible or above the objective cutoff: no further action
2157                  here. Break from the evaluation loop and assume the node will be purged
2158                  by the caller.
2159                  * One side below cutoff: Install the branch (i.e., fix the variable). Break
2160                  from the evaluation loop and assume the node will be reoptimised by the
2161                  caller.
2162                */
2163                // reset
2164                choice[i].possibleBranch->resetNumberBranchesLeft();
2165                if (choice[i].upMovement < 1.0e100) {
2166                    if (choice[i].downMovement < 1.0e100) {
2167                        // feasible - no action
2168                    } else {
2169                        // up feasible, down infeasible
2170                        anyAction = -1;
2171                        //printf("Down infeasible for choice %d sequence %d\n",i,
2172                        // model->object(choice[i].objectNumber)->columnNumber());
2173                        if (!solveAll) {
2174                            choice[i].possibleBranch->way(1);
2175                            choice[i].possibleBranch->branch();
2176                            break;
2177                        } else {
2178                            choice[i].fix = 1;
2179                        }
2180                    }
2181                } else {
2182                    if (choice[i].downMovement < 1.0e100) {
2183                        // down feasible, up infeasible
2184                        anyAction = -1;
2185                        //printf("Up infeasible for choice %d sequence %d\n",i,
2186                        // model->object(choice[i].objectNumber)->columnNumber());
2187                        if (!solveAll) {
2188                            choice[i].possibleBranch->way(-1);
2189                            choice[i].possibleBranch->branch();
2190                            break;
2191                        } else {
2192                            choice[i].fix = -1;
2193                        }
2194                    } else {
2195                        // neither side feasible
2196                        anyAction = -2;
2197                        //printf("Both infeasible for choice %d sequence %d\n",i,
2198                        // model->object(choice[i].objectNumber)->columnNumber());
2199                        break;
2200                    }
2201                }
2202                bool hitMaxTime = ( CoinCpuTime() - model->getDblParam(CbcModel::CbcStartSeconds) >
2203                                    model->getDblParam(CbcModel::CbcMaximumSeconds));
2204                if (hitMaxTime) {
2205                    numberStrong = i + 1;
2206                    break;
2207                }
2208            }
2209            if (!clp) {
2210                // Delete the snapshot
2211                solver->unmarkHotStart();
2212            } else {
2213                delete [] newLower;
2214                delete [] newUpper;
2215                delete [] outputStuff;
2216                int i;
2217                for (i = 0; i < 2*numberStrong; i++)
2218                    delete [] outputSolution[i];
2219                delete [] outputSolution;
2220            }
2221            solver->setIntParam(OsiMaxNumIterationHotStart, saveLimit);
2222            // restore basis
2223            solver->setWarmStart(ws);
2224            // Unless infeasible we will carry on
2225            // But we could fix anyway
2226            if (anyAction == -1 && solveAll) {
2227                // apply and take off
2228                for (i = 0 ; i < numberStrong ; i++) {
2229                    if (choice[i].fix) {
2230                        choice[i].possibleBranch->way(choice[i].fix) ;
2231                        choice[i].possibleBranch->branch() ;
2232                    }
2233                }
2234                bool feasible = true;
2235                if (checkFeasibility) {
2236                    // check branching did not make infeasible
2237                    int iColumn;
2238                    int numberColumns = solver->getNumCols();
2239                    const double * columnLower = solver->getColLower();
2240                    const double * columnUpper = solver->getColUpper();
2241                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2242                        if (columnLower[iColumn] > columnUpper[iColumn] + 1.0e-5)
2243                            feasible = false;
2244                    }
2245                }
2246                if (feasible) {
2247                    // can do quick optimality check
2248                    int easy = 2;
2249                    solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, &easy) ;
2250                    solver->resolve() ;
2251                    solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
2252                    feasible = solver->isProvenOptimal();
2253                }
2254                if (feasible) {
2255                    memcpy(saveSolution, solver->getColSolution(), numberColumns*sizeof(double));
2256                    model->reserveCurrentSolution(saveSolution);
2257                    memcpy(saveLower, solver->getColLower(), numberColumns*sizeof(double));
2258                    memcpy(saveUpper, solver->getColUpper(), numberColumns*sizeof(double));
2259                    // Clean up all candidates whih are fixed
2260                    int numberLeft = 0;
2261                    for (i = 0 ; i < numberStrong ; i++) {
2262                        CbcStrongInfo thisChoice = choice[i];
2263                        choice[i].possibleBranch = NULL;
2264                        const OsiObject * object = model->object(thisChoice.objectNumber);
2265                        int preferredWay;
2266                        double infeasibility = object->infeasibility(&usefulInfo, preferredWay);
2267                        if (!infeasibility) {
2268                            // take out
2269                            delete thisChoice.possibleBranch;
2270                        } else {
2271                            choice[numberLeft++] = thisChoice;
2272                        }
2273                    }
2274                    numberStrong = numberLeft;
2275                    for (; i < maximumStrong; i++) {
2276                        delete choice[i].possibleBranch;
2277                        choice[i].possibleBranch = NULL;
2278                    }
2279                    // If all fixed then round again
2280                    if (!numberLeft) {
2281                        finished = false;
2282                        numberStrong = 0;
2283                        saveNumberStrong = 0;
2284                        maximumStrong = 1;
2285                    } else {
2286                        anyAction = 0;
2287                    }
2288                    // If these two uncommented then different action
2289                    anyAction = -1;
2290                    finished = true;
2291                    //printf("some fixed but continuing %d left\n",numberLeft);
2292                } else {
2293                    anyAction = -2; // say infeasible
2294                }
2295            }
2296            delete ws;
2297            int numberNodes = model->getNodeCount();
2298            // update number of strong iterations etc
2299            model->incrementStrongInfo(numberStrongDone, numberStrongIterations,
2300                                       anyAction == -2 ? 0 : numberStrongInfeasible, anyAction == -2);
2301
2302            /*
2303              anyAction >= 0 indicates that strong branching didn't produce any monotone
2304              variables. Sift through the candidates for the best one.
2305
2306              QUERY: Setting numberNodes looks to be a distributed noop. numberNodes is
2307              local to this code block. Perhaps should be numberNodes_ from model?
2308              Unclear what this calculation is doing.
2309            */
2310            if (anyAction >= 0) {
2311
2312                // get average cost per iteration and assume stopped ones
2313                // would stop after 50% more iterations at average cost??? !!! ???
2314                double averageCostPerIteration = 0.0;
2315                double totalNumberIterations = 1.0;
2316                int smallestNumberInfeasibilities = COIN_INT_MAX;
2317                for (i = 0; i < numberStrong; i++) {
2318                    totalNumberIterations += choice[i].numItersDown +
2319                                             choice[i].numItersUp ;
2320                    averageCostPerIteration += choice[i].downMovement +
2321                                               choice[i].upMovement;
2322                    smallestNumberInfeasibilities =
2323                        CoinMin(CoinMin(choice[i].numIntInfeasDown ,
2324                                        choice[i].numIntInfeasUp ),
2325                                smallestNumberInfeasibilities);
2326                }
2327                //if (smallestNumberInfeasibilities>=numberIntegerInfeasibilities)
2328                //numberNodes=1000000; // switch off search for better solution
2329                numberNodes = 1000000; // switch off anyway
2330                averageCostPerIteration /= totalNumberIterations;
2331                // all feasible - choose best bet
2332
2333                // New method does all at once so it can be more sophisticated
2334                // in deciding how to balance actions.
2335                // But it does need arrays
2336                double * changeUp = new double [numberStrong];
2337                int * numberInfeasibilitiesUp = new int [numberStrong];
2338                double * changeDown = new double [numberStrong];
2339                int * numberInfeasibilitiesDown = new int [numberStrong];
2340                CbcBranchingObject ** objects = new CbcBranchingObject * [ numberStrong];
2341                for (i = 0 ; i < numberStrong ; i++) {
2342                    int iColumn = choice[i].possibleBranch->variable() ;
2343                    model->messageHandler()->message(CBC_STRONG, *model->messagesPointer())
2344                    << i << iColumn
2345                    << choice[i].downMovement << choice[i].numIntInfeasDown
2346                    << choice[i].upMovement << choice[i].numIntInfeasUp
2347                    << choice[i].possibleBranch->value()
2348                    << CoinMessageEol;
2349                    changeUp[i] = choice[i].upMovement;
2350                    numberInfeasibilitiesUp[i] = choice[i].numIntInfeasUp;
2351                    changeDown[i] = choice[i].downMovement;
2352                    numberInfeasibilitiesDown[i] = choice[i].numIntInfeasDown;
2353                    objects[i] = choice[i].possibleBranch;
2354                }
2355                int whichObject = decision->bestBranch(objects, numberStrong, numberUnsatisfied_,
2356                                                       changeUp, numberInfeasibilitiesUp,
2357                                                       changeDown, numberInfeasibilitiesDown,
2358                                                       objectiveValue_);
2359                // move branching object and make sure it will not be deleted
2360                if (whichObject >= 0) {
2361                    branch_ = objects[whichObject];
2362                    if (model->messageHandler()->logLevel() > 3)
2363                        printf("Choosing column %d\n", choice[whichObject].possibleBranch->variable()) ;
2364                    choice[whichObject].possibleBranch = NULL;
2365                }
2366                delete [] changeUp;
2367                delete [] numberInfeasibilitiesUp;
2368                delete [] changeDown;
2369                delete [] numberInfeasibilitiesDown;
2370                delete [] objects;
2371            }
2372#     ifdef COIN_HAS_CLP
2373            if (osiclp && !allNormal) {
2374                // back to normal
2375                osiclp->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
2376            }
2377#     endif
2378        }
2379        /*
2380          Simple branching. Probably just one, but we may have got here
2381          because of an odd branch e.g. a cut
2382        */
2383        else {
2384            // not strong
2385            // C) create branching object
2386            branch_ = choice[bestChoice].possibleBranch;
2387            choice[bestChoice].possibleBranch = NULL;
2388        }
2389    }
2390    // Set guessed solution value
2391    guessedObjectiveValue_ = objectiveValue_ + estimatedDegradation;
2392    /*
2393      Cleanup, then we're outta here.
2394    */
2395    if (!model->branchingMethod() || dynamicBranchingObject)
2396        delete decision;
2397
2398    for (i = 0; i < maximumStrong; i++)
2399        delete choice[i].possibleBranch;
2400    delete [] choice;
2401    delete [] saveLower;
2402    delete [] saveUpper;
2403
2404    // restore solution
2405    solver->setColSolution(saveSolution);
2406    delete [] saveSolution;
2407# ifdef COIN_HAS_CLP
2408    if (osiclp)
2409        osiclp->setSpecialOptions(saveClpOptions);
2410# endif
2411    return anyAction;
2412}
2413
2414/*
2415  Version for dynamic pseudo costs.
2416
2417  **** For now just return if anything odd
2418  later allow even if odd
2419
2420  The routine scans through the object list of the model looking for objects
2421  that indicate infeasibility. It tests each object using strong branching
2422  and selects the one with the least objective degradation.  A corresponding
2423  branching object is left attached to lastNode.
2424  This version gives preference in evaluation to variables which
2425  have not been evaluated many times.  It also uses numberStrong
2426  to say give up if last few tries have not changed incumbent.
2427  See Achterberg, Koch and Martin.
2428
2429  If strong branching is disabled, a candidate object is chosen essentially
2430  at random (whatever object ends up in pos'n 0 of the candidate array).
2431
2432  If a branching candidate is found to be monotone, bounds are set to fix the
2433  variable and the routine immediately returns (the caller is expected to
2434  reoptimize).
2435
2436  If a branching candidate is found to result in infeasibility in both
2437  directions, the routine immediately returns an indication of infeasibility.
2438
2439  Returns:  0   both branch directions are feasible
2440  -1    branching variable is monotone
2441  -2    infeasible
2442  -3   Use another method
2443
2444  For now just fix on objective from strong branching.
2445*/
2446
2447int CbcNode::chooseDynamicBranch (CbcModel *model, CbcNode *lastNode,
2448                                  OsiSolverBranch * & /*branches*/,
2449                                  int numberPassesLeft)
2450
2451{
2452    if (lastNode)
2453        depth_ = lastNode->depth_ + 1;
2454    else
2455        depth_ = 0;
2456    // Go to other choose if hot start
2457    if (model->hotstartSolution() &&
2458            (((model->moreSpecialOptions()&1024) == 0) || false))
2459        return -3;
2460    delete branch_;
2461    branch_ = NULL;
2462    OsiSolverInterface * solver = model->solver();
2463    // get information on solver type
2464    const OsiAuxInfo * auxInfo = solver->getAuxiliaryInfo();
2465    const OsiBabSolver * auxiliaryInfo = dynamic_cast<const OsiBabSolver *> (auxInfo);
2466    if (!auxiliaryInfo) {
2467        // use one from CbcModel
2468        auxiliaryInfo = model->solverCharacteristics();
2469    }
2470    int numberObjects = model->numberObjects();
2471    // If very odd set of objects then use older chooseBranch
2472    bool useOldWay = false;
2473    // point to useful information
2474    OsiBranchingInformation usefulInfo = model->usefulInformation();
2475    if (numberObjects > model->numberIntegers()) {
2476        for (int i = model->numberIntegers(); i < numberObjects; i++) {
2477            OsiObject * object = model->modifiableObject(i);
2478            CbcObject * obj =   dynamic_cast <CbcObject *>(object) ;
2479            if (!obj || !obj->optionalObject()) {
2480                int preferredWay;
2481                double infeasibility = object->infeasibility(&usefulInfo, preferredWay);
2482                if (infeasibility) {
2483                    useOldWay = true;
2484                    break;
2485                }
2486            }
2487        }
2488    }
2489    if ((model->specialOptions()&128) != 0)
2490        useOldWay = false; // allow
2491    // For now return if not simple
2492    if (useOldWay)
2493        return -3;
2494    // Modify useful info
2495    usefulInfo.depth_ = depth_;
2496    if ((model->specialOptions()&128) != 0) {
2497        // SOS - shadow prices
2498        int numberRows = solver->getNumRows();
2499        const double * pi = usefulInfo.pi_;
2500        double sumPi = 0.0;
2501        for (int i = 0; i < numberRows; i++)
2502            sumPi += fabs(pi[i]);
2503        sumPi /= static_cast<double> (numberRows);
2504        // and scale back
2505        sumPi *= 0.01;
2506        usefulInfo.defaultDual_ = sumPi; // switch on
2507        int numberColumns = solver->getNumCols();
2508        int size = CoinMax(numberColumns, 2 * numberRows);
2509        usefulInfo.usefulRegion_ = new double [size];
2510        CoinZeroN(usefulInfo.usefulRegion_, size);
2511        usefulInfo.indexRegion_ = new int [size];
2512        // pi may change
2513        usefulInfo.pi_ = CoinCopyOfArray(usefulInfo.pi_, numberRows);
2514    }
2515    assert (auxiliaryInfo);
2516    double cutoff = model->getCutoff();
2517    const double * lower = solver->getColLower();
2518    const double * upper = solver->getColUpper();
2519    // See if user thinks infeasible
2520    int anyAction = model->problemFeasibility()->feasible(model, 0);
2521    if (anyAction) {
2522        // will return -2 if infeasible , 0 if treat as integer
2523        return anyAction - 1;
2524    }
2525    int i;
2526    int saveStateOfSearch = model->stateOfSearch() % 10;
2527    int numberStrong = model->numberStrong();
2528    /* Ranging is switched off.
2529       The idea is that you can find out the effect of one iteration
2530       on each unsatisfied variable cheaply.  Then use this
2531       if you have not got much else to go on.
2532    */
2533    //#define RANGING
2534#ifdef RANGING
2535    // must have clp
2536#ifndef COIN_HAS_CLP
2537#  warning("Ranging switched off as not Clp");
2538#undef RANGING
2539#endif
2540    // Pass number
2541    int kPass = 0;
2542    int numberRows = solver->getNumRows();
2543#endif
2544    int numberColumns = model->getNumCols();
2545    double * saveUpper = new double[numberColumns];
2546    double * saveLower = new double[numberColumns];
2547    for (i = 0; i < numberColumns; i++) {
2548        saveLower[i] = lower[i];
2549        saveUpper[i] = upper[i];
2550    }
2551
2552    // Save solution in case heuristics need good solution later
2553
2554    double * saveSolution = new double[numberColumns];
2555    memcpy(saveSolution, solver->getColSolution(), numberColumns*sizeof(double));
2556    model->reserveCurrentSolution(saveSolution);
2557    const double * hotstartSolution = model->hotstartSolution();
2558    const int * hotstartPriorities = model->hotstartPriorities();
2559    double integerTolerance =
2560        model->getDblParam(CbcModel::CbcIntegerTolerance);
2561    if (hotstartSolution) {
2562        if ((model->moreSpecialOptions()&1024) != 0) {
2563            int nBad = 0;
2564            int nUnsat = 0;
2565            int nDiff = 0;
2566            for (int i = 0; i < numberObjects; i++) {
2567                OsiObject * object = model->modifiableObject(i);
2568                const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
2569                if (thisOne) {
2570                    int iColumn = thisOne->columnNumber();
2571                    double targetValue = hotstartSolution[iColumn];
2572                    double value = saveSolution[iColumn];
2573                    if (fabs(value - floor(value + 0.5)) > 1.0e-6) {
2574                        nUnsat++;
2575#ifdef CLP_INVESTIGATE
2576                        printf("H %d is %g target %g\n", iColumn, value, targetValue);
2577#endif
2578                    } else if (fabs(targetValue - value) > 1.0e-6) {
2579                        nDiff++;
2580                    }
2581                    if (targetValue < saveLower[iColumn] ||
2582                            targetValue > saveUpper[iColumn]) {
2583#ifdef CLP_INVESTIGATE
2584                        printf("%d has target %g and current bounds %g and %g\n",
2585                               iColumn, targetValue, saveLower[iColumn], saveUpper[iColumn]);
2586#endif
2587                        nBad++;
2588                    }
2589                }
2590            }
2591#ifdef CLP_INVESTIGATE
2592            printf("Hot %d unsatisfied, %d outside limits, %d different\n",
2593                   nUnsat, nBad, nDiff);
2594#endif
2595            if (nBad) {
2596                // switch off as not possible
2597                hotstartSolution = NULL;
2598                model->setHotstartSolution(NULL, NULL);
2599                usefulInfo.hotstartSolution_ = NULL;
2600            }
2601        }
2602    }
2603    /*
2604      Get a branching decision object. Use the default dynamic decision criteria unless
2605      the user has loaded a decision method into the model.
2606    */
2607    CbcBranchDecision *decision = model->branchingMethod();
2608    if (!decision)
2609        decision = new CbcBranchDynamicDecision();
2610    int xMark = 0;
2611    // Get arrays to sort
2612    double * sort = new double[numberObjects];
2613    int * whichObject = new int[numberObjects];
2614#ifdef RANGING
2615    int xPen = 0;
2616    int * objectMark = new int[2*numberObjects+1];
2617#endif
2618    // Arrays with movements
2619    double * upEstimate = new double[numberObjects];
2620    double * downEstimate = new double[numberObjects];
2621    double estimatedDegradation = 0.0;
2622    int numberNodes = model->getNodeCount();
2623    int saveLogLevel = model->logLevel();
2624#if 0
2625    if ((numberNodes % 500) == 0) {
2626        model->setLogLevel(6);
2627        // Get average up and down costs
2628        double averageUp = 0.0;
2629        double averageDown = 0.0;
2630        int numberUp = 0;
2631        int numberDown = 0;
2632        int i;
2633        for ( i = 0; i < numberObjects; i++) {
2634            OsiObject * object = model->modifiableObject(i);
2635            CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2636                dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2637            assert(dynamicObject);
2638            int  numberUp2 = 0;
2639            int numberDown2 = 0;
2640            double up = 0.0;
2641            double down = 0.0;
2642            if (dynamicObject->numberTimesUp()) {
2643                numberUp++;
2644                averageUp += dynamicObject->upDynamicPseudoCost();
2645                numberUp2 += dynamicObject->numberTimesUp();
2646                up = dynamicObject->upDynamicPseudoCost();
2647            }
2648            if (dynamicObject->numberTimesDown()) {
2649                numberDown++;
2650                averageDown += dynamicObject->downDynamicPseudoCost();
2651                numberDown2 += dynamicObject->numberTimesDown();
2652                down = dynamicObject->downDynamicPseudoCost();
2653            }
2654            if (numberUp2 || numberDown2)
2655                printf("col %d - up %d times cost %g, - down %d times cost %g\n",
2656                       dynamicObject->columnNumber(), numberUp2, up, numberDown2, down);
2657        }
2658        if (numberUp)
2659            averageUp /= static_cast<double> (numberUp);
2660        else
2661            averageUp = 1.0;
2662        if (numberDown)
2663            averageDown /= static_cast<double> (numberDown);
2664        else
2665            averageDown = 1.0;
2666        printf("total - up %d vars average %g, - down %d vars average %g\n",
2667               numberUp, averageUp, numberDown, averageDown);
2668    }
2669#endif
2670    int numberBeforeTrust = model->numberBeforeTrust();
2671    // May go round twice if strong branching fixes all local candidates
2672    bool finished = false;
2673    int numberToFix = 0;
2674# ifdef COIN_HAS_CLP
2675    OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
2676    int saveClpOptions = 0;
2677    if (osiclp) {
2678        // for faster hot start
2679        saveClpOptions = osiclp->specialOptions();
2680        osiclp->setSpecialOptions(saveClpOptions | 8192);
2681    }
2682# else
2683    OsiSolverInterface *osiclp = NULL ;
2684# endif
2685    //const CglTreeProbingInfo * probingInfo = NULL; //model->probingInfo();
2686    // Old code left in with DEPRECATED_STRATEGY
2687    assert (model->searchStrategy() == -1 ||
2688            model->searchStrategy() == 1 ||
2689            model->searchStrategy() == 2);
2690#ifdef DEPRECATED_STRATEGY
2691    int saveSearchStrategy2 = model->searchStrategy();
2692#endif
2693    // Get average up and down costs
2694    {
2695        double averageUp = 0.0;
2696        double averageDown = 0.0;
2697        int numberUp = 0;
2698        int numberDown = 0;
2699        int i;
2700        for ( i = 0; i < numberObjects; i++) {
2701            OsiObject * object = model->modifiableObject(i);
2702            CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2703                dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2704            if (dynamicObject) {
2705                if (dynamicObject->numberTimesUp()) {
2706                    numberUp++;
2707                    averageUp += dynamicObject->upDynamicPseudoCost();
2708                }
2709                if (dynamicObject->numberTimesDown()) {
2710                    numberDown++;
2711                    averageDown += dynamicObject->downDynamicPseudoCost();
2712                }
2713            }
2714        }
2715        if (numberUp)
2716            averageUp /= static_cast<double> (numberUp);
2717        else
2718            averageUp = 1.0;
2719        if (numberDown)
2720            averageDown /= static_cast<double> (numberDown);
2721        else
2722            averageDown = 1.0;
2723        for ( i = 0; i < numberObjects; i++) {
2724            OsiObject * object = model->modifiableObject(i);
2725            CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2726                dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2727            if (dynamicObject) {
2728                if (!dynamicObject->numberTimesUp())
2729                    dynamicObject->setUpDynamicPseudoCost(averageUp);
2730                if (!dynamicObject->numberTimesDown())
2731                    dynamicObject->setDownDynamicPseudoCost(averageDown);
2732            }
2733        }
2734    }
2735    /*
2736      1 strong
2737      2 no strong
2738      3 strong just before solution
2739      4 no strong just before solution
2740      5 strong first time or before solution
2741      6 strong first time
2742    */
2743    int useShadow = model->moreSpecialOptions() & 7;
2744    if (useShadow > 2) {
2745        if (model->getSolutionCount()) {
2746            if (numberNodes || useShadow < 5) {
2747                useShadow = 0;
2748                // zap pseudo shadow prices
2749                model->pseudoShadow(-1);
2750                // and switch off
2751                model->setMoreSpecialOptions(model->moreSpecialOptions()&(~1023));
2752            } else {
2753                useShadow = 1;
2754            }
2755        } else if (useShadow < 5) {
2756            useShadow -= 2;
2757        } else {
2758            useShadow = 1;
2759        }
2760    }
2761    if (useShadow) {
2762        // pseudo shadow prices
2763        model->pseudoShadow((model->moreSpecialOptions() >> 3)&63);
2764    }
2765#ifdef DEPRECATED_STRATEGY
2766    { // in for tabbing
2767    } else if (saveSearchStrategy2 < 1999) {
2768        // pseudo shadow prices
2769        model->pseudoShadow(NULL, NULL);
2770    } else if (saveSearchStrategy2 < 2999) {
2771        // leave old ones
2772    } else if (saveSearchStrategy2 < 3999) {
2773        // pseudo shadow prices at root
2774        if (!numberNodes)
2775            model->pseudoShadow(NULL, NULL);
2776    } else {
2777        abort();
2778    }
2779    if (saveSearchStrategy2 >= 0)
2780        saveSearchStrategy2 = saveSearchStrategy2 % 1000;
2781    if (saveSearchStrategy2 == 999)
2782        saveSearchStrategy2 = -1;
2783    int saveSearchStrategy = saveSearchStrategy2 < 99 ? saveSearchStrategy2 : saveSearchStrategy2 - 100;
2784#endif //DEPRECATED_STRATEGY
2785    int numberNotTrusted = 0;
2786    int numberStrongDone = 0;
2787    int numberUnfinished = 0;
2788    int numberStrongInfeasible = 0;
2789    int numberStrongIterations = 0;
2790    // so we can save lots of stuff
2791    CbcStrongInfo choice;
2792    CbcDynamicPseudoCostBranchingObject * choiceObject = NULL;
2793    if (model->allDynamic()) {
2794        CbcSimpleIntegerDynamicPseudoCost * object = NULL;
2795        choiceObject = new CbcDynamicPseudoCostBranchingObject(model, 0, -1, 0.5, object);
2796    }
2797    choice.possibleBranch = choiceObject;
2798    numberPassesLeft = CoinMax(numberPassesLeft, 2);
2799    while (!finished) {
2800        numberPassesLeft--;
2801        finished = true;
2802        decision->initialize(model);
2803        // Some objects may compute an estimate of best solution from here
2804        estimatedDegradation = 0.0;
2805        numberToFix = 0;
2806        int numberToDo = 0;
2807        int iBestNot = -1;
2808        int iBestGot = -1;
2809        double best = 0.0;
2810        numberNotTrusted = 0;
2811        numberStrongDone = 0;
2812        numberUnfinished = 0;
2813        numberStrongInfeasible = 0;
2814        numberStrongIterations = 0;
2815#ifdef RANGING
2816        int * which = objectMark + numberObjects + 1;
2817        int neededPenalties;
2818        int optionalPenalties;
2819#endif
2820        // We may go round this loop three times (only if we think we have solution)
2821        for (int iPass = 0; iPass < 3; iPass++) {
2822
2823            // Some objects may compute an estimate of best solution from here
2824            estimatedDegradation = 0.0;
2825            numberUnsatisfied_ = 0;
2826            // initialize sum of "infeasibilities"
2827            sumInfeasibilities_ = 0.0;
2828            int bestPriority = COIN_INT_MAX;
2829#if 0
2830            int number01 = 0;
2831            const cliqueEntry * entry = NULL;
2832            const int * toZero = NULL;
2833            const int * toOne = NULL;
2834            const int * backward = NULL;
2835            int numberUnsatisProbed = 0;
2836            int numberUnsatisNotProbed = 0; // 0-1
2837            if (probingInfo) {
2838                number01 = probingInfo->numberIntegers();
2839                entry = probingInfo->fixEntries();
2840                toZero = probingInfo->toZero();
2841                toOne = probingInfo->toOne();
2842                backward = probingInfo->backward();
2843                if (!toZero[number01] || number01 < numberObjects || true) {
2844                    // no info
2845                    probingInfo = NULL;
2846                }
2847            }
2848#endif
2849            /*
2850              Scan for branching objects that indicate infeasibility. Choose candidates
2851              using priority as the first criteria, then integer infeasibility.
2852
2853              The algorithm is to fill the array with a set of good candidates (by
2854              infeasibility) with priority bestPriority.  Finding a candidate with
2855              priority better (less) than bestPriority flushes the choice array. (This
2856              serves as initialization when the first candidate is found.)
2857
2858            */
2859            numberToDo = 0;
2860#ifdef RANGING
2861            neededPenalties = 0;
2862            optionalPenalties = numberObjects;
2863#endif
2864            iBestNot = -1;
2865            double bestNot = 0.0;
2866            iBestGot = -1;
2867            best = 0.0;
2868            /* Problem type as set by user or found by analysis.  This will be extended
2869            0 - not known
2870            1 - Set partitioning <=
2871            2 - Set partitioning ==
2872            3 - Set covering
2873            4 - all +- 1 or all +1 and odd
2874            */
2875            int problemType = model->problemType();
2876            bool canDoOneHot = false;
2877            for (i = 0; i < numberObjects; i++) {
2878                OsiObject * object = model->modifiableObject(i);
2879                CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2880                    dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2881                int preferredWay;
2882                double infeasibility = object->infeasibility(&usefulInfo, preferredWay);
2883                int priorityLevel = object->priority();
2884                if (hotstartSolution) {
2885                    // we are doing hot start
2886                    const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
2887                    if (thisOne) {
2888                        int iColumn = thisOne->columnNumber();
2889                        bool canDoThisHot = true;
2890                        double targetValue = hotstartSolution[iColumn];
2891                        if (saveUpper[iColumn] > saveLower[iColumn]) {
2892                            double value = saveSolution[iColumn];
2893                            if (hotstartPriorities)
2894                                priorityLevel = hotstartPriorities[iColumn];
2895                            //double originalLower = thisOne->originalLower();
2896                            //double originalUpper = thisOne->originalUpper();
2897                            // switch off if not possible
2898                            if (targetValue >= saveLower[iColumn] && targetValue <= saveUpper[iColumn]) {
2899                                /* priority outranks rest always if negative
2900                                   otherwise can be downgraded if at correct level.
2901                                   Infeasibility may be increased to choose 1.0 values first.
2902                                   choose one near wanted value
2903                                */
2904                                if (fabs(value - targetValue) > integerTolerance) {
2905                                    //if (infeasibility>0.01)
2906                                    //infeasibility = fabs(1.0e6-fabs(value-targetValue));
2907                                    //else
2908                                    infeasibility = fabs(value - targetValue);
2909                                    //if (targetValue==1.0)
2910                                    //infeasibility += 1.0;
2911                                    if (value > targetValue) {
2912                                        preferredWay = -1;
2913                                    } else {
2914                                        preferredWay = 1;
2915                                    }
2916                                    priorityLevel = CoinAbs(priorityLevel);
2917                                } else if (priorityLevel < 0) {
2918                                    priorityLevel = CoinAbs(priorityLevel);
2919                                    if (targetValue == saveLower[iColumn]) {
2920                                        infeasibility = integerTolerance + 1.0e-12;
2921                                        preferredWay = -1;
2922                                    } else if (targetValue == saveUpper[iColumn]) {
2923                                        infeasibility = integerTolerance + 1.0e-12;
2924                                        preferredWay = 1;
2925                                    } else {
2926                                        // can't
2927                                        priorityLevel += 10000000;
2928                                        canDoThisHot = false;
2929                                    }
2930                                } else {
2931                                    priorityLevel += 10000000;
2932                                    canDoThisHot = false;
2933                                }
2934                            } else {
2935                                // switch off if not possible
2936                                canDoThisHot = false;
2937                            }
2938                            if (canDoThisHot)
2939                                canDoOneHot = true;
2940                        } else if (targetValue < saveLower[iColumn] || targetValue > saveUpper[iColumn]) {
2941                        }
2942                    } else {
2943                        priorityLevel += 10000000;
2944                    }
2945                }
2946#define ZERO_ONE 0
2947#define ZERO_FAKE 1.0e20;
2948#if ZERO_ONE==1
2949                // branch on 0-1 first (temp)
2950                if (fabs(saveSolution[dynamicObject->columnNumber()]) < 1.0)
2951                    priorityLevel--;
2952#endif
2953#if ZERO_ONE==2
2954                if (fabs(saveSolution[dynamicObject->columnNumber()]) < 1.0)
2955                    infeasibility *= ZERO_FAKE;
2956#endif
2957                if (infeasibility) {
2958                    int iColumn = numberColumns + i;
2959                    bool gotDown = false;
2960                    int numberThisDown = 0;
2961                    bool gotUp = false;
2962                    int numberThisUp = 0;
2963                    double downGuess = object->downEstimate();
2964                    double upGuess = object->upEstimate();
2965                    if (dynamicObject) {
2966                        // Use this object's numberBeforeTrust
2967                        int numberBeforeTrust = dynamicObject->numberBeforeTrust();
2968                        iColumn = dynamicObject->columnNumber();
2969                        gotDown = false;
2970                        numberThisDown = dynamicObject->numberTimesDown();
2971                        if (numberThisDown >= numberBeforeTrust)
2972                            gotDown = true;
2973                        gotUp = false;
2974                        numberThisUp = dynamicObject->numberTimesUp();
2975                        if (numberThisUp >= numberBeforeTrust)
2976                            gotUp = true;
2977                        if (!depth_ && false) {
2978                            // try closest to 0.5
2979                            double part = saveSolution[iColumn] - floor(saveSolution[iColumn]);
2980                            infeasibility = fabs(0.5 - part);
2981                        }
2982                        if (problemType > 0 && problemType < 4 && false) {
2983                            // try closest to 0.5
2984                            double part = saveSolution[iColumn] - floor(saveSolution[iColumn]);
2985                            infeasibility = 0.5 - fabs(0.5 - part);
2986                        }
2987#if 0
2988                        if (probingInfo) {
2989                            int iSeq = backward[iColumn];
2990                            assert (iSeq >= 0);
2991                            infeasibility = 1.0 + (toZero[iSeq+1] - toZero[iSeq]) +
2992                                            5.0 * CoinMin(toOne[iSeq] - toZero[iSeq], toZero[iSeq+1] - toOne[iSeq]);
2993                            if (toZero[iSeq+1] > toZero[iSeq]) {
2994                                numberUnsatisProbed++;
2995                            } else {
2996                                numberUnsatisNotProbed++;
2997                            }
2998                        }
2999#endif
3000                    } else {
3001                        // see if SOS
3002                        CbcSOS * sosObject =
3003                            dynamic_cast <CbcSOS *>(object) ;
3004                        if (sosObject) {
3005                            gotDown = false;
3006                            numberThisDown = sosObject->numberTimesDown();
3007                            if (numberThisDown >= numberBeforeTrust)
3008                                gotDown = true;
3009                            gotUp = false;
3010                            numberThisUp = sosObject->numberTimesUp();
3011                            if (numberThisUp >= numberBeforeTrust)
3012                                gotUp = true;
3013                        } else {
3014                            gotDown = true;
3015                            numberThisDown = 999999;
3016                            downGuess = 1.0e20;
3017                            gotUp = true;
3018                            numberThisUp = 999999;
3019                            upGuess = 1.0e20;
3020                            numberPassesLeft = 0;
3021                        }
3022                    }
3023                    // Increase estimated degradation to solution
3024                    estimatedDegradation += CoinMin(downGuess, upGuess);
3025                    downEstimate[i] = downGuess;
3026                    upEstimate[i] = upGuess;
3027                    numberUnsatisfied_++;
3028                    sumInfeasibilities_ += infeasibility;
3029                    // Better priority? Flush choices.
3030                    if (priorityLevel < bestPriority) {
3031                        numberToDo = 0;
3032                        bestPriority = priorityLevel;
3033                        iBestGot = -1;
3034                        best = 0.0;
3035                        numberNotTrusted = 0;
3036#ifdef RANGING
3037                        neededPenalties = 0;
3038                        optionalPenalties = numberObjects;
3039#endif
3040                    } else if (priorityLevel > bestPriority) {
3041                        continue;
3042                    }
3043                    if (!gotUp || !gotDown)
3044                        numberNotTrusted++;
3045                    // Check for suitability based on infeasibility.
3046                    if ((gotDown && gotUp) && numberStrong > 0) {
3047                        sort[numberToDo] = -infeasibility;
3048                        if (infeasibility > best) {
3049                            best = infeasibility;
3050                            iBestGot = numberToDo;
3051                        }
3052#ifdef RANGING
3053                        if (dynamicObject) {
3054                            objectMark[--optionalPenalties] = numberToDo;
3055                            which[optionalPenalties] = iColumn;
3056                        }
3057#endif
3058                    } else {
3059#ifdef RANGING
3060                        if (dynamicObject) {
3061                            objectMark[neededPenalties] = numberToDo;
3062                            which[neededPenalties++] = iColumn;
3063                        }
3064#endif
3065                        sort[numberToDo] = -10.0 * infeasibility;
3066                        if (!(numberThisUp + numberThisDown))
3067                            sort[numberToDo] *= 100.0; // make even more likely
3068                        if (iColumn < numberColumns) {
3069                            double part = saveSolution[iColumn] - floor(saveSolution[iColumn]);
3070                            if (1.0 - fabs(part - 0.5) > bestNot) {
3071                                iBestNot = numberToDo;
3072                                bestNot = 1.0 - fabs(part - 0.5);
3073                            }
3074                        } else {
3075                            // SOS
3076                            if (-sort[numberToDo] > bestNot) {
3077                                iBestNot = numberToDo;
3078                                bestNot = -sort[numberToDo];
3079                            }
3080                        }
3081                    }
3082                    if (model->messageHandler()->logLevel() > 3) {
3083                        printf("%d (%d) down %d %g up %d %g - infeas %g - sort %g solution %g\n",
3084                               i, iColumn, numberThisDown, object->downEstimate(), numberThisUp, object->upEstimate(),
3085                               infeasibility, sort[numberToDo], saveSolution[iColumn]);
3086                    }
3087                    whichObject[numberToDo++] = i;
3088                } else {
3089                    // for debug
3090                    downEstimate[i] = -1.0;
3091                    upEstimate[i] = -1.0;
3092                }
3093            }
3094            if (numberUnsatisfied_) {
3095                //if (probingInfo&&false)
3096                //printf("nunsat %d, %d probed, %d other 0-1\n",numberUnsatisfied_,
3097                // numberUnsatisProbed,numberUnsatisNotProbed);
3098                // some infeasibilities - go to next steps
3099                if (!canDoOneHot && hotstartSolution) {
3100                    // switch off as not possible
3101                    hotstartSolution = NULL;
3102                    model->setHotstartSolution(NULL, NULL);
3103                    usefulInfo.hotstartSolution_ = NULL;
3104                }
3105                break;
3106            } else if (!iPass) {
3107                // may just need resolve
3108                model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3109                //double newObjValue = solver->getObjSense()*solver->getObjValue();
3110                //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3111                if (!solver->isProvenOptimal()) {
3112                    // infeasible
3113                    anyAction = -2;
3114                    break;
3115                }
3116            } else if (iPass == 1) {
3117                // looks like a solution - get paranoid
3118                bool roundAgain = false;
3119                // get basis
3120                CoinWarmStartBasis * ws = dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
3121                if (!ws)
3122                    break;
3123                double tolerance;
3124                solver->getDblParam(OsiPrimalTolerance, tolerance);
3125                for (i = 0; i < numberColumns; i++) {
3126                    double value = saveSolution[i];
3127                    if (value < lower[i] - tolerance) {
3128                        saveSolution[i] = lower[i];
3129                        roundAgain = true;
3130                        ws->setStructStatus(i, CoinWarmStartBasis::atLowerBound);
3131                    } else if (value > upper[i] + tolerance) {
3132                        saveSolution[i] = upper[i];
3133                        roundAgain = true;
3134                        ws->setStructStatus(i, CoinWarmStartBasis::atUpperBound);
3135                    }
3136                }
3137                if (roundAgain) {
3138                    // restore basis
3139                    solver->setWarmStart(ws);
3140                    solver->setColSolution(saveSolution);
3141                    delete ws;
3142                    bool takeHint;
3143                    OsiHintStrength strength;
3144                    solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
3145                    solver->setHintParam(OsiDoDualInResolve, false, OsiHintDo) ;
3146                    model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3147                    //double newObjValue = solver->getObjSense()*solver->getObjValue();
3148                    //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3149                    solver->setHintParam(OsiDoDualInResolve, takeHint, strength) ;
3150                    if (!solver->isProvenOptimal()) {
3151                        // infeasible
3152                        anyAction = -2;
3153                        break;
3154                    }
3155                } else {
3156                    delete ws;
3157                    break;
3158                }
3159            }
3160        }
3161        if (anyAction == -2) {
3162            break;
3163        }
3164        // skip if solution
3165        if (!numberUnsatisfied_)
3166            break;
3167        int skipAll = (numberNotTrusted == 0 || numberToDo == 1) ? 1 : 0;
3168        bool doneHotStart = false;
3169        //DEPRECATED_STRATEGYint searchStrategy = saveSearchStrategy>=0 ? (saveSearchStrategy%10) : -1;
3170        int searchStrategy = model->searchStrategy();
3171        // But adjust depending on ratio of iterations
3172        if (searchStrategy > 0) {
3173            if (numberBeforeTrust >= 5 && numberBeforeTrust <= 10) {
3174                if (searchStrategy != 2) {
3175                    assert (searchStrategy == 1);
3176                    if (depth_ > 5) {
3177                        int numberIterations = model->getIterationCount();
3178                        int numberStrongIterations = model->numberStrongIterations();
3179                        if (numberStrongIterations > numberIterations + 10000) {
3180                            searchStrategy = 2;
3181                            skipAll = 1;
3182                        } else if (numberStrongIterations*4 + 1000 < numberIterations) {
3183                            searchStrategy = 3;
3184                            skipAll = 0;
3185                        }
3186                    } else {
3187                        searchStrategy = 3;
3188                        skipAll = 0;
3189                    }
3190                }
3191            }
3192        }
3193        // worth trying if too many times
3194        // Save basis
3195        CoinWarmStart * ws = NULL;
3196        // save limit
3197        int saveLimit = 0;
3198        solver->getIntParam(OsiMaxNumIterationHotStart, saveLimit);
3199        if (!numberPassesLeft)
3200            skipAll = 1;
3201        if (!skipAll) {
3202            ws = solver->getWarmStart();
3203            int limit = 100;
3204            if (!saveStateOfSearch && saveLimit < limit && saveLimit == 100)
3205                solver->setIntParam(OsiMaxNumIterationHotStart, limit);
3206        }
3207        // Say which one will be best
3208        int whichChoice = 0;
3209        int bestChoice;
3210        if (iBestGot >= 0)
3211            bestChoice = iBestGot;
3212        else
3213            bestChoice = iBestNot;
3214        assert (bestChoice >= 0);
3215        // If we have hit max time don't do strong branching
3216        bool hitMaxTime = ( CoinCpuTime() - model->getDblParam(CbcModel::CbcStartSeconds) >
3217                            model->getDblParam(CbcModel::CbcMaximumSeconds));
3218        // also give up if we are looping round too much
3219        if (hitMaxTime || numberPassesLeft <= 0 || useShadow == 2) {
3220            int iObject = whichObject[bestChoice];
3221            OsiObject * object = model->modifiableObject(iObject);
3222            int preferredWay;
3223            object->infeasibility(&usefulInfo, preferredWay);
3224            CbcObject * obj =
3225                dynamic_cast <CbcObject *>(object) ;
3226            assert (obj);
3227            branch_ = obj->createCbcBranch(solver, &usefulInfo, preferredWay);
3228            {
3229                CbcBranchingObject * branchObj =
3230                    dynamic_cast <CbcBranchingObject *>(branch_) ;
3231                assert (branchObj);
3232                branchObj->way(preferredWay);
3233            }
3234            delete ws;
3235            ws = NULL;
3236            break;
3237        } else {
3238            // say do fast
3239            int easy = 1;
3240            if (!skipAll)
3241                solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, &easy) ;
3242            int iDo;
3243#ifdef RANGING
3244            bool useRanging = model->allDynamic() && !skipAll;
3245            if (useRanging) {
3246                double currentObjective = solver->getObjValue() * solver->getObjSense();
3247                double gap = cutoff - currentObjective;
3248                // relax a bit
3249                gap *= 1.0000001;
3250                gap = CoinMax(1.0e-5, gap);
3251                // off penalties if too much
3252                double needed = neededPenalties;
3253                needed *= numberRows;
3254                if (numberNodes) {
3255                    if (needed > 1.0e6) {
3256                        neededPenalties = 0;
3257                    } else if (gap < 1.0e5) {
3258                        // maybe allow some not needed
3259                        int extra = static_cast<int> ((1.0e6 - needed) / numberRows);
3260                        int nStored = numberObjects - optionalPenalties;
3261                        extra = CoinMin(extra, nStored);
3262                        for (int i = 0; i < extra; i++) {
3263                            objectMark[neededPenalties] = objectMark[optionalPenalties+i];
3264                            which[neededPenalties++] = which[optionalPenalties+i];;
3265                        }
3266                    }
3267                }
3268                if (osiclp && neededPenalties) {
3269                    assert (!doneHotStart);
3270                    xPen += neededPenalties;
3271                    which--;
3272                    which[0] = neededPenalties;
3273                    osiclp->passInRanges(which);
3274                    // Mark hot start and get ranges
3275                    if (kPass) {
3276                        // until can work out why solution can go funny
3277                        int save = osiclp->specialOptions();
3278                        osiclp->setSpecialOptions(save | 256);
3279                        solver->markHotStart();
3280                        osiclp->setSpecialOptions(save);
3281                    } else {
3282                        solver->markHotStart();
3283                    }
3284                    doneHotStart = true;
3285                    xMark++;
3286                    kPass++;
3287                    osiclp->passInRanges(NULL);
3288                    const double * downCost = osiclp->upRange();
3289                    const double * upCost = osiclp->downRange();
3290                    bool problemFeasible = true;
3291                    int numberFixed = 0;
3292                    for (int i = 0; i < neededPenalties; i++) {
3293                        int j = objectMark[i];
3294                        int iObject = whichObject[j];
3295                        OsiObject * object = model->modifiableObject(iObject);
3296                        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3297                            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3298                        // Use this object's numberBeforeTrust
3299                        int numberBeforeTrust = dynamicObject->numberBeforeTrust();
3300                        int iSequence = dynamicObject->columnNumber();
3301                        double value = saveSolution[iSequence];
3302                        value -= floor(value);
3303                        double upPenalty = CoinMin(upCost[i], 1.0e110) * (1.0 - value);
3304                        double downPenalty = CoinMin(downCost[i], 1.0e110) * value;
3305                        int numberThisDown = dynamicObject->numberTimesDown();
3306                        int numberThisUp = dynamicObject->numberTimesUp();
3307                        if (!numberBeforeTrust) {
3308                            // override
3309                            downEstimate[iObject] = downPenalty;
3310                            upEstimate[iObject] = upPenalty;
3311                            double min1 = CoinMin(downEstimate[iObject],
3312                                                  upEstimate[iObject]);
3313                            double max1 = CoinMax(downEstimate[iObject],
3314                                                  upEstimate[iObject]);
3315                            min1 = 0.8 * min1 + 0.2 * max1;
3316                            sort[j] = - min1;
3317                        } else if (numberThisDown < numberBeforeTrust ||
3318                                   numberThisUp < numberBeforeTrust) {
3319                            double invTrust = 1.0 / static_cast<double> (numberBeforeTrust);
3320                            if (numberThisDown < numberBeforeTrust) {
3321                                double fraction = numberThisDown * invTrust;
3322                                downEstimate[iObject] = fraction * downEstimate[iObject] + (1.0 - fraction) * downPenalty;
3323                            }
3324                            if (numberThisUp < numberBeforeTrust) {
3325                                double fraction = numberThisUp * invTrust;
3326                                upEstimate[iObject] = fraction * upEstimate[iObject] + (1.0 - fraction) * upPenalty;
3327                            }
3328                            double min1 = CoinMin(downEstimate[iObject],
3329                                                  upEstimate[iObject]);
3330                            double max1 = CoinMax(downEstimate[iObject],
3331                                                  upEstimate[iObject]);
3332                            min1 = 0.8 * min1 + 0.2 * max1;
3333                            min1 *= 10.0;
3334                            if (!(numberThisDown + numberThisUp))
3335                                min1 *= 100.0;
3336                            sort[j] = - min1;
3337                        }
3338                        if (CoinMax(downPenalty, upPenalty) > gap) {
3339                            printf("gap %g object %d has down range %g, up %g\n",
3340                                   gap, i, downPenalty, upPenalty);
3341                            //sort[j] -= 1.0e50; // make more likely to be chosen
3342                            int number;
3343                            if (downPenalty > gap) {
3344                                number = dynamicObject->numberTimesDown();
3345                                if (upPenalty > gap)
3346                                    problemFeasible = false;
3347                                CbcBranchingObject * branch = dynamicObject->createCbcBranch(solver, &usefulInfo, 1);
3348                                //branch->fix(solver,saveLower,saveUpper,1);
3349                                delete branch;
3350                            } else {
3351                                number = dynamicObject->numberTimesUp();
3352                                CbcBranchingObject * branch = dynamicObject->createCbcBranch(solver, &usefulInfo, 1);
3353                                //branch->fix(solver,saveLower,saveUpper,-1);
3354                                delete branch;
3355                            }
3356                            if (number >= numberBeforeTrust)
3357                                dynamicObject->setNumberBeforeTrust(number + 1);
3358                            numberFixed++;
3359                        }
3360                        if (!numberNodes)
3361                            printf("%d pen down ps %g -> %g up ps %g -> %g\n",
3362                                   iObject, downPenalty, downPenalty, upPenalty, upPenalty);
3363                    }
3364                    if (numberFixed && problemFeasible) {
3365                        assert(doneHotStart);
3366                        solver->unmarkHotStart();
3367                        model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3368                        //double newObjValue = solver->getObjSense()*solver->getObjValue();
3369                        //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3370                        solver->markHotStart();
3371                        problemFeasible = solver->isProvenOptimal();
3372                    }
3373                    if (!problemFeasible) {
3374                        fprintf(stdout, "both ways infeas on ranging - code needed\n");
3375                        anyAction = -2;
3376                        if (!choiceObject) {
3377                            delete choice.possibleBranch;
3378                            choice.possibleBranch = NULL;
3379                        }
3380                        //printf("Both infeasible for choice %d sequence %d\n",i,
3381                        // model->object(choice.objectNumber)->columnNumber());
3382                        // Delete the snapshot
3383                        solver->unmarkHotStart();
3384                        // back to normal
3385                        solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
3386                        // restore basis
3387                        solver->setWarmStart(ws);
3388                        doneHotStart = false;
3389                        delete ws;
3390                        ws = NULL;
3391                        break;
3392                    }
3393                }
3394            }
3395#endif          /* RANGING */
3396            {
3397                int numberIterations = model->getIterationCount();
3398                //numberIterations += (model->numberExtraIterations()>>2);
3399                const int * strongInfo = model->strongInfo();
3400                //int numberDone = strongInfo[0]-strongInfo[3];
3401                int numberFixed = strongInfo[1] - strongInfo[4];
3402                int numberInfeasible = strongInfo[2] - strongInfo[5];
3403                assert (!strongInfo[3]);
3404                assert (!strongInfo[4]);
3405                assert (!strongInfo[5]);
3406                int numberStrongIterations = model->numberStrongIterations();
3407                int numberRows = solver->getNumRows();
3408                if (numberStrongIterations > numberIterations + CoinMin(100, 10*numberRows) && depth_ >= 4 && numberNodes > 100) {
3409                    if (20*numberInfeasible + 4*numberFixed < numberNodes) {
3410                        // Say never do
3411                        skipAll = -1;
3412                    }
3413                }
3414            }
3415            // make sure best will be first
3416            if (iBestGot >= 0)
3417                sort[iBestGot] = -COIN_DBL_MAX;
3418            // Actions 0 - exit for repeat, 1 resolve and try old choice,2 exit for continue
3419            if (anyAction)
3420                numberToDo = 0; // skip as we will be trying again
3421            // Sort
3422            CoinSort_2(sort, sort + numberToDo, whichObject);
3423            // Change in objective opposite infeasible
3424            double worstFeasible = 0.0;
3425            // Just first if strong off
3426            if (!numberStrong)
3427                numberToDo = CoinMin(numberToDo, 1);
3428            if (searchStrategy == 2)
3429                numberToDo = CoinMin(numberToDo, 20);
3430            iDo = 0;
3431            int saveLimit2;
3432            solver->getIntParam(OsiMaxNumIterationHotStart, saveLimit2);
3433            bool doQuickly = false; // numberToDo>2*numberStrong;
3434            if (searchStrategy == 2)
3435                doQuickly = true;
3436            int numberTest = numberNotTrusted > 0 ? numberStrong : (numberStrong + 1) / 2;
3437            if (searchStrategy == 3) {
3438                // Previously decided we need strong
3439                doQuickly = false;
3440                numberTest = numberStrong;
3441            }
3442            // Try nearly always off
3443            if (skipAll >= 0) {
3444                if (searchStrategy < 2) {
3445                    //if ((numberNodes%20)!=0) {
3446                    if ((model->specialOptions()&8) == 0) {
3447                        numberTest = 0;
3448                        doQuickly = true;
3449                    }
3450                    //} else {
3451                    //doQuickly=false;
3452                    //numberTest=2*numberStrong;
3453                    //skipAll=0;
3454                    //}
3455                }
3456            } else {
3457                // Just take first
3458                doQuickly = true;
3459                numberTest = 1;
3460            }
3461            int testDepth = (skipAll >= 0) ? 8 : 4;
3462            if (depth_ < testDepth && numberStrong) {
3463                if (searchStrategy != 2) {
3464                    doQuickly = false;
3465                    int numberRows = solver->getNumRows();
3466                    // whether to do this or not is important - think
3467                    if (numberRows < 300 || numberRows + numberColumns < 2500) {
3468                        if (depth_ < 7)
3469                            numberStrong = CoinMin(3 * numberStrong, numberToDo);
3470                        if (!depth_)
3471                            numberStrong = CoinMin(6 * numberStrong, numberToDo);
3472                    }
3473                    numberTest = numberStrong;
3474                    skipAll = 0;
3475                }
3476            }
3477            // Do at least 5 strong
3478            if (numberColumns < 1000 && (depth_ < 15 || numberNodes < 1000000))
3479                numberTest = CoinMax(numberTest, 5);
3480            if ((model->specialOptions()&8) == 0) {
3481                if (skipAll) {
3482                    numberTest = 0;
3483                    doQuickly = true;
3484                }
3485            } else {
3486                // do 5 as strong is fixing
3487                numberTest = CoinMax(numberTest, 5);
3488            }
3489            // see if switched off
3490            if (skipAll < 0) {
3491                numberTest = 0;
3492                doQuickly = true;
3493            }
3494            int realMaxHotIterations = 999999;
3495            if (skipAll < 0)
3496                numberToDo = 1;
3497#ifdef DO_ALL_AT_ROOT
3498            if (!numberNodes) {
3499                printf("DOX %d test %d unsat %d - nobj %d\n",
3500                       numberToDo, numberTest, numberUnsatisfied_,
3501                       numberObjects);
3502                numberTest = numberToDo;
3503            }
3504#endif
3505            for ( iDo = 0; iDo < numberToDo; iDo++) {
3506                int iObject = whichObject[iDo];
3507                OsiObject * object = model->modifiableObject(iObject);
3508                CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3509                    dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3510                int iColumn = dynamicObject ? dynamicObject->columnNumber() : numberColumns + iObject;
3511                int preferredWay;
3512                double infeasibility = object->infeasibility(&usefulInfo, preferredWay);
3513                // may have become feasible
3514                if (!infeasibility)
3515                    continue;
3516#ifndef NDEBUG
3517                if (iColumn < numberColumns) {
3518                    const double * solution = model->testSolution();
3519                    assert (saveSolution[iColumn] == solution[iColumn]);
3520                }
3521#endif
3522                CbcSimpleInteger * obj =
3523                    dynamic_cast <CbcSimpleInteger *>(object) ;
3524                if (obj) {
3525                    if (choiceObject) {
3526                        obj->fillCreateBranch(choiceObject, &usefulInfo, preferredWay);
3527                        choiceObject->setObject(dynamicObject);
3528                    } else {
3529                        choice.possibleBranch = obj->createCbcBranch(solver, &usefulInfo, preferredWay);
3530                    }
3531                } else {
3532                    CbcObject * obj =
3533                        dynamic_cast <CbcObject *>(object) ;
3534                    assert (obj);
3535                    choice.possibleBranch = obj->createCbcBranch(solver, &usefulInfo, preferredWay);
3536                }
3537                // Save which object it was
3538                choice.objectNumber = iObject;
3539                choice.numIntInfeasUp = numberUnsatisfied_;
3540                choice.numIntInfeasDown = numberUnsatisfied_;
3541                choice.upMovement = upEstimate[iObject];
3542                choice.downMovement = downEstimate[iObject];
3543                assert (choice.upMovement >= 0.0);
3544                assert (choice.downMovement >= 0.0);
3545                choice.fix = 0; // say not fixed
3546                // see if can skip strong branching
3547                int canSkip = choice.possibleBranch->fillStrongInfo(choice);
3548                if ((numberTest <= 0 || skipAll)) {
3549                    if (iDo > 20) {
3550#ifdef DO_ALL_AT_ROOT
3551                        if (!numberNodes)
3552                            printf("DOY test %d - iDo %d\n",
3553                                   numberTest, iDo);
3554#endif
3555                        if (!choiceObject) {
3556                            delete choice.possibleBranch;
3557                            choice.possibleBranch = NULL;
3558                        }
3559                        break; // give up anyway
3560                    }
3561                }
3562                if (model->messageHandler()->logLevel() > 3 && numberBeforeTrust && dynamicObject)
3563                    dynamicObject->print(1, choice.possibleBranch->value());
3564                if (skipAll < 0)
3565                    canSkip = true;
3566                if (!canSkip) {
3567                    if (!doneHotStart) {
3568                        // Mark hot start
3569                        doneHotStart = true;
3570                        solver->markHotStart();
3571                        xMark++;
3572                    }
3573                    numberTest--;
3574                    // just do a few
3575                    if (searchStrategy == 2)
3576                        solver->setIntParam(OsiMaxNumIterationHotStart, 10);
3577                    double objectiveChange ;
3578                    double newObjectiveValue = 1.0e100;
3579                    int j;
3580                    // status is 0 finished, 1 infeasible and other
3581                    int iStatus;
3582                    /*
3583                      Try the down direction first. (Specify the initial branching alternative as
3584                      down with a call to way(-1). Each subsequent call to branch() performs the
3585                      specified branch and advances the branch object state to the next branch
3586                      alternative.)
3587                    */
3588                    choice.possibleBranch->way(-1) ;
3589                    choice.possibleBranch->branch() ;
3590                    solver->solveFromHotStart() ;
3591                    bool needHotStartUpdate = false;
3592                    numberStrongDone++;
3593                    numberStrongIterations += solver->getIterationCount();
3594                    /*
3595                      We now have an estimate of objective degradation that we can use for strong
3596                      branching. If we're over the cutoff, the variable is monotone up.
3597                      If we actually made it to optimality, check for a solution, and if we have
3598                      a good one, call setBestSolution to process it. Note that this may reduce the
3599                      cutoff, so we check again to see if we can declare this variable monotone.
3600                    */
3601                    if (solver->isProvenOptimal())
3602                        iStatus = 0; // optimal
3603                    else if (solver->isIterationLimitReached()
3604                             && !solver->isDualObjectiveLimitReached())
3605                        iStatus = 2; // unknown
3606                    else
3607                        iStatus = 1; // infeasible
3608                    if (iStatus != 2 && solver->getIterationCount() >
3609                            realMaxHotIterations)
3610                        numberUnfinished++;
3611                    newObjectiveValue = solver->getObjSense() * solver->getObjValue();
3612                    choice.numItersDown = solver->getIterationCount();
3613                    objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_, 0.0);
3614                    // Update branching information if wanted
3615                    CbcBranchingObject * cbcobj = dynamic_cast<CbcBranchingObject *> (choice.possibleBranch);
3616                    if (cbcobj) {
3617                        CbcObject * object = cbcobj->object();
3618                        assert (object) ;
3619                        CbcObjectUpdateData update = object->createUpdateInformation(solver, this, cbcobj);
3620                        update.objectNumber_ = choice.objectNumber;
3621                        model->addUpdateInformation(update);
3622                    } else {
3623                        decision->updateInformation( solver, this);
3624                    }
3625                    if (!iStatus) {
3626                        choice.finishedDown = true ;
3627                        if (newObjectiveValue >= cutoff) {
3628                            objectiveChange = 1.0e100; // say infeasible
3629                            numberStrongInfeasible++;
3630                        } else {
3631                            //#define TIGHTEN_BOUNDS
3632#ifdef TIGHTEN_BOUNDS
3633                            // Can we tighten bounds?
3634                            if (iColumn < numberColumns && cutoff < 1.0e20
3635                                    && objectiveChange > 1.0e-5) {
3636                                double value = saveSolution[iColumn];
3637                                double down = value - floor(value);
3638                                double changePer = objectiveChange / (down + 1.0e-7);
3639                                double distance = (cutoff - objectiveValue_) / changePer;
3640                                distance += 1.0e-3;
3641                                if (distance < 5.0) {
3642                                    double newLower = ceil(value - distance);
3643                                    if (newLower > saveLower[iColumn]) {
3644                                        //printf("Could increase lower bound on %d from %g to %g\n",
3645                                        //   iColumn,saveLower[iColumn],newLower);
3646                                        saveLower[iColumn] = newLower;
3647                                        solver->setColLower(iColumn, newLower);
3648                                    }
3649                                }
3650                            }
3651#endif
3652                            // See if integer solution
3653                            if (model->feasibleSolution(choice.numIntInfeasDown,
3654                                                        choice.numObjInfeasDown)
3655                                    && model->problemFeasibility()->feasible(model, -1) >= 0) {
3656                                if (auxiliaryInfo->solutionAddsCuts()) {
3657                                    needHotStartUpdate = true;
3658                                    solver->unmarkHotStart();
3659                                }
3660                                model->setBestSolution(CBC_STRONGSOL,
3661                                                       newObjectiveValue,
3662                                                       solver->getColSolution()) ;
3663                                if (needHotStartUpdate) {
3664                                    model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3665                                    newObjectiveValue = solver->getObjSense() * solver->getObjValue();
3666                                    //objectiveValue_ = CoinMax(objectiveValue_,newObjectiveValue);
3667                                    objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_, 0.0);
3668                                    model->feasibleSolution(choice.numIntInfeasDown,
3669                                                            choice.numObjInfeasDown);
3670                                }
3671                                model->setLastHeuristic(NULL);
3672                                model->incrementUsed(solver->getColSolution());
3673                                cutoff = model->getCutoff();
3674                                if (newObjectiveValue >= cutoff) {      //  *new* cutoff
3675                                    objectiveChange = 1.0e100 ;
3676                                    numberStrongInfeasible++;
3677                                }
3678                            }
3679                        }
3680                    } else if (iStatus == 1) {
3681                        objectiveChange = 1.0e100 ;
3682                        numberStrongInfeasible++;
3683                    } else {
3684                        // Can't say much as we did not finish
3685                        choice.finishedDown = false ;
3686                        numberUnfinished++;
3687                    }
3688                    choice.downMovement = objectiveChange ;
3689
3690                    // restore bounds
3691                    for ( j = 0; j < numberColumns; j++) {
3692                        if (saveLower[j] != lower[j])
3693                            solver->setColLower(j, saveLower[j]);
3694                        if (saveUpper[j] != upper[j])
3695                            solver->setColUpper(j, saveUpper[j]);
3696                    }
3697                    if (needHotStartUpdate) {
3698                        needHotStartUpdate = false;
3699                        model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3700                        //double newObjValue = solver->getObjSense()*solver->getObjValue();
3701                        //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3702                        //we may again have an integer feasible solution
3703                        int numberIntegerInfeasibilities;
3704                        int numberObjectInfeasibilities;
3705                        if (model->feasibleSolution(
3706                                    numberIntegerInfeasibilities,
3707                                    numberObjectInfeasibilities)) {
3708#ifdef BONMIN
3709                            //In this case node has become integer feasible, let us exit the loop
3710                            std::cout << "Node has become integer feasible" << std::endl;
3711                            numberUnsatisfied_ = 0;
3712                            break;
3713#endif
3714                            double objValue = solver->getObjValue();
3715                            model->setBestSolution(CBC_STRONGSOL,
3716                                                   objValue,
3717                                                   solver->getColSolution()) ;
3718                            model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3719                            //double newObjValue = solver->getObjSense()*solver->getObjValue();
3720                            //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3721                            cutoff = model->getCutoff();
3722                        }
3723                        solver->markHotStart();
3724                    }
3725#ifdef DO_ALL_AT_ROOT
3726                    if (!numberNodes)
3727                        printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
3728                               choice.objectNumber, iStatus, newObjectiveValue, choice.numItersDown,
3729                               choice.downMovement, choice.finishedDown, choice.numIntInfeasDown,
3730                               choice.numObjInfeasDown);
3731#endif
3732
3733                    // repeat the whole exercise, forcing the variable up
3734                    choice.possibleBranch->branch();
3735                    solver->solveFromHotStart() ;
3736                    numberStrongDone++;
3737                    numberStrongIterations += solver->getIterationCount();
3738                    /*
3739                      We now have an estimate of objective degradation that we can use for strong
3740                      branching. If we're over the cutoff, the variable is monotone up.
3741                      If we actually made it to optimality, check for a solution, and if we have
3742                      a good one, call setBestSolution to process it. Note that this may reduce the
3743                      cutoff, so we check again to see if we can declare this variable monotone.
3744                    */
3745                    if (solver->isProvenOptimal())
3746                        iStatus = 0; // optimal
3747                    else if (solver->isIterationLimitReached()
3748                             && !solver->isDualObjectiveLimitReached())
3749                        iStatus = 2; // unknown
3750                    else
3751                        iStatus = 1; // infeasible
3752                    if (iStatus != 2 && solver->getIterationCount() >
3753                            realMaxHotIterations)
3754                        numberUnfinished++;
3755                    newObjectiveValue = solver->getObjSense() * solver->getObjValue();
3756                    choice.numItersUp = solver->getIterationCount();
3757                    objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_, 0.0);
3758                    // Update branching information if wanted
3759                    cbcobj = dynamic_cast<CbcBranchingObject *> (choice.possibleBranch);
3760                    if (cbcobj) {
3761                        CbcObject * object = cbcobj->object();
3762                        assert (object) ;
3763                        CbcObjectUpdateData update = object->createUpdateInformation(solver, this, cbcobj);
3764                        update.objectNumber_ = choice.objectNumber;
3765                        model->addUpdateInformation(update);
3766                    } else {
3767                        decision->updateInformation( solver, this);
3768                    }
3769                    if (!iStatus) {
3770                        choice.finishedUp = true ;
3771                        if (newObjectiveValue >= cutoff) {
3772                            objectiveChange = 1.0e100; // say infeasible
3773                            numberStrongInfeasible++;
3774                        } else {
3775#ifdef TIGHTEN_BOUNDS
3776                            // Can we tighten bounds?
3777                            if (iColumn < numberColumns && cutoff < 1.0e20
3778                                    && objectiveChange > 1.0e-5) {
3779                                double value = saveSolution[iColumn];
3780                                double up = ceil(value) - value;
3781                                double changePer = objectiveChange / (up + 1.0e-7);
3782                                double distance = (cutoff - objectiveValue_) / changePer;
3783                                distance += 1.0e-3;
3784                                if (distance < 5.0) {
3785                                    double newUpper = floor(value + distance);
3786                                    if (newUpper < saveUpper[iColumn]) {
3787                                        //printf("Could decrease upper bound on %d from %g to %g\n",
3788                                        //   iColumn,saveUpper[iColumn],newUpper);
3789                                        saveUpper[iColumn] = newUpper;
3790                                        solver->setColUpper(iColumn, newUpper);
3791                                    }
3792                                }
3793                            }
3794#endif
3795                            // See if integer solution
3796                            if (model->feasibleSolution(choice.numIntInfeasUp,
3797                                                        choice.numObjInfeasUp)
3798                                    && model->problemFeasibility()->feasible(model, -1) >= 0) {
3799#ifdef BONMIN
3800                                std::cout << "Node has become integer feasible" << std::endl;
3801                                numberUnsatisfied_ = 0;
3802                                break;
3803#endif
3804                                if (auxiliaryInfo->solutionAddsCuts()) {
3805                                    needHotStartUpdate = true;
3806                                    solver->unmarkHotStart();
3807                                }
3808                                model->setBestSolution(CBC_STRONGSOL,
3809                                                       newObjectiveValue,
3810                                                       solver->getColSolution()) ;
3811                                if (needHotStartUpdate) {
3812                                    model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3813                                    newObjectiveValue = solver->getObjSense() * solver->getObjValue();
3814                                    //objectiveValue_ = CoinMax(objectiveValue_,newObjectiveValue);
3815                                    objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_, 0.0);
3816                                    model->feasibleSolution(choice.numIntInfeasDown,
3817                                                            choice.numObjInfeasDown);
3818                                }
3819                                model->setLastHeuristic(NULL);
3820                                model->incrementUsed(solver->getColSolution());
3821                                cutoff = model->getCutoff();
3822                                if (newObjectiveValue >= cutoff) {      //  *new* cutoff
3823                                    objectiveChange = 1.0e100 ;
3824                                    numberStrongInfeasible++;
3825                                }
3826                            }
3827                        }
3828                    } else if (iStatus == 1) {
3829                        objectiveChange = 1.0e100 ;
3830                        numberStrongInfeasible++;
3831                    } else {
3832                        // Can't say much as we did not finish
3833                        choice.finishedUp = false ;
3834                        numberUnfinished++;
3835                    }
3836                    choice.upMovement = objectiveChange ;
3837
3838                    // restore bounds
3839                    for ( j = 0; j < numberColumns; j++) {
3840                        if (saveLower[j] != lower[j])
3841                            solver->setColLower(j, saveLower[j]);
3842                        if (saveUpper[j] != upper[j])
3843                            solver->setColUpper(j, saveUpper[j]);
3844                    }
3845                    if (needHotStartUpdate) {
3846                        needHotStartUpdate = false;
3847                        model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3848                        //double newObjValue = solver->getObjSense()*solver->getObjValue();
3849                        //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3850                        //we may again have an integer feasible solution
3851                        int numberIntegerInfeasibilities;
3852                        int numberObjectInfeasibilities;
3853                        if (model->feasibleSolution(
3854                                    numberIntegerInfeasibilities,
3855                                    numberObjectInfeasibilities)) {
3856                            double objValue = solver->getObjValue();
3857                            model->setBestSolution(CBC_STRONGSOL,
3858                                                   objValue,
3859                                                   solver->getColSolution()) ;
3860                            model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3861                            //double newObjValue = solver->getObjSense()*solver->getObjValue();
3862                            //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3863                            cutoff = model->getCutoff();
3864                        }
3865                        solver->markHotStart();
3866                    }
3867
3868#ifdef DO_ALL_AT_ROOT
3869                    if (!numberNodes)
3870                        printf("Up on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
3871                               choice.objectNumber, iStatus, newObjectiveValue, choice.numItersUp,
3872                               choice.upMovement, choice.finishedUp, choice.numIntInfeasUp,
3873                               choice.numObjInfeasUp);
3874#endif
3875                }
3876
3877                solver->setIntParam(OsiMaxNumIterationHotStart, saveLimit2);
3878                /*
3879                  End of evaluation for this candidate variable. Possibilities are:
3880                  * Both sides below cutoff; this variable is a candidate for branching.
3881                  * Both sides infeasible or above the objective cutoff: no further action
3882                  here. Break from the evaluation loop and assume the node will be purged
3883                  by the caller.
3884                  * One side below cutoff: Install the branch (i.e., fix the variable). Break
3885                  from the evaluation loop and assume the node will be reoptimised by the
3886                  caller.
3887                */
3888                // reset
3889                choice.possibleBranch->resetNumberBranchesLeft();
3890                if (choice.upMovement < 1.0e100) {
3891                    if (choice.downMovement < 1.0e100) {
3892                        // In case solution coming in was odd
3893                        choice.upMovement = CoinMax(0.0, choice.upMovement);
3894                        choice.downMovement = CoinMax(0.0, choice.downMovement);
3895#if ZERO_ONE==2
3896                        // branch on 0-1 first (temp)
3897                        if (fabs(choice.possibleBranch->value()) < 1.0) {
3898                            choice.upMovement *= ZERO_FAKE;
3899                            choice.downMovement *= ZERO_FAKE;
3900                        }
3901#endif
3902                        // feasible - see which best
3903                        if (!canSkip) {
3904                            if (model->messageHandler()->logLevel() > 3)
3905                                printf("sort %g downest %g upest %g ", sort[iDo], downEstimate[iObject],
3906                                       upEstimate[iObject]);
3907                            model->messageHandler()->message(CBC_STRONG, *model->messagesPointer())
3908                            << iObject << iColumn
3909                            << choice.downMovement << choice.numIntInfeasDown
3910                            << choice.upMovement << choice.numIntInfeasUp
3911                            << choice.possibleBranch->value()
3912                            << CoinMessageEol;
3913                        }
3914                        int betterWay;
3915                        {
3916                            CbcBranchingObject * branchObj =
3917                                dynamic_cast <CbcBranchingObject *>(branch_) ;
3918                            if (branch_)
3919                                assert (branchObj);
3920                            betterWay = decision->betterBranch(choice.possibleBranch,
3921                                                               branchObj,
3922                                                               choice.upMovement,
3923                                                               choice.numIntInfeasUp ,
3924                                                               choice.downMovement,
3925                                                               choice.numIntInfeasDown );
3926                        }
3927                        if (betterWay) {
3928                            // C) create branching object
3929                            if (choiceObject) {
3930                                delete branch_;
3931                                branch_ = choice.possibleBranch->clone();
3932                            } else {
3933                                delete branch_;
3934                                branch_ = choice.possibleBranch;
3935                                choice.possibleBranch = NULL;
3936                            }
3937                            {
3938                                CbcBranchingObject * branchObj =
3939                                    dynamic_cast <CbcBranchingObject *>(branch_) ;
3940                                assert (branchObj);
3941                                //branchObj->way(preferredWay);
3942                                branchObj->way(betterWay);
3943                            }
3944                            bestChoice = choice.objectNumber;
3945                            whichChoice = iDo;
3946                            if (numberStrong <= 1) {
3947                                delete ws;
3948                                ws = NULL;
3949                                break;
3950                            }
3951                        } else {
3952                            if (!choiceObject) {
3953                                delete choice.possibleBranch;
3954                                choice.possibleBranch = NULL;
3955                            }
3956                            if (iDo >= 2*numberStrong) {
3957                                delete ws;
3958                                ws = NULL;
3959                                break;
3960                            }
3961                            if (!dynamicObject || dynamicObject->numberTimesUp() > 1) {
3962                                if (iDo - whichChoice >= numberStrong) {
3963                                    if (!choiceObject) {
3964                                        delete choice.possibleBranch;
3965                                        choice.possibleBranch = NULL;
3966                                    }
3967                                    break; // give up
3968                                }
3969                            } else {
3970                                if (iDo - whichChoice >= 2*numberStrong) {
3971                                    delete ws;
3972                                    ws = NULL;
3973                                    if (!choiceObject) {
3974                                        delete choice.possibleBranch;
3975                                        choice.possibleBranch = NULL;
3976                                    }
3977                                    break; // give up
3978                                }
3979                            }
3980                        }
3981                    } else {
3982                        // up feasible, down infeasible
3983                        anyAction = -1;
3984                        worstFeasible = CoinMax(worstFeasible, choice.upMovement);
3985                        model->messageHandler()->message(CBC_STRONG, *model->messagesPointer())
3986                        << iObject << iColumn
3987                        << choice.downMovement << choice.numIntInfeasDown
3988                        << choice.upMovement << choice.numIntInfeasUp
3989                        << choice.possibleBranch->value()
3990                        << CoinMessageEol;
3991                        //printf("Down infeasible for choice %d sequence %d\n",i,
3992                        // model->object(choice.objectNumber)->columnNumber());
3993                        choice.fix = 1;
3994                        numberToFix++;
3995                        choice.possibleBranch->fix(solver, saveLower, saveUpper, 1);
3996                        if (!choiceObject) {
3997                            choice.possibleBranch = NULL;
3998                        } else {
3999                            //choiceObject = new CbcDynamicPseudoCostBranchingObject(*choiceObject);
4000                            choice.possibleBranch = choiceObject;
4001                        }
4002                        assert(doneHotStart);
4003                        solver->unmarkHotStart();
4004                        model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
4005                        //double newObjValue = solver->getObjSense()*solver->getObjValue();
4006                        //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
4007                        solver->markHotStart();
4008                        // may be infeasible (if other way stopped on iterations)
4009                        if (!solver->isProvenOptimal()) {
4010                            // neither side feasible
4011                            anyAction = -2;
4012                            if (!choiceObject) {
4013                                delete choice.possibleBranch;
4014                                choice.possibleBranch = NULL;
4015                            }
4016                            //printf("Both infeasible for choice %d sequence %d\n",i,
4017                            // model->object(choice.objectNumber)->columnNumber());
4018                            delete ws;
4019                            ws = NULL;
4020                            break;
4021                        }
4022                    }
4023                } else {
4024                    if (choice.downMovement < 1.0e100) {
4025                        // down feasible, up infeasible
4026                        anyAction = -1;
4027                        worstFeasible = CoinMax(worstFeasible, choice.downMovement);
4028                        model->messageHandler()->message(CBC_STRONG, *model->messagesPointer())
4029                        << iObject << iColumn
4030                        << choice.downMovement << choice.numIntInfeasDown
4031                        << choice.upMovement << choice.numIntInfeasUp
4032                        << choice.possibleBranch->value()
4033                        << CoinMessageEol;
4034                        choice.fix = -1;
4035                        numberToFix++;
4036                        choice.possibleBranch->fix(solver, saveLower, saveUpper, -1);
4037                        if (!choiceObject) {
4038                            choice.possibleBranch = NULL;
4039                        } else {
4040                            //choiceObject = new CbcDynamicPseudoCostBranchingObject(*choiceObject);
4041                            choice.possibleBranch = choiceObject;
4042                        }
4043                        assert(doneHotStart);
4044                        solver->unmarkHotStart();
4045                        model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
4046                        //double newObjValue = solver->getObjSense()*solver->getObjValue();
4047                        //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
4048                        solver->markHotStart();
4049                        // may be infeasible (if other way stopped on iterations)
4050                        if (!solver->isProvenOptimal()) {
4051                            // neither side feasible
4052                            anyAction = -2;
4053                            if (!choiceObject) {
4054                                delete choice.possibleBranch;
4055                                choice.possibleBranch = NULL;
4056                            }
4057                            delete ws;
4058                            ws = NULL;
4059                            break;
4060                        }
4061                    } else {
4062                        // neither side feasible
4063                        anyAction = -2;
4064                        if (!choiceObject) {
4065                            delete choice.possibleBranch;
4066                            choice.possibleBranch = NULL;
4067                        }
4068                        delete ws;
4069                        ws = NULL;
4070                        break;
4071                    }
4072                }
4073                // Check max time
4074                hitMaxTime = ( CoinCpuTime() - model->getDblParam(CbcModel::CbcStartSeconds) >
4075                               model->getDblParam(CbcModel::CbcMaximumSeconds));
4076                if (hitMaxTime) {
4077                    // make sure rest are fast
4078                    doQuickly = true;
4079                    for ( int jDo = iDo + 1; jDo < numberToDo; jDo++) {
4080                        int iObject = whichObject[iDo];
4081                        OsiObject * object = model->modifiableObject(iObject);
4082                        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
4083                            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
4084                        if (dynamicObject)
4085                            dynamicObject->setNumberBeforeTrust(0);
4086                    }
4087                    numberTest = 0;
4088                }
4089                if (!choiceObject) {
4090                    delete choice.possibleBranch;
4091                }
4092            }
4093            if (model->messageHandler()->logLevel() > 3) {
4094                if (anyAction == -2) {
4095                    printf("infeasible\n");
4096                } else if (anyAction == -1) {
4097                    printf("%d fixed AND choosing %d iDo %d iChosenWhen %d numberToDo %d\n", numberToFix, bestChoice,
4098                           iDo, whichChoice, numberToDo);
4099                } else {
4100                    int iObject = whichObject[whichChoice];
4101                    OsiObject * object = model->modifiableObject(iObject);
4102                    CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
4103                        dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
4104                    if (dynamicObject) {
4105                        int iColumn = dynamicObject->columnNumber();
4106                        printf("choosing %d (column %d) iChosenWhen %d numberToDo %d\n", bestChoice,
4107                               iColumn, whichChoice, numberToDo);
4108                    }
4109                }
4110            }
4111            if (doneHotStart) {
4112                // Delete the snapshot
4113                solver->unmarkHotStart();
4114                // back to normal
4115                solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
4116                // restore basis
4117                solver->setWarmStart(ws);
4118            }
4119            solver->setIntParam(OsiMaxNumIterationHotStart, saveLimit);
4120            // Unless infeasible we will carry on
4121            // But we could fix anyway
4122            if (numberToFix && !hitMaxTime) {
4123                if (anyAction != -2) {
4124                    // apply and take off
4125                    bool feasible = true;
4126                    // can do quick optimality check
4127                    int easy = 2;
4128                    solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, &easy) ;
4129                    model->resolve(NULL, 11, saveSolution, saveLower, saveUpper) ;
4130                    //double newObjValue = solver->getObjSense()*solver->getObjValue();
4131                    //objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
4132                    solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
4133                    feasible = solver->isProvenOptimal();
4134                    if (feasible) {
4135                        anyAction = 0;
4136                        // See if candidate still possible
4137                        if (branch_) {
4138                            const OsiObject * object = model->object(bestChoice);
4139                            int preferredWay;
4140                            double infeasibility = object->infeasibility(&usefulInfo, preferredWay);
4141                            if (!infeasibility) {
4142                                // take out
4143                                delete branch_;
4144                                branch_ = NULL;
4145                            } else {
4146                                CbcBranchingObject * branchObj =
4147                                    dynamic_cast <CbcBranchingObject *>(branch_) ;
4148                                assert (branchObj);
4149                                branchObj->way(preferredWay);
4150                            }
4151                        }
4152                    } else {
4153                        anyAction = -2;
4154                        finished = true;
4155                    }
4156                }
4157            }
4158            // If  fixed then round again
4159            if (!branch_ && anyAction != -2) {
4160                finished = false;
4161            }
4162            delete ws;
4163        }
4164    }
4165    // update number of strong iterations etc
4166    model->incrementStrongInfo(numberStrongDone, numberStrongIterations,
4167                               anyAction == -2 ? 0 : numberToFix, anyAction == -2);
4168    if (model->searchStrategy() == -1) {
4169#ifndef COIN_DEVELOP
4170        if (solver->messageHandler()->logLevel() > 1)
4171#endif
4172            printf("%d strong, %d iters, %d inf, %d not finished, %d not trusted\n",
4173                   numberStrongDone, numberStrongIterations, numberStrongInfeasible, numberUnfinished,
4174                   numberNotTrusted);
4175        // decide what to do
4176        int strategy = 1;
4177        if (((numberUnfinished*4 > numberStrongDone &&
4178                numberStrongInfeasible*40 < numberStrongDone) ||
4179                numberStrongInfeasible < 0) && model->numberStrong() < 10 && model->numberBeforeTrust() <= 20 && model->numberObjects() > CoinMax(1000, solver->getNumRows())) {
4180            strategy = 2;
4181#ifdef COIN_DEVELOP
4182            //if (model->logLevel()>1)
4183            printf("going to strategy 2\n");
4184#endif
4185            // Weaken
4186            model->setNumberStrong(2);
4187            model->setNumberBeforeTrust(1);
4188            model->synchronizeNumberBeforeTrust();
4189        }
4190        if (numberNodes)
4191            strategy = 1;  // should only happen after hot start
4192        model->setSearchStrategy(strategy);
4193    } else if (numberStrongDone) {
4194        //printf("%d strongB, %d iters, %d inf, %d not finished, %d not trusted\n",
4195        //   numberStrongDone,numberStrongIterations,numberStrongInfeasible,numberUnfinished,
4196        //   numberNotTrusted);
4197    }
4198    if (model->searchStrategy() == 1 && numberNodes > 500 && numberNodes < -510) {
4199#ifndef COIN_DEVELOP
4200        if (solver->messageHandler()->logLevel() > 1)
4201#endif
4202            printf("after %d nodes - %d strong, %d iters, %d inf, %d not finished, %d not trusted\n",
4203                   numberNodes, numberStrongDone, numberStrongIterations, numberStrongInfeasible, numberUnfinished,
4204                   numberNotTrusted);
4205        // decide what to do
4206        if (numberUnfinished*10 > numberStrongDone + 1 ||
4207                !numberStrongInfeasible) {
4208            printf("going to strategy 2\n");
4209            // Weaken
4210            model->setNumberStrong(2);
4211            model->setNumberBeforeTrust(1);
4212            model->synchronizeNumberBeforeTrust();
4213            model->setSearchStrategy(2);
4214        }
4215    }
4216    // Set guessed solution value
4217    guessedObjectiveValue_ = objectiveValue_ + estimatedDegradation;
4218
4219    /*
4220      Cleanup, then we're finished
4221    */
4222    if (!model->branchingMethod())
4223        delete decision;
4224
4225    delete choiceObject;
4226    delete [] sort;
4227    delete [] whichObject;
4228#ifdef RANGING
4229    delete [] objectMark;
4230#endif
4231    delete [] saveLower;
4232    delete [] saveUpper;
4233    delete [] upEstimate;
4234    delete [] downEstimate;
4235# ifdef COIN_HAS_CLP
4236    if (osiclp) {
4237        osiclp->setSpecialOptions(saveClpOptions);
4238    }
4239# endif
4240    // restore solution
4241    solver->setColSolution(saveSolution);
4242    model->reserveCurrentSolution(saveSolution);
4243    delete [] saveSolution;
4244    model->setStateOfSearch(saveStateOfSearch);
4245    model->setLogLevel(saveLogLevel);
4246    // delete extra regions
4247    if (usefulInfo.usefulRegion_) {
4248        delete [] usefulInfo.usefulRegion_;
4249        delete [] usefulInfo.indexRegion_;
4250        delete [] usefulInfo.pi_;
4251        usefulInfo.usefulRegion_ = NULL;
4252        usefulInfo.indexRegion_ = NULL;
4253        usefulInfo.pi_ = NULL;
4254    }
4255    useShadow = model->moreSpecialOptions() & 7;
4256    if ((useShadow == 5 && model->getSolutionCount()) || useShadow == 6) {
4257        // zap pseudo shadow prices
4258        model->pseudoShadow(-1);
4259        // and switch off
4260        model->setMoreSpecialOptions(model->moreSpecialOptions()&(~1023));
4261    }
4262    return anyAction;
4263}
4264int CbcNode::analyze (CbcModel *model, double * results)
4265{
4266    int i;
4267    int numberIterationsAllowed = model->numberAnalyzeIterations();
4268    OsiSolverInterface * solver = model->solver();
4269    objectiveValue_ = solver->getObjSense() * solver->getObjValue();
4270    double cutoff = model->getCutoff();
4271    const double * lower = solver->getColLower();
4272    const double * upper = solver->getColUpper();
4273    const double * dj = solver->getReducedCost();
4274    int numberObjects = model->numberObjects();
4275    int numberColumns = model->getNumCols();
4276    // Initialize arrays
4277    int numberIntegers = model->numberIntegers();
4278    int * back = new int[numberColumns];
4279    const int * integerVariable = model->integerVariable();
4280    for (i = 0; i < numberColumns; i++)
4281        back[i] = -1;
4282    // What results is
4283    double * newLower = results;
4284    double * objLower = newLower + numberIntegers;
4285    double * newUpper = objLower + numberIntegers;
4286    double * objUpper = newUpper + numberIntegers;
4287    for (i = 0; i < numberIntegers; i++) {
4288        int iColumn = integerVariable[i];
4289        back[iColumn] = i;
4290        newLower[i] = 0.0;
4291        objLower[i] = -COIN_DBL_MAX;
4292        newUpper[i] = 0.0;
4293        objUpper[i] = -COIN_DBL_MAX;
4294    }
4295    double * saveUpper = new double[numberColumns];
4296    double * saveLower = new double[numberColumns];
4297    int anyAction = 0;
4298    // Save solution in case heuristics need good solution later
4299
4300    double * saveSolution = new double[numberColumns];
4301    memcpy(saveSolution, solver->getColSolution(), numberColumns*sizeof(double));
4302    model->reserveCurrentSolution(saveSolution);
4303    for (i = 0; i < numberColumns; i++) {
4304        saveLower[i] = lower[i];
4305        saveUpper[i] = upper[i];
4306    }
4307    // Get arrays to sort
4308    double * sort = new double[numberObjects];
4309    int * whichObject = new int[numberObjects];
4310    int numberToFix = 0;
4311    int numberToDo = 0;
4312    double integerTolerance =
4313        model->getDblParam(CbcModel::CbcIntegerTolerance);
4314    // point to useful information
4315    OsiBranchingInformation usefulInfo = model->usefulInformation();
4316    // and modify
4317    usefulInfo.depth_ = depth_;
4318
4319    // compute current state
4320    int numberObjectInfeasibilities; // just odd ones
4321    int numberIntegerInfeasibilities;
4322    model->feasibleSolution(
4323        numberIntegerInfeasibilities,
4324        numberObjectInfeasibilities);
4325# ifdef COIN_HAS_CLP
4326    OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
4327    int saveClpOptions = 0;
4328    bool fastIterations = (model->specialOptions() & 8) != 0;
4329    if (osiclp && fastIterations) {
4330        // for faster hot start
4331        saveClpOptions = osiclp->specialOptions();
4332        osiclp->setSpecialOptions(saveClpOptions | 8192);
4333    }
4334# else
4335    bool fastIterations = false ;
4336# endif
4337    /*
4338      Scan for branching objects that indicate infeasibility. Choose candidates
4339      using priority as the first criteria, then integer infeasibility.
4340
4341      The algorithm is to fill the array with a set of good candidates (by
4342      infeasibility) with priority bestPriority.  Finding a candidate with
4343      priority better (less) than bestPriority flushes the choice array. (This
4344      serves as initialization when the first candidate is found.)
4345
4346    */
4347    numberToDo = 0;
4348    for (i = 0; i < numberObjects; i++) {
4349        OsiObject * object = model->modifiableObject(i);
4350        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
4351            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
4352        if (!dynamicObject)
4353            continue;
4354        int preferredWay;
4355        double infeasibility = object->infeasibility(&usefulInfo, preferredWay);
4356        int iColumn = dynamicObject->columnNumber();
4357        if (saveUpper[iColumn] == saveLower[iColumn])
4358            continue;
4359        if (infeasibility)
4360            sort[numberToDo] = -1.0e10 - infeasibility;
4361        else
4362            sort[numberToDo] = -fabs(dj[iColumn]);
4363        whichObject[numberToDo++] = i;
4364    }
4365    // Save basis
4366    CoinWarmStart * ws = solver->getWarmStart();
4367    int saveLimit;
4368    solver->getIntParam(OsiMaxNumIterationHotStart, saveLimit);
4369    int targetIterations = CoinMax(500, numberIterationsAllowed / numberObjects);
4370    if (saveLimit < targetIterations)
4371        solver->setIntParam(OsiMaxNumIterationHotStart, targetIterations);
4372    // Mark hot start
4373    solver->markHotStart();
4374    // Sort
4375    CoinSort_2(sort, sort + numberToDo, whichObject);
4376    double * currentSolution = model->currentSolution();
4377    double objMin = 1.0e50;
4378    double objMax = -1.0e50;
4379    bool needResolve = false;
4380    int iDo;
4381    for (iDo = 0; iDo < numberToDo; iDo++) {
4382        CbcStrongInfo choice;
4383        int iObject = whichObject[iDo];
4384        OsiObject * object = model->modifiableObject(iObject);
4385        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
4386            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
4387        if (!dynamicObject)
4388            continue;
4389        int iColumn = dynamicObject->columnNumber();
4390        int preferredWay;
4391        object->infeasibility(&usefulInfo, preferredWay);
4392        double value = currentSolution[iColumn];
4393        double nearest = floor(value + 0.5);
4394        double lowerValue = floor(value);
4395        bool satisfied = false;
4396        if (fabs(value - nearest) <= integerTolerance || value < saveLower[iColumn] || value > saveUpper[iColumn]) {
4397            satisfied = true;
4398            double newValue;
4399            if (nearest < saveUpper[iColumn]) {
4400                newValue = nearest + 1.0001 * integerTolerance;
4401                lowerValue = nearest;
4402            } else {
4403                newValue = nearest - 1.0001 * integerTolerance;
4404                lowerValue = nearest - 1;
4405            }
4406            currentSolution[iColumn] = newValue;
4407        }
4408        double upperValue = lowerValue + 1.0;
4409        //CbcSimpleInteger * obj =
4410        //dynamic_cast <CbcSimpleInteger *>(object) ;
4411        //if (obj) {
4412        //choice.possibleBranch=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
4413        //} else {
4414        CbcObject * obj =
4415            dynamic_cast <CbcObject *>(object) ;
4416        assert (obj);
4417        choice.possibleBranch = obj->createCbcBranch(solver, &usefulInfo, preferredWay);
4418        //}
4419        currentSolution[iColumn] = value;
4420        // Save which object it was
4421        choice.objectNumber = iObject;
4422        choice.numIntInfeasUp = numberUnsatisfied_;
4423        choice.numIntInfeasDown = numberUnsatisfied_;
4424        choice.downMovement = 0.0;
4425        choice.upMovement = 0.0;
4426        choice.numItersDown = 0;
4427        choice.numItersUp = 0;
4428        choice.fix = 0; // say not fixed
4429        double objectiveChange ;
4430        double newObjectiveValue = 1.0e100;
4431        int j;
4432        // status is 0 finished, 1 infeasible and other
4433        int iStatus;
4434        /*
4435          Try the down direction first. (Specify the initial branching alternative as
4436          down with a call to way(-1). Each subsequent call to branch() performs the
4437          specified branch and advances the branch object state to the next branch
4438          alternative.)
4439        */
4440        choice.possibleBranch->way(-1) ;
4441        choice.possibleBranch->branch() ;
4442        if (fabs(value - lowerValue) > integerTolerance) {
4443            solver->solveFromHotStart() ;
4444            /*
4445              We now have an estimate of objective degradation that we can use for strong
4446              branching. If we're over the cutoff, the variable is monotone up.
4447              If we actually made it to optimality, check for a solution, and if we have
4448              a good one, call setBestSolution to process it. Note that this may reduce the
4449              cutoff, so we check again to see if we can declare this variable monotone.
4450            */
4451            if (solver->isProvenOptimal())
4452                iStatus = 0; // optimal
4453            else if (solver->isIterationLimitReached()
4454                     && !solver->isDualObjectiveLimitReached())
4455                iStatus = 2; // unknown
4456            else
4457                iStatus = 1; // infeasible
4458            newObjectiveValue = solver->getObjSense() * solver->getObjValue();
4459            choice.numItersDown = solver->getIterationCount();
4460            numberIterationsAllowed -= choice.numItersDown;
4461            objectiveChange = newObjectiveValue  - objectiveValue_;
4462            if (!iStatus) {
4463                choice.finishedDown = true ;
4464                if (newObjectiveValue >= cutoff) {
4465                    objectiveChange = 1.0e100; // say infeasible
4466                } else {
4467                    // See if integer solution
4468                    if (model->feasibleSolution(choice.numIntInfeasDown,
4469                                                choice.numObjInfeasDown)
4470                            && model->problemFeasibility()->feasible(model, -1) >= 0) {
4471                        model->setBestSolution(CBC_STRONGSOL,
4472                                               newObjectiveValue,
4473                                               solver->getColSolution()) ;
4474                        model->setLastHeuristic(NULL);
4475                        model->incrementUsed(solver->getColSolution());
4476                        cutoff = model->getCutoff();
4477                        if (newObjectiveValue >= cutoff)        //  *new* cutoff
4478                            objectiveChange = 1.0e100 ;
4479                    }
4480                }
4481            } else if (iStatus == 1) {
4482                objectiveChange = 1.0e100 ;
4483            } else {
4484                // Can't say much as we did not finish
4485                choice.finishedDown = false ;
4486            }
4487            choice.downMovement = objectiveChange ;
4488        }
4489        // restore bounds
4490        for ( j = 0; j < numberColumns; j++) {
4491            if (saveLower[j] != lower[j])
4492                solver->setColLower(j, saveLower[j]);
4493            if (saveUpper[j] != upper[j])
4494                solver->setColUpper(j, saveUpper[j]);
4495        }
4496        // repeat the whole exercise, forcing the variable up
4497        choice.possibleBranch->branch();
4498        if (fabs(value - upperValue) > integerTolerance) {
4499            solver->solveFromHotStart() ;
4500            /*
4501              We now have an estimate of objective degradation that we can use for strong
4502              branching. If we're over the cutoff, the variable is monotone up.
4503              If we actually made it to optimality, check for a solution, and if we have
4504              a good one, call setBestSolution to process it. Note that this may reduce the
4505              cutoff, so we check again to see if we can declare this variable monotone.
4506            */
4507            if (solver->isProvenOptimal())
4508                iStatus = 0; // optimal
4509            else if (solver->isIterationLimitReached()
4510                     && !solver->isDualObjectiveLimitReached())
4511                iStatus = 2; // unknown
4512            else
4513                iStatus = 1; // infeasible
4514            newObjectiveValue = solver->getObjSense() * solver->getObjValue();
4515            choice.numItersUp = solver->getIterationCount();
4516            numberIterationsAllowed -= choice.numItersUp;
4517            objectiveChange = newObjectiveValue  - objectiveValue_;
4518            if (!iStatus) {
4519                choice.finishedUp = true ;
4520                if (newObjectiveValue >= cutoff) {
4521                    objectiveChange = 1.0e100; // say infeasible
4522                } else {
4523                    // See if integer solution
4524                    if (model->feasibleSolution(choice.numIntInfeasUp,
4525                                                choice.numObjInfeasUp)
4526                            && model->problemFeasibility()->feasible(model, -1) >= 0) {
4527                        model->setBestSolution(CBC_STRONGSOL,
4528                                               newObjectiveValue,
4529                                               solver->getColSolution()) ;
4530                        model->setLastHeuristic(NULL);
4531                        model->incrementUsed(solver->getColSolution());
4532                        cutoff = model->getCutoff();
4533                        if (newObjectiveValue >= cutoff)        //  *new* cutoff
4534                            objectiveChange = 1.0e100 ;
4535                    }
4536                }
4537            } else if (iStatus == 1) {
4538                objectiveChange = 1.0e100 ;
4539            } else {
4540                // Can't say much as we did not finish
4541                choice.finishedUp = false ;
4542            }
4543            choice.upMovement = objectiveChange ;
4544
4545            // restore bounds
4546            for ( j = 0; j < numberColumns; j++) {
4547                if (saveLower[j] != lower[j])
4548                    solver->setColLower(j, saveLower[j]);
4549                if (saveUpper[j] != upper[j])
4550                    solver->setColUpper(j, saveUpper[j]);
4551            }
4552        }
4553        // If objective goes above certain amount we can set bound
4554        int jInt = back[iColumn];
4555        newLower[jInt] = upperValue;
4556        if (choice.finishedDown)
4557            objLower[jInt] = choice.downMovement + objectiveValue_;
4558        else
4559            objLower[jInt] = objectiveValue_;
4560        newUpper[jInt] = lowerValue;
4561        if (choice.finishedUp)
4562            objUpper[jInt] = choice.upMovement + objectiveValue_;
4563        else
4564            objUpper[jInt] = objectiveValue_;
4565        objMin = CoinMin(CoinMin(objLower[jInt], objUpper[jInt]), objMin);
4566        /*
4567          End of evaluation for this candidate variable. Possibilities are:
4568          * Both sides below cutoff; this variable is a candidate for branching.
4569          * Both sides infeasible or above the objective cutoff: no further action
4570          here. Break from the evaluation loop and assume the node will be purged
4571          by the caller.
4572          * One side below cutoff: Install the branch (i.e., fix the variable). Break
4573          from the evaluation loop and assume the node will be reoptimised by the
4574          caller.
4575        */
4576        if (choice.upMovement < 1.0e100) {
4577            if (choice.downMovement < 1.0e100) {
4578                objMax = CoinMax(CoinMax(objLower[jInt], objUpper[jInt]), objMax);
4579                // In case solution coming in was odd
4580                choice.upMovement = CoinMax(0.0, choice.upMovement);
4581                choice.downMovement = CoinMax(0.0, choice.downMovement);
4582                // feasible -
4583                model->messageHandler()->message(CBC_STRONG, *model->messagesPointer())
4584                << iObject << iColumn
4585                << choice.downMovement << choice.numIntInfeasDown
4586                << choice.upMovement << choice.numIntInfeasUp
4587                << value
4588                << CoinMessageEol;
4589            } else {
4590                // up feasible, down infeasible
4591                anyAction = -1;
4592                if (!satisfied)
4593                    needResolve = true;
4594                choice.fix = 1;
4595                numberToFix++;
4596                saveLower[iColumn] = upperValue;
4597                solver->setColLower(iColumn, upperValue);
4598            }
4599        } else {
4600            if (choice.downMovement < 1.0e100) {
4601                // down feasible, up infeasible
4602                anyAction = -1;
4603                if (!satisfied)
4604                    needResolve = true;
4605                choice.fix = -1;
4606                numberToFix++;
4607                saveUpper[iColumn] = lowerValue;
4608                solver->setColUpper(iColumn, lowerValue);
4609            } else {
4610                // neither side feasible
4611                anyAction = -2;
4612                printf("Both infeasible for choice %d sequence %d\n", i,
4613                       model->object(choice.objectNumber)->columnNumber());
4614                delete ws;
4615                ws = NULL;
4616                //solver->writeMps("bad");
4617                numberToFix = -1;
4618                delete choice.possibleBranch;
4619                choice.possibleBranch = NULL;
4620                break;
4621            }
4622        }
4623        delete choice.possibleBranch;
4624        if (numberIterationsAllowed <= 0)
4625            break;
4626        //printf("obj %d, col %d, down %g up %g value %g\n",iObject,iColumn,
4627        //     choice.downMovement,choice.upMovement,value);
4628    }
4629    printf("Best possible solution %g, can fix more if solution of %g found - looked at %d variables in %d iterations\n",
4630           objMin, objMax, iDo, model->numberAnalyzeIterations() - numberIterationsAllowed);
4631    model->setNumberAnalyzeIterations(numberIterationsAllowed);
4632    // Delete the snapshot
4633    solver->unmarkHotStart();
4634    // back to normal
4635    solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
4636    solver->setIntParam(OsiMaxNumIterationHotStart, saveLimit);
4637    // restore basis
4638    solver->setWarmStart(ws);
4639    delete ws;
4640
4641    delete [] sort;
4642    delete [] whichObject;
4643    delete [] saveLower;
4644    delete [] saveUpper;
4645    delete [] back;
4646    // restore solution
4647    solver->setColSolution(saveSolution);
4648# ifdef COIN_HAS_CLP
4649    if (osiclp)
4650        osiclp->setSpecialOptions(saveClpOptions);
4651# endif
4652    model->reserveCurrentSolution(saveSolution);
4653    delete [] saveSolution;
4654    if (needResolve)
4655        solver->resolve();
4656    return numberToFix;
4657}
4658
4659
4660CbcNode::CbcNode(const CbcNode & rhs)
4661        : CoinTreeNode(rhs)
4662{
4663#ifdef CHECK_NODE
4664    printf("CbcNode %x Constructor from rhs %x\n", this, &rhs);
4665#endif
4666    if (rhs.nodeInfo_)
4667        nodeInfo_ = rhs.nodeInfo_->clone();
4668    else
4669        nodeInfo_ = NULL;
4670    objectiveValue_ = rhs.objectiveValue_;
4671    guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
4672    sumInfeasibilities_ = rhs.sumInfeasibilities_;
4673    if (rhs.branch_)
4674        branch_ = rhs.branch_->clone();
4675    else
4676        branch_ = NULL;
4677    depth_ = rhs.depth_;
4678    numberUnsatisfied_ = rhs.numberUnsatisfied_;
4679    nodeNumber_ = rhs.nodeNumber_;
4680    state_ = rhs.state_;
4681    if (nodeInfo_)
4682        assert ((state_&2) != 0);
4683    else
4684        assert ((state_&2) == 0);
4685}
4686
4687CbcNode &
4688CbcNode::operator=(const CbcNode & rhs)
4689{
4690    if (this != &rhs) {
4691        delete nodeInfo_;
4692        if (rhs.nodeInfo_)
4693            nodeInfo_ = rhs.nodeInfo_->clone();
4694        else
4695            nodeInfo_ = NULL;
4696        objectiveValue_ = rhs.objectiveValue_;
4697        guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
4698        sumInfeasibilities_ = rhs.sumInfeasibilities_;
4699        if (rhs.branch_)
4700            branch_ = rhs.branch_->clone();
4701        else
4702            branch_ = NULL,
4703                      depth_ = rhs.depth_;
4704        numberUnsatisfied_ = rhs.numberUnsatisfied_;
4705        nodeNumber_ = rhs.nodeNumber_;
4706        state_ = rhs.state_;
4707        if (nodeInfo_)
4708            assert ((state_&2) != 0);
4709        else
4710            assert ((state_&2) == 0);
4711    }
4712    return *this;
4713}
4714CbcNode::~CbcNode ()
4715{
4716#ifdef CHECK_NODE
4717    if (nodeInfo_) {
4718        printf("CbcNode %x Destructor nodeInfo %x (%d)\n",
4719               this, nodeInfo_, nodeInfo_->numberPointingToThis());
4720        //assert(nodeInfo_->numberPointingToThis()>=0);
4721    } else {
4722        printf("CbcNode %x Destructor nodeInfo %x (?)\n",
4723               this, nodeInfo_);
4724    }
4725#endif
4726    if (nodeInfo_) {
4727        // was if (nodeInfo_&&(state_&2)!=0) {
4728        nodeInfo_->nullOwner();
4729        int numberToDelete = nodeInfo_->numberBranchesLeft();
4730        //    CbcNodeInfo * parent = nodeInfo_->parent();
4731        //assert (nodeInfo_->numberPointingToThis()>0);
4732        if (nodeInfo_->decrement(numberToDelete) == 0 || (state_&2) == 0) {
4733            if ((state_&2) == 0)
4734                nodeInfo_->nullParent();
4735            delete nodeInfo_;
4736        } else {
4737            //printf("node %x nodeinfo %x parent %x\n",this,nodeInfo_,nodeInfo_->parent());
4738            // anyway decrement parent
4739            //if (parent)
4740            ///parent->decrement(1);
4741        }
4742    }
4743    delete branch_;
4744}
4745// Decrement  active cut counts
4746void
4747CbcNode::decrementCuts(int change)
4748{
4749    if (nodeInfo_)
4750        assert ((state_&2) != 0);
4751    else
4752        assert ((state_&2) == 0);
4753    if (nodeInfo_) {
4754        nodeInfo_->decrementCuts(change);
4755    }
4756}
4757void
4758CbcNode::decrementParentCuts(CbcModel * model, int change)
4759{
4760    if (nodeInfo_)
4761        assert ((state_&2) != 0);
4762    else
4763        assert ((state_&2) == 0);
4764    if (nodeInfo_) {
4765        nodeInfo_->decrementParentCuts(model, change);
4766    }
4767}
4768
4769/*
4770  Initialize reference counts (numberPointingToThis, numberBranchesLeft_)
4771  in the attached nodeInfo_.
4772*/
4773void
4774CbcNode::initializeInfo()
4775{
4776    assert(nodeInfo_ && branch_) ;
4777    nodeInfo_->initializeInfo(branch_->numberBranches());
4778    assert ((state_&2) != 0);
4779    assert (nodeInfo_->numberBranchesLeft() ==
4780            branch_->numberBranchesLeft());
4781}
4782// Nulls out node info
4783void
4784CbcNode::nullNodeInfo()
4785{
4786    nodeInfo_ = NULL;
4787    // say not active
4788    state_ &= ~2;
4789}
4790
4791int
4792CbcNode::branch(OsiSolverInterface * solver)
4793{
4794    double changeInGuessed;
4795    assert (nodeInfo_->numberBranchesLeft() ==
4796            branch_->numberBranchesLeft());
4797    if (!solver)
4798        changeInGuessed = branch_->branch();
4799    else
4800        changeInGuessed = branch_->branch(solver);
4801    guessedObjectiveValue_ += changeInGuessed;
4802    //#define PRINTIT
4803#ifdef PRINTIT
4804    int numberLeft = nodeInfo_->numberBranchesLeft();
4805    CbcNodeInfo * parent = nodeInfo_->parent();
4806    int parentNodeNumber = -1;
4807    CbcBranchingObject * object1 =
4808        dynamic_cast<CbcBranchingObject *>(branch_) ;
4809    //OsiObject * object = object1->
4810    //int sequence = object->columnNumber);
4811    int id = -1;
4812    double value = 0.0;
4813    if (object1) {
4814        id = object1->variable();
4815        value = object1->value();
4816    }
4817    printf("id %d value %g objvalue %g\n", id, value, objectiveValue_);
4818    if (parent)
4819        parentNodeNumber = parent->nodeNumber();
4820    printf("Node number %d, %s, way %d, depth %d, parent node number %d\n",
4821           nodeInfo_->nodeNumber(), (numberLeft == 2) ? "leftBranch" : "rightBranch",
4822           way(), depth_, parentNodeNumber);
4823    assert (parentNodeNumber != nodeInfo_->nodeNumber());
4824#endif
4825    return nodeInfo_->branchedOn();
4826}
4827/* Active arm of the attached OsiBranchingObject.
4828
4829   In the simplest instance, coded -1 for the down arm of the branch, +1 for
4830   the up arm. But see OsiBranchingObject::way()
4831   Use nodeInfo--.numberBranchesLeft_ to see how active
4832*/
4833int
4834CbcNode::way() const
4835{
4836    if (branch_) {
4837        CbcBranchingObject * obj =
4838            dynamic_cast <CbcBranchingObject *>(branch_) ;
4839        if (obj) {
4840            return obj->way();
4841        } else {
4842            OsiTwoWayBranchingObject * obj2 =
4843                dynamic_cast <OsiTwoWayBranchingObject *>(branch_) ;
4844            assert (obj2);
4845            return obj2->way();
4846        }
4847    } else {
4848        return 0;
4849    }
4850}
4851/* Create a branching object for the node
4852
4853   The routine scans the object list of the model and selects a set of
4854   unsatisfied objects as candidates for branching. The candidates are
4855   evaluated, and an appropriate branch object is installed.
4856
4857   The numberPassesLeft is decremented to stop fixing one variable each time
4858   and going on and on (e.g. for stock cutting, air crew scheduling)
4859
4860   If evaluation determines that an object is monotone or infeasible,
4861   the routine returns immediately. In the case of a monotone object,
4862   the branch object has already been called to modify the model.
4863
4864   Return value:
4865   <ul>
4866   <li>  0: A branching object has been installed
4867   <li> -1: A monotone object was discovered
4868   <li> -2: An infeasible object was discovered
4869   </ul>
4870   Branch state:
4871   <ul>
4872   <li> -1: start
4873   <li> -1: A monotone object was discovered
4874   <li> -2: An infeasible object was discovered
4875   </ul>
4876*/
4877int
4878CbcNode::chooseOsiBranch (CbcModel * model,
4879                          CbcNode * lastNode,
4880                          OsiBranchingInformation * usefulInfo,
4881                          int branchState)
4882{
4883    int returnStatus = 0;
4884    if (lastNode)
4885        depth_ = lastNode->depth_ + 1;
4886    else
4887        depth_ = 0;
4888    OsiSolverInterface * solver = model->solver();
4889    objectiveValue_ = solver->getObjValue() * solver->getObjSense();
4890    usefulInfo->objectiveValue_ = objectiveValue_;
4891    usefulInfo->depth_ = depth_;
4892    const double * saveInfoSol = usefulInfo->solution_;
4893    double * saveSolution = new double[solver->getNumCols()];
4894    memcpy(saveSolution, solver->getColSolution(), solver->getNumCols()*sizeof(double));
4895    usefulInfo->solution_ = saveSolution;
4896    OsiChooseVariable * choose = model->branchingMethod()->chooseMethod();
4897    int numberUnsatisfied = -1;
4898    if (branchState < 0) {
4899        // initialize
4900        // initialize sum of "infeasibilities"
4901        sumInfeasibilities_ = 0.0;
4902        numberUnsatisfied = choose->setupList(usefulInfo, true);
4903        numberUnsatisfied_ = numberUnsatisfied;
4904        branchState = 0;
4905        if (numberUnsatisfied_ < 0) {
4906            // infeasible
4907            delete [] saveSolution;
4908            return -2;
4909        }
4910    }
4911    // unset best
4912    int best = -1;
4913    choose->setBestObjectIndex(-1);
4914    if (numberUnsatisfied) {
4915        if (branchState > 0 || !choose->numberOnList()) {
4916            // we need to return at once - don't do strong branching or anything
4917            if (choose->numberOnList() || !choose->numberStrong()) {
4918                best = choose->candidates()[0];
4919                choose->setBestObjectIndex(best);
4920            } else {
4921                // nothing on list - need to try again - keep any solution
4922                numberUnsatisfied = choose->setupList(usefulInfo, false);
4923                numberUnsatisfied_ = numberUnsatisfied;
4924                if (numberUnsatisfied) {
4925                    best = choose->candidates()[0];
4926                    choose->setBestObjectIndex(best);
4927                }
4928            }
4929        } else {
4930            // carry on with strong branching or whatever
4931            int returnCode = choose->chooseVariable(solver, usefulInfo, true);
4932            // update number of strong iterations etc
4933            model->incrementStrongInfo(choose->numberStrongDone(), choose->numberStrongIterations(),
4934                                       returnCode == -1 ? 0 : choose->numberStrongFixed(), returnCode == -1);
4935            if (returnCode > 1) {
4936                // has fixed some
4937                returnStatus = -1;
4938            } else if (returnCode == -1) {
4939                // infeasible
4940                returnStatus = -2;
4941            } else if (returnCode == 0) {
4942                // normal
4943                returnStatus = 0;
4944                numberUnsatisfied = 1;
4945            } else {
4946                // ones on list satisfied - double check
4947                numberUnsatisfied = choose->setupList(usefulInfo, false);
4948                numberUnsatisfied_ = numberUnsatisfied;
4949                if (numberUnsatisfied) {
4950                    best = choose->candidates()[0];
4951                    choose->setBestObjectIndex(best);
4952                }
4953            }
4954        }
4955    }
4956    delete branch_;
4957    branch_ = NULL;
4958    guessedObjectiveValue_ = COIN_DBL_MAX;//objectiveValue_; // for now
4959    if (!returnStatus) {
4960        if (numberUnsatisfied) {
4961            // create branching object
4962            const OsiObject * obj = model->solver()->object(choose->bestObjectIndex());
4963            //const OsiSolverInterface * solver = usefulInfo->solver_;
4964            branch_ = obj->createBranch(model->solver(), usefulInfo, obj->whichWay());
4965        }
4966    }
4967    usefulInfo->solution_ = saveInfoSol;
4968    delete [] saveSolution;
4969    // may have got solution
4970    if (choose->goodSolution()
4971            && model->problemFeasibility()->feasible(model, -1) >= 0) {
4972        // yes
4973        double objValue = choose->goodObjectiveValue();
4974        model->setBestSolution(CBC_STRONGSOL,
4975                               objValue,
4976                               choose->goodSolution()) ;
4977        model->setLastHeuristic(NULL);
4978        model->incrementUsed(choose->goodSolution());
4979        choose->clearGoodSolution();
4980    }
4981    return returnStatus;
4982}
4983int
4984CbcNode::chooseClpBranch (CbcModel * model,
4985                          CbcNode * lastNode)
4986{
4987    assert(lastNode);
4988    depth_ = lastNode->depth_ + 1;
4989    delete branch_;
4990    branch_ = NULL;
4991    OsiSolverInterface * solver = model->solver();
4992    //double saveObjectiveValue = solver->getObjValue();
4993    //double objectiveValue = CoinMax(solver->getObjSense()*saveObjectiveValue,objectiveValue_);
4994    const double * lower = solver->getColLower();
4995    const double * upper = solver->getColUpper();
4996    // point to useful information
4997    OsiBranchingInformation usefulInfo = model->usefulInformation();
4998    // and modify
4999    usefulInfo.depth_ = depth_;
5000    int i;
5001    //bool beforeSolution = model->getSolutionCount()==0;
5002    int numberObjects = model->numberObjects();
5003    int numberColumns = model->getNumCols();
5004    double * saveUpper = new double[numberColumns];
5005    double * saveLower = new double[numberColumns];
5006
5007    // Save solution in case heuristics need good solution later
5008
5009    double * saveSolution = new double[numberColumns];
5010    memcpy(saveSolution, solver->getColSolution(), numberColumns*sizeof(double));
5011    model->reserveCurrentSolution(saveSolution);
5012    for (i = 0; i < numberColumns; i++) {
5013        saveLower[i] = lower[i];
5014        saveUpper[i] = upper[i];
5015    }
5016    // Save basis
5017    CoinWarmStart * ws = solver->getWarmStart();
5018    numberUnsatisfied_ = 0;
5019    // initialize sum of "infeasibilities"
5020    sumInfeasibilities_ = 0.0;
5021    // Note looks as if off end (hidden one)
5022    OsiObject * object = model->modifiableObject(numberObjects);
5023    CbcGeneralDepth * thisOne = dynamic_cast <CbcGeneralDepth *> (object);
5024    assert (thisOne);
5025    OsiClpSolverInterface * clpSolver
5026    = dynamic_cast<OsiClpSolverInterface *> (solver);
5027    assert (clpSolver);
5028    ClpSimplex * simplex = clpSolver->getModelPtr();
5029    int preferredWay;
5030    double infeasibility = object->infeasibility(&usefulInfo, preferredWay);
5031    if (thisOne->whichSolution() >= 0) {
5032        ClpNode * nodeInfo = thisOne->nodeInfo(thisOne->whichSolution());
5033        nodeInfo->applyNode(simplex, 2);
5034        int saveLogLevel = simplex->logLevel();
5035        simplex->setLogLevel(0);
5036        simplex->dual();
5037        simplex->setLogLevel(saveLogLevel);
5038        double cutoff = model->getCutoff();
5039        bool goodSolution = true;
5040        if (simplex->status()) {
5041            //simplex->writeMps("bad7.mps",2);
5042            if (nodeInfo->objectiveValue() > cutoff - 1.0e-2)
5043                goodSolution = false;
5044            else
5045                assert (!simplex->status());
5046        }
5047        if (goodSolution) {
5048            double newObjectiveValue = solver->getObjSense() * solver->getObjValue();
5049            // See if integer solution
5050            int numInf;
5051            int numInf2;
5052            bool gotSol = model->feasibleSolution(numInf, numInf2);
5053            if (!gotSol) {
5054                printf("numinf %d\n", numInf);
5055                double * sol = simplex->primalColumnSolution();
5056                for (int i = 0; i < numberColumns; i++) {
5057                    if (simplex->isInteger(i)) {
5058                        double value = floor(sol[i] + 0.5);
5059                        if (fabs(value - sol[i]) > 1.0e-7) {
5060                            printf("%d value %g\n", i, sol[i]);
5061                            if (fabs(value - sol[i]) < 1.0e-3) {
5062                                sol[i] = value;
5063                            }
5064                        }
5065                    }
5066                }
5067                simplex->writeMps("bad8.mps", 2);
5068                bool gotSol = model->feasibleSolution(numInf, numInf2);
5069                if (!gotSol)
5070                    assert (gotSol);
5071            }
5072            model->setBestSolution(CBC_STRONGSOL,
5073                                   newObjectiveValue,
5074                                   solver->getColSolution()) ;
5075            model->setLastHeuristic(NULL);
5076            model->incrementUsed(solver->getColSolution());
5077        }
5078    }
5079    // restore bounds
5080    {
5081        for (int j = 0; j < numberColumns; j++) {
5082            if (saveLower[j] != lower[j])
5083                solver->setColLower(j, saveLower[j]);
5084            if (saveUpper[j] != upper[j])
5085                solver->setColUpper(j, saveUpper[j]);
5086        }
5087    }
5088    // restore basis
5089    solver->setWarmStart(ws);
5090    delete ws;
5091    int anyAction;
5092    //#define CHECK_PATH
5093#ifdef CHECK_PATH
5094    extern int gotGoodNode_Z;
5095    if (gotGoodNode_Z >= 0)
5096        printf("good node %d %g\n", gotGoodNode_Z, infeasibility);
5097#endif
5098    if (infeasibility > 0.0) {
5099        if (infeasibility == COIN_DBL_MAX) {
5100            anyAction = -2; // infeasible
5101        } else {
5102            branch_ = thisOne->createCbcBranch(solver, &usefulInfo, preferredWay);
5103            // Set to firat one (and change when re-pushing)
5104            CbcGeneralBranchingObject * branch =
5105                dynamic_cast <CbcGeneralBranchingObject *> (branch_);
5106            branch->state(objectiveValue_, sumInfeasibilities_,
5107                          numberUnsatisfied_, 0);
5108            branch->setNode(this);
5109            anyAction = 0;
5110        }
5111    } else {
5112        anyAction = -1;
5113    }
5114#ifdef CHECK_PATH
5115    gotGoodNode_Z = -1;
5116#endif
5117    // Set guessed solution value
5118    guessedObjectiveValue_ = objectiveValue_ + 1.0e-5;
5119    delete [] saveLower;
5120    delete [] saveUpper;
5121
5122    // restore solution
5123    solver->setColSolution(saveSolution);
5124    delete [] saveSolution;
5125    return anyAction;
5126}
5127/* Double checks in case node can change its mind!
5128   Returns objective value
5129   Can change objective etc */
5130double
5131CbcNode::checkIsCutoff(double cutoff)
5132{
5133    branch_->checkIsCutoff(cutoff);
5134    return objectiveValue_;
5135}
Note: See TracBrowser for help on using the repository browser.