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

Last change on this file since 1514 was 1514, checked in by forrest, 8 years ago

fix minor bugs

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