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

Last change on this file since 2366 was 2366, checked in by forrest, 12 months ago

clean up when second check shows infeasible

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