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

Last change on this file since 2357 was 2357, checked in by forrest, 3 months ago

changes for SC

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 249.6 KB
Line 
1/* $Id: CbcNode.cpp 2357 2018-01-19 15:12:47Z 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                            for (int i = 0; i < numberObjects; i++) {
2386                                OsiObject * object = model->modifiableObject(i);
2387                                CbcSimpleIntegerDynamicPseudoCost * obj =
2388                                    dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2389                                assert (obj);
2390                                int iColumn = obj->columnNumber();
2391                                if (iColumn == whichAway) {
2392                                    int preferredWay = (saveSolution[iColumn] > solution[iColumn])
2393                                                       ? -1 : +1;
2394                                    usefulInfo.integerTolerance_ = 0.0;
2395                                    branch_ = obj->createCbcBranch(solver, &usefulInfo, preferredWay);
2396                                    break;
2397                                }
2398                            }
2399                            anyAction = 0;
2400                            break;
2401                        }
2402                    }
2403                    delete [] solution;
2404                }
2405            } else if (iPass == 1) {
2406                // looks like a solution - get paranoid
2407                bool roundAgain = false;
2408                // get basis
2409                CoinWarmStartBasis * ws = dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
2410                if (!ws)
2411                    break;
2412                double tolerance;
2413                solver->getDblParam(OsiPrimalTolerance, tolerance);
2414                for (i = 0; i < numberColumns; i++) {
2415                    double value = saveSolution[i];
2416                    if (value < lower[i] - tolerance) {
2417                        saveSolution[i] = lower[i];
2418                        roundAgain = true;
2419                        ws->setStructStatus(i, CoinWarmStartBasis::atLowerBound);
2420                    } else if (value > upper[i] + tolerance) {
2421                        saveSolution[i] = upper[i];
2422                        roundAgain = true;
2423                        ws->setStructStatus(i, CoinWarmStartBasis::atUpperBound);
2424                    }
2425                }
2426                if (roundAgain) {
2427                    // restore basis
2428                    solver->setWarmStart(ws);
2429                    solver->setColSolution(saveSolution);
2430                    delete ws;
2431                    bool takeHint;
2432                    OsiHintStrength strength;
2433                    solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
2434                    solver->setHintParam(OsiDoDualInResolve, false, OsiHintDo) ;
2435                    model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
2436                    double newObjValue = solver->getObjSense()*solver->getObjValue();
2437                    objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
2438                    solver->setHintParam(OsiDoDualInResolve, takeHint, strength) ;
2439                    if (!solver->isProvenOptimal()) {
2440                        // infeasible
2441                        anyAction = -2;
2442                        break;
2443                    }
2444                } else {
2445                    delete ws;
2446                    break;
2447                }
2448            }
2449        }
2450        if (anyAction == -2) {
2451            break;
2452        }
2453        // skip if solution
2454        if (!numberUnsatisfied_)
2455            break;
2456        int skipAll = (numberNotTrusted == 0 || numberToDo == 1) ? 1 : 0;
2457        bool doneHotStart = false;
2458        //DEPRECATED_STRATEGYint searchStrategy = saveSearchStrategy>=0 ? (saveSearchStrategy%10) : -1;
2459        int searchStrategy = model->searchStrategy();
2460        // But adjust depending on ratio of iterations
2461        if (searchStrategy > 0) {
2462          if (numberBeforeTrust >= /*5*/ 10 && numberBeforeTrust <= 10) {
2463                if (searchStrategy != 2) {
2464                    assert (searchStrategy == 1);
2465                    if (depth_ > 5) {
2466                        int numberIterations = model->getIterationCount();
2467                        int numberStrongIterations = model->numberStrongIterations();
2468                        if (numberStrongIterations > numberIterations + 10000) {
2469                            searchStrategy = 2;
2470                            skipAll = 1;
2471                        } else if (numberStrongIterations*4 + 1000 < numberIterations) {
2472                            searchStrategy = 3;
2473                            skipAll = 0;
2474                        }
2475                    } else {
2476                        searchStrategy = 3;
2477                        skipAll = 0;
2478                    }
2479                }
2480            }
2481        }
2482        // worth trying if too many times
2483        // Save basis
2484        CoinWarmStart * ws = NULL;
2485        // save limit
2486        int saveLimit = 0;
2487        solver->getIntParam(OsiMaxNumIterationHotStart, saveLimit);
2488        if (!numberPassesLeft)
2489            skipAll = 1;
2490        if (!skipAll) {
2491            ws = solver->getWarmStart();
2492            int limit = 100;
2493            if (!saveStateOfSearch && saveLimit < limit && saveLimit == 100)
2494                solver->setIntParam(OsiMaxNumIterationHotStart, limit);
2495        }
2496        // Say which one will be best
2497        int whichChoice = 0;
2498        int bestChoice;
2499        if (iBestGot >= 0)
2500            bestChoice = iBestGot;
2501        else
2502            bestChoice = iBestNot;
2503        assert (bestChoice >= 0);
2504        // If we have hit max time don't do strong branching
2505        bool hitMaxTime = (model->getCurrentSeconds() >
2506                            model->getDblParam(CbcModel::CbcMaximumSeconds));
2507        // also give up if we are looping round too much
2508        if (hitMaxTime || numberPassesLeft <= 0 || useShadow == 2) {
2509            int iObject = whichObject[bestChoice];
2510            OsiObject * object = model->modifiableObject(iObject);
2511            int preferredWay;
2512            object->infeasibility(&usefulInfo, preferredWay);
2513            CbcObject * obj =
2514                dynamic_cast <CbcObject *>(object) ;
2515            assert (obj);
2516            branch_ = obj->createCbcBranch(solver, &usefulInfo, preferredWay);
2517            {
2518                CbcBranchingObject * branchObj =
2519                    dynamic_cast <CbcBranchingObject *>(branch_) ;
2520                assert (branchObj);
2521                branchObj->way(preferredWay);
2522            }
2523            delete ws;
2524            ws = NULL;
2525            break;
2526        } else {
2527            // say do fast
2528            int easy = 1;
2529            if (!skipAll)
2530                solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, &easy) ;
2531            int iDo;
2532#define RESET_BOUNDS
2533#ifdef RANGING
2534            bool useRanging = model->allDynamic() && !skipAll;
2535            if (useRanging) {
2536                double currentObjective = solver->getObjValue() * solver->getObjSense();
2537                double gap = cutoff - currentObjective;
2538                // relax a bit
2539                gap *= 1.0000001;
2540                gap = CoinMax(1.0e-5, gap);
2541                // off penalties if too much
2542                double needed = neededPenalties;
2543                needed *= numberRows;
2544                if (numberNodes) {
2545                    if (needed > 1.0e6) {
2546                        neededPenalties = 0;
2547                    } else if (gap < 1.0e5) {
2548                        // maybe allow some not needed
2549                        int extra = static_cast<int> ((1.0e6 - needed) / numberRows);
2550                        int nStored = numberObjects - optionalPenalties;
2551                        extra = CoinMin(extra, nStored);
2552                        for (int i = 0; i < extra; i++) {
2553                            objectMark[neededPenalties] = objectMark[optionalPenalties+i];
2554                            which[neededPenalties++] = which[optionalPenalties+i];;
2555                        }
2556                    }
2557                }
2558                if (osiclp && neededPenalties) {
2559                    assert (!doneHotStart);
2560                    xPen += neededPenalties;
2561                    which--;
2562                    which[0] = neededPenalties;
2563                    osiclp->passInRanges(which);
2564                    // Mark hot start and get ranges
2565                    if (kPass) {
2566                        // until can work out why solution can go funny
2567                        int save = osiclp->specialOptions();
2568                        osiclp->setSpecialOptions(save | 256);
2569                        solver->markHotStart();
2570#ifdef RESET_BOUNDS
2571                        memcpy(saveLower,solver->getColLower(),solver->getNumCols()*sizeof(double));
2572                        memcpy(saveUpper,solver->getColUpper(),solver->getNumCols()*sizeof(double));
2573#endif
2574                        osiclp->setSpecialOptions(save);
2575                    } else {
2576                        solver->markHotStart();
2577#ifdef RESET_BOUNDS
2578                        memcpy(saveLower,solver->getColLower(),solver->getNumCols()*sizeof(double));
2579                        memcpy(saveUpper,solver->getColUpper(),solver->getNumCols()*sizeof(double));
2580#endif
2581                    }
2582                    doneHotStart = true;
2583                    xMark++;
2584                    kPass++;
2585                    osiclp->passInRanges(NULL);
2586                    const double * downCost = osiclp->upRange();
2587                    const double * upCost = osiclp->downRange();
2588                    bool problemFeasible = true;
2589                    int numberFixed = 0;
2590                    for (int i = 0; i < neededPenalties; i++) {
2591                        int j = objectMark[i];
2592                        int iObject = whichObject[j];
2593                        OsiObject * object = model->modifiableObject(iObject);
2594                        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2595                            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2596                        // Use this object's numberBeforeTrust
2597                        int numberBeforeTrustThis = dynamicObject->numberBeforeTrust();
2598                        int iSequence = dynamicObject->columnNumber();
2599                        double value = saveSolution[iSequence];
2600                        value -= floor(value);
2601                        double upPenalty = CoinMin(upCost[i], 1.0e110) * (1.0 - value);
2602                        double downPenalty = CoinMin(downCost[i], 1.0e110) * value;
2603                        int numberThisDown = dynamicObject->numberTimesDown();
2604                        int numberThisUp = dynamicObject->numberTimesUp();
2605                        if (!numberBeforeTrustThis) {
2606                            // override
2607                            downEstimate[iObject] = downPenalty;
2608                            upEstimate[iObject] = upPenalty;
2609                            double min1 = CoinMin(downEstimate[iObject],
2610                                                  upEstimate[iObject]);
2611                            double max1 = CoinMax(downEstimate[iObject],
2612                                                  upEstimate[iObject]);
2613                            min1 = 0.8 * min1 + 0.2 * max1;
2614                            sort[j] = - min1;
2615                        } else if (numberThisDown < numberBeforeTrustThis ||
2616                                   numberThisUp < numberBeforeTrustThis) {
2617                            double invTrust = 1.0 / static_cast<double> (numberBeforeTrustThis);
2618                            if (numberThisDown < numberBeforeTrustThis) {
2619                                double fraction = numberThisDown * invTrust;
2620                                downEstimate[iObject] = fraction * downEstimate[iObject] + (1.0 - fraction) * downPenalty;
2621                            }
2622                            if (numberThisUp < numberBeforeTrustThis) {
2623                                double fraction = numberThisUp * invTrust;
2624                                upEstimate[iObject] = fraction * upEstimate[iObject] + (1.0 - fraction) * upPenalty;
2625                            }
2626                            double min1 = CoinMin(downEstimate[iObject],
2627                                                  upEstimate[iObject]);
2628                            double max1 = CoinMax(downEstimate[iObject],
2629                                                  upEstimate[iObject]);
2630                            min1 = 0.8 * min1 + 0.2 * max1;
2631                            min1 *= 10.0;
2632                            if (!(numberThisDown + numberThisUp))
2633                                min1 *= 100.0;
2634                            sort[j] = - min1;
2635                        }
2636                        // seems unreliable
2637                        if (false&&CoinMax(downPenalty, upPenalty) > gap) {
2638                            COIN_DETAIL_PRINT(printf("gap %g object %d has down range %g, up %g\n",
2639                                                     gap, i, downPenalty, upPenalty));
2640                            //sort[j] -= 1.0e50; // make more likely to be chosen
2641                            int number;
2642                            if (downPenalty > gap) {
2643                                number = dynamicObject->numberTimesDown();
2644                                if (upPenalty > gap)
2645                                    problemFeasible = false;
2646                                CbcBranchingObject * branch = dynamicObject->createCbcBranch(solver, &usefulInfo, 1);
2647                                //branch->fix(solver,saveLower,saveUpper,1);
2648                                delete branch;
2649                            } else {
2650                                number = dynamicObject->numberTimesUp();
2651                                CbcBranchingObject * branch = dynamicObject->createCbcBranch(solver, &usefulInfo, 1);
2652                                //branch->fix(solver,saveLower,saveUpper,-1);
2653                                delete branch;
2654                            }
2655                            if (number >= numberBeforeTrustThis)
2656                              dynamicObject->setNumberBeforeTrust(CoinMin(number + 1,5*numberBeforeTrust));
2657                            numberFixed++;
2658                        }
2659                        if (!numberNodes)
2660                            COIN_DETAIL_PRINT(printf("%d pen down ps %g -> %g up ps %g -> %g\n",
2661                                                     iObject, downPenalty, downPenalty, upPenalty, upPenalty));
2662                    }
2663                    if (numberFixed && problemFeasible) {
2664                        assert(doneHotStart);
2665                        solver->unmarkHotStart();
2666                        model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
2667#ifdef CHECK_DEBUGGER_PATH
2668                        if ((model->specialOptions()&1) != 0 && onOptimalPath) {
2669                          const OsiRowCutDebugger *debugger = solver->getRowCutDebugger() ;
2670                          if (!debugger) {
2671                            printf("Strong branching down on %d went off optimal path\n",iObject);
2672                            abort();
2673                          }
2674                        }
2675#endif
2676                        double newObjValue = solver->getObjSense()*solver->getObjValue();
2677                        objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
2678                        solver->markHotStart();
2679#ifdef RESET_BOUNDS
2680                        memcpy(saveLower,solver->getColLower(),solver->getNumCols()*sizeof(double));
2681                        memcpy(saveUpper,solver->getColUpper(),solver->getNumCols()*sizeof(double));
2682#endif
2683                        problemFeasible = solver->isProvenOptimal();
2684                    }
2685                    if (!problemFeasible) {
2686                      COIN_DETAIL_PRINT(fprintf(stdout, "both ways infeas on ranging - code needed\n"));
2687                        anyAction = -2;
2688                        if (!choiceObject) {
2689                            delete choice.possibleBranch;
2690                            choice.possibleBranch = NULL;
2691                        }
2692                        //printf("Both infeasible for choice %d sequence %d\n",i,
2693                        // model->object(choice.objectNumber)->columnNumber());
2694                        // Delete the snapshot
2695                        solver->unmarkHotStart();
2696                        // back to normal
2697                        solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
2698                        // restore basis
2699                        solver->setWarmStart(ws);
2700                        doneHotStart = false;
2701                        delete ws;
2702                        ws = NULL;
2703                        break;
2704                    }
2705                }
2706            }
2707#endif          /* RANGING */
2708            {
2709                int numberIterations = model->getIterationCount();
2710                //numberIterations += (model->numberExtraIterations()>>2);
2711                const int * strongInfo = model->strongInfo();
2712                //int numberDone = strongInfo[0]-strongInfo[3];
2713                int numberFixed = strongInfo[1] - strongInfo[4];
2714                int numberInfeasible = strongInfo[2] - strongInfo[5];
2715                assert (!strongInfo[3]);
2716                assert (!strongInfo[4]);
2717                assert (!strongInfo[5]);
2718                int numberStrongIterations = model->numberStrongIterations();
2719                int numberRows = solver->getNumRows();
2720                if (numberStrongIterations > numberIterations + CoinMin(100, 10*numberRows) && depth_ >= 4 && numberNodes > 100) {
2721                    if (20*numberInfeasible + 4*numberFixed < numberNodes) {
2722                        // Say never do
2723                        if (numberBeforeTrust == 10)
2724                          skipAll = -1;
2725                    }
2726                }
2727            }
2728            // make sure best will be first
2729            if (iBestGot >= 0)
2730                sort[iBestGot] = -COIN_DBL_MAX;
2731            // Actions 0 - exit for repeat, 1 resolve and try old choice,2 exit for continue
2732            if (anyAction)
2733                numberToDo = 0; // skip as we will be trying again
2734            // Sort
2735            CoinSort_2(sort, sort + numberToDo, whichObject);
2736            // Change in objective opposite infeasible
2737            double worstFeasible = 0.0;
2738            // Just first if strong off
2739            if (!numberStrong)
2740                numberToDo = CoinMin(numberToDo, 1);
2741            if (searchStrategy == 2)
2742                numberToDo = CoinMin(numberToDo, 20);
2743            iDo = 0;
2744            int saveLimit2;
2745            solver->getIntParam(OsiMaxNumIterationHotStart, saveLimit2);
2746            int numberTest = numberNotTrusted > 0 ? numberStrong : (numberStrong + 1) / 2;
2747            if (searchStrategy == 3) {
2748                // Previously decided we need strong
2749                numberTest = numberStrong;
2750            }
2751            // Try nearly always off
2752            if (skipAll >= 0) {
2753                if (searchStrategy < 2) {
2754                    //if ((numberNodes%20)!=0) {
2755                    if ((model->specialOptions()&8) == 0) {
2756                        numberTest = 0;
2757                    }
2758                    //} else {
2759                    //numberTest=2*numberStrong;
2760                    //skipAll=0;
2761                    //}
2762                }
2763            } else {
2764                // Just take first
2765                numberTest = 1;
2766            }
2767            int testDepth = (skipAll >= 0) ? 8 : 4;
2768            if (depth_ < testDepth && numberStrong) {
2769                if (searchStrategy != 2) {
2770                    int numberRows = solver->getNumRows();
2771                    // whether to do this or not is important - think
2772                    if (numberRows < 300 || numberRows + numberColumns < 2500) {
2773                        if (depth_ < 7)
2774                            numberStrong = CoinMin(3 * numberStrong, numberToDo);
2775                        if (!depth_)
2776                            numberStrong = CoinMin(6 * numberStrong, numberToDo);
2777                    }
2778                    numberTest = numberStrong;
2779                    skipAll = 0;
2780                }
2781            }
2782            // Do at least 5 strong
2783            if (numberColumns < 1000 && (depth_ < 15 || numberNodes < 1000000))
2784                numberTest = CoinMax(numberTest, 5);
2785            if ((model->specialOptions()&8) == 0) {
2786                if (skipAll) {
2787                    numberTest = 0;
2788                }
2789            } else {
2790                // do 5 as strong is fixing
2791                numberTest = CoinMax(numberTest, 5);
2792            }
2793            // see if switched off
2794            if (skipAll < 0) {
2795                numberTest = 0;
2796            }
2797            int realMaxHotIterations = 999999;
2798            if (skipAll < 0)
2799                numberToDo = 1;
2800            strongType=0;
2801#ifdef DO_ALL_AT_ROOT
2802            if (model->strongStrategy()) {
2803              int iStrategy=model->strongStrategy();
2804              int kDepth = iStrategy/100;
2805              if (kDepth)
2806                iStrategy -= 100*kDepth;
2807              else
2808                kDepth=5;
2809              double objValue = solver->getObjSense()*solver->getObjValue();
2810              double bestPossible = model->getBestPossibleObjValue();
2811              bestPossible += 1.0e-7*(1.0+fabs(bestPossible));
2812              int jStrategy = iStrategy/10;
2813              if (jStrategy) {
2814                if ((jStrategy&1)!=0&&!depth_)
2815                  strongType=2;
2816                else if ((jStrategy&2)!=0&&depth_<=kDepth)
2817                  strongType=2;
2818                else if ((jStrategy&4)!=0&&objValue<bestPossible)
2819                  strongType=2;
2820                iStrategy-=10*jStrategy;
2821              }
2822              if (!strongType) {
2823                if ((iStrategy&1)!=0&&!depth_)
2824                  strongType=1;
2825                else if ((iStrategy&2)!=0&&depth_<=kDepth)
2826                  strongType=1;
2827                else if ((iStrategy&4)!=0&&objValue<bestPossible)
2828                  strongType=1;
2829              }
2830              saveNumberToDo=numberToDo;
2831              if (strongType==2) {
2832                // add in satisfied
2833                const int * integerVariable = model->integerVariable();
2834                int numberIntegers = model->numberIntegers();
2835                if (numberIntegers==numberObjects) {
2836                  numberToDo=0;
2837                  for (int i=0;i<numberIntegers;i++) {
2838                    int iColumn=integerVariable[i];
2839                    if (saveUpper[iColumn]>saveLower[iColumn]) {
2840                      whichObject [numberToDo++]=i;
2841                    }
2842                  }
2843                  saveSatisfiedVariables=numberToDo-saveNumberToDo;
2844                } else {
2845                  strongType=1;
2846                }
2847              }
2848              if (strongType) {
2849                numberTest = numberToDo;
2850                numberStrong=numberToDo;
2851                skipAll=0;
2852                searchStrategy=0;
2853                solver->setIntParam(OsiMaxNumIterationHotStart, 100000);
2854                //printf("Strong branching type %d\n",strongType);
2855              }
2856            }
2857#endif
2858#ifdef COIN_HAS_NTY
2859            const int * orbits = NULL;
2860#endif
2861#ifdef COIN_HAS_NTY
2862            if (orbitOption==2 /* was >1*/) {
2863              CbcSymmetry * symmetryInfo = model->symmetryInfo();
2864              CbcNodeInfo * infoX = lastNode ? lastNode->nodeInfo() : NULL;
2865              bool worthTrying = false;
2866              if (infoX) {
2867                CbcNodeInfo * info = infoX;
2868                for (int i=0;i<NTY_BAD_DEPTH;i++) {
2869                  if (!info->parent()) {
2870                    worthTrying = true;
2871                    break;
2872                  }
2873                  info = info->parent();
2874                  if (info->symmetryWorked()) {
2875                    worthTrying = true;
2876                    break;
2877                  }
2878                }
2879              } else {
2880                worthTrying=true;
2881              }
2882              if (symmetryInfo && worthTrying) {
2883                symmetryInfo->ChangeBounds(solver->getColLower(),
2884                                            solver->getColUpper(),
2885                                            solver->getNumCols(),false);
2886                symmetryInfo->Compute_Symmetry();
2887                symmetryInfo->fillOrbits();
2888                orbits = symmetryInfo->whichOrbit();
2889                int iColumn=-1;
2890                if (orbits && symmetryInfo->numberUsefulOrbits()) {
2891                  bool doBranch=true;
2892                  int numberUsefulOrbits = symmetryInfo->numberUsefulOrbits();
2893                  if (numberUsefulOrbits<2) {
2894                    assert (numberUsefulOrbits);
2895                    double largest=-1.0;
2896                    for (int i=0;i<numberColumns;i++) {
2897                      if (orbits[i]>=0) {
2898                        if (saveSolution[i]>largest) {
2899                          largest=saveSolution[i];
2900                          iColumn=i;
2901                        }
2902                      }
2903                    }
2904                  } else {
2905#if COIN_HAS_NTY2 == 1
2906                    // take largest
2907                    int iOrbit=symmetryInfo->largestOrbit(solver->getColLower(),
2908                                                          solver->getColUpper());
2909                    double largest=-1.0;
2910                    for (int i=0;i<numberColumns;i++) {
2911                      if (orbits[i]==iOrbit) {
2912                        if (saveSolution[i]>largest) {
2913                          largest=saveSolution[i];
2914                          iColumn=i;
2915                        }
2916                      }
2917                    }
2918#endif
2919                    if (orbitOption==2) {
2920                      // strong
2921                      int nDo=0;
2922                      const double * lower = solver->getColLower();
2923                      const double * upper = solver->getColUpper();
2924                      const int * integerVariable = model->integerVariable();
2925                      for (int iOrbit = 0; iOrbit < numberUsefulOrbits; iOrbit++) {
2926                        double distance=1.0;
2927                        int iColumn = -1;
2928                        int numberIntegers = model->numberIntegers();
2929                        for (int j=0;j<numberIntegers;j++) {
2930                          int i=integerVariable[j];
2931                          if (orbits[i]==iOrbit &&lower[i]==0.0&&upper[i]==1.0) {
2932                            double away = fabs(saveSolution[i]-0.5);
2933                            if (away<distance&&away<0.4999) {
2934                              distance=away;
2935                              iColumn=j;
2936                            }
2937                          }
2938                        }
2939                        if (iColumn>=0) 
2940                          whichObject[nDo++]=iColumn;
2941                      }
2942                      if (nDo)
2943                        numberToDo=nDo;
2944                      doBranch=false;
2945                    } else if (orbitOption==3) {
2946                      // subset
2947                      int nDo=0;
2948                      for (int iDo = 0; iDo < numberToDo; iDo++) {
2949                        int iObject = whichObject[iDo];
2950                        OsiObject * object = model->modifiableObject(iObject);
2951                        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2952                          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2953                        int iColumn = dynamicObject ? dynamicObject->columnNumber() : -1;
2954                        if (iColumn<0||orbits[iColumn]>=0)
2955                          whichObject[nDo++]=whichObject[iDo];
2956                      }
2957                      assert(nDo);
2958                      //printf("nDo %d\n",nDo);
2959                      numberToDo=nDo;
2960                      doBranch=false;
2961                      /* need NULL as if two in same orbit and strong branching fixes
2962                         then we may be in trouble.
2963                         Strong option should be OK as only one in set done.
2964                       */ 
2965                      orbits=NULL;
2966                    }
2967                  }
2968                  if(doBranch) {
2969                    orbitOption=0;
2970                    branch_ = new CbcOrbitalBranchingObject(model,iColumn,1,0,NULL);
2971                    if (infoX)
2972                      infoX->setSymmetryWorked();
2973                    numberToDo=0;
2974                  }
2975                }
2976              }
2977            }
2978#endif
2979            for ( iDo = 0; iDo < numberToDo; iDo++) {
2980                int iObject = whichObject[iDo];
2981                OsiObject * object = model->modifiableObject(iObject);
2982                CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
2983                    dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
2984                int iColumn = dynamicObject ? dynamicObject->columnNumber() : numberColumns + iObject;
2985                int preferredWay;
2986                double infeasibility = object->infeasibility(&usefulInfo, preferredWay);
2987                bool feasibleSolution=false;
2988                double predictedChange=0.0;
2989                // may have become feasible
2990                if (!infeasibility) {
2991                  if(strongType!=2||solver->getColLower()[iColumn]==solver->getColUpper()[iColumn])
2992                    continue;
2993                }
2994#ifndef NDEBUG
2995                if (iColumn < numberColumns) {
2996                    const double * solution = model->testSolution();
2997                    assert (saveSolution[iColumn] == solution[iColumn]);
2998                }
2999#endif
3000                CbcSimpleInteger * obj =
3001                    dynamic_cast <CbcSimpleInteger *>(object) ;
3002                if (obj) {
3003                    if (choiceObject) {
3004                        obj->fillCreateBranch(choiceObject, &usefulInfo, preferredWay);
3005                        choiceObject->setObject(dynamicObject);
3006                    } else {
3007                        choice.possibleBranch = obj->createCbcBranch(solver, &usefulInfo, preferredWay);
3008                    }
3009                } else {
3010                    CbcObject * obj =
3011                        dynamic_cast <CbcObject *>(object) ;
3012                    assert (obj);
3013                    choice.possibleBranch = obj->createCbcBranch(solver, &usefulInfo, preferredWay);
3014                }
3015                // Save which object it was
3016                choice.objectNumber = iObject;
3017                choice.numIntInfeasUp = numberUnsatisfied_;
3018                choice.numIntInfeasDown = numberUnsatisfied_;
3019                if (strongType!=2) {
3020                  choice.upMovement = upEstimate[iObject];
3021                  choice.downMovement = downEstimate[iObject];
3022                } else {
3023                  choice.upMovement = 0.1;
3024                  choice.downMovement = 0.1;
3025                }
3026                assert (choice.upMovement >= 0.0);
3027                assert (choice.downMovement >= 0.0);
3028                choice.fix = 0; // say not fixed
3029                // see if can skip strong branching
3030                int canSkip = choice.possibleBranch->fillStrongInfo(choice);
3031                if ((numberTest <= 0 || skipAll)) {
3032                    if (iDo > 20) {
3033                        if (!choiceObject) {
3034                            delete choice.possibleBranch;
3035                            choice.possibleBranch = NULL;
3036                        }
3037                        break; // give up anyway
3038                    }
3039                }
3040                if (model->messageHandler()->logLevel() > 3 && numberBeforeTrust && dynamicObject)
3041                    dynamicObject->print(1, choice.possibleBranch->value());
3042                if (strongType)
3043                  canSkip=0;
3044                if (skipAll < 0)
3045                    canSkip = 1;
3046                if (!canSkip) {
3047                    if (!doneHotStart) {
3048                        // Mark hot start
3049                        doneHotStart = true;
3050                        solver->markHotStart();
3051#ifdef RESET_BOUNDS
3052                        memcpy(saveLower,solver->getColLower(),solver->getNumCols()*sizeof(double));
3053                        memcpy(saveUpper,solver->getColUpper(),solver->getNumCols()*sizeof(double));
3054#endif
3055                        if (!solver->isProvenOptimal()) {
3056                          skipAll=-2;
3057                          canSkip = 1;
3058                        }
3059                        xMark++;
3060                    }
3061                }
3062                if (!canSkip) {
3063                    numberTest--;
3064                    // just do a few
3065                    if (searchStrategy == 2)
3066                        solver->setIntParam(OsiMaxNumIterationHotStart, 10);
3067                    double objectiveChange ;
3068                    double newObjectiveValue = 1.0e100;
3069                    int j;
3070#ifdef COIN_HAS_CLP
3071                    int saveMaxHotIts=0;
3072                    int saveOsiClpOptions=0;
3073                    if (osiclp && goToEndInStrongBranching) {
3074                      /* How dogged to be in strong branching
3075                         0 - default
3076                         1 - go to end on first time
3077                         2 - always go to end
3078                      */
3079                      osiclp->getIntParam(OsiMaxNumIterationHotStart, saveMaxHotIts);
3080                      saveOsiClpOptions=osiclp->specialOptions();
3081                      if (goToEndInStrongBranching==2 ||
3082                          (dynamicObject &&
3083                           dynamicObject->numberTimesBranched()==0)) {
3084                        if (osiclp->getNumRows()<200||goToEndInStrongBranching==2) {
3085                        osiclp->setIntParam(OsiMaxNumIterationHotStart, 
3086                                            10*(osiclp->getNumRows()+numberColumns));
3087                        osiclp->setSpecialOptions(saveOsiClpOptions & (~32));
3088                        }
3089                      }
3090                    }
3091#endif
3092                    // status is 0 finished, 1 infeasible and other
3093                    int iStatus;
3094                    /*
3095                      Try the down direction first. (Specify the initial branching alternative as
3096                      down with a call to way(-1). Each subsequent call to branch() performs the
3097                      specified branch and advances the branch object state to the next branch
3098                      alternative.)
3099                    */
3100                    choice.possibleBranch->way(-1) ;
3101                    predictedChange = choice.possibleBranch->branch() ;
3102#ifdef COIN_HAS_NTY
3103                    if (orbits) {
3104                      // can fix all in orbit
3105                      int fixOrbit = orbits[iObject];
3106                      if (fixOrbit>=0) {
3107                        //printf("fixing all in orbit %d for column %d\n",fixOrbit,iObject);
3108                        for (int i=0;i<numberColumns;i++) {
3109                          if (orbits[i]==fixOrbit)
3110                            solver->setColUpper(i,0.0);
3111                        }
3112                      }
3113                    }
3114#endif
3115                    solver->solveFromHotStart() ;
3116                    bool needHotStartUpdate = false;
3117                    numberStrongDone++;
3118                    numberStrongIterations += solver->getIterationCount();
3119                    /*
3120                      We now have an estimate of objective degradation that we can use for strong
3121                      branching. If we're over the cutoff, the variable is monotone up.
3122                      If we actually made it to optimality, check for a solution, and if we have
3123                      a good one, call setBestSolution to process it. Note that this may reduce the
3124                      cutoff, so we check again to see if we can declare this variable monotone.
3125                    */
3126                    if (solver->isProvenOptimal())
3127                        iStatus = 0; // optimal
3128                    else if (solver->isIterationLimitReached() 
3129                             && !solver->isDualObjectiveLimitReached()) {
3130                        iStatus = 2; // unknown
3131                    } else {
3132                        iStatus = 1; // infeasible
3133#ifdef CONFLICT_CUTS
3134#undef CONFLICT_CUTS
3135                        //#define CONFLICT_CUTS 2
3136#endif
3137#ifdef CONFLICT_CUTS
3138# ifdef COIN_HAS_CLP
3139                        if (osiclp&&(model->moreSpecialOptions()&4194304)!=0) {
3140                          const CbcFullNodeInfo * topOfTree =
3141                            model->topOfTree();
3142                          if (topOfTree) {
3143#if CONFLICT_CUTS==2
3144                            OsiRowCut * cut = osiclp->smallModelCut(topOfTree->lower(),
3145                                                                    topOfTree->upper(),
3146                                                                    model->numberRowsAtContinuous(),
3147                                                                    model->whichGenerator());
3148#else
3149                            OsiRowCut * cut = osiclp->modelCut(topOfTree->lower(),
3150                                                               topOfTree->upper(),
3151                                                               model->numberRowsAtContinuous(),
3152                                                               model->whichGenerator(),0);
3153#endif
3154                            if (cut) {
3155                              if (model->messageHandler()->logLevel() > 1)
3156                                printf("Conflict cut found in strong branching (%d elements)\n",
3157                                       cut->row().getNumElements());
3158                              //cut->print();
3159                              if ((model->specialOptions()&1) != 0) {
3160                                const OsiRowCutDebugger *debugger = model->continuousSolver()->getRowCutDebugger() ;
3161                                if (debugger) {
3162                                  if (debugger->invalidCut(*cut)) {
3163                                    model->continuousSolver()->applyRowCuts(1,cut);
3164                                    model->continuousSolver()->writeMps("bad");
3165                                  }
3166                                  CoinAssert (!debugger->invalidCut(*cut));
3167                                }
3168                              }
3169                              model->makeGlobalCut(cut) ;
3170                            }
3171                          }
3172                        }
3173#endif
3174#endif
3175                    }
3176                    // say infeasible if branch says so
3177                    if (predictedChange==COIN_DBL_MAX)
3178                      iStatus=1;
3179                    if (iStatus != 2 && solver->getIterationCount() >
3180                            realMaxHotIterations)
3181                        numberUnfinished++;
3182                    newObjectiveValue = solver->getObjSense() * solver->getObjValue();
3183                    choice.numItersDown = solver->getIterationCount();
3184                    objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_, 0.0);
3185                    // Update branching information if wanted
3186                    CbcBranchingObject * cbcobj = dynamic_cast<CbcBranchingObject *> (choice.possibleBranch);
3187                    if (cbcobj) {
3188                        CbcObject * object = cbcobj->object();
3189                        assert (object) ;
3190                        CbcObjectUpdateData update = object->createUpdateInformation(solver, this, cbcobj);
3191                        update.objectNumber_ = choice.objectNumber;
3192                        model->addUpdateInformation(update);
3193                    } else {
3194                        decision->updateInformation( solver, this);
3195                    }
3196                    if (!iStatus) {
3197                        choice.finishedDown = true ;
3198                        if (newObjectiveValue >= cutoff) {
3199                            objectiveChange = 1.0e100; // say infeasible
3200                            numberStrongInfeasible++;
3201                        } else {
3202#define CBCNODE_TIGHTEN_BOUNDS
3203#ifdef CBCNODE_TIGHTEN_BOUNDS
3204                            // Can we tighten bounds?
3205                            if (iColumn < numberColumns && cutoff < 1.0e20
3206                                    && objectiveChange > 1.0e-5) {
3207                                double value = saveSolution[iColumn];
3208                                double down = value - floor(value-integerTolerance);
3209                                double changePer = objectiveChange / (down + 1.0e-7);
3210                                double distance = (cutoff - objectiveValue_) / changePer;
3211                                distance += 1.0e-3;
3212                                if (distance < 5.0) {
3213                                    double newLower = ceil(value - distance);
3214                                    if (newLower > saveLower[iColumn]) {
3215                                        //printf("Could increase lower bound on %d from %g to %g\n",
3216                                        //   iColumn,saveLower[iColumn],newLower);
3217                                        saveLower[iColumn] = newLower;
3218                                        solver->setColLower(iColumn, newLower);
3219                                    }
3220                                }
3221                            }
3222#endif
3223                            // See if integer solution
3224                            feasibleSolution = 
3225                              model->feasibleSolution(choice.numIntInfeasDown,
3226                                                      choice.numObjInfeasDown);
3227                            if (feasibleSolution
3228                                    && model->problemFeasibility()->feasible(model, -1) >= 0) {
3229                                if (auxiliaryInfo->solutionAddsCuts()) {
3230                                    needHotStartUpdate = true;
3231                                    solver->unmarkHotStart();
3232                                }
3233                                model->setLogLevel(saveLogLevel);
3234                                model->setBestSolution(CBC_STRONGSOL,
3235                                                       newObjectiveValue,
3236                                                       solver->getColSolution()) ;
3237                                if (needHotStartUpdate) {
3238                                    model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3239                                    newObjectiveValue = solver->getObjSense() * solver->getObjValue();
3240                                    objectiveValue_ = CoinMax(objectiveValue_,newObjectiveValue);
3241                                    objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_, 0.0);
3242                                    model->feasibleSolution(choice.numIntInfeasDown,
3243                                                            choice.numObjInfeasDown);
3244                                }
3245                                model->setLastHeuristic(NULL);
3246                                model->incrementUsed(solver->getColSolution());
3247                                cutoff = model->getCutoff();
3248                                if (newObjectiveValue >= cutoff) {      //  *new* cutoff
3249                                    objectiveChange = 1.0e100 ;
3250                                    numberStrongInfeasible++;
3251                                }
3252                            }
3253                        }
3254                    } else if (iStatus == 1) {
3255                        choice.finishedDown = true ;
3256                        objectiveChange = COIN_DBL_MAX ;
3257                        numberStrongInfeasible++;
3258                    } else {
3259                        // Can't say much as we did not finish
3260                        choice.finishedDown = false ;
3261                        numberUnfinished++;
3262                    }
3263                    choice.downMovement = objectiveChange ;
3264
3265                    // restore bounds
3266                    for ( j = 0; j < numberColumns; j++) {
3267                        if (saveLower[j] != lower[j])
3268                            solver->setColLower(j, saveLower[j]);
3269                        if (saveUpper[j] != upper[j])
3270                            solver->setColUpper(j, saveUpper[j]);
3271                    }
3272                    if (needHotStartUpdate) {
3273                        needHotStartUpdate = false;
3274                        model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3275#ifdef CHECK_DEBUGGER_PATH
3276                        if ((model->specialOptions()&1) != 0 && onOptimalPath) {
3277                          const OsiRowCutDebugger *debugger = solver->getRowCutDebugger() ;
3278                          if (!debugger) {
3279                            printf("Strong branching down on %d went off optimal path\n",iObject);
3280                            model->solver()->writeMps("query");
3281                            abort();
3282                          }
3283                        }
3284#endif
3285                        double newObjValue = solver->getObjSense()*solver->getObjValue();
3286                        objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3287                        //we may again have an integer feasible solution
3288                        int numberIntegerInfeasibilities;
3289                        int numberObjectInfeasibilities;
3290                        if (model->feasibleSolution(
3291                                    numberIntegerInfeasibilities,
3292                                    numberObjectInfeasibilities)) {
3293#ifdef BONMIN
3294                            //In this case node has become integer feasible, let us exit the loop
3295                            std::cout << "Node has become integer feasible" << std::endl;
3296                            numberUnsatisfied_ = 0;
3297                            break;
3298#endif
3299                            double objValue = solver->getObjValue();
3300                            model->setLogLevel(saveLogLevel);
3301                            model->setBestSolution(CBC_STRONGSOL,
3302                                                   objValue,
3303                                                   solver->getColSolution()) ;
3304                            model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3305                            double newObjValue = solver->getObjSense()*solver->getObjValue();
3306                            objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3307                            cutoff = model->getCutoff();
3308                        }
3309                        solver->markHotStart();
3310#ifdef RESET_BOUNDS
3311                        memcpy(saveLower,solver->getColLower(),solver->getNumCols()*sizeof(double));
3312                        memcpy(saveUpper,solver->getColUpper(),solver->getNumCols()*sizeof(double));
3313#endif
3314                        if (!solver->isProvenOptimal()) {
3315                          skipAll=-2;
3316                          canSkip = 1;
3317                        }
3318                        xMark++;
3319                    }
3320#if 0 //def DO_ALL_AT_ROOT
3321                    if (strongType)
3322                        printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
3323                               choice.objectNumber, iStatus, newObjectiveValue, choice.numItersDown,
3324                               choice.downMovement, choice.finishedDown, choice.numIntInfeasDown,
3325                               choice.numObjInfeasDown);
3326#endif
3327
3328                    // repeat the whole exercise, forcing the variable up
3329                    predictedChange=choice.possibleBranch->branch();
3330                    solver->solveFromHotStart() ;
3331#ifdef COIN_HAS_CLP
3332                    if (osiclp && goToEndInStrongBranching) {
3333                      osiclp->setIntParam(OsiMaxNumIterationHotStart, saveMaxHotIts);
3334                      osiclp->setSpecialOptions(saveOsiClpOptions);
3335                    }
3336#endif
3337                    numberStrongDone++;
3338                    numberStrongIterations += solver->getIterationCount();
3339                    /*
3340                      We now have an estimate of objective degradation that we can use for strong
3341                      branching. If we're over the cutoff, the variable is monotone up.
3342                      If we actually made it to optimality, check for a solution, and if we have
3343                      a good one, call setBestSolution to process it. Note that this may reduce the
3344                      cutoff, so we check again to see if we can declare this variable monotone.
3345                    */
3346                    if (solver->isProvenOptimal())
3347                        iStatus = 0; // optimal
3348                    else if (solver->isIterationLimitReached()
3349                             && !solver->isDualObjectiveLimitReached()) {
3350                        iStatus = 2; // unknown
3351                    } else {
3352                        iStatus = 1; // infeasible
3353#ifdef CONFLICT_CUTS
3354# ifdef COIN_HAS_CLP
3355                        if (osiclp&&(model->moreSpecialOptions()&4194304)!=0) {
3356                          const CbcFullNodeInfo * topOfTree =
3357                            model->topOfTree();
3358                          if (topOfTree) {
3359#if CONFLICT_CUTS==2
3360                            OsiRowCut * cut = osiclp->smallModelCut(topOfTree->lower(),
3361                                                                    topOfTree->upper(),
3362                                                                    model->numberRowsAtContinuous(),
3363                                                                    model->whichGenerator());
3364#else
3365                            OsiRowCut * cut = osiclp->modelCut(topOfTree->lower(),
3366                                                               topOfTree->upper(),
3367                                                               model->numberRowsAtContinuous(),
3368                                                               model->whichGenerator(),0);
3369#endif
3370                            if (cut) {
3371                              //printf("XXXXXX found conflict cut in strong branching\n");
3372                              //cut->print();
3373                              if ((model->specialOptions()&1) != 0) {
3374                                const OsiRowCutDebugger *debugger = model->continuousSolver()->getRowCutDebugger() ;
3375                                if (debugger) {
3376                                  if (debugger->invalidCut(*cut)) {
3377                                    model->continuousSolver()->applyRowCuts(1,cut);
3378                                    model->continuousSolver()->writeMps("bad");
3379                                  }
3380                                  CoinAssert (!debugger->invalidCut(*cut));
3381                                }
3382                              }
3383                              model->makeGlobalCut(cut) ;
3384                            }
3385                          }
3386                        }
3387#endif
3388#endif
3389                    }
3390                    // say infeasible if branch says so
3391                    if (predictedChange==COIN_DBL_MAX)
3392                      iStatus=1;
3393                    if (iStatus != 2 && solver->getIterationCount() >
3394                            realMaxHotIterations)
3395                        numberUnfinished++;
3396                    newObjectiveValue = solver->getObjSense() * solver->getObjValue();
3397                    choice.numItersUp = solver->getIterationCount();
3398                    objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_, 0.0);
3399                    // Update branching information if wanted
3400                    cbcobj = dynamic_cast<CbcBranchingObject *> (choice.possibleBranch);
3401                    if (cbcobj) {
3402                        CbcObject * object = cbcobj->object();
3403                        assert (object) ;
3404                        CbcObjectUpdateData update = object->createUpdateInformation(solver, this, cbcobj);
3405                        update.objectNumber_ = choice.objectNumber;
3406                        model->addUpdateInformation(update);
3407                    } else {
3408                        decision->updateInformation( solver, this);
3409                    }
3410                    if (!iStatus) {
3411                        choice.finishedUp = true ;
3412                        if (newObjectiveValue >= cutoff) {
3413                            objectiveChange = 1.0e100; // say infeasible
3414                            numberStrongInfeasible++;
3415                        } else {
3416#ifdef CBCNODE_TIGHTEN_BOUNDS
3417                            // Can we tighten bounds?
3418                            if (iColumn < numberColumns && cutoff < 1.0e20
3419                                    && objectiveChange > 1.0e-5) {
3420                                double value = saveSolution[iColumn];
3421                                double up = ceil(value+integerTolerance) - value;
3422                                double changePer = objectiveChange / (up + 1.0e-7);
3423                                double distance = (cutoff - objectiveValue_) / changePer;
3424                                distance += 1.0e-3;
3425                                if (distance < 5.0) {
3426                                    double newUpper = floor(value + distance);
3427                                    if (newUpper < saveUpper[iColumn]) {
3428                                        //printf("Could decrease upper bound on %d from %g to %g\n",
3429                                        //   iColumn,saveUpper[iColumn],newUpper);
3430                                        saveUpper[iColumn] = newUpper;
3431                                        solver->setColUpper(iColumn, newUpper);
3432                                    }
3433                                }
3434                            }
3435#endif
3436                            // See if integer solution
3437                            feasibleSolution = 
3438                              model->feasibleSolution(choice.numIntInfeasUp,
3439                                                      choice.numObjInfeasUp);
3440                            if (feasibleSolution
3441                                    && model->problemFeasibility()->feasible(model, -1) >= 0) {
3442#ifdef BONMIN
3443                                std::cout << "Node has become integer feasible" << std::endl;
3444                                numberUnsatisfied_ = 0;
3445                                break;
3446#endif
3447                                if (auxiliaryInfo->solutionAddsCuts()) {
3448                                    needHotStartUpdate = true;
3449                                    solver->unmarkHotStart();
3450                                }
3451                                model->setLogLevel(saveLogLevel);
3452                                model->setBestSolution(CBC_STRONGSOL,
3453                                                       newObjectiveValue,
3454                                                       solver->getColSolution()) ;
3455                                if (choice.finishedDown) {
3456                                  double cutoff = model->getCutoff();
3457                                  double downObj = objectiveValue_
3458                                    + choice.downMovement ;
3459                                  if (downObj >= cutoff) {     
3460                                    choice.downMovement = 1.0e100 ;
3461                                    numberStrongInfeasible++;
3462                                }
3463
3464                                }
3465                                if (needHotStartUpdate) {
3466                                    model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3467#ifdef CHECK_DEBUGGER_PATH
3468                                    if ((model->specialOptions()&1) != 0 && onOptimalPath) {
3469                                      const OsiRowCutDebugger *debugger = solver->getRowCutDebugger() ;
3470                                      if (!debugger) {
3471                                        printf("Strong branching up on %d went off optimal path\n",iObject);
3472                                        abort();
3473                                      }
3474                                    }
3475#endif
3476                                    newObjectiveValue = solver->getObjSense() * solver->getObjValue();
3477                                    objectiveValue_ = CoinMax(objectiveValue_,newObjectiveValue);
3478                                    objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_, 0.0);
3479                                    model->feasibleSolution(choice.numIntInfeasDown,
3480                                                            choice.numObjInfeasDown);
3481                                }
3482                                model->setLastHeuristic(NULL);
3483                                model->incrementUsed(solver->getColSolution());
3484                                cutoff = model->getCutoff();
3485                                if (newObjectiveValue >= cutoff) {      //  *new* cutoff
3486                                    objectiveChange = 1.0e100 ;
3487                                    numberStrongInfeasible++;
3488                                }
3489                            }
3490                        }
3491                    } else if (iStatus == 1) {
3492                        choice.finishedUp = true ;
3493                        objectiveChange = COIN_DBL_MAX;
3494                        numberStrongInfeasible++;
3495                    } else {
3496                        // Can't say much as we did not finish
3497                        choice.finishedUp = false ;
3498                        numberUnfinished++;
3499                    }
3500                    choice.upMovement = objectiveChange ;
3501
3502                    // restore bounds
3503                    for ( j = 0; j < numberColumns; j++) {
3504                        if (saveLower[j] != lower[j])
3505                            solver->setColLower(j, saveLower[j]);
3506                        if (saveUpper[j] != upper[j])
3507                            solver->setColUpper(j, saveUpper[j]);
3508                    }
3509                    if (needHotStartUpdate) {
3510                        needHotStartUpdate = false;
3511                        model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3512#ifdef CHECK_DEBUGGER_PATH
3513                        if ((model->specialOptions()&1) != 0 && onOptimalPath) {
3514                          const OsiRowCutDebugger *debugger = solver->getRowCutDebugger() ;
3515                          if (!debugger) {
3516                            printf("Strong branching up on %d went off optimal path\n",iObject);
3517                            abort();
3518                          }
3519                        }
3520#endif
3521                        double newObjValue = solver->getObjSense()*solver->getObjValue();
3522                        objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3523                        //we may again have an integer feasible solution
3524                        int numberIntegerInfeasibilities;
3525                        int numberObjectInfeasibilities;
3526                        if (model->feasibleSolution(
3527                                    numberIntegerInfeasibilities,
3528                                    numberObjectInfeasibilities)) {
3529                            double objValue = solver->getObjValue();
3530                            model->setLogLevel(saveLogLevel);
3531                            model->setBestSolution(CBC_STRONGSOL,
3532                                                   objValue,
3533                                                   solver->getColSolution()) ;
3534                            model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3535#ifdef CHECK_DEBUGGER_PATH
3536                            if ((model->specialOptions()&1) != 0 && onOptimalPath) {
3537                              const OsiRowCutDebugger *debugger = solver->getRowCutDebugger() ;
3538                              if (!debugger) {
3539                                printf("Strong branching up on %d went off optimal path\n",iObject);
3540                                abort();
3541                              }
3542                            }
3543#endif
3544                            double newObjValue = solver->getObjSense()*solver->getObjValue();
3545                            objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3546                            cutoff = model->getCutoff();
3547                        }
3548                        solver->markHotStart();
3549#ifdef RESET_BOUNDS
3550                        memcpy(saveLower,solver->getColLower(),solver->getNumCols()*sizeof(double));
3551                        memcpy(saveUpper,solver->getColUpper(),solver->getNumCols()*sizeof(double));
3552#endif
3553                        if (!solver->isProvenOptimal()) {
3554                          skipAll=-2;
3555                          canSkip = 1;
3556                        }
3557                        xMark++;
3558                    }
3559
3560#if 0 //def DO_ALL_AT_ROOT
3561                    if (strongType)
3562                        printf("Up on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
3563                               choice.objectNumber, iStatus, newObjectiveValue, choice.numItersUp,
3564                               choice.upMovement, choice.finishedUp, choice.numIntInfeasUp,
3565                               choice.numObjInfeasUp);
3566#endif
3567                }
3568
3569                solver->setIntParam(OsiMaxNumIterationHotStart, saveLimit2);
3570                /*
3571                  End of evaluation for this candidate variable. Possibilities are:
3572                  * Both sides below cutoff; this variable is a candidate for branching.
3573                  * Both sides infeasible or above the objective cutoff: no further action
3574                  here. Break from the evaluation loop and assume the node will be purged
3575                  by the caller.
3576                  * One side below cutoff: Install the branch (i.e., fix the variable). Break
3577                  from the evaluation loop and assume the node will be reoptimised by the
3578                  caller.
3579                */
3580                // reset
3581                choice.possibleBranch->resetNumberBranchesLeft();
3582                if (choice.upMovement < 1.0e100) {
3583                    if (choice.downMovement < 1.0e100) {
3584                        // In case solution coming in was odd
3585                        choice.upMovement = CoinMax(0.0, choice.upMovement);
3586                        choice.downMovement = CoinMax(0.0, choice.downMovement);
3587#if ZERO_ONE==2
3588                        // branch on 0-1 first (temp)
3589                        if (fabs(choice.possibleBranch->value()) < 1.0) {
3590                            choice.upMovement *= ZERO_FAKE;
3591                            choice.downMovement *= ZERO_FAKE;
3592                        }
3593#endif
3594                        // feasible - see which best
3595                        if (!canSkip) {
3596                            if (model->messageHandler()->logLevel() > 3)
3597                                printf("sort %g downest %g upest %g ", sort[iDo], downEstimate[iObject],
3598                                       upEstimate[iObject]);
3599                            model->messageHandler()->message(CBC_STRONG, *model->messagesPointer())
3600                            << iObject << iColumn
3601                            << choice.downMovement << choice.numIntInfeasDown
3602                            << choice.upMovement << choice.numIntInfeasUp
3603                            << choice.possibleBranch->value()
3604                            << CoinMessageEol;
3605                        }
3606                        int betterWay=0;
3607                        // If was feasible (extra strong branching) skip
3608                        if (infeasibility) {
3609                            CbcBranchingObject * branchObj =
3610                                dynamic_cast <CbcBranchingObject *>(branch_) ;
3611                            if (branch_)
3612                                assert (branchObj);
3613                            betterWay = decision->betterBranch(choice.possibleBranch,
3614                                                               branchObj,
3615                                                               choice.upMovement,
3616                                                               choice.numIntInfeasUp ,
3617                                                               choice.downMovement,
3618                                                               choice.numIntInfeasDown );
3619                        }
3620                        if (betterWay) {
3621                            // C) create branching object
3622                            if (choiceObject) {
3623                                delete branch_;
3624                                branch_ = choice.possibleBranch->clone();
3625                            } else {
3626                                delete branch_;
3627                                branch_ = choice.possibleBranch;
3628                                choice.possibleBranch = NULL;
3629                            }
3630                            {
3631                                CbcBranchingObject * branchObj =
3632                                    dynamic_cast <CbcBranchingObject *>(branch_) ;
3633                                assert (branchObj);
3634                                //branchObj->way(preferredWay);
3635                                branchObj->way(betterWay);
3636                            }
3637                            bestChoice = choice.objectNumber;
3638                            whichChoice = iDo;
3639                            if (numberStrong <= 1) {
3640                                delete ws;
3641                                ws = NULL;
3642                                break;
3643                            }
3644                        } else {
3645                            if (!choiceObject) {
3646                                delete choice.possibleBranch;
3647                                choice.possibleBranch = NULL;
3648                            }
3649                            if (iDo >= 2*numberStrong) {
3650                                delete ws;
3651                                ws = NULL;
3652                                break;
3653                            }
3654                            if (!dynamicObject || dynamicObject->numberTimesUp() > 1) {
3655                                if (iDo - whichChoice >= numberStrong) {
3656                                    if (!choiceObject) {
3657                                        delete choice.possibleBranch;
3658                                        choice.possibleBranch = NULL;
3659                                    }
3660                                    break; // give up
3661                                }
3662                            } else {
3663                                if (iDo - whichChoice >= 2*numberStrong) {
3664                                    delete ws;
3665                                    ws = NULL;
3666                                    if (!choiceObject) {
3667                                        delete choice.possibleBranch;
3668                                        choice.possibleBranch = NULL;
3669                                    }
3670                                    break; // give up
3671                                }
3672                            }
3673                        }
3674                    } else {
3675                        // up feasible, down infeasible
3676                        anyAction = -1;
3677                        worstFeasible = CoinMax(worstFeasible, choice.upMovement);
3678                        model->messageHandler()->message(CBC_STRONG, *model->messagesPointer())
3679                        << iObject << iColumn
3680                        << choice.downMovement << choice.numIntInfeasDown
3681                        << choice.upMovement << choice.numIntInfeasUp
3682                        << choice.possibleBranch->value()
3683                        << CoinMessageEol;
3684                        //printf("Down infeasible for choice %d sequence %d\n",i,
3685                        // model->object(choice.objectNumber)->columnNumber());
3686                        choice.fix = 1;
3687                        numberToFix++;
3688                        choice.possibleBranch->fix(solver, saveLower, saveUpper, 1);
3689                        if (!choiceObject) {
3690                            delete choice.possibleBranch;
3691                            choice.possibleBranch = NULL;
3692                        } else {
3693                            //choiceObject = new CbcDynamicPseudoCostBranchingObject(*choiceObject);
3694                            choice.possibleBranch = choiceObject;
3695                        }
3696                        assert(doneHotStart);
3697                        solver->unmarkHotStart();
3698                        model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3699#ifdef CHECK_DEBUGGER_PATH
3700                        if ((model->specialOptions()&1) != 0 && onOptimalPath) {
3701                          const OsiRowCutDebugger *debugger = solver->getRowCutDebugger() ;
3702                          if (!debugger) {
3703                            printf("Strong branching down on %d went off optimal path\n",iObject);
3704                            abort();
3705                          }
3706                        }
3707#endif
3708                        double newObjValue = solver->getObjSense()*solver->getObjValue();
3709                        objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3710                        bool goneInfeasible = (!solver->isProvenOptimal()||solver->isDualObjectiveLimitReached());
3711                        solver->markHotStart();
3712#ifdef RESET_BOUNDS
3713                        memcpy(saveLower,solver->getColLower(),solver->getNumCols()*sizeof(double));
3714                        memcpy(saveUpper,solver->getColUpper(),solver->getNumCols()*sizeof(double));
3715#endif
3716                        if (!solver->isProvenOptimal()) {
3717                          skipAll=-2;
3718                          canSkip = 1;
3719                        }
3720                        xMark++;
3721                        // may be infeasible (if other way stopped on iterations)
3722                        if (goneInfeasible) {
3723                            // neither side feasible
3724                            anyAction = -2;
3725                            if (!choiceObject) {
3726                                delete choice.possibleBranch;
3727                                choice.possibleBranch = NULL;
3728                            }
3729                            //printf("Both infeasible for choice %d sequence %d\n",i,
3730                            // model->object(choice.objectNumber)->columnNumber());
3731                            delete ws;
3732                            ws = NULL;
3733                            break;
3734                        }
3735                    }
3736                } else {
3737                    if (choice.downMovement < 1.0e100) {
3738                        // down feasible, up infeasible
3739                        anyAction = -1;
3740                        worstFeasible = CoinMax(worstFeasible, choice.downMovement);
3741                        model->messageHandler()->message(CBC_STRONG, *model->messagesPointer())
3742                        << iObject << iColumn
3743                        << choice.downMovement << choice.numIntInfeasDown
3744                        << choice.upMovement << choice.numIntInfeasUp
3745                        << choice.possibleBranch->value()
3746                        << CoinMessageEol;
3747                        choice.fix = -1;
3748                        numberToFix++;
3749                        choice.possibleBranch->fix(solver, saveLower, saveUpper, -1);
3750                        if (!choiceObject) {
3751                            delete choice.possibleBranch;
3752                            choice.possibleBranch = NULL;
3753                        } else {
3754                            //choiceObject = new CbcDynamicPseudoCostBranchingObject(*choiceObject);
3755                            choice.possibleBranch = choiceObject;
3756                        }
3757                        assert(doneHotStart);
3758                        solver->unmarkHotStart();
3759                        model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3760#ifdef CHECK_DEBUGGER_PATH
3761                        if ((model->specialOptions()&1) != 0 && onOptimalPath) {
3762                          const OsiRowCutDebugger *debugger = solver->getRowCutDebugger() ;
3763                          if (!debugger) {
3764                            printf("Strong branching down on %d went off optimal path\n",iObject);
3765                            solver->writeMps("query");
3766                            abort();
3767                          }
3768                        }
3769#endif
3770                        double newObjValue = solver->getObjSense()*solver->getObjValue();
3771                        objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3772                        bool goneInfeasible = (!solver->isProvenOptimal()||solver->isDualObjectiveLimitReached());
3773                        solver->markHotStart();
3774#ifdef RESET_BOUNDS
3775                        memcpy(saveLower,solver->getColLower(),solver->getNumCols()*sizeof(double));
3776                        memcpy(saveUpper,solver->getColUpper(),solver->getNumCols()*sizeof(double));
3777#endif
3778                        if (!solver->isProvenOptimal()) {
3779                          skipAll=-2;
3780                          canSkip = 1;
3781                        }
3782                        xMark++;
3783                        // may be infeasible (if other way stopped on iterations)
3784                        if (goneInfeasible) {
3785                            // neither side feasible
3786                            anyAction = -2;
3787                            if (!choiceObject) {
3788                                delete choice.possibleBranch;
3789                                choice.possibleBranch = NULL;
3790                            }
3791                            delete ws;
3792                            ws = NULL;
3793                            break;
3794                        }
3795                    } else {
3796                        // neither side feasible
3797                        anyAction = -2;
3798                        if (!choiceObject) {
3799                            delete choice.possibleBranch;
3800                            choice.possibleBranch = NULL;
3801                        }
3802                        delete ws;
3803                        ws = NULL;
3804                        break;
3805                    }
3806                }
3807                // Check max time
3808                hitMaxTime = (model->getCurrentSeconds() >
3809                               model->getDblParam(CbcModel::CbcMaximumSeconds));
3810                if (hitMaxTime) {
3811                    // make sure rest are fast
3812                    for ( int jDo = iDo + 1; jDo < numberToDo; jDo++) {
3813                        int iObject = whichObject[iDo];
3814                        OsiObject * object = model->modifiableObject(iObject);
3815                        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3816                            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3817                        if (dynamicObject)
3818                            dynamicObject->setNumberBeforeTrust(0);
3819                    }
3820                    numberTest = 0;
3821                }
3822                if (!choiceObject) {
3823                    delete choice.possibleBranch;
3824                }
3825            }
3826            if (model->messageHandler()->logLevel() > 3) {
3827                if (anyAction == -2) {
3828                    printf("infeasible\n");
3829                } else if (anyAction == -1) {
3830                    printf("%d fixed AND choosing %d iDo %d iChosenWhen %d numberToDo %d\n", numberToFix, bestChoice,
3831                           iDo, whichChoice, numberToDo);
3832                } else {
3833                    int iObject = whichObject[whichChoice];
3834                    OsiObject * object = model->modifiableObject(iObject);
3835                    CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3836                        dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3837                    if (dynamicObject) {
3838                        int iColumn = dynamicObject->columnNumber();
3839                        printf("choosing %d (column %d) iChosenWhen %d numberToDo %d\n", bestChoice,
3840                               iColumn, whichChoice, numberToDo);
3841                    }
3842                }
3843            }
3844            if (doneHotStart) {
3845                // Delete the snapshot
3846                solver->unmarkHotStart();
3847                // back to normal
3848                solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
3849                // restore basis
3850                solver->setWarmStart(ws);
3851            }
3852            solver->setIntParam(OsiMaxNumIterationHotStart, saveLimit);
3853            // Unless infeasible we will carry on
3854            // But we could fix anyway
3855            if (numberToFix && !hitMaxTime) {
3856                if (anyAction != -2) {
3857                    // apply and take off
3858                    bool feasible = true;
3859                    // can do quick optimality check
3860                    int easy = 2;
3861                    solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, &easy) ;
3862                    model->resolve(NULL, 11, saveSolution, saveLower, saveUpper) ;
3863#ifdef CHECK_DEBUGGER_PATH
3864                    if ((model->specialOptions()&1) != 0 && onOptimalPath) {
3865                      const OsiRowCutDebugger *debugger = solver->getRowCutDebugger() ;
3866                      if (!debugger) {
3867                        printf("Strong branching went off optimal path\n");
3868                        abort();
3869                      }
3870                    }
3871#endif
3872                    double newObjValue = solver->getObjSense()*solver->getObjValue();
3873                    objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3874                    solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
3875                    feasible = solver->isProvenOptimal();
3876                    if (feasible) {
3877                        anyAction = 0;
3878                    } else {
3879                        anyAction = -2;
3880                        finished = true;
3881                    }
3882                }
3883            }
3884            // If  fixed then round again
3885            // See if candidate still possible
3886            if (branch_) {
3887                 const OsiObject * object = model->object(bestChoice);
3888                 double infeasibility = object->checkInfeasibility(&usefulInfo);
3889                 if (!infeasibility) {
3890                   // take out
3891                   delete branch_;
3892                   branch_ = NULL;
3893                 } else {
3894                   // get preferred way
3895                   int preferredWay;
3896                   object->infeasibility(&usefulInfo, preferredWay);
3897                   CbcBranchingObject * branchObj =
3898                     dynamic_cast <CbcBranchingObject *>(branch_) ;
3899                   assert (branchObj);
3900                   branchObj->way(preferredWay);
3901#ifdef CBCNODE_TIGHTEN_BOUNDS
3902                   bool fixed = branchObj->tighten(solver);
3903                   if (fixed) {
3904                     //printf("Variable now fixed!\n");
3905                     // take out
3906                     delete branch_;
3907                     branch_ = NULL;
3908                   }
3909#endif
3910                 }
3911            }
3912            if (!branch_ && anyAction != -2 && !hitMaxTime) {
3913                finished = false;
3914            }
3915            delete ws;
3916        }
3917    }
3918    // update number of strong iterations etc
3919    model->incrementStrongInfo(numberStrongDone, numberStrongIterations,
3920                               anyAction == -2 ? 0 : numberToFix, anyAction == -2);
3921    if (model->searchStrategy() == -1) {
3922#ifndef COIN_DEVELOP
3923        if (solver->messageHandler()->logLevel() > 1)
3924#endif
3925            printf("%d strong, %d iters, %d inf, %d not finished, %d not trusted\n",
3926                   numberStrongDone, numberStrongIterations, numberStrongInfeasible, numberUnfinished,
3927                   numberNotTrusted);
3928        // decide what to do
3929        int strategy = 1;
3930        if (((numberUnfinished*4 > numberStrongDone &&
3931                numberStrongInfeasible*40 < numberStrongDone) ||
3932                numberStrongInfeasible < 0) && model->numberStrong() < 10 && model->numberBeforeTrust() <= 20 && model->numberObjects() > CoinMax(1000, solver->getNumRows())) {
3933            strategy = 2;
3934#ifdef COIN_DEVELOP
3935            //if (model->logLevel()>1)
3936            printf("going to strategy 2\n");
3937#endif
3938            // Weaken
3939            model->setNumberStrong(2);
3940            model->setNumberBeforeTrust(1);
3941            model->synchronizeNumberBeforeTrust();
3942        }
3943        if (numberNodes)
3944            strategy = 1;  // should only happen after hot start
3945        model->setSearchStrategy(strategy);
3946    } else if (numberStrongDone) {
3947        //printf("%d strongB, %d iters, %d inf, %d not finished, %d not trusted\n",
3948        //   numberStrongDone,numberStrongIterations,numberStrongInfeasible,numberUnfinished,
3949        //   numberNotTrusted);
3950    }
3951    if (model->searchStrategy() == 1 && numberNodes > 500 && numberNodes < -510) {
3952#ifndef COIN_DEVELOP
3953        if (solver->messageHandler()->logLevel() > 1)
3954#endif
3955            printf("after %d nodes - %d strong, %d iters, %d inf, %d not finished, %d not trusted\n",
3956                   numberNodes, numberStrongDone, numberStrongIterations, numberStrongInfeasible, numberUnfinished,
3957                   numberNotTrusted);
3958        // decide what to do
3959        if (numberUnfinished*10 > numberStrongDone + 1 ||
3960                !numberStrongInfeasible) {
3961          COIN_DETAIL_PRINT(printf("going to strategy 2\n"));
3962            // Weaken
3963            model->setNumberStrong(2);
3964            model->setNumberBeforeTrust(1);
3965            model->synchronizeNumberBeforeTrust();
3966            model->setSearchStrategy(2);
3967        }
3968    }
3969    if (numberUnfinished*10 < numberStrongDone &&
3970        model->numberStrongIterations()*20 < model->getIterationCount()&&
3971                                !auxiliaryInfo->solutionAddsCuts()) {
3972        //printf("increasing trust\n");
3973        model->synchronizeNumberBeforeTrust(2);
3974    }
3975
3976    // Set guessed solution value
3977    guessedObjectiveValue_ = objectiveValue_ + estimatedDegradation;
3978    int kColumn=-1;
3979    if (branch_) {
3980      CbcObject * obj = (dynamic_cast<CbcBranchingObject *>(branch_))->object(); 
3981      CbcSimpleInteger * branchObj =
3982        dynamic_cast <CbcSimpleInteger *>(obj) ;
3983      if (branchObj) {
3984        kColumn=branchObj->columnNumber();
3985      }
3986    }
3987#ifdef COIN_HAS_NTY
3988    if (orbitOption&&kColumn>=0) {
3989      CbcSymmetry * symmetryInfo = model->symmetryInfo();
3990      CbcNodeInfo * infoX = lastNode ? lastNode->nodeInfo() : NULL;
3991      bool worthTrying = false;
3992      if (infoX) {
3993        CbcNodeInfo * info = infoX;
3994        for (int i=0;i<NTY_BAD_DEPTH;i++) {
3995          if (!info->parent()) {
3996            worthTrying = true;
3997            break;
3998          }
3999          info = info->parent();
4000          if (info->symmetryWorked()) {
4001            worthTrying = true;
4002            break;
4003          }
4004        }
4005      } else {
4006        worthTrying=true;
4007      }
4008      if (orbitOption==3&&depth_>5)
4009        worthTrying=false;
4010      if (symmetryInfo && worthTrying) {
4011        if ((orbitOption&1)==1) {
4012          symmetryInfo->ChangeBounds(solver->getColLower(),
4013                                     solver->getColUpper(),
4014                                     solver->getNumCols(),false);
4015          symmetryInfo->Compute_Symmetry();
4016          symmetryInfo->fillOrbits();
4017        }
4018        const int * orbits = symmetryInfo->whichOrbit();
4019        if (orbits && orbits[kColumn]>=0) {
4020          int numberUsefulOrbits = symmetryInfo->numberUsefulOrbits();
4021          if (solver->messageHandler()->logLevel() > 1)
4022            printf("Orbital Branching on %d - way %d n %d\n",kColumn,way(),numberUsefulOrbits);
4023          if (numberUsefulOrbits<1000||orbitOption==3) {
4024            delete branch_;
4025            branch_ = new CbcOrbitalBranchingObject(model,kColumn,1,0,NULL);
4026            if (infoX)
4027              infoX->setSymmetryWorked();
4028          }
4029        }
4030      }
4031    }
4032#endif
4033    if (model->logLevel()>1)
4034    printf ("Node %d depth %d unsatisfied %d sum %g obj %g guess %g branching on %d\n",
4035            model->getNodeCount(),depth_,numberUnsatisfied_,
4036            sumInfeasibilities_,objectiveValue_,guessedObjectiveValue_,
4037            kColumn);
4038#ifdef DO_ALL_AT_ROOT
4039    if (strongType) {
4040      char general[200];
4041      if (strongType==1)
4042        sprintf(general,"Strong branching on all %d unsatisfied, %d iterations (depth %d)\n",
4043                saveNumberToDo,numberStrongIterations,depth_);
4044      else
4045        sprintf(general,"Strong branching on all %d unfixed variables (%d unsatisfied), %d iterations (depth %d)\n",
4046                saveNumberToDo+saveSatisfiedVariables,saveNumberToDo,numberStrongIterations,depth_);
4047      model->messageHandler()->message(CBC_FPUMP2,model->messages())
4048        << general << CoinMessageEol ;
4049    }
4050#endif
4051#ifdef DEBUG_SOLUTION
4052    if(onOptimalPath&&anyAction==-2) {
4053      printf("Gone off optimal path in CbcNode\n");
4054      assert(!onOptimalPath||anyAction!=-2);
4055    }
4056#endif
4057    /*
4058      Cleanup, then we're finished
4059    */
4060    if (!model->branchingMethod())
4061        delete decision;
4062
4063    delete choiceObject;
4064    delete [] sort;
4065    delete [] whichObject;
4066#ifdef RANGING
4067    delete [] objectMark;
4068#endif
4069    delete [] saveLower;
4070    delete [] saveUpper;
4071    delete [] upEstimate;
4072    delete [] downEstimate;
4073# ifdef COIN_HAS_CLP
4074    if (osiclp) {
4075        osiclp->setSpecialOptions(saveClpOptions);
4076    }
4077# endif
4078    // restore solution
4079    solver->setColSolution(saveSolution);
4080    model->reserveCurrentSolution(saveSolution);
4081    delete [] saveSolution;
4082    model->setStateOfSearch(saveStateOfSearch);
4083    model->setLogLevel(saveLogLevel);
4084    // delete extra regions
4085    if (usefulInfo.usefulRegion_) {
4086        delete [] usefulInfo.usefulRegion_;
4087        delete [] usefulInfo.indexRegion_;
4088        delete [] usefulInfo.pi_;
4089        usefulInfo.usefulRegion_ = NULL;
4090        usefulInfo.indexRegion_ = NULL;
4091        usefulInfo.pi_ = NULL;
4092    }
4093    useShadow = model->moreSpecialOptions() & 7;
4094    if ((useShadow == 5 && model->getSolutionCount()) || useShadow == 6) {
4095        // zap pseudo shadow prices
4096        model->pseudoShadow(-1);
4097        // and switch off
4098        model->setMoreSpecialOptions(model->moreSpecialOptions()&(~1023));
4099    }
4100    return anyAction;
4101}
4102// 0 is down, 1 is up
4103typedef struct {
4104  double initialValue; // initial value
4105  double upLowerBound; // Lower bound when going up
4106  double downUpperBound; // Upper bound when going down
4107  double movement[2]; // cost  (and initial away from feasible)
4108  double sumModified[2]; // Sum of integer changes
4109  int modified[2]; // Number integers changed
4110  int numIntInfeas[2]; // without odd ones
4111  int numObjInfeas[2]; // just odd ones
4112  bool finished[2]; // true if solver finished
4113  int numIters[2]; // number of iterations in solver (-1 if never solved)
4114  double * integerSolution; // output if thinks integer solution
4115# ifdef COIN_HAS_CLP
4116  ClpDualRowSteepest * steepest;
4117#endif
4118  int columnNumber; // Which column it is
4119} StrongInfo;
4120typedef struct {
4121  double integerTolerance;
4122  double * originalSolution;
4123  CoinWarmStart * ws;
4124  double * newObjective;
4125# ifdef COIN_HAS_CLP
4126  ClpDualRowSteepest * dualRowPivot;
4127  ClpPrimalColumnPivot * primalColumnPivot;
4128# endif
4129  int * back;
4130  int solveType;
4131} StrongStaticInfo;
4132typedef struct {
4133  StrongStaticInfo *staticInfo;
4134  StrongInfo * choice;
4135  OsiSolverInterface * solver;
4136  double * tempSolution;
4137  CoinWarmStart * tempBasis;
4138  int whichChoice;
4139} StrongBundle;
4140/* return 1 if possible solution (for solveType 100 if infeasible)
4141   2 set if down was infeasible
4142   4 set if up was infeasible
4143 */
4144int solveAnalyze(void * info) {
4145  StrongBundle * bundle = reinterpret_cast<StrongBundle *>(info);
4146  StrongInfo * choice = bundle->choice;
4147  StrongStaticInfo * staticInfo = bundle->staticInfo;
4148  OsiSolverInterface * solver = bundle->solver;
4149  int solveType = staticInfo->solveType;
4150  if (solveType==77) {
4151    return 0;
4152  }
4153  const double * saveSolution = staticInfo->originalSolution;
4154  int iColumn = choice->columnNumber;
4155  const int * back = staticInfo->back;
4156  double newObjectiveValue = 1.0e100;
4157  double integerTolerance = staticInfo->integerTolerance;
4158  double bestSolutionValue=COIN_DBL_MAX;
4159  int returnStatus=0;
4160  // status is 0 finished, 1 infeasible and other
4161  int iStatus;
4162  /*
4163    Try the down direction first. (Specify the initial branching alternative as
4164    down with a call to way(-1). Each subsequent call to branch() performs the
4165    specified branch and advances the branch object state to the next branch
4166    alternative.)
4167  */
4168  for (int iWay=0;iWay<2;iWay++) {
4169    if (choice->numIters[iWay]==0) {
4170      int numberColumns=solver->getNumCols();
4171      if (solveType!=100) {
4172        double saveBound;
4173        if (iWay==0) {
4174          saveBound = solver->getColUpper()[iColumn];
4175          solver->setColUpper(iColumn,choice->downUpperBound);
4176        } else {
4177          saveBound = solver->getColLower()[iColumn];
4178          solver->setColLower(iColumn,choice->upLowerBound);
4179        }
4180        if ((solveType&2)==0) {
4181          solver->solveFromHotStart() ;
4182        } else {
4183          // restore basis
4184          solver->setWarmStart(staticInfo->ws);
4185# ifdef COIN_HAS_CLP
4186          if (staticInfo->dualRowPivot) {
4187            OsiClpSolverInterface * osiclp = dynamic_cast<OsiClpSolverInterface *>(solver);
4188            ClpSimplex * simplex = osiclp->getModelPtr();
4189            simplex->setDualRowPivotAlgorithm(*staticInfo->dualRowPivot);
4190            //simplex->dualRowPivot()->saveWeights(simplex,4);
4191            simplex->setWhatsChanged(ALL_SAME_EXCEPT_COLUMN_BOUNDS);
4192            simplex->dual(0,5);
4193          } else {
4194#endif
4195          solver->resolve();
4196# ifdef COIN_HAS_CLP
4197          }
4198#endif
4199        }
4200        if (iWay==0)
4201          solver->setColUpper(iColumn,saveBound);
4202        else
4203          solver->setColLower(iColumn,saveBound);
4204        /*
4205          We now have an estimate of objective degradation that we can use for strong
4206          branching. If we're over the cutoff, the variable is monotone up.
4207          If we actually made it to optimality, check for a solution, and if we have
4208          a good one, call setBestSolution to process it. Note that this may reduce the
4209          cutoff, so we check again to see if we can declare this variable monotone.
4210        */
4211        if (solver->isProvenOptimal()) {
4212          iStatus = 0; // optimal
4213        } else if (solver->isIterationLimitReached()
4214                   && !solver->isDualObjectiveLimitReached()) {
4215          iStatus = 2; // unknown
4216        } else {
4217          iStatus = 1; // infeasible
4218        }
4219        newObjectiveValue = solver->getObjSense() * solver->getObjValue();
4220        choice->numIters[iWay] = solver->getIterationCount();
4221        // Look at interaction
4222        const double * thisSolution = solver->getColSolution();
4223        int numberModified=0;
4224        double sumModified=0.0;
4225        int numberInfeas=0;
4226        for (int i=0;i<numberColumns;i++) {
4227          if (back[i]>=0) {
4228            double value = thisSolution[i];
4229            if (iColumn!=i) {
4230              double difference = fabs(saveSolution[i]-value);
4231              if (difference>integerTolerance) {
4232                numberModified++;
4233                sumModified += difference;
4234              }
4235            }
4236            if (fabs(value-floor(value+0.5))>integerTolerance)
4237              numberInfeas++;;
4238          }
4239        }
4240        choice->numIntInfeas[iWay]=numberInfeas;
4241        choice->sumModified[iWay] = sumModified;
4242        choice->modified[iWay] = numberModified;
4243        if (!iStatus) {
4244          choice->finished[iWay] = true ;
4245          if (!numberInfeas) {
4246            returnStatus=1;
4247            if (!choice->integerSolution) {
4248              bestSolutionValue=newObjectiveValue;
4249              choice->integerSolution=CoinCopyOfArray(thisSolution,numberColumns);;
4250            } else if (bestSolutionValue>newObjectiveValue) {
4251              memcpy(choice->integerSolution,thisSolution,numberColumns*sizeof(double));
4252            }
4253          }
4254        } else if (iStatus == 1) {
4255          newObjectiveValue = 1.0e100 ;
4256        } else {
4257          // Can't say much as we did not finish
4258          choice->finished[iWay] = false ;
4259        }
4260        choice->movement[iWay] = newObjectiveValue ;
4261      } else {
4262# ifdef COIN_HAS_CLP
4263        OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
4264        ClpSimplex * simplex = osiclp ? osiclp->getModelPtr() : NULL;
4265#endif
4266        // doing continuous and general integer
4267        solver->setColSolution(staticInfo->originalSolution);
4268        solver->setWarmStart(staticInfo->ws);
4269        double saveBound;
4270        double newBound;
4271        if (iWay==0) {
4272          saveBound=solver->getColUpper()[iColumn];
4273          solver->setColUpper(iColumn,choice->downUpperBound);
4274          newBound=choice->downUpperBound;
4275        } else {
4276          saveBound=solver->getColLower()[iColumn];
4277          solver->setColLower(iColumn,choice->upLowerBound);
4278          newBound=choice->upLowerBound;
4279        }
4280# if 0 //def COIN_HAS_CLP
4281        if (simplex) {
4282          // set solution to new bound (if basic will be recomputed)
4283          simplex->primalColumnSolution()[iColumn]=newBound;
4284        }
4285#endif
4286        solver->setHintParam(OsiDoDualInResolve, true, OsiHintDo) ;
4287#define PRINT_ANALYZE  0
4288#if PRINT_ANALYZE>0
4289        osiclp->getModelPtr()->setLogLevel(1);
4290        solver->setHintParam(OsiDoReducePrint, false, OsiHintTry);
4291#endif
4292        solver->resolve();
4293        if (iWay==0) {
4294#if PRINT_ANALYZE>0
4295          printf("column %d down original %g <= %g <= %g upper now %g - result %s\n",
4296                 iColumn,solver->getColLower()[iColumn],
4297                 staticInfo->originalSolution[iColumn],saveBound,
4298                 newBound,solver->isProvenOptimal() ? "ok" : "infeas");
4299#endif
4300          solver->setColUpper(iColumn,saveBound);
4301        } else {
4302#if PRINT_ANALYZE>0
4303          printf("column %d up original %g <= %g <= %g lower now %g - result %s\n",
4304                 iColumn,saveBound,staticInfo->originalSolution[iColumn],
4305                 solver->getColUpper()[iColumn],
4306                 newBound,solver->isProvenOptimal() ? "ok" : "infeas");
4307#endif
4308          solver->setColLower(iColumn,saveBound);
4309        }
4310        choice->numIters[iWay] = solver->getIterationCount();
4311        if (solver->isProvenOptimal()) {
4312          //printf("Way %d - all way %d iterations - column %d\n",
4313          //     iWay,solver->getIterationCount(),iColumn);
4314          // can go all way
4315          choice->movement[iWay] = newBound;
4316        } else {
4317          // zero objective
4318          double offset;
4319          solver->getDblParam(OsiObjOffset,offset);
4320          solver->setDblParam(OsiObjOffset, 0.0);
4321          solver->setObjective(staticInfo->newObjective+numberColumns);
4322          if (iWay==0) {
4323            solver->setObjCoeff(iColumn,1.0);
4324          } else {
4325            solver->setObjCoeff(iColumn,-1.0);
4326          }
4327          solver->setColSolution(staticInfo->originalSolution);
4328          solver->setWarmStart(staticInfo->ws);
4329          solver->setHintParam(OsiDoDualInResolve, false, OsiHintDo) ;
4330          solver->resolve();
4331          //printf("Way %d - first solve %d iterations, second %d - column %d\n",
4332          //     iWay,choice->numIters[iWay],solver->getIterationCount(),iColumn);
4333          choice->movement[iWay] = solver->getColSolution()[iColumn];
4334          choice->numIters[iWay] += solver->getIterationCount();
4335#if PRINT_ANALYZE>0
4336          if (iWay==0) {
4337            printf("column %d down can get to %g - result %s\n",
4338                 iColumn,solver->getColSolution()[iColumn],solver->isProvenOptimal() ? "ok" : "infeas");
4339          } else {
4340            printf("column %d up can get to %g - result %s\n",
4341                 iColumn,solver->getColSolution()[iColumn],solver->isProvenOptimal() ? "ok" : "infeas");
4342          }
4343#endif
4344          // reset objective
4345          solver->setDblParam(OsiObjOffset, offset);
4346          solver->setObjective(staticInfo->newObjective);
4347          if (!solver->isProvenOptimal()) {
4348# ifdef COIN_HAS_CLP
4349            OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
4350            ClpSimplex * simplex = osiclp->getModelPtr();
4351            double sum = simplex->sumPrimalInfeasibilities();
4352            sum /= static_cast<double>(simplex->numberPrimalInfeasibilities());
4353            if (sum>1.0e-3) {
4354#endif
4355            choice->modified[0]=1;
4356            returnStatus=1;
4357            solver->writeMps("bad","mps");
4358            abort();
4359# ifdef COIN_HAS_CLP
4360            }
4361#endif
4362          }
4363        }
4364        //solver->setObjCoeff(iColumn,0.0);
4365      }
4366    }
4367  }
4368  return returnStatus;
4369}
4370#ifdef THREADS_IN_ANALYZE
4371void * cbc_parallelManager(void * stuff)
4372{
4373  CoinPthreadStuff * driver = reinterpret_cast<CoinPthreadStuff *>(stuff);
4374  int whichThread=driver->whichThread();
4375  CoinThreadInfo * threadInfo = driver->threadInfoPointer(whichThread);
4376  threadInfo->status=-1;
4377  int * which = threadInfo->stuff;
4378  pthread_barrier_wait(driver->barrierPointer());
4379#if 0
4380  int status=-1;
4381  while (status!=100)
4382    status=timedWait(driver,1000,2);
4383  pthread_cond_signal(driver->conditionPointer(1));
4384  pthread_mutex_unlock(driver->mutexPointer(1,whichThread));
4385#endif
4386  // so now mutex_ is locked
4387  int whichLocked=0;
4388  while (true) {
4389    pthread_mutex_t * mutexPointer = driver->mutexPointer(whichLocked,whichThread);
4390    // wait
4391    //printf("Child waiting for %d - status %d %d %d\n",
4392    //     whichLocked,lockedX[0],lockedX[1],lockedX[2]);
4393#ifdef DETAIL_THREAD
4394    printf("thread %d about to lock mutex %d\n",whichThread,whichLocked);
4395#endif
4396    pthread_mutex_lock (mutexPointer);
4397    whichLocked++;
4398    if (whichLocked==3)
4399      whichLocked=0;
4400    int unLock=whichLocked+1;
4401    if (unLock==3)
4402      unLock=0;
4403    //printf("child pointer %p status %d\n",threadInfo,threadInfo->status);
4404    assert(threadInfo->status>=0);
4405    if (threadInfo->status==1000)
4406      pthread_exit(NULL);
4407    int type=threadInfo->status;
4408    int & returnCode=which[0];
4409    int iPass=which[1];
4410    //CoinIndexedVector * array;
4411    //double dummy;
4412    switch(type) {
4413      // dummy
4414    case 0:
4415      break;
4416    case 1:
4417      returnCode=solveAnalyze(threadInfo->extraInfo);
4418      threadInfo->stuff[3]=0;
4419      break;
4420    case 100:
4421      // initialization
4422      break;
4423    }
4424    threadInfo->status= (type!=1) ? -1 : -2;
4425#ifdef DETAIL_THREAD
4426    printf("thread %d about to unlock mutex %d\n",whichThread,unLock);
4427#endif
4428    pthread_mutex_unlock (driver->mutexPointer(unLock,whichThread));
4429  }
4430}
4431#endif
4432int CbcNode::analyze (CbcModel *model, double * results)
4433{
4434#define COIN_DETAIL
4435  int i;
4436  int numberIterationsAllowed = model->numberAnalyzeIterations();
4437  int numberColumns = model->getNumCols();
4438  int numberRows = model->getNumRows();
4439  int numberObjects = model->numberObjects();
4440  int numberIntegers = model->numberIntegers();
4441  int numberLookIntegers=0;
4442  int highestPriority=COIN_INT_MAX;
4443  int * back = new int[numberColumns];
4444  const int * integerVariable = model->integerVariable();
4445  for (i = 0; i < numberIntegers; i++) {
4446    highestPriority = CoinMin(highestPriority,model->modifiableObject(i)->priority());
4447  }
4448  for (i = 0; i < numberColumns; i++) 
4449    back[i] = -1;
4450  for (i = 0; i < numberIntegers; i++) {
4451    int iColumn = integerVariable[i];
4452    back[iColumn] = i;
4453    if (model->modifiableObject(i)->priority()==highestPriority) {
4454      numberLookIntegers++;
4455    } else {
4456      back[iColumn] = i+numberColumns;
4457    }
4458  }
4459  /*
4460    0 - just look
4461    0 (1) bit - use to set priorities
4462    1 (2) bit - look at bounds on all variables and more iterations
4463    2 (4) bit - do threaded (if parallelMode()==1 then not repeatable if any fixed)
4464    3 (8) bit -
4465    4 (16) bit - do even if m*n>1,000,000
4466    5 (32) bit - printing time
4467    6 (64) bit - save mps file
4468  */
4469  int solveType;
4470  char general[200];
4471  if (numberIterationsAllowed>0) {
4472    solveType = 0;
4473  } else {
4474    solveType = - numberIterationsAllowed;
4475    if ((solveType&16)==0) {
4476      double size=numberRows;
4477      size*=numberLookIntegers;
4478      if (size>1000000) {
4479        if ((solveType&32)!=0)
4480          model->messageHandler()->message(CBC_GENERAL, *model->messagesPointer())
4481            << "Skipping analyze as problem too large"
4482            << CoinMessageEol;
4483        return 0;
4484      }
4485    }
4486    sprintf(general,"Analyze options %d",solveType);
4487    model->messageHandler()->message(CBC_GENERAL, *model->messagesPointer())
4488      << general
4489      << CoinMessageEol;
4490    if ((solveType&1)!=0)
4491      model->messageHandler()->message(CBC_GENERAL, *model->messagesPointer())
4492        << "Using to set priorities (probably bad idea)"
4493        << CoinMessageEol;
4494    if ((solveType&2)!=0)
4495      model->messageHandler()->message(CBC_GENERAL, *model->messagesPointer())
4496        << "Use more iterations and look at continuous/general integer variables"
4497        << CoinMessageEol;
4498    if ((solveType&4)!=0)
4499      model->messageHandler()->message(CBC_GENERAL, *model->messagesPointer())
4500        << "Use threads"
4501        << CoinMessageEol;
4502    if ((solveType&32)!=0)
4503      model->messageHandler()->message(CBC_GENERAL, *model->messagesPointer())
4504        << "32 switches on more printing, (16 bit allows large problems)"
4505        << CoinMessageEol;
4506  }
4507  OsiSolverInterface * solver = model->solver();
4508  objectiveValue_ = solver->getObjSense() * solver->getObjValue();
4509  const double * lower = solver->getColLower();
4510  const double * upper = solver->getColUpper();
4511  const double * dj = solver->getReducedCost();
4512  // What results is
4513  double * newLower = results;
4514  double * objLower = newLower + numberIntegers;
4515  double * newUpper = objLower + numberIntegers;
4516  double * objUpper = newUpper + numberIntegers;
4517  double * interAction = objUpper + numberIntegers;
4518  for (i = 0; i < numberIntegers; i++) {
4519    int iColumn = integerVariable[i];
4520    newLower[i] = lower[iColumn];
4521    objLower[i] = -COIN_DBL_MAX;
4522    newUpper[i] = upper[iColumn];
4523    objUpper[i] = -COIN_DBL_MAX;
4524    interAction[i] = 0.0;
4525  }
4526  double * objMovement=new double[2*numberIntegers];
4527  memset(objMovement,0,2*numberIntegers*sizeof(double));
4528  double * saveUpper = new double[numberColumns];
4529  double * saveLower = new double[numberColumns];
4530  // Save solution in case heuristics need good solution later
4531 
4532  double * saveSolution = new double[numberColumns];
4533  memcpy(saveSolution, solver->getColSolution(), numberColumns*sizeof(double));
4534  model->reserveCurrentSolution(saveSolution);
4535  for (i = 0; i < numberColumns; i++) {
4536    saveLower[i] = lower[i];
4537    saveUpper[i] = upper[i];
4538  }
4539  // Get arrays to sort
4540  double * sort = new double[numberObjects];
4541  int * whichObject = new int[numberObjects];
4542  int numberToFix = 0;
4543  int numberToDo = 0;
4544  double integerTolerance =
4545    model->getDblParam(CbcModel::CbcIntegerTolerance);
4546  // point to useful information
4547  OsiBranchingInformation usefulInfo = model->usefulInformation();
4548  // and modify
4549  usefulInfo.depth_ = depth_;
4550 
4551  // compute current state
4552  int numberObjectInfeasibilities; // just odd ones
4553  int numberIntegerInfeasibilities;
4554  model->feasibleSolution(
4555                          numberIntegerInfeasibilities,
4556                          numberObjectInfeasibilities);
4557  if (solveType) {
4558    if ((solveType&2)==0)
4559      numberIterationsAllowed=200*numberIntegerInfeasibilities;
4560    else
4561      numberIterationsAllowed=COIN_INT_MAX;
4562  }
4563  int saveAllowed=numberIterationsAllowed;
4564# ifdef COIN_HAS_CLP
4565  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
4566  int saveClpOptions = 0;
4567  bool fastIterations = (model->specialOptions() & 8) != 0;
4568  if (osiclp) {
4569    saveClpOptions = osiclp->specialOptions();
4570    // for faster hot start
4571    if (fastIterations) 
4572      osiclp->setSpecialOptions(saveClpOptions | 8192);
4573    else
4574      osiclp->setSpecialOptions(saveClpOptions | 2048); // switch off crunch
4575  }
4576# else
4577  bool fastIterations = false ;
4578# endif
4579  /*
4580    Scan for branching objects that indicate infeasibility.
4581   
4582    The algorithm is to fill the array with a set of good candidates (by
4583    infeasibility).
4584   
4585  */
4586  numberToDo = 0;
4587  for (i = 0; i < numberObjects; i++) {
4588    OsiObject * object = model->modifiableObject(i);
4589    CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
4590      dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
4591    if (!dynamicObject)
4592      continue;
4593    if (dynamicObject->priority()!=highestPriority)
4594      continue;
4595    double infeasibility = object->checkInfeasibility(&usefulInfo);
4596    int iColumn = dynamicObject->columnNumber();
4597    if (saveUpper[iColumn] == saveLower[iColumn])
4598      continue;
4599    if (infeasibility)
4600      sort[numberToDo] = -1.0e10 - infeasibility;
4601    else
4602      sort[numberToDo] = -fabs(dj[iColumn]);
4603    whichObject[numberToDo++] = i;
4604  }
4605  // Save basis
4606  CoinWarmStart * ws = solver->getWarmStart();
4607  int saveLimit;
4608  solver->getIntParam(OsiMaxNumIterationHotStart, saveLimit);
4609  int targetIterations = CoinMax(500, numberIterationsAllowed / numberObjects);
4610  if (saveLimit < targetIterations)
4611    solver->setIntParam(OsiMaxNumIterationHotStart, targetIterations);
4612  if ((solveType&2)==0) {
4613    // Mark hot start
4614    solver->markHotStart();
4615  }
4616  solver->setHintParam(OsiDoDualInResolve, true, OsiHintDo) ;
4617  // Sort
4618  CoinSort_2(sort, sort + numberToDo, whichObject);
4619  double * currentSolution = model->currentSolution();
4620  double objMin = 1.0e50;
4621  double objMax = -1.0e50;
4622  bool needResolve = false;
4623  int maxChoices=1;
4624  int currentChoice=0;
4625  int numberThreads=0;
4626  bool doAtEnd=false;
4627  if (model->parallelMode() && (solveType&4)!=0) {
4628    numberThreads=model->getNumberThreads();
4629    sprintf(general,"Using %d threads in analysis\n",numberThreads);
4630    model->messageHandler()->message(CBC_GENERAL, *model->messagesPointer())
4631      << general
4632      << CoinMessageEol;
4633    if (model->parallelMode()==1) {
4634      maxChoices=numberThreads;
4635    } else {
4636      maxChoices = numberToDo;
4637      if ((solveType&2)!=0) 
4638        maxChoices = numberColumns;
4639      doAtEnd=true;
4640    }
4641  }
4642  StrongInfo * choices = new StrongInfo[maxChoices];
4643  StrongStaticInfo staticInfo;
4644  int numberBundles = CoinMax(1,numberThreads);
4645  StrongBundle * bundles = new StrongBundle[numberBundles];
4646  /*
4647    0 - available - no need to look at results
4648    1 - not available
4649    2 - available - need to look at results
4650  */
4651#ifndef NUMBER_THREADS
4652#define NUMBER_THREADS 4
4653#endif
4654  int status[NUMBER_THREADS];
4655  memset(status,0,sizeof(status));
4656  memset(&staticInfo,0,sizeof(staticInfo));
4657  staticInfo.solveType = solveType;
4658  staticInfo.originalSolution=saveSolution;
4659  staticInfo.back=back;
4660  staticInfo.ws=ws;
4661  staticInfo.integerTolerance=integerTolerance;
4662  double time1 = model->getCurrentSeconds();
4663#define DO_STEEPEST_SERIAL 1
4664# ifdef COIN_HAS_CLP
4665  if (osiclp&&(solveType&2)!=0&&(!numberThreads||DO_STEEPEST_SERIAL)) {
4666    ClpSimplex * simplex = osiclp->getModelPtr();
4667    simplex->setLogLevel(0);
4668    simplex->dual(0,1);
4669    ClpDualRowPivot * dualRowPivot=simplex->dualRowPivot();
4670    ClpDualRowSteepest * steep = dynamic_cast<ClpDualRowSteepest *>(dualRowPivot);
4671    if (steep) {
4672      staticInfo.dualRowPivot=new ClpDualRowSteepest (*steep);
4673      staticInfo.dualRowPivot->setMode(1); // full steepest edge
4674      simplex->spareIntArray_[0]=0;
4675      simplex->spareIntArray_[1]=numberRows;
4676      staticInfo.dualRowPivot->saveWeights(simplex,7);
4677    }
4678  }
4679#endif
4680  for (int i=0;i<numberBundles;i++) {
4681    memset(bundles+i,0,sizeof(StrongBundle));
4682    bundles[i].staticInfo=&staticInfo;
4683  }
4684#if defined (THREADS_IN_ANALYZE) && defined (COIN_HAS_CLP)
4685#define USE_STRONG_THREADS
4686  CoinPthreadStuff threadInfo(numberThreads,cbc_parallelManager);
4687  int threadNeedsRefreshing[NUMBER_THREADS];
4688  for (int i=0;i<numberThreads;i++) {
4689    threadInfo.threadInfo_[i].extraInfo2 = solver->clone();
4690    threadNeedsRefreshing[i]=0;
4691  }
4692# ifdef COIN_HAS_CLP
4693  int numberSteepThreads=0;
4694  int step=numberThreads ? (numberRows+numberThreads-1)/numberThreads : 0;
4695  int first=0;
4696  for (int i=0;i<numberThreads;i++) {
4697    if (osiclp&&(solveType&2)!=0&&!DO_STEEPEST_SERIAL) {
4698      OsiSolverInterface * solver=
4699        reinterpret_cast<OsiSolverInterface *>(threadInfo.threadInfo_[i].extraInfo2);
4700      OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
4701      ClpSimplex * simplex = osiclp->getModelPtr();
4702      simplex->setLogLevel(0);
4703      simplex->dual(0,1);
4704      ClpDualRowPivot * dualRowPivot=simplex->dualRowPivot();
4705      ClpDualRowSteepest * steep = dynamic_cast<ClpDualRowSteepest *>(dualRowPivot);
4706      if (steep) {
4707        numberSteepThreads=numberThreads;
4708        ClpDualRowSteepest * dualRowPivot=new ClpDualRowSteepest (*steep);
4709        dualRowPivot->setMode(1); // full steepest edge
4710        simplex->spareIntArray_[0]=0;
4711        simplex->spareIntArray_[1]=numberRows;
4712        simplex->spareIntArray_[0]=first;
4713        simplex->spareIntArray_[1]=CoinMin(first+step,numberRows);
4714        first += step;
4715        if (i==0)
4716          staticInfo.dualRowPivot=dualRowPivot;
4717        choices[i].steepest=dualRowPivot;
4718        dualRowPivot->saveWeights(simplex,7);
4719      }
4720    }
4721  }
4722  if (numberSteepThreads&&false) {
4723    int numberDone=0;
4724    int iDo=0;
4725    staticInfo.solveType = 200;
4726    while (numberDone<numberSteepThreads) {
4727      int iThread;
4728      threadInfo.waitParallelTask(1,iThread,iDo<numberToDo);
4729      int threadStatus = 1+threadInfo.threadInfo_[iThread].status;
4730      iThread=iThread;
4731      if (threadStatus==0&&iDo<numberSteepThreads) {
4732        StrongInfo & choice = choices[iThread];
4733        StrongBundle & bundle = bundles[iThread];
4734        bundle.whichChoice=iThread;
4735        memset(&choice,0,sizeof(StrongInfo));
4736        iDo++; //started this one
4737        bundle.choice=&choice;
4738        bundle.solver = solver;
4739        bundle.solver=reinterpret_cast<OsiSolverInterface *>(threadInfo.threadInfo_[iThread].extraInfo2);
4740        threadStatus=0;
4741#ifdef DETAIL_THREAD
4742        printf("Starting steep task on thread %d\n",
4743               choice.iThread);
4744#endif
4745        threadInfo.startParallelTask(1,iThread,&bundle);
4746      }
4747      if (!threadStatus) {
4748#ifdef _MSC_VER
4749        Sleep(1);
4750#else
4751        usleep(1000);
4752#endif
4753        continue;
4754      }
4755      if (threadStatus) {
4756        numberDone++;
4757        // say available
4758        threadInfo.sayIdle(iThread);
4759      }
4760      staticInfo.solveType = solveType;
4761    }
4762    OsiSolverInterface * solver0=
4763      reinterpret_cast<OsiSolverInterface *>(threadInfo.threadInfo_[0].extraInfo2);
4764    CoinIndexedVector * savedWeights0 = staticInfo.dualRowPivot->savedWeights();
4765    int * index0 = savedWeights0->getIndices();
4766    double * weight0 = savedWeights0->denseVector();
4767    int step=(numberRows+numberSteepThreads-1)/numberSteepThreads;
4768    int first=step;
4769    //memset(weight0+first,0,(numberRows-first)*sizeof(double));
4770    for (int i=1;i<numberSteepThreads;i++) {
4771      int n=CoinMin(step,numberRows-first);
4772      CoinIndexedVector * savedWeights = choices[i].steepest->savedWeights();
4773      int * index = savedWeights->getIndices();
4774      double * weight = savedWeights->denseVector();
4775      memcpy(index0+first,index+first,n*sizeof(int));
4776      memcpy(weight0+first,weight+first,n*sizeof(double));
4777      first += step;
4778      delete choices[i].steepest;
4779      choices[i].steepest=NULL;
4780    }
4781    //for (int j=0;j<numberRows;j++)
4782    //weight0[j]=1.0;
4783  }
4784#endif
4785#endif
4786  double bestSolutionValue = model->getMinimizationObjValue();
4787  double * bestSolution = NULL;
4788  double cutoff;
4789  solver->getDblParam(OsiDualObjectiveLimit,cutoff);
4790  double maxMovement = 2.0*(cutoff-objectiveValue_)+1.0e-6;
4791  /*
4792    Now calculate the cost forcing the variable up and down.
4793  */
4794  int iDo=0;
4795  int iDone=-1;
4796  int numberDone=0;
4797  int iThread=0;
4798  int threadStatus=0;
4799  int whenPrint = (numberToDo+9)/10;
4800  while (numberDone<numberToDo) {
4801    if ((solveType&32)!=0&&numberDone==whenPrint) {
4802      whenPrint += (numberToDo+9)/10;
4803      sprintf(general,"%d variables looked at - %d changed by initial strong branching (%.2f seconds - %d iterations)",numberDone,numberToFix,model->getCurrentSeconds()-time1,saveAllowed-numberIterationsAllowed);
4804      model->messageHandler()->message(CBC_GENERAL, *model->messagesPointer())
4805        << general
4806        << CoinMessageEol;
4807    }
4808#ifdef USE_STRONG_THREADS
4809    if (numberThreads) {
4810      threadInfo.waitParallelTask(1,iThread,iDo<numberToDo);
4811      threadStatus = 1+threadInfo.threadInfo_[iThread].status;
4812      if (!doAtEnd)
4813        currentChoice=iThread;
4814      if (threadNeedsRefreshing[iThread]) { 
4815        OsiSolverInterface * solver=
4816          reinterpret_cast<OsiSolverInterface *>(threadInfo.threadInfo_[iThread].extraInfo2);
4817        if ((threadNeedsRefreshing[iThread]&1)!=0) 
4818          solver->setColLower(saveLower);
4819        if ((threadNeedsRefreshing[iThread]&2)!=0) 
4820          solver->setColUpper(saveUpper);
4821        threadNeedsRefreshing[iThread]=0;
4822      }
4823    }
4824#endif
4825    if (threadStatus==0&&iDo<numberToDo) {
4826      StrongInfo & choice = choices[currentChoice];
4827      StrongBundle & bundle = bundles[iThread];
4828      bundle.whichChoice=currentChoice;
4829      memset(&choice,0,sizeof(StrongInfo));
4830      int iObject = whichObject[iDo];
4831      iDo++; //started this one
4832      OsiObject * object = model->modifiableObject(iObject);
4833      CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
4834        dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
4835      int iColumn = dynamicObject->columnNumber();
4836      double value = currentSolution[iColumn];
4837      double nearest = floor(value + 0.5);
4838      double lowerValue = floor(value);
4839      bool satisfied = false;
4840      if (fabs(value - nearest) <= integerTolerance ||
4841          value < saveLower[iColumn] || value > saveUpper[iColumn]) {
4842        satisfied = true;
4843        if (nearest < saveUpper[iColumn]) {
4844          lowerValue = nearest;
4845        } else {
4846          lowerValue = nearest - 1;
4847        }
4848      }
4849      double upperValue = lowerValue + 1.0;
4850      // Save which object it was
4851      choice.columnNumber = iColumn;
4852      choice.initialValue=value;
4853      choice.upLowerBound=upperValue;
4854      choice.downUpperBound=lowerValue;
4855      choice.numIntInfeas[1] = numberUnsatisfied_;
4856      choice.numIntInfeas[0] = numberUnsatisfied_;
4857      choice.movement[0] = 0.0;
4858      choice.movement[1] = 0.0;
4859      choice.numIters[0] = 0;
4860      choice.numIters[1] = 0;
4861      if (fabs(value - lowerValue) <= integerTolerance) 
4862        choice.numIters[0]=-1; // mark as not done
4863      if (fabs(value - upperValue) <= integerTolerance) 
4864        choice.numIters[1]=-1; // mark as not done
4865      bundle.choice=&choice;
4866      bundle.solver = solver;
4867#ifdef USE_STRONG_THREADS
4868      if (numberThreads) {
4869        bundle.solver=reinterpret_cast<OsiSolverInterface *>(threadInfo.threadInfo_[iThread].extraInfo2);
4870        threadStatus=0;
4871#ifdef DETAIL_THREAD
4872        printf("Starting task for column %d on thread %d\n",
4873               choice.columnNumber,iThread);
4874#endif
4875        threadInfo.startParallelTask(1,iThread,&bundle);
4876      } else {
4877#endif
4878        threadStatus=2;
4879        solveAnalyze(&bundle);
4880#ifdef USE_STRONG_THREADS
4881      }
4882#endif
4883    }
4884    if (!threadStatus) {
4885#ifdef _MSC_VER
4886        Sleep(1);
4887#else
4888        usleep(1000);
4889#endif
4890      continue;
4891    }
4892    if (threadStatus) {
4893      int whichChoice = bundles[iThread].whichChoice; 
4894      StrongInfo & choice = choices[whichChoice];
4895      int iColumn=choice.columnNumber;
4896      if (choice.integerSolution) {
4897        double * foundSolution = choice.integerSolution;
4898        solver->setColSolution(foundSolution);
4899        // See if integer solution
4900        int numberInfeas=0;
4901        int numberOddInfeas=0;
4902        if (model->feasibleSolution(numberInfeas,numberOddInfeas)
4903            && model->problemFeasibility()->feasible(model, -1) >= 0) {
4904          double newObjectiveValue;
4905          solver->getDblParam(OsiObjOffset,newObjectiveValue);
4906          newObjectiveValue=-newObjectiveValue;
4907          const double * cost = solver->getObjCoefficients();
4908          for ( int i = 0 ; i < numberColumns ; i++ )
4909            newObjectiveValue += cost[i] * foundSolution[i];
4910          if (newObjectiveValue<bestSolutionValue) {
4911            if (doAtEnd) {
4912              if (!bestSolution)
4913                bestSolution = CoinCopyOfArray(foundSolution,numberColumns);
4914              else
4915                memcpy(bestSolution,foundSolution,numberColumns*sizeof(double));
4916              bestSolutionValue=newObjectiveValue;
4917            } else {
4918              model->setBestSolution(CBC_STRONGSOL,
4919                                     newObjectiveValue,
4920                                     foundSolution) ;
4921              model->setLastHeuristic(NULL);
4922              model->incrementUsed(solver->getColSolution());
4923              bestSolutionValue = model->getMinimizationObjValue();
4924            }
4925          }
4926        }
4927        delete [] foundSolution;
4928      }
4929      for (int iWay=0;iWay<2;iWay++) {
4930        numberIterationsAllowed -= choice.numIters[iWay];
4931        choice.movement[iWay] -= objectiveValue_;
4932      }
4933      // If objective goes above certain amount we can set bound
4934      int jInt = back[iColumn];
4935      OsiObject * object = model->modifiableObject(jInt);
4936      CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
4937        dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
4938      if (dynamicObject) {
4939        if (choice.numIters[0]>=0) {
4940          dynamicObject->addToSumDownCost(CoinMin(choice.movement[0],maxMovement));
4941          dynamicObject->addToSumDownChange(choice.initialValue-choice.downUpperBound);
4942        }
4943        if (choice.numIters[1]>=0) {
4944          dynamicObject->addToSumUpCost(CoinMin(choice.movement[1],maxMovement));
4945          dynamicObject->addToSumUpChange(choice.upLowerBound-choice.initialValue);
4946        }
4947      }
4948      newLower[jInt] = choice.upLowerBound;
4949      if (choice.finished[0])
4950        objLower[jInt] = choice.movement[0] + objectiveValue_;
4951      else
4952        objLower[jInt] = objectiveValue_;
4953      newUpper[jInt] = choice.downUpperBound;
4954      if (choice.finished[1])
4955        objUpper[jInt] = choice.movement[1] + objectiveValue_;
4956      else
4957        objUpper[jInt] = objectiveValue_;
4958      objMin = CoinMin(CoinMin(objLower[jInt], objUpper[jInt]), objMin);
4959      objMovement[2*jInt]=choice.movement[0];
4960      objMovement[2*jInt+1]=choice.movement[1];
4961      double sumModified = choice.modified[0] + choice.modified[1] +
4962        1.0e-15*(choice.sumModified[0]+choice.sumModified[1]);
4963      if (choice.numIters[0]>=0&&choice.numIters[1]>=0)
4964        sumModified *= 0.6;
4965      interAction[jInt] = sumModified;
4966      /*
4967        End of evaluation for this candidate variable. Possibilities are:
4968        * Both sides below cutoff; this variable is a candidate for branching.
4969        * Both sides infeasible or above the objective cutoff: no further action
4970        here. Break from the evaluation loop and assume the node will be purged
4971        by the caller.
4972        * One side below cutoff: Install the branch (i.e., fix the variable). Break
4973        from the evaluation loop and assume the node will be reoptimised by the
4974        caller.
4975      */
4976      threadStatus=0;
4977      currentChoice++;
4978      numberDone++;
4979#ifdef USE_STRONG_THREADS
4980      // say available
4981      if (numberThreads) {
4982        threadInfo.sayIdle(iThread);
4983      }
4984#endif
4985      if (doAtEnd)
4986        continue;
4987      if (choice.movement[1] < 1.0e100) {
4988        if (choice.movement[0] < 1.0e100) {
4989          objMax = CoinMax(CoinMax(objLower[jInt], objUpper[jInt]), objMax);
4990          // In case solution coming in was odd
4991          choice.movement[1] = CoinMax(0.0, choice.movement[1]);
4992          choice.movement[0] = CoinMax(0.0, choice.movement[0]);
4993          // feasible -
4994          model->messageHandler()->message(CBC_STRONG, *model->messagesPointer())
4995            << iColumn << iColumn
4996            << choice.movement[0] << choice.numIntInfeas[0]
4997            << choice.movement[1] << choice.numIntInfeas[1]
4998            << choice.initialValue
4999            << CoinMessageEol;
5000        } else {
5001          // up feasible, down infeasible
5002          needResolve = true;
5003          numberToFix++;
5004          saveLower[iColumn] = choice.upLowerBound;
5005          solver->setColLower(iColumn, choice.upLowerBound);
5006#ifdef USE_STRONG_THREADS
5007          for (int i=0;i<numberThreads;i++) {
5008            threadNeedsRefreshing[i] |= 1;
5009          }
5010#endif
5011        }
5012      } else {
5013        if (choice.movement[0] < 1.0e100) {
5014          // down feasible, up infeasible
5015          needResolve = true;
5016          numberToFix++;
5017          saveUpper[iColumn] = choice.downUpperBound;
5018          solver->setColUpper(iColumn, choice.downUpperBound);
5019#ifdef USE_STRONG_THREADS
5020          for (int i=0;i<numberThreads;i++) {
5021            threadNeedsRefreshing[i] |= 2;
5022          }
5023#endif
5024        } else {
5025          // neither side feasible
5026          COIN_DETAIL_PRINT(printf("Both infeasible for choice %d sequence %d\n", i,
5027                                   model->object(choice.objectNumber)->columnNumber()));
5028          //solver->writeMps("bad");
5029          numberToFix = -1;
5030          break;
5031        }
5032      }
5033      if (numberIterationsAllowed <= 0)
5034        break;
5035      if (currentChoice==maxChoices)
5036        currentChoice=0;
5037    }
5038    //printf("obj %d, col %d, down %g up %g value %g\n",iObject,iColumn,
5039    //     choice.downMovement,choice.upMovement,value);
5040  }
5041  // Do at end if deterministic
5042  if (doAtEnd) {
5043    if (bestSolution) {
5044      model->setBestSolution(CBC_STRONGSOL,
5045                             bestSolutionValue,
5046                             bestSolution) ;
5047      model->setLastHeuristic(NULL);
5048      model->incrementUsed(solver->getColSolution());
5049      delete [] bestSolution;
5050    }
5051    for (int iDo = 0; iDo < numberLookIntegers; iDo++) {
5052      StrongInfo & choice = choices[iDo];
5053      int iColumn = choice.columnNumber;
5054      int iObject = iColumn;
5055      int jInt = back[iColumn];
5056      double value = choice.initialValue;
5057      double lowerValue = choice.downUpperBound;
5058      double upperValue = choice.upLowerBound;
5059      if (choice.movement[1] < 1.0e100) {
5060        if (choice.movement[0] < 1.0e100) {
5061          objMax = CoinMax(CoinMax(objLower[jInt], objUpper[jInt]), objMax);
5062          // In case solution coming in was odd
5063          choice.movement[1] = CoinMax(0.0, choice.movement[1]);
5064          choice.movement[0] = CoinMax(0.0, choice.movement[0]);
5065          // feasible -
5066          model->messageHandler()->message(CBC_STRONG, *model->messagesPointer())
5067            << iObject << iColumn
5068            << choice.movement[0] << choice.numIntInfeas[0]
5069            << choice.movement[1] << choice.numIntInfeas[1]
5070            << value
5071            << CoinMessageEol;
5072        } else {
5073          // up feasible, down infeasible
5074          numberToFix++;
5075          saveLower[iColumn] = upperValue;
5076          solver->setColLower(iColumn, upperValue);
5077        }
5078      } else {
5079        if (choice.movement[0] < 1.0e100) {
5080          // down feasible, up infeasible
5081          needResolve = true;
5082          numberToFix++;
5083          saveUpper[iColumn] = lowerValue;
5084          solver->setColUpper(iColumn, lowerValue);
5085        } else {
5086          // neither side feasible
5087          COIN_DETAIL_PRINT(printf("Both infeasible for choice %d sequence %d\n", i,
5088                                   model->object(choice.objectNumber)->columnNumber()));
5089          //solver->writeMps("bad");
5090          numberToFix = -1;
5091          break;
5092        }
5093      }
5094    }
5095  }
5096  if (false) {
5097    const double * lower = solver->getColLower();
5098    const double * upper = solver->getColUpper();
5099    for (int i=0;i<numberColumns;i++) {
5100      if (lower[i]!=saveLower[i]||upper[i]!=saveUpper[i]) {
5101        printf("%d changed- saved %g,%g now %g,%g\n",i,saveLower[i],saveUpper[i],
5102               lower[i],upper[i]);
5103      }
5104    }
5105  }
5106  COIN_DETAIL_PRINT(printf("Best possible solution %g, can fix more if solution of %g found - looked at %d variables in %d iterations\n",
5107                           objMin, objMax, iDo, model->numberAnalyzeIterations() - numberIterationsAllowed));
5108  model->setNumberAnalyzeIterations(numberIterationsAllowed);
5109  if (numberToFix>0) {
5110    sprintf(general,"%d variable bounds modified by initial strong branching (%.2f seconds - %d iterations)",numberToFix,model->getCurrentSeconds()-time1,saveAllowed-numberIterationsAllowed);
5111  } else if (numberToFix<0) {
5112    sprintf(general,"initial strong branching found to be infeasible (%.2f seconds - %d iterations)",model->getCurrentSeconds()-time1,saveAllowed-numberIterationsAllowed);
5113  } else if ((solveType&32)!=0) {
5114    sprintf(general,"No variables fixed by initial strong branching (%.2f seconds - %d iterations)",model->getCurrentSeconds()-time1,saveAllowed-numberIterationsAllowed);
5115  } else {
5116    general[0]='\0';
5117  }
5118  if (general[0]!='\0')
5119    model->messageHandler()->message(CBC_GENERAL, *model->messagesPointer())
5120      << general
5121      << CoinMessageEol;
5122  double smallestEffect=COIN_DBL_MAX;
5123  double largestEffect=0.0;
5124  for (i = 0; i < numberIntegers; i++) {
5125    int iColumn=integerVariable[i];
5126    if (back[iColumn]>=numberColumns)
5127      continue;
5128    smallestEffect = CoinMin(smallestEffect,interAction[i]);
5129    largestEffect = CoinMax(largestEffect,interAction[i]);
5130  }
5131  double groupValue[11];
5132  int groupCounts[11]={0,0,0,0,0,0,0,0,0,0,0};
5133  groupValue[10]=largestEffect;
5134  for (int i=0;i<10;i++)
5135    groupValue[i]=smallestEffect+i*0.1*(largestEffect-smallestEffect);
5136  sprintf(general,"Looked at %d integer variables - smallest interaction %g",
5137          numberLookIntegers,smallestEffect);
5138  model->messageHandler()->message((solveType&32)==0 ? CBC_FPUMP2 : CBC_FPUMP1, 
5139                                    *model->messagesPointer())
5140    << general
5141    << CoinMessageEol;
5142  for (int i = 0; i < numberIntegers; i++) {
5143    int iColumn=integerVariable[i];
5144    if (back[iColumn]>=numberColumns)
5145      continue;
5146    double value = interAction[i];
5147    int j;
5148    for (j=0;j<11;j++) {
5149      if (value<=groupValue[j]||j==10)
5150        break;
5151    }
5152    groupCounts[j]++;
5153  }
5154  general[0]='\0';
5155  for (int i=0;i<11;i++)
5156    sprintf(general+strlen(general),"%d <= %g ",groupCounts[i],groupValue[i]);
5157  model->messageHandler()->message((solveType&32)==0 ? CBC_FPUMP2 : CBC_FPUMP1, 
5158                                    *model->messagesPointer())
5159    << general
5160    << CoinMessageEol;
5161  smallestEffect=COIN_DBL_MAX;
5162  largestEffect=0.0;
516