source: stable/2.7/Cbc/src/CbcNode.cpp @ 1577

Last change on this file since 1577 was 1573, checked in by lou, 9 years ago

Change to EPL license notice.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 197.6 KB
Line 
1/* $Id: CbcNode.cpp 1573 2011-01-05 01:12:36Z lou $ */
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 = ( CoinCpuTime() - model->getDblParam(CbcModel::CbcStartSeconds) >
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 = ( CoinCpuTime() - model->getDblParam(CbcModel::CbcStartSeconds) >
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 = ( CoinCpuTime() - model->getDblParam(CbcModel::CbcStartSeconds) >
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 TIGHTEN_BOUNDS
2864#ifdef 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 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()) {
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()) {
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 = ( CoinCpuTime() - model->getDblParam(CbcModel::CbcStartSeconds) >
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                        // See if candidate still possible
3386                        if (branch_) {
3387                            const OsiObject * object = model->object(bestChoice);
3388                            int preferredWay;
3389                            double infeasibility = object->infeasibility(&usefulInfo, preferredWay);
3390                            if (!infeasibility) {
3391                                // take out
3392                                delete branch_;
3393                                branch_ = NULL;
3394                            } else {
3395                                CbcBranchingObject * branchObj =
3396                                    dynamic_cast <CbcBranchingObject *>(branch_) ;
3397                                assert (branchObj);
3398                                branchObj->way(preferredWay);
3399                            }
3400                        }
3401                    } else {
3402                        anyAction = -2;
3403                        finished = true;
3404                    }
3405                }
3406            }
3407            // If  fixed then round again
3408            if (!branch_ && anyAction != -2) {
3409                finished = false;
3410            }
3411            delete ws;
3412        }
3413    }
3414    // update number of strong iterations etc
3415    model->incrementStrongInfo(numberStrongDone, numberStrongIterations,
3416                               anyAction == -2 ? 0 : numberToFix, anyAction == -2);
3417    if (model->searchStrategy() == -1) {
3418#ifndef COIN_DEVELOP
3419        if (solver->messageHandler()->logLevel() > 1)
3420#endif
3421            printf("%d strong, %d iters, %d inf, %d not finished, %d not trusted\n",
3422                   numberStrongDone, numberStrongIterations, numberStrongInfeasible, numberUnfinished,
3423                   numberNotTrusted);
3424        // decide what to do
3425        int strategy = 1;
3426        if (((numberUnfinished*4 > numberStrongDone &&
3427                numberStrongInfeasible*40 < numberStrongDone) ||
3428                numberStrongInfeasible < 0) && model->numberStrong() < 10 && model->numberBeforeTrust() <= 20 && model->numberObjects() > CoinMax(1000, solver->getNumRows())) {
3429            strategy = 2;
3430#ifdef COIN_DEVELOP
3431            //if (model->logLevel()>1)
3432            printf("going to strategy 2\n");
3433#endif
3434            // Weaken
3435            model->setNumberStrong(2);
3436            model->setNumberBeforeTrust(1);
3437            model->synchronizeNumberBeforeTrust();
3438        }
3439        if (numberNodes)
3440            strategy = 1;  // should only happen after hot start
3441        model->setSearchStrategy(strategy);
3442    } else if (numberStrongDone) {
3443        //printf("%d strongB, %d iters, %d inf, %d not finished, %d not trusted\n",
3444        //   numberStrongDone,numberStrongIterations,numberStrongInfeasible,numberUnfinished,
3445        //   numberNotTrusted);
3446    }
3447    if (model->searchStrategy() == 1 && numberNodes > 500 && numberNodes < -510) {
3448#ifndef COIN_DEVELOP
3449        if (solver->messageHandler()->logLevel() > 1)
3450#endif
3451            printf("after %d nodes - %d strong, %d iters, %d inf, %d not finished, %d not trusted\n",
3452                   numberNodes, numberStrongDone, numberStrongIterations, numberStrongInfeasible, numberUnfinished,
3453                   numberNotTrusted);
3454        // decide what to do
3455        if (numberUnfinished*10 > numberStrongDone + 1 ||
3456                !numberStrongInfeasible) {
3457            printf("going to strategy 2\n");
3458            // Weaken
3459            model->setNumberStrong(2);
3460            model->setNumberBeforeTrust(1);
3461            model->synchronizeNumberBeforeTrust();
3462            model->setSearchStrategy(2);
3463        }
3464    }
3465    if (numberUnfinished*10 < numberStrongDone &&
3466            numberStrongIterations*20 < model->getIterationCount()) {
3467        //printf("increasing trust\n");
3468        model->synchronizeNumberBeforeTrust(2);
3469    }
3470
3471    // Set guessed solution value
3472    guessedObjectiveValue_ = objectiveValue_ + estimatedDegradation;
3473
3474    /*
3475      Cleanup, then we're finished
3476    */
3477    if (!model->branchingMethod())
3478        delete decision;
3479
3480    delete choiceObject;
3481    delete [] sort;
3482    delete [] whichObject;
3483#ifdef RANGING
3484    delete [] objectMark;
3485#endif
3486    delete [] saveLower;
3487    delete [] saveUpper;
3488    delete [] upEstimate;
3489    delete [] downEstimate;
3490# ifdef COIN_HAS_CLP
3491    if (osiclp) {
3492        osiclp->setSpecialOptions(saveClpOptions);
3493    }
3494# endif
3495    // restore solution
3496    solver->setColSolution(saveSolution);
3497    model->reserveCurrentSolution(saveSolution);
3498    delete [] saveSolution;
3499    model->setStateOfSearch(saveStateOfSearch);
3500    model->setLogLevel(saveLogLevel);
3501    // delete extra regions
3502    if (usefulInfo.usefulRegion_) {
3503        delete [] usefulInfo.usefulRegion_;
3504        delete [] usefulInfo.indexRegion_;
3505        delete [] usefulInfo.pi_;
3506        usefulInfo.usefulRegion_ = NULL;
3507        usefulInfo.indexRegion_ = NULL;
3508        usefulInfo.pi_ = NULL;
3509    }
3510    useShadow = model->moreSpecialOptions() & 7;
3511    if ((useShadow == 5 && model->getSolutionCount()) || useShadow == 6) {
3512        // zap pseudo shadow prices
3513        model->pseudoShadow(-1);
3514        // and switch off
3515        model->setMoreSpecialOptions(model->moreSpecialOptions()&(~1023));
3516    }
3517    return anyAction;
3518}
3519int CbcNode::analyze (CbcModel *model, double * results)
3520{
3521    int i;
3522    int numberIterationsAllowed = model->numberAnalyzeIterations();
3523    OsiSolverInterface * solver = model->solver();
3524    objectiveValue_ = solver->getObjSense() * solver->getObjValue();
3525    double cutoff = model->getCutoff();
3526    const double * lower = solver->getColLower();
3527    const double * upper = solver->getColUpper();
3528    const double * dj = solver->getReducedCost();
3529    int numberObjects = model->numberObjects();
3530    int numberColumns = model->getNumCols();
3531    // Initialize arrays
3532    int numberIntegers = model->numberIntegers();
3533    int * back = new int[numberColumns];
3534    const int * integerVariable = model->integerVariable();
3535    for (i = 0; i < numberColumns; i++)
3536        back[i] = -1;
3537    // What results is
3538    double * newLower = results;
3539    double * objLower = newLower + numberIntegers;
3540    double * newUpper = objLower + numberIntegers;
3541    double * objUpper = newUpper + numberIntegers;
3542    for (i = 0; i < numberIntegers; i++) {
3543        int iColumn = integerVariable[i];
3544        back[iColumn] = i;
3545        newLower[i] = 0.0;
3546        objLower[i] = -COIN_DBL_MAX;
3547        newUpper[i] = 0.0;
3548        objUpper[i] = -COIN_DBL_MAX;
3549    }
3550    double * saveUpper = new double[numberColumns];
3551    double * saveLower = new double[numberColumns];
3552    int anyAction = 0;
3553    // Save solution in case heuristics need good solution later
3554
3555    double * saveSolution = new double[numberColumns];
3556    memcpy(saveSolution, solver->getColSolution(), numberColumns*sizeof(double));
3557    model->reserveCurrentSolution(saveSolution);
3558    for (i = 0; i < numberColumns; i++) {
3559        saveLower[i] = lower[i];
3560        saveUpper[i] = upper[i];
3561    }
3562    // Get arrays to sort
3563    double * sort = new double[numberObjects];
3564    int * whichObject = new int[numberObjects];
3565    int numberToFix = 0;
3566    int numberToDo = 0;
3567    double integerTolerance =
3568        model->getDblParam(CbcModel::CbcIntegerTolerance);
3569    // point to useful information
3570    OsiBranchingInformation usefulInfo = model->usefulInformation();
3571    // and modify
3572    usefulInfo.depth_ = depth_;
3573
3574    // compute current state
3575    int numberObjectInfeasibilities; // just odd ones
3576    int numberIntegerInfeasibilities;
3577    model->feasibleSolution(
3578        numberIntegerInfeasibilities,
3579        numberObjectInfeasibilities);
3580# ifdef COIN_HAS_CLP
3581    OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
3582    int saveClpOptions = 0;
3583    bool fastIterations = (model->specialOptions() & 8) != 0;
3584    if (osiclp && fastIterations) {
3585        // for faster hot start
3586        saveClpOptions = osiclp->specialOptions();
3587        osiclp->setSpecialOptions(saveClpOptions | 8192);
3588    }
3589# else
3590    bool fastIterations = false ;
3591# endif
3592    /*
3593      Scan for branching objects that indicate infeasibility. Choose candidates
3594      using priority as the first criteria, then integer infeasibility.
3595
3596      The algorithm is to fill the array with a set of good candidates (by
3597      infeasibility) with priority bestPriority.  Finding a candidate with
3598      priority better (less) than bestPriority flushes the choice array. (This
3599      serves as initialization when the first candidate is found.)
3600
3601    */
3602    numberToDo = 0;
3603    for (i = 0; i < numberObjects; i++) {
3604        OsiObject * object = model->modifiableObject(i);
3605        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3606            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3607        if (!dynamicObject)
3608            continue;
3609        int preferredWay;
3610        double infeasibility = object->infeasibility(&usefulInfo, preferredWay);
3611        int iColumn = dynamicObject->columnNumber();
3612        if (saveUpper[iColumn] == saveLower[iColumn])
3613            continue;
3614        if (infeasibility)
3615            sort[numberToDo] = -1.0e10 - infeasibility;
3616        else
3617            sort[numberToDo] = -fabs(dj[iColumn]);
3618        whichObject[numberToDo++] = i;
3619    }
3620    // Save basis
3621    CoinWarmStart * ws = solver->getWarmStart();
3622    int saveLimit;
3623    solver->getIntParam(OsiMaxNumIterationHotStart, saveLimit);
3624    int targetIterations = CoinMax(500, numberIterationsAllowed / numberObjects);
3625    if (saveLimit < targetIterations)
3626        solver->setIntParam(OsiMaxNumIterationHotStart, targetIterations);
3627    // Mark hot start
3628    solver->markHotStart();
3629    // Sort
3630    CoinSort_2(sort, sort + numberToDo, whichObject);
3631    double * currentSolution = model->currentSolution();
3632    double objMin = 1.0e50;
3633    double objMax = -1.0e50;
3634    bool needResolve = false;
3635    /*
3636      Now calculate the cost forcing the variable up and down.
3637    */
3638    int iDo;
3639    for (iDo = 0; iDo < numberToDo; iDo++) {
3640        CbcStrongInfo choice;
3641        int iObject = whichObject[iDo];
3642        OsiObject * object = model->modifiableObject(iObject);
3643        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3644            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3645        if (!dynamicObject)
3646            continue;
3647        int iColumn = dynamicObject->columnNumber();
3648        int preferredWay;
3649        /*
3650          Update the information held in the object.
3651        */
3652        object->infeasibility(&usefulInfo, preferredWay);
3653        double value = currentSolution[iColumn];
3654        double nearest = floor(value + 0.5);
3655        double lowerValue = floor(value);
3656        bool satisfied = false;
3657        if (fabs(value - nearest) <= integerTolerance || value < saveLower[iColumn] || value > saveUpper[iColumn]) {
3658            satisfied = true;
3659            double newValue;
3660            if (nearest < saveUpper[iColumn]) {
3661                newValue = nearest + 1.0001 * integerTolerance;
3662                lowerValue = nearest;
3663            } else {
3664                newValue = nearest - 1.0001 * integerTolerance;
3665                lowerValue = nearest - 1;
3666            }
3667            currentSolution[iColumn] = newValue;
3668        }
3669        double upperValue = lowerValue + 1.0;
3670        //CbcSimpleInteger * obj =
3671        //dynamic_cast <CbcSimpleInteger *>(object) ;
3672        //if (obj) {
3673        //choice.possibleBranch=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
3674        //} else {
3675        CbcObject * obj =
3676            dynamic_cast <CbcObject *>(object) ;
3677        assert (obj);
3678        choice.possibleBranch = obj->createCbcBranch(solver, &usefulInfo, preferredWay);
3679        //}
3680        currentSolution[iColumn] = value;
3681        // Save which object it was
3682        choice.objectNumber = iObject;
3683        choice.numIntInfeasUp = numberUnsatisfied_;
3684        choice.numIntInfeasDown = numberUnsatisfied_;
3685        choice.downMovement = 0.0;
3686        choice.upMovement = 0.0;
3687        choice.numItersDown = 0;
3688        choice.numItersUp = 0;
3689        choice.fix = 0; // say not fixed
3690        double objectiveChange ;
3691        double newObjectiveValue = 1.0e100;
3692        int j;
3693        // status is 0 finished, 1 infeasible and other
3694        int iStatus;
3695        /*
3696          Try the down direction first. (Specify the initial branching alternative as
3697          down with a call to way(-1). Each subsequent call to branch() performs the
3698          specified branch and advances the branch object state to the next branch
3699          alternative.)
3700        */
3701        choice.possibleBranch->way(-1) ;
3702        choice.possibleBranch->branch() ;
3703        if (fabs(value - lowerValue) > integerTolerance) {
3704            solver->solveFromHotStart() ;
3705            /*
3706              We now have an estimate of objective degradation that we can use for strong
3707              branching. If we're over the cutoff, the variable is monotone up.
3708              If we actually made it to optimality, check for a solution, and if we have
3709              a good one, call setBestSolution to process it. Note that this may reduce the
3710              cutoff, so we check again to see if we can declare this variable monotone.
3711            */
3712            if (solver->isProvenOptimal())
3713                iStatus = 0; // optimal
3714            else if (solver->isIterationLimitReached()
3715                     && !solver->isDualObjectiveLimitReached())
3716                iStatus = 2; // unknown
3717            else
3718                iStatus = 1; // infeasible
3719            newObjectiveValue = solver->getObjSense() * solver->getObjValue();
3720            choice.numItersDown = solver->getIterationCount();
3721            numberIterationsAllowed -= choice.numItersDown;
3722            objectiveChange = newObjectiveValue  - objectiveValue_;
3723            if (!iStatus) {
3724                choice.finishedDown = true ;
3725                if (newObjectiveValue >= cutoff) {
3726                    objectiveChange = 1.0e100; // say infeasible
3727                } else {
3728                    // See if integer solution
3729                    if (model->feasibleSolution(choice.numIntInfeasDown,
3730                                                choice.numObjInfeasDown)
3731                            && model->problemFeasibility()->feasible(model, -1) >= 0) {
3732                        model->setBestSolution(CBC_STRONGSOL,
3733                                               newObjectiveValue,
3734                                               solver->getColSolution()) ;
3735                        model->setLastHeuristic(NULL);
3736                        model->incrementUsed(solver->getColSolution());
3737                        cutoff = model->getCutoff();
3738                        if (newObjectiveValue >= cutoff)        //  *new* cutoff
3739                            objectiveChange = 1.0e100 ;
3740                    }
3741                }
3742            } else if (iStatus == 1) {
3743                objectiveChange = 1.0e100 ;
3744            } else {
3745                // Can't say much as we did not finish
3746                choice.finishedDown = false ;
3747            }
3748            choice.downMovement = objectiveChange ;
3749        }
3750        // restore bounds
3751        for ( j = 0; j < numberColumns; j++) {
3752            if (saveLower[j] != lower[j])
3753                solver->setColLower(j, saveLower[j]);
3754            if (saveUpper[j] != upper[j])
3755                solver->setColUpper(j, saveUpper[j]);
3756        }
3757        // repeat the whole exercise, forcing the variable up
3758        choice.possibleBranch->branch();
3759        if (fabs(value - upperValue) > integerTolerance) {
3760            solver->solveFromHotStart() ;
3761            /*
3762              We now have an estimate of objective degradation that we can use for strong
3763              branching. If we're over the cutoff, the variable is monotone up.
3764              If we actually made it to optimality, check for a solution, and if we have
3765              a good one, call setBestSolution to process it. Note that this may reduce the
3766              cutoff, so we check again to see if we can declare this variable monotone.
3767            */
3768            if (solver->isProvenOptimal())
3769                iStatus = 0; // optimal
3770            else if (solver->isIterationLimitReached()
3771                     && !solver->isDualObjectiveLimitReached())
3772                iStatus = 2; // unknown
3773            else
3774                iStatus = 1; // infeasible
3775            newObjectiveValue = solver->getObjSense() * solver->getObjValue();
3776            choice.numItersUp = solver->getIterationCount();
3777            numberIterationsAllowed -= choice.numItersUp;
3778            objectiveChange = newObjectiveValue  - objectiveValue_;
3779            if (!iStatus) {
3780                choice.finishedUp = true ;
3781                if (newObjectiveValue >= cutoff) {
3782                    objectiveChange = 1.0e100; // say infeasible
3783                } else {
3784                    // See if integer solution
3785                    if (model->feasibleSolution(choice.numIntInfeasUp,
3786                                                choice.numObjInfeasUp)
3787                            && model->problemFeasibility()->feasible(model, -1) >= 0) {
3788                        model->setBestSolution(CBC_STRONGSOL,
3789                                               newObjectiveValue,
3790                                               solver->getColSolution()) ;
3791                        model->setLastHeuristic(NULL);
3792                        model->incrementUsed(solver->getColSolution());
3793                        cutoff = model->getCutoff();
3794                        if (newObjectiveValue >= cutoff)        //  *new* cutoff
3795                            objectiveChange = 1.0e100 ;
3796                    }
3797                }
3798            } else if (iStatus == 1) {
3799                objectiveChange = 1.0e100 ;
3800            } else {
3801                // Can't say much as we did not finish
3802                choice.finishedUp = false ;
3803            }
3804            choice.upMovement = objectiveChange ;
3805
3806            // restore bounds
3807            for ( j = 0; j < numberColumns; j++) {
3808                if (saveLower[j] != lower[j])
3809                    solver->setColLower(j, saveLower[j]);
3810                if (saveUpper[j] != upper[j])
3811                    solver->setColUpper(j, saveUpper[j]);
3812            }
3813        }
3814        // If objective goes above certain amount we can set bound
3815        int jInt = back[iColumn];
3816        newLower[jInt] = upperValue;
3817        if (choice.finishedDown)
3818            objLower[jInt] = choice.downMovement + objectiveValue_;
3819        else
3820            objLower[jInt] = objectiveValue_;
3821        newUpper[jInt] = lowerValue;
3822        if (choice.finishedUp)
3823            objUpper[jInt] = choice.upMovement + objectiveValue_;
3824        else
3825            objUpper[jInt] = objectiveValue_;
3826        objMin = CoinMin(CoinMin(objLower[jInt], objUpper[jInt]), objMin);
3827        /*
3828          End of evaluation for this candidate variable. Possibilities are:
3829          * Both sides below cutoff; this variable is a candidate for branching.
3830          * Both sides infeasible or above the objective cutoff: no further action
3831          here. Break from the evaluation loop and assume the node will be purged
3832          by the caller.
3833          * One side below cutoff: Install the branch (i.e., fix the variable). Break
3834          from the evaluation loop and assume the node will be reoptimised by the
3835          caller.
3836        */
3837        if (choice.upMovement < 1.0e100) {
3838            if (choice.downMovement < 1.0e100) {
3839                objMax = CoinMax(CoinMax(objLower[jInt], objUpper[jInt]), objMax);
3840                // In case solution coming in was odd
3841                choice.upMovement = CoinMax(0.0, choice.upMovement);
3842                choice.downMovement = CoinMax(0.0, choice.downMovement);
3843                // feasible -
3844                model->messageHandler()->message(CBC_STRONG, *model->messagesPointer())
3845                << iObject << iColumn
3846                << choice.downMovement << choice.numIntInfeasDown
3847                << choice.upMovement << choice.numIntInfeasUp
3848                << value
3849                << CoinMessageEol;
3850            } else {
3851                // up feasible, down infeasible
3852                anyAction = -1;
3853                if (!satisfied)
3854                    needResolve = true;
3855                choice.fix = 1;
3856                numberToFix++;
3857                saveLower[iColumn] = upperValue;
3858                solver->setColLower(iColumn, upperValue);
3859            }
3860        } else {
3861            if (choice.downMovement < 1.0e100) {
3862                // down feasible, up infeasible
3863                anyAction = -1;
3864                if (!satisfied)
3865                    needResolve = true;
3866                choice.fix = -1;
3867                numberToFix++;
3868                saveUpper[iColumn] = lowerValue;
3869                solver->setColUpper(iColumn, lowerValue);
3870            } else {
3871                // neither side feasible
3872                anyAction = -2;
3873                printf("Both infeasible for choice %d sequence %d\n", i,
3874                       model->object(choice.objectNumber)->columnNumber());
3875                delete ws;
3876                ws = NULL;
3877                //solver->writeMps("bad");
3878                numberToFix = -1;
3879                delete choice.possibleBranch;
3880                choice.possibleBranch = NULL;
3881                break;
3882            }
3883        }
3884        delete choice.possibleBranch;
3885        if (numberIterationsAllowed <= 0)
3886            break;
3887        //printf("obj %d, col %d, down %g up %g value %g\n",iObject,iColumn,
3888        //     choice.downMovement,choice.upMovement,value);
3889    }
3890    printf("Best possible solution %g, can fix more if solution of %g found - looked at %d variables in %d iterations\n",
3891           objMin, objMax, iDo, model->numberAnalyzeIterations() - numberIterationsAllowed);
3892    model->setNumberAnalyzeIterations(numberIterationsAllowed);
3893    // Delete the snapshot
3894    solver->unmarkHotStart();
3895    // back to normal
3896    solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
3897    solver->setIntParam(OsiMaxNumIterationHotStart, saveLimit);
3898    // restore basis
3899    solver->setWarmStart(ws);
3900    delete ws;
3901
3902    delete [] sort;
3903    delete [] whichObject;
3904    delete [] saveLower;
3905    delete [] saveUpper;
3906    delete [] back;
3907    // restore solution
3908    solver->setColSolution(saveSolution);
3909# ifdef COIN_HAS_CLP
3910    if (osiclp)
3911        osiclp->setSpecialOptions(saveClpOptions);
3912# endif
3913    model->reserveCurrentSolution(saveSolution);
3914    delete [] saveSolution;
3915    if (needResolve)
3916        solver->resolve();
3917    return numberToFix;
3918}
3919
3920
3921CbcNode::CbcNode(const CbcNode & rhs)
3922        : CoinTreeNode(rhs)
3923{
3924#ifdef CHECK_NODE
3925    printf("CbcNode %x Constructor from rhs %x\n", this, &rhs);
3926#endif
3927    if (rhs.nodeInfo_)
3928        nodeInfo_ = rhs.nodeInfo_->clone();
3929    else
3930        nodeInfo_ = NULL;
3931    objectiveValue_ = rhs.objectiveValue_;
3932    guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
3933    sumInfeasibilities_ = rhs.sumInfeasibilities_;
3934    if (rhs.branch_)
3935        branch_ = rhs.branch_->clone();
3936    else
3937        branch_ = NULL;
3938    depth_ = rhs.depth_;
3939    numberUnsatisfied_ = rhs.numberUnsatisfied_;
3940    nodeNumber_ = rhs.nodeNumber_;
3941    state_ = rhs.state_;
3942    if (nodeInfo_)
3943        assert ((state_&2) != 0);
3944    else
3945        assert ((state_&2) == 0);
3946}
3947
3948CbcNode &
3949CbcNode::operator=(const CbcNode & rhs)
3950{
3951    if (this != &rhs) {
3952        delete nodeInfo_;
3953        if (rhs.nodeInfo_)
3954            nodeInfo_ = rhs.nodeInfo_->clone();
3955        else
3956            nodeInfo_ = NULL;
3957        objectiveValue_ = rhs.objectiveValue_;
3958        guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
3959        sumInfeasibilities_ = rhs.sumInfeasibilities_;
3960        if (rhs.branch_)
3961            branch_ = rhs.branch_->clone();
3962        else
3963            branch_ = NULL,
3964                      depth_ = rhs.depth_;
3965        numberUnsatisfied_ = rhs.numberUnsatisfied_;
3966        nodeNumber_ = rhs.nodeNumber_;
3967        state_ = rhs.state_;
3968        if (nodeInfo_)
3969            assert ((state_&2) != 0);
3970        else
3971            assert ((state_&2) == 0);
3972    }
3973    return *this;
3974}
3975CbcNode::~CbcNode ()
3976{
3977#ifdef CHECK_NODE
3978    if (nodeInfo_) {
3979        printf("CbcNode %x Destructor nodeInfo %x (%d)\n",
3980               this, nodeInfo_, nodeInfo_->numberPointingToThis());
3981        //assert(nodeInfo_->numberPointingToThis()>=0);
3982    } else {
3983        printf("CbcNode %x Destructor nodeInfo %x (?)\n",
3984               this, nodeInfo_);
3985    }
3986#endif
3987    if (nodeInfo_) {
3988        // was if (nodeInfo_&&(state_&2)!=0) {
3989        nodeInfo_->nullOwner();
3990        int numberToDelete = nodeInfo_->numberBranchesLeft();
3991        //    CbcNodeInfo * parent = nodeInfo_->parent();
3992        //assert (nodeInfo_->numberPointingToThis()>0);
3993        if (nodeInfo_->decrement(numberToDelete) == 0 || (state_&2) == 0) {
3994            if ((state_&2) == 0)
3995                nodeInfo_->nullParent();
3996            delete nodeInfo_;
3997        } else {
3998            //printf("node %x nodeinfo %x parent %x\n",this,nodeInfo_,nodeInfo_->parent());
3999            // anyway decrement parent
4000            //if (parent)
4001            ///parent->decrement(1);
4002        }
4003    }
4004    delete branch_;
4005}
4006// Decrement  active cut counts
4007void
4008CbcNode::decrementCuts(int change)
4009{
4010    if (nodeInfo_)
4011        assert ((state_&2) != 0);
4012    else
4013        assert ((state_&2) == 0);
4014    if (nodeInfo_) {
4015        nodeInfo_->decrementCuts(change);
4016    }
4017}
4018void
4019CbcNode::decrementParentCuts(CbcModel * model, int change)
4020{
4021    if (nodeInfo_)
4022        assert ((state_&2) != 0);
4023    else
4024        assert ((state_&2) == 0);
4025    if (nodeInfo_) {
4026        nodeInfo_->decrementParentCuts(model, change);
4027    }
4028}
4029
4030/*
4031  Initialize reference counts (numberPointingToThis, numberBranchesLeft_)
4032  in the attached nodeInfo_.
4033*/
4034void
4035CbcNode::initializeInfo()
4036{
4037    assert(nodeInfo_ && branch_) ;
4038    nodeInfo_->initializeInfo(branch_->numberBranches());
4039    assert ((state_&2) != 0);
4040    assert (nodeInfo_->numberBranchesLeft() ==
4041            branch_->numberBranchesLeft());
4042}
4043// Nulls out node info
4044void
4045CbcNode::nullNodeInfo()
4046{
4047    nodeInfo_ = NULL;
4048    // say not active
4049    state_ &= ~2;
4050}
4051
4052int
4053CbcNode::branch(OsiSolverInterface * solver)
4054{
4055    double changeInGuessed;
4056    assert (nodeInfo_->numberBranchesLeft() ==
4057            branch_->numberBranchesLeft());
4058    if (!solver)
4059        changeInGuessed = branch_->branch();
4060    else
4061        changeInGuessed = branch_->branch(solver);
4062    guessedObjectiveValue_ += changeInGuessed;
4063    //#define PRINTIT
4064#ifdef PRINTIT
4065    int numberLeft = nodeInfo_->numberBranchesLeft();
4066    CbcNodeInfo * parent = nodeInfo_->parent();
4067    int parentNodeNumber = -1;
4068    CbcBranchingObject * object1 =
4069        dynamic_cast<CbcBranchingObject *>(branch_) ;
4070    //OsiObject * object = object1->
4071    //int sequence = object->columnNumber);
4072    int id = -1;
4073    double value = 0.0;
4074    if (object1) {
4075        id = object1->variable();
4076        value = object1->value();
4077    }
4078    printf("id %d value %g objvalue %g\n", id, value, objectiveValue_);
4079    if (parent)
4080        parentNodeNumber = parent->nodeNumber();
4081    printf("Node number %d, %s, way %d, depth %d, parent node number %d\n",
4082           nodeInfo_->nodeNumber(), (numberLeft == 2) ? "leftBranch" : "rightBranch",
4083           way(), depth_, parentNodeNumber);
4084    assert (parentNodeNumber != nodeInfo_->nodeNumber());
4085#endif
4086    return nodeInfo_->branchedOn();
4087}
4088/* Active arm of the attached OsiBranchingObject.
4089
4090   In the simplest instance, coded -1 for the down arm of the branch, +1 for
4091   the up arm. But see OsiBranchingObject::way()
4092   Use nodeInfo--.numberBranchesLeft_ to see how active
4093
4094   Except that there is no OsiBranchingObject::way(), and this'll fail in any
4095   event because we have various OsiXXXBranchingObjects which aren't descended
4096   from CbcBranchingObjects. I think branchIndex() is the appropriate
4097   equivalent, but could be wrong. (lh, 061220)
4098
4099   071212: I'm finally getting back to cbc-generic and rescuing a lot of my
4100   annotation from branches/devel (which was killed in summer). I'm going to
4101   put back an assert(obj) just to see what happens. It's still present as of
4102   the most recent change to CbcNode (r833).
4103
4104   080104: Yep, we can arrive here with an OsiBranchingObject. Removed the
4105   assert, it's served its purpose.
4106
4107   080226: John finally noticed this problem and added a way() method to the
4108   OsiBranchingObject hierarchy. Removing my workaround.
4109
4110*/
4111int
4112CbcNode::way() const
4113{
4114    if (branch_) {
4115        CbcBranchingObject * obj =
4116            dynamic_cast <CbcBranchingObject *>(branch_) ;
4117        if (obj) {
4118            return obj->way();
4119        } else {
4120            OsiTwoWayBranchingObject * obj2 =
4121                dynamic_cast <OsiTwoWayBranchingObject *>(branch_) ;
4122            assert (obj2);
4123            return obj2->way();
4124        }
4125    } else {
4126        return 0;
4127    }
4128}
4129/* Create a branching object for the node
4130
4131   The routine scans the object list of the model and selects a set of
4132   unsatisfied objects as candidates for branching. The candidates are
4133   evaluated, and an appropriate branch object is installed.
4134
4135   The numberPassesLeft is decremented to stop fixing one variable each time
4136   and going on and on (e.g. for stock cutting, air crew scheduling)
4137
4138   If evaluation determines that an object is monotone or infeasible,
4139   the routine returns immediately. In the case of a monotone object,
4140   the branch object has already been called to modify the model.
4141
4142   Return value:
4143   <ul>
4144   <li>  0: A branching object has been installed
4145   <li> -1: A monotone object was discovered
4146   <li> -2: An infeasible object was discovered
4147   </ul>
4148   Branch state:
4149   <ul>
4150   <li> -1: start
4151   <li> -1: A monotone object was discovered
4152   <li> -2: An infeasible object was discovered
4153   </ul>
4154*/
4155int
4156CbcNode::chooseOsiBranch (CbcModel * model,
4157                          CbcNode * lastNode,
4158                          OsiBranchingInformation * usefulInfo,
4159                          int branchState)
4160{
4161    int returnStatus = 0;
4162    if (lastNode)
4163        depth_ = lastNode->depth_ + 1;
4164    else
4165        depth_ = 0;
4166    OsiSolverInterface * solver = model->solver();
4167    objectiveValue_ = solver->getObjValue() * solver->getObjSense();
4168    usefulInfo->objectiveValue_ = objectiveValue_;
4169    usefulInfo->depth_ = depth_;
4170    const double * saveInfoSol = usefulInfo->solution_;
4171    double * saveSolution = new double[solver->getNumCols()];
4172    memcpy(saveSolution, solver->getColSolution(), solver->getNumCols()*sizeof(double));
4173    usefulInfo->solution_ = saveSolution;
4174    OsiChooseVariable * choose = model->branchingMethod()->chooseMethod();
4175    int numberUnsatisfied = -1;
4176    if (branchState < 0) {
4177        // initialize
4178        // initialize sum of "infeasibilities"
4179        sumInfeasibilities_ = 0.0;
4180        numberUnsatisfied = choose->setupList(usefulInfo, true);
4181        numberUnsatisfied_ = numberUnsatisfied;
4182        branchState = 0;
4183        if (numberUnsatisfied_ < 0) {
4184            // infeasible
4185            delete [] saveSolution;
4186            return -2;
4187        }
4188    }
4189    // unset best
4190    int best = -1;
4191    choose->setBestObjectIndex(-1);
4192    if (numberUnsatisfied) {
4193        if (branchState > 0 || !choose->numberOnList()) {
4194            // we need to return at once - don't do strong branching or anything
4195            if (choose->numberOnList() || !choose->numberStrong()) {
4196                best = choose->candidates()[0];
4197                choose->setBestObjectIndex(best);
4198            } else {
4199                // nothing on list - need to try again - keep any solution
4200                numberUnsatisfied = choose->setupList(usefulInfo, false);
4201                numberUnsatisfied_ = numberUnsatisfied;
4202                if (numberUnsatisfied) {
4203                    best = choose->candidates()[0];
4204                    choose->setBestObjectIndex(best);
4205                }
4206            }
4207        } else {
4208            // carry on with strong branching or whatever
4209            int returnCode = choose->chooseVariable(solver, usefulInfo, true);
4210            // update number of strong iterations etc
4211            model->incrementStrongInfo(choose->numberStrongDone(), choose->numberStrongIterations(),
4212                                       returnCode == -1 ? 0 : choose->numberStrongFixed(), returnCode == -1);
4213            if (returnCode > 1) {
4214                // has fixed some
4215                returnStatus = -1;
4216            } else if (returnCode == -1) {
4217                // infeasible
4218                returnStatus = -2;
4219            } else if (returnCode == 0) {
4220                // normal
4221                returnStatus = 0;
4222                numberUnsatisfied = 1;
4223            } else {
4224                // ones on list satisfied - double check
4225                numberUnsatisfied = choose->setupList(usefulInfo, false);
4226                numberUnsatisfied_ = numberUnsatisfied;
4227                if (numberUnsatisfied) {
4228                    best = choose->candidates()[0];
4229                    choose->setBestObjectIndex(best);
4230                }
4231            }
4232        }
4233    }
4234    delete branch_;
4235    branch_ = NULL;
4236    guessedObjectiveValue_ = COIN_DBL_MAX;//objectiveValue_; // for now
4237    if (!returnStatus) {
4238        if (numberUnsatisfied) {
4239            // create branching object
4240            const OsiObject * obj = model->solver()->object(choose->bestObjectIndex());
4241            //const OsiSolverInterface * solver = usefulInfo->solver_;
4242            branch_ = obj->createBranch(model->solver(), usefulInfo, obj->whichWay());
4243        }
4244    }
4245    usefulInfo->solution_ = saveInfoSol;
4246    delete [] saveSolution;
4247    // may have got solution
4248    if (choose->goodSolution()
4249            && model->problemFeasibility()->feasible(model, -1) >= 0) {
4250        // yes
4251        double objValue = choose->goodObjectiveValue();
4252        model->setBestSolution(CBC_STRONGSOL,
4253                               objValue,
4254                               choose->goodSolution()) ;
4255        model->setLastHeuristic(NULL);
4256        model->incrementUsed(choose->goodSolution());
4257        choose->clearGoodSolution();
4258    }
4259    return returnStatus;
4260}
4261int
4262CbcNode::chooseClpBranch (CbcModel * model,
4263                          CbcNode * lastNode)
4264{
4265    assert(lastNode);
4266    depth_ = lastNode->depth_ + 1;
4267    delete branch_;
4268    branch_ = NULL;
4269    OsiSolverInterface * solver = model->solver();
4270    //double saveObjectiveValue = solver->getObjValue();
4271    //double objectiveValue = CoinMax(solver->getObjSense()*saveObjectiveValue,objectiveValue_);
4272    const double * lower = solver->getColLower();
4273    const double * upper = solver->getColUpper();
4274    // point to useful information
4275    OsiBranchingInformation usefulInfo = model->usefulInformation();
4276    // and modify
4277    usefulInfo.depth_ = depth_;
4278    int i;
4279    //bool beforeSolution = model->getSolutionCount()==0;
4280    int numberObjects = model->numberObjects();
4281    int numberColumns = model->getNumCols();
4282    double * saveUpper = new double[numberColumns];
4283    double * saveLower = new double[numberColumns];
4284
4285    // Save solution in case heuristics need good solution later
4286
4287    double * saveSolution = new double[numberColumns];
4288    memcpy(saveSolution, solver->getColSolution(), numberColumns*sizeof(double));
4289    model->reserveCurrentSolution(saveSolution);
4290    for (i = 0; i < numberColumns; i++) {
4291        saveLower[i] = lower[i];
4292        saveUpper[i] = upper[i];
4293    }
4294    // Save basis
4295    CoinWarmStart * ws = solver->getWarmStart();
4296    numberUnsatisfied_ = 0;
4297    // initialize sum of "infeasibilities"
4298    sumInfeasibilities_ = 0.0;
4299    // Note looks as if off end (hidden one)
4300    OsiObject * object = model->modifiableObject(numberObjects);
4301    CbcGeneralDepth * thisOne = dynamic_cast <CbcGeneralDepth *> (object);
4302    assert (thisOne);
4303    OsiClpSolverInterface * clpSolver
4304    = dynamic_cast<OsiClpSolverInterface *> (solver);
4305    assert (clpSolver);
4306    ClpSimplex * simplex = clpSolver->getModelPtr();
4307    int preferredWay;
4308    double infeasibility = object->infeasibility(&usefulInfo, preferredWay);
4309    if (thisOne->whichSolution() >= 0) {
4310        ClpNode * nodeInfo = thisOne->nodeInfo(thisOne->whichSolution());
4311        nodeInfo->applyNode(simplex, 2);
4312        int saveLogLevel = simplex->logLevel();
4313        simplex->setLogLevel(0);
4314        simplex->dual();
4315        simplex->setLogLevel(saveLogLevel);
4316        double cutoff = model->getCutoff();
4317        bool goodSolution = true;
4318        if (simplex->status()) {
4319            //simplex->writeMps("bad7.mps",2);
4320            if (nodeInfo->objectiveValue() > cutoff - 1.0e-2)
4321                goodSolution = false;
4322            else
4323                assert (!simplex->status());
4324        }
4325        if (goodSolution) {
4326            double newObjectiveValue = solver->getObjSense() * solver->getObjValue();
4327            // See if integer solution
4328            int numInf;
4329            int numInf2;
4330            bool gotSol = model->feasibleSolution(numInf, numInf2);
4331            if (!gotSol) {
4332                printf("numinf %d\n", numInf);
4333                double * sol = simplex->primalColumnSolution();
4334                for (int i = 0; i < numberColumns; i++) {
4335                    if (simplex->isInteger(i)) {
4336                        double value = floor(sol[i] + 0.5);
4337                        if (fabs(value - sol[i]) > 1.0e-7) {
4338                            printf("%d value %g\n", i, sol[i]);
4339                            if (fabs(value - sol[i]) < 1.0e-3) {
4340                                sol[i] = value;
4341                            }
4342                        }
4343                    }
4344                }
4345                simplex->writeMps("bad8.mps", 2);
4346                bool gotSol = model->feasibleSolution(numInf, numInf2);
4347                if (!gotSol)
4348                    assert (gotSol);
4349            }
4350            model->setBestSolution(CBC_STRONGSOL,
4351                                   newObjectiveValue,
4352                                   solver->getColSolution()) ;
4353            model->setLastHeuristic(NULL);
4354            model->incrementUsed(solver->getColSolution());
4355        }
4356    }
4357    // restore bounds
4358    {
4359        for (int j = 0; j < numberColumns; j++) {
4360            if (saveLower[j] != lower[j])
4361                solver->setColLower(j, saveLower[j]);
4362            if (saveUpper[j] != upper[j])
4363                solver->setColUpper(j, saveUpper[j]);
4364        }
4365    }
4366    // restore basis
4367    solver->setWarmStart(ws);
4368    delete ws;
4369    int anyAction;
4370    //#define CHECK_PATH
4371#ifdef CHECK_PATH
4372    extern int gotGoodNode_Z;
4373    if (gotGoodNode_Z >= 0)
4374        printf("good node %d %g\n", gotGoodNode_Z, infeasibility);
4375#endif
4376    if (infeasibility > 0.0) {
4377        if (infeasibility == COIN_DBL_MAX) {
4378            anyAction = -2; // infeasible
4379        } else {
4380            branch_ = thisOne->createCbcBranch(solver, &usefulInfo, preferredWay);
4381            // Set to firat one (and change when re-pushing)
4382            CbcGeneralBranchingObject * branch =
4383                dynamic_cast <CbcGeneralBranchingObject *> (branch_);
4384            branch->state(objectiveValue_, sumInfeasibilities_,
4385                          numberUnsatisfied_, 0);
4386            branch->setNode(this);
4387            anyAction = 0;
4388        }
4389    } else {
4390        anyAction = -1;
4391    }
4392#ifdef CHECK_PATH
4393    gotGoodNode_Z = -1;
4394#endif
4395    // Set guessed solution value
4396    guessedObjectiveValue_ = objectiveValue_ + 1.0e-5;
4397    delete [] saveLower;
4398    delete [] saveUpper;
4399
4400    // restore solution
4401    solver->setColSolution(saveSolution);
4402    delete [] saveSolution;
4403    return anyAction;
4404}
4405/* Double checks in case node can change its mind!
4406   Returns objective value
4407   Can change objective etc */
4408double
4409CbcNode::checkIsCutoff(double cutoff)
4410{
4411    branch_->checkIsCutoff(cutoff);
4412    return objectiveValue_;
4413}
4414
Note: See TracBrowser for help on using the repository browser.