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

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

First try at orbital branching

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 210.1 KB
Line 
1/* $Id: CbcNode.cpp 2048 2014-07-16 09:29:16Z 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                        int preferredWay;
1418                        double infeasibility = object->infeasibility(&usefulInfo, preferredWay);
1419                        if (!infeasibility) {
1420                            // take out
1421                            delete thisChoice.possibleBranch;
1422                        } else {
1423                            choice[numberLeft++] = thisChoice;
1424                        }
1425                    }
1426                    numberStrong = numberLeft;
1427                    for (; i < maximumStrong; i++) {
1428                        delete choice[i].possibleBranch;
1429                        choice[i].possibleBranch = NULL;
1430                    }
1431                    // If all fixed then round again
1432                    if (!numberLeft) {
1433                        finished = false;
1434                        numberStrong = 0;
1435                        saveNumberStrong = 0;
1436                        maximumStrong = 1;
1437                    } else {
1438                        anyAction = 0;
1439                    }
1440                    // If these two uncommented then different action
1441                    anyAction = -1;
1442                    finished = true;
1443                    //printf("some fixed but continuing %d left\n",numberLeft);
1444                } else {
1445                    anyAction = -2; // say infeasible
1446                }
1447            }
1448            delete ws;
1449            //int numberNodes = model->getNodeCount();
1450            // update number of strong iterations etc
1451            model->incrementStrongInfo(numberStrongDone, numberStrongIterations,
1452                                       anyAction == -2 ? 0 : numberStrongInfeasible, anyAction == -2);
1453
1454            /*
1455              anyAction >= 0 indicates that strong branching didn't produce any monotone
1456              variables. Sift through the candidates for the best one.
1457
1458              QUERY: Setting numberNodes looks to be a distributed noop. numberNodes is
1459              local to this code block. Perhaps should be numberNodes_ from model?
1460              Unclear what this calculation is doing.
1461            */
1462            if (anyAction >= 0) {
1463
1464                // get average cost per iteration and assume stopped ones
1465                // would stop after 50% more iterations at average cost??? !!! ???
1466                double averageCostPerIteration = 0.0;
1467                double totalNumberIterations = 1.0;
1468                int smallestNumberInfeasibilities = COIN_INT_MAX;
1469                for (i = 0; i < numberStrong; i++) {
1470                    totalNumberIterations += choice[i].numItersDown +
1471                                             choice[i].numItersUp ;
1472                    averageCostPerIteration += choice[i].downMovement +
1473                                               choice[i].upMovement;
1474                    smallestNumberInfeasibilities =
1475                        CoinMin(CoinMin(choice[i].numIntInfeasDown ,
1476                                        choice[i].numIntInfeasUp ),
1477                                smallestNumberInfeasibilities);
1478                }
1479                //if (smallestNumberInfeasibilities>=numberIntegerInfeasibilities)
1480                //numberNodes=1000000; // switch off search for better solution
1481                averageCostPerIteration /= totalNumberIterations;
1482                // all feasible - choose best bet
1483
1484                // New method does all at once so it can be more sophisticated
1485                // in deciding how to balance actions.
1486                // But it does need arrays
1487                double * changeUp = new double [numberStrong];
1488                int * numberInfeasibilitiesUp = new int [numberStrong];
1489                double * changeDown = new double [numberStrong];
1490                int * numberInfeasibilitiesDown = new int [numberStrong];
1491                CbcBranchingObject ** objects = new CbcBranchingObject * [ numberStrong];
1492                for (i = 0 ; i < numberStrong ; i++) {
1493                    int iColumn = choice[i].possibleBranch->variable() ;
1494                    model->messageHandler()->message(CBC_STRONG, *model->messagesPointer())
1495                    << i << iColumn
1496                    << choice[i].downMovement << choice[i].numIntInfeasDown
1497                    << choice[i].upMovement << choice[i].numIntInfeasUp
1498                    << choice[i].possibleBranch->value()
1499                    << CoinMessageEol;
1500                    changeUp[i] = choice[i].upMovement;
1501                    numberInfeasibilitiesUp[i] = choice[i].numIntInfeasUp;
1502                    changeDown[i] = choice[i].downMovement;
1503                    numberInfeasibilitiesDown[i] = choice[i].numIntInfeasDown;
1504                    objects[i] = choice[i].possibleBranch;
1505                }
1506                int whichObject = decision->bestBranch(objects, numberStrong, numberUnsatisfied_,
1507                                                       changeUp, numberInfeasibilitiesUp,
1508                                                       changeDown, numberInfeasibilitiesDown,
1509                                                       objectiveValue_);
1510                // move branching object and make sure it will not be deleted
1511                if (whichObject >= 0) {
1512                    branch_ = objects[whichObject];
1513                    if (model->messageHandler()->logLevel() > 3)
1514                        printf("Choosing column %d\n", choice[whichObject].possibleBranch->variable()) ;
1515                    choice[whichObject].possibleBranch = NULL;
1516                }
1517                delete [] changeUp;
1518                delete [] numberInfeasibilitiesUp;
1519                delete [] changeDown;
1520                delete [] numberInfeasibilitiesDown;
1521                delete [] objects;
1522            }
1523#     ifdef COIN_HAS_CLP
1524            if (osiclp && !allNormal) {
1525                // back to normal
1526                osiclp->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
1527            }
1528#     endif
1529        }
1530        /*
1531          Simple branching. Probably just one, but we may have got here
1532          because of an odd branch e.g. a cut
1533        */
1534        else {
1535            // not strong
1536            // C) create branching object
1537            branch_ = choice[bestChoice].possibleBranch;
1538            choice[bestChoice].possibleBranch = NULL;
1539        }
1540    }
1541    // Set guessed solution value
1542    guessedObjectiveValue_ = objectiveValue_ + estimatedDegradation;
1543    //printf ("Node %d depth %d unsatisfied %d sum %g obj %g guess %g\n",
1544    //      model->getNodeCount(),depth_,numberUnsatisfied_,
1545    //      sumInfeasibilities_,objectiveValue_,guessedObjectiveValue_);
1546    /*
1547      Cleanup, then we're outta here.
1548    */
1549    if (!model->branchingMethod() || dynamicBranchingObject)
1550        delete decision;
1551
1552    for (i = 0; i < maximumStrong; i++)
1553        delete choice[i].possibleBranch;
1554    delete [] choice;
1555    delete [] saveLower;
1556    delete [] saveUpper;
1557
1558    // restore solution
1559    solver->setColSolution(saveSolution);
1560    delete [] saveSolution;
1561# ifdef COIN_HAS_CLP
1562    if (osiclp)
1563        osiclp->setSpecialOptions(saveClpOptions);
1564# endif
1565    return anyAction;
1566}
1567
1568/*
1569  Version for dynamic pseudo costs.
1570
1571  **** For now just return if anything odd
1572  later allow even if odd
1573
1574  The routine scans through the object list of the model looking for objects
1575  that indicate infeasibility. It tests each object using strong branching
1576  and selects the one with the least objective degradation.  A corresponding
1577  branching object is left attached to lastNode.
1578  This version gives preference in evaluation to variables which
1579  have not been evaluated many times.  It also uses numberStrong
1580  to say give up if last few tries have not changed incumbent.
1581  See Achterberg, Koch and Martin.
1582
1583  If strong branching is disabled, a candidate object is chosen essentially
1584  at random (whatever object ends up in pos'n 0 of the candidate array).
1585
1586  If a branching candidate is found to be monotone, bounds are set to fix the
1587  variable and the routine immediately returns (the caller is expected to
1588  reoptimize).
1589
1590  If a branching candidate is found to result in infeasibility in both
1591  directions, the routine immediately returns an indication of infeasibility.
1592
1593  Returns:  0   both branch directions are feasible
1594  -1    branching variable is monotone
1595  -2    infeasible
1596  -3   Use another method
1597
1598  For now just fix on objective from strong branching.
1599*/
1600
1601int CbcNode::chooseDynamicBranch (CbcModel *model, CbcNode *lastNode,
1602                                  OsiSolverBranch * & /*branches*/,
1603                                  int numberPassesLeft)
1604
1605{
1606    if (lastNode)
1607        depth_ = lastNode->depth_ + 1;
1608    else
1609        depth_ = 0;
1610    // Go to other choose if hot start
1611    if (model->hotstartSolution() &&
1612            (((model->moreSpecialOptions()&1024) == 0) || false))
1613        return -3;
1614    delete branch_;
1615    branch_ = NULL;
1616    OsiSolverInterface * solver = model->solver();
1617    // get information on solver type
1618    const OsiAuxInfo * auxInfo = solver->getAuxiliaryInfo();
1619    const OsiBabSolver * auxiliaryInfo = dynamic_cast<const OsiBabSolver *> (auxInfo);
1620    if (!auxiliaryInfo) {
1621        // use one from CbcModel
1622        auxiliaryInfo = model->solverCharacteristics();
1623    }
1624    int numberObjects = model->numberObjects();
1625    // If very odd set of objects then use older chooseBranch
1626    bool useOldWay = false;
1627    // point to useful information
1628    OsiBranchingInformation usefulInfo = model->usefulInformation();
1629    if (numberObjects > model->numberIntegers()) {
1630        for (int i = model->numberIntegers(); i < numberObjects; i++) {
1631            OsiObject * object = model->modifiableObject(i);
1632            CbcObject * obj =   dynamic_cast <CbcObject *>(object) ;
1633            if (!obj || !obj->optionalObject()) {
1634                int preferredWay;
1635                double infeasibility = object->infeasibility(&usefulInfo, preferredWay);
1636                if (infeasibility) {
1637                    useOldWay = true;
1638                    break;
1639                }
1640            } else {
1641              obj->initializeForBranching(model);
1642            }
1643        }
1644    }
1645    if ((model->specialOptions()&128) != 0)
1646        useOldWay = false; // allow
1647    // For now return if not simple
1648    if (useOldWay)
1649        return -3;
1650    // Modify useful info
1651    usefulInfo.depth_ = depth_;
1652    if ((model->specialOptions()&128) != 0) {
1653        // SOS - shadow prices
1654        int numberRows = solver->getNumRows();
1655        const double * pi = usefulInfo.pi_;
1656        double sumPi = 0.0;
1657        for (int i = 0; i < numberRows; i++)
1658            sumPi += fabs(pi[i]);
1659        sumPi /= static_cast<double> (numberRows);
1660        // and scale back
1661        sumPi *= 0.01;
1662        usefulInfo.defaultDual_ = sumPi; // switch on
1663        int numberColumns = solver->getNumCols();
1664        int size = CoinMax(numberColumns, 2 * numberRows);
1665        usefulInfo.usefulRegion_ = new double [size];
1666        CoinZeroN(usefulInfo.usefulRegion_, size);
1667        usefulInfo.indexRegion_ = new int [size];
1668        // pi may change
1669        usefulInfo.pi_ = CoinCopyOfArray(usefulInfo.pi_, numberRows);
1670    }
1671    assert (auxiliaryInfo);
1672    double cutoff = model->getCutoff();
1673    const double * lower = solver->getColLower();
1674    const double * upper = solver->getColUpper();
1675    // See if user thinks infeasible
1676    int anyAction = model->problemFeasibility()->feasible(model, 0);
1677    if (anyAction) {
1678        // will return -2 if infeasible , 0 if treat as integer
1679        return anyAction - 1;
1680    }
1681    int i;
1682    int saveStateOfSearch = model->stateOfSearch() % 10;
1683    int numberStrong = model->numberStrong();
1684    /* Ranging is switched off.
1685       The idea is that you can find out the effect of one iteration
1686       on each unsatisfied variable cheaply.  Then use this
1687       if you have not got much else to go on.
1688    */
1689    //#define RANGING
1690#ifdef RANGING
1691    // must have clp
1692#ifndef COIN_HAS_CLP
1693#  warning("Ranging switched off as not Clp");
1694#undef RANGING
1695#endif
1696    // Pass number
1697    int kPass = 0;
1698    int numberRows = solver->getNumRows();
1699#endif
1700    int numberColumns = model->getNumCols();
1701    double * saveUpper = new double[numberColumns];
1702    double * saveLower = new double[numberColumns];
1703    for (i = 0; i < numberColumns; i++) {
1704        saveLower[i] = lower[i];
1705        saveUpper[i] = upper[i];
1706    }
1707
1708    // Save solution in case heuristics need good solution later
1709
1710    double * saveSolution = new double[numberColumns];
1711    memcpy(saveSolution, solver->getColSolution(), numberColumns*sizeof(double));
1712    model->reserveCurrentSolution(saveSolution);
1713    const double * hotstartSolution = model->hotstartSolution();
1714    const int * hotstartPriorities = model->hotstartPriorities();
1715    double integerTolerance =
1716        model->getDblParam(CbcModel::CbcIntegerTolerance);
1717    if (hotstartSolution) {
1718        if ((model->moreSpecialOptions()&1024) != 0) {
1719            int nBad = 0;
1720            int nUnsat = 0;
1721            int nDiff = 0;
1722            for (int i = 0; i < numberObjects; i++) {
1723                OsiObject * object = model->modifiableObject(i);
1724                const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
1725                if (thisOne) {
1726                    int iColumn = thisOne->columnNumber();
1727                    double targetValue = hotstartSolution[iColumn];
1728                    double value = saveSolution[iColumn];
1729                    if (fabs(value - floor(value + 0.5)) > 1.0e-6) {
1730                        nUnsat++;
1731#ifdef CLP_INVESTIGATE
1732                        printf("H %d is %g target %g\n", iColumn, value, targetValue);
1733#endif
1734                    } else if (fabs(targetValue - value) > 1.0e-6) {
1735                        nDiff++;
1736                    }
1737                    if (targetValue < saveLower[iColumn] ||
1738                            targetValue > saveUpper[iColumn]) {
1739#ifdef CLP_INVESTIGATE
1740                        printf("%d has target %g and current bounds %g and %g\n",
1741                               iColumn, targetValue, saveLower[iColumn], saveUpper[iColumn]);
1742#endif
1743                        nBad++;
1744                    }
1745                }
1746            }
1747#ifdef CLP_INVESTIGATE
1748            printf("Hot %d unsatisfied, %d outside limits, %d different\n",
1749                   nUnsat, nBad, nDiff);
1750#endif
1751            if (nBad) {
1752                // switch off as not possible
1753                hotstartSolution = NULL;
1754                model->setHotstartSolution(NULL, NULL);
1755                usefulInfo.hotstartSolution_ = NULL;
1756            }
1757        }
1758    }
1759    /*
1760      Get a branching decision object. Use the default dynamic decision criteria unless
1761      the user has loaded a decision method into the model.
1762    */
1763    CbcBranchDecision *decision = model->branchingMethod();
1764    if (!decision)
1765        decision = new CbcBranchDynamicDecision();
1766    int xMark = 0;
1767    // Get arrays to sort
1768    double * sort = new double[numberObjects];
1769    int * whichObject = new int[numberObjects];
1770#ifdef RANGING
1771    int xPen = 0;
1772    int * objectMark = new int[2*numberObjects+1];
1773#endif
1774    // Arrays with movements
1775    double * upEstimate = new double[numberObjects];
1776    double * downEstimate = new double[numberObjects];
1777    double estimatedDegradation = 0.0;
1778    int numberNodes = model->getNodeCount();
1779    int saveLogLevel = model->logLevel();
1780#ifdef JJF_ZERO
1781    if ((numberNodes % 500) == 0) {
1782        model->setLogLevel(6);
1783        // Get average up and down costs
1784        double averageUp = 0.0;
1785        double averageDown = 0.0;
1786        int numberUp = 0;
1787        int numberDown = 0;
1788        int i;
1789        for ( i = 0; i < numberObjects; i++) {
1790            OsiObject * object = model->modifiableObject(i);
1791            CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
1792                dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
1793            assert(dynamicObject);
1794            int  numberUp2 = 0;
1795            int numberDown2 = 0;
1796            double up = 0.0;
1797            double down = 0.0;
1798            if (dynamicObject->numberTimesUp()) {
1799                numberUp++;
1800                averageUp += dynamicObject->upDynamicPseudoCost();
1801                numberUp2 += dynamicObject->numberTimesUp();
1802                up = dynamicObject->upDynamicPseudoCost();
1803            }
1804            if (dynamicObject->numberTimesDown()) {
1805                numberDown++;
1806                averageDown += dynamicObject->downDynamicPseudoCost();
1807                numberDown2 += dynamicObject->numberTimesDown();
1808                down = dynamicObject->downDynamicPseudoCost();
1809            }
1810            if (numberUp2 || numberDown2)
1811                printf("col %d - up %d times cost %g, - down %d times cost %g\n",
1812                       dynamicObject->columnNumber(), numberUp2, up, numberDown2, down);
1813        }
1814        if (numberUp)
1815            averageUp /= static_cast<double> (numberUp);
1816        else
1817            averageUp = 1.0;
1818        if (numberDown)
1819            averageDown /= static_cast<double> (numberDown);
1820        else
1821            averageDown = 1.0;
1822        printf("total - up %d vars average %g, - down %d vars average %g\n",
1823               numberUp, averageUp, numberDown, averageDown);
1824    }
1825#endif
1826    int numberBeforeTrust = model->numberBeforeTrust();
1827    // May go round twice if strong branching fixes all local candidates
1828    bool finished = false;
1829    int numberToFix = 0;
1830# ifdef COIN_HAS_CLP
1831    OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
1832    int saveClpOptions = 0;
1833    if (osiclp) {
1834        // for faster hot start
1835        saveClpOptions = osiclp->specialOptions();
1836        osiclp->setSpecialOptions(saveClpOptions | 8192);
1837    }
1838# else
1839    OsiSolverInterface *osiclp = NULL ;
1840# endif
1841    //const CglTreeProbingInfo * probingInfo = NULL; //model->probingInfo();
1842    // Old code left in with DEPRECATED_STRATEGY
1843    assert (model->searchStrategy() == -1 ||
1844            model->searchStrategy() == 1 ||
1845            model->searchStrategy() == 2);
1846#ifdef DEPRECATED_STRATEGY
1847    int saveSearchStrategy2 = model->searchStrategy();
1848#endif
1849    // Get average up and down costs
1850    {
1851        double averageUp = 0.0;
1852        double averageDown = 0.0;
1853        int numberUp = 0;
1854        int numberDown = 0;
1855        int i;
1856        for ( i = 0; i < numberObjects; i++) {
1857            OsiObject * object = model->modifiableObject(i);
1858            CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
1859                dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
1860            if (dynamicObject) {
1861                if (dynamicObject->numberTimesUp()) {
1862                    numberUp++;
1863                    averageUp += dynamicObject->upDynamicPseudoCost();
1864                }
1865                if (dynamicObject->numberTimesDown()) {
1866                    numberDown++;
1867                    averageDown += dynamicObject->downDynamicPseudoCost();
1868                }
1869            }
1870        }
1871        if (numberUp)
1872            averageUp /= static_cast<double> (numberUp);
1873        else
1874            averageUp = 1.0;
1875        if (numberDown)
1876            averageDown /= static_cast<double> (numberDown);
1877        else
1878            averageDown = 1.0;
1879        for ( i = 0; i < numberObjects; i++) {
1880            OsiObject * object = model->modifiableObject(i);
1881            CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
1882                dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
1883            if (dynamicObject) {
1884                if (!dynamicObject->numberTimesUp())
1885                    dynamicObject->setUpDynamicPseudoCost(averageUp);
1886                if (!dynamicObject->numberTimesDown())
1887                    dynamicObject->setDownDynamicPseudoCost(averageDown);
1888            }
1889        }
1890    }
1891    /*
1892      1 strong
1893      2 no strong
1894      3 strong just before solution
1895      4 no strong just before solution
1896      5 strong first time or before solution
1897      6 strong first time
1898    */
1899    int useShadow = model->moreSpecialOptions() & 7;
1900    if (useShadow > 2) {
1901        if (model->getSolutionCount()) {
1902            if (numberNodes || useShadow < 5) {
1903                useShadow = 0;
1904                // zap pseudo shadow prices
1905                model->pseudoShadow(-1);
1906                // and switch off
1907                model->setMoreSpecialOptions(model->moreSpecialOptions()&(~1023));
1908            } else {
1909                useShadow = 1;
1910            }
1911        } else if (useShadow < 5) {
1912            useShadow -= 2;
1913        } else {
1914            useShadow = 1;
1915        }
1916    }
1917    if (useShadow) {
1918        // pseudo shadow prices
1919        model->pseudoShadow((model->moreSpecialOptions() >> 3)&63);
1920    }
1921#ifdef DEPRECATED_STRATEGY
1922    { // in for tabbing
1923    } else if (saveSearchStrategy2 < 1999) {
1924        // pseudo shadow prices
1925        model->pseudoShadow(NULL, NULL);
1926    } else if (saveSearchStrategy2 < 2999) {
1927        // leave old ones
1928    } else if (saveSearchStrategy2 < 3999) {
1929        // pseudo shadow prices at root
1930        if (!numberNodes)
1931            model->pseudoShadow(NULL, NULL);
1932    } else {
1933        abort();
1934    }
1935    if (saveSearchStrategy2 >= 0)
1936        saveSearchStrategy2 = saveSearchStrategy2 % 1000;
1937    if (saveSearchStrategy2 == 999)
1938        saveSearchStrategy2 = -1;
1939    int saveSearchStrategy = saveSearchStrategy2 < 99 ? saveSearchStrategy2 : saveSearchStrategy2 - 100;
1940#endif //DEPRECATED_STRATEGY
1941    int numberNotTrusted = 0;
1942    int numberStrongDone = 0;
1943    int numberUnfinished = 0;
1944    int numberStrongInfeasible = 0;
1945    int numberStrongIterations = 0;
1946    int strongType=0;
1947#define DO_ALL_AT_ROOT
1948#ifdef DO_ALL_AT_ROOT
1949    int saveSatisfiedVariables=0;
1950    int saveNumberToDo=0;
1951#endif
1952    // so we can save lots of stuff
1953    CbcStrongInfo choice;
1954    CbcDynamicPseudoCostBranchingObject * choiceObject = NULL;
1955    if (model->allDynamic()) {
1956        CbcSimpleIntegerDynamicPseudoCost * object = NULL;
1957        choiceObject = new CbcDynamicPseudoCostBranchingObject(model, 0, -1, 0.5, object);
1958    }
1959    choice.possibleBranch = choiceObject;
1960    numberPassesLeft = CoinMax(numberPassesLeft, 2);
1961#ifdef COIN_HAS_NTY
1962    // 1 after, 2 strong, 3 subset
1963    int orbitOption = (model->moreSpecialOptions2()&(128|256))>>7;
1964#endif
1965    //#define DEBUG_SOLUTION
1966#ifdef DEBUG_SOLUTION
1967    bool onOptimalPath=false;
1968    if ((model->specialOptions()&1) != 0) {
1969      const OsiRowCutDebugger *debugger = model->continuousSolver()->getRowCutDebugger() ;
1970      if (debugger) {
1971        const OsiRowCutDebugger *debugger2 = model->solver()->getRowCutDebugger() ;
1972        printf("On optimal in CbcNode %s\n",debugger2 ? "" : "but bad cuts");
1973        onOptimalPath=true;
1974      }
1975    }
1976#endif
1977    while (!finished) {
1978        numberPassesLeft--;
1979        finished = true;
1980        decision->initialize(model);
1981        // Some objects may compute an estimate of best solution from here
1982        estimatedDegradation = 0.0;
1983        numberToFix = 0;
1984        int numberToDo = 0;
1985        int iBestNot = -1;
1986        int iBestGot = -1;
1987        double best = 0.0;
1988        numberNotTrusted = 0;
1989        numberStrongDone = 0;
1990        numberUnfinished = 0;
1991        numberStrongInfeasible = 0;
1992        numberStrongIterations = 0;
1993#ifdef RANGING
1994        int * which = objectMark + numberObjects + 1;
1995        int neededPenalties;
1996        int optionalPenalties;
1997#endif
1998        // We may go round this loop three times (only if we think we have solution)
1999        for (int iPass = 0; iPass < 3; iPass++) {
2000
2001            // Some objects may compute an estimate of best solution from here
2002            estimatedDegradation = 0.0;
2003            numberUnsatisfied_ = 0;
2004            // initialize sum of "infeasibilities"
2005            sumInfeasibilities_ = 0.0;
2006            int bestPriority = COIN_INT_MAX;
2007#ifdef JJF_ZERO
2008            int number01 = 0;
2009            const cliqueEntry * entry = NULL;
2010            const int * toZero = NULL;
2011            const int * toOne = NULL;
2012            const int * backward = NULL;
2013            int numberUnsatisProbed = 0;
2014            int numberUnsatisNotProbed = 0; // 0-1
2015            if (probingInfo) {
2016                number01 = probingInfo->numberIntegers();
2017                entry = probingInfo->fixEntries();
2018                toZero = probingInfo->toZero();
2019                toOne = probingInfo->toOne();
2020                backward = probingInfo->backward();
2021                if (!toZero[number01] || number01 < numberObjects || true) {
2022                    // no info
2023                    probingInfo = NULL;
2024                }
2025            }
2026#endif
2027            /*
2028              Scan for branching objects that indicate infeasibility. Choose candidates
2029              using priority as the first criteria, then integer infeasibility.
2030
2031              The algorithm is to fill the array with a set of good candidates (by
2032              infeasibility) with priority bestPriority.  Finding a candidate with
2033              priority better (less) than bestPriority flushes the choice array. (This
2034              serves as initialization when the first candidate is found.)
2035
2036            */
2037            numberToDo = 0;
2038#ifdef RANGING
2039            neededPenalties = 0;
2040            optionalPenalties = numberObjects;
2041#endif
2042            iBestNot = -1;
2043            double bestNot = 0.0;
2044            iBestGot = -1;
2045            best = 0.0;
2046            /* Problem type as set by user or found by analysis.  This will be extended
2047            0 - not known
2048            1 - Set partitioning <=
2049            2 - Set partitioning ==
2050            3 - Set covering
2051            4 - all +- 1 or all +1 and odd
2052            */
2053            int problemType = model->problemType();
2054            bool canDoOneHot = false;
2055            for (i = 0; i < numberObjects; i++) {
2056                OsiObject * object = model->modifiableObject(i);
2057                CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2058                    dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2059                int preferredWay;
2060                double infeasibility = object->infeasibility(&usefulInfo, preferredWay);
2061                int priorityLevel = object->priority();
2062                if (hotstartSolution) {
2063                    // we are doing hot start
2064                    const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
2065                    if (thisOne) {
2066                        int iColumn = thisOne->columnNumber();
2067                        bool canDoThisHot = true;
2068                        double targetValue = hotstartSolution[iColumn];
2069                        if (saveUpper[iColumn] > saveLower[iColumn]) {
2070                            double value = saveSolution[iColumn];
2071                            if (hotstartPriorities)
2072                                priorityLevel = hotstartPriorities[iColumn];
2073                            //double originalLower = thisOne->originalLower();
2074                            //double originalUpper = thisOne->originalUpper();
2075                            // switch off if not possible
2076                            if (targetValue >= saveLower[iColumn] && targetValue <= saveUpper[iColumn]) {
2077                                /* priority outranks rest always if negative
2078                                   otherwise can be downgraded if at correct level.
2079                                   Infeasibility may be increased to choose 1.0 values first.
2080                                   choose one near wanted value
2081                                */
2082                                if (fabs(value - targetValue) > integerTolerance) {
2083                                    //if (infeasibility>0.01)
2084                                    //infeasibility = fabs(1.0e6-fabs(value-targetValue));
2085                                    //else
2086                                    infeasibility = fabs(value - targetValue);
2087                                    //if (targetValue==1.0)
2088                                    //infeasibility += 1.0;
2089                                    if (value > targetValue) {
2090                                        preferredWay = -1;
2091                                    } else {
2092                                        preferredWay = 1;
2093                                    }
2094                                    priorityLevel = CoinAbs(priorityLevel);
2095                                } else if (priorityLevel < 0) {
2096                                    priorityLevel = CoinAbs(priorityLevel);
2097                                    if (targetValue == saveLower[iColumn]) {
2098                                        infeasibility = integerTolerance + 1.0e-12;
2099                                        preferredWay = -1;
2100                                    } else if (targetValue == saveUpper[iColumn]) {
2101                                        infeasibility = integerTolerance + 1.0e-12;
2102                                        preferredWay = 1;
2103                                    } else {
2104                                        // can't
2105                                        priorityLevel += 10000000;
2106                                        canDoThisHot = false;
2107                                    }
2108                                } else {
2109                                    priorityLevel += 10000000;
2110                                    canDoThisHot = false;
2111                                }
2112                            } else {
2113                                // switch off if not possible
2114                                canDoThisHot = false;
2115                            }
2116                            if (canDoThisHot)
2117                                canDoOneHot = true;
2118                        } else if (targetValue < saveLower[iColumn] || targetValue > saveUpper[iColumn]) {
2119                        }
2120                    } else {
2121                        priorityLevel += 10000000;
2122                    }
2123                }
2124#define ZERO_ONE 0
2125#define ZERO_FAKE 1.0e20;
2126#if ZERO_ONE==1
2127                // branch on 0-1 first (temp)
2128                if (fabs(saveSolution[dynamicObject->columnNumber()]) < 1.0)
2129                    priorityLevel--;
2130#endif
2131#if ZERO_ONE==2
2132                if (fabs(saveSolution[dynamicObject->columnNumber()]) < 1.0)
2133                    infeasibility *= ZERO_FAKE;
2134#endif
2135                if (infeasibility) {
2136                    int iColumn = numberColumns + i;
2137                    bool gotDown = false;
2138                    int numberThisDown = 0;
2139                    bool gotUp = false;
2140                    int numberThisUp = 0;
2141                    double downGuess = object->downEstimate();
2142                    double upGuess = object->upEstimate();
2143                    if (dynamicObject) {
2144                        // Use this object's numberBeforeTrust
2145                        int numberBeforeTrustThis = dynamicObject->numberBeforeTrust();
2146                        iColumn = dynamicObject->columnNumber();
2147                        gotDown = false;
2148                        numberThisDown = dynamicObject->numberTimesDown();
2149                        if (numberThisDown >= numberBeforeTrustThis)
2150                            gotDown = true;
2151                        gotUp = false;
2152                        numberThisUp = dynamicObject->numberTimesUp();
2153                        if (numberThisUp >= numberBeforeTrustThis)
2154                            gotUp = true;
2155                        if (!depth_ && false) {
2156                            // try closest to 0.5
2157                            double part = saveSolution[iColumn] - floor(saveSolution[iColumn]);
2158                            infeasibility = fabs(0.5 - part);
2159                        }
2160                        if (problemType > 0 && problemType < 4 && false) {
2161                            // try closest to 0.5
2162                            double part = saveSolution[iColumn] - floor(saveSolution[iColumn]);
2163                            infeasibility = 0.5 - fabs(0.5 - part);
2164                        }
2165#ifdef JJF_ZERO
2166                        if (probingInfo) {
2167                            int iSeq = backward[iColumn];
2168                            assert (iSeq >= 0);
2169                            infeasibility = 1.0 + (toZero[iSeq+1] - toZero[iSeq]) +
2170                                            5.0 * CoinMin(toOne[iSeq] - toZero[iSeq], toZero[iSeq+1] - toOne[iSeq]);
2171                            if (toZero[iSeq+1] > toZero[iSeq]) {
2172                                numberUnsatisProbed++;
2173                            } else {
2174                                numberUnsatisNotProbed++;
2175                            }
2176                        }
2177#endif
2178                    } else {
2179                        // see if SOS
2180                        CbcSOS * sosObject =
2181                            dynamic_cast <CbcSOS *>(object) ;
2182                        if (sosObject) {
2183                            gotDown = false;
2184                            numberThisDown = sosObject->numberTimesDown();
2185                            if (numberThisDown >= numberBeforeTrust)
2186                                gotDown = true;
2187                            gotUp = false;
2188                            numberThisUp = sosObject->numberTimesUp();
2189                            if (numberThisUp >= numberBeforeTrust)
2190                                gotUp = true;
2191                        } else {
2192                            gotDown = true;
2193                            numberThisDown = 999999;
2194                            downGuess = 1.0e20;
2195                            gotUp = true;
2196                            numberThisUp = 999999;
2197                            upGuess = 1.0e20;
2198                            numberPassesLeft = 0;
2199                        }
2200                    }
2201                    // Increase estimated degradation to solution
2202                    estimatedDegradation += CoinMin(downGuess, upGuess);
2203                    downEstimate[i] = downGuess;
2204                    upEstimate[i] = upGuess;
2205                    numberUnsatisfied_++;
2206                    sumInfeasibilities_ += infeasibility;
2207                    // Better priority? Flush choices.
2208                    if (priorityLevel < bestPriority) {
2209                        numberToDo = 0;
2210                        bestPriority = priorityLevel;
2211                        iBestGot = -1;
2212                        best = 0.0;
2213                        numberNotTrusted = 0;
2214#ifdef RANGING
2215                        neededPenalties = 0;
2216                        optionalPenalties = numberObjects;
2217#endif
2218                    } else if (priorityLevel > bestPriority) {
2219                        continue;
2220                    }
2221                    if (!gotUp || !gotDown)
2222                        numberNotTrusted++;
2223                    // Check for suitability based on infeasibility.
2224                    if ((gotDown && gotUp) && numberStrong > 0) {
2225                        sort[numberToDo] = -infeasibility;
2226                        if (infeasibility > best) {
2227                            best = infeasibility;
2228                            iBestGot = numberToDo;
2229                        }
2230#ifdef RANGING
2231                        if (dynamicObject) {
2232                            objectMark[--optionalPenalties] = numberToDo;
2233                            which[optionalPenalties] = iColumn;
2234                        }
2235#endif
2236                    } else {
2237#ifdef RANGING
2238                        if (dynamicObject) {
2239                            objectMark[neededPenalties] = numberToDo;
2240                            which[neededPenalties++] = iColumn;
2241                        }
2242#endif
2243                        sort[numberToDo] = -10.0 * infeasibility;
2244                        if (!(numberThisUp + numberThisDown))
2245                            sort[numberToDo] *= 100.0; // make even more likely
2246                        if (iColumn < numberColumns) {
2247                            double part = saveSolution[iColumn] - floor(saveSolution[iColumn]);
2248                            if (1.0 - fabs(part - 0.5) > bestNot) {
2249                                iBestNot = numberToDo;
2250                                bestNot = 1.0 - fabs(part - 0.5);
2251                            }
2252                        } else {
2253                            // SOS
2254                            if (-sort[numberToDo] > bestNot) {
2255                                iBestNot = numberToDo;
2256                                bestNot = -sort[numberToDo];
2257                            }
2258                        }
2259                    }
2260                    if (model->messageHandler()->logLevel() > 3) {
2261                        printf("%d (%d) down %d %g up %d %g - infeas %g - sort %g solution %g\n",
2262                               i, iColumn, numberThisDown, object->downEstimate(), numberThisUp, object->upEstimate(),
2263                               infeasibility, sort[numberToDo], saveSolution[iColumn]);
2264                    }
2265                    whichObject[numberToDo++] = i;
2266                } else {
2267                    // for debug
2268                    downEstimate[i] = -1.0;
2269                    upEstimate[i] = -1.0;
2270                }
2271            }
2272            if (numberUnsatisfied_) {
2273                //if (probingInfo&&false)
2274                //printf("nunsat %d, %d probed, %d other 0-1\n",numberUnsatisfied_,
2275                // numberUnsatisProbed,numberUnsatisNotProbed);
2276                // some infeasibilities - go to next steps
2277                if (!canDoOneHot && hotstartSolution) {
2278                    // switch off as not possible
2279                    hotstartSolution = NULL;
2280                    model->setHotstartSolution(NULL, NULL);
2281                    usefulInfo.hotstartSolution_ = NULL;
2282                }
2283                break;
2284            } else if (!iPass) {
2285                // may just need resolve
2286                model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
2287                double newObjValue = solver->getObjSense()*solver->getObjValue();
2288                objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
2289                if (!solver->isProvenOptimal()) {
2290                    // infeasible
2291                    anyAction = -2;
2292                    break;
2293                }
2294                // Double check looks OK - just look at rows with all integers
2295                if (model->allDynamic()) {
2296                    double * solution = CoinCopyOfArray(saveSolution, numberColumns);
2297                    for (int i = 0; i < numberColumns; i++) {
2298                        if (model->isInteger(i))
2299                            solution[i] = floor(solution[i] + 0.5);
2300                    }
2301                    int numberRows = solver->getNumRows();
2302                    double * rowActivity = new double [numberRows];
2303                    CoinZeroN(rowActivity, numberRows);
2304                    solver->getMatrixByCol()->times(solution, rowActivity);
2305                    //const double * element = model->solver()->getMatrixByCol()->getElements();
2306                    const int * row = model->solver()->getMatrixByCol()->getIndices();
2307                    const CoinBigIndex * columnStart = model->solver()->getMatrixByCol()->getVectorStarts();
2308                    const int * columnLength = model->solver()->getMatrixByCol()->getVectorLengths();
2309                    int nFree = 0;
2310                    int nFreeNon = 0;
2311                    int nFixedNon = 0;
2312                    double mostAway = 0.0;
2313                    int whichAway = -1;
2314                    const double * columnLower = solver->getColLower();
2315                    const double * columnUpper = solver->getColUpper();
2316                    for (int i = 0; i < numberColumns; i++) {
2317                        if (!model->isInteger(i)) {
2318                            // mark rows as flexible
2319                            CoinBigIndex start = columnStart[i];
2320                            CoinBigIndex end = start + columnLength[i];
2321                            for (CoinBigIndex j = start; j < end; j++) {
2322                                int iRow = row[j];
2323                                rowActivity[iRow] = COIN_DBL_MAX;
2324                            }
2325                        } else if (columnLower[i] < columnUpper[i]) {
2326                            if (fabs(solution[i] - saveSolution[i]) > 
2327                                integerTolerance) {
2328                                nFreeNon++;
2329                                if (fabs(solution[i] - saveSolution[i]) > mostAway) {
2330                                    mostAway = fabs(solution[i] - saveSolution[i]);
2331                                    whichAway = i;
2332                                }
2333                            } else {
2334                                nFree++;
2335                            }
2336                        } else if (solution[i] != saveSolution[i]) {
2337                            nFixedNon++;
2338                        }
2339                    }
2340                    const double * lower = solver->getRowLower();
2341                    const double * upper = solver->getRowUpper();
2342                    bool satisfied = true;
2343                    for (int i = 0; i < numberRows; i++) {
2344                        double value = rowActivity[i];
2345                        if (value != COIN_DBL_MAX) {
2346                            if (value > upper[i] + 1.0e-5 || value < lower[i] - 1.0e-5) {
2347                                satisfied = false;
2348                            }
2349                        }
2350                    }
2351                    delete [] rowActivity;
2352                    delete [] solution;
2353                    if (!satisfied) {
2354#ifdef CLP_INVESTIGATE
2355                        printf("%d free ok %d free off target %d fixed off target\n",
2356                               nFree, nFreeNon, nFixedNon);
2357#endif
2358                        if (nFreeNon) {
2359                            // try branching on these
2360                            delete branch_;
2361                            for (int i = 0; i < numberObjects; i++) {
2362                                OsiObject * object = model->modifiableObject(i);
2363                                CbcSimpleIntegerDynamicPseudoCost * obj =
2364                                    dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2365                                assert (obj);
2366                                int iColumn = obj->columnNumber();
2367                                if (iColumn == whichAway) {
2368                                    int preferredWay = (saveSolution[iColumn] > solution[iColumn])
2369                                                       ? -1 : +1;
2370                                    usefulInfo.integerTolerance_ = 0.0;
2371                                    branch_ = obj->createCbcBranch(solver, &usefulInfo, preferredWay);
2372                                    break;
2373                                }
2374                            }
2375                            anyAction = 0;
2376                            break;
2377                        }
2378                    }
2379                }
2380            } else if (iPass == 1) {
2381                // looks like a solution - get paranoid
2382                bool roundAgain = false;
2383                // get basis
2384                CoinWarmStartBasis * ws = dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
2385                if (!ws)
2386                    break;
2387                double tolerance;
2388                solver->getDblParam(OsiPrimalTolerance, tolerance);
2389                for (i = 0; i < numberColumns; i++) {
2390                    double value = saveSolution[i];
2391                    if (value < lower[i] - tolerance) {
2392                        saveSolution[i] = lower[i];
2393                        roundAgain = true;
2394                        ws->setStructStatus(i, CoinWarmStartBasis::atLowerBound);
2395                    } else if (value > upper[i] + tolerance) {
2396                        saveSolution[i] = upper[i];
2397                        roundAgain = true;
2398                        ws->setStructStatus(i, CoinWarmStartBasis::atUpperBound);
2399                    }
2400                }
2401                if (roundAgain) {
2402                    // restore basis
2403                    solver->setWarmStart(ws);
2404                    solver->setColSolution(saveSolution);
2405                    delete ws;
2406                    bool takeHint;
2407                    OsiHintStrength strength;
2408                    solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
2409                    solver->setHintParam(OsiDoDualInResolve, false, OsiHintDo) ;
2410                    model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
2411                    double newObjValue = solver->getObjSense()*solver->getObjValue();
2412                    objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
2413                    solver->setHintParam(OsiDoDualInResolve, takeHint, strength) ;
2414                    if (!solver->isProvenOptimal()) {
2415                        // infeasible
2416                        anyAction = -2;
2417                        break;
2418                    }
2419                } else {
2420                    delete ws;
2421                    break;
2422                }
2423            }
2424        }
2425        if (anyAction == -2) {
2426            break;
2427        }
2428        // skip if solution
2429        if (!numberUnsatisfied_)
2430            break;
2431        int skipAll = (numberNotTrusted == 0 || numberToDo == 1) ? 1 : 0;
2432        bool doneHotStart = false;
2433        //DEPRECATED_STRATEGYint searchStrategy = saveSearchStrategy>=0 ? (saveSearchStrategy%10) : -1;
2434        int searchStrategy = model->searchStrategy();
2435        // But adjust depending on ratio of iterations
2436        if (searchStrategy > 0) {
2437          if (numberBeforeTrust >= /*5*/ 10 && numberBeforeTrust <= 10) {
2438                if (searchStrategy != 2) {
2439                    assert (searchStrategy == 1);
2440                    if (depth_ > 5) {
2441                        int numberIterations = model->getIterationCount();
2442                        int numberStrongIterations = model->numberStrongIterations();
2443                        if (numberStrongIterations > numberIterations + 10000) {
2444                            searchStrategy = 2;
2445                            skipAll = 1;
2446                        } else if (numberStrongIterations*4 + 1000 < numberIterations) {
2447                            searchStrategy = 3;
2448                            skipAll = 0;
2449                        }
2450                    } else {
2451                        searchStrategy = 3;
2452                        skipAll = 0;
2453                    }
2454                }
2455            }
2456        }
2457        // worth trying if too many times
2458        // Save basis
2459        CoinWarmStart * ws = NULL;
2460        // save limit
2461        int saveLimit = 0;
2462        solver->getIntParam(OsiMaxNumIterationHotStart, saveLimit);
2463        if (!numberPassesLeft)
2464            skipAll = 1;
2465        if (!skipAll) {
2466            ws = solver->getWarmStart();
2467            int limit = 100;
2468            if (!saveStateOfSearch && saveLimit < limit && saveLimit == 100)
2469                solver->setIntParam(OsiMaxNumIterationHotStart, limit);
2470        }
2471        // Say which one will be best
2472        int whichChoice = 0;
2473        int bestChoice;
2474        if (iBestGot >= 0)
2475            bestChoice = iBestGot;
2476        else
2477            bestChoice = iBestNot;
2478        assert (bestChoice >= 0);
2479        // If we have hit max time don't do strong branching
2480        bool hitMaxTime = (model->getCurrentSeconds() >
2481                            model->getDblParam(CbcModel::CbcMaximumSeconds));
2482        // also give up if we are looping round too much
2483        if (hitMaxTime || numberPassesLeft <= 0 || useShadow == 2) {
2484            int iObject = whichObject[bestChoice];
2485            OsiObject * object = model->modifiableObject(iObject);
2486            int preferredWay;
2487            object->infeasibility(&usefulInfo, preferredWay);
2488            CbcObject * obj =
2489                dynamic_cast <CbcObject *>(object) ;
2490            assert (obj);
2491            branch_ = obj->createCbcBranch(solver, &usefulInfo, preferredWay);
2492            {
2493                CbcBranchingObject * branchObj =
2494                    dynamic_cast <CbcBranchingObject *>(branch_) ;
2495                assert (branchObj);
2496                branchObj->way(preferredWay);
2497            }
2498            delete ws;
2499            ws = NULL;
2500            break;
2501        } else {
2502            // say do fast
2503            int easy = 1;
2504            if (!skipAll)
2505                solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, &easy) ;
2506            int iDo;
2507#define RESET_BOUNDS
2508#ifdef RANGING
2509            bool useRanging = model->allDynamic() && !skipAll;
2510            if (useRanging) {
2511                double currentObjective = solver->getObjValue() * solver->getObjSense();
2512                double gap = cutoff - currentObjective;
2513                // relax a bit
2514                gap *= 1.0000001;
2515                gap = CoinMax(1.0e-5, gap);
2516                // off penalties if too much
2517                double needed = neededPenalties;
2518                needed *= numberRows;
2519                if (numberNodes) {
2520                    if (needed > 1.0e6) {
2521                        neededPenalties = 0;
2522                    } else if (gap < 1.0e5) {
2523                        // maybe allow some not needed
2524                        int extra = static_cast<int> ((1.0e6 - needed) / numberRows);
2525                        int nStored = numberObjects - optionalPenalties;
2526                        extra = CoinMin(extra, nStored);
2527                        for (int i = 0; i < extra; i++) {
2528                            objectMark[neededPenalties] = objectMark[optionalPenalties+i];
2529                            which[neededPenalties++] = which[optionalPenalties+i];;
2530                        }
2531                    }
2532                }
2533                if (osiclp && neededPenalties) {
2534                    assert (!doneHotStart);
2535                    xPen += neededPenalties;
2536                    which--;
2537                    which[0] = neededPenalties;
2538                    osiclp->passInRanges(which);
2539                    // Mark hot start and get ranges
2540                    if (kPass) {
2541                        // until can work out why solution can go funny
2542                        int save = osiclp->specialOptions();
2543                        osiclp->setSpecialOptions(save | 256);
2544                        solver->markHotStart();
2545#ifdef RESET_BOUNDS
2546                        memcpy(saveLower,solver->getColLower(),solver->getNumCols()*sizeof(double));
2547                        memcpy(saveUpper,solver->getColUpper(),solver->getNumCols()*sizeof(double));
2548#endif
2549                        osiclp->setSpecialOptions(save);
2550                    } else {
2551                        solver->markHotStart();
2552#ifdef RESET_BOUNDS
2553                        memcpy(saveLower,solver->getColLower(),solver->getNumCols()*sizeof(double));
2554                        memcpy(saveUpper,solver->getColUpper(),solver->getNumCols()*sizeof(double));
2555#endif
2556                    }
2557                    doneHotStart = true;
2558                    xMark++;
2559                    kPass++;
2560                    osiclp->passInRanges(NULL);
2561                    const double * downCost = osiclp->upRange();
2562                    const double * upCost = osiclp->downRange();
2563                    bool problemFeasible = true;
2564                    int numberFixed = 0;
2565                    for (int i = 0; i < neededPenalties; i++) {
2566                        int j = objectMark[i];
2567                        int iObject = whichObject[j];
2568                        OsiObject * object = model->modifiableObject(iObject);
2569                        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2570                            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2571                        // Use this object's numberBeforeTrust
2572                        int numberBeforeTrustThis = dynamicObject->numberBeforeTrust();
2573                        int iSequence = dynamicObject->columnNumber();
2574                        double value = saveSolution[iSequence];
2575                        value -= floor(value);
2576                        double upPenalty = CoinMin(upCost[i], 1.0e110) * (1.0 - value);
2577                        double downPenalty = CoinMin(downCost[i], 1.0e110) * value;
2578                        int numberThisDown = dynamicObject->numberTimesDown();
2579                        int numberThisUp = dynamicObject->numberTimesUp();
2580                        if (!numberBeforeTrustThis) {
2581                            // override
2582                            downEstimate[iObject] = downPenalty;
2583                            upEstimate[iObject] = upPenalty;
2584                            double min1 = CoinMin(downEstimate[iObject],
2585                                                  upEstimate[iObject]);
2586                            double max1 = CoinMax(downEstimate[iObject],
2587                                                  upEstimate[iObject]);
2588                            min1 = 0.8 * min1 + 0.2 * max1;
2589                            sort[j] = - min1;
2590                        } else if (numberThisDown < numberBeforeTrustThis ||
2591                                   numberThisUp < numberBeforeTrustThis) {
2592                            double invTrust = 1.0 / static_cast<double> (numberBeforeTrustThis);
2593                            if (numberThisDown < numberBeforeTrustThis) {
2594                                double fraction = numberThisDown * invTrust;
2595                                downEstimate[iObject] = fraction * downEstimate[iObject] + (1.0 - fraction) * downPenalty;
2596                            }
2597                            if (numberThisUp < numberBeforeTrustThis) {
2598                                double fraction = numberThisUp * invTrust;
2599                                upEstimate[iObject] = fraction * upEstimate[iObject] + (1.0 - fraction) * upPenalty;
2600                            }
2601                            double min1 = CoinMin(downEstimate[iObject],
2602                                                  upEstimate[iObject]);
2603                            double max1 = CoinMax(downEstimate[iObject],
2604                                                  upEstimate[iObject]);
2605                            min1 = 0.8 * min1 + 0.2 * max1;
2606                            min1 *= 10.0;
2607                            if (!(numberThisDown + numberThisUp))
2608                                min1 *= 100.0;
2609                            sort[j] = - min1;
2610                        }
2611                        // seems unreliable
2612                        if (false&&CoinMax(downPenalty, upPenalty) > gap) {
2613                            COIN_DETAIL_PRINT(printf("gap %g object %d has down range %g, up %g\n",
2614                                                     gap, i, downPenalty, upPenalty));
2615                            //sort[j] -= 1.0e50; // make more likely to be chosen
2616                            int number;
2617                            if (downPenalty > gap) {
2618                                number = dynamicObject->numberTimesDown();
2619                                if (upPenalty > gap)
2620                                    problemFeasible = false;
2621                                CbcBranchingObject * branch = dynamicObject->createCbcBranch(solver, &usefulInfo, 1);
2622                                //branch->fix(solver,saveLower,saveUpper,1);
2623                                delete branch;
2624                            } else {
2625                                number = dynamicObject->numberTimesUp();
2626                                CbcBranchingObject * branch = dynamicObject->createCbcBranch(solver, &usefulInfo, 1);
2627                                //branch->fix(solver,saveLower,saveUpper,-1);
2628                                delete branch;
2629                            }
2630                            if (number >= numberBeforeTrustThis)
2631                              dynamicObject->setNumberBeforeTrust(CoinMin(number + 1,5*numberBeforeTrust));
2632                            numberFixed++;
2633                        }
2634                        if (!numberNodes)
2635                            COIN_DETAIL_PRINT(printf("%d pen down ps %g -> %g up ps %g -> %g\n",
2636                                                     iObject, downPenalty, downPenalty, upPenalty, upPenalty));
2637                    }
2638                    if (numberFixed && problemFeasible) {
2639                        assert(doneHotStart);
2640                        solver->unmarkHotStart();
2641                        model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
2642                        double newObjValue = solver->getObjSense()*solver->getObjValue();
2643                        objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
2644                        solver->markHotStart();
2645#ifdef RESET_BOUNDS
2646                        memcpy(saveLower,solver->getColLower(),solver->getNumCols()*sizeof(double));
2647                        memcpy(saveUpper,solver->getColUpper(),solver->getNumCols()*sizeof(double));
2648#endif
2649                        problemFeasible = solver->isProvenOptimal();
2650                    }
2651                    if (!problemFeasible) {
2652                      COIN_DETAIL_PRINT(fprintf(stdout, "both ways infeas on ranging - code needed\n"));
2653                        anyAction = -2;
2654                        if (!choiceObject) {
2655                            delete choice.possibleBranch;
2656                            choice.possibleBranch = NULL;
2657                        }
2658                        //printf("Both infeasible for choice %d sequence %d\n",i,
2659                        // model->object(choice.objectNumber)->columnNumber());
2660                        // Delete the snapshot
2661                        solver->unmarkHotStart();
2662                        // back to normal
2663                        solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
2664                        // restore basis
2665                        solver->setWarmStart(ws);
2666                        doneHotStart = false;
2667                        delete ws;
2668                        ws = NULL;
2669                        break;
2670                    }
2671                }
2672            }
2673#endif          /* RANGING */
2674            {
2675                int numberIterations = model->getIterationCount();
2676                //numberIterations += (model->numberExtraIterations()>>2);
2677                const int * strongInfo = model->strongInfo();
2678                //int numberDone = strongInfo[0]-strongInfo[3];
2679                int numberFixed = strongInfo[1] - strongInfo[4];
2680                int numberInfeasible = strongInfo[2] - strongInfo[5];
2681                assert (!strongInfo[3]);
2682                assert (!strongInfo[4]);
2683                assert (!strongInfo[5]);
2684                int numberStrongIterations = model->numberStrongIterations();
2685                int numberRows = solver->getNumRows();
2686                if (numberStrongIterations > numberIterations + CoinMin(100, 10*numberRows) && depth_ >= 4 && numberNodes > 100) {
2687                    if (20*numberInfeasible + 4*numberFixed < numberNodes) {
2688                        // Say never do
2689                        if (numberBeforeTrust == 5)
2690                          skipAll = -1;
2691                    }
2692                }
2693            }
2694            // make sure best will be first
2695            if (iBestGot >= 0)
2696                sort[iBestGot] = -COIN_DBL_MAX;
2697            // Actions 0 - exit for repeat, 1 resolve and try old choice,2 exit for continue
2698            if (anyAction)
2699                numberToDo = 0; // skip as we will be trying again
2700            // Sort
2701            CoinSort_2(sort, sort + numberToDo, whichObject);
2702            // Change in objective opposite infeasible
2703            double worstFeasible = 0.0;
2704            // Just first if strong off
2705            if (!numberStrong)
2706                numberToDo = CoinMin(numberToDo, 1);
2707            if (searchStrategy == 2)
2708                numberToDo = CoinMin(numberToDo, 20);
2709            iDo = 0;
2710            int saveLimit2;
2711            solver->getIntParam(OsiMaxNumIterationHotStart, saveLimit2);
2712            int numberTest = numberNotTrusted > 0 ? numberStrong : (numberStrong + 1) / 2;
2713            if (searchStrategy == 3) {
2714                // Previously decided we need strong
2715                numberTest = numberStrong;
2716            }
2717            // Try nearly always off
2718            if (skipAll >= 0) {
2719                if (searchStrategy < 2) {
2720                    //if ((numberNodes%20)!=0) {
2721                    if ((model->specialOptions()&8) == 0) {
2722                        numberTest = 0;
2723                    }
2724                    //} else {
2725                    //numberTest=2*numberStrong;
2726                    //skipAll=0;
2727                    //}
2728                }
2729            } else {
2730                // Just take first
2731                numberTest = 1;
2732            }
2733            int testDepth = (skipAll >= 0) ? 8 : 4;
2734            if (depth_ < testDepth && numberStrong) {
2735                if (searchStrategy != 2) {
2736                    int numberRows = solver->getNumRows();
2737                    // whether to do this or not is important - think
2738                    if (numberRows < 300 || numberRows + numberColumns < 2500) {
2739                        if (depth_ < 7)
2740                            numberStrong = CoinMin(3 * numberStrong, numberToDo);
2741                        if (!depth_)
2742                            numberStrong = CoinMin(6 * numberStrong, numberToDo);
2743                    }
2744                    numberTest = numberStrong;
2745                    skipAll = 0;
2746                }
2747            }
2748            // Do at least 5 strong
2749            if (numberColumns < 1000 && (depth_ < 15 || numberNodes < 1000000))
2750                numberTest = CoinMax(numberTest, 5);
2751            if ((model->specialOptions()&8) == 0) {
2752                if (skipAll) {
2753                    numberTest = 0;
2754                }
2755            } else {
2756                // do 5 as strong is fixing
2757                numberTest = CoinMax(numberTest, 5);
2758            }
2759            // see if switched off
2760            if (skipAll < 0) {
2761                numberTest = 0;
2762            }
2763            int realMaxHotIterations = 999999;
2764            if (skipAll < 0)
2765                numberToDo = 1;
2766            strongType=0;
2767#ifdef DO_ALL_AT_ROOT
2768            if (model->strongStrategy()) {
2769              int iStrategy=model->strongStrategy();
2770              int kDepth = iStrategy/100;
2771              if (kDepth)
2772                iStrategy -= 100*kDepth;
2773              else
2774                kDepth=5;
2775              double objValue = solver->getObjSense()*solver->getObjValue();
2776              double bestPossible = model->getBestPossibleObjValue();
2777              bestPossible += 1.0e-7*(1.0+fabs(bestPossible));
2778              int jStrategy = iStrategy/10;
2779              if (jStrategy) {
2780                if ((jStrategy&1)!=0&&!depth_)
2781                  strongType=2;
2782                else if ((jStrategy&2)!=0&&depth_<=kDepth)
2783                  strongType=2;
2784                else if ((jStrategy&4)!=0&&objValue<bestPossible)
2785                  strongType=2;
2786                iStrategy-=10*jStrategy;
2787              }
2788              if (!strongType) {
2789                if ((iStrategy&1)!=0&&!depth_)
2790                  strongType=1;
2791                else if ((iStrategy&2)!=0&&depth_<=kDepth)
2792                  strongType=1;
2793                else if ((iStrategy&4)!=0&&objValue<bestPossible)
2794                  strongType=1;
2795              }
2796              saveNumberToDo=numberToDo;
2797              if (strongType==2) {
2798                // add in satisfied
2799                const int * integerVariable = model->integerVariable();
2800                int numberIntegers = model->numberIntegers();
2801                if (numberIntegers==numberObjects) {
2802                  numberToDo=0;
2803                  for (int i=0;i<numberIntegers;i++) {
2804                    int iColumn=integerVariable[i];
2805                    if (saveUpper[iColumn]>saveLower[iColumn]) {
2806                      whichObject [numberToDo++]=i;
2807                    }
2808                  }
2809                  saveSatisfiedVariables=numberToDo-saveNumberToDo;
2810                } else {
2811                  strongType=1;
2812                }
2813              }
2814              if (strongType) {
2815                numberTest = numberToDo;
2816                numberStrong=numberToDo;
2817                skipAll=0;
2818                searchStrategy=0;
2819                solver->setIntParam(OsiMaxNumIterationHotStart, 100000);
2820                //printf("Strong branching type %d\n",strongType);
2821              }
2822            }
2823#endif
2824#ifdef COIN_HAS_NTY
2825            const int * orbits = NULL;
2826#endif
2827#ifdef COIN_HAS_NTY
2828            if (orbitOption>1) {
2829              CbcSymmetry * symmetryInfo = model->symmetryInfo();
2830              CbcNodeInfo * infoX = lastNode ? lastNode->nodeInfo() : NULL;
2831              bool worthTrying = false;
2832              if (infoX) {
2833                CbcNodeInfo * info = infoX;
2834                for (int i=0;i<NTY_BAD_DEPTH;i++) {
2835                  if (!info->parent()) {
2836                    worthTrying = true;
2837                    break;
2838                  }
2839                  info = info->parent();
2840                  if (info->symmetryWorked()) {
2841                    worthTrying = true;
2842                    break;
2843                  }
2844                }
2845              } else {
2846                worthTrying=true;
2847              }
2848              if (symmetryInfo && worthTrying) {
2849                symmetryInfo->ChangeBounds(solver->getColLower(),
2850                                            solver->getColUpper(),
2851                                            solver->getNumCols(),false);
2852                symmetryInfo->Compute_Symmetry();
2853                symmetryInfo->fillOrbits();
2854                orbits = symmetryInfo->whichOrbit();
2855                int iColumn=-1;
2856                if (orbits && symmetryInfo->numberUsefulObjects()) {
2857                  bool doBranch=true;
2858                  int numberUsefulOrbits = symmetryInfo->numberUsefulObjects();
2859                  if (numberUsefulOrbits<2) {
2860                    assert (numberUsefulOrbits);
2861                    double largest=-1.0;
2862                    for (int i=0;i<numberColumns;i++) {
2863                      if (orbits[i]>=0) {
2864                        if (saveSolution[i]>largest) {
2865                          largest=saveSolution[i];
2866                          iColumn=i;
2867                        }
2868                      }
2869                    }
2870                  } else {
2871#if COIN_HAS_NTY2 == 1
2872                    // take largest
2873                    int iOrbit=symmetryInfo->largestOrbit(solver->getColLower(),
2874                                                          solver->getColUpper());
2875                    double largest=-1.0;
2876                    for (int i=0;i<numberColumns;i++) {
2877                      if (orbits[i]==iOrbit) {
2878                        if (saveSolution[i]>largest) {
2879                          largest=saveSolution[i];
2880                          iColumn=i;
2881                        }
2882                      }
2883                    }
2884#endif
2885                    if (orbitOption==2) {
2886                      // strong
2887                      int nDo=0;
2888                      const double * lower = solver->getColLower();
2889                      const double * upper = solver->getColUpper();
2890                      for (int iOrbit = 0; iOrbit < numberUsefulOrbits; iOrbit++) {
2891                        double distance=1.0;
2892                        int iColumn = -1;
2893                        for (int i=0;i<numberColumns;i++) {
2894                          if (orbits[i]==iOrbit &&lower[i]==0.0&&upper[i]==1.0) {
2895                            double away = fabs(saveSolution[i]-0.5);
2896                            if (away<distance&&away<0.4999) {
2897                              distance=away;
2898                              iColumn=i;
2899                            }
2900                          }
2901                        }
2902                        if (iColumn>=0) 
2903                          whichObject[nDo++]=iColumn;
2904                      }
2905                      if (nDo)
2906                        numberToDo=nDo;
2907                      doBranch=false;
2908                    } else if (orbitOption==3) {
2909                      // subset
2910                      int nDo=0;
2911                      for (int iDo = 0; iDo < numberToDo; iDo++) {
2912                        int iObject = whichObject[iDo];
2913                        OsiObject * object = model->modifiableObject(iObject);
2914                        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2915                          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2916                        int iColumn = dynamicObject ? dynamicObject->columnNumber() : -1;
2917                        if (iColumn<0||orbits[iColumn]>=0)
2918                          whichObject[nDo++]=whichObject[iDo];
2919                      }
2920                      assert(nDo);
2921                      //printf("nDo %d\n",nDo);
2922                      numberToDo=nDo;
2923                      doBranch=false;
2924                      /* need NULL as if two in same orbit and strong branching fixes
2925                         then we may be in trouble.
2926                         Strong option should be OK as only one in set done.
2927                       */ 
2928                      orbits=NULL;
2929                    }
2930                  }
2931                  if(doBranch) {
2932                    orbitOption=0;
2933                    branch_ = new CbcOrbitalBranchingObject(model,iColumn,1,0,NULL);
2934                    if (infoX)
2935                      infoX->setSymmetryWorked();
2936                    numberToDo=0;
2937                  }
2938                }
2939              }
2940            }
2941#endif
2942            for ( iDo = 0; iDo < numberToDo; iDo++) {
2943                int iObject = whichObject[iDo];
2944                OsiObject * object = model->modifiableObject(iObject);
2945                CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2946                    dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2947                int iColumn = dynamicObject ? dynamicObject->columnNumber() : numberColumns + iObject;
2948                int preferredWay;
2949                double infeasibility = object->infeasibility(&usefulInfo, preferredWay);
2950                bool feasibleSolution=false;
2951                double predictedChange=0.0;
2952                // may have become feasible
2953                if (!infeasibility) {
2954                  if(strongType!=2||solver->getColLower()[iColumn]==solver->getColUpper()[iColumn])
2955                    continue;
2956                }
2957#ifndef NDEBUG
2958                if (iColumn < numberColumns) {
2959                    const double * solution = model->testSolution();
2960                    assert (saveSolution[iColumn] == solution[iColumn]);
2961                }
2962#endif
2963                CbcSimpleInteger * obj =
2964                    dynamic_cast <CbcSimpleInteger *>(object) ;
2965                if (obj) {
2966                    if (choiceObject) {
2967                        obj->fillCreateBranch(choiceObject, &usefulInfo, preferredWay);
2968                        choiceObject->setObject(dynamicObject);
2969                    } else {
2970                        choice.possibleBranch = obj->createCbcBranch(solver, &usefulInfo, preferredWay);
2971                    }
2972                } else {
2973                    CbcObject * obj =
2974                        dynamic_cast <CbcObject *>(object) ;
2975                    assert (obj);
2976                    choice.possibleBranch = obj->createCbcBranch(solver, &usefulInfo, preferredWay);
2977                }
2978                // Save which object it was
2979                choice.objectNumber = iObject;
2980                choice.numIntInfeasUp = numberUnsatisfied_;
2981                choice.numIntInfeasDown = numberUnsatisfied_;
2982                if (strongType!=2) {
2983                  choice.upMovement = upEstimate[iObject];
2984                  choice.downMovement = downEstimate[iObject];
2985                } else {
2986                  choice.upMovement = 0.1;
2987                  choice.downMovement = 0.1;
2988                }
2989                assert (choice.upMovement >= 0.0);
2990                assert (choice.downMovement >= 0.0);
2991                choice.fix = 0; // say not fixed
2992                // see if can skip strong branching
2993                int canSkip = choice.possibleBranch->fillStrongInfo(choice);
2994                if ((numberTest <= 0 || skipAll)) {
2995                    if (iDo > 20) {
2996                        if (!choiceObject) {
2997                            delete choice.possibleBranch;
2998                            choice.possibleBranch = NULL;
2999                        }
3000                        break; // give up anyway
3001                    }
3002                }
3003                if (model->messageHandler()->logLevel() > 3 && numberBeforeTrust && dynamicObject)
3004                    dynamicObject->print(1, choice.possibleBranch->value());
3005                if (strongType)
3006                  canSkip=0;
3007                if (skipAll < 0)
3008                    canSkip = 1;
3009                if (!canSkip) {
3010                    if (!doneHotStart) {
3011                        // Mark hot start
3012                        doneHotStart = true;
3013                        solver->markHotStart();
3014#ifdef RESET_BOUNDS
3015                        memcpy(saveLower,solver->getColLower(),solver->getNumCols()*sizeof(double));
3016                        memcpy(saveUpper,solver->getColUpper(),solver->getNumCols()*sizeof(double));
3017#endif
3018                        if (!solver->isProvenOptimal()) {
3019                          skipAll=-2;
3020                          canSkip = 1;
3021                        }
3022                        xMark++;
3023                    }
3024                }
3025                if (!canSkip) {
3026                    numberTest--;
3027                    // just do a few
3028                    if (searchStrategy == 2)
3029                        solver->setIntParam(OsiMaxNumIterationHotStart, 10);
3030                    double objectiveChange ;
3031                    double newObjectiveValue = 1.0e100;
3032                    int j;
3033                    // status is 0 finished, 1 infeasible and other
3034                    int iStatus;
3035                    /*
3036                      Try the down direction first. (Specify the initial branching alternative as
3037                      down with a call to way(-1). Each subsequent call to branch() performs the
3038                      specified branch and advances the branch object state to the next branch
3039                      alternative.)
3040                    */
3041                    choice.possibleBranch->way(-1) ;
3042                    predictedChange = choice.possibleBranch->branch() ;
3043#ifdef COIN_HAS_NTY
3044                    if (orbits) {
3045                      // can fix all in orbit
3046                      int fixOrbit = orbits[iObject];
3047                      if (fixOrbit>=0) {
3048                        //printf("fixing all in orbit %d for column %d\n",fixOrbit,iObject);
3049                        for (int i=0;i<numberColumns;i++) {
3050                          if (orbits[i]==fixOrbit)
3051                            solver->setColUpper(i,0.0);
3052                        }
3053                      }
3054                    }
3055#endif
3056                    solver->solveFromHotStart() ;
3057                    bool needHotStartUpdate = false;
3058                    numberStrongDone++;
3059                    numberStrongIterations += solver->getIterationCount();
3060                    /*
3061                      We now have an estimate of objective degradation that we can use for strong
3062                      branching. If we're over the cutoff, the variable is monotone up.
3063                      If we actually made it to optimality, check for a solution, and if we have
3064                      a good one, call setBestSolution to process it. Note that this may reduce the
3065                      cutoff, so we check again to see if we can declare this variable monotone.
3066                    */
3067                    if (solver->isProvenOptimal())
3068                        iStatus = 0; // optimal
3069                    else if (solver->isIterationLimitReached() 
3070                             && !solver->isDualObjectiveLimitReached()) {
3071                        iStatus = 2; // unknown
3072                    } else {
3073                        iStatus = 1; // infeasible
3074#ifdef CONFLICT_CUTS
3075# ifdef COIN_HAS_CLP
3076                        if (osiclp&&(model->moreSpecialOptions()&4194304)!=0) {
3077                          const CbcFullNodeInfo * topOfTree =
3078                            model->topOfTree();
3079                          if (topOfTree) {
3080                            OsiRowCut * cut = osiclp->smallModelCut(topOfTree->lower(),
3081                                                                    topOfTree->upper(),
3082                                                                    model->numberRowsAtContinuous(),
3083                                                                    model->whichGenerator());
3084                            if (cut) {
3085                              printf("XXXXXX found conflict cut in strong branching\n");
3086                              //cut->print();
3087                              if ((model->specialOptions()&1) != 0) {
3088                                const OsiRowCutDebugger *debugger = model->continuousSolver()->getRowCutDebugger() ;
3089                                if (debugger) {
3090                                  if (debugger->invalidCut(*cut)) {
3091                                    model->continuousSolver()->applyRowCuts(1,cut);
3092                                    model->continuousSolver()->writeMps("bad");
3093                                  }
3094                                  CoinAssert (!debugger->invalidCut(*cut));
3095                                }
3096                              }
3097                              model->makeGlobalCut(cut) ;
3098                            }
3099                          }
3100                        }
3101#endif
3102#endif
3103                    }
3104                    // say infeasible if branch says so
3105                    if (predictedChange==COIN_DBL_MAX)
3106                      iStatus=1;
3107                    if (iStatus != 2 && solver->getIterationCount() >
3108                            realMaxHotIterations)
3109                        numberUnfinished++;
3110                    newObjectiveValue = solver->getObjSense() * solver->getObjValue();
3111                    choice.numItersDown = solver->getIterationCount();
3112                    objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_, 0.0);
3113                    // Update branching information if wanted
3114                    CbcBranchingObject * cbcobj = dynamic_cast<CbcBranchingObject *> (choice.possibleBranch);
3115                    if (cbcobj) {
3116                        CbcObject * object = cbcobj->object();
3117                        assert (object) ;
3118                        CbcObjectUpdateData update = object->createUpdateInformation(solver, this, cbcobj);
3119                        update.objectNumber_ = choice.objectNumber;
3120                        model->addUpdateInformation(update);
3121                    } else {
3122                        decision->updateInformation( solver, this);
3123                    }
3124                    if (!iStatus) {
3125                        choice.finishedDown = true ;
3126                        if (newObjectiveValue >= cutoff) {
3127                            objectiveChange = 1.0e100; // say infeasible
3128                            numberStrongInfeasible++;
3129                        } else {
3130#define CBCNODE_TIGHTEN_BOUNDS
3131#ifdef CBCNODE_TIGHTEN_BOUNDS
3132                            // Can we tighten bounds?
3133                            if (iColumn < numberColumns && cutoff < 1.0e20
3134                                    && objectiveChange > 1.0e-5) {
3135                                double value = saveSolution[iColumn];
3136                                double down = value - floor(value-integerTolerance);
3137                                double changePer = objectiveChange / (down + 1.0e-7);
3138                                double distance = (cutoff - objectiveValue_) / changePer;
3139                                distance += 1.0e-3;
3140                                if (distance < 5.0) {
3141                                    double newLower = ceil(value - distance);
3142                                    if (newLower > saveLower[iColumn]) {
3143                                        //printf("Could increase lower bound on %d from %g to %g\n",
3144                                        //   iColumn,saveLower[iColumn],newLower);
3145                                        saveLower[iColumn] = newLower;
3146                                        solver->setColLower(iColumn, newLower);
3147                                    }
3148                                }
3149                            }
3150#endif
3151                            // See if integer solution
3152                            feasibleSolution = 
3153                              model->feasibleSolution(choice.numIntInfeasDown,
3154                                                      choice.numObjInfeasDown);
3155                            if (feasibleSolution
3156                                    && model->problemFeasibility()->feasible(model, -1) >= 0) {
3157                                if (auxiliaryInfo->solutionAddsCuts()) {
3158                                    needHotStartUpdate = true;
3159                                    solver->unmarkHotStart();
3160                                }
3161                                model->setLogLevel(saveLogLevel);
3162                                model->setBestSolution(CBC_STRONGSOL,
3163                                                       newObjectiveValue,
3164                                                       solver->getColSolution()) ;
3165                                if (needHotStartUpdate) {
3166                                    model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3167                                    newObjectiveValue = solver->getObjSense() * solver->getObjValue();
3168                                    objectiveValue_ = CoinMax(objectiveValue_,newObjectiveValue);
3169                                    objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_, 0.0);
3170                                    model->feasibleSolution(choice.numIntInfeasDown,
3171                                                            choice.numObjInfeasDown);
3172                                }
3173                                model->setLastHeuristic(NULL);
3174                                model->incrementUsed(solver->getColSolution());
3175                                cutoff = model->getCutoff();
3176                                if (newObjectiveValue >= cutoff) {      //  *new* cutoff
3177                                    objectiveChange = 1.0e100 ;
3178                                    numberStrongInfeasible++;
3179                                }
3180                            }
3181                        }
3182                    } else if (iStatus == 1) {
3183                        objectiveChange = 1.0e100 ;
3184                        numberStrongInfeasible++;
3185                    } else {
3186                        // Can't say much as we did not finish
3187                        choice.finishedDown = false ;
3188                        numberUnfinished++;
3189                    }
3190                    choice.downMovement = objectiveChange ;
3191
3192                    // restore bounds
3193                    for ( j = 0; j < numberColumns; j++) {
3194                        if (saveLower[j] != lower[j])
3195                            solver->setColLower(j, saveLower[j]);
3196                        if (saveUpper[j] != upper[j])
3197                            solver->setColUpper(j, saveUpper[j]);
3198                    }
3199                    if (needHotStartUpdate) {
3200                        needHotStartUpdate = false;
3201                        model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3202                        double newObjValue = solver->getObjSense()*solver->getObjValue();
3203                        objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3204                        //we may again have an integer feasible solution
3205                        int numberIntegerInfeasibilities;
3206                        int numberObjectInfeasibilities;
3207                        if (model->feasibleSolution(
3208                                    numberIntegerInfeasibilities,
3209                                    numberObjectInfeasibilities)) {
3210#ifdef BONMIN
3211                            //In this case node has become integer feasible, let us exit the loop
3212                            std::cout << "Node has become integer feasible" << std::endl;
3213                            numberUnsatisfied_ = 0;
3214                            break;
3215#endif
3216                            double objValue = solver->getObjValue();
3217                            model->setLogLevel(saveLogLevel);
3218                            model->setBestSolution(CBC_STRONGSOL,
3219                                                   objValue,
3220                                                   solver->getColSolution()) ;
3221                            model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3222                            double newObjValue = solver->getObjSense()*solver->getObjValue();
3223                            objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3224                            cutoff = model->getCutoff();
3225                        }
3226                        solver->markHotStart();
3227#ifdef RESET_BOUNDS
3228                        memcpy(saveLower,solver->getColLower(),solver->getNumCols()*sizeof(double));
3229                        memcpy(saveUpper,solver->getColUpper(),solver->getNumCols()*sizeof(double));
3230#endif
3231                        if (!solver->isProvenOptimal()) {
3232                          skipAll=-2;
3233                          canSkip = 1;
3234                        }
3235                        xMark++;
3236                    }
3237#if 0 //def DO_ALL_AT_ROOT
3238                    if (strongType)
3239                        printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
3240                               choice.objectNumber, iStatus, newObjectiveValue, choice.numItersDown,
3241                               choice.downMovement, choice.finishedDown, choice.numIntInfeasDown,
3242                               choice.numObjInfeasDown);
3243#endif
3244
3245                    // repeat the whole exercise, forcing the variable up
3246                    predictedChange=choice.possibleBranch->branch();
3247                    solver->solveFromHotStart() ;
3248                    numberStrongDone++;
3249                    numberStrongIterations += solver->getIterationCount();
3250                    /*
3251                      We now have an estimate of objective degradation that we can use for strong
3252                      branching. If we're over the cutoff, the variable is monotone up.
3253                      If we actually made it to optimality, check for a solution, and if we have
3254                      a good one, call setBestSolution to process it. Note that this may reduce the
3255                      cutoff, so we check again to see if we can declare this variable monotone.
3256                    */
3257                    if (solver->isProvenOptimal())
3258                        iStatus = 0; // optimal
3259                    else if (solver->isIterationLimitReached()
3260                             && !solver->isDualObjectiveLimitReached()) {
3261                        iStatus = 2; // unknown
3262                    } else {
3263                        iStatus = 1; // infeasible
3264#ifdef CONFLICT_CUTS
3265# ifdef COIN_HAS_CLP
3266                        if (osiclp&&(model->moreSpecialOptions()&4194304)!=0) {
3267                          const CbcFullNodeInfo * topOfTree =
3268                            model->topOfTree();
3269                          if (topOfTree) {
3270                            OsiRowCut * cut = osiclp->smallModelCut(topOfTree->lower(),
3271                                                                    topOfTree->upper(),
3272                                                                    model->numberRowsAtContinuous(),
3273                                                                    model->whichGenerator());
3274                            if (cut) {
3275                              printf("XXXXXX found conflict cut in strong branching\n");
3276                              //cut->print();
3277                              if ((model->specialOptions()&1) != 0) {
3278                                const OsiRowCutDebugger *debugger = model->continuousSolver()->getRowCutDebugger() ;
3279                                if (debugger) {
3280                                  if (debugger->invalidCut(*cut)) {
3281                                    model->continuousSolver()->applyRowCuts(1,cut);
3282                                    model->continuousSolver()->writeMps("bad");
3283                                  }
3284                                  CoinAssert (!debugger->invalidCut(*cut));
3285                                }
3286                              }
3287                              model->makeGlobalCut(cut) ;
3288                            }
3289                          }
3290                        }
3291#endif
3292#endif
3293                    }
3294                    // say infeasible if branch says so
3295                    if (predictedChange==COIN_DBL_MAX)
3296                      iStatus=1;
3297                    if (iStatus != 2 && solver->getIterationCount() >
3298                            realMaxHotIterations)
3299                        numberUnfinished++;
3300                    newObjectiveValue = solver->getObjSense() * solver->getObjValue();
3301                    choice.numItersUp = solver->getIterationCount();
3302                    objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_, 0.0);
3303                    // Update branching information if wanted
3304                    cbcobj = dynamic_cast<CbcBranchingObject *> (choice.possibleBranch);
3305                    if (cbcobj) {
3306                        CbcObject * object = cbcobj->object();
3307                        assert (object) ;
3308                        CbcObjectUpdateData update = object->createUpdateInformation(solver, this, cbcobj);
3309                        update.objectNumber_ = choice.objectNumber;
3310                        model->addUpdateInformation(update);
3311                    } else {
3312                        decision->updateInformation( solver, this);
3313                    }
3314                    if (!iStatus) {
3315                        choice.finishedUp = true ;
3316                        if (newObjectiveValue >= cutoff) {
3317                            objectiveChange = 1.0e100; // say infeasible
3318                            numberStrongInfeasible++;
3319                        } else {
3320#ifdef CBCNODE_TIGHTEN_BOUNDS
3321                            // Can we tighten bounds?
3322                            if (iColumn < numberColumns && cutoff < 1.0e20
3323                                    && objectiveChange > 1.0e-5) {
3324                                double value = saveSolution[iColumn];
3325                                double up = ceil(value+integerTolerance) - value;
3326                                double changePer = objectiveChange / (up + 1.0e-7);
3327                                double distance = (cutoff - objectiveValue_) / changePer;
3328                                distance += 1.0e-3;
3329                                if (distance < 5.0) {
3330                                    double newUpper = floor(value + distance);
3331                                    if (newUpper < saveUpper[iColumn]) {
3332                                        //printf("Could decrease upper bound on %d from %g to %g\n",
3333                                        //   iColumn,saveUpper[iColumn],newUpper);
3334                                        saveUpper[iColumn] = newUpper;
3335                                        solver->setColUpper(iColumn, newUpper);
3336                                    }
3337                                }
3338                            }
3339#endif
3340                            // See if integer solution
3341                            feasibleSolution = 
3342                              model->feasibleSolution(choice.numIntInfeasUp,
3343                                                      choice.numObjInfeasUp);
3344                            if (feasibleSolution
3345                                    && model->problemFeasibility()->feasible(model, -1) >= 0) {
3346#ifdef BONMIN
3347                                std::cout << "Node has become integer feasible" << std::endl;
3348                                numberUnsatisfied_ = 0;
3349                                break;
3350#endif
3351                                if (auxiliaryInfo->solutionAddsCuts()) {
3352                                    needHotStartUpdate = true;
3353                                    solver->unmarkHotStart();
3354                                }
3355                                model->setLogLevel(saveLogLevel);
3356                                model->setBestSolution(CBC_STRONGSOL,
3357                                                       newObjectiveValue,
3358                                                       solver->getColSolution()) ;
3359                                if (choice.finishedDown) {
3360                                  double cutoff = model->getCutoff();
3361                                  double downObj = objectiveValue_
3362                                    + choice.downMovement ;
3363                                  if (downObj >= cutoff) {     
3364                                    choice.downMovement = 1.0e100 ;
3365                                    numberStrongInfeasible++;
3366                                }
3367
3368                                }
3369                                if (needHotStartUpdate) {
3370                                    model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3371                                    newObjectiveValue = solver->getObjSense() * solver->getObjValue();
3372                                    objectiveValue_ = CoinMax(objectiveValue_,newObjectiveValue);
3373                                    objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_, 0.0);
3374                                    model->feasibleSolution(choice.numIntInfeasDown,
3375                                                            choice.numObjInfeasDown);
3376                                }
3377                                model->setLastHeuristic(NULL);
3378                                model->incrementUsed(solver->getColSolution());
3379                                cutoff = model->getCutoff();
3380                                if (newObjectiveValue >= cutoff) {      //  *new* cutoff
3381                                    objectiveChange = 1.0e100 ;
3382                                    numberStrongInfeasible++;
3383                                }
3384                            }
3385                        }
3386                    } else if (iStatus == 1) {
3387                        objectiveChange = 1.0e100 ;
3388                        numberStrongInfeasible++;
3389                    } else {
3390                        // Can't say much as we did not finish
3391                        choice.finishedUp = false ;
3392                        numberUnfinished++;
3393                    }
3394                    choice.upMovement = objectiveChange ;
3395
3396                    // restore bounds
3397                    for ( j = 0; j < numberColumns; j++) {
3398                        if (saveLower[j] != lower[j])
3399                            solver->setColLower(j, saveLower[j]);
3400                        if (saveUpper[j] != upper[j])
3401                            solver->setColUpper(j, saveUpper[j]);
3402                    }
3403                    if (needHotStartUpdate) {
3404                        needHotStartUpdate = false;
3405                        model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3406                        double newObjValue = solver->getObjSense()*solver->getObjValue();
3407                        objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3408                        //we may again have an integer feasible solution
3409                        int numberIntegerInfeasibilities;
3410                        int numberObjectInfeasibilities;
3411                        if (model->feasibleSolution(
3412                                    numberIntegerInfeasibilities,
3413                                    numberObjectInfeasibilities)) {
3414                            double objValue = solver->getObjValue();
3415                            model->setLogLevel(saveLogLevel);
3416                            model->setBestSolution(CBC_STRONGSOL,
3417                                                   objValue,
3418                                                   solver->getColSolution()) ;
3419                            model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3420                            double newObjValue = solver->getObjSense()*solver->getObjValue();
3421                            objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3422                            cutoff = model->getCutoff();
3423                        }
3424                        solver->markHotStart();
3425#ifdef RESET_BOUNDS
3426                        memcpy(saveLower,solver->getColLower(),solver->getNumCols()*sizeof(double));
3427                        memcpy(saveUpper,solver->getColUpper(),solver->getNumCols()*sizeof(double));
3428#endif
3429                        if (!solver->isProvenOptimal()) {
3430                          skipAll=-2;
3431                          canSkip = 1;
3432                        }
3433                        xMark++;
3434                    }
3435
3436#if 0 //def DO_ALL_AT_ROOT
3437                    if (strongType)
3438                        printf("Up on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
3439                               choice.objectNumber, iStatus, newObjectiveValue, choice.numItersUp,
3440                               choice.upMovement, choice.finishedUp, choice.numIntInfeasUp,
3441                               choice.numObjInfeasUp);
3442#endif
3443                }
3444
3445                solver->setIntParam(OsiMaxNumIterationHotStart, saveLimit2);
3446                /*
3447                  End of evaluation for this candidate variable. Possibilities are:
3448                  * Both sides below cutoff; this variable is a candidate for branching.
3449                  * Both sides infeasible or above the objective cutoff: no further action
3450                  here. Break from the evaluation loop and assume the node will be purged
3451                  by the caller.
3452                  * One side below cutoff: Install the branch (i.e., fix the variable). Break
3453                  from the evaluation loop and assume the node will be reoptimised by the
3454                  caller.
3455                */
3456                // reset
3457                choice.possibleBranch->resetNumberBranchesLeft();
3458                if (choice.upMovement < 1.0e100) {
3459                    if (choice.downMovement < 1.0e100) {
3460                        // In case solution coming in was odd
3461                        choice.upMovement = CoinMax(0.0, choice.upMovement);
3462                        choice.downMovement = CoinMax(0.0, choice.downMovement);
3463#if ZERO_ONE==2
3464                        // branch on 0-1 first (temp)
3465                        if (fabs(choice.possibleBranch->value()) < 1.0) {
3466                            choice.upMovement *= ZERO_FAKE;
3467                            choice.downMovement *= ZERO_FAKE;
3468                        }
3469#endif
3470                        // feasible - see which best
3471                        if (!canSkip) {
3472                            if (model->messageHandler()->logLevel() > 3)
3473                                printf("sort %g downest %g upest %g ", sort[iDo], downEstimate[iObject],
3474                                       upEstimate[iObject]);
3475                            model->messageHandler()->message(CBC_STRONG, *model->messagesPointer())
3476                            << iObject << iColumn
3477                            << choice.downMovement << choice.numIntInfeasDown
3478                            << choice.upMovement << choice.numIntInfeasUp
3479                            << choice.possibleBranch->value()
3480                            << CoinMessageEol;
3481                        }
3482                        int betterWay=0;
3483                        // If was feasible (extra strong branching) skip
3484                        if (infeasibility) {
3485                            CbcBranchingObject * branchObj =
3486                                dynamic_cast <CbcBranchingObject *>(branch_) ;
3487                            if (branch_)
3488                                assert (branchObj);
3489                            betterWay = decision->betterBranch(choice.possibleBranch,
3490                                                               branchObj,
3491                                                               choice.upMovement,
3492                                                               choice.numIntInfeasUp ,
3493                                                               choice.downMovement,
3494                                                               choice.numIntInfeasDown );
3495                        }
3496                        if (betterWay) {
3497                            // C) create branching object
3498                            if (choiceObject) {
3499                                delete branch_;
3500                                branch_ = choice.possibleBranch->clone();
3501                            } else {
3502                                delete branch_;
3503                                branch_ = choice.possibleBranch;
3504                                choice.possibleBranch = NULL;
3505                            }
3506                            {
3507                                CbcBranchingObject * branchObj =
3508                                    dynamic_cast <CbcBranchingObject *>(branch_) ;
3509                                assert (branchObj);
3510                                //branchObj->way(preferredWay);
3511                                branchObj->way(betterWay);
3512                            }
3513                            bestChoice = choice.objectNumber;
3514                            whichChoice = iDo;
3515                            if (numberStrong <= 1) {
3516                                delete ws;
3517                                ws = NULL;
3518                                break;
3519                            }
3520                        } else {
3521                            if (!choiceObject) {
3522                                delete choice.possibleBranch;
3523                                choice.possibleBranch = NULL;
3524                            }
3525                            if (iDo >= 2*numberStrong) {
3526                                delete ws;
3527                                ws = NULL;
3528                                break;
3529                            }
3530                            if (!dynamicObject || dynamicObject->numberTimesUp() > 1) {
3531                                if (iDo - whichChoice >= numberStrong) {
3532                                    if (!choiceObject) {
3533                                        delete choice.possibleBranch;
3534                                        choice.possibleBranch = NULL;
3535                                    }
3536                                    break; // give up
3537                                }
3538                            } else {
3539                                if (iDo - whichChoice >= 2*numberStrong) {
3540                                    delete ws;
3541                                    ws = NULL;
3542                                    if (!choiceObject) {
3543                                        delete choice.possibleBranch;
3544                                        choice.possibleBranch = NULL;
3545                                    }
3546                                    break; // give up
3547                                }
3548                            }
3549                        }
3550                    } else {
3551                        // up feasible, down infeasible
3552                        anyAction = -1;
3553                        worstFeasible = CoinMax(worstFeasible, choice.upMovement);
3554                        model->messageHandler()->message(CBC_STRONG, *model->messagesPointer())
3555                        << iObject << iColumn
3556                        << choice.downMovement << choice.numIntInfeasDown
3557                        << choice.upMovement << choice.numIntInfeasUp
3558                        << choice.possibleBranch->value()
3559                        << CoinMessageEol;
3560                        //printf("Down infeasible for choice %d sequence %d\n",i,
3561                        // model->object(choice.objectNumber)->columnNumber());
3562                        choice.fix = 1;
3563                        numberToFix++;
3564                        choice.possibleBranch->fix(solver, saveLower, saveUpper, 1);
3565                        if (!choiceObject) {
3566                            delete choice.possibleBranch;
3567                            choice.possibleBranch = NULL;
3568                        } else {
3569                            //choiceObject = new CbcDynamicPseudoCostBranchingObject(*choiceObject);
3570                            choice.possibleBranch = choiceObject;
3571                        }
3572                        assert(doneHotStart);
3573                        solver->unmarkHotStart();
3574                        model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3575                        double newObjValue = solver->getObjSense()*solver->getObjValue();
3576                        objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3577                        bool goneInfeasible = (!solver->isProvenOptimal()||solver->isDualObjectiveLimitReached());
3578                        solver->markHotStart();
3579#ifdef RESET_BOUNDS
3580                        memcpy(saveLower,solver->getColLower(),solver->getNumCols()*sizeof(double));
3581                        memcpy(saveUpper,solver->getColUpper(),solver->getNumCols()*sizeof(double));
3582#endif
3583                        if (!solver->isProvenOptimal()) {
3584                          skipAll=-2;
3585                          canSkip = 1;
3586                        }
3587                        xMark++;
3588                        // may be infeasible (if other way stopped on iterations)
3589                        if (goneInfeasible) {
3590                            // neither side feasible
3591                            anyAction = -2;
3592                            if (!choiceObject) {
3593                                delete choice.possibleBranch;
3594                                choice.possibleBranch = NULL;
3595                            }
3596                            //printf("Both infeasible for choice %d sequence %d\n",i,
3597                            // model->object(choice.objectNumber)->columnNumber());
3598                            delete ws;
3599                            ws = NULL;
3600                            break;
3601                        }
3602                    }
3603                } else {
3604                    if (choice.downMovement < 1.0e100) {
3605                        // down feasible, up infeasible
3606                        anyAction = -1;
3607                        worstFeasible = CoinMax(worstFeasible, choice.downMovement);
3608                        model->messageHandler()->message(CBC_STRONG, *model->messagesPointer())
3609                        << iObject << iColumn
3610                        << choice.downMovement << choice.numIntInfeasDown
3611                        << choice.upMovement << choice.numIntInfeasUp
3612                        << choice.possibleBranch->value()
3613                        << CoinMessageEol;
3614                        choice.fix = -1;
3615                        numberToFix++;
3616                        choice.possibleBranch->fix(solver, saveLower, saveUpper, -1);
3617                        if (!choiceObject) {
3618                            delete choice.possibleBranch;
3619                            choice.possibleBranch = NULL;
3620                        } else {
3621                            //choiceObject = new CbcDynamicPseudoCostBranchingObject(*choiceObject);
3622                            choice.possibleBranch = choiceObject;
3623                        }
3624                        assert(doneHotStart);
3625                        solver->unmarkHotStart();
3626                        model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3627                        double newObjValue = solver->getObjSense()*solver->getObjValue();
3628                        objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3629                        bool goneInfeasible = (!solver->isProvenOptimal()||solver->isDualObjectiveLimitReached());
3630                        solver->markHotStart();
3631#ifdef RESET_BOUNDS
3632                        memcpy(saveLower,solver->getColLower(),solver->getNumCols()*sizeof(double));
3633                        memcpy(saveUpper,solver->getColUpper(),solver->getNumCols()*sizeof(double));
3634#endif
3635                        if (!solver->isProvenOptimal()) {
3636                          skipAll=-2;
3637                          canSkip = 1;
3638                        }
3639                        xMark++;
3640                        // may be infeasible (if other way stopped on iterations)
3641                        if (goneInfeasible) {
3642                            // neither side feasible
3643                            anyAction = -2;
3644                            if (!choiceObject) {
3645                                delete choice.possibleBranch;
3646                                choice.possibleBranch = NULL;
3647                            }
3648                            delete ws;
3649                            ws = NULL;
3650                            break;
3651                        }
3652                    } else {
3653                        // neither side feasible
3654                        anyAction = -2;
3655                        if (!choiceObject) {
3656                            delete choice.possibleBranch;
3657                            choice.possibleBranch = NULL;
3658                        }
3659                        delete ws;
3660                        ws = NULL;
3661                        break;
3662                    }
3663                }
3664                // Check max time
3665                hitMaxTime = (model->getCurrentSeconds() >
3666                               model->getDblParam(CbcModel::CbcMaximumSeconds));
3667                if (hitMaxTime) {
3668                    // make sure rest are fast
3669                    for ( int jDo = iDo + 1; jDo < numberToDo; jDo++) {
3670                        int iObject = whichObject[iDo];
3671                        OsiObject * object = model->modifiableObject(iObject);
3672                        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3673                            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3674                        if (dynamicObject)
3675                            dynamicObject->setNumberBeforeTrust(0);
3676                    }
3677                    numberTest = 0;
3678                }
3679                if (!choiceObject) {
3680                    delete choice.possibleBranch;
3681                }
3682            }
3683            if (model->messageHandler()->logLevel() > 3) {
3684                if (anyAction == -2) {
3685                    printf("infeasible\n");
3686                } else if (anyAction == -1) {
3687                    printf("%d fixed AND choosing %d iDo %d iChosenWhen %d numberToDo %d\n", numberToFix, bestChoice,
3688                           iDo, whichChoice, numberToDo);
3689                } else {
3690                    int iObject = whichObject[whichChoice];
3691                    OsiObject * object = model->modifiableObject(iObject);
3692                    CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3693                        dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3694                    if (dynamicObject) {
3695                        int iColumn = dynamicObject->columnNumber();
3696                        printf("choosing %d (column %d) iChosenWhen %d numberToDo %d\n", bestChoice,
3697                               iColumn, whichChoice, numberToDo);
3698                    }
3699                }
3700            }
3701            if (doneHotStart) {
3702                // Delete the snapshot
3703                solver->unmarkHotStart();
3704                // back to normal
3705                solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
3706                // restore basis
3707                solver->setWarmStart(ws);
3708            }
3709            solver->setIntParam(OsiMaxNumIterationHotStart, saveLimit);
3710            // Unless infeasible we will carry on
3711            // But we could fix anyway
3712            if (numberToFix && !hitMaxTime) {
3713                if (anyAction != -2) {
3714                    // apply and take off
3715                    bool feasible = true;
3716                    // can do quick optimality check
3717                    int easy = 2;
3718                    solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, &easy) ;
3719                    model->resolve(NULL, 11, saveSolution, saveLower, saveUpper) ;
3720                    double newObjValue = solver->getObjSense()*solver->getObjValue();
3721                    objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3722                    solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
3723                    feasible = solver->isProvenOptimal();
3724                    if (feasible) {
3725                        anyAction = 0;
3726                    } else {
3727                        anyAction = -2;
3728                        finished = true;
3729                    }
3730                }
3731            }
3732            // If  fixed then round again
3733            // See if candidate still possible
3734            if (branch_) {
3735                 const OsiObject * object = model->object(bestChoice);
3736                 int preferredWay;
3737                 double infeasibility = object->infeasibility(&usefulInfo, preferredWay);
3738                 if (!infeasibility) {
3739                   // take out
3740                   delete branch_;
3741                   branch_ = NULL;
3742                 } else {
3743                   CbcBranchingObject * branchObj =
3744                     dynamic_cast <CbcBranchingObject *>(branch_) ;
3745                   assert (branchObj);
3746                   branchObj->way(preferredWay);
3747#ifdef CBCNODE_TIGHTEN_BOUNDS
3748                   bool fixed = branchObj->tighten(solver);
3749                   if (fixed) {
3750                     //printf("Variable now fixed!\n");
3751                     // take out
3752                     delete branch_;
3753                     branch_ = NULL;
3754                   }
3755#endif
3756                 }
3757            }
3758            if (!branch_ && anyAction != -2 && !hitMaxTime) {
3759                finished = false;
3760            }
3761            delete ws;
3762        }
3763    }
3764    // update number of strong iterations etc
3765    model->incrementStrongInfo(numberStrongDone, numberStrongIterations,
3766                               anyAction == -2 ? 0 : numberToFix, anyAction == -2);
3767    if (model->searchStrategy() == -1) {
3768#ifndef COIN_DEVELOP
3769        if (solver->messageHandler()->logLevel() > 1)
3770#endif
3771            printf("%d strong, %d iters, %d inf, %d not finished, %d not trusted\n",
3772                   numberStrongDone, numberStrongIterations, numberStrongInfeasible, numberUnfinished,
3773                   numberNotTrusted);
3774        // decide what to do
3775        int strategy = 1;
3776        if (((numberUnfinished*4 > numberStrongDone &&
3777                numberStrongInfeasible*40 < numberStrongDone) ||
3778                numberStrongInfeasible < 0) && model->numberStrong() < 10 && model->numberBeforeTrust() <= 20 && model->numberObjects() > CoinMax(1000, solver->getNumRows())) {
3779            strategy = 2;
3780#ifdef COIN_DEVELOP
3781            //if (model->logLevel()>1)
3782            printf("going to strategy 2\n");
3783#endif
3784            // Weaken
3785            model->setNumberStrong(2);
3786            model->setNumberBeforeTrust(1);
3787            model->synchronizeNumberBeforeTrust();
3788        }
3789        if (numberNodes)
3790            strategy = 1;  // should only happen after hot start
3791        model->setSearchStrategy(strategy);
3792    } else if (numberStrongDone) {
3793        //printf("%d strongB, %d iters, %d inf, %d not finished, %d not trusted\n",
3794        //   numberStrongDone,numberStrongIterations,numberStrongInfeasible,numberUnfinished,
3795        //   numberNotTrusted);
3796    }
3797    if (model->searchStrategy() == 1 && numberNodes > 500 && numberNodes < -510) {
3798#ifndef COIN_DEVELOP
3799        if (solver->messageHandler()->logLevel() > 1)
3800#endif
3801            printf("after %d nodes - %d strong, %d iters, %d inf, %d not finished, %d not trusted\n",
3802                   numberNodes, numberStrongDone, numberStrongIterations, numberStrongInfeasible, numberUnfinished,
3803                   numberNotTrusted);
3804        // decide what to do
3805        if (numberUnfinished*10 > numberStrongDone + 1 ||
3806                !numberStrongInfeasible) {
3807          COIN_DETAIL_PRINT(printf("going to strategy 2\n"));
3808            // Weaken
3809            model->setNumberStrong(2);
3810            model->setNumberBeforeTrust(1);
3811            model->synchronizeNumberBeforeTrust();
3812            model->setSearchStrategy(2);
3813        }
3814    }
3815    if (numberUnfinished*10 < numberStrongDone &&
3816        model->numberStrongIterations()*20 < model->getIterationCount()&&
3817                                !auxiliaryInfo->solutionAddsCuts()) {
3818        //printf("increasing trust\n");
3819        model->synchronizeNumberBeforeTrust(2);
3820    }
3821
3822    // Set guessed solution value
3823    guessedObjectiveValue_ = objectiveValue_ + estimatedDegradation;
3824    int kColumn=-1;
3825    if (branch_) {
3826      CbcObject * obj = (dynamic_cast<CbcBranchingObject *>(branch_))->object(); 
3827      CbcSimpleInteger * branchObj =
3828        dynamic_cast <CbcSimpleInteger *>(obj) ;
3829      if (branchObj) {
3830        kColumn=branchObj->columnNumber();
3831      }
3832    }
3833#ifdef COIN_HAS_NTY
3834    if (orbitOption&&kColumn>=0) {
3835      CbcSymmetry * symmetryInfo = model->symmetryInfo();
3836      CbcNodeInfo * infoX = lastNode ? lastNode->nodeInfo() : NULL;
3837      bool worthTrying = false;
3838      if (infoX) {
3839        CbcNodeInfo * info = infoX;
3840        for (int i=0;i<NTY_BAD_DEPTH;i++) {
3841          if (!info->parent()) {
3842            worthTrying = true;
3843            break;
3844          }
3845          info = info->parent();
3846          if (info->symmetryWorked()) {
3847            worthTrying = true;
3848            break;
3849          }
3850        }
3851      } else {
3852        worthTrying=true;
3853      }
3854      if (symmetryInfo && worthTrying) {
3855        if (orbitOption==1) {
3856          symmetryInfo->ChangeBounds(solver->getColLower(),
3857                                     solver->getColUpper(),
3858                                     solver->getNumCols(),false);
3859          symmetryInfo->Compute_Symmetry();
3860          symmetryInfo->fillOrbits();
3861        }
3862        const int * orbits = symmetryInfo->whichOrbit();
3863        if (orbits && orbits[kColumn]>=0) {
3864          int numberUsefulOrbits = symmetryInfo->numberUsefulObjects();
3865          if (numberUsefulOrbits<1000) {
3866            delete branch_;
3867            branch_ = new CbcOrbitalBranchingObject(model,kColumn,1,0,NULL);
3868            if (infoX)
3869              infoX->setSymmetryWorked();
3870          }
3871        }
3872      }
3873    }
3874#endif
3875    if (model->logLevel()>1)
3876    printf ("Node %d depth %d unsatisfied %d sum %g obj %g guess %g branching on %d\n",
3877            model->getNodeCount(),depth_,numberUnsatisfied_,
3878            sumInfeasibilities_,objectiveValue_,guessedObjectiveValue_,
3879            kColumn);
3880#ifdef DO_ALL_AT_ROOT
3881    if (strongType) {
3882      char general[200];
3883      if (strongType==1)
3884        sprintf(general,"Strong branching on all %d unsatisfied, %d iterations (depth %d)\n",
3885                saveNumberToDo,numberStrongIterations,depth_);
3886      else
3887        sprintf(general,"Strong branching on all %d unfixed variables (%d unsatisfied), %d iterations (depth %d)\n",
3888                saveNumberToDo+saveSatisfiedVariables,saveNumberToDo,numberStrongIterations,depth_);
3889      model->messageHandler()->message(CBC_FPUMP2,model->messages())
3890        << general << CoinMessageEol ;
3891    }
3892#endif
3893#ifdef DEBUG_SOLUTION
3894    if(onOptimalPath&&anyAction==-2) {
3895      printf("Gone off optimal path in CbcNode\n");
3896      assert(!onOptimalPath||anyAction!=-2);
3897    }
3898#endif
3899    /*
3900      Cleanup, then we're finished
3901    */
3902    if (!model->branchingMethod())
3903        delete decision;
3904
3905    delete choiceObject;
3906    delete [] sort;
3907    delete [] whichObject;
3908#ifdef RANGING
3909    delete [] objectMark;
3910#endif
3911    delete [] saveLower;
3912    delete [] saveUpper;
3913    delete [] upEstimate;
3914    delete [] downEstimate;
3915# ifdef COIN_HAS_CLP
3916    if (osiclp) {
3917        osiclp->setSpecialOptions(saveClpOptions);
3918    }
3919# endif
3920    // restore solution
3921    solver->setColSolution(saveSolution);
3922    model->reserveCurrentSolution(saveSolution);
3923    delete [] saveSolution;
3924    model->setStateOfSearch(saveStateOfSearch);
3925    model->setLogLevel(saveLogLevel);
3926    // delete extra regions
3927    if (usefulInfo.usefulRegion_) {
3928        delete [] usefulInfo.usefulRegion_;
3929        delete [] usefulInfo.indexRegion_;
3930        delete [] usefulInfo.pi_;
3931        usefulInfo.usefulRegion_ = NULL;
3932        usefulInfo.indexRegion_ = NULL;
3933        usefulInfo.pi_ = NULL;
3934    }
3935    useShadow = model->moreSpecialOptions() & 7;
3936    if ((useShadow == 5 && model->getSolutionCount()) || useShadow == 6) {
3937        // zap pseudo shadow prices
3938        model->pseudoShadow(-1);
3939        // and switch off
3940        model->setMoreSpecialOptions(model->moreSpecialOptions()&(~1023));
3941    }
3942    return anyAction;
3943}
3944int CbcNode::analyze (CbcModel *model, double * results)
3945{
3946    int i;
3947    int numberIterationsAllowed = model->numberAnalyzeIterations();
3948    OsiSolverInterface * solver = model->solver();
3949    objectiveValue_ = solver->getObjSense() * solver->getObjValue();
3950    double cutoff = model->getCutoff();
3951    const double * lower = solver->getColLower();
3952    const double * upper = solver->getColUpper();
3953    const double * dj = solver->getReducedCost();
3954    int numberObjects = model->numberObjects();
3955    int numberColumns = model->getNumCols();
3956    // Initialize arrays
3957    int numberIntegers = model->numberIntegers();
3958    int * back = new int[numberColumns];
3959    const int * integerVariable = model->integerVariable();
3960    for (i = 0; i < numberColumns; i++)
3961        back[i] = -1;
3962    // What results is
3963    double * newLower = results;
3964    double * objLower = newLower + numberIntegers;
3965    double * newUpper = objLower + numberIntegers;
3966    double * objUpper = newUpper + numberIntegers;
3967    for (i = 0; i < numberIntegers; i++) {
3968        int iColumn = integerVariable[i];
3969        back[iColumn] = i;
3970        newLower[i] = 0.0;
3971        objLower[i] = -COIN_DBL_MAX;
3972        newUpper[i] = 0.0;
3973        objUpper[i] = -COIN_DBL_MAX;
3974    }
3975    double * saveUpper = new double[numberColumns];
3976    double * saveLower = new double[numberColumns];
3977    // Save solution in case heuristics need good solution later
3978
3979    double * saveSolution = new double[numberColumns];
3980    memcpy(saveSolution, solver->getColSolution(), numberColumns*sizeof(double));
3981    model->reserveCurrentSolution(saveSolution);
3982    for (i = 0; i < numberColumns; i++) {
3983        saveLower[i] = lower[i];
3984        saveUpper[i] = upper[i];
3985    }
3986    // Get arrays to sort
3987    double * sort = new double[numberObjects];
3988    int * whichObject = new int[numberObjects];
3989    int numberToFix = 0;
3990    int numberToDo = 0;
3991    double integerTolerance =
3992        model->getDblParam(CbcModel::CbcIntegerTolerance);
3993    // point to useful information
3994    OsiBranchingInformation usefulInfo = model->usefulInformation();
3995    // and modify
3996    usefulInfo.depth_ = depth_;
3997
3998    // compute current state
3999    int numberObjectInfeasibilities; // just odd ones
4000    int numberIntegerInfeasibilities;
4001    model->feasibleSolution(
4002        numberIntegerInfeasibilities,
4003        numberObjectInfeasibilities);
4004# ifdef COIN_HAS_CLP
4005    OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
4006    int saveClpOptions = 0;
4007    bool fastIterations = (model->specialOptions() & 8) != 0;
4008    if (osiclp && fastIterations) {
4009        // for faster hot start
4010        saveClpOptions = osiclp->specialOptions();
4011        osiclp->setSpecialOptions(saveClpOptions | 8192);
4012    }
4013# else
4014    bool fastIterations = false ;
4015# endif
4016    /*
4017      Scan for branching objects that indicate infeasibility. Choose candidates
4018      using priority as the first criteria, then integer infeasibility.
4019
4020      The algorithm is to fill the array with a set of good candidates (by
4021      infeasibility) with priority bestPriority.  Finding a candidate with
4022      priority better (less) than bestPriority flushes the choice array. (This
4023      serves as initialization when the first candidate is found.)
4024
4025    */
4026    numberToDo = 0;
4027    for (i = 0; i < numberObjects; i++) {
4028        OsiObject * object = model->modifiableObject(i);
4029        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
4030            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
4031        if (!dynamicObject)
4032            continue;
4033        int preferredWay;
4034        double infeasibility = object->infeasibility(&usefulInfo, preferredWay);
4035        int iColumn = dynamicObject->columnNumber();
4036        if (saveUpper[iColumn] == saveLower[iColumn])
4037            continue;
4038        if (infeasibility)
4039            sort[numberToDo] = -1.0e10 - infeasibility;
4040        else
4041            sort[numberToDo] = -fabs(dj[iColumn]);
4042        whichObject[numberToDo++] = i;
4043    }
4044    // Save basis
4045    CoinWarmStart * ws = solver->getWarmStart();
4046    int saveLimit;
4047    solver->getIntParam(OsiMaxNumIterationHotStart, saveLimit);
4048    int targetIterations = CoinMax(500, numberIterationsAllowed / numberObjects);
4049    if (saveLimit < targetIterations)
4050        solver->setIntParam(OsiMaxNumIterationHotStart, targetIterations);
4051    // Mark hot start
4052    solver->markHotStart();
4053    // Sort
4054    CoinSort_2(sort, sort + numberToDo, whichObject);
4055    double * currentSolution = model->currentSolution();
4056    double objMin = 1.0e50;
4057    double objMax = -1.0e50;
4058    bool needResolve = false;
4059    /*
4060      Now calculate the cost forcing the variable up and down.
4061    */
4062    int iDo;
4063    for (iDo = 0; iDo < numberToDo; iDo++) {
4064        CbcStrongInfo choice;
4065        int iObject = whichObject[iDo];
4066        OsiObject * object = model->modifiableObject(iObject);
4067        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
4068            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
4069        if (!dynamicObject)
4070            continue;
4071        int iColumn = dynamicObject->columnNumber();
4072        int preferredWay;
4073        /*
4074          Update the information held in the object.
4075        */
4076        object->infeasibility(&usefulInfo, preferredWay);
4077        double value = currentSolution[iColumn];
4078        double nearest = floor(value + 0.5);
4079        double lowerValue = floor(value);
4080        bool satisfied = false;
4081        if (fabs(value - nearest) <= integerTolerance || value < saveLower[iColumn] || value > saveUpper[iColumn]) {
4082            satisfied = true;
4083            double newValue;
4084            if (nearest < saveUpper[iColumn]) {
4085                newValue = nearest + 1.0001 * integerTolerance;
4086                lowerValue = nearest;
4087            } else {
4088                newValue = nearest - 1.0001 * integerTolerance;
4089                lowerValue = nearest - 1;
4090            }
4091            currentSolution[iColumn] = newValue;
4092        }
4093        double upperValue = lowerValue + 1.0;
4094        //CbcSimpleInteger * obj =
4095        //dynamic_cast <CbcSimpleInteger *>(object) ;
4096        //if (obj) {
4097        //choice.possibleBranch=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
4098        //} else {
4099        CbcObject * obj =
4100            dynamic_cast <CbcObject *>(object) ;
4101        assert (obj);
4102        choice.possibleBranch = obj->createCbcBranch(solver, &usefulInfo, preferredWay);
4103        //}
4104        currentSolution[iColumn] = value;
4105        // Save which object it was
4106        choice.objectNumber = iObject;
4107        choice.numIntInfeasUp = numberUnsatisfied_;
4108        choice.numIntInfeasDown = numberUnsatisfied_;
4109        choice.downMovement = 0.0;
4110        choice.upMovement = 0.0;
4111        choice.numItersDown = 0;
4112        choice.numItersUp = 0;
4113        choice.fix = 0; // say not fixed
4114        double objectiveChange ;
4115        double newObjectiveValue = 1.0e100;
4116        int j;
4117        // status is 0 finished, 1 infeasible and other
4118        int iStatus;
4119        /*
4120          Try the down direction first. (Specify the initial branching alternative as
4121          down with a call to way(-1). Each subsequent call to branch() performs the
4122          specified branch and advances the branch object state to the next branch
4123          alternative.)
4124        */
4125        choice.possibleBranch->way(-1) ;
4126        choice.possibleBranch->branch() ;
4127        if (fabs(value - lowerValue) > integerTolerance) {
4128            solver->solveFromHotStart() ;
4129            /*
4130              We now have an estimate of objective degradation that we can use for strong
4131              branching. If we're over the cutoff, the variable is monotone up.
4132              If we actually made it to optimality, check for a solution, and if we have
4133              a good one, call setBestSolution to process it. Note that this may reduce the
4134              cutoff, so we check again to see if we can declare this variable monotone.
4135            */
4136            if (solver->isProvenOptimal())
4137                iStatus = 0; // optimal
4138            else if (solver->isIterationLimitReached()
4139                     && !solver->isDualObjectiveLimitReached())
4140                iStatus = 2; // unknown
4141            else
4142                iStatus = 1; // infeasible
4143            newObjectiveValue = solver->getObjSense() * solver->getObjValue();
4144            choice.numItersDown = solver->getIterationCount();
4145            numberIterationsAllowed -= choice.numItersDown;
4146            objectiveChange = newObjectiveValue  - objectiveValue_;
4147            if (!iStatus) {
4148                choice.finishedDown = true ;
4149                if (newObjectiveValue >= cutoff) {
4150                    objectiveChange = 1.0e100; // say infeasible
4151                } else {
4152                    // See if integer solution
4153                    if (model->feasibleSolution(choice.numIntInfeasDown,
4154                                                choice.numObjInfeasDown)
4155                            && model->problemFeasibility()->feasible(model, -1) >= 0) {
4156                        model->setBestSolution(CBC_STRONGSOL,
4157                                               newObjectiveValue,
4158                                               solver->getColSolution()) ;
4159                        model->setLastHeuristic(NULL);
4160                        model->incrementUsed(solver->getColSolution());
4161                        cutoff = model->getCutoff();
4162                        if (newObjectiveValue >= cutoff)        //  *new* cutoff
4163                            objectiveChange = 1.0e100 ;
4164                    }
4165                }
4166            } else if (iStatus == 1) {
4167                objectiveChange = 1.0e100 ;
4168            } else {
4169                // Can't say much as we did not finish
4170                choice.finishedDown = false ;
4171            }
4172            choice.downMovement = objectiveChange ;
4173        }
4174        // restore bounds
4175        for ( j = 0; j < numberColumns; j++) {
4176            if (saveLower[j] != lower[j])
4177                solver->setColLower(j, saveLower[j]);
4178            if (saveUpper[j] != upper[j])
4179                solver->setColUpper(j, saveUpper[j]);
4180        }
4181        // repeat the whole exercise, forcing the variable up
4182        choice.possibleBranch->branch();
4183        if (fabs(value - upperValue) > integerTolerance) {
4184            solver->solveFromHotStart() ;
4185            /*
4186              We now have an estimate of objective degradation that we can use for strong
4187              branching. If we're over the cutoff, the variable is monotone up.
4188              If we actually made it to optimality, check for a solution, and if we have
4189              a good one, call setBestSolution to process it. Note that this may reduce the
4190              cutoff, so we check again to see if we can declare this variable monotone.
4191            */
4192            if (solver->isProvenOptimal())
4193                iStatus = 0; // optimal
4194            else if (solver->isIterationLimitReached()
4195                     && !solver->isDualObjectiveLimitReached())
4196                iStatus = 2; // unknown
4197            else
4198                iStatus = 1; // infeasible
4199            newObjectiveValue = solver->getObjSense() * solver->getObjValue();
4200            choice.numItersUp = solver->getIterationCount();
4201            numberIterationsAllowed -= choice.numItersUp;
4202            objectiveChange = newObjectiveValue  - objectiveValue_;
4203            if (!iStatus) {
4204                choice.finishedUp = true ;
4205                if (newObjectiveValue >= cutoff) {
4206                    objectiveChange = 1.0e100; // say infeasible
4207                } else {
4208                    // See if integer solution
4209                    if (model->feasibleSolution(choice.numIntInfeasUp,
4210                                                choice.numObjInfeasUp)
4211                            && model->problemFeasibility()->feasible(model, -1) >= 0) {
4212                        model->setBestSolution(CBC_STRONGSOL,
4213                                               newObjectiveValue,
4214                                               solver->getColSolution()) ;
4215                        model->setLastHeuristic(NULL);
4216                        model->incrementUsed(solver->getColSolution());
4217                        cutoff = model->getCutoff();
4218                        if (newObjectiveValue >= cutoff)        //  *new* cutoff
4219                            objectiveChange = 1.0e100 ;
4220                    }
4221                }
4222            } else if (iStatus == 1) {
4223                objectiveChange = 1.0e100 ;
4224            } else {
4225                // Can't say much as we did not finish
4226                choice.finishedUp = false ;
4227            }
4228            choice.upMovement = objectiveChange ;
4229
4230            // restore bounds
4231            for ( j = 0; j < numberColumns; j++) {
4232                if (saveLower[j] != lower[j])
4233                    solver->setColLower(j, saveLower[j]);
4234                if (saveUpper[j] != upper[j])
4235                    solver->setColUpper(j, saveUpper[j]);
4236            }
4237        }
4238        // If objective goes above certain amount we can set bound
4239        int jInt = back[iColumn];
4240        newLower[jInt] = upperValue;
4241        if (choice.finishedDown)
4242            objLower[jInt] = choice.downMovement + objectiveValue_;
4243        else
4244            objLower[jInt] = objectiveValue_;
4245        newUpper[jInt] = lowerValue;
4246        if (choice.finishedUp)
4247            objUpper[jInt] = choice.upMovement + objectiveValue_;
4248        else
4249            objUpper[jInt] = objectiveValue_;
4250        objMin = CoinMin(CoinMin(objLower[jInt], objUpper[jInt]), objMin);
4251        /*
4252          End of evaluation for this candidate variable. Possibilities are:
4253          * Both sides below cutoff; this variable is a candidate for branching.
4254          * Both sides infeasible or above the objective cutoff: no further action
4255          here. Break from the evaluation loop and assume the node will be purged
4256          by the caller.
4257          * One side below cutoff: Install the branch (i.e., fix the variable). Break
4258          from the evaluation loop and assume the node will be reoptimised by the
4259          caller.
4260        */
4261        if (choice.upMovement < 1.0e100) {
4262            if (choice.downMovement < 1.0e100) {
4263                objMax = CoinMax(CoinMax(objLower[jInt], objUpper[jInt]), objMax);
4264                // In case solution coming in was odd
4265                choice.upMovement = CoinMax(0.0, choice.upMovement);
4266                choice.downMovement = CoinMax(0.0, choice.downMovement);
4267                // feasible -
4268                model->messageHandler()->message(CBC_STRONG, *model->messagesPointer())
4269                << iObject << iColumn
4270                << choice.downMovement << choice.numIntInfeasDown
4271                << choice.upMovement << choice.numIntInfeasUp
4272                << value
4273                << CoinMessageEol;
4274            } else {
4275                // up feasible, down infeasible
4276                if (!satisfied)
4277                    needResolve = true;
4278                choice.fix = 1;
4279                numberToFix++;
4280                saveLower[iColumn] = upperValue;
4281                solver->setColLower(iColumn, upperValue);
4282            }
4283        } else {
4284            if (choice.downMovement < 1.0e100) {
4285                // down feasible, up infeasible
4286                if (!satisfied)
4287                    needResolve = true;
4288                choice.fix = -1;
4289                numberToFix++;
4290                saveUpper[iColumn] = lowerValue;
4291                solver->setColUpper(iColumn, lowerValue);
4292            } else {
4293                // neither side feasible
4294                COIN_DETAIL_PRINT(printf("Both infeasible for choice %d sequence %d\n", i,
4295                                         model->object(choice.objectNumber)->columnNumber()));
4296                delete ws;
4297                ws = NULL;
4298                //solver->writeMps("bad");
4299                numberToFix = -1;
4300                delete choice.possibleBranch;
4301                choice.possibleBranch = NULL;
4302                break;
4303            }
4304        }
4305        delete choice.possibleBranch;
4306        if (numberIterationsAllowed <= 0)
4307            break;
4308        //printf("obj %d, col %d, down %g up %g value %g\n",iObject,iColumn,
4309        //     choice.downMovement,choice.upMovement,value);
4310    }
4311    COIN_DETAIL_PRINT(printf("Best possible solution %g, can fix more if solution of %g found - looked at %d variables in %d iterations\n",
4312                             objMin, objMax, iDo, model->numberAnalyzeIterations() - numberIterationsAllowed));
4313    model->setNumberAnalyzeIterations(numberIterationsAllowed);
4314    // Delete the snapshot
4315    solver->unmarkHotStart();
4316    // back to normal
4317    solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
4318    solver->setIntParam(OsiMaxNumIterationHotStart, saveLimit);
4319    // restore basis
4320    solver->setWarmStart(ws);
4321    delete ws;
4322
4323    delete [] sort;
4324    delete [] whichObject;
4325    delete [] saveLower;
4326    delete [] saveUpper;
4327    delete [] back;
4328    // restore solution
4329    solver->setColSolution(saveSolution);
4330# ifdef COIN_HAS_CLP
4331    if (osiclp)
4332        osiclp->setSpecialOptions(saveClpOptions);
4333# endif
4334    model->reserveCurrentSolution(saveSolution);
4335    delete [] saveSolution;
4336    if (needResolve)
4337        solver->resolve();
4338    return numberToFix;
4339}
4340
4341
4342CbcNode::CbcNode(const CbcNode & rhs)
4343        : CoinTreeNode(rhs)
4344{
4345#ifdef CHECK_NODE
4346    printf("CbcNode %x Constructor from rhs %x\n", this, &rhs);
4347#endif
4348    if (rhs.nodeInfo_)
4349        nodeInfo_ = rhs.nodeInfo_->clone();
4350    else
4351        nodeInfo_ = NULL;
4352    objectiveValue_ = rhs.objectiveValue_;
4353    guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
4354    sumInfeasibilities_ = rhs.sumInfeasibilities_;
4355    if (rhs.branch_)
4356        branch_ = rhs.branch_->clone();
4357    else
4358        branch_ = NULL;
4359    depth_ = rhs.depth_;
4360    numberUnsatisfied_ = rhs.numberUnsatisfied_;
4361    nodeNumber_ = rhs.nodeNumber_;
4362    state_ = rhs.state_;
4363    if (nodeInfo_)
4364        assert ((state_&2) != 0);
4365    else
4366        assert ((state_&2) == 0);
4367}
4368
4369CbcNode &
4370CbcNode::operator=(const CbcNode & rhs)
4371{
4372    if (this != &rhs) {
4373        delete nodeInfo_;
4374        if (rhs.nodeInfo_)
4375            nodeInfo_ = rhs.nodeInfo_->clone();
4376        else
4377            nodeInfo_ = NULL;
4378        objectiveValue_ = rhs.objectiveValue_;
4379        guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
4380        sumInfeasibilities_ = rhs.sumInfeasibilities_;
4381        if (rhs.branch_)
4382            branch_ = rhs.branch_->clone();
4383        else
4384            branch_ = NULL,
4385                      depth_ = rhs.depth_;
4386        numberUnsatisfied_ = rhs.numberUnsatisfied_;
4387        nodeNumber_ = rhs.nodeNumber_;
4388        state_ = rhs.state_;
4389        if (nodeInfo_)
4390            assert ((state_&2) != 0);
4391        else
4392            assert ((state_&2) == 0);
4393    }
4394    return *this;
4395}
4396CbcNode::~CbcNode ()
4397{
4398#ifdef CHECK_NODE
4399    if (nodeInfo_) {
4400        printf("CbcNode %x Destructor nodeInfo %x (%d)\n",
4401               this, nodeInfo_, nodeInfo_->numberPointingToThis());
4402        //assert(nodeInfo_->numberPointingToThis()>=0);
4403    } else {
4404        printf("CbcNode %x Destructor nodeInfo %x (?)\n",
4405               this, nodeInfo_);
4406    }
4407#endif
4408    if (nodeInfo_) {
4409        // was if (nodeInfo_&&(state_&2)!=0) {
4410        nodeInfo_->nullOwner();
4411        int numberToDelete = nodeInfo_->numberBranchesLeft();
4412        //    CbcNodeInfo * parent = nodeInfo_->parent();
4413        //assert (nodeInfo_->numberPointingToThis()>0);
4414        if (nodeInfo_->decrement(numberToDelete) == 0 || (state_&2) == 0) {
4415            if ((state_&2) == 0)
4416                nodeInfo_->nullParent();
4417            delete nodeInfo_;
4418        } else {
4419            //printf("node %x nodeinfo %x parent %x\n",this,nodeInfo_,nodeInfo_->parent());
4420            // anyway decrement parent
4421            //if (parent)
4422            ///parent->decrement(1);
4423        }
4424    }
4425    delete branch_;
4426}
4427// Decrement  active cut counts
4428void
4429CbcNode::decrementCuts(int change)
4430{
4431    if (nodeInfo_)
4432        assert ((state_&2) != 0);
4433    else
4434        assert ((state_&2) == 0);
4435    if (nodeInfo_) {
4436        nodeInfo_->decrementCuts(change);
4437    }
4438}
4439void
4440CbcNode::decrementParentCuts(CbcModel * model, int change)
4441{
4442    if (nodeInfo_)
4443        assert ((state_&2) != 0);
4444    else
4445        assert ((state_&2) == 0);
4446    if (nodeInfo_) {
4447        nodeInfo_->decrementParentCuts(model, change);
4448    }
4449}
4450
4451/*
4452  Initialize reference counts (numberPointingToThis, numberBranchesLeft_)
4453  in the attached nodeInfo_.
4454*/
4455void
4456CbcNode::initializeInfo()
4457{
4458    assert(nodeInfo_ && branch_) ;
4459    nodeInfo_->initializeInfo(branch_->numberBranches());
4460    assert ((state_&2) != 0);
4461    assert (nodeInfo_->numberBranchesLeft() ==
4462            branch_->numberBranchesLeft());
4463}
4464// Nulls out node info
4465void
4466CbcNode::nullNodeInfo()
4467{
4468    nodeInfo_ = NULL;
4469    // say not active
4470    state_ &= ~2;
4471}
4472
4473int
4474CbcNode::branch(OsiSolverInterface * solver)
4475{
4476    double changeInGuessed;
4477    assert (nodeInfo_->numberBranchesLeft() ==
4478            branch_->numberBranchesLeft());
4479    if (!solver)
4480        changeInGuessed = branch_->branch();
4481    else
4482        changeInGuessed = branch_->branch(solver);
4483    guessedObjectiveValue_ += changeInGuessed;
4484    //#define PRINTIT
4485#ifdef PRINTIT
4486    int numberLeft = nodeInfo_->numberBranchesLeft();
4487    CbcNodeInfo * parent = nodeInfo_->parent();
4488    int parentNodeNumber = -1;
4489    CbcBranchingObject * object1 =
4490        dynamic_cast<CbcBranchingObject *>(branch_) ;
4491    //OsiObject * object = object1->
4492    //int sequence = object->columnNumber);
4493    int id = -1;
4494    double value = 0.0;
4495    if (object1) {
4496        id = object1->variable();
4497        value = object1->value();
4498    }
4499    printf("id %d value %g objvalue %g\n", id, value, objectiveValue_);
4500    if (parent)
4501        parentNodeNumber = parent->nodeNumber();
4502    printf("Node number %d, %s, way %d, depth %d, parent node number %d\n",
4503           nodeInfo_->nodeNumber(), (numberLeft == 2) ? "leftBranch" : "rightBranch",
4504           way(), depth_, parentNodeNumber);
4505    assert (parentNodeNumber != nodeInfo_->nodeNumber());
4506#endif
4507    return nodeInfo_->branchedOn();
4508}
4509/* Active arm of the attached OsiBranchingObject.
4510
4511   In the simplest instance, coded -1 for the down arm of the branch, +1 for
4512   the up arm. But see OsiBranchingObject::way()
4513   Use nodeInfo--.numberBranchesLeft_ to see how active
4514
4515   Except that there is no OsiBranchingObject::way(), and this'll fail in any
4516   event because we have various OsiXXXBranchingObjects which aren't descended
4517   from CbcBranchingObjects. I think branchIndex() is the appropriate
4518   equivalent, but could be wrong. (lh, 061220)
4519
4520   071212: I'm finally getting back to cbc-generic and rescuing a lot of my
4521   annotation from branches/devel (which was killed in summer). I'm going to
4522   put back an assert(obj) just to see what happens. It's still present as of
4523   the most recent change to CbcNode (r833).
4524
4525   080104: Yep, we can arrive here with an OsiBranchingObject. Removed the
4526   assert, it's served its purpose.
4527
4528   080226: John finally noticed this problem and added a way() method to the
4529   OsiBranchingObject hierarchy. Removing my workaround.
4530
4531*/
4532int
4533CbcNode::way() const
4534{
4535    if (branch_) {
4536        CbcBranchingObject * obj =
4537            dynamic_cast <CbcBranchingObject *>(branch_) ;
4538        if (obj) {
4539            return obj->way();
4540        } else {
4541            OsiTwoWayBranchingObject * obj2 =
4542                dynamic_cast <OsiTwoWayBranchingObject *>(branch_) ;
4543            assert (obj2);
4544            return obj2->way();
4545        }
4546    } else {
4547        return 0;
4548    }
4549}
4550/* Create a branching object for the node
4551
4552   The routine scans the object list of the model and selects a set of
4553   unsatisfied objects as candidates for branching. The candidates are
4554   evaluated, and an appropriate branch object is installed.
4555
4556   The numberPassesLeft is decremented to stop fixing one variable each time
4557   and going on and on (e.g. for stock cutting, air crew scheduling)
4558
4559   If evaluation determines that an object is monotone or infeasible,
4560   the routine returns immediately. In the case of a monotone object,
4561   the branch object has already been called to modify the model.
4562
4563   Return value:
4564   <ul>
4565   <li>  0: A branching object has been installed
4566   <li> -1: A monotone object was discovered
4567   <li> -2: An infeasible object was discovered
4568   </ul>
4569   Branch state:
4570   <ul>
4571   <li> -1: start
4572   <li> -1: A monotone object was discovered
4573   <li> -2: An infeasible object was discovered
4574   </ul>
4575*/
4576int
4577CbcNode::chooseOsiBranch (CbcModel * model,
4578                          CbcNode * lastNode,
4579                          OsiBranchingInformation * usefulInfo,
4580                          int branchState)
4581{
4582    int returnStatus = 0;
4583    if (lastNode)
4584        depth_ = lastNode->depth_ + 1;
4585    else
4586        depth_ = 0;
4587    OsiSolverInterface * solver = model->solver();
4588    objectiveValue_ = solver->getObjValue() * solver->getObjSense();
4589    usefulInfo->objectiveValue_ = objectiveValue_;
4590    usefulInfo->depth_ = depth_;
4591    const double * saveInfoSol = usefulInfo->solution_;
4592    double * saveSolution = new double[solver->getNumCols()];
4593    memcpy(saveSolution, solver->getColSolution(), solver->getNumCols()*sizeof(double));
4594    usefulInfo->solution_ = saveSolution;
4595    OsiChooseVariable * choose = model->branchingMethod()->chooseMethod();
4596    int numberUnsatisfied = -1;
4597    if (branchState < 0) {
4598        // initialize
4599        // initialize sum of "infeasibilities"
4600        sumInfeasibilities_ = 0.0;
4601        numberUnsatisfied = choose->setupList(usefulInfo, true);
4602        numberUnsatisfied_ = numberUnsatisfied;
4603        branchState = 0;
4604        if (numberUnsatisfied_ < 0) {
4605            // infeasible
4606            delete [] saveSolution;
4607            return -2;
4608        }
4609    }
4610    // unset best
4611    int best = -1;
4612    choose->setBestObjectIndex(-1);
4613    if (numberUnsatisfied) {
4614        if (branchState > 0 || !choose->numberOnList()) {
4615            // we need to return at once - don't do strong branching or anything
4616            if (choose->numberOnList() || !choose->numberStrong()) {
4617                best = choose->candidates()[0];
4618                choose->setBestObjectIndex(best);
4619            } else {
4620                // nothing on list - need to try again - keep any solution
4621                numberUnsatisfied = choose->setupList(usefulInfo, false);
4622                numberUnsatisfied_ = numberUnsatisfied;
4623                if (numberUnsatisfied) {
4624                    best = choose->candidates()[0];
4625                    choose->setBestObjectIndex(best);
4626                }
4627            }
4628        } else {
4629            // carry on with strong branching or whatever
4630            int returnCode = choose->chooseVariable(solver, usefulInfo, true);
4631            // update number of strong iterations etc
4632            model->incrementStrongInfo(choose->numberStrongDone(), choose->numberStrongIterations(),
4633                                       returnCode == -1 ? 0 : choose->numberStrongFixed(), returnCode == -1);
4634            if (returnCode > 1) {
4635                // has fixed some
4636                returnStatus = -1;
4637            } else if (returnCode == -1) {
4638                // infeasible
4639                returnStatus = -2;
4640            } else if (returnCode == 0) {
4641                // normal
4642                returnStatus = 0;
4643                numberUnsatisfied = 1;
4644            } else {
4645                // ones on list satisfied - double check
4646                numberUnsatisfied = choose->setupList(usefulInfo, false);
4647                numberUnsatisfied_ = numberUnsatisfied;
4648                if (numberUnsatisfied) {
4649                    best = choose->candidates()[0];
4650                    choose->setBestObjectIndex(best);
4651                }
4652            }
4653        }
4654    }
4655    delete branch_;
4656    branch_ = NULL;
4657    guessedObjectiveValue_ = COIN_DBL_MAX;//objectiveValue_; // for now
4658    if (!returnStatus) {
4659        if (numberUnsatisfied) {
4660            // create branching object
4661            const OsiObject * obj = model->solver()->object(choose->bestObjectIndex());
4662            //const OsiSolverInterface * solver = usefulInfo->solver_;
4663            branch_ = obj->createBranch(model->solver(), usefulInfo, obj->whichWay());
4664        }
4665    }
4666    usefulInfo->solution_ = saveInfoSol;
4667    delete [] saveSolution;
4668    // may have got solution
4669    if (choose->goodSolution()
4670            && model->problemFeasibility()->feasible(model, -1) >= 0) {
4671        // yes
4672        double objValue = choose->goodObjectiveValue();
4673        model->setBestSolution(CBC_STRONGSOL,
4674                               objValue,
4675                               choose->goodSolution()) ;
4676        model->setLastHeuristic(NULL);
4677        model->incrementUsed(choose->goodSolution());
4678        choose->clearGoodSolution();
4679    }
4680    return returnStatus;
4681}
4682int
4683CbcNode::chooseClpBranch (CbcModel * model,
4684                          CbcNode * lastNode)
4685{
4686    assert(lastNode);
4687    depth_ = lastNode->depth_ + 1;
4688    delete branch_;
4689    branch_ = NULL;
4690    OsiSolverInterface * solver = model->solver();
4691    //double saveObjectiveValue = solver->getObjValue();
4692    //double objectiveValue = CoinMax(solver->getObjSense()*saveObjectiveValue,objectiveValue_);
4693    const double * lower = solver->getColLower();
4694    const double * upper = solver->getColUpper();
4695    // point to useful information
4696    OsiBranchingInformation usefulInfo = model->usefulInformation();
4697    // and modify
4698    usefulInfo.depth_ = depth_;
4699    int i;
4700    //bool beforeSolution = model->getSolutionCount()==0;
4701    int numberObjects = model->numberObjects();
4702    int numberColumns = model->getNumCols();
4703    double * saveUpper = new double[numberColumns];
4704    double * saveLower = new double[numberColumns];
4705
4706    // Save solution in case heuristics need good solution later
4707
4708    double * saveSolution = new double[numberColumns];
4709    memcpy(saveSolution, solver->getColSolution(), numberColumns*sizeof(double));
4710    model->reserveCurrentSolution(saveSolution);
4711    for (i = 0; i < numberColumns; i++) {
4712        saveLower[i] = lower[i];
4713        saveUpper[i] = upper[i];
4714    }
4715    // Save basis
4716    CoinWarmStart * ws = solver->getWarmStart();
4717    numberUnsatisfied_ = 0;
4718    // initialize sum of "infeasibilities"
4719    sumInfeasibilities_ = 0.0;
4720    // Note looks as if off end (hidden one)
4721    OsiObject * object = model->modifiableObject(numberObjects);
4722    CbcGeneralDepth * thisOne = dynamic_cast <CbcGeneralDepth *> (object);
4723    assert (thisOne);
4724    OsiClpSolverInterface * clpSolver
4725    = dynamic_cast<OsiClpSolverInterface *> (solver);
4726    assert (clpSolver);
4727    ClpSimplex * simplex = clpSolver->getModelPtr();
4728    int preferredWay;
4729    double infeasibility = object->infeasibility(&usefulInfo, preferredWay);
4730    if (thisOne->whichSolution() >= 0) {
4731        ClpNode * nodeInfo=NULL;
4732        if ((model->moreSpecialOptions()&33554432)==0) {
4733          nodeInfo = thisOne->nodeInfo(thisOne->whichSolution());
4734          nodeInfo->applyNode(simplex, 2);
4735        } else {
4736          // from diving
4737          CbcSubProblem ** nodes = reinterpret_cast<CbcSubProblem **>
4738            (model->temporaryPointer());
4739          assert (nodes);
4740          int numberDo=thisOne->numberNodes()-1;
4741          for (int iNode=0;iNode<numberDo;iNode++)
4742            nodes[iNode]->apply(solver,1);
4743          nodes[numberDo]->apply(solver,9+16);
4744        }
4745        int saveLogLevel = simplex->logLevel();
4746        simplex->setLogLevel(0);
4747        simplex->dual();
4748        simplex->setLogLevel(saveLogLevel);
4749        double cutoff = model->getCutoff();
4750        bool goodSolution = true;
4751        if (simplex->status()) {
4752            //simplex->writeMps("bad7.mps",2);
4753            if (nodeInfo) {
4754              if (nodeInfo->objectiveValue() > cutoff - 1.0e-2)
4755                goodSolution = false;
4756              else
4757                assert (!simplex->status());
4758            } else {
4759              // debug diving
4760              assert (!simplex->status());
4761            }
4762        }
4763        if (goodSolution) {
4764            double newObjectiveValue = solver->getObjSense() * solver->getObjValue();
4765            // See if integer solution
4766            int numInf;
4767            int numInf2;
4768            bool gotSol = model->feasibleSolution(numInf, numInf2);
4769            if (!gotSol) {
4770              COIN_DETAIL_PRINT(printf("numinf %d\n", numInf));
4771                double * sol = simplex->primalColumnSolution();
4772                for (int i = 0; i < numberColumns; i++) {
4773                    if (simplex->isInteger(i)) {
4774                        double value = floor(sol[i] + 0.5);
4775                        if (fabs(value - sol[i]) > 1.0e-7) {
4776                          COIN_DETAIL_PRINT(printf("%d value %g\n", i, sol[i]));
4777                            if (fabs(value - sol[i]) < 1.0e-3) {
4778                                sol[i] = value;
4779                            }
4780                        }
4781                    }
4782                }
4783                simplex->writeMps("bad8.mps", 2);
4784                bool gotSol = model->feasibleSolution(numInf, numInf2);
4785                if (!gotSol)
4786                    assert (gotSol);
4787            }
4788            model->setBestSolution(CBC_STRONGSOL,
4789                                   newObjectiveValue,
4790                                   solver->getColSolution()) ;
4791            model->setLastHeuristic(NULL);
4792            model->incrementUsed(solver->getColSolution());
4793        }
4794    }
4795    // restore bounds
4796    {
4797        for (int j = 0; j < numberColumns; j++) {
4798            if (saveLower[j] != lower[j])
4799                solver->setColLower(j, saveLower[j]);
4800            if (saveUpper[j] != upper[j])
4801                solver->setColUpper(j, saveUpper[j]);
4802        }
4803    }
4804    // restore basis
4805    solver->setWarmStart(ws);
4806    delete ws;
4807    int anyAction;
4808    //#define CHECK_PATH
4809#ifdef CHECK_PATH
4810    extern int gotGoodNode_Z;
4811    if (gotGoodNode_Z >= 0)
4812        printf("good node %d %g\n", gotGoodNode_Z, infeasibility);
4813#endif
4814    if (infeasibility > 0.0) {
4815        if (infeasibility == COIN_DBL_MAX) {
4816            anyAction = -2; // infeasible
4817        } else {
4818            branch_ = thisOne->createCbcBranch(solver, &usefulInfo, preferredWay);
4819            if (branch_) {
4820              // Set to first one (and change when re-pushing)
4821              CbcGeneralBranchingObject * branch =
4822                dynamic_cast <CbcGeneralBranchingObject *> (branch_);
4823              branch->state(objectiveValue_, sumInfeasibilities_,
4824                            numberUnsatisfied_, 0);
4825              branch->setNode(this);
4826              anyAction = 0;
4827            } else {
4828              anyAction = -2; // mark as infeasible
4829            }
4830        }
4831    } else {
4832        anyAction = -1;
4833    }
4834#ifdef CHECK_PATH
4835    gotGoodNode_Z = -1;
4836#endif
4837    // Set guessed solution value
4838    guessedObjectiveValue_ = objectiveValue_ + 1.0e-5;
4839    delete [] saveLower;
4840    delete [] saveUpper;
4841
4842    // restore solution
4843    solver->setColSolution(saveSolution);
4844    delete [] saveSolution;
4845    return anyAction;
4846}
4847/* Double checks in case node can change its mind!
4848   Returns objective value
4849   Can change objective etc */
4850double
4851CbcNode::checkIsCutoff(double cutoff)
4852{
4853    branch_->checkIsCutoff(cutoff);
4854    return objectiveValue_;
4855}
4856
Note: See TracBrowser for help on using the repository browser.