source: trunk/Cbc/src/CbcNode.cpp @ 2092

Last change on this file since 2092 was 2092, checked in by forrest, 5 years ago

symmetry and diving

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