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

Last change on this file since 2103 was 2103, checked in by tkr, 3 years ago

Fix for missing usleep with MSVC. See r2004 in Clp

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