source: releases/2.7.1/Cbc/src/CbcNode.cpp @ 1995

Last change on this file since 1995 was 1675, checked in by stefan, 8 years ago

sync with trunk rev1674

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 197.5 KB
Line 
1/* $Id: CbcNode.cpp 1675 2011-06-19 17:23:14Z tkr $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#if defined(_MSC_VER)
7// Turn off compiler warning about long names
8#  pragma warning(disable:4786)
9#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 = (model->getCurrentSeconds() >
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 = (model->getCurrentSeconds() >
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 = (model->getCurrentSeconds() >
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                            COIN_DETAIL_PRINT(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                            COIN_DETAIL_PRINT(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                      COIN_DETAIL_PRINT(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                        if (!solver->isProvenOptimal()) {
2803                          skipAll=-2;
2804                          canSkip = true;
2805                        }
2806                        xMark++;
2807                    }
2808                }
2809                if (!canSkip) {
2810                    numberTest--;
2811                    // just do a few
2812                    if (searchStrategy == 2)
2813                        solver->setIntParam(OsiMaxNumIterationHotStart, 10);
2814                    double objectiveChange ;
2815                    double newObjectiveValue = 1.0e100;
2816                    int j;
2817                    // status is 0 finished, 1 infeasible and other
2818                    int iStatus;
2819                    /*
2820                      Try the down direction first. (Specify the initial branching alternative as
2821                      down with a call to way(-1). Each subsequent call to branch() performs the
2822                      specified branch and advances the branch object state to the next branch
2823                      alternative.)
2824                    */
2825                    choice.possibleBranch->way(-1) ;
2826                    choice.possibleBranch->branch() ;
2827                    solver->solveFromHotStart() ;
2828                    bool needHotStartUpdate = false;
2829                    numberStrongDone++;
2830                    numberStrongIterations += solver->getIterationCount();
2831                    /*
2832                      We now have an estimate of objective degradation that we can use for strong
2833                      branching. If we're over the cutoff, the variable is monotone up.
2834                      If we actually made it to optimality, check for a solution, and if we have
2835                      a good one, call setBestSolution to process it. Note that this may reduce the
2836                      cutoff, so we check again to see if we can declare this variable monotone.
2837                    */
2838                    if (solver->isProvenOptimal())
2839                        iStatus = 0; // optimal
2840                    else if (solver->isIterationLimitReached() 
2841                             && !solver->isDualObjectiveLimitReached()) {
2842                        iStatus = 2; // unknown
2843                    } else {
2844                        iStatus = 1; // infeasible
2845                    }
2846                    if (iStatus != 2 && solver->getIterationCount() >
2847                            realMaxHotIterations)
2848                        numberUnfinished++;
2849                    newObjectiveValue = solver->getObjSense() * solver->getObjValue();
2850                    choice.numItersDown = solver->getIterationCount();
2851                    objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_, 0.0);
2852                    // Update branching information if wanted
2853                    CbcBranchingObject * cbcobj = dynamic_cast<CbcBranchingObject *> (choice.possibleBranch);
2854                    if (cbcobj) {
2855                        CbcObject * object = cbcobj->object();
2856                        assert (object) ;
2857                        CbcObjectUpdateData update = object->createUpdateInformation(solver, this, cbcobj);
2858                        update.objectNumber_ = choice.objectNumber;
2859                        model->addUpdateInformation(update);
2860                    } else {
2861                        decision->updateInformation( solver, this);
2862                    }
2863                    if (!iStatus) {
2864                        choice.finishedDown = true ;
2865                        if (newObjectiveValue >= cutoff) {
2866                            objectiveChange = 1.0e100; // say infeasible
2867                            numberStrongInfeasible++;
2868                        } else {
2869#define CBCNODE_TIGHTEN_BOUNDS
2870#ifdef CBCNODE_TIGHTEN_BOUNDS
2871                            // Can we tighten bounds?
2872                            if (iColumn < numberColumns && cutoff < 1.0e20
2873                                    && objectiveChange > 1.0e-5) {
2874                                double value = saveSolution[iColumn];
2875                                double down = value - floor(value);
2876                                double changePer = objectiveChange / (down + 1.0e-7);
2877                                double distance = (cutoff - objectiveValue_) / changePer;
2878                                distance += 1.0e-3;
2879                                if (distance < 5.0) {
2880                                    double newLower = ceil(value - distance);
2881                                    if (newLower > saveLower[iColumn]) {
2882                                        //printf("Could increase lower bound on %d from %g to %g\n",
2883                                        //   iColumn,saveLower[iColumn],newLower);
2884                                        saveLower[iColumn] = newLower;
2885                                        solver->setColLower(iColumn, newLower);
2886                                    }
2887                                }
2888                            }
2889#endif
2890                            // See if integer solution
2891                            if (model->feasibleSolution(choice.numIntInfeasDown,
2892                                                        choice.numObjInfeasDown)
2893                                    && model->problemFeasibility()->feasible(model, -1) >= 0) {
2894                                if (auxiliaryInfo->solutionAddsCuts()) {
2895                                    needHotStartUpdate = true;
2896                                    solver->unmarkHotStart();
2897                                }
2898                                model->setLogLevel(saveLogLevel);
2899                                model->setBestSolution(CBC_STRONGSOL,
2900                                                       newObjectiveValue,
2901                                                       solver->getColSolution()) ;
2902                                if (needHotStartUpdate) {
2903                                    model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
2904                                    newObjectiveValue = solver->getObjSense() * solver->getObjValue();
2905                                    objectiveValue_ = CoinMax(objectiveValue_,newObjectiveValue);
2906                                    objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_, 0.0);
2907                                    model->feasibleSolution(choice.numIntInfeasDown,
2908                                                            choice.numObjInfeasDown);
2909                                }
2910                                model->setLastHeuristic(NULL);
2911                                model->incrementUsed(solver->getColSolution());
2912                                cutoff = model->getCutoff();
2913                                if (newObjectiveValue >= cutoff) {      //  *new* cutoff
2914                                    objectiveChange = 1.0e100 ;
2915                                    numberStrongInfeasible++;
2916                                }
2917                            }
2918                        }
2919                    } else if (iStatus == 1) {
2920                        objectiveChange = 1.0e100 ;
2921                        numberStrongInfeasible++;
2922                    } else {
2923                        // Can't say much as we did not finish
2924                        choice.finishedDown = false ;
2925                        numberUnfinished++;
2926                    }
2927                    choice.downMovement = objectiveChange ;
2928
2929                    // restore bounds
2930                    for ( j = 0; j < numberColumns; j++) {
2931                        if (saveLower[j] != lower[j])
2932                            solver->setColLower(j, saveLower[j]);
2933                        if (saveUpper[j] != upper[j])
2934                            solver->setColUpper(j, saveUpper[j]);
2935                    }
2936                    if (needHotStartUpdate) {
2937                        needHotStartUpdate = false;
2938                        model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
2939                        double newObjValue = solver->getObjSense()*solver->getObjValue();
2940                        objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
2941                        //we may again have an integer feasible solution
2942                        int numberIntegerInfeasibilities;
2943                        int numberObjectInfeasibilities;
2944                        if (model->feasibleSolution(
2945                                    numberIntegerInfeasibilities,
2946                                    numberObjectInfeasibilities)) {
2947#ifdef BONMIN
2948                            //In this case node has become integer feasible, let us exit the loop
2949                            std::cout << "Node has become integer feasible" << std::endl;
2950                            numberUnsatisfied_ = 0;
2951                            break;
2952#endif
2953                            double objValue = solver->getObjValue();
2954                            model->setLogLevel(saveLogLevel);
2955                            model->setBestSolution(CBC_STRONGSOL,
2956                                                   objValue,
2957                                                   solver->getColSolution()) ;
2958                            model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
2959                            double newObjValue = solver->getObjSense()*solver->getObjValue();
2960                            objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
2961                            cutoff = model->getCutoff();
2962                        }
2963                        solver->markHotStart();
2964                    }
2965#ifdef DO_ALL_AT_ROOT
2966                    if (!numberNodes)
2967                        printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
2968                               choice.objectNumber, iStatus, newObjectiveValue, choice.numItersDown,
2969                               choice.downMovement, choice.finishedDown, choice.numIntInfeasDown,
2970                               choice.numObjInfeasDown);
2971#endif
2972
2973                    // repeat the whole exercise, forcing the variable up
2974                    choice.possibleBranch->branch();
2975                    solver->solveFromHotStart() ;
2976                    numberStrongDone++;
2977                    numberStrongIterations += solver->getIterationCount();
2978                    /*
2979                      We now have an estimate of objective degradation that we can use for strong
2980                      branching. If we're over the cutoff, the variable is monotone up.
2981                      If we actually made it to optimality, check for a solution, and if we have
2982                      a good one, call setBestSolution to process it. Note that this may reduce the
2983                      cutoff, so we check again to see if we can declare this variable monotone.
2984                    */
2985                    if (solver->isProvenOptimal())
2986                        iStatus = 0; // optimal
2987                    else if (solver->isIterationLimitReached()
2988                             && !solver->isDualObjectiveLimitReached()) {
2989                        iStatus = 2; // unknown
2990                    } else {
2991                        iStatus = 1; // infeasible
2992                    }
2993                    if (iStatus != 2 && solver->getIterationCount() >
2994                            realMaxHotIterations)
2995                        numberUnfinished++;
2996                    newObjectiveValue = solver->getObjSense() * solver->getObjValue();
2997                    choice.numItersUp = solver->getIterationCount();
2998                    objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_, 0.0);
2999                    // Update branching information if wanted
3000                    cbcobj = dynamic_cast<CbcBranchingObject *> (choice.possibleBranch);
3001                    if (cbcobj) {
3002                        CbcObject * object = cbcobj->object();
3003                        assert (object) ;
3004                        CbcObjectUpdateData update = object->createUpdateInformation(solver, this, cbcobj);
3005                        update.objectNumber_ = choice.objectNumber;
3006                        model->addUpdateInformation(update);
3007                    } else {
3008                        decision->updateInformation( solver, this);
3009                    }
3010                    if (!iStatus) {
3011                        choice.finishedUp = true ;
3012                        if (newObjectiveValue >= cutoff) {
3013                            objectiveChange = 1.0e100; // say infeasible
3014                            numberStrongInfeasible++;
3015                        } else {
3016#ifdef CBCNODE_TIGHTEN_BOUNDS
3017                            // Can we tighten bounds?
3018                            if (iColumn < numberColumns && cutoff < 1.0e20
3019                                    && objectiveChange > 1.0e-5) {
3020                                double value = saveSolution[iColumn];
3021                                double up = ceil(value) - value;
3022                                double changePer = objectiveChange / (up + 1.0e-7);
3023                                double distance = (cutoff - objectiveValue_) / changePer;
3024                                distance += 1.0e-3;
3025                                if (distance < 5.0) {
3026                                    double newUpper = floor(value + distance);
3027                                    if (newUpper < saveUpper[iColumn]) {
3028                                        //printf("Could decrease upper bound on %d from %g to %g\n",
3029                                        //   iColumn,saveUpper[iColumn],newUpper);
3030                                        saveUpper[iColumn] = newUpper;
3031                                        solver->setColUpper(iColumn, newUpper);
3032                                    }
3033                                }
3034                            }
3035#endif
3036                            // See if integer solution
3037                            if (model->feasibleSolution(choice.numIntInfeasUp,
3038                                                        choice.numObjInfeasUp)
3039                                    && model->problemFeasibility()->feasible(model, -1) >= 0) {
3040#ifdef BONMIN
3041                                std::cout << "Node has become integer feasible" << std::endl;
3042                                numberUnsatisfied_ = 0;
3043                                break;
3044#endif
3045                                if (auxiliaryInfo->solutionAddsCuts()) {
3046                                    needHotStartUpdate = true;
3047                                    solver->unmarkHotStart();
3048                                }
3049                                model->setLogLevel(saveLogLevel);
3050                                model->setBestSolution(CBC_STRONGSOL,
3051                                                       newObjectiveValue,
3052                                                       solver->getColSolution()) ;
3053                                if (choice.finishedDown) {
3054                                  double cutoff = model->getCutoff();
3055                                  double downObj = objectiveValue_
3056                                    + choice.downMovement ;
3057                                  if (downObj >= cutoff) {     
3058                                    choice.downMovement = 1.0e100 ;
3059                                    numberStrongInfeasible++;
3060                                }
3061
3062                                }
3063                                if (needHotStartUpdate) {
3064                                    model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3065                                    newObjectiveValue = solver->getObjSense() * solver->getObjValue();
3066                                    objectiveValue_ = CoinMax(objectiveValue_,newObjectiveValue);
3067                                    objectiveChange = CoinMax(newObjectiveValue  - objectiveValue_, 0.0);
3068                                    model->feasibleSolution(choice.numIntInfeasDown,
3069                                                            choice.numObjInfeasDown);
3070                                }
3071                                model->setLastHeuristic(NULL);
3072                                model->incrementUsed(solver->getColSolution());
3073                                cutoff = model->getCutoff();
3074                                if (newObjectiveValue >= cutoff) {      //  *new* cutoff
3075                                    objectiveChange = 1.0e100 ;
3076                                    numberStrongInfeasible++;
3077                                }
3078                            }
3079                        }
3080                    } else if (iStatus == 1) {
3081                        objectiveChange = 1.0e100 ;
3082                        numberStrongInfeasible++;
3083                    } else {
3084                        // Can't say much as we did not finish
3085                        choice.finishedUp = false ;
3086                        numberUnfinished++;
3087                    }
3088                    choice.upMovement = objectiveChange ;
3089
3090                    // restore bounds
3091                    for ( j = 0; j < numberColumns; j++) {
3092                        if (saveLower[j] != lower[j])
3093                            solver->setColLower(j, saveLower[j]);
3094                        if (saveUpper[j] != upper[j])
3095                            solver->setColUpper(j, saveUpper[j]);
3096                    }
3097                    if (needHotStartUpdate) {
3098                        needHotStartUpdate = false;
3099                        model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3100                        double newObjValue = solver->getObjSense()*solver->getObjValue();
3101                        objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3102                        //we may again have an integer feasible solution
3103                        int numberIntegerInfeasibilities;
3104                        int numberObjectInfeasibilities;
3105                        if (model->feasibleSolution(
3106                                    numberIntegerInfeasibilities,
3107                                    numberObjectInfeasibilities)) {
3108                            double objValue = solver->getObjValue();
3109                            model->setLogLevel(saveLogLevel);
3110                            model->setBestSolution(CBC_STRONGSOL,
3111                                                   objValue,
3112                                                   solver->getColSolution()) ;
3113                            model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3114                            double newObjValue = solver->getObjSense()*solver->getObjValue();
3115                            objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3116                            cutoff = model->getCutoff();
3117                        }
3118                        solver->markHotStart();
3119                    }
3120
3121#ifdef DO_ALL_AT_ROOT
3122                    if (!numberNodes)
3123                        printf("Up on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
3124                               choice.objectNumber, iStatus, newObjectiveValue, choice.numItersUp,
3125                               choice.upMovement, choice.finishedUp, choice.numIntInfeasUp,
3126                               choice.numObjInfeasUp);
3127#endif
3128                }
3129
3130                solver->setIntParam(OsiMaxNumIterationHotStart, saveLimit2);
3131                /*
3132                  End of evaluation for this candidate variable. Possibilities are:
3133                  * Both sides below cutoff; this variable is a candidate for branching.
3134                  * Both sides infeasible or above the objective cutoff: no further action
3135                  here. Break from the evaluation loop and assume the node will be purged
3136                  by the caller.
3137                  * One side below cutoff: Install the branch (i.e., fix the variable). Break
3138                  from the evaluation loop and assume the node will be reoptimised by the
3139                  caller.
3140                */
3141                // reset
3142                choice.possibleBranch->resetNumberBranchesLeft();
3143                if (choice.upMovement < 1.0e100) {
3144                    if (choice.downMovement < 1.0e100) {
3145                        // In case solution coming in was odd
3146                        choice.upMovement = CoinMax(0.0, choice.upMovement);
3147                        choice.downMovement = CoinMax(0.0, choice.downMovement);
3148#if ZERO_ONE==2
3149                        // branch on 0-1 first (temp)
3150                        if (fabs(choice.possibleBranch->value()) < 1.0) {
3151                            choice.upMovement *= ZERO_FAKE;
3152                            choice.downMovement *= ZERO_FAKE;
3153                        }
3154#endif
3155                        // feasible - see which best
3156                        if (!canSkip) {
3157                            if (model->messageHandler()->logLevel() > 3)
3158                                printf("sort %g downest %g upest %g ", sort[iDo], downEstimate[iObject],
3159                                       upEstimate[iObject]);
3160                            model->messageHandler()->message(CBC_STRONG, *model->messagesPointer())
3161                            << iObject << iColumn
3162                            << choice.downMovement << choice.numIntInfeasDown
3163                            << choice.upMovement << choice.numIntInfeasUp
3164                            << choice.possibleBranch->value()
3165                            << CoinMessageEol;
3166                        }
3167                        int betterWay;
3168                        {
3169                            CbcBranchingObject * branchObj =
3170                                dynamic_cast <CbcBranchingObject *>(branch_) ;
3171                            if (branch_)
3172                                assert (branchObj);
3173                            betterWay = decision->betterBranch(choice.possibleBranch,
3174                                                               branchObj,
3175                                                               choice.upMovement,
3176                                                               choice.numIntInfeasUp ,
3177                                                               choice.downMovement,
3178                                                               choice.numIntInfeasDown );
3179                        }
3180                        if (betterWay) {
3181                            // C) create branching object
3182                            if (choiceObject) {
3183                                delete branch_;
3184                                branch_ = choice.possibleBranch->clone();
3185                            } else {
3186                                delete branch_;
3187                                branch_ = choice.possibleBranch;
3188                                choice.possibleBranch = NULL;
3189                            }
3190                            {
3191                                CbcBranchingObject * branchObj =
3192                                    dynamic_cast <CbcBranchingObject *>(branch_) ;
3193                                assert (branchObj);
3194                                //branchObj->way(preferredWay);
3195                                branchObj->way(betterWay);
3196                            }
3197                            bestChoice = choice.objectNumber;
3198                            whichChoice = iDo;
3199                            if (numberStrong <= 1) {
3200                                delete ws;
3201                                ws = NULL;
3202                                break;
3203                            }
3204                        } else {
3205                            if (!choiceObject) {
3206                                delete choice.possibleBranch;
3207                                choice.possibleBranch = NULL;
3208                            }
3209                            if (iDo >= 2*numberStrong) {
3210                                delete ws;
3211                                ws = NULL;
3212                                break;
3213                            }
3214                            if (!dynamicObject || dynamicObject->numberTimesUp() > 1) {
3215                                if (iDo - whichChoice >= numberStrong) {
3216                                    if (!choiceObject) {
3217                                        delete choice.possibleBranch;
3218                                        choice.possibleBranch = NULL;
3219                                    }
3220                                    break; // give up
3221                                }
3222                            } else {
3223                                if (iDo - whichChoice >= 2*numberStrong) {
3224                                    delete ws;
3225                                    ws = NULL;
3226                                    if (!choiceObject) {
3227                                        delete choice.possibleBranch;
3228                                        choice.possibleBranch = NULL;
3229                                    }
3230                                    break; // give up
3231                                }
3232                            }
3233                        }
3234                    } else {
3235                        // up feasible, down infeasible
3236                        anyAction = -1;
3237                        worstFeasible = CoinMax(worstFeasible, choice.upMovement);
3238                        model->messageHandler()->message(CBC_STRONG, *model->messagesPointer())
3239                        << iObject << iColumn
3240                        << choice.downMovement << choice.numIntInfeasDown
3241                        << choice.upMovement << choice.numIntInfeasUp
3242                        << choice.possibleBranch->value()
3243                        << CoinMessageEol;
3244                        //printf("Down infeasible for choice %d sequence %d\n",i,
3245                        // model->object(choice.objectNumber)->columnNumber());
3246                        choice.fix = 1;
3247                        numberToFix++;
3248                        choice.possibleBranch->fix(solver, saveLower, saveUpper, 1);
3249                        if (!choiceObject) {
3250                            delete choice.possibleBranch;
3251                            choice.possibleBranch = NULL;
3252                        } else {
3253                            //choiceObject = new CbcDynamicPseudoCostBranchingObject(*choiceObject);
3254                            choice.possibleBranch = choiceObject;
3255                        }
3256                        assert(doneHotStart);
3257                        solver->unmarkHotStart();
3258                        model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3259                        double newObjValue = solver->getObjSense()*solver->getObjValue();
3260                        objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3261                        solver->markHotStart();
3262                        // may be infeasible (if other way stopped on iterations)
3263                        if (!solver->isProvenOptimal()||solver->isDualObjectiveLimitReached()) {
3264                            // neither side feasible
3265                            anyAction = -2;
3266                            if (!choiceObject) {
3267                                delete choice.possibleBranch;
3268                                choice.possibleBranch = NULL;
3269                            }
3270                            //printf("Both infeasible for choice %d sequence %d\n",i,
3271                            // model->object(choice.objectNumber)->columnNumber());
3272                            delete ws;
3273                            ws = NULL;
3274                            break;
3275                        }
3276                    }
3277                } else {
3278                    if (choice.downMovement < 1.0e100) {
3279                        // down feasible, up infeasible
3280                        anyAction = -1;
3281                        worstFeasible = CoinMax(worstFeasible, choice.downMovement);
3282                        model->messageHandler()->message(CBC_STRONG, *model->messagesPointer())
3283                        << iObject << iColumn
3284                        << choice.downMovement << choice.numIntInfeasDown
3285                        << choice.upMovement << choice.numIntInfeasUp
3286                        << choice.possibleBranch->value()
3287                        << CoinMessageEol;
3288                        choice.fix = -1;
3289                        numberToFix++;
3290                        choice.possibleBranch->fix(solver, saveLower, saveUpper, -1);
3291                        if (!choiceObject) {
3292                            delete choice.possibleBranch;
3293                            choice.possibleBranch = NULL;
3294                        } else {
3295                            //choiceObject = new CbcDynamicPseudoCostBranchingObject(*choiceObject);
3296                            choice.possibleBranch = choiceObject;
3297                        }
3298                        assert(doneHotStart);
3299                        solver->unmarkHotStart();
3300                        model->resolve(NULL, 11, saveSolution, saveLower, saveUpper);
3301                        double newObjValue = solver->getObjSense()*solver->getObjValue();
3302                        objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3303                        solver->markHotStart();
3304                        // may be infeasible (if other way stopped on iterations)
3305                        if (!solver->isProvenOptimal()||solver->isDualObjectiveLimitReached()) {
3306                            // neither side feasible
3307                            anyAction = -2;
3308                            if (!choiceObject) {
3309                                delete choice.possibleBranch;
3310                                choice.possibleBranch = NULL;
3311                            }
3312                            delete ws;
3313                            ws = NULL;
3314                            break;
3315                        }
3316                    } else {
3317                        // neither side feasible
3318                        anyAction = -2;
3319                        if (!choiceObject) {
3320                            delete choice.possibleBranch;
3321                            choice.possibleBranch = NULL;
3322                        }
3323                        delete ws;
3324                        ws = NULL;
3325                        break;
3326                    }
3327                }
3328                // Check max time
3329                hitMaxTime = (model->getCurrentSeconds() >
3330                               model->getDblParam(CbcModel::CbcMaximumSeconds));
3331                if (hitMaxTime) {
3332                    // make sure rest are fast
3333                    doQuickly = true;
3334                    for ( int jDo = iDo + 1; jDo < numberToDo; jDo++) {
3335                        int iObject = whichObject[iDo];
3336                        OsiObject * object = model->modifiableObject(iObject);
3337                        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3338                            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3339                        if (dynamicObject)
3340                            dynamicObject->setNumberBeforeTrust(0);
3341                    }
3342                    numberTest = 0;
3343                }
3344                if (!choiceObject) {
3345                    delete choice.possibleBranch;
3346                }
3347            }
3348            if (model->messageHandler()->logLevel() > 3) {
3349                if (anyAction == -2) {
3350                    printf("infeasible\n");
3351                } else if (anyAction == -1) {
3352                    printf("%d fixed AND choosing %d iDo %d iChosenWhen %d numberToDo %d\n", numberToFix, bestChoice,
3353                           iDo, whichChoice, numberToDo);
3354                } else {
3355                    int iObject = whichObject[whichChoice];
3356                    OsiObject * object = model->modifiableObject(iObject);
3357                    CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3358                        dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3359                    if (dynamicObject) {
3360                        int iColumn = dynamicObject->columnNumber();
3361                        printf("choosing %d (column %d) iChosenWhen %d numberToDo %d\n", bestChoice,
3362                               iColumn, whichChoice, numberToDo);
3363                    }
3364                }
3365            }
3366            if (doneHotStart) {
3367                // Delete the snapshot
3368                solver->unmarkHotStart();
3369                // back to normal
3370                solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
3371                // restore basis
3372                solver->setWarmStart(ws);
3373            }
3374            solver->setIntParam(OsiMaxNumIterationHotStart, saveLimit);
3375            // Unless infeasible we will carry on
3376            // But we could fix anyway
3377            if (numberToFix && !hitMaxTime) {
3378                if (anyAction != -2) {
3379                    // apply and take off
3380                    bool feasible = true;
3381                    // can do quick optimality check
3382                    int easy = 2;
3383                    solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, &easy) ;
3384                    model->resolve(NULL, 11, saveSolution, saveLower, saveUpper) ;
3385                    double newObjValue = solver->getObjSense()*solver->getObjValue();
3386                    objectiveValue_ = CoinMax(objectiveValue_,newObjValue);
3387                    solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
3388                    feasible = solver->isProvenOptimal();
3389                    if (feasible) {
3390                        anyAction = 0;
3391                    } else {
3392                        anyAction = -2;
3393                        finished = true;
3394                    }
3395                }
3396            }
3397            // If  fixed then round again
3398            // See if candidate still possible
3399            if (branch_) {
3400                 const OsiObject * object = model->object(bestChoice);
3401                 int preferredWay;
3402                 double infeasibility = object->infeasibility(&usefulInfo, preferredWay);
3403                 if (!infeasibility) {
3404                   // take out
3405                   delete branch_;
3406                   branch_ = NULL;
3407                 } else {
3408                   CbcBranchingObject * branchObj =
3409                     dynamic_cast <CbcBranchingObject *>(branch_) ;
3410                   assert (branchObj);
3411                   branchObj->way(preferredWay);
3412#ifdef CBCNODE_TIGHTEN_BOUNDS
3413                   bool fixed = branchObj->tighten(solver);
3414                   if (fixed) {
3415                     //printf("Variable now fixed!\n");
3416                     // take out
3417                     delete branch_;
3418                     branch_ = NULL;
3419                   }
3420#endif
3421                 }
3422            }
3423            if (!branch_ && anyAction != -2 && !hitMaxTime) {
3424                finished = false;
3425            }
3426            delete ws;
3427        }
3428    }
3429    // update number of strong iterations etc
3430    model->incrementStrongInfo(numberStrongDone, numberStrongIterations,
3431                               anyAction == -2 ? 0 : numberToFix, anyAction == -2);
3432    if (model->searchStrategy() == -1) {
3433#ifndef COIN_DEVELOP
3434        if (solver->messageHandler()->logLevel() > 1)
3435#endif
3436            printf("%d strong, %d iters, %d inf, %d not finished, %d not trusted\n",
3437                   numberStrongDone, numberStrongIterations, numberStrongInfeasible, numberUnfinished,
3438                   numberNotTrusted);
3439        // decide what to do
3440        int strategy = 1;
3441        if (((numberUnfinished*4 > numberStrongDone &&
3442                numberStrongInfeasible*40 < numberStrongDone) ||
3443                numberStrongInfeasible < 0) && model->numberStrong() < 10 && model->numberBeforeTrust() <= 20 && model->numberObjects() > CoinMax(1000, solver->getNumRows())) {
3444            strategy = 2;
3445#ifdef COIN_DEVELOP
3446            //if (model->logLevel()>1)
3447            printf("going to strategy 2\n");
3448#endif
3449            // Weaken
3450            model->setNumberStrong(2);
3451            model->setNumberBeforeTrust(1);
3452            model->synchronizeNumberBeforeTrust();
3453        }
3454        if (numberNodes)
3455            strategy = 1;  // should only happen after hot start
3456        model->setSearchStrategy(strategy);
3457    } else if (numberStrongDone) {
3458        //printf("%d strongB, %d iters, %d inf, %d not finished, %d not trusted\n",
3459        //   numberStrongDone,numberStrongIterations,numberStrongInfeasible,numberUnfinished,
3460        //   numberNotTrusted);
3461    }
3462    if (model->searchStrategy() == 1 && numberNodes > 500 && numberNodes < -510) {
3463#ifndef COIN_DEVELOP
3464        if (solver->messageHandler()->logLevel() > 1)
3465#endif
3466            printf("after %d nodes - %d strong, %d iters, %d inf, %d not finished, %d not trusted\n",
3467                   numberNodes, numberStrongDone, numberStrongIterations, numberStrongInfeasible, numberUnfinished,
3468                   numberNotTrusted);
3469        // decide what to do
3470        if (numberUnfinished*10 > numberStrongDone + 1 ||
3471                !numberStrongInfeasible) {
3472          COIN_DETAIL_PRINT(printf("going to strategy 2\n"));
3473            // Weaken
3474            model->setNumberStrong(2);
3475            model->setNumberBeforeTrust(1);
3476            model->synchronizeNumberBeforeTrust();
3477            model->setSearchStrategy(2);
3478        }
3479    }
3480    if (numberUnfinished*10 < numberStrongDone &&
3481            numberStrongIterations*20 < model->getIterationCount()&&
3482                                !auxiliaryInfo->solutionAddsCuts()) {
3483        //printf("increasing trust\n");
3484        model->synchronizeNumberBeforeTrust(2);
3485    }
3486
3487    // Set guessed solution value
3488    guessedObjectiveValue_ = objectiveValue_ + estimatedDegradation;
3489
3490    /*
3491      Cleanup, then we're finished
3492    */
3493    if (!model->branchingMethod())
3494        delete decision;
3495
3496    delete choiceObject;
3497    delete [] sort;
3498    delete [] whichObject;
3499#ifdef RANGING
3500    delete [] objectMark;
3501#endif
3502    delete [] saveLower;
3503    delete [] saveUpper;
3504    delete [] upEstimate;
3505    delete [] downEstimate;
3506# ifdef COIN_HAS_CLP
3507    if (osiclp) {
3508        osiclp->setSpecialOptions(saveClpOptions);
3509    }
3510# endif
3511    // restore solution
3512    solver->setColSolution(saveSolution);
3513    model->reserveCurrentSolution(saveSolution);
3514    delete [] saveSolution;
3515    model->setStateOfSearch(saveStateOfSearch);
3516    model->setLogLevel(saveLogLevel);
3517    // delete extra regions
3518    if (usefulInfo.usefulRegion_) {
3519        delete [] usefulInfo.usefulRegion_;
3520        delete [] usefulInfo.indexRegion_;
3521        delete [] usefulInfo.pi_;
3522        usefulInfo.usefulRegion_ = NULL;
3523        usefulInfo.indexRegion_ = NULL;
3524        usefulInfo.pi_ = NULL;
3525    }
3526    useShadow = model->moreSpecialOptions() & 7;
3527    if ((useShadow == 5 && model->getSolutionCount()) || useShadow == 6) {
3528        // zap pseudo shadow prices
3529        model->pseudoShadow(-1);
3530        // and switch off
3531        model->setMoreSpecialOptions(model->moreSpecialOptions()&(~1023));
3532    }
3533    return anyAction;
3534}
3535int CbcNode::analyze (CbcModel *model, double * results)
3536{
3537    int i;
3538    int numberIterationsAllowed = model->numberAnalyzeIterations();
3539    OsiSolverInterface * solver = model->solver();
3540    objectiveValue_ = solver->getObjSense() * solver->getObjValue();
3541    double cutoff = model->getCutoff();
3542    const double * lower = solver->getColLower();
3543    const double * upper = solver->getColUpper();
3544    const double * dj = solver->getReducedCost();
3545    int numberObjects = model->numberObjects();
3546    int numberColumns = model->getNumCols();
3547    // Initialize arrays
3548    int numberIntegers = model->numberIntegers();
3549    int * back = new int[numberColumns];
3550    const int * integerVariable = model->integerVariable();
3551    for (i = 0; i < numberColumns; i++)
3552        back[i] = -1;
3553    // What results is
3554    double * newLower = results;
3555    double * objLower = newLower + numberIntegers;
3556    double * newUpper = objLower + numberIntegers;
3557    double * objUpper = newUpper + numberIntegers;
3558    for (i = 0; i < numberIntegers; i++) {
3559        int iColumn = integerVariable[i];
3560        back[iColumn] = i;
3561        newLower[i] = 0.0;
3562        objLower[i] = -COIN_DBL_MAX;
3563        newUpper[i] = 0.0;
3564        objUpper[i] = -COIN_DBL_MAX;
3565    }
3566    double * saveUpper = new double[numberColumns];
3567    double * saveLower = new double[numberColumns];
3568    int anyAction = 0;
3569    // Save solution in case heuristics need good solution later
3570
3571    double * saveSolution = new double[numberColumns];
3572    memcpy(saveSolution, solver->getColSolution(), numberColumns*sizeof(double));
3573    model->reserveCurrentSolution(saveSolution);
3574    for (i = 0; i < numberColumns; i++) {
3575        saveLower[i] = lower[i];
3576        saveUpper[i] = upper[i];
3577    }
3578    // Get arrays to sort
3579    double * sort = new double[numberObjects];
3580    int * whichObject = new int[numberObjects];
3581    int numberToFix = 0;
3582    int numberToDo = 0;
3583    double integerTolerance =
3584        model->getDblParam(CbcModel::CbcIntegerTolerance);
3585    // point to useful information
3586    OsiBranchingInformation usefulInfo = model->usefulInformation();
3587    // and modify
3588    usefulInfo.depth_ = depth_;
3589
3590    // compute current state
3591    int numberObjectInfeasibilities; // just odd ones
3592    int numberIntegerInfeasibilities;
3593    model->feasibleSolution(
3594        numberIntegerInfeasibilities,
3595        numberObjectInfeasibilities);
3596# ifdef COIN_HAS_CLP
3597    OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
3598    int saveClpOptions = 0;
3599    bool fastIterations = (model->specialOptions() & 8) != 0;
3600    if (osiclp && fastIterations) {
3601        // for faster hot start
3602        saveClpOptions = osiclp->specialOptions();
3603        osiclp->setSpecialOptions(saveClpOptions | 8192);
3604    }
3605# else
3606    bool fastIterations = false ;
3607# endif
3608    /*
3609      Scan for branching objects that indicate infeasibility. Choose candidates
3610      using priority as the first criteria, then integer infeasibility.
3611
3612      The algorithm is to fill the array with a set of good candidates (by
3613      infeasibility) with priority bestPriority.  Finding a candidate with
3614      priority better (less) than bestPriority flushes the choice array. (This
3615      serves as initialization when the first candidate is found.)
3616
3617    */
3618    numberToDo = 0;
3619    for (i = 0; i < numberObjects; i++) {
3620        OsiObject * object = model->modifiableObject(i);
3621        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3622            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3623        if (!dynamicObject)
3624            continue;
3625        int preferredWay;
3626        double infeasibility = object->infeasibility(&usefulInfo, preferredWay);
3627        int iColumn = dynamicObject->columnNumber();
3628        if (saveUpper[iColumn] == saveLower[iColumn])
3629            continue;
3630        if (infeasibility)
3631            sort[numberToDo] = -1.0e10 - infeasibility;
3632        else
3633            sort[numberToDo] = -fabs(dj[iColumn]);
3634        whichObject[numberToDo++] = i;
3635    }
3636    // Save basis
3637    CoinWarmStart * ws = solver->getWarmStart();
3638    int saveLimit;
3639    solver->getIntParam(OsiMaxNumIterationHotStart, saveLimit);
3640    int targetIterations = CoinMax(500, numberIterationsAllowed / numberObjects);
3641    if (saveLimit < targetIterations)
3642        solver->setIntParam(OsiMaxNumIterationHotStart, targetIterations);
3643    // Mark hot start
3644    solver->markHotStart();
3645    // Sort
3646    CoinSort_2(sort, sort + numberToDo, whichObject);
3647    double * currentSolution = model->currentSolution();
3648    double objMin = 1.0e50;
3649    double objMax = -1.0e50;
3650    bool needResolve = false;
3651    /*
3652      Now calculate the cost forcing the variable up and down.
3653    */
3654    int iDo;
3655    for (iDo = 0; iDo < numberToDo; iDo++) {
3656        CbcStrongInfo choice;
3657        int iObject = whichObject[iDo];
3658        OsiObject * object = model->modifiableObject(iObject);
3659        CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
3660            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
3661        if (!dynamicObject)
3662            continue;
3663        int iColumn = dynamicObject->columnNumber();
3664        int preferredWay;
3665        /*
3666          Update the information held in the object.
3667        */
3668        object->infeasibility(&usefulInfo, preferredWay);
3669        double value = currentSolution[iColumn];
3670        double nearest = floor(value + 0.5);
3671        double lowerValue = floor(value);
3672        bool satisfied = false;
3673        if (fabs(value - nearest) <= integerTolerance || value < saveLower[iColumn] || value > saveUpper[iColumn]) {
3674            satisfied = true;
3675            double newValue;
3676            if (nearest < saveUpper[iColumn]) {
3677                newValue = nearest + 1.0001 * integerTolerance;
3678                lowerValue = nearest;
3679            } else {
3680                newValue = nearest - 1.0001 * integerTolerance;
3681                lowerValue = nearest - 1;
3682            }
3683            currentSolution[iColumn] = newValue;
3684        }
3685        double upperValue = lowerValue + 1.0;
3686        //CbcSimpleInteger * obj =
3687        //dynamic_cast <CbcSimpleInteger *>(object) ;
3688        //if (obj) {
3689        //choice.possibleBranch=obj->createCbcBranch(solver,&usefulInfo,preferredWay);
3690        //} else {
3691        CbcObject * obj =
3692            dynamic_cast <CbcObject *>(object) ;
3693        assert (obj);
3694        choice.possibleBranch = obj->createCbcBranch(solver, &usefulInfo, preferredWay);
3695        //}
3696        currentSolution[iColumn] = value;
3697        // Save which object it was
3698        choice.objectNumber = iObject;
3699        choice.numIntInfeasUp = numberUnsatisfied_;
3700        choice.numIntInfeasDown = numberUnsatisfied_;
3701        choice.downMovement = 0.0;
3702        choice.upMovement = 0.0;
3703        choice.numItersDown = 0;
3704        choice.numItersUp = 0;
3705        choice.fix = 0; // say not fixed
3706        double objectiveChange ;
3707        double newObjectiveValue = 1.0e100;
3708        int j;
3709        // status is 0 finished, 1 infeasible and other
3710        int iStatus;
3711        /*
3712          Try the down direction first. (Specify the initial branching alternative as
3713          down with a call to way(-1). Each subsequent call to branch() performs the
3714          specified branch and advances the branch object state to the next branch
3715          alternative.)
3716        */
3717        choice.possibleBranch->way(-1) ;
3718        choice.possibleBranch->branch() ;
3719        if (fabs(value - lowerValue) > integerTolerance) {
3720            solver->solveFromHotStart() ;
3721            /*
3722              We now have an estimate of objective degradation that we can use for strong
3723              branching. If we're over the cutoff, the variable is monotone up.
3724              If we actually made it to optimality, check for a solution, and if we have
3725              a good one, call setBestSolution to process it. Note that this may reduce the
3726              cutoff, so we check again to see if we can declare this variable monotone.
3727            */
3728            if (solver->isProvenOptimal())
3729                iStatus = 0; // optimal
3730            else if (solver->isIterationLimitReached()
3731                     && !solver->isDualObjectiveLimitReached())
3732                iStatus = 2; // unknown
3733            else
3734                iStatus = 1; // infeasible
3735            newObjectiveValue = solver->getObjSense() * solver->getObjValue();
3736            choice.numItersDown = solver->getIterationCount();
3737            numberIterationsAllowed -= choice.numItersDown;
3738            objectiveChange = newObjectiveValue  - objectiveValue_;
3739            if (!iStatus) {
3740                choice.finishedDown = true ;
3741                if (newObjectiveValue >= cutoff) {
3742                    objectiveChange = 1.0e100; // say infeasible
3743                } else {
3744                    // See if integer solution
3745                    if (model->feasibleSolution(choice.numIntInfeasDown,
3746                                                choice.numObjInfeasDown)
3747                            && model->problemFeasibility()->feasible(model, -1) >= 0) {
3748                        model->setBestSolution(CBC_STRONGSOL,
3749                                               newObjectiveValue,
3750                                               solver->getColSolution()) ;
3751                        model->setLastHeuristic(NULL);
3752                        model->incrementUsed(solver->getColSolution());
3753                        cutoff = model->getCutoff();
3754                        if (newObjectiveValue >= cutoff)        //  *new* cutoff
3755                            objectiveChange = 1.0e100 ;
3756                    }
3757                }
3758            } else if (iStatus == 1) {
3759                objectiveChange = 1.0e100 ;
3760            } else {
3761                // Can't say much as we did not finish
3762                choice.finishedDown = false ;
3763            }
3764            choice.downMovement = objectiveChange ;
3765        }
3766        // restore bounds
3767        for ( j = 0; j < numberColumns; j++) {
3768            if (saveLower[j] != lower[j])
3769                solver->setColLower(j, saveLower[j]);
3770            if (saveUpper[j] != upper[j])
3771                solver->setColUpper(j, saveUpper[j]);
3772        }
3773        // repeat the whole exercise, forcing the variable up
3774        choice.possibleBranch->branch();
3775        if (fabs(value - upperValue) > integerTolerance) {
3776            solver->solveFromHotStart() ;
3777            /*
3778              We now have an estimate of objective degradation that we can use for strong
3779              branching. If we're over the cutoff, the variable is monotone up.
3780              If we actually made it to optimality, check for a solution, and if we have
3781              a good one, call setBestSolution to process it. Note that this may reduce the
3782              cutoff, so we check again to see if we can declare this variable monotone.
3783            */
3784            if (solver->isProvenOptimal())
3785                iStatus = 0; // optimal
3786            else if (solver->isIterationLimitReached()
3787                     && !solver->isDualObjectiveLimitReached())
3788                iStatus = 2; // unknown
3789            else
3790                iStatus = 1; // infeasible
3791            newObjectiveValue = solver->getObjSense() * solver->getObjValue();
3792            choice.numItersUp = solver->getIterationCount();
3793            numberIterationsAllowed -= choice.numItersUp;
3794            objectiveChange = newObjectiveValue  - objectiveValue_;
3795            if (!iStatus) {
3796                choice.finishedUp = true ;
3797                if (newObjectiveValue >= cutoff) {
3798                    objectiveChange = 1.0e100; // say infeasible
3799                } else {
3800                    // See if integer solution
3801                    if (model->feasibleSolution(choice.numIntInfeasUp,
3802                                                choice.numObjInfeasUp)
3803                            && model->problemFeasibility()->feasible(model, -1) >= 0) {
3804                        model->setBestSolution(CBC_STRONGSOL,
3805                                               newObjectiveValue,
3806                                               solver->getColSolution()) ;
3807                        model->setLastHeuristic(NULL);
3808                        model->incrementUsed(solver->getColSolution());
3809                        cutoff = model->getCutoff();
3810                        if (newObjectiveValue >= cutoff)        //  *new* cutoff
3811                            objectiveChange = 1.0e100 ;
3812                    }
3813                }
3814            } else if (iStatus == 1) {
3815                objectiveChange = 1.0e100 ;
3816            } else {
3817                // Can't say much as we did not finish
3818                choice.finishedUp = false ;
3819            }
3820            choice.upMovement = objectiveChange ;
3821
3822            // restore bounds
3823            for ( j = 0; j < numberColumns; j++) {
3824                if (saveLower[j] != lower[j])
3825                    solver->setColLower(j, saveLower[j]);
3826                if (saveUpper[j] != upper[j])
3827                    solver->setColUpper(j, saveUpper[j]);
3828            }
3829        }
3830        // If objective goes above certain amount we can set bound
3831        int jInt = back[iColumn];
3832        newLower[jInt] = upperValue;
3833        if (choice.finishedDown)
3834            objLower[jInt] = choice.downMovement + objectiveValue_;
3835        else
3836            objLower[jInt] = objectiveValue_;
3837        newUpper[jInt] = lowerValue;
3838        if (choice.finishedUp)
3839            objUpper[jInt] = choice.upMovement + objectiveValue_;
3840        else
3841            objUpper[jInt] = objectiveValue_;
3842        objMin = CoinMin(CoinMin(objLower[jInt], objUpper[jInt]), objMin);
3843        /*
3844          End of evaluation for this candidate variable. Possibilities are:
3845          * Both sides below cutoff; this variable is a candidate for branching.
3846          * Both sides infeasible or above the objective cutoff: no further action
3847          here. Break from the evaluation loop and assume the node will be purged
3848          by the caller.
3849          * One side below cutoff: Install the branch (i.e., fix the variable). Break
3850          from the evaluation loop and assume the node will be reoptimised by the
3851          caller.
3852        */
3853        if (choice.upMovement < 1.0e100) {
3854            if (choice.downMovement < 1.0e100) {
3855                objMax = CoinMax(CoinMax(objLower[jInt], objUpper[jInt]), objMax);
3856                // In case solution coming in was odd
3857                choice.upMovement = CoinMax(0.0, choice.upMovement);
3858                choice.downMovement = CoinMax(0.0, choice.downMovement);
3859                // feasible -
3860                model->messageHandler()->message(CBC_STRONG, *model->messagesPointer())
3861                << iObject << iColumn
3862                << choice.downMovement << choice.numIntInfeasDown
3863                << choice.upMovement << choice.numIntInfeasUp
3864                << value
3865                << CoinMessageEol;
3866            } else {
3867                // up feasible, down infeasible
3868                anyAction = -1;
3869                if (!satisfied)
3870                    needResolve = true;
3871                choice.fix = 1;
3872                numberToFix++;
3873                saveLower[iColumn] = upperValue;
3874                solver->setColLower(iColumn, upperValue);
3875            }
3876        } else {
3877            if (choice.downMovement < 1.0e100) {
3878                // down feasible, up infeasible
3879                anyAction = -1;
3880                if (!satisfied)
3881                    needResolve = true;
3882                choice.fix = -1;
3883                numberToFix++;
3884                saveUpper[iColumn] = lowerValue;
3885                solver->setColUpper(iColumn, lowerValue);
3886            } else {
3887                // neither side feasible
3888                anyAction = -2;
3889                COIN_DETAIL_PRINT(printf("Both infeasible for choice %d sequence %d\n", i,
3890                                         model->object(choice.objectNumber)->columnNumber()));
3891                delete ws;
3892                ws = NULL;
3893                //solver->writeMps("bad");
3894                numberToFix = -1;
3895                delete choice.possibleBranch;
3896                choice.possibleBranch = NULL;
3897                break;
3898            }
3899        }
3900        delete choice.possibleBranch;
3901        if (numberIterationsAllowed <= 0)
3902            break;
3903        //printf("obj %d, col %d, down %g up %g value %g\n",iObject,iColumn,
3904        //     choice.downMovement,choice.upMovement,value);
3905    }
3906    COIN_DETAIL_PRINT(printf("Best possible solution %g, can fix more if solution of %g found - looked at %d variables in %d iterations\n",
3907                             objMin, objMax, iDo, model->numberAnalyzeIterations() - numberIterationsAllowed));
3908    model->setNumberAnalyzeIterations(numberIterationsAllowed);
3909    // Delete the snapshot
3910    solver->unmarkHotStart();
3911    // back to normal
3912    solver->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
3913    solver->setIntParam(OsiMaxNumIterationHotStart, saveLimit);
3914    // restore basis
3915    solver->setWarmStart(ws);
3916    delete ws;
3917
3918    delete [] sort;
3919    delete [] whichObject;
3920    delete [] saveLower;
3921    delete [] saveUpper;
3922    delete [] back;
3923    // restore solution
3924    solver->setColSolution(saveSolution);
3925# ifdef COIN_HAS_CLP
3926    if (osiclp)
3927        osiclp->setSpecialOptions(saveClpOptions);
3928# endif
3929    model->reserveCurrentSolution(saveSolution);
3930    delete [] saveSolution;
3931    if (needResolve)
3932        solver->resolve();
3933    return numberToFix;
3934}
3935
3936
3937CbcNode::CbcNode(const CbcNode & rhs)
3938        : CoinTreeNode(rhs)
3939{
3940#ifdef CHECK_NODE
3941    printf("CbcNode %x Constructor from rhs %x\n", this, &rhs);
3942#endif
3943    if (rhs.nodeInfo_)
3944        nodeInfo_ = rhs.nodeInfo_->clone();
3945    else
3946        nodeInfo_ = NULL;
3947    objectiveValue_ = rhs.objectiveValue_;
3948    guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
3949    sumInfeasibilities_ = rhs.sumInfeasibilities_;
3950    if (rhs.branch_)
3951        branch_ = rhs.branch_->clone();
3952    else
3953        branch_ = NULL;
3954    depth_ = rhs.depth_;
3955    numberUnsatisfied_ = rhs.numberUnsatisfied_;
3956    nodeNumber_ = rhs.nodeNumber_;
3957    state_ = rhs.state_;
3958    if (nodeInfo_)
3959        assert ((state_&2) != 0);
3960    else
3961        assert ((state_&2) == 0);
3962}
3963
3964CbcNode &
3965CbcNode::operator=(const CbcNode & rhs)
3966{
3967    if (this != &rhs) {
3968        delete nodeInfo_;
3969        if (rhs.nodeInfo_)
3970            nodeInfo_ = rhs.nodeInfo_->clone();
3971        else
3972            nodeInfo_ = NULL;
3973        objectiveValue_ = rhs.objectiveValue_;
3974        guessedObjectiveValue_ = rhs.guessedObjectiveValue_;
3975        sumInfeasibilities_ = rhs.sumInfeasibilities_;
3976        if (rhs.branch_)
3977            branch_ = rhs.branch_->clone();
3978        else
3979            branch_ = NULL,
3980                      depth_ = rhs.depth_;
3981        numberUnsatisfied_ = rhs.numberUnsatisfied_;
3982        nodeNumber_ = rhs.nodeNumber_;
3983        state_ = rhs.state_;
3984        if (nodeInfo_)
3985            assert ((state_&2) != 0);
3986        else
3987            assert ((state_&2) == 0);
3988    }
3989    return *this;
3990}
3991CbcNode::~CbcNode ()
3992{
3993#ifdef CHECK_NODE
3994    if (nodeInfo_) {
3995        printf("CbcNode %x Destructor nodeInfo %x (%d)\n",
3996               this, nodeInfo_, nodeInfo_->numberPointingToThis());
3997        //assert(nodeInfo_->numberPointingToThis()>=0);
3998    } else {
3999        printf("CbcNode %x Destructor nodeInfo %x (?)\n",
4000               this, nodeInfo_);
4001    }
4002#endif
4003    if (nodeInfo_) {
4004        // was if (nodeInfo_&&(state_&2)!=0) {
4005        nodeInfo_->nullOwner();
4006        int numberToDelete = nodeInfo_->numberBranchesLeft();
4007        //    CbcNodeInfo * parent = nodeInfo_->parent();
4008        //assert (nodeInfo_->numberPointingToThis()>0);
4009        if (nodeInfo_->decrement(numberToDelete) == 0 || (state_&2) == 0) {
4010            if ((state_&2) == 0)
4011                nodeInfo_->nullParent();
4012            delete nodeInfo_;
4013        } else {
4014            //printf("node %x nodeinfo %x parent %x\n",this,nodeInfo_,nodeInfo_->parent());
4015            // anyway decrement parent
4016            //if (parent)
4017            ///parent->decrement(1);
4018        }
4019    }
4020    delete branch_;
4021}
4022// Decrement  active cut counts
4023void
4024CbcNode::decrementCuts(int change)
4025{
4026    if (nodeInfo_)
4027        assert ((state_&2) != 0);
4028    else
4029        assert ((state_&2) == 0);
4030    if (nodeInfo_) {
4031        nodeInfo_->decrementCuts(change);
4032    }
4033}
4034void
4035CbcNode::decrementParentCuts(CbcModel * model, int change)
4036{
4037    if (nodeInfo_)
4038        assert ((state_&2) != 0);
4039    else
4040        assert ((state_&2) == 0);
4041    if (nodeInfo_) {
4042        nodeInfo_->decrementParentCuts(model, change);
4043    }
4044}
4045
4046/*
4047  Initialize reference counts (numberPointingToThis, numberBranchesLeft_)
4048  in the attached nodeInfo_.
4049*/
4050void
4051CbcNode::initializeInfo()
4052{
4053    assert(nodeInfo_ && branch_) ;
4054    nodeInfo_->initializeInfo(branch_->numberBranches());
4055    assert ((state_&2) != 0);
4056    assert (nodeInfo_->numberBranchesLeft() ==
4057            branch_->numberBranchesLeft());
4058}
4059// Nulls out node info
4060void
4061CbcNode::nullNodeInfo()
4062{
4063    nodeInfo_ = NULL;
4064    // say not active
4065    state_ &= ~2;
4066}
4067
4068int
4069CbcNode::branch(OsiSolverInterface * solver)
4070{
4071    double changeInGuessed;
4072    assert (nodeInfo_->numberBranchesLeft() ==
4073            branch_->numberBranchesLeft());
4074    if (!solver)
4075        changeInGuessed = branch_->branch();
4076    else
4077        changeInGuessed = branch_->branch(solver);
4078    guessedObjectiveValue_ += changeInGuessed;
4079    //#define PRINTIT
4080#ifdef PRINTIT
4081    int numberLeft = nodeInfo_->numberBranchesLeft();
4082    CbcNodeInfo * parent = nodeInfo_->parent();
4083    int parentNodeNumber = -1;
4084    CbcBranchingObject * object1 =
4085        dynamic_cast<CbcBranchingObject *>(branch_) ;
4086    //OsiObject * object = object1->
4087    //int sequence = object->columnNumber);
4088    int id = -1;
4089    double value = 0.0;
4090    if (object1) {
4091        id = object1->variable();
4092        value = object1->value();
4093    }
4094    printf("id %d value %g objvalue %g\n", id, value, objectiveValue_);
4095    if (parent)
4096        parentNodeNumber = parent->nodeNumber();
4097    printf("Node number %d, %s, way %d, depth %d, parent node number %d\n",
4098           nodeInfo_->nodeNumber(), (numberLeft == 2) ? "leftBranch" : "rightBranch",
4099           way(), depth_, parentNodeNumber);
4100    assert (parentNodeNumber != nodeInfo_->nodeNumber());
4101#endif
4102    return nodeInfo_->branchedOn();
4103}
4104/* Active arm of the attached OsiBranchingObject.
4105
4106   In the simplest instance, coded -1 for the down arm of the branch, +1 for
4107   the up arm. But see OsiBranchingObject::way()
4108   Use nodeInfo--.numberBranchesLeft_ to see how active
4109
4110   Except that there is no OsiBranchingObject::way(), and this'll fail in any
4111   event because we have various OsiXXXBranchingObjects which aren't descended
4112   from CbcBranchingObjects. I think branchIndex() is the appropriate
4113   equivalent, but could be wrong. (lh, 061220)
4114
4115   071212: I'm finally getting back to cbc-generic and rescuing a lot of my
4116   annotation from branches/devel (which was killed in summer). I'm going to
4117   put back an assert(obj) just to see what happens. It's still present as of
4118   the most recent change to CbcNode (r833).
4119
4120   080104: Yep, we can arrive here with an OsiBranchingObject. Removed the
4121   assert, it's served its purpose.
4122
4123   080226: John finally noticed this problem and added a way() method to the
4124   OsiBranchingObject hierarchy. Removing my workaround.
4125
4126*/
4127int
4128CbcNode::way() const
4129{
4130    if (branch_) {
4131        CbcBranchingObject * obj =
4132            dynamic_cast <CbcBranchingObject *>(branch_) ;
4133        if (obj) {
4134            return obj->way();
4135        } else {
4136            OsiTwoWayBranchingObject * obj2 =
4137                dynamic_cast <OsiTwoWayBranchingObject *>(branch_) ;
4138            assert (obj2);
4139            return obj2->way();
4140        }
4141    } else {
4142        return 0;
4143    }
4144}
4145/* Create a branching object for the node
4146
4147   The routine scans the object list of the model and selects a set of
4148   unsatisfied objects as candidates for branching. The candidates are
4149   evaluated, and an appropriate branch object is installed.
4150
4151   The numberPassesLeft is decremented to stop fixing one variable each time
4152   and going on and on (e.g. for stock cutting, air crew scheduling)
4153
4154   If evaluation determines that an object is monotone or infeasible,
4155   the routine returns immediately. In the case of a monotone object,
4156   the branch object has already been called to modify the model.
4157
4158   Return value:
4159   <ul>
4160   <li>  0: A branching object has been installed
4161   <li> -1: A monotone object was discovered
4162   <li> -2: An infeasible object was discovered
4163   </ul>
4164   Branch state:
4165   <ul>
4166   <li> -1: start
4167   <li> -1: A monotone object was discovered
4168   <li> -2: An infeasible object was discovered
4169   </ul>
4170*/
4171int
4172CbcNode::chooseOsiBranch (CbcModel * model,
4173                          CbcNode * lastNode,
4174                          OsiBranchingInformation * usefulInfo,
4175                          int branchState)
4176{
4177    int returnStatus = 0;
4178    if (lastNode)
4179        depth_ = lastNode->depth_ + 1;
4180    else
4181        depth_ = 0;
4182    OsiSolverInterface * solver = model->solver();
4183    objectiveValue_ = solver->getObjValue() * solver->getObjSense();
4184    usefulInfo->objectiveValue_ = objectiveValue_;
4185    usefulInfo->depth_ = depth_;
4186    const double * saveInfoSol = usefulInfo->solution_;
4187    double * saveSolution = new double[solver->getNumCols()];
4188    memcpy(saveSolution, solver->getColSolution(), solver->getNumCols()*sizeof(double));
4189    usefulInfo->solution_ = saveSolution;
4190    OsiChooseVariable * choose = model->branchingMethod()->chooseMethod();
4191    int numberUnsatisfied = -1;
4192    if (branchState < 0) {
4193        // initialize
4194        // initialize sum of "infeasibilities"
4195        sumInfeasibilities_ = 0.0;
4196        numberUnsatisfied = choose->setupList(usefulInfo, true);
4197        numberUnsatisfied_ = numberUnsatisfied;
4198        branchState = 0;
4199        if (numberUnsatisfied_ < 0) {
4200            // infeasible
4201            delete [] saveSolution;
4202            return -2;
4203        }
4204    }
4205    // unset best
4206    int best = -1;
4207    choose->setBestObjectIndex(-1);
4208    if (numberUnsatisfied) {
4209        if (branchState > 0 || !choose->numberOnList()) {
4210            // we need to return at once - don't do strong branching or anything
4211            if (choose->numberOnList() || !choose->numberStrong()) {
4212                best = choose->candidates()[0];
4213                choose->setBestObjectIndex(best);
4214            } else {
4215                // nothing on list - need to try again - keep any solution
4216                numberUnsatisfied = choose->setupList(usefulInfo, false);
4217                numberUnsatisfied_ = numberUnsatisfied;
4218                if (numberUnsatisfied) {
4219                    best = choose->candidates()[0];
4220                    choose->setBestObjectIndex(best);
4221                }
4222            }
4223        } else {
4224            // carry on with strong branching or whatever
4225            int returnCode = choose->chooseVariable(solver, usefulInfo, true);
4226            // update number of strong iterations etc
4227            model->incrementStrongInfo(choose->numberStrongDone(), choose->numberStrongIterations(),
4228                                       returnCode == -1 ? 0 : choose->numberStrongFixed(), returnCode == -1);
4229            if (returnCode > 1) {
4230                // has fixed some
4231                returnStatus = -1;
4232            } else if (returnCode == -1) {
4233                // infeasible
4234                returnStatus = -2;
4235            } else if (returnCode == 0) {
4236                // normal
4237                returnStatus = 0;
4238                numberUnsatisfied = 1;
4239            } else {
4240                // ones on list satisfied - double check
4241                numberUnsatisfied = choose->setupList(usefulInfo, false);
4242                numberUnsatisfied_ = numberUnsatisfied;
4243                if (numberUnsatisfied) {
4244                    best = choose->candidates()[0];
4245                    choose->setBestObjectIndex(best);
4246                }
4247            }
4248        }
4249    }
4250    delete branch_;
4251    branch_ = NULL;
4252    guessedObjectiveValue_ = COIN_DBL_MAX;//objectiveValue_; // for now
4253    if (!returnStatus) {
4254        if (numberUnsatisfied) {
4255            // create branching object
4256            const OsiObject * obj = model->solver()->object(choose->bestObjectIndex());
4257            //const OsiSolverInterface * solver = usefulInfo->solver_;
4258            branch_ = obj->createBranch(model->solver(), usefulInfo, obj->whichWay());
4259        }
4260    }
4261    usefulInfo->solution_ = saveInfoSol;
4262    delete [] saveSolution;
4263    // may have got solution
4264    if (choose->goodSolution()
4265            && model->problemFeasibility()->feasible(model, -1) >= 0) {
4266        // yes
4267        double objValue = choose->goodObjectiveValue();
4268        model->setBestSolution(CBC_STRONGSOL,
4269                               objValue,
4270                               choose->goodSolution()) ;
4271        model->setLastHeuristic(NULL);
4272        model->incrementUsed(choose->goodSolution());
4273        choose->clearGoodSolution();
4274    }
4275    return returnStatus;
4276}
4277int
4278CbcNode::chooseClpBranch (CbcModel * model,
4279                          CbcNode * lastNode)
4280{
4281    assert(lastNode);
4282    depth_ = lastNode->depth_ + 1;
4283    delete branch_;
4284    branch_ = NULL;
4285    OsiSolverInterface * solver = model->solver();
4286    //double saveObjectiveValue = solver->getObjValue();
4287    //double objectiveValue = CoinMax(solver->getObjSense()*saveObjectiveValue,objectiveValue_);
4288    const double * lower = solver->getColLower();
4289    const double * upper = solver->getColUpper();
4290    // point to useful information
4291    OsiBranchingInformation usefulInfo = model->usefulInformation();
4292    // and modify
4293    usefulInfo.depth_ = depth_;
4294    int i;
4295    //bool beforeSolution = model->getSolutionCount()==0;
4296    int numberObjects = model->numberObjects();
4297    int numberColumns = model->getNumCols();
4298    double * saveUpper = new double[numberColumns];
4299    double * saveLower = new double[numberColumns];
4300
4301    // Save solution in case heuristics need good solution later
4302
4303    double * saveSolution = new double[numberColumns];
4304    memcpy(saveSolution, solver->getColSolution(), numberColumns*sizeof(double));
4305    model->reserveCurrentSolution(saveSolution);
4306    for (i = 0; i < numberColumns; i++) {
4307        saveLower[i] = lower[i];
4308        saveUpper[i] = upper[i];
4309    }
4310    // Save basis
4311    CoinWarmStart * ws = solver->getWarmStart();
4312    numberUnsatisfied_ = 0;
4313    // initialize sum of "infeasibilities"
4314    sumInfeasibilities_ = 0.0;
4315    // Note looks as if off end (hidden one)
4316    OsiObject * object = model->modifiableObject(numberObjects);
4317    CbcGeneralDepth * thisOne = dynamic_cast <CbcGeneralDepth *> (object);
4318    assert (thisOne);
4319    OsiClpSolverInterface * clpSolver
4320    = dynamic_cast<OsiClpSolverInterface *> (solver);
4321    assert (clpSolver);
4322    ClpSimplex * simplex = clpSolver->getModelPtr();
4323    int preferredWay;
4324    double infeasibility = object->infeasibility(&usefulInfo, preferredWay);
4325    if (thisOne->whichSolution() >= 0) {
4326        ClpNode * nodeInfo = thisOne->nodeInfo(thisOne->whichSolution());
4327        nodeInfo->applyNode(simplex, 2);
4328        int saveLogLevel = simplex->logLevel();
4329        simplex->setLogLevel(0);
4330        simplex->dual();
4331        simplex->setLogLevel(saveLogLevel);
4332        double cutoff = model->getCutoff();
4333        bool goodSolution = true;
4334        if (simplex->status()) {
4335            //simplex->writeMps("bad7.mps",2);
4336            if (nodeInfo->objectiveValue() > cutoff - 1.0e-2)
4337                goodSolution = false;
4338            else
4339                assert (!simplex->status());
4340        }
4341        if (goodSolution) {
4342            double newObjectiveValue = solver->getObjSense() * solver->getObjValue();
4343            // See if integer solution
4344            int numInf;
4345            int numInf2;
4346            bool gotSol = model->feasibleSolution(numInf, numInf2);
4347            if (!gotSol) {
4348              COIN_DETAIL_PRINT(printf("numinf %d\n", numInf));
4349                double * sol = simplex->primalColumnSolution();
4350                for (int i = 0; i < numberColumns; i++) {
4351                    if (simplex->isInteger(i)) {
4352                        double value = floor(sol[i] + 0.5);
4353                        if (fabs(value - sol[i]) > 1.0e-7) {
4354                          COIN_DETAIL_PRINT(printf("%d value %g\n", i, sol[i]));
4355                            if (fabs(value - sol[i]) < 1.0e-3) {
4356                                sol[i] = value;
4357                            }
4358                        }
4359                    }
4360                }
4361                simplex->writeMps("bad8.mps", 2);
4362                bool gotSol = model->feasibleSolution(numInf, numInf2);
4363                if (!gotSol)
4364                    assert (gotSol);
4365            }
4366            model->setBestSolution(CBC_STRONGSOL,
4367                                   newObjectiveValue,
4368                                   solver->getColSolution()) ;
4369            model->setLastHeuristic(NULL);
4370            model->incrementUsed(solver->getColSolution());
4371        }
4372    }
4373    // restore bounds
4374    {
4375        for (int j = 0; j < numberColumns; j++) {
4376            if (saveLower[j] != lower[j])
4377                solver->setColLower(j, saveLower[j]);
4378            if (saveUpper[j] != upper[j])
4379                solver->setColUpper(j, saveUpper[j]);
4380        }
4381    }
4382    // restore basis
4383    solver->setWarmStart(ws);
4384    delete ws;
4385    int anyAction;
4386    //#define CHECK_PATH
4387#ifdef CHECK_PATH
4388    extern int gotGoodNode_Z;
4389    if (gotGoodNode_Z >= 0)
4390        printf("good node %d %g\n", gotGoodNode_Z, infeasibility);
4391#endif
4392    if (infeasibility > 0.0) {
4393        if (infeasibility == COIN_DBL_MAX) {
4394            anyAction = -2; // infeasible
4395        } else {
4396            branch_ = thisOne->createCbcBranch(solver, &usefulInfo, preferredWay);
4397            // Set to firat one (and change when re-pushing)
4398            CbcGeneralBranchingObject * branch =
4399                dynamic_cast <CbcGeneralBranchingObject *> (branch_);
4400            branch->state(objectiveValue_, sumInfeasibilities_,
4401                          numberUnsatisfied_, 0);
4402            branch->setNode(this);
4403            anyAction = 0;
4404        }
4405    } else {
4406        anyAction = -1;
4407    }
4408#ifdef CHECK_PATH
4409    gotGoodNode_Z = -1;
4410#endif
4411    // Set guessed solution value
4412    guessedObjectiveValue_ = objectiveValue_ + 1.0e-5;
4413    delete [] saveLower;
4414    delete [] saveUpper;
4415
4416    // restore solution
4417    solver->setColSolution(saveSolution);
4418    delete [] saveSolution;
4419    return anyAction;
4420}
4421/* Double checks in case node can change its mind!
4422   Returns objective value
4423   Can change objective etc */
4424double
4425CbcNode::checkIsCutoff(double cutoff)
4426{
4427    branch_->checkIsCutoff(cutoff);
4428    return objectiveValue_;
4429}
4430
Note: See TracBrowser for help on using the repository browser.