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

Last change on this file since 2100 was 2100, checked in by forrest, 5 years ago

fix from stable + fix in analyze

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