source: branches/sandbox/Cbc/src/CbcNode.cpp @ 1315

Last change on this file since 1315 was 1315, checked in by forrest, 10 years ago

final changes before cleaning

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