source: stable/2.8/Cbc/src/CbcNode.cpp @ 1870

Last change on this file since 1870 was 1839, checked in by forrest, 6 years ago

multiple root solvers, stronger strong branching and cutoff as constraint

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