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

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

first attempt at cleaning up threads

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