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

Last change on this file since 1826 was 1802, checked in by forrest, 7 years ago

add Proximity heuristic (Fischetti and Monaci) - shouldn't break anything

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