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

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

fix bad bound message

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