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

Last change on this file since 1659 was 1659, checked in by forrest, 8 years ago

fix for bonmin

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