source: trunk/Cbc/src/CbcModel.cpp @ 2014

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

fix infeas status

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 739.1 KB
Line 
1/* $Id: CbcModel.cpp 2014 2014-03-09 17:42:30Z forrest $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#if defined(_MSC_VER)
7// Turn off compiler warning about long names
8#  pragma warning(disable:4786)
9#endif
10
11#include "CbcConfig.h"
12
13#include <string>
14//#define CBC_DEBUG 1
15//#define CHECK_CUT_COUNTS
16//#define CHECK_NODE
17//#define CHECK_NODE_FULL
18//#define NODE_LOG
19//#define GLOBAL_CUTS_JUST_POINTERS
20#ifdef CGL_DEBUG_GOMORY
21extern int gomory_try;
22#endif
23#include <cassert>
24#include <cmath>
25#include <cfloat>
26#ifdef COIN_HAS_CLP
27// include Presolve from Clp
28#include "ClpPresolve.hpp"
29#include "OsiClpSolverInterface.hpp"
30#include "ClpNode.hpp"
31#include "ClpDualRowDantzig.hpp"
32#include "ClpSimplexPrimal.hpp"
33#endif
34
35#include "CbcEventHandler.hpp"
36
37#include "OsiSolverInterface.hpp"
38#include "OsiAuxInfo.hpp"
39#include "OsiSolverBranch.hpp"
40#include "OsiChooseVariable.hpp"
41#include "CoinWarmStartBasis.hpp"
42#include "CoinPackedMatrix.hpp"
43#include "CoinHelperFunctions.hpp"
44#include "CbcBranchActual.hpp"
45#include "CbcBranchDynamic.hpp"
46#include "CbcHeuristic.hpp"
47#include "CbcHeuristicFPump.hpp"
48#include "CbcHeuristicRINS.hpp"
49#include "CbcHeuristicDive.hpp"
50#include "CbcModel.hpp"
51#include "CbcTreeLocal.hpp"
52#include "CbcStatistics.hpp"
53#include "CbcStrategy.hpp"
54#include "CbcMessage.hpp"
55#include "OsiRowCut.hpp"
56#include "OsiColCut.hpp"
57#include "OsiRowCutDebugger.hpp"
58#include "OsiCuts.hpp"
59#include "CbcCountRowCut.hpp"
60#include "CbcCutGenerator.hpp"
61#include "CbcFeasibilityBase.hpp"
62#include "CbcFathom.hpp"
63#include "CbcFullNodeInfo.hpp"
64// include Probing
65#include "CglProbing.hpp"
66#include "CglGomory.hpp"
67#include "CglTwomir.hpp"
68// include preprocessing
69#include "CglPreProcess.hpp"
70#include "CglDuplicateRow.hpp"
71#include "CglStored.hpp"
72#include "CglClique.hpp"
73#include "CglKnapsackCover.hpp"
74
75#include "CoinTime.hpp"
76#include "CoinMpsIO.hpp"
77
78#include "CbcCompareActual.hpp"
79#include "CbcTree.hpp"
80// This may be dummy
81#include "CbcThread.hpp"
82/* Various functions local to CbcModel.cpp */
83
84static void * doRootCbcThread(void * voidInfo);
85
86namespace {
87
88//-------------------------------------------------------------------
89// Returns the greatest common denominator of two
90// positive integers, a and b, found using Euclid's algorithm
91//-------------------------------------------------------------------
92static int gcd(int a, int b)
93{
94    int remainder = -1;
95    // make sure a<=b (will always remain so)
96    if (a > b) {
97        // Swap a and b
98        int temp = a;
99        a = b;
100        b = temp;
101    }
102    // if zero then gcd is nonzero (zero may occur in rhs of packed)
103    if (!a) {
104        if (b) {
105            return b;
106        } else {
107            printf("**** gcd given two zeros!!\n");
108            abort();
109        }
110    }
111    while (remainder) {
112        remainder = b % a;
113        b = a;
114        a = remainder;
115    }
116    return b;
117}
118
119
120
121#ifdef CHECK_NODE_FULL
122
123/*
124  Routine to verify that tree linkage is correct. The invariant that is tested
125  is
126
127  reference count = (number of actual references) + (number of branches left)
128
129  The routine builds a set of paired arrays, info and count, by traversing the
130  tree. Each CbcNodeInfo is recorded in info, and the number of times it is
131  referenced (via the parent field) is recorded in count. Then a final check is
132  made to see if the numberPointingToThis_ field agrees.
133*/
134
135void verifyTreeNodes (const CbcTree * branchingTree, const CbcModel &model)
136
137{
138    if (model.getNodeCount() == 661) return;
139    printf("*** CHECKING tree after %d nodes\n", model.getNodeCount()) ;
140
141    int j ;
142    int nNodes = branchingTree->size() ;
143# define MAXINFO 1000
144    int *count = new int [MAXINFO] ;
145    CbcNodeInfo **info = new CbcNodeInfo*[MAXINFO] ;
146    int nInfo = 0 ;
147    /*
148      Collect all CbcNodeInfo objects in info, by starting from each live node and
149      traversing back to the root. Nodes in the live set should have unexplored
150      branches remaining.
151
152      TODO: The `while (nodeInfo)' loop could be made to break on reaching a
153        common ancester (nodeInfo is found in info[k]). Alternatively, the
154        check could change to signal an error if nodeInfo is not found above a
155        common ancestor.
156    */
157    for (j = 0 ; j < nNodes ; j++) {
158        CbcNode *node = branchingTree->nodePointer(j) ;
159        if (!node)
160            continue;
161        CbcNodeInfo *nodeInfo = node->nodeInfo() ;
162        int change = node->nodeInfo()->numberBranchesLeft() ;
163        assert(change) ;
164        while (nodeInfo) {
165            int k ;
166            for (k = 0 ; k < nInfo ; k++) {
167                if (nodeInfo == info[k]) break ;
168            }
169            if (k == nInfo) {
170                assert(nInfo < MAXINFO) ;
171                nInfo++ ;
172                info[k] = nodeInfo ;
173                count[k] = 0 ;
174            }
175            nodeInfo = nodeInfo->parent() ;
176        }
177    }
178    /*
179      Walk the info array. For each nodeInfo, look up its parent in info and
180      increment the corresponding count.
181    */
182    for (j = 0 ; j < nInfo ; j++) {
183        CbcNodeInfo *nodeInfo = info[j] ;
184        nodeInfo = nodeInfo->parent() ;
185        if (nodeInfo) {
186            int k ;
187            for (k = 0 ; k < nInfo ; k++) {
188                if (nodeInfo == info[k]) break ;
189            }
190            assert (k < nInfo) ;
191            count[k]++ ;
192        }
193    }
194    /*
195      Walk the info array one more time and check that the invariant holds. The
196      number of references (numberPointingToThis()) should equal the sum of the
197      number of actual references (held in count[]) plus the number of potential
198      references (unexplored branches, numberBranchesLeft()).
199    */
200    for (j = 0; j < nInfo; j++) {
201        CbcNodeInfo * nodeInfo = info[j] ;
202        if (nodeInfo) {
203            int k ;
204            for (k = 0; k < nInfo; k++)
205                if (nodeInfo == info[k])
206                    break ;
207            printf("Nodeinfo %x - %d left, %d count\n",
208                   nodeInfo,
209                   nodeInfo->numberBranchesLeft(),
210                   nodeInfo->numberPointingToThis()) ;
211            assert(nodeInfo->numberPointingToThis() ==
212                   count[k] + nodeInfo->numberBranchesLeft()) ;
213        }
214    }
215
216    delete [] count ;
217    delete [] info ;
218
219    return ;
220}
221
222#endif  /* CHECK_NODE_FULL */
223
224
225
226#ifdef CHECK_CUT_COUNTS
227
228/*
229  Routine to verify that cut reference counts are correct.
230*/
231void verifyCutCounts (const CbcTree * branchingTree, CbcModel &model)
232
233{
234    printf("*** CHECKING cuts after %d nodes\n", model.getNodeCount()) ;
235
236    int j ;
237    int nNodes = branchingTree->size() ;
238
239    /*
240      cut.tempNumber_ exists for the purpose of doing this verification. Clear it
241      in all cuts. We traverse the tree by starting from each live node and working
242      back to the root. At each CbcNodeInfo, check for cuts.
243    */
244    for (j = 0 ; j < nNodes ; j++) {
245        CbcNode *node = branchingTree->nodePointer(j) ;
246        CbcNodeInfo * nodeInfo = node->nodeInfo() ;
247        assert (node->nodeInfo()->numberBranchesLeft()) ;
248        while (nodeInfo) {
249            int k ;
250            for (k = 0 ; k < nodeInfo->numberCuts() ; k++) {
251                CbcCountRowCut *cut = nodeInfo->cuts()[k] ;
252                if (cut) cut->tempNumber_ = 0;
253            }
254            nodeInfo = nodeInfo->parent() ;
255        }
256    }
257    /*
258      Walk the live set again, this time collecting the list of cuts in use at each
259      node. addCuts1 will collect the cuts in model.addedCuts_. Take into account
260      that when we recreate the basis for a node, we compress out the slack cuts.
261    */
262    for (j = 0 ; j < nNodes ; j++) {
263        CoinWarmStartBasis *debugws = model.getEmptyBasis() ;
264        CbcNode *node = branchingTree->nodePointer(j) ;
265        CbcNodeInfo *nodeInfo = node->nodeInfo();
266        int change = node->nodeInfo()->numberBranchesLeft() ;
267        printf("Node %d %x (info %x) var %d way %d obj %g", j, node,
268               node->nodeInfo(), node->columnNumber(), node->way(),
269               node->objectiveValue()) ;
270
271        model.addCuts1(node, debugws) ;
272
273        int i ;
274        int numberRowsAtContinuous = model.numberRowsAtContinuous() ;
275        CbcCountRowCut **addedCuts = model.addedCuts() ;
276        for (i = 0 ; i < model.currentNumberCuts() ; i++) {
277            CoinWarmStartBasis::Status status =
278                debugws->getArtifStatus(i + numberRowsAtContinuous) ;
279            if (status != CoinWarmStartBasis::basic && addedCuts[i]) {
280                addedCuts[i]->tempNumber_ += change ;
281            }
282        }
283
284        while (nodeInfo) {
285            nodeInfo = nodeInfo->parent() ;
286            if (nodeInfo) printf(" -> %x", nodeInfo);
287        }
288        printf("\n") ;
289        delete debugws ;
290    }
291    /*
292      The moment of truth: We've tallied up the references by direct scan of the  search tree. Check for agreement with the count in the cut.
293
294      TODO: Rewrite to check and print mismatch only when tempNumber_ == 0?
295    */
296    for (j = 0 ; j < nNodes ; j++) {
297        CbcNode *node = branchingTree->nodePointer(j) ;
298        CbcNodeInfo *nodeInfo = node->nodeInfo();
299        while (nodeInfo) {
300            int k ;
301            for (k = 0 ; k < nodeInfo->numberCuts() ; k++) {
302                CbcCountRowCut *cut = nodeInfo->cuts()[k] ;
303                if (cut && cut->tempNumber_ >= 0) {
304                    if (cut->tempNumber_ != cut->numberPointingToThis())
305                        printf("mismatch %x %d %x %d %d\n", nodeInfo, k,
306                               cut, cut->tempNumber_, cut->numberPointingToThis()) ;
307                    else
308                        printf("   match %x %d %x %d %d\n", nodeInfo, k,
309                               cut, cut->tempNumber_, cut->numberPointingToThis()) ;
310                    cut->tempNumber_ = -1 ;
311                }
312            }
313            nodeInfo = nodeInfo->parent() ;
314        }
315    }
316
317    return ;
318}
319
320#endif /* CHECK_CUT_COUNTS */
321
322
323#ifdef CHECK_CUT_SIZE
324
325/*
326  Routine to verify that cut reference counts are correct.
327*/
328void verifyCutSize (const CbcTree * branchingTree, CbcModel &model)
329{
330
331    int j ;
332    int nNodes = branchingTree->size() ;
333    int totalCuts = 0;
334
335    /*
336      cut.tempNumber_ exists for the purpose of doing this verification. Clear it
337      in all cuts. We traverse the tree by starting from each live node and working
338      back to the root. At each CbcNodeInfo, check for cuts.
339    */
340    for (j = 0 ; j < nNodes ; j++) {
341        CbcNode *node = branchingTree->nodePointer(j) ;
342        CbcNodeInfo * nodeInfo = node->nodeInfo() ;
343        assert (node->nodeInfo()->numberBranchesLeft()) ;
344        while (nodeInfo) {
345            totalCuts += nodeInfo->numberCuts();
346            nodeInfo = nodeInfo->parent() ;
347        }
348    }
349    printf("*** CHECKING cuts (size) after %d nodes - %d cuts\n", model.getNodeCount(), totalCuts) ;
350    return ;
351}
352
353#endif /* CHECK_CUT_SIZE */
354
355}
356
357/* End unnamed namespace for CbcModel.cpp */
358
359
360void
361CbcModel::analyzeObjective ()
362/*
363  Try to find a minimum change in the objective function. The first scan
364  checks that there are no continuous variables with non-zero coefficients,
365  and grabs the largest objective coefficient associated with an unfixed
366  integer variable. The second scan attempts to scale up the objective
367  coefficients to a point where they are sufficiently close to integer that
368  we can pretend they are integer, and calculate a gcd over the coefficients
369  of interest. This will be the minimum increment for the scaled coefficients.
370  The final action is to scale the increment back for the original coefficients
371  and install it, if it's better than the existing value.
372
373  John's note: We could do better than this.
374
375  John's second note - apologies for changing s to z
376*/
377{
378    const double *objective = getObjCoefficients() ;
379    const double *lower = getColLower() ;
380    const double *upper = getColUpper() ;
381    /*
382      Scan continuous and integer variables to see if continuous
383      are cover or network with integral rhs.
384    */
385    double continuousMultiplier = 1.0;
386    double * coeffMultiplier = NULL;
387    double largestObj = 0.0;
388    double smallestObj = COIN_DBL_MAX;
389    {
390        const double *rowLower = getRowLower() ;
391        const double *rowUpper = getRowUpper() ;
392        int numberRows = solver_->getNumRows() ;
393        double * rhs = new double [numberRows];
394        memset(rhs, 0, numberRows*sizeof(double));
395        int iColumn;
396        int numberColumns = solver_->getNumCols() ;
397        // Column copy of matrix
398        int problemType = -1;
399        const double * element = solver_->getMatrixByCol()->getElements();
400        const int * row = solver_->getMatrixByCol()->getIndices();
401        const CoinBigIndex * columnStart = solver_->getMatrixByCol()->getVectorStarts();
402        const int * columnLength = solver_->getMatrixByCol()->getVectorLengths();
403        int numberInteger = 0;
404        int numberIntegerObj = 0;
405        int numberGeneralIntegerObj = 0;
406        int numberIntegerWeight = 0;
407        int numberContinuousObj = 0;
408        double cost = COIN_DBL_MAX;
409        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
410            if (upper[iColumn] == lower[iColumn]) {
411                CoinBigIndex start = columnStart[iColumn];
412                CoinBigIndex end = start + columnLength[iColumn];
413                for (CoinBigIndex j = start; j < end; j++) {
414                    int iRow = row[j];
415                    rhs[iRow] += lower[iColumn] * element[j];
416                }
417            } else {
418                double objValue = objective[iColumn];
419                if (solver_->isInteger(iColumn))
420                    numberInteger++;
421                if (objValue) {
422                    if (!solver_->isInteger(iColumn)) {
423                        numberContinuousObj++;
424                    } else {
425                        largestObj = CoinMax(largestObj, fabs(objValue));
426                        smallestObj = CoinMin(smallestObj, fabs(objValue));
427                        numberIntegerObj++;
428                        if (cost == COIN_DBL_MAX)
429                            cost = objValue;
430                        else if (cost != objValue)
431                            cost = -COIN_DBL_MAX;
432                        int gap = static_cast<int> (upper[iColumn] - lower[iColumn]);
433                        if (gap > 1) {
434                            numberGeneralIntegerObj++;
435                            numberIntegerWeight += gap;
436                        }
437                    }
438                }
439            }
440        }
441        int iType = 0;
442        if (!numberContinuousObj && numberIntegerObj <= 5 && numberIntegerWeight <= 100 &&
443                numberIntegerObj*3 < numberObjects_ && !parentModel_ && solver_->getNumRows() > 100)
444          iType = 1 + 4 + ((moreSpecialOptions_&536870912)==0) ? 2 : 0;
445        else if (!numberContinuousObj && numberIntegerObj <= 100 &&
446                 numberIntegerObj*5 < numberObjects_ && numberIntegerWeight <= 100 &&
447                 !parentModel_ &&
448                 solver_->getNumRows() > 100 && cost != -COIN_DBL_MAX)
449          iType = 4 + ((moreSpecialOptions_&536870912)==0) ? 2 : 0;
450        else if (!numberContinuousObj && numberIntegerObj <= 100 &&
451                 numberIntegerObj*5 < numberObjects_ &&
452                 !parentModel_ &&
453                 solver_->getNumRows() > 100 && cost != -COIN_DBL_MAX)
454            iType = 8;
455        int iTest = getMaximumNodes();
456        if (iTest >= 987654320 && iTest < 987654330 && numberObjects_ && !parentModel_) {
457            iType = iTest - 987654320;
458            printf("Testing %d integer variables out of %d objects (%d integer) have cost of %g - %d continuous\n",
459                   numberIntegerObj, numberObjects_, numberInteger, cost, numberContinuousObj);
460            if (iType == 9)
461                exit(77);
462            if (numberContinuousObj)
463                iType = 0;
464        }
465
466        //if (!numberContinuousObj&&(numberIntegerObj<=5||cost!=-COIN_DBL_MAX)&&
467        //numberIntegerObj*3<numberObjects_&&!parentModel_&&solver_->getNumRows()>100) {
468        if (iType) {
469            /*
470            A) put high priority on (if none)
471            B) create artificial objective (if clp)
472            */
473            int iPriority = -1;
474            for (int i = 0; i < numberObjects_; i++) {
475                int k = object_[i]->priority();
476                if (iPriority == -1)
477                    iPriority = k;
478                else if (iPriority != k)
479                    iPriority = -2;
480            }
481            bool branchOnSatisfied = ((iType & 1) != 0);
482            bool createFake = ((iType & 2) != 0);
483            bool randomCost = ((iType & 4) != 0);
484            if (iPriority >= 0) {
485                char general[200];
486                if (cost == -COIN_DBL_MAX) {
487                    sprintf(general, "%d integer variables out of %d objects (%d integer) have costs - high priority",
488                            numberIntegerObj, numberObjects_, numberInteger);
489                } else if (cost == COIN_DBL_MAX) {
490                    sprintf(general, "No integer variables out of %d objects (%d integer) have costs",
491                            numberObjects_, numberInteger);
492                    branchOnSatisfied = false;
493                } else {
494                    sprintf(general, "%d integer variables out of %d objects (%d integer) have cost of %g - high priority",
495                            numberIntegerObj, numberObjects_, numberInteger, cost);
496                }
497                messageHandler()->message(CBC_GENERAL,
498                                          messages())
499                << general << CoinMessageEol ;
500                sprintf(general, "branch on satisfied %c create fake objective %c random cost %c",
501                        branchOnSatisfied ? 'Y' : 'N',
502                        createFake ? 'Y' : 'N',
503                        randomCost ? 'Y' : 'N');
504                messageHandler()->message(CBC_GENERAL,
505                                          messages())
506                << general << CoinMessageEol ;
507                // switch off clp type branching
508                // no ? fastNodeDepth_ = -1;
509                int highPriority = (branchOnSatisfied) ? -999 : 100;
510                for (int i = 0; i < numberObjects_; i++) {
511                    CbcSimpleInteger * thisOne = dynamic_cast <CbcSimpleInteger *> (object_[i]);
512                    object_[i]->setPriority(1000);
513                    if (thisOne) {
514                        int iColumn = thisOne->columnNumber();
515                        if (objective[iColumn])
516                            thisOne->setPriority(highPriority);
517                    }
518                }
519            }
520#ifdef COIN_HAS_CLP
521            OsiClpSolverInterface * clpSolver
522            = dynamic_cast<OsiClpSolverInterface *> (solver_);
523            if (clpSolver && createFake) {
524                // Create artificial objective to be used when all else fixed
525                int numberColumns = clpSolver->getNumCols();
526                double * fakeObj = new double [numberColumns];
527                // Column copy
528                const CoinPackedMatrix  * matrixByCol = clpSolver->getMatrixByCol();
529                //const double * element = matrixByCol.getElements();
530                //const int * row = matrixByCol.getIndices();
531                //const CoinBigIndex * columnStart = matrixByCol.getVectorStarts();
532                const int * columnLength = matrixByCol->getVectorLengths();
533                const double * solution = clpSolver->getColSolution();
534#ifdef JJF_ZERO
535                int nAtBound = 0;
536                for (int i = 0; i < numberColumns; i++) {
537                    double lowerValue = lower[i];
538                    double upperValue = upper[i];
539                    if (clpSolver->isInteger(i)) {
540                        double lowerValue = lower[i];
541                        double upperValue = upper[i];
542                        double value = solution[i];
543                        if (value < lowerValue + 1.0e-6 ||
544                                value > upperValue - 1.0e-6)
545                            nAtBound++;
546                    }
547                }
548#endif
549                /*
550                  Generate a random objective function for problems where the given objective
551                  function is not terribly useful. (Nearly feasible, single integer variable,
552                  that sort of thing.
553                */
554                CoinDrand48(true, 1234567);
555                for (int i = 0; i < numberColumns; i++) {
556                    double lowerValue = lower[i];
557                    double upperValue = upper[i];
558                    double value = (randomCost) ? ceil((CoinDrand48() + 0.5) * 1000)
559                                   : i + 1 + columnLength[i] * 1000;
560                    value *= 0.001;
561                    //value += columnLength[i];
562                    if (lowerValue > -1.0e5 || upperValue < 1.0e5) {
563                        if (fabs(lowerValue) > fabs(upperValue))
564                            value = - value;
565                        if (clpSolver->isInteger(i)) {
566                            double solValue = solution[i];
567                            // Better to add in 0.5 or 1.0??
568                            if (solValue < lowerValue + 1.0e-6)
569                                value = fabs(value) + 0.5; //fabs(value*1.5);
570                            else if (solValue > upperValue - 1.0e-6)
571                                value = -fabs(value) - 0.5; //-fabs(value*1.5);
572                        }
573                    } else {
574                        value = 0.0;
575                    }
576                    fakeObj[i] = value;
577                }
578                // pass to solver
579                clpSolver->setFakeObjective(fakeObj);
580                delete [] fakeObj;
581            }
582#endif
583        } else if (largestObj < smallestObj*5.0 && !parentModel_ &&
584                   !numberContinuousObj &&
585                   !numberGeneralIntegerObj &&
586                   numberIntegerObj*2 < numberColumns) {
587            // up priorities on costed
588            int iPriority = -1;
589            for (int i = 0; i < numberObjects_; i++) {
590                int k = object_[i]->priority();
591                if (iPriority == -1)
592                    iPriority = k;
593                else if (iPriority != k)
594                    iPriority = -2;
595            }
596            if (iPriority >= 100) {
597#ifdef CLP_INVESTIGATE
598                printf("Setting variables with obj to high priority\n");
599#endif
600                for (int i = 0; i < numberObjects_; i++) {
601                    CbcSimpleInteger * obj =
602                        dynamic_cast <CbcSimpleInteger *>(object_[i]) ;
603                    if (obj) {
604                        int iColumn = obj->columnNumber();
605                        if (objective[iColumn])
606                            object_[i]->setPriority(iPriority - 1);
607                    }
608                }
609            }
610        }
611        int iRow;
612        for (iRow = 0; iRow < numberRows; iRow++) {
613            if (rowLower[iRow] > -1.0e20 &&
614                    fabs(rowLower[iRow] - rhs[iRow] - floor(rowLower[iRow] - rhs[iRow] + 0.5)) > 1.0e-10) {
615                continuousMultiplier = 0.0;
616                break;
617            }
618            if (rowUpper[iRow] < 1.0e20 &&
619                    fabs(rowUpper[iRow] - rhs[iRow] - floor(rowUpper[iRow] - rhs[iRow] + 0.5)) > 1.0e-10) {
620                continuousMultiplier = 0.0;
621                break;
622            }
623            // set rhs to limiting value
624            if (rowLower[iRow] != rowUpper[iRow]) {
625                if (rowLower[iRow] > -1.0e20) {
626                    if (rowUpper[iRow] < 1.0e20) {
627                        // no good
628                        continuousMultiplier = 0.0;
629                        break;
630                    } else {
631                        rhs[iRow] = rowLower[iRow] - rhs[iRow];
632                        if (problemType < 0)
633                            problemType = 3; // set cover
634                        else if (problemType != 3)
635                            problemType = 4;
636                    }
637                } else {
638                    rhs[iRow] = rowUpper[iRow] - rhs[iRow];
639                    if (problemType < 0)
640                        problemType = 1; // set partitioning <=
641                    else if (problemType != 1)
642                        problemType = 4;
643                }
644            } else {
645                rhs[iRow] = rowUpper[iRow] - rhs[iRow];
646                if (problemType < 0)
647                    problemType = 3; // set partitioning ==
648                else if (problemType != 2)
649                    problemType = 2;
650            }
651            if (fabs(rhs[iRow] - 1.0) > 1.0e-12)
652                problemType = 4;
653        }
654        if (continuousMultiplier) {
655            // 1 network, 2 cover, 4 negative cover
656            int possible = 7;
657            bool unitRhs = true;
658            // See which rows could be set cover
659            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
660                if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
661                    CoinBigIndex start = columnStart[iColumn];
662                    CoinBigIndex end = start + columnLength[iColumn];
663                    for (CoinBigIndex j = start; j < end; j++) {
664                        double value = element[j];
665                        if (value == 1.0) {
666                        } else if (value == -1.0) {
667                            rhs[row[j]] = -0.5;
668                        } else {
669                            rhs[row[j]] = -COIN_DBL_MAX;
670                        }
671                    }
672                }
673            }
674            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
675                if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
676                    if (!isInteger(iColumn)) {
677                        CoinBigIndex start = columnStart[iColumn];
678                        CoinBigIndex end = start + columnLength[iColumn];
679                        double rhsValue = 0.0;
680                        // 1 all ones, -1 all -1s, 2 all +- 1, 3 no good
681                        int type = 0;
682                        for (CoinBigIndex j = start; j < end; j++) {
683                            double value = element[j];
684                            if (fabs(value) != 1.0) {
685                                type = 3;
686                                break;
687                            } else if (value == 1.0) {
688                                if (!type)
689                                    type = 1;
690                                else if (type != 1)
691                                    type = 2;
692                            } else {
693                                if (!type)
694                                    type = -1;
695                                else if (type != -1)
696                                    type = 2;
697                            }
698                            int iRow = row[j];
699                            if (rhs[iRow] == -COIN_DBL_MAX) {
700                                type = 3;
701                                break;
702                            } else if (rhs[iRow] == -0.5) {
703                                // different values
704                                unitRhs = false;
705                            } else if (rhsValue) {
706                                if (rhsValue != rhs[iRow])
707                                    unitRhs = false;
708                            } else {
709                                rhsValue = rhs[iRow];
710                            }
711                        }
712                        // if no elements OK
713                        if (type == 3) {
714                            // no good
715                            possible = 0;
716                            break;
717                        } else if (type == 2) {
718                            if (end - start > 2) {
719                                // no good
720                                possible = 0;
721                                break;
722                            } else {
723                                // only network
724                                possible &= 1;
725                                if (!possible)
726                                    break;
727                            }
728                        } else if (type == 1) {
729                            // only cover
730                            possible &= 2;
731                            if (!possible)
732                                break;
733                        } else if (type == -1) {
734                            // only negative cover
735                            possible &= 4;
736                            if (!possible)
737                                break;
738                        }
739                    }
740                }
741            }
742            if ((possible == 2 || possible == 4) && !unitRhs) {
743#if COIN_DEVELOP>1
744                printf("XXXXXX Continuous all +1 but different rhs\n");
745#endif
746                possible = 0;
747            }
748            // may be all integer
749            if (possible != 7) {
750                if (!possible)
751                    continuousMultiplier = 0.0;
752                else if (possible == 1)
753                    continuousMultiplier = 1.0;
754                else
755                    continuousMultiplier = 0.0; // 0.5 was incorrect;
756#if COIN_DEVELOP>1
757                if (continuousMultiplier)
758                    printf("XXXXXX multiplier of %g\n", continuousMultiplier);
759#endif
760                if (continuousMultiplier == 0.5) {
761                    coeffMultiplier = new double [numberColumns];
762                    bool allOne = true;
763                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
764                        coeffMultiplier[iColumn] = 1.0;
765                        if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
766                            if (!isInteger(iColumn)) {
767                                CoinBigIndex start = columnStart[iColumn];
768                                int iRow = row[start];
769                                double value = rhs[iRow];
770                                assert (value >= 0.0);
771                                if (value != 0.0 && value != 1.0)
772                                    allOne = false;
773                                coeffMultiplier[iColumn] = 0.5 * value;
774                            }
775                        }
776                    }
777                    if (allOne) {
778                        // back to old way
779                        delete [] coeffMultiplier;
780                        coeffMultiplier = NULL;
781                    }
782                }
783            } else {
784                // all integer
785                problemType_ = problemType;
786#if COIN_DEVELOP>1
787                printf("Problem type is %d\n", problemType_);
788#endif
789            }
790        }
791
792        // But try again
793        if (continuousMultiplier < 1.0) {
794            memset(rhs, 0, numberRows*sizeof(double));
795            int * count = new int [numberRows];
796            memset(count, 0, numberRows*sizeof(int));
797            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
798                CoinBigIndex start = columnStart[iColumn];
799                CoinBigIndex end = start + columnLength[iColumn];
800                if (upper[iColumn] == lower[iColumn]) {
801                    for (CoinBigIndex j = start; j < end; j++) {
802                        int iRow = row[j];
803                        rhs[iRow] += lower[iColumn] * element[j];
804                    }
805                } else if (solver_->isInteger(iColumn)) {
806                    for (CoinBigIndex j = start; j < end; j++) {
807                        int iRow = row[j];
808                        if (fabs(element[j] - floor(element[j] + 0.5)) > 1.0e-10)
809                            rhs[iRow]  = COIN_DBL_MAX;
810                    }
811                } else {
812                    for (CoinBigIndex j = start; j < end; j++) {
813                        int iRow = row[j];
814                        count[iRow]++;
815                        if (fabs(element[j]) != 1.0)
816                            rhs[iRow]  = COIN_DBL_MAX;
817                    }
818                }
819            }
820            // now look at continuous
821            bool allGood = true;
822            double direction = solver_->getObjSense() ;
823            int numberObj = 0;
824            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
825                if (upper[iColumn] > lower[iColumn]) {
826                    double objValue = objective[iColumn] * direction;
827                    if (objValue && !solver_->isInteger(iColumn)) {
828                        numberObj++;
829                        CoinBigIndex start = columnStart[iColumn];
830                        CoinBigIndex end = start + columnLength[iColumn];
831                        if (objValue > 0.0) {
832                            // wants to be as low as possible
833                            if (lower[iColumn] < -1.0e10 || fabs(lower[iColumn] - floor(lower[iColumn] + 0.5)) > 1.0e-10) {
834                                allGood = false;
835                                break;
836                            } else if (upper[iColumn] < 1.0e10 && fabs(upper[iColumn] - floor(upper[iColumn] + 0.5)) > 1.0e-10) {
837                                allGood = false;
838                                break;
839                            }
840                            bool singletonRow = true;
841                            bool equality = false;
842                            for (CoinBigIndex j = start; j < end; j++) {
843                                int iRow = row[j];
844                                if (count[iRow] > 1)
845                                    singletonRow = false;
846                                else if (rowLower[iRow] == rowUpper[iRow])
847                                    equality = true;
848                                double rhsValue = rhs[iRow];
849                                double lowerValue = rowLower[iRow];
850                                double upperValue = rowUpper[iRow];
851                                if (rhsValue < 1.0e20) {
852                                    if (lowerValue > -1.0e20)
853                                        lowerValue -= rhsValue;
854                                    if (upperValue < 1.0e20)
855                                        upperValue -= rhsValue;
856                                }
857                                if (fabs(rhsValue) > 1.0e20 || fabs(rhsValue - floor(rhsValue + 0.5)) > 1.0e-10
858                                        || fabs(element[j]) != 1.0) {
859                                    // no good
860                                    allGood = false;
861                                    break;
862                                }
863                                if (element[j] > 0.0) {
864                                    if (lowerValue > -1.0e20 && fabs(lowerValue - floor(lowerValue + 0.5)) > 1.0e-10) {
865                                        // no good
866                                        allGood = false;
867                                        break;
868                                    }
869                                } else {
870                                    if (upperValue < 1.0e20 && fabs(upperValue - floor(upperValue + 0.5)) > 1.0e-10) {
871                                        // no good
872                                        allGood = false;
873                                        break;
874                                    }
875                                }
876                            }
877                            if (!singletonRow && end > start + 1 && !equality)
878                                allGood = false;
879                            if (!allGood)
880                                break;
881                        } else {
882                            // wants to be as high as possible
883                            if (upper[iColumn] > 1.0e10 || fabs(upper[iColumn] - floor(upper[iColumn] + 0.5)) > 1.0e-10) {
884                                allGood = false;
885                                break;
886                            } else if (lower[iColumn] > -1.0e10 && fabs(lower[iColumn] - floor(lower[iColumn] + 0.5)) > 1.0e-10) {
887                                allGood = false;
888                                break;
889                            }
890                            bool singletonRow = true;
891                            bool equality = false;
892                            for (CoinBigIndex j = start; j < end; j++) {
893                                int iRow = row[j];
894                                if (count[iRow] > 1)
895                                    singletonRow = false;
896                                else if (rowLower[iRow] == rowUpper[iRow])
897                                    equality = true;
898                                double rhsValue = rhs[iRow];
899                                double lowerValue = rowLower[iRow];
900                                double upperValue = rowUpper[iRow];
901                                if (rhsValue < 1.0e20) {
902                                    if (lowerValue > -1.0e20)
903                                        lowerValue -= rhsValue;
904                                    if (upperValue < 1.0e20)
905                                        upperValue -= rhsValue;
906                                }
907                                if (fabs(rhsValue) > 1.0e20 || fabs(rhsValue - floor(rhsValue + 0.5)) > 1.0e-10
908                                        || fabs(element[j]) != 1.0) {
909                                    // no good
910                                    allGood = false;
911                                    break;
912                                }
913                                if (element[j] < 0.0) {
914                                    if (lowerValue > -1.0e20 && fabs(lowerValue - floor(lowerValue + 0.5)) > 1.0e-10) {
915                                        // no good
916                                        allGood = false;
917                                        break;
918                                    }
919                                } else {
920                                    if (upperValue < 1.0e20 && fabs(upperValue - floor(upperValue + 0.5)) > 1.0e-10) {
921                                        // no good
922                                        allGood = false;
923                                        break;
924                                    }
925                                }
926                            }
927                            if (!singletonRow && end > start + 1 && !equality)
928                                allGood = false;
929                            if (!allGood)
930                                break;
931                        }
932                    }
933                }
934            }
935            delete [] count;
936            if (allGood) {
937#if COIN_DEVELOP>1
938                if (numberObj)
939                    printf("YYYY analysis says all continuous with costs will be integer\n");
940#endif
941                continuousMultiplier = 1.0;
942            }
943        }
944        delete [] rhs;
945    }
946    /*
947      Take a first scan to see if there are unfixed continuous variables in the
948      objective.  If so, the minimum objective change could be arbitrarily small.
949      Also pick off the maximum coefficient of an unfixed integer variable.
950
951      If the objective is found to contain only integer variables, set the
952      fathoming discipline to strict.
953    */
954    double maximumCost = 0.0 ;
955    //double trueIncrement=0.0;
956    int iColumn ;
957    int numberColumns = getNumCols() ;
958    double scaleFactor = 1.0; // due to rhs etc
959    /*
960      Original model did not have integer bounds.
961    */
962    if ((specialOptions_&65536) == 0) {
963        /* be on safe side (later look carefully as may be able to
964           to get 0.5 say if bounds are multiples of 0.5 */
965        for (iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
966            if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
967                double value;
968                value = fabs(lower[iColumn]);
969                if (floor(value + 0.5) != value) {
970                    scaleFactor = CoinMin(scaleFactor, 0.5);
971                    if (floor(2.0*value + 0.5) != 2.0*value) {
972                        scaleFactor = CoinMin(scaleFactor, 0.25);
973                        if (floor(4.0*value + 0.5) != 4.0*value) {
974                            scaleFactor = 0.0;
975                        }
976                    }
977                }
978                value = fabs(upper[iColumn]);
979                if (floor(value + 0.5) != value) {
980                    scaleFactor = CoinMin(scaleFactor, 0.5);
981                    if (floor(2.0*value + 0.5) != 2.0*value) {
982                        scaleFactor = CoinMin(scaleFactor, 0.25);
983                        if (floor(4.0*value + 0.5) != 4.0*value) {
984                            scaleFactor = 0.0;
985                        }
986                    }
987                }
988            }
989        }
990    }
991    bool possibleMultiple = continuousMultiplier != 0.0 && scaleFactor != 0.0 ;
992    if (possibleMultiple) {
993        for (iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
994            if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
995                maximumCost = CoinMax(maximumCost, fabs(objective[iColumn])) ;
996            }
997        }
998    }
999    setIntParam(CbcModel::CbcFathomDiscipline, possibleMultiple) ;
1000    /*
1001      If a nontrivial increment is possible, try and figure it out. We're looking
1002      for gcd(c<j>) for all c<j> that are coefficients of unfixed integer
1003      variables. Since the c<j> might not be integers, try and inflate them
1004      sufficiently that they look like integers (and we'll deflate the gcd
1005      later).
1006
1007      2520.0 is used as it is a nice multiple of 2,3,5,7
1008    */
1009    if (possibleMultiple && maximumCost) {
1010        int increment = 0 ;
1011        double multiplier = 2520.0 ;
1012        while (10.0*multiplier*maximumCost < 1.0e8)
1013            multiplier *= 10.0 ;
1014        int bigIntegers = 0; // Count of large costs which are integer
1015        for (iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
1016            if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
1017                double objValue = fabs(objective[iColumn]);
1018                if (!isInteger(iColumn)) {
1019                    if (!coeffMultiplier)
1020                        objValue *= continuousMultiplier;
1021                    else
1022                        objValue *= coeffMultiplier[iColumn];
1023                }
1024                if (objValue) {
1025                    double value = objValue * multiplier ;
1026                    if (value < 2.1e9) {
1027                        int nearest = static_cast<int> (floor(value + 0.5)) ;
1028                        if (fabs(value - floor(value + 0.5)) > 1.0e-8) {
1029                            increment = 0 ;
1030                            break ;
1031                        } else if (!increment) {
1032                            increment = nearest ;
1033                        } else {
1034                            increment = gcd(increment, nearest) ;
1035                        }
1036                    } else {
1037                        // large value - may still be multiple of 1.0
1038                        if (fabs(objValue - floor(objValue + 0.5)) > 1.0e-8) {
1039                            increment = 0;
1040                            break;
1041                        } else {
1042                            bigIntegers++;
1043                        }
1044                    }
1045                }
1046            }
1047        }
1048        delete [] coeffMultiplier;
1049        /*
1050          If the increment beats the current value for objective change, install it.
1051        */
1052        if (increment) {
1053            double value = increment ;
1054            double cutoff = getDblParam(CbcModel::CbcCutoffIncrement) ;
1055            if (bigIntegers) {
1056                // allow for 1.0
1057                increment = gcd(increment, static_cast<int> (multiplier));
1058                value = increment;
1059            }
1060            value /= multiplier ;
1061            value *= scaleFactor;
1062            //trueIncrement=CoinMax(cutoff,value);;
1063            if (value*0.999 > cutoff) {
1064                messageHandler()->message(CBC_INTEGERINCREMENT,
1065                                          messages())
1066                << value << CoinMessageEol ;
1067                setDblParam(CbcModel::CbcCutoffIncrement, CoinMax(value*0.999,value-1.0e-4)) ;
1068            }
1069        }
1070    }
1071
1072    return ;
1073}
1074
1075/*
1076saveModel called (carved out of) BranchandBound
1077*/
1078void CbcModel::saveModel(OsiSolverInterface * saveSolver, double * checkCutoffForRestart, bool * feasible)
1079{
1080    if (saveSolver && (specialOptions_&32768) != 0) {
1081        // See if worth trying reduction
1082        *checkCutoffForRestart = getCutoff();
1083        bool tryNewSearch = solverCharacteristics_->reducedCostsAccurate() &&
1084                            (*checkCutoffForRestart < 1.0e20);
1085        int numberColumns = getNumCols();
1086        if (tryNewSearch) {
1087#ifdef CLP_INVESTIGATE
1088            printf("after %d nodes, cutoff %g - looking\n",
1089                   numberNodes_, getCutoff());
1090#endif
1091            saveSolver->resolve();
1092            double direction = saveSolver->getObjSense() ;
1093            double gap = *checkCutoffForRestart - saveSolver->getObjValue() * direction ;
1094            double tolerance;
1095            saveSolver->getDblParam(OsiDualTolerance, tolerance) ;
1096            if (gap <= 0.0)
1097                gap = tolerance;
1098            gap += 100.0 * tolerance;
1099            double integerTolerance = getDblParam(CbcIntegerTolerance) ;
1100
1101            const double *lower = saveSolver->getColLower() ;
1102            const double *upper = saveSolver->getColUpper() ;
1103            const double *solution = saveSolver->getColSolution() ;
1104            const double *reducedCost = saveSolver->getReducedCost() ;
1105
1106            int numberFixed = 0 ;
1107            int numberFixed2 = 0;
1108            for (int i = 0 ; i < numberIntegers_ ; i++) {
1109                int iColumn = integerVariable_[i] ;
1110                double djValue = direction * reducedCost[iColumn] ;
1111                if (upper[iColumn] - lower[iColumn] > integerTolerance) {
1112                    if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) {
1113                        saveSolver->setColUpper(iColumn, lower[iColumn]) ;
1114                        numberFixed++ ;
1115                    } else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) {
1116                        saveSolver->setColLower(iColumn, upper[iColumn]) ;
1117                        numberFixed++ ;
1118                    }
1119                } else {
1120                    numberFixed2++;
1121                }
1122            }
1123#ifdef COIN_DEVELOP
1124            /*
1125              We're debugging. (specialOptions 1)
1126            */
1127            if ((specialOptions_&1) != 0) {
1128                const OsiRowCutDebugger *debugger = saveSolver->getRowCutDebugger() ;
1129                if (debugger) {
1130                    printf("Contains optimal\n") ;
1131                    OsiSolverInterface * temp = saveSolver->clone();
1132                    const double * solution = debugger->optimalSolution();
1133                    const double *lower = temp->getColLower() ;
1134                    const double *upper = temp->getColUpper() ;
1135                    int n = temp->getNumCols();
1136                    for (int i = 0; i < n; i++) {
1137                        if (temp->isInteger(i)) {
1138                            double value = floor(solution[i] + 0.5);
1139                            assert (value >= lower[i] && value <= upper[i]);
1140                            temp->setColLower(i, value);
1141                            temp->setColUpper(i, value);
1142                        }
1143                    }
1144                    temp->writeMps("reduced_fix");
1145                    delete temp;
1146                    saveSolver->writeMps("reduced");
1147                } else {
1148                    abort();
1149                }
1150            }
1151            printf("Restart could fix %d integers (%d already fixed)\n",
1152                   numberFixed + numberFixed2, numberFixed2);
1153#endif
1154            numberFixed += numberFixed2;
1155            if (numberFixed*20 < numberColumns)
1156                tryNewSearch = false;
1157        }
1158        if (tryNewSearch) {
1159            // back to solver without cuts?
1160            OsiSolverInterface * solver2 = continuousSolver_->clone();
1161            const double *lower = saveSolver->getColLower() ;
1162            const double *upper = saveSolver->getColUpper() ;
1163            for (int i = 0 ; i < numberIntegers_ ; i++) {
1164                int iColumn = integerVariable_[i] ;
1165                solver2->setColLower(iColumn, lower[iColumn]);
1166                solver2->setColUpper(iColumn, upper[iColumn]);
1167            }
1168            // swap
1169            delete saveSolver;
1170            saveSolver = solver2;
1171            double * newSolution = new double[numberColumns];
1172            double objectiveValue = *checkCutoffForRestart;
1173            CbcSerendipity heuristic(*this);
1174            if (bestSolution_)
1175                heuristic.setInputSolution(bestSolution_, bestObjective_);
1176            heuristic.setFractionSmall(0.9);
1177            heuristic.setFeasibilityPumpOptions(1008013);
1178            // Use numberNodes to say how many are original rows
1179            heuristic.setNumberNodes(continuousSolver_->getNumRows());
1180#ifdef COIN_DEVELOP
1181            if (continuousSolver_->getNumRows() <
1182                    saveSolver->getNumRows())
1183                printf("%d rows added ZZZZZ\n",
1184                       solver_->getNumRows() - continuousSolver_->getNumRows());
1185#endif
1186            int returnCode = heuristic.smallBranchAndBound(saveSolver,
1187                             -1, newSolution,
1188                             objectiveValue,
1189                             *checkCutoffForRestart, "Reduce");
1190            if (returnCode < 0) {
1191#ifdef COIN_DEVELOP
1192                printf("Restart - not small enough to do search after fixing\n");
1193#endif
1194                delete [] newSolution;
1195            } else {
1196                if ((returnCode&1) != 0) {
1197                    // increment number of solutions so other heuristics can test
1198                    numberSolutions_++;
1199                    numberHeuristicSolutions_++;
1200                    lastHeuristic_ = NULL;
1201                    setBestSolution(CBC_ROUNDING, objectiveValue, newSolution) ;
1202                }
1203                delete [] newSolution;
1204                *feasible = false; // stop search
1205            }
1206#if 0 // probably not needed def CBC_THREAD
1207            if (master_) {
1208                lockThread();
1209                if (parallelMode() > 0) {
1210                    while (master_->waitForThreadsInTree(0)) {
1211                        lockThread();
1212                        double dummyBest;
1213                        tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
1214                        //unlockThread();
1215                    }
1216                }
1217                master_->waitForThreadsInTree(2);
1218                delete master_;
1219                master_ = NULL;
1220                masterThread_ = NULL;
1221            }
1222#endif
1223        }
1224    }
1225}
1226/*
1227Adds integers, called from BranchandBound()
1228*/
1229void CbcModel::AddIntegers()
1230{
1231    int numberColumns = continuousSolver_->getNumCols();
1232    int numberRows = continuousSolver_->getNumRows();
1233    int * del = new int [CoinMax(numberColumns, numberRows)];
1234    int * original = new int [numberColumns];
1235    char * possibleRow = new char [numberRows];
1236    {
1237        const CoinPackedMatrix * rowCopy = continuousSolver_->getMatrixByRow();
1238        const int * column = rowCopy->getIndices();
1239        const int * rowLength = rowCopy->getVectorLengths();
1240        const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
1241        const double * rowLower = continuousSolver_->getRowLower();
1242        const double * rowUpper = continuousSolver_->getRowUpper();
1243        const double * element = rowCopy->getElements();
1244        for (int i = 0; i < numberRows; i++) {
1245            int nLeft = 0;
1246            bool possible = false;
1247            if (rowLower[i] < -1.0e20) {
1248                double value = rowUpper[i];
1249                if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1250                    possible = true;
1251            } else if (rowUpper[i] > 1.0e20) {
1252                double value = rowLower[i];
1253                if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1254                    possible = true;
1255            } else {
1256                double value = rowUpper[i];
1257                if (rowLower[i] == rowUpper[i] &&
1258                        fabs(value - floor(value + 0.5)) < 1.0e-8)
1259                    possible = true;
1260            }
1261            double allSame = (possible) ? 0.0 : -1.0;
1262            for (CoinBigIndex j = rowStart[i];
1263                    j < rowStart[i] + rowLength[i]; j++) {
1264                int iColumn = column[j];
1265                if (continuousSolver_->isInteger(iColumn)) {
1266                    if (fabs(element[j]) != 1.0)
1267                        possible = false;
1268                } else {
1269                    nLeft++;
1270                    if (!allSame) {
1271                      allSame = fabs(element[j]);
1272                    } else if (allSame>0.0) {
1273                      if (allSame!=fabs(element[j]))
1274                        allSame = -1.0;
1275                    }
1276                }
1277            }
1278            if (nLeft == rowLength[i] && allSame > 0.0)
1279                possibleRow[i] = 2;
1280            else if (possible || !nLeft)
1281                possibleRow[i] = 1;
1282            else
1283                possibleRow[i] = 0;
1284        }
1285    }
1286    int nDel = 0;
1287    for (int i = 0; i < numberColumns; i++) {
1288        original[i] = i;
1289        if (continuousSolver_->isInteger(i))
1290            del[nDel++] = i;
1291    }
1292    int nExtra = 0;
1293    OsiSolverInterface * copy1 = continuousSolver_->clone();
1294    int nPass = 0;
1295    while (nDel && nPass < 10) {
1296        nPass++;
1297        OsiSolverInterface * copy2 = copy1->clone();
1298        int nLeft = 0;
1299        for (int i = 0; i < nDel; i++)
1300            original[del[i]] = -1;
1301        for (int i = 0; i < numberColumns; i++) {
1302            int kOrig = original[i];
1303            if (kOrig >= 0)
1304                original[nLeft++] = kOrig;
1305        }
1306        assert (nLeft == numberColumns - nDel);
1307        copy2->deleteCols(nDel, del);
1308        numberColumns = copy2->getNumCols();
1309        const CoinPackedMatrix * rowCopy = copy2->getMatrixByRow();
1310        numberRows = rowCopy->getNumRows();
1311        const int * column = rowCopy->getIndices();
1312        const int * rowLength = rowCopy->getVectorLengths();
1313        const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
1314        const double * rowLower = copy2->getRowLower();
1315        const double * rowUpper = copy2->getRowUpper();
1316        const double * element = rowCopy->getElements();
1317        const CoinPackedMatrix * columnCopy = copy2->getMatrixByCol();
1318        const int * columnLength = columnCopy->getVectorLengths();
1319        nDel = 0;
1320        // Could do gcd stuff on ones with costs
1321        for (int i = 0; i < numberRows; i++) {
1322            if (!rowLength[i]) {
1323                del[nDel++] = i;
1324            } else if (possibleRow[i]) {
1325                if (rowLength[i] == 1) {
1326                    int k = rowStart[i];
1327                    int iColumn = column[k];
1328                    if (!copy2->isInteger(iColumn)) {
1329                        double mult = 1.0 / fabs(element[k]);
1330                        if (rowLower[i] < -1.0e20) {
1331                            // treat rhs as multiple of 1 unless elements all same
1332                            double value = ((possibleRow[i]==2) ? rowUpper[i] : 1.0) * mult;
1333                            if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
1334                                del[nDel++] = i;
1335                                if (columnLength[iColumn] == 1) {
1336                                    copy2->setInteger(iColumn);
1337                                    int kOrig = original[iColumn];
1338                                    setOptionalInteger(kOrig);
1339                                }
1340                            }
1341                        } else if (rowUpper[i] > 1.0e20) {
1342                            // treat rhs as multiple of 1 unless elements all same
1343                            double value = ((possibleRow[i]==2) ? rowLower[i] : 1.0) * mult;
1344                            if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
1345                                del[nDel++] = i;
1346                                if (columnLength[iColumn] == 1) {
1347                                    copy2->setInteger(iColumn);
1348                                    int kOrig = original[iColumn];
1349                                    setOptionalInteger(kOrig);
1350                                }
1351                            }
1352                        } else {
1353                            // treat rhs as multiple of 1 unless elements all same
1354                            double value = ((possibleRow[i]==2) ? rowUpper[i] : 1.0) * mult;
1355                            if (rowLower[i] == rowUpper[i] &&
1356                                    fabs(value - floor(value + 0.5)) < 1.0e-8) {
1357                                del[nDel++] = i;
1358                                copy2->setInteger(iColumn);
1359                                int kOrig = original[iColumn];
1360                                setOptionalInteger(kOrig);
1361                            }
1362                        }
1363                    }
1364                } else {
1365                    // only if all singletons
1366                    bool possible = false;
1367                    if (rowLower[i] < -1.0e20) {
1368                        double value = rowUpper[i];
1369                        if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1370                            possible = true;
1371                    } else if (rowUpper[i] > 1.0e20) {
1372                        double value = rowLower[i];
1373                        if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1374                            possible = true;
1375                    } else {
1376                        double value = rowUpper[i];
1377                        if (rowLower[i] == rowUpper[i] &&
1378                                fabs(value - floor(value + 0.5)) < 1.0e-8)
1379                            possible = true;
1380                    }
1381                    if (possible) {
1382                        for (CoinBigIndex j = rowStart[i];
1383                                j < rowStart[i] + rowLength[i]; j++) {
1384                            int iColumn = column[j];
1385                            if (columnLength[iColumn] != 1 || fabs(element[j]) != 1.0) {
1386                                possible = false;
1387                                break;
1388                            }
1389                        }
1390                        if (possible) {
1391                            for (CoinBigIndex j = rowStart[i];
1392                                    j < rowStart[i] + rowLength[i]; j++) {
1393                                int iColumn = column[j];
1394                                if (!copy2->isInteger(iColumn)) {
1395                                    copy2->setInteger(iColumn);
1396                                    int kOrig = original[iColumn];
1397                                    setOptionalInteger(kOrig);
1398                                }
1399                            }
1400                            del[nDel++] = i;
1401                        }
1402                    }
1403                }
1404            }
1405        }
1406        if (nDel) {
1407            copy2->deleteRows(nDel, del);
1408            // pack down possible
1409            int n=0;
1410            for (int i=0;i<nDel;i++)
1411              possibleRow[del[i]]=-1;
1412            for (int i=0;i<numberRows;i++) {
1413              if (possibleRow[i]>=0) 
1414                possibleRow[n++]=possibleRow[i];
1415            }
1416        }
1417        if (nDel != numberRows) {
1418            nDel = 0;
1419            for (int i = 0; i < numberColumns; i++) {
1420                if (copy2->isInteger(i)) {
1421                    del[nDel++] = i;
1422                    nExtra++;
1423                }
1424            }
1425        } else {
1426            nDel = 0;
1427        }
1428        delete copy1;
1429        copy1 = copy2->clone();
1430        delete copy2;
1431    }
1432    // See if what's left is a network
1433    bool couldBeNetwork = false;
1434    if (copy1->getNumRows() && copy1->getNumCols()) {
1435#ifdef COIN_HAS_CLP
1436        OsiClpSolverInterface * clpSolver
1437        = dynamic_cast<OsiClpSolverInterface *> (copy1);
1438        if (false && clpSolver) {
1439            numberRows = clpSolver->getNumRows();
1440            char * rotate = new char[numberRows];
1441            int n = clpSolver->getModelPtr()->findNetwork(rotate, 1.0);
1442            delete [] rotate;
1443#ifdef CLP_INVESTIGATE
1444            printf("INTA network %d rows out of %d\n", n, numberRows);
1445#endif
1446            if (CoinAbs(n) == numberRows) {
1447                couldBeNetwork = true;
1448                for (int i = 0; i < numberRows; i++) {
1449                    if (!possibleRow[i]) {
1450                        couldBeNetwork = false;
1451#ifdef CLP_INVESTIGATE
1452                        printf("but row %d is bad\n", i);
1453#endif
1454                        break;
1455                    }
1456                }
1457            }
1458        } else
1459#endif
1460        {
1461            numberColumns = copy1->getNumCols();
1462            numberRows = copy1->getNumRows();
1463            const double * rowLower = copy1->getRowLower();
1464            const double * rowUpper = copy1->getRowUpper();
1465            couldBeNetwork = true;
1466            for (int i = 0; i < numberRows; i++) {
1467                if (rowLower[i] > -1.0e20 && fabs(rowLower[i] - floor(rowLower[i] + 0.5)) > 1.0e-12) {
1468                    couldBeNetwork = false;
1469                    break;
1470                }
1471                if (rowUpper[i] < 1.0e20 && fabs(rowUpper[i] - floor(rowUpper[i] + 0.5)) > 1.0e-12) {
1472                    couldBeNetwork = false;
1473                    break;
1474                }
1475                if (possibleRow[i]==0) {
1476                    couldBeNetwork = false;
1477                    break;
1478                }
1479            }
1480            if (couldBeNetwork) {
1481                const CoinPackedMatrix  * matrixByCol = copy1->getMatrixByCol();
1482                const double * element = matrixByCol->getElements();
1483                //const int * row = matrixByCol->getIndices();
1484                const CoinBigIndex * columnStart = matrixByCol->getVectorStarts();
1485                const int * columnLength = matrixByCol->getVectorLengths();
1486                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1487                    CoinBigIndex start = columnStart[iColumn];
1488                    CoinBigIndex end = start + columnLength[iColumn];
1489                    if (end > start + 2) {
1490                        couldBeNetwork = false;
1491                        break;
1492                    }
1493                    int type = 0;
1494                    for (CoinBigIndex j = start; j < end; j++) {
1495                        double value = element[j];
1496                        if (fabs(value) != 1.0) {
1497                            couldBeNetwork = false;
1498                            break;
1499                        } else if (value == 1.0) {
1500                            if ((type&1) == 0)
1501                                type |= 1;
1502                            else
1503                                type = 7;
1504                        } else if (value == -1.0) {
1505                            if ((type&2) == 0)
1506                                type |= 2;
1507                            else
1508                                type = 7;
1509                        }
1510                    }
1511                    if (type > 3) {
1512                        couldBeNetwork = false;
1513                        break;
1514                    }
1515                }
1516            }
1517        }
1518    }
1519    if (couldBeNetwork) {
1520        for (int i = 0; i < numberColumns; i++)
1521            setOptionalInteger(original[i]);
1522    }
1523    if (nExtra || couldBeNetwork) {
1524        numberColumns = copy1->getNumCols();
1525        numberRows = copy1->getNumRows();
1526        if (!numberColumns || !numberRows) {
1527            int numberColumns = solver_->getNumCols();
1528            for (int i = 0; i < numberColumns; i++)
1529                assert(solver_->isInteger(i));
1530        }
1531#ifdef CLP_INVESTIGATE
1532        if (couldBeNetwork || nExtra)
1533            printf("INTA %d extra integers, %d left%s\n", nExtra,
1534                   numberColumns,
1535                   couldBeNetwork ? ", all network" : "");
1536#endif
1537        findIntegers(true, 2);
1538        convertToDynamic();
1539    }
1540#ifdef CLP_INVESTIGATE
1541    if (!couldBeNetwork && copy1->getNumCols() &&
1542            copy1->getNumRows()) {
1543        printf("INTA %d rows and %d columns remain\n",
1544               copy1->getNumRows(), copy1->getNumCols());
1545        if (copy1->getNumCols() < 200) {
1546            copy1->writeMps("moreint");
1547            printf("INTA Written remainder to moreint.mps.gz %d rows %d cols\n",
1548                   copy1->getNumRows(), copy1->getNumCols());
1549        }
1550    }
1551#endif
1552    delete copy1;
1553    delete [] del;
1554    delete [] original;
1555    delete [] possibleRow;
1556    // double check increment
1557    analyzeObjective();
1558}
1559/**
1560  \todo
1561  Normally, it looks like we enter here from command dispatch in the main
1562  routine, after calling the solver for an initial solution
1563  (CbcModel::initialSolve, which simply calls the solver's initialSolve
1564  routine.) The first thing we do is call resolve. Presumably there are
1565  circumstances where this is nontrivial? There's also a call from
1566  CbcModel::originalModel (tied up with integer presolve), which should be
1567  checked.
1568
1569*/
1570
1571/*
1572  The overall flow can be divided into three stages:
1573    * Prep: Check that the lp relaxation remains feasible at the root. If so,
1574      do all the setup for B&C.
1575    * Process the root node: Generate cuts, apply heuristics, and in general do
1576      the best we can to resolve the problem without B&C.
1577    * Do B&C search until we hit a limit or exhaust the search tree.
1578
1579  Keep in mind that in general there is no node in the search tree that
1580  corresponds to the active subproblem. The active subproblem is represented
1581  by the current state of the model,  of the solver, and of the constraint
1582  system held by the solver.
1583*/
1584#ifdef COIN_HAS_CPX
1585#include "OsiCpxSolverInterface.hpp"
1586#include "cplex.h"
1587#endif
1588void CbcModel::branchAndBound(int doStatistics)
1589
1590{
1591  if (!parentModel_)
1592    /*
1593      Capture a time stamp before we start (unless set).
1594    */
1595    if (!dblParam_[CbcStartSeconds]) {
1596      if (!useElapsedTime())
1597        dblParam_[CbcStartSeconds] = CoinCpuTime();
1598      else
1599        dblParam_[CbcStartSeconds] = CoinGetTimeOfDay();
1600    }
1601    dblParam_[CbcSmallestChange] = COIN_DBL_MAX;
1602    dblParam_[CbcSumChange] = 0.0;
1603    dblParam_[CbcLargestChange] = 0.0;
1604    intParam_[CbcNumberBranches] = 0;
1605    // Force minimization !!!!
1606    bool flipObjective = (solver_->getObjSense()<0.0);
1607    if (flipObjective)
1608      flipModel();
1609    dblParam_[CbcOptimizationDirection] = 1.0; // was solver_->getObjSense();
1610    strongInfo_[0] = 0;
1611    strongInfo_[1] = 0;
1612    strongInfo_[2] = 0;
1613    strongInfo_[3] = 0;
1614    strongInfo_[4] = 0;
1615    strongInfo_[5] = 0;
1616    strongInfo_[6] = 0;
1617    numberStrongIterations_ = 0;
1618    currentNode_ = NULL;
1619    // See if should do cuts old way
1620    if (parallelMode() < 0) {
1621        specialOptions_ |= 4096 + 8192;
1622    } else if (parallelMode() > 0) {
1623        specialOptions_ |= 4096;
1624    }
1625    int saveMoreSpecialOptions = moreSpecialOptions_;
1626    if (dynamic_cast<CbcTreeLocal *> (tree_))
1627        specialOptions_ |= 4096 + 8192;
1628#ifdef COIN_HAS_CLP
1629    {
1630        OsiClpSolverInterface * clpSolver
1631        = dynamic_cast<OsiClpSolverInterface *> (solver_);
1632        if (clpSolver) {
1633            // pass in disaster handler
1634            CbcDisasterHandler handler(this);
1635            clpSolver->passInDisasterHandler(&handler);
1636            // Initialise solvers seed (unless users says not)
1637            if ((specialOptions_&4194304)==0)
1638              clpSolver->getModelPtr()->setRandomSeed(1234567);
1639#ifdef JJF_ZERO
1640            // reduce factorization frequency
1641            int frequency = clpSolver->getModelPtr()->factorizationFrequency();
1642            clpSolver->getModelPtr()->setFactorizationFrequency(CoinMin(frequency, 120));
1643#endif
1644        }
1645    }
1646#endif
1647    // original solver (only set if pre-processing)
1648    OsiSolverInterface * originalSolver = NULL;
1649    int numberOriginalObjects = numberObjects_;
1650    OsiObject ** originalObject = NULL;
1651    // Save whether there were any objects
1652    bool noObjects = (numberObjects_ == 0);
1653    // Set up strategies
1654    /*
1655      See if the user has supplied a strategy object and deal with it if present.
1656      The call to setupOther will set numberStrong_ and numberBeforeTrust_, and
1657      perform integer preprocessing, if requested.
1658
1659      We need to hang on to a pointer to solver_. setupOther will assign a
1660      preprocessed solver to model, but will instruct assignSolver not to trash the
1661      existing one.
1662    */
1663    if (strategy_) {
1664        // May do preprocessing
1665        originalSolver = solver_;
1666        strategy_->setupOther(*this);
1667        if (strategy_->preProcessState()) {
1668            // pre-processing done
1669            if (strategy_->preProcessState() < 0) {
1670                // infeasible (or unbounded)
1671                status_ = 0 ;
1672                if (!solver_->isProvenDualInfeasible()) {
1673                  handler_->message(CBC_INFEAS, messages_) << CoinMessageEol ;
1674                  secondaryStatus_ = 1;
1675                } else {
1676                  handler_->message(CBC_UNBOUNDED, 
1677                                    messages_) << CoinMessageEol ;
1678                  secondaryStatus_ = 7;
1679                }
1680                originalContinuousObjective_ = COIN_DBL_MAX;
1681                if (flipObjective)
1682                  flipModel();
1683                return ;
1684            } else if (numberObjects_ && object_) {
1685                numberOriginalObjects = numberObjects_;
1686                // redo sequence
1687                numberIntegers_ = 0;
1688                int numberColumns = getNumCols();
1689                int nOrig = originalSolver->getNumCols();
1690                CglPreProcess * process = strategy_->process();
1691                assert (process);
1692                const int * originalColumns = process->originalColumns();
1693                // allow for cliques etc
1694                nOrig = CoinMax(nOrig, originalColumns[numberColumns-1] + 1);
1695                // try and redo debugger
1696                OsiRowCutDebugger * debugger = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
1697                if (debugger) {
1698                  if (numberColumns<=debugger->numberColumns())
1699                    debugger->redoSolution(numberColumns, originalColumns);
1700                  else
1701                    debugger=NULL; // no idea how to handle (SOS?)
1702                }
1703                // User-provided solution might have been best. Synchronise.
1704                if (bestSolution_) {
1705                    // need to redo - in case no better found in BAB
1706                    // just get integer part right
1707                    for (int i = 0; i < numberColumns; i++) {
1708                        int jColumn = originalColumns[i];
1709                        bestSolution_[i] = bestSolution_[jColumn];
1710                    }
1711                }
1712                originalObject = object_;
1713                // object number or -1
1714                int * temp = new int[nOrig];
1715                int iColumn;
1716                for (iColumn = 0; iColumn < nOrig; iColumn++)
1717                    temp[iColumn] = -1;
1718                int iObject;
1719                int nNonInt = 0;
1720                for (iObject = 0; iObject < numberOriginalObjects; iObject++) {
1721                    iColumn = originalObject[iObject]->columnNumber();
1722                    if (iColumn < 0) {
1723                        nNonInt++;
1724                    } else {
1725                        temp[iColumn] = iObject;
1726                    }
1727                }
1728                int numberNewIntegers = 0;
1729                int numberOldIntegers = 0;
1730                int numberOldOther = 0;
1731                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1732                    int jColumn = originalColumns[iColumn];
1733                    if (temp[jColumn] >= 0) {
1734                        int iObject = temp[jColumn];
1735                        CbcSimpleInteger * obj =
1736                            dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ;
1737                        if (obj)
1738                            numberOldIntegers++;
1739                        else
1740                            numberOldOther++;
1741                    } else if (isInteger(iColumn)) {
1742                        numberNewIntegers++;
1743                    }
1744                }
1745                /*
1746                  Allocate an array to hold the indices of the integer variables.
1747                  Make a large enough array for all objects
1748                */
1749                numberObjects_ = numberNewIntegers + numberOldIntegers + numberOldOther + nNonInt;
1750                object_ = new OsiObject * [numberObjects_];
1751                delete [] integerVariable_;
1752                integerVariable_ = new int [numberNewIntegers+numberOldIntegers];
1753                /*
1754                  Walk the variables again, filling in the indices and creating objects for
1755                  the integer variables. Initially, the objects hold the index and upper &
1756                  lower bounds.
1757                */
1758                numberIntegers_ = 0;
1759                int n = originalColumns[numberColumns-1] + 1;
1760                int * backward = new int[n];
1761                int i;
1762                for ( i = 0; i < n; i++)
1763                    backward[i] = -1;
1764                for (i = 0; i < numberColumns; i++)
1765                    backward[originalColumns[i]] = i;
1766                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1767                    int jColumn = originalColumns[iColumn];
1768                    if (temp[jColumn] >= 0) {
1769                        int iObject = temp[jColumn];
1770                        CbcSimpleInteger * obj =
1771                            dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ;
1772                        if (obj) {
1773                            object_[numberIntegers_] = originalObject[iObject]->clone();
1774                            // redo ids etc
1775                            //object_[numberIntegers_]->resetSequenceEtc(numberColumns,originalColumns);
1776                            object_[numberIntegers_]->resetSequenceEtc(numberColumns, backward);
1777                            integerVariable_[numberIntegers_++] = iColumn;
1778                        }
1779                    } else if (isInteger(iColumn)) {
1780                        object_[numberIntegers_] =
1781                            new CbcSimpleInteger(this, iColumn);
1782                        integerVariable_[numberIntegers_++] = iColumn;
1783                    }
1784                }
1785                delete [] backward;
1786                numberObjects_ = numberIntegers_;
1787                // Now append other column stuff
1788                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1789                    int jColumn = originalColumns[iColumn];
1790                    if (temp[jColumn] >= 0) {
1791                        int iObject = temp[jColumn];
1792                        CbcSimpleInteger * obj =
1793                            dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ;
1794                        if (!obj) {
1795                            object_[numberObjects_] = originalObject[iObject]->clone();
1796                            // redo ids etc
1797                            CbcObject * obj =
1798                                dynamic_cast <CbcObject *>(object_[numberObjects_]) ;
1799                            assert (obj);
1800                            obj->redoSequenceEtc(this, numberColumns, originalColumns);
1801                            numberObjects_++;
1802                        }
1803                    }
1804                }
1805                // now append non column stuff
1806                for (iObject = 0; iObject < numberOriginalObjects; iObject++) {
1807                    iColumn = originalObject[iObject]->columnNumber();
1808                    if (iColumn < 0) {
1809                        // already has column numbers changed
1810                        object_[numberObjects_] = originalObject[iObject]->clone();
1811#ifdef JJF_ZERO
1812                        // redo ids etc
1813                        CbcObject * obj =
1814                            dynamic_cast <CbcObject *>(object_[numberObjects_]) ;
1815                        assert (obj);
1816                        obj->redoSequenceEtc(this, numberColumns, originalColumns);
1817#endif
1818                        numberObjects_++;
1819                    }
1820                }
1821                delete [] temp;
1822                if (!numberObjects_)
1823                    handler_->message(CBC_NOINT, messages_) << CoinMessageEol ;
1824            } else {
1825                int numberColumns = getNumCols();
1826                CglPreProcess * process = strategy_->process();
1827                assert (process);
1828                const int * originalColumns = process->originalColumns();
1829                // try and redo debugger
1830                OsiRowCutDebugger * debugger = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
1831                if (debugger)
1832                    debugger->redoSolution(numberColumns, originalColumns);
1833            }
1834        } else {
1835            //no preprocessing
1836            originalSolver = NULL;
1837        }
1838        strategy_->setupCutGenerators(*this);
1839        strategy_->setupHeuristics(*this);
1840        // Set strategy print level to models
1841        strategy_->setupPrinting(*this, handler_->logLevel());
1842    }
1843    eventHappened_ = false;
1844    CbcEventHandler *eventHandler = getEventHandler() ;
1845    if (eventHandler)
1846        eventHandler->setModel(this);
1847#define CLIQUE_ANALYSIS
1848#ifdef CLIQUE_ANALYSIS
1849    // set up for probing
1850    // If we're doing clever stuff with cliques, additional info here.
1851    if (!parentModel_)
1852        probingInfo_ = new CglTreeProbingInfo(solver_);
1853    else
1854        probingInfo_ = NULL;
1855#else
1856    probingInfo_ = NULL;
1857#endif
1858
1859    // Try for dominated columns
1860    if ((specialOptions_&64) != 0) {
1861        CglDuplicateRow dupcuts(solver_);
1862        dupcuts.setMode(2);
1863        CglStored * storedCuts = dupcuts.outDuplicates(solver_);
1864        if (storedCuts) {
1865          COIN_DETAIL_PRINT(printf("adding dup cuts\n"));
1866            addCutGenerator(storedCuts, 1, "StoredCuts from dominated",
1867                            true, false, false, -200);
1868        }
1869    }
1870    if (!nodeCompare_)
1871        nodeCompare_ = new CbcCompareDefault();;
1872    // See if hot start wanted
1873    CbcCompareBase * saveCompare = NULL;
1874    // User supplied hotstart. Adapt for preprocessing.
1875    if (hotstartSolution_) {
1876        if (strategy_ && strategy_->preProcessState() > 0) {
1877            CglPreProcess * process = strategy_->process();
1878            assert (process);
1879            int n = solver_->getNumCols();
1880            const int * originalColumns = process->originalColumns();
1881            // columns should be in order ... but
1882            double * tempS = new double[n];
1883            for (int i = 0; i < n; i++) {
1884                int iColumn = originalColumns[i];
1885                tempS[i] = hotstartSolution_[iColumn];
1886            }
1887            delete [] hotstartSolution_;
1888            hotstartSolution_ = tempS;
1889            if (hotstartPriorities_) {
1890                int * tempP = new int [n];
1891                for (int i = 0; i < n; i++) {
1892                    int iColumn = originalColumns[i];
1893                    tempP[i] = hotstartPriorities_[iColumn];
1894                }
1895                delete [] hotstartPriorities_;
1896                hotstartPriorities_ = tempP;
1897            }
1898        }
1899        saveCompare = nodeCompare_;
1900        // depth first
1901        nodeCompare_ = new CbcCompareDepth();
1902    }
1903    if (!problemFeasibility_)
1904        problemFeasibility_ = new CbcFeasibilityBase();
1905# ifdef CBC_DEBUG
1906    std::string problemName ;
1907    solver_->getStrParam(OsiProbName, problemName) ;
1908    printf("Problem name - %s\n", problemName.c_str()) ;
1909    solver_->setHintParam(OsiDoReducePrint, false, OsiHintDo, 0) ;
1910# endif
1911    /*
1912      Assume we're done, and see if we're proven wrong.
1913    */
1914    status_ = 0 ;
1915    secondaryStatus_ = 0;
1916    phase_ = 0;
1917    /*
1918      Scan the variables, noting the integer variables. Create an
1919      CbcSimpleInteger object for each integer variable.
1920    */
1921    findIntegers(false) ;
1922    // Say not dynamic pseudo costs
1923    ownership_ &= ~0x40000000;
1924    // If dynamic pseudo costs then do
1925    if (numberBeforeTrust_)
1926        convertToDynamic();
1927    // Set up char array to say if integer (speed)
1928    delete [] integerInfo_;
1929    {
1930        int n = solver_->getNumCols();
1931        integerInfo_ = new char [n];
1932        for (int i = 0; i < n; i++) {
1933            if (solver_->isInteger(i))
1934                integerInfo_[i] = 1;
1935            else
1936                integerInfo_[i] = 0;
1937        }
1938    }
1939    if (preferredWay_) {
1940        // set all unset ones
1941        for (int iObject = 0 ; iObject < numberObjects_ ; iObject++) {
1942            CbcObject * obj =
1943                dynamic_cast <CbcObject *>(object_[iObject]) ;
1944            if (obj && !obj->preferredWay())
1945                obj->setPreferredWay(preferredWay_);
1946        }
1947    }
1948    /*
1949      Ensure that objects on the lists of OsiObjects, heuristics, and cut
1950      generators attached to this model all refer to this model.
1951    */
1952    synchronizeModel() ;
1953    if (!solverCharacteristics_) {
1954        OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
1955        if (solverCharacteristics) {
1956            solverCharacteristics_ = solverCharacteristics;
1957        } else {
1958            // replace in solver
1959            OsiBabSolver defaultC;
1960            solver_->setAuxiliaryInfo(&defaultC);
1961            solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
1962        }
1963    }
1964
1965    solverCharacteristics_->setSolver(solver_);
1966    // Set so we can tell we are in initial phase in resolve
1967    continuousObjective_ = -COIN_DBL_MAX ;
1968    /*
1969      Solve the relaxation.
1970
1971      Apparently there are circumstances where this will be non-trivial --- i.e.,
1972      we've done something since initialSolve that's trashed the solution to the
1973      continuous relaxation.
1974    */
1975    /* Tell solver we are in Branch and Cut
1976       Could use last parameter for subtle differences */
1977    solver_->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
1978#ifdef COIN_HAS_CLP
1979    {
1980        OsiClpSolverInterface * clpSolver
1981        = dynamic_cast<OsiClpSolverInterface *> (solver_);
1982        if (clpSolver) {
1983            ClpSimplex * clpSimplex = clpSolver->getModelPtr();
1984            if ((specialOptions_&32) == 0) {
1985                // take off names
1986                clpSimplex->dropNames();
1987            }
1988            // no crunch if mostly continuous
1989            if ((clpSolver->specialOptions()&(1 + 8)) != (1 + 8)) {
1990                int numberColumns = solver_->getNumCols();
1991                if (numberColumns > 1000 && numberIntegers_*4 < numberColumns)
1992                    clpSolver->setSpecialOptions(clpSolver->specialOptions()&(~1));
1993            }
1994            //#define NO_CRUNCH
1995#ifdef NO_CRUNCH
1996            printf("TEMP switching off crunch\n");
1997            int iOpt = clpSolver->specialOptions();
1998            iOpt &= ~1;
1999            iOpt |= 65536;
2000            clpSolver->setSpecialOptions(iOpt);
2001#endif
2002        }
2003    }
2004#endif
2005    bool feasible;
2006    numberSolves_ = 0 ;
2007    // If NLP then we assume already solved outside branchAndbound
2008    if (!solverCharacteristics_->solverType() || solverCharacteristics_->solverType() == 4) {
2009        feasible = resolve(NULL, 0) != 0 ;
2010    } else {
2011        // pick up given status
2012        feasible = (solver_->isProvenOptimal() &&
2013                    !solver_->isDualObjectiveLimitReached()) ;
2014    }
2015    if (problemFeasibility_->feasible(this, 0) < 0) {
2016        feasible = false; // pretend infeasible
2017    }
2018    numberSavedSolutions_ = 0;
2019    int saveNumberStrong = numberStrong_;
2020    int saveNumberBeforeTrust = numberBeforeTrust_;
2021    /*
2022      If the linear relaxation of the root is infeasible, bail out now. Otherwise,
2023      continue with processing the root node.
2024    */
2025    if (!feasible) {
2026        status_ = 0 ;
2027        if (!solver_->isProvenDualInfeasible()) {
2028            handler_->message(CBC_INFEAS, messages_) << CoinMessageEol ;
2029            secondaryStatus_ = 1;
2030        } else {
2031            handler_->message(CBC_UNBOUNDED, messages_) << CoinMessageEol ;
2032            secondaryStatus_ = 7;
2033        }
2034        originalContinuousObjective_ = COIN_DBL_MAX;
2035        solverCharacteristics_ = NULL;
2036        if (flipObjective)
2037          flipModel();
2038        return ;
2039    } else if (!numberObjects_) { 
2040        // nothing to do
2041        // Undo preprocessing performed during BaB.
2042        if (strategy_ && strategy_->preProcessState() > 0) {
2043          // undo preprocessing
2044          CglPreProcess * process = strategy_->process();
2045          assert (process);
2046          int n = originalSolver->getNumCols();
2047          if (bestSolution_) {
2048            delete [] bestSolution_;
2049            bestSolution_ = new double [n];
2050            process->postProcess(*solver_);
2051          }
2052          strategy_->deletePreProcess();
2053          // Solution now back in originalSolver
2054          delete solver_;
2055          solver_ = originalSolver;
2056          if (bestSolution_) {
2057            bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
2058            memcpy(bestSolution_, solver_->getColSolution(), n*sizeof(double));
2059          }
2060          // put back original objects if there were any
2061          if (originalObject) {
2062            int iColumn;
2063            assert (ownObjects_);
2064            for (iColumn = 0; iColumn < numberObjects_; iColumn++)
2065              delete object_[iColumn];
2066            delete [] object_;
2067            numberObjects_ = numberOriginalObjects;
2068            object_ = originalObject;
2069            delete [] integerVariable_;
2070            numberIntegers_ = 0;
2071            for (iColumn = 0; iColumn < n; iColumn++) {
2072              if (solver_->isInteger(iColumn))
2073                numberIntegers_++;
2074            }
2075            integerVariable_ = new int[numberIntegers_];
2076            numberIntegers_ = 0;
2077            for (iColumn = 0; iColumn < n; iColumn++) {
2078            if (solver_->isInteger(iColumn))
2079              integerVariable_[numberIntegers_++] = iColumn;
2080            }
2081          }
2082        }
2083        if (flipObjective)
2084          flipModel();
2085        solverCharacteristics_ = NULL;
2086        bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
2087        int numberColumns = solver_->getNumCols();
2088        delete [] bestSolution_;
2089        bestSolution_ = new double[numberColumns];
2090        CoinCopyN(solver_->getColSolution(), numberColumns, bestSolution_);
2091        return ;
2092    }
2093    /*
2094      See if we're using the Osi side of the branching hierarchy. If so, either
2095      convert existing CbcObjects to OsiObjects, or generate them fresh. In the
2096      first case, CbcModel owns the objects on the object_ list. In the second
2097      case, the solver holds the objects and object_ simply points to the
2098      solver's list.
2099
2100      080417 The conversion code here (the block protected by `if (obj)') cannot
2101      possibly be correct. On the Osi side, descent is OsiObject -> OsiObject2 ->
2102      all other Osi object classes. On the Cbc side, it's OsiObject -> CbcObject
2103      -> all other Cbc object classes. It's structurally impossible for any Osi
2104      object to descend from CbcObject. The only thing I can see is that this is
2105      really dead code, and object detection is now handled from the Osi side.
2106    */
2107    // Convert to Osi if wanted
2108    //OsiBranchingInformation * persistentInfo = NULL;
2109    if (branchingMethod_ && branchingMethod_->chooseMethod()) {
2110        //persistentInfo = new OsiBranchingInformation(solver_);
2111        if (numberOriginalObjects) {
2112            for (int iObject = 0 ; iObject < numberObjects_ ; iObject++) {
2113                CbcObject * obj =
2114                    dynamic_cast <CbcObject *>(object_[iObject]) ;
2115                if (obj) {
2116                    CbcSimpleInteger * obj2 =
2117                        dynamic_cast <CbcSimpleInteger *>(obj) ;
2118                    if (obj2) {
2119                        // back to Osi land
2120                        object_[iObject] = obj2->osiObject();
2121                        delete obj;
2122                    } else {
2123                        OsiSimpleInteger * obj3 =
2124                            dynamic_cast <OsiSimpleInteger *>(obj) ;
2125                        if (!obj3) {
2126                            OsiSOS * obj4 =
2127                                dynamic_cast <OsiSOS *>(obj) ;
2128                            if (!obj4) {
2129                                CbcSOS * obj5 =
2130                                    dynamic_cast <CbcSOS *>(obj) ;
2131                                if (obj5) {
2132                                    // back to Osi land
2133                                    object_[iObject] = obj5->osiObject(solver_);
2134                                } else {
2135                                    printf("Code up CbcObject type in Osi land\n");
2136                                    abort();
2137                                }
2138                            }
2139                        }
2140                    }
2141                }
2142            }
2143            // and add to solver
2144            //if (!solver_->numberObjects()) {
2145            solver_->addObjects(numberObjects_, object_);
2146            //} else {
2147            //if (solver_->numberObjects()!=numberOriginalObjects) {
2148            //printf("should have trapped that solver has objects before\n");
2149            //abort();
2150            //}
2151            //}
2152        } else {
2153            /*
2154              As of 080104, findIntegersAndSOS is misleading --- the default OSI
2155              implementation finds only integers.
2156            */
2157            // do from solver
2158            deleteObjects(false);
2159            solver_->findIntegersAndSOS(false);
2160            numberObjects_ = solver_->numberObjects();
2161            object_ = solver_->objects();
2162            ownObjects_ = false;
2163        }
2164        branchingMethod_->chooseMethod()->setSolver(solver_);
2165    }
2166    // take off heuristics if have to (some do not work with SOS, for example)
2167    // object should know what's safe.
2168    {
2169        int numberOdd = 0;
2170        int numberSOS = 0;
2171        for (int i = 0; i < numberObjects_; i++) {
2172            if (!object_[i]->canDoHeuristics())
2173                numberOdd++;
2174            CbcSOS * obj =
2175                dynamic_cast <CbcSOS *>(object_[i]) ;
2176            if (obj)
2177                numberSOS++;
2178        }
2179        if (numberOdd) {
2180            if (numberHeuristics_) {
2181                int k = 0;
2182                for (int i = 0; i < numberHeuristics_; i++) {
2183                    if (!heuristic_[i]->canDealWithOdd())
2184                        delete heuristic_[i];
2185                    else
2186                        heuristic_[k++] = heuristic_[i];
2187                }
2188                if (!k) {
2189                    delete [] heuristic_;
2190                    heuristic_ = NULL;
2191                }
2192                numberHeuristics_ = k;
2193                handler_->message(CBC_HEURISTICS_OFF, messages_) << numberOdd << CoinMessageEol ;
2194            }
2195            // If odd switch off AddIntegers
2196            specialOptions_ &= ~65536;
2197        } else if (numberSOS) {
2198            specialOptions_ |= 128; // say can do SOS in dynamic mode
2199            // switch off fast nodes for now
2200            fastNodeDepth_ = -1;
2201            moreSpecialOptions_ &= ~33554432; // no diving
2202        }
2203        if (numberThreads_ > 0) {
2204            // switch off fast nodes for now
2205            fastNodeDepth_ = -1;
2206        }
2207    }
2208    // Save objective (just so user can access it)
2209    originalContinuousObjective_ = solver_->getObjValue()* solver_->getObjSense();
2210    bestPossibleObjective_ = originalContinuousObjective_;
2211    sumChangeObjective1_ = 0.0;
2212    sumChangeObjective2_ = 0.0;
2213    /*
2214      OsiRowCutDebugger knows an optimal answer for a subset of MIP problems.
2215      Assuming it recognises the problem, when called upon it will check a cut to
2216      see if it cuts off the optimal answer.
2217    */
2218    // If debugger exists set specialOptions_ bit
2219    if (solver_->getRowCutDebuggerAlways()) {
2220        specialOptions_ |= 1;
2221    }
2222
2223# ifdef CBC_DEBUG
2224    if ((specialOptions_&1) == 0)
2225        solver_->activateRowCutDebugger(problemName.c_str()) ;
2226    if (solver_->getRowCutDebuggerAlways())
2227        specialOptions_ |= 1;
2228# endif
2229
2230    /*
2231      Begin setup to process a feasible root node.
2232    */
2233    bestObjective_ = CoinMin(bestObjective_, 1.0e50) ;
2234    if (!bestSolution_) {
2235        numberSolutions_ = 0 ;
2236        numberHeuristicSolutions_ = 0 ;
2237    }
2238    stateOfSearch_ = 0;
2239    // Everything is minimization
2240    {
2241        // needed to sync cutoffs
2242        double value ;
2243        solver_->getDblParam(OsiDualObjectiveLimit, value) ;
2244        dblParam_[CbcCurrentCutoff] = value * solver_->getObjSense();
2245    }
2246    double cutoff = getCutoff() ;
2247    double direction = solver_->getObjSense() ;
2248    dblParam_[CbcOptimizationDirection] = direction;
2249    if (cutoff < 1.0e20 && direction < 0.0)
2250        messageHandler()->message(CBC_CUTOFF_WARNING1,
2251                                  messages())
2252        << cutoff << -cutoff << CoinMessageEol ;
2253    if (cutoff > bestObjective_)
2254        cutoff = bestObjective_ ;
2255    setCutoff(cutoff) ;
2256    /*
2257      We probably already have a current solution, but just in case ...
2258    */
2259    int numberColumns = getNumCols() ;
2260    if (!currentSolution_)
2261        currentSolution_ = new double[numberColumns] ;
2262    testSolution_ = currentSolution_;
2263    /*
2264      Create a copy of the solver, thus capturing the original (root node)
2265      constraint system (aka the continuous system).
2266    */
2267    continuousSolver_ = solver_->clone() ;
2268
2269    // add cutoff as constraint if wanted
2270    if (cutoffRowNumber_==-2) {
2271      if (!parentModel_) {
2272        int numberColumns=solver_->getNumCols();
2273        double * obj = CoinCopyOfArray(solver_->getObjCoefficients(),numberColumns);
2274        int * indices = new int [numberColumns];
2275        int n=0;
2276        for (int i=0;i<numberColumns;i++) {
2277          if (obj[i]) {
2278            indices[n]=i;
2279            obj[n++]=obj[i];
2280          }
2281        }
2282        if (n) {
2283          double cutoff=getCutoff();
2284          // relax a little bit
2285          cutoff += 1.0e-4;
2286          double offset;
2287          solver_->getDblParam(OsiObjOffset, offset);
2288          cutoffRowNumber_ = solver_->getNumRows();
2289          solver_->addRow(n,indices,obj,-COIN_DBL_MAX,CoinMin(cutoff,1.0e25)+offset);
2290        } else {
2291          // no objective!
2292          cutoffRowNumber_ = -1;
2293        }
2294        delete [] indices;
2295        delete [] obj;
2296      } else {
2297        // switch off
2298        cutoffRowNumber_ = -1;
2299      }
2300    }
2301    numberRowsAtContinuous_ = getNumRows() ;
2302    solver_->saveBaseModel();
2303    /*
2304      Check the objective to see if we can deduce a nontrivial increment. If
2305      it's better than the current value for CbcCutoffIncrement, it'll be
2306      installed.
2307    */
2308    if (solverCharacteristics_->reducedCostsAccurate())
2309        analyzeObjective() ;
2310    {
2311        // may be able to change cutoff now
2312        double cutoff = getCutoff();
2313        double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
2314        if (cutoff > bestObjective_ - increment) {
2315            cutoff = bestObjective_ - increment ;
2316            setCutoff(cutoff) ;
2317        }
2318    }
2319#ifdef COIN_HAS_CLP
2320    // Possible save of pivot method
2321    ClpDualRowPivot * savePivotMethod = NULL;
2322    {
2323        // pass tolerance and increment to solver
2324        OsiClpSolverInterface * clpSolver
2325        = dynamic_cast<OsiClpSolverInterface *> (solver_);
2326        if (clpSolver)
2327            clpSolver->setStuff(getIntegerTolerance(), getCutoffIncrement());
2328#ifdef CLP_RESOLVE
2329        if ((moreSpecialOptions_&1048576)!=0&&!parentModel_&&clpSolver) {
2330          resolveClp(clpSolver,0);
2331        }
2332#endif
2333    }
2334#endif
2335    /*
2336      Set up for cut generation. addedCuts_ holds the cuts which are relevant for
2337      the active subproblem. whichGenerator will be used to record the generator
2338      that produced a given cut.
2339    */
2340    maximumWhich_ = 1000 ;
2341    delete [] whichGenerator_;
2342    whichGenerator_ = new int[maximumWhich_] ;
2343    memset(whichGenerator_, 0, maximumWhich_*sizeof(int));
2344    maximumNumberCuts_ = 0 ;
2345    currentNumberCuts_ = 0 ;
2346    delete [] addedCuts_ ;
2347    addedCuts_ = NULL ;
2348    OsiObject ** saveObjects = NULL;
2349    maximumRows_ = numberRowsAtContinuous_;
2350    currentDepth_ = 0;
2351    workingBasis_.resize(maximumRows_, numberColumns);
2352    /*
2353      Set up an empty heap and associated data structures to hold the live set
2354      (problems which require further exploration).
2355    */
2356    CbcCompareDefault * compareActual
2357    = dynamic_cast<CbcCompareDefault *> (nodeCompare_);
2358    if (compareActual) {
2359        compareActual->setBestPossible(direction*solver_->getObjValue());
2360        compareActual->setCutoff(getCutoff());
2361#ifdef JJF_ZERO
2362        if (false && !numberThreads_ && !parentModel_) {
2363            printf("CbcTreeArray ? threads ? parentArray\n");
2364            // Setup new style tree
2365            delete tree_;
2366            tree_ = new CbcTreeArray();
2367        }
2368#endif
2369    }
2370    tree_->setComparison(*nodeCompare_) ;
2371    /*
2372      Used to record the path from a node to the root of the search tree, so that
2373      we can then traverse from the root to the node when restoring a subproblem.
2374    */
2375    maximumDepth_ = 10 ;
2376    delete [] walkback_ ;
2377    walkback_ = new CbcNodeInfo * [maximumDepth_] ;
2378    lastDepth_ = 0;
2379    delete [] lastNodeInfo_ ;
2380    lastNodeInfo_ = new CbcNodeInfo * [maximumDepth_] ;
2381    delete [] lastNumberCuts_ ;
2382    lastNumberCuts_ = new int [maximumDepth_] ;
2383    maximumCuts_ = 100;
2384    lastNumberCuts2_ = 0;
2385    delete [] lastCut_;
2386    lastCut_ = new const OsiRowCut * [maximumCuts_];
2387    /*
2388      Used to generate bound edits for CbcPartialNodeInfo.
2389    */
2390    double * lowerBefore = new double [numberColumns] ;
2391    double * upperBefore = new double [numberColumns] ;
2392    /*
2393    Set up to run heuristics and generate cuts at the root node. The heavy
2394    lifting is hidden inside the calls to doHeuristicsAtRoot and solveWithCuts.
2395
2396    To start, tell cut generators they can be a bit more aggressive at the
2397    root node.
2398
2399    QUESTION: phase_ = 0 is documented as `initial solve', phase = 1 as `solve
2400        with cuts at root'. Is phase_ = 1 the correct indication when
2401        doHeurisiticsAtRoot is called to run heuristics outside of the main
2402        cut / heurisitc / reoptimise loop in solveWithCuts?
2403
2404      Generate cuts at the root node and reoptimise. solveWithCuts does the heavy
2405      lifting. It will iterate a generate/reoptimise loop (including reduced cost
2406      fixing) until no cuts are generated, the change in objective falls off,  or
2407      the limit on the number of rounds of cut generation is exceeded.
2408
2409      At the end of all this, any cuts will be recorded in cuts and also
2410      installed in the solver's constraint system. We'll have reoptimised, and
2411      removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
2412      adjusted accordingly).
2413
2414      Tell cut generators they can be a bit more aggressive at root node
2415
2416      TODO: Why don't we make a copy of the solution after solveWithCuts?
2417      TODO: If numberUnsatisfied == 0, don't we have a solution?
2418    */
2419    phase_ = 1;
2420    int iCutGenerator;
2421    for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
2422        // If parallel switch off global cuts
2423        if (numberThreads_) {
2424            generator_[iCutGenerator]->setGlobalCuts(false);
2425            generator_[iCutGenerator]->setGlobalCutsAtRoot(false);
2426        }
2427        CglCutGenerator * generator = generator_[iCutGenerator]->generator();
2428        generator->setAggressiveness(generator->getAggressiveness() + 100);
2429        if (!generator->canDoGlobalCuts())
2430          generator->setGlobalCuts(false);
2431    }
2432    OsiCuts cuts ;
2433    int anyAction = -1 ;
2434    numberOldActiveCuts_ = 0 ;
2435    numberNewCuts_ = 0 ;
2436    // Array to mark solution
2437    delete [] usedInSolution_;
2438    usedInSolution_ = new int[numberColumns];
2439    CoinZeroN(usedInSolution_, numberColumns);
2440    /*
2441      For printing totals and for CbcNode (numberNodes_)
2442    */
2443    numberIterations_ = 0 ;
2444    numberNodes_ = 0 ;
2445    numberNodes2_ = 0 ;
2446    maximumStatistics_ = 0;
2447    maximumDepthActual_ = 0;
2448    numberDJFixed_ = 0.0;
2449    if (!parentModel_) {
2450        if ((specialOptions_&262144) != 0) {
2451            // create empty stored cuts
2452            //storedRowCuts_ = new CglStored(solver_->getNumCols());
2453        } else if ((specialOptions_&524288) != 0 && storedRowCuts_) {
2454            // tighten and set best solution
2455            // A) tight bounds on integer variables
2456            /*
2457                storedRowCuts_ are coming in from outside, probably for nonlinear.
2458              John was unsure about origin.
2459            */
2460            const double * lower = solver_->getColLower();
2461            const double * upper = solver_->getColUpper();
2462            const double * tightLower = storedRowCuts_->tightLower();
2463            const double * tightUpper = storedRowCuts_->tightUpper();
2464            int nTightened = 0;
2465            for (int i = 0; i < numberIntegers_; i++) {
2466                int iColumn = integerVariable_[i];
2467                if (tightLower[iColumn] > lower[iColumn]) {
2468                    nTightened++;
2469                    solver_->setColLower(iColumn, tightLower[iColumn]);
2470                }
2471                if (tightUpper[iColumn] < upper[iColumn]) {
2472                    nTightened++;
2473                    solver_->setColUpper(iColumn, tightUpper[iColumn]);
2474                }
2475            }
2476            if (nTightened)
2477              COIN_DETAIL_PRINT(printf("%d tightened by alternate cuts\n", nTightened));
2478            if (storedRowCuts_->bestObjective() < bestObjective_) {
2479                // B) best solution
2480                double objValue = storedRowCuts_->bestObjective();
2481                setBestSolution(CBC_SOLUTION, objValue,
2482                                storedRowCuts_->bestSolution()) ;
2483                // Do heuristics
2484                // Allow RINS
2485                for (int i = 0; i < numberHeuristics_; i++) {
2486                    CbcHeuristicRINS * rins
2487                    = dynamic_cast<CbcHeuristicRINS *> (heuristic_[i]);
2488                    if (rins) {
2489                        rins->setLastNode(-100);
2490                    }
2491                }
2492            }
2493        }
2494    }
2495#ifdef SWITCH_VARIABLES
2496    // see if any switching variables
2497    if (numberIntegers_<solver_->getNumCols())
2498      findSwitching();
2499#endif
2500    /*
2501      Run heuristics at the root. This is the only opportunity to run FPump; it
2502      will be removed from the heuristics list by doHeuristicsAtRoot.
2503    */
2504    // See if multiple runs wanted
2505    CbcModel ** rootModels=NULL;
2506    if (!parentModel_&&multipleRootTries_%100) {
2507      double rootTimeCpu=CoinCpuTime();
2508      double startTimeRoot=CoinGetTimeOfDay();
2509      int numberRootThreads=1;
2510      /* undocumented fine tuning
2511         aabbcc where cc is number of tries
2512         bb if nonzero is number of threads
2513         aa if nonzero just do heuristics
2514      */
2515      int numberModels = multipleRootTries_%100;
2516#ifdef CBC_THREAD
2517      numberRootThreads = (multipleRootTries_/100)%100;
2518      if (!numberRootThreads)
2519        numberRootThreads=numberModels;
2520#endif
2521      int otherOptions = (multipleRootTries_/10000)%100;
2522      rootModels = new CbcModel * [numberModels];
2523      unsigned int newSeed = randomSeed_;
2524      if (newSeed==0) {
2525        double time = fabs(CoinGetTimeOfDay());
2526        while (time>=COIN_INT_MAX)
2527          time *= 0.5;
2528        newSeed = static_cast<unsigned int>(time);
2529      } else if (newSeed<0) {
2530        newSeed = 123456789;
2531#ifdef COIN_HAS_CLP
2532        OsiClpSolverInterface * clpSolver
2533          = dynamic_cast<OsiClpSolverInterface *> (solver_);
2534        if (clpSolver) {
2535          newSeed += clpSolver->getModelPtr()->randomNumberGenerator()->getSeed();
2536        }
2537#endif
2538      }
2539      CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis *> (solver_->getEmptyWarmStart());
2540      int numberIterations=0;
2541      for (int i=0;i<numberModels;i++) { 
2542        rootModels[i]=new CbcModel(*this);
2543        rootModels[i]->setNumberThreads(0);
2544        rootModels[i]->setMaximumNodes(otherOptions ? -1 : 0);
2545        rootModels[i]->setRandomSeed(newSeed+10000000*i);
2546        rootModels[i]->randomNumberGenerator()->setSeed(newSeed+50000000*i);
2547        rootModels[i]->setMultipleRootTries(0);
2548        // use seed
2549        rootModels[i]->setSpecialOptions(specialOptions_ |(4194304|8388608));
2550        rootModels[i]->setMoreSpecialOptions(moreSpecialOptions_ & 
2551                                             (~134217728));
2552        rootModels[i]->solver_->setWarmStart(basis);
2553#ifdef COIN_HAS_CLP
2554        OsiClpSolverInterface * clpSolver
2555          = dynamic_cast<OsiClpSolverInterface *> (rootModels[i]->solver_);
2556        if (clpSolver) {
2557          ClpSimplex * simplex = clpSolver->getModelPtr();
2558          if (defaultHandler_)
2559            simplex->setDefaultMessageHandler();
2560          simplex->setRandomSeed(newSeed+20000000*i);
2561          simplex->allSlackBasis();
2562          int logLevel=simplex->logLevel();
2563          if (logLevel==1)
2564            simplex->setLogLevel(0);
2565          if (i!=0) {
2566            double random = simplex->randomNumberGenerator()->randomDouble();
2567            int bias = static_cast<int>(random*(numberIterations/4));
2568            simplex->setMaximumIterations(numberIterations/2+bias);
2569            simplex->primal();
2570            simplex->setMaximumIterations(COIN_INT_MAX);
2571            simplex->dual();
2572          } else {
2573            simplex->primal();
2574            numberIterations=simplex->numberIterations();
2575          }
2576          simplex->setLogLevel(logLevel);
2577          clpSolver->setWarmStart(NULL);
2578        }
2579#endif
2580        for (int j=0;j<numberHeuristics_;j++)
2581          rootModels[i]->heuristic_[j]->setSeed(rootModels[i]->heuristic_[j]->getSeed()+100000000*i);
2582        for (int j=0;j<numberCutGenerators_;j++)
2583          rootModels[i]->generator_[j]->generator()->refreshSolver(rootModels[i]->solver_);
2584      }
2585      delete basis;
2586#ifdef CBC_THREAD
2587      if (numberRootThreads==1) {
2588#endif
2589        for (int iModel=0;iModel<numberModels;iModel++) {
2590          doRootCbcThread(rootModels[iModel]);
2591          // see if solved at root node
2592          if (rootModels[iModel]->getMaximumNodes()) {
2593            feasible=false;
2594            break;
2595          }
2596        }
2597#ifdef CBC_THREAD
2598      } else {
2599        Coin_pthread_t * threadId = new Coin_pthread_t [numberRootThreads];
2600        for (int kModel=0;kModel<numberModels;kModel+=numberRootThreads) {
2601          bool finished=false;
2602          for (int iModel=kModel;iModel<CoinMin(numberModels,kModel+numberRootThreads);iModel++) {
2603            pthread_create(&(threadId[iModel-kModel].thr), NULL, 
2604                           doRootCbcThread,
2605                           rootModels[iModel]);
2606          }
2607          // wait
2608          for (int iModel=kModel;iModel<CoinMin(numberModels,kModel+numberRootThreads);iModel++) {
2609            pthread_join(threadId[iModel-kModel].thr, NULL);
2610          }
2611          // see if solved at root node
2612          for (int iModel=kModel;iModel<CoinMin(numberModels,kModel+numberRootThreads);iModel++) {
2613            if (rootModels[iModel]->getMaximumNodes())
2614              finished=true;
2615          }
2616          if (finished) {
2617            feasible=false;
2618            break;
2619          }
2620        }
2621        delete [] threadId;
2622      }
2623#endif
2624      // sort solutions
2625      int * which = new int [numberModels];
2626      double * value = new double [numberModels];
2627      int numberSolutions=0;
2628      for (int iModel=0;iModel<numberModels;iModel++) {
2629        if (rootModels[iModel]->bestSolution()) {
2630          which[numberSolutions]=iModel;
2631          value[numberSolutions++]=
2632            -rootModels[iModel]->getMinimizationObjValue();
2633        }
2634      }
2635      char general[100];
2636      rootTimeCpu=CoinCpuTime()-rootTimeCpu;
2637      if (numberRootThreads==1)
2638        sprintf(general,"Multiple root solvers took a total of %.2f seconds\n",
2639                rootTimeCpu);
2640      else
2641        sprintf(general,"Multiple root solvers took a total of %.2f seconds (%.2f elapsed)\n",
2642                rootTimeCpu,CoinGetTimeOfDay()-startTimeRoot);
2643      messageHandler()->message(CBC_GENERAL,
2644                                messages())
2645        << general << CoinMessageEol ;
2646      CoinSort_2(value,value+numberSolutions,which);
2647      // to get name
2648      CbcHeuristicRINS dummyHeuristic;
2649      dummyHeuristic.setHeuristicName("Multiple root solvers");
2650      lastHeuristic_=&dummyHeuristic;
2651      for (int i=0;i<numberSolutions;i++) {
2652        double objValue = - value[i];
2653        if (objValue<getCutoff()) {
2654          int iModel=which[i];
2655          setBestSolution(CBC_ROUNDING,objValue,
2656                          rootModels[iModel]->bestSolution());
2657        }
2658      }
2659      lastHeuristic_=NULL;
2660      delete [] which;
2661      delete [] value;
2662    }
2663    // Do heuristics
2664    if (numberObjects_&&!rootModels)
2665        doHeuristicsAtRoot();
2666    if (solverCharacteristics_->solutionAddsCuts()) {
2667      // With some heuristics solver needs a resolve here
2668      solver_->resolve();
2669      if(!isProvenOptimal()){
2670        solver_->initialSolve();
2671      }
2672    }
2673    /*
2674      Grepping through the code, it would appear that this is a command line
2675      debugging hook.  There's no obvious place in the code where this is set to
2676      a negative value.
2677
2678      User hook, says John.
2679    */
2680    if ( intParam_[CbcMaxNumNode] < 0
2681      ||numberSolutions_>=getMaximumSolutions())
2682      eventHappened_ = true; // stop as fast as possible
2683    stoppedOnGap_ = false ;
2684    // See if can stop on gap
2685    bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
2686    if(canStopOnGap()) {
2687        if (bestPossibleObjective_ < getCutoff())
2688            stoppedOnGap_ = true ;
2689        feasible = false;
2690        //eventHappened_=true; // stop as fast as possible
2691    }
2692    /*
2693      Set up for statistics collection, if requested. Standard values are
2694      documented in CbcModel.hpp. The magic number 100 will trigger a dump of
2695      CbcSimpleIntegerDynamicPseudoCost objects (no others). Looks like another
2696      command line debugging hook.
2697    */
2698    statistics_ = NULL;
2699    // Do on switch
2700    if (doStatistics > 0 && doStatistics <= 100) {
2701        maximumStatistics_ = 10000;
2702        statistics_ = new CbcStatistics * [maximumStatistics_];
2703        memset(statistics_, 0, maximumStatistics_*sizeof(CbcStatistics *));
2704    }
2705    // See if we can add integers
2706    if (noObjects && numberIntegers_ < solver_->getNumCols() && (specialOptions_&65536) != 0 && !parentModel_)
2707        AddIntegers();
2708    /*
2709      Do an initial round of cut generation for the root node. Depending on the
2710      type of underlying solver, we may want to do this even if the initial query
2711      to the objects indicates they're satisfied.
2712
2713      solveWithCuts does the heavy lifting. It will iterate a generate/reoptimise
2714      loop (including reduced cost fixing) until no cuts are generated, the
2715      change in objective falls off,  or the limit on the number of rounds of cut
2716      generation is exceeded.
2717
2718      At the end of all this, any cuts will be recorded in cuts and also
2719      installed in the solver's constraint system. We'll have reoptimised, and
2720      removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
2721      adjusted accordingly).
2722    */
2723    int iObject ;
2724    int preferredWay ;
2725    int numberUnsatisfied = 0 ;
2726    delete [] currentSolution_;
2727    currentSolution_ = new double [numberColumns];
2728    testSolution_ = currentSolution_;
2729    memcpy(currentSolution_, solver_->getColSolution(),
2730           numberColumns*sizeof(double)) ;
2731    // point to useful information
2732    OsiBranchingInformation usefulInfo = usefulInformation();
2733
2734    for (iObject = 0 ; iObject < numberObjects_ ; iObject++) {
2735        double infeasibility =
2736            object_[iObject]->infeasibility(&usefulInfo, preferredWay) ;
2737        if (infeasibility ) numberUnsatisfied++ ;
2738    }
2739    // replace solverType
2740    double * tightBounds = NULL;
2741    if (solverCharacteristics_->tryCuts())  {
2742
2743        if (numberUnsatisfied)   {
2744            // User event
2745            if (!eventHappened_ && feasible) {
2746                if (rootModels) {
2747                  // for fixings
2748                  int numberColumns=solver_->getNumCols();
2749                  tightBounds = new double [2*numberColumns];
2750                  {
2751                    const double * lower = solver_->getColLower();
2752                    const double * upper = solver_->getColUpper();
2753                    for (int i=0;i<numberColumns;i++) {
2754                      tightBounds[2*i+0]=lower[i];
2755                      tightBounds[2*i+1]=upper[i];
2756                    }
2757                  }
2758                  int numberModels = multipleRootTries_%100;
2759                  const OsiSolverInterface ** solvers = new 
2760                    const OsiSolverInterface * [numberModels];
2761                  int numberRows=continuousSolver_->getNumRows();
2762                  int maxCuts=0;
2763                  for (int i=0;i<numberModels;i++) {
2764                    solvers[i]=rootModels[i]->solver();
2765                    const double * lower = solvers[i]->getColLower();
2766                    const double * upper = solvers[i]->getColUpper();
2767                    for (int j=0;j<numberColumns;j++) {
2768                      tightBounds[2*j+0]=CoinMax(lower[j],tightBounds[2*j+0]);
2769                      tightBounds[2*j+1]=CoinMin(upper[j],tightBounds[2*j+1]);
2770                    }
2771                    int numberRows2=solvers[i]->getNumRows();
2772                    assert (numberRows2>=numberRows);
2773                    maxCuts += numberRows2-numberRows;
2774                    // accumulate statistics
2775                    for (int j=0;j<numberCutGenerators_;j++) {
2776                      generator_[j]->addStatistics(rootModels[i]->cutGenerator(j));
2777                    }
2778                  }
2779                  for (int j=0;j<numberCutGenerators_;j++) {
2780                    generator_[j]->scaleBackStatistics(numberModels);
2781                  }
2782                  //CbcRowCuts rowCut(maxCuts);
2783                  const OsiRowCutDebugger *debugger = NULL;
2784                  if ((specialOptions_&1) != 0) 
2785                    debugger = solver_->getRowCutDebugger() ;
2786                  for (int iModel=0;iModel<numberModels;iModel++) {
2787                    int numberRows2=solvers[iModel]->getNumRows();
2788                    const CoinPackedMatrix * rowCopy = solvers[iModel]->getMatrixByRow();
2789                    const int * rowLength = rowCopy->getVectorLengths(); 
2790                    const double * elements = rowCopy->getElements();
2791                    const int * column = rowCopy->getIndices();
2792                    const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
2793                    const double * rowLower = solvers[iModel]->getRowLower();
2794                    const double * rowUpper = solvers[iModel]->getRowUpper();
2795                    for (int iRow=numberRows;iRow<numberRows2;iRow++) {
2796                      OsiRowCut rc;
2797                      rc.setLb(rowLower[iRow]);
2798                      rc.setUb(rowUpper[iRow]);
2799                      CoinBigIndex start = rowStart[iRow];
2800                      rc.setRow(rowLength[iRow],column+start,elements+start,false);
2801                      if (debugger)
2802                        CoinAssert (!debugger->invalidCut(rc));
2803                      globalCuts_.addCutIfNotDuplicate(rc);
2804                    }
2805                    //int cutsAdded=globalCuts_.numberCuts()-numberCuts;
2806                    //numberCuts += cutsAdded;
2807                    //printf("Model %d gave %d cuts (out of %d possible)\n",
2808                    //     iModel,cutsAdded,numberRows2-numberRows);
2809                  }
2810                  // normally replace global cuts
2811                  //if (!globalCuts_.())
2812                  //globalCuts_=rowCutrowCut.addCuts(globalCuts_);
2813                  //rowCut.addCuts(globalCuts_);
2814                  int nTightened=0;
2815                  bool feasible=true;
2816                  {
2817                    double tolerance=1.0e-5;
2818                    const double * lower = solver_->getColLower();
2819                    const double * upper = solver_->getColUpper();
2820                    for (int i=0;i<numberColumns;i++) {
2821                      if (tightBounds[2*i+0]>tightBounds[2*i+1]) {
2822                        feasible=false;
2823                        printf("Bad bounds on %d %g,%g was %g,%g\n",
2824                               i,tightBounds[2*i+0],tightBounds[2*i+1],
2825                               lower[i],upper[i]);
2826                      }
2827                      //int k=0;
2828                      double oldLower=lower[i];
2829                      double oldUpper=upper[i];
2830                      if (tightBounds[2*i+0]>oldLower+tolerance) {
2831                        nTightened++;
2832                        //k++;
2833                        solver_->setColLower(i,tightBounds[2*i+0]);
2834                      }
2835                      if (tightBounds[2*i+1]<oldUpper-tolerance) {
2836                        nTightened++;
2837                        //k++;
2838                        solver_->setColUpper(i,tightBounds[2*i+1]);
2839                      }
2840                      //if (k)
2841                      //printf("new bounds on %d %g,%g was %g,%g\n",
2842                      //       i,tightBounds[2*i+0],tightBounds[2*i+1],
2843                      //       oldLower,oldUpper);
2844                    }
2845                    if (!feasible)
2846                      abort(); // deal with later
2847                  }
2848                  delete [] tightBounds;
2849                  tightBounds=NULL;
2850                  char printBuffer[200];
2851                  sprintf(printBuffer,"%d solvers added %d different cuts out of pool of %d",
2852                          numberModels,globalCuts_.sizeRowCuts(),maxCuts);
2853                  messageHandler()->message(CBC_GENERAL, messages())
2854                    << printBuffer << CoinMessageEol ;
2855                  if (nTightened) {
2856                    sprintf(printBuffer,"%d bounds were tightened",
2857                          nTightened);
2858                    messageHandler()->message(CBC_GENERAL, messages())
2859                      << printBuffer << CoinMessageEol ;
2860                  }
2861                  delete [] solvers;
2862                }
2863                if (!parentModel_&&(moreSpecialOptions_&67108864) != 0) {
2864                  // load cuts from file
2865                  FILE * fp = fopen("global.cuts","rb");
2866                  if (fp) {
2867                    size_t nRead;
2868                    int numberColumns=solver_->getNumCols();
2869                    int numCols;
2870                    nRead = fread(&numCols, sizeof(int), 1, fp);
2871                    if (nRead != 1)
2872                      throw("Error in fread");
2873                    if (numberColumns!=numCols) {
2874                      printf("Mismatch on columns %d %d\n",numberColumns,numCols);
2875                      fclose(fp);
2876                    } else {
2877                      // If rootModel just do some
2878                      double threshold=-1.0;
2879                      if (!multipleRootTries_)
2880                        threshold=0.5;
2881                      int initialCuts=0;
2882                      int initialGlobal = globalCuts_.sizeRowCuts();
2883                      double * elements = new double [numberColumns+2];
2884                      int * indices = new int [numberColumns];
2885                      int numberEntries=1;
2886                      while (numberEntries>0) {
2887                        nRead = fread(&numberEntries, sizeof(int), 1, fp);
2888                        if (nRead != 1)
2889                          throw("Error in fread");
2890                        double randomNumber=randomNumberGenerator_.randomDouble();
2891                        if (numberEntries>0) {
2892                          initialCuts++;
2893                          nRead = fread(elements, sizeof(double), numberEntries+2, fp);
2894                          if (nRead != static_cast<size_t>(numberEntries+2))
2895                            throw("Error in fread");
2896                          nRead = fread(indices, sizeof(int), numberEntries, fp);
2897                          if (nRead != static_cast<size_t>(numberEntries))
2898                            throw("Error in fread");
2899                          if (randomNumber>threshold) {
2900                            OsiRowCut rc;
2901                            rc.setLb(elements[numberEntries]);
2902                            rc.setUb(elements[numberEntries+1]);
2903                            rc.setRow(numberEntries,indices,elements,
2904                                      false);
2905                            rc.setGloballyValidAsInteger(2);
2906                            globalCuts_.addCutIfNotDuplicate(rc) ;
2907                          } 
2908                        }
2909                      }
2910                      fclose(fp);
2911                      // fixes
2912                      int nTightened=0;
2913                      fp = fopen("global.fix","rb");
2914                      if (fp) {
2915                        nRead = fread(indices, sizeof(int), 2, fp);
2916                        if (nRead != 2)
2917                          throw("Error in fread");
2918                        if (numberColumns!=indices[0]) {
2919                          printf("Mismatch on columns %d %d\n",numberColumns,
2920                                 indices[0]);
2921                        } else {
2922                          indices[0]=1;
2923                          while (indices[0]>=0) {
2924                            nRead = fread(indices, sizeof(int), 2, fp);
2925                            if (nRead != 2)
2926                              throw("Error in fread");
2927                            int iColumn=indices[0];
2928                            if (iColumn>=0) {
2929                              nTightened++;
2930                              nRead = fread(elements, sizeof(double), 4, fp);
2931                              if (nRead != 4)
2932                                throw("Error in fread");
2933                              solver_->setColLower(iColumn,elements[0]);
2934                              solver_->setColUpper(iColumn,elements[1]);
2935                            } 
2936                          }
2937                        }
2938                      }
2939                      if (fp)
2940                        fclose(fp);
2941                      char printBuffer[200];
2942                      sprintf(printBuffer,"%d cuts read in of which %d were unique, %d bounds tightened",
2943                             initialCuts,
2944                             globalCuts_.sizeRowCuts()-initialGlobal,nTightened); 
2945                      messageHandler()->message(CBC_GENERAL, messages())
2946                        << printBuffer << CoinMessageEol ;
2947                      delete [] elements;
2948                      delete [] indices;
2949                    }
2950                  }
2951                }
2952                feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
2953                                         NULL);
2954                if (multipleRootTries_&&
2955                    (moreSpecialOptions_&134217728)!=0) {
2956                  FILE * fp=NULL;
2957                  size_t nRead;
2958                  int numberColumns=solver_->getNumCols();
2959                  int initialCuts=0;
2960                  if ((moreSpecialOptions_&134217728)!=0) {
2961                    // append so go down to end
2962                    fp = fopen("global.cuts","r+b");
2963                    if (fp) {
2964                      int numCols;
2965                      nRead = fread(&numCols, sizeof(int), 1, fp);
2966                      if (nRead != 1)
2967                        throw("Error in fread");
2968                      if (numberColumns!=numCols) {
2969                        printf("Mismatch on columns %d %d\n",numberColumns,numCols);
2970                        fclose(fp);
2971                        fp=NULL;
2972                      }
2973                    }
2974                  }
2975                  double * elements = new double [numberColumns+2];
2976                  int * indices = new int [numberColumns];
2977                  if (fp) {
2978                    int numberEntries=1;
2979                    while (numberEntries>0) {
2980                      fpos_t position;
2981                      fgetpos(fp, &position);
2982                      nRead = fread(&numberEntries, sizeof(int), 1, fp);
2983                      if (nRead != 1)
2984                        throw("Error in fread");
2985                      if (numberEntries>0) {
2986                        initialCuts++;
2987                        nRead = fread(elements, sizeof(double), numberEntries+2, fp);
2988                        if (nRead != static_cast<size_t>(numberEntries+2))
2989                          throw("Error in fread");
2990                        nRead = fread(indices, sizeof(int), numberEntries, fp);
2991                        if (nRead != static_cast<size_t>(numberEntries))
2992                          throw("Error in fread");
2993                      } else {
2994                        // end
2995                        fsetpos(fp, &position);
2996                      }
2997                    }
2998                  } else {
2999                    fp = fopen("global.cuts","wb");
3000                    size_t nWrite;
3001                    nWrite=fwrite(&numberColumns,sizeof(int),1,fp);
3002                    if (nWrite != 1)
3003                      throw("Error in fwrite");
3004                  }
3005                  size_t nWrite;
3006                  // now append binding cuts
3007                  int numberC=continuousSolver_->getNumRows();
3008                  int numberRows=solver_->getNumRows();
3009                  printf("Saving %d cuts (up from %d)\n",
3010                         initialCuts+numberRows-numberC,initialCuts);
3011                  const double * rowLower = solver_->getRowLower();
3012                  const double * rowUpper = solver_->getRowUpper();
3013                  // Row copy
3014                  CoinPackedMatrix matrixByRow(*solver_->getMatrixByRow());
3015                  const double * elementByRow = matrixByRow.getElements();
3016                  const int * column = matrixByRow.getIndices();
3017                  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
3018                  const int * rowLength = matrixByRow.getVectorLengths();
3019                  for (int iRow=numberC;iRow<numberRows;iRow++) {
3020                    int n=rowLength[iRow];
3021                    assert (n);
3022                    CoinBigIndex start=rowStart[iRow];
3023                    memcpy(elements,elementByRow+start,n*sizeof(double));
3024                    memcpy(indices,column+start,n*sizeof(int));
3025                    elements[n]=rowLower[iRow];
3026                    elements[n+1]=rowUpper[iRow];
3027                    nWrite=fwrite(&n,sizeof(int),1,fp);
3028                    if (nWrite != 1)
3029                      throw("Error in fwrite");
3030                    nWrite=fwrite(elements,sizeof(double),n+2,fp);
3031                    if (nWrite != static_cast<size_t>(n+2))
3032                      throw("Error in fwrite");
3033                    nWrite=fwrite(indices,sizeof(int),n,fp);
3034                    if (nWrite != static_cast<size_t>(n))
3035                      throw("Error in fwrite");
3036                  }
3037                  // eof marker
3038                  int eofMarker=-1;
3039                  nWrite=fwrite(&eofMarker,sizeof(int),1,fp);
3040                  if (nWrite != 1)
3041                    throw("Error in fwrite");
3042                  fclose(fp);
3043                  // do tighter bounds (? later extra to original columns)
3044                  int nTightened=0;
3045                  const double * lower = solver_->getColLower();
3046                  const double * upper = solver_->getColUpper();
3047                  const double * originalLower = continuousSolver_->getColLower();
3048                  const double * originalUpper = continuousSolver_->getColUpper();
3049                  double tolerance=1.0e-5;
3050                  for (int i=0;i<numberColumns;i++) {
3051                    if (lower[i]>originalLower[i]+tolerance) {
3052                      nTightened++;
3053                    }
3054                    if (upper[i]<originalUpper[i]-tolerance) {
3055                      nTightened++;
3056                    }
3057                  }
3058                  if (nTightened) {
3059                    fp = fopen("global.fix","wb");
3060                    size_t nWrite;
3061                    indices[0]=numberColumns;
3062                    if (originalColumns_)
3063                      indices[1]=COIN_INT_MAX;
3064                    else
3065                      indices[1]=-1;
3066                    nWrite=fwrite(indices,sizeof(int),2,fp);
3067                    if (nWrite != 2)
3068                      throw("Error in fwrite");
3069                    for (int i=0;i<numberColumns;i++) {
3070                      int nTightened=0;
3071                      if (lower[i]>originalLower[i]+tolerance) {
3072                        nTightened++;
3073                      }
3074                      if (upper[i]<originalUpper[i]-tolerance) {
3075                        nTightened++;
3076                      }
3077                      if (nTightened) {
3078                        indices[0]=i;
3079                        if (originalColumns_)
3080                          indices[1]=originalColumns_[i];
3081                        elements[0]=lower[i];
3082                        elements[1]=upper[i];
3083                        elements[2]=originalLower[i];
3084                        elements[3]=originalUpper[i];
3085                        nWrite=fwrite(indices,sizeof(int),2,fp);
3086                        if (nWrite != 2)
3087                          throw("Error in fwrite");
3088                        nWrite=fwrite(elements,sizeof(double),4,fp);
3089                        if (nWrite != 4)
3090                          throw("Error in fwrite");
3091                      }
3092                    }
3093                    // eof marker
3094                    indices[0]=-1;
3095                    nWrite=fwrite(indices,sizeof(int),2,fp);
3096                    if (nWrite != 2)
3097                      throw("Error in fwrite");
3098                    fclose(fp);
3099                  }
3100                  delete [] elements;
3101                  delete [] indices;
3102                }
3103                if ((specialOptions_&524288) != 0 && !parentModel_
3104                        && storedRowCuts_) {
3105                    if (feasible) {
3106                        /* pick up stuff and try again
3107                        add cuts, maybe keep around
3108                        do best solution and if so new heuristics
3109                        obviously tighten bounds
3110                        */
3111                        // A and B probably done on entry
3112                        // A) tight bounds on integer variables
3113                        const double * lower = solver_->getColLower();
3114                        const double * upper = solver_->getColUpper();
3115                        const double * tightLower = storedRowCuts_->tightLower();
3116                        const double * tightUpper = storedRowCuts_->tightUpper();
3117                        int nTightened = 0;
3118                        for (int i = 0; i < numberIntegers_; i++) {
3119                            int iColumn = integerVariable_[i];
3120                            if (tightLower[iColumn] > lower[iColumn]) {
3121                                nTightened++;
3122                                solver_->setColLower(iColumn, tightLower[iColumn]);
3123                            }
3124                            if (tightUpper[iColumn] < upper[iColumn]) {
3125                                nTightened++;
3126                                solver_->setColUpper(iColumn, tightUpper[iColumn]);
3127                            }
3128                        }
3129                        if (nTightened)
3130                          COIN_DETAIL_PRINT(printf("%d tightened by alternate cuts\n", nTightened));
3131                        if (storedRowCuts_->bestObjective() < bestObjective_) {
3132                            // B) best solution
3133                            double objValue = storedRowCuts_->bestObjective();
3134                            setBestSolution(CBC_SOLUTION, objValue,
3135                                            storedRowCuts_->bestSolution()) ;
3136                            // Do heuristics
3137                            // Allow RINS
3138                            for (int i = 0; i < numberHeuristics_; i++) {
3139                                CbcHeuristicRINS * rins
3140                                = dynamic_cast<CbcHeuristicRINS *> (heuristic_[i]);
3141                                if (rins) {
3142                                    rins->setLastNode(-100);
3143                                }
3144                            }
3145                            doHeuristicsAtRoot();
3146                        }
3147#ifdef JJF_ZERO
3148                        int nCuts = storedRowCuts_->sizeRowCuts();
3149                        // add to global list
3150                        for (int i = 0; i < nCuts; i++) {
3151                            OsiRowCut newCut(*storedRowCuts_->rowCutPointer(i));
3152                            newCut.setGloballyValidAsInteger(2);
3153                            newCut.mutableRow().setTestForDuplicateIndex(false);
3154                            globalCuts_.insert(newCut) ;
3155                        }
3156#else
3157                        addCutGenerator(storedRowCuts_, -99, "Stored from previous run",
3158                                        true, false, false, -200);
3159#endif
3160                        // Set cuts as active
3161                        delete [] addedCuts_ ;
3162                        maximumNumberCuts_ = cuts.sizeRowCuts();
3163                        if (maximumNumberCuts_) {
3164                            addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
3165                        } else {
3166                            addedCuts_ = NULL;
3167                        }
3168                        for (int i = 0; i < maximumNumberCuts_; i++)
3169                            addedCuts_[i] = new CbcCountRowCut(*cuts.rowCutPtr(i),
3170                                                               NULL, -1, -1, 2);
3171                        COIN_DETAIL_PRINT(printf("size %d\n", cuts.sizeRowCuts()));
3172                        cuts = OsiCuts();
3173                        currentNumberCuts_ = maximumNumberCuts_;
3174                        feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
3175                                                 NULL);
3176                        for (int i = 0; i < maximumNumberCuts_; i++)
3177                            delete addedCuts_[i];
3178                    }
3179                    delete storedRowCuts_;
3180                    storedRowCuts_ = NULL;
3181                }
3182            } else {
3183                feasible = false;
3184            }
3185        }       else if (solverCharacteristics_->solutionAddsCuts() ||
3186                   solverCharacteristics_->alwaysTryCutsAtRootNode()) {
3187            // may generate cuts and turn the solution
3188            //to an infeasible one
3189            feasible = solveWithCuts(cuts, 2,
3190                                     NULL);
3191        }
3192    }
3193    if (rootModels) {
3194      int numberModels = multipleRootTries_%100;
3195      for (int i=0;i<numberModels;i++)
3196        delete rootModels[i];
3197      delete [] rootModels;
3198    }
3199    // check extra info on feasibility
3200    if (!solverCharacteristics_->mipFeasible())
3201        feasible = false;
3202    // If max nodes==0 - don't do strong branching
3203    if (!getMaximumNodes()) {
3204      if (feasible)
3205        feasible=false;
3206      else
3207        setMaximumNodes(1); //allow to stop on success
3208    }
3209    topOfTree_=NULL;
3210#ifdef CLP_RESOLVE
3211    if ((moreSpecialOptions_&2097152)!=0&&!parentModel_&&feasible) {
3212      OsiClpSolverInterface * clpSolver
3213        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3214      if (clpSolver)
3215        resolveClp(clpSolver,0);
3216    }
3217#endif
3218    // make cut generators less aggressive
3219    for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
3220        CglCutGenerator * generator = generator_[iCutGenerator]->generator();
3221        generator->setAggressiveness(generator->getAggressiveness() - 100);
3222    }
3223    currentNumberCuts_ = numberNewCuts_ ;
3224    if (solverCharacteristics_->solutionAddsCuts()) {
3225      // With some heuristics solver needs a resolve here (don't know if this is bug in heuristics)
3226      solver_->resolve();
3227      if(!isProvenOptimal()){
3228        solver_->initialSolve();
3229      }
3230    }
3231    // See if can stop on gap
3232    bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
3233    if(canStopOnGap()) {
3234        if (bestPossibleObjective_ < getCutoff())
3235            stoppedOnGap_ = true ;
3236        feasible = false;
3237    }
3238    // User event
3239    if (eventHappened_)
3240        feasible = false;
3241#if defined(COIN_HAS_CLP)&&defined(COIN_HAS_CPX)
3242    /*
3243      This is the notion of using Cbc stuff to get going, then calling cplex to
3244      finish off.
3245    */
3246    if (feasible && (specialOptions_&16384) != 0 && fastNodeDepth_ == -2 && !parentModel_) {
3247        // Use Cplex to do search!
3248        double time1 = CoinCpuTime();
3249        OsiClpSolverInterface * clpSolver
3250        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3251        OsiCpxSolverInterface cpxSolver;
3252        double direction = clpSolver->getObjSense();
3253        cpxSolver.setObjSense(direction);
3254        // load up cplex
3255        const CoinPackedMatrix * matrix = continuousSolver_->getMatrixByCol();
3256        const double * rowLower = continuousSolver_->getRowLower();
3257        const double * rowUpper = continuousSolver_->getRowUpper();
3258        const double * columnLower = continuousSolver_->getColLower();
3259        const double * columnUpper = continuousSolver_->getColUpper();
3260        const double * objective = continuousSolver_->getObjCoefficients();
3261        cpxSolver.loadProblem(*matrix, columnLower, columnUpper,
3262                              objective, rowLower, rowUpper);
3263        double * setSol = new double [numberIntegers_];
3264        int * setVar = new int [numberIntegers_];
3265        // cplex doesn't know about objective offset
3266        double offset = clpSolver->getModelPtr()->objectiveOffset();
3267        for (int i = 0; i < numberIntegers_; i++) {
3268            int iColumn = integerVariable_[i];
3269            cpxSolver.setInteger(iColumn);
3270            if (bestSolution_) {
3271                setSol[i] = bestSolution_[iColumn];
3272                setVar[i] = iColumn;
3273            }
3274        }
3275        CPXENVptr env = cpxSolver.getEnvironmentPtr();
3276        CPXLPptr lpPtr = cpxSolver.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL);
3277        cpxSolver.switchToMIP();
3278        if (bestSolution_) {
3279            CPXcopymipstart(env, lpPtr, numberIntegers_, setVar, setSol);
3280        }
3281        if (clpSolver->getNumRows() > continuousSolver_->getNumRows() && false) {
3282            // add cuts
3283            const CoinPackedMatrix * matrix = clpSolver->getMatrixByRow();
3284            const double * rhs = clpSolver->getRightHandSide();
3285            const char * rowSense = clpSolver->getRowSense();
3286            const double * elementByRow = matrix->getElements();
3287            const int * column = matrix->getIndices();
3288            const CoinBigIndex * rowStart = matrix->getVectorStarts();
3289            const int * rowLength = matrix->getVectorLengths();
3290            int nStart = continuousSolver_->getNumRows();
3291            int nRows = clpSolver->getNumRows();
3292            int size = rowStart[nRows-1] + rowLength[nRows-1] -
3293                       rowStart[nStart];
3294            int nAdd = 0;
3295            double * rmatval = new double [size];
3296            int * rmatind = new int [size];
3297            int * rmatbeg = new int [nRows-nStart+1];
3298            size = 0;
3299            rmatbeg[0] = 0;
3300            for (int i = nStart; i < nRows; i++) {
3301                for (int k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
3302                    rmatind[size] = column[k];
3303                    rmatval[size++] = elementByRow[k];
3304                }
3305                nAdd++;
3306                rmatbeg[nAdd] = size;
3307            }
3308            CPXaddlazyconstraints(env, lpPtr, nAdd, size,
3309                                  rhs, rowSense, rmatbeg,
3310                                  rmatind, rmatval, NULL);
3311            CPXsetintparam( env, CPX_PARAM_REDUCE,
3312                            // CPX_PREREDUCE_NOPRIMALORDUAL (0)
3313                            CPX_PREREDUCE_PRIMALONLY);
3314        }
3315        if (getCutoff() < 1.0e50) {
3316            double useCutoff = getCutoff() + offset;
3317            if (bestObjective_ < 1.0e50)
3318                useCutoff = bestObjective_ + offset + 1.0e-7;
3319            cpxSolver.setDblParam(OsiDualObjectiveLimit, useCutoff*
3320                                  direction);
3321            if ( direction > 0.0 )
3322                CPXsetdblparam( env, CPX_PARAM_CUTUP, useCutoff ) ; // min
3323            else
3324                CPXsetdblparam( env, CPX_PARAM_CUTLO, useCutoff ) ; // max
3325        }
3326        CPXsetdblparam(env, CPX_PARAM_EPGAP, dblParam_[CbcAllowableFractionGap]);
3327        delete [] setSol;
3328        delete [] setVar;
3329        char printBuffer[200];
3330        if (offset) {
3331            sprintf(printBuffer, "Add %g to all Cplex messages for true objective",
3332                    -offset);
3333            messageHandler()->message(CBC_GENERAL, messages())
3334            << printBuffer << CoinMessageEol ;
3335            cpxSolver.setDblParam(OsiObjOffset, offset);
3336        }
3337        cpxSolver.branchAndBound();
3338        double timeTaken = CoinCpuTime() - time1;
3339        sprintf(printBuffer, "Cplex took %g seconds",
3340                timeTaken);
3341        messageHandler()->message(CBC_GENERAL, messages())
3342        << printBuffer << CoinMessageEol ;
3343        numberExtraNodes_ = CPXgetnodecnt(env, lpPtr);
3344        numberExtraIterations_ = CPXgetmipitcnt(env, lpPtr);
3345        double value = cpxSolver.getObjValue() * direction;
3346        if (cpxSolver.isProvenOptimal() && value <= getCutoff()) {
3347            feasible = true;
3348            clpSolver->setWarmStart(NULL);
3349            // try and do solution
3350            double * newSolution =
3351                CoinCopyOfArray(cpxSolver.getColSolution(),
3352                                getNumCols());
3353            setBestSolution(CBC_STRONGSOL, value, newSolution) ;
3354            delete [] newSolution;
3355        }
3356        feasible = false;
3357    }
3358#endif
3359    if (!parentModel_&&(moreSpecialOptions_&268435456) != 0) {
3360      // try idiotic idea
3361      CbcObject * obj = new CbcIdiotBranch(this);
3362      obj->setPriority(1); // temp
3363      addObjects(1, &obj);
3364      delete obj;
3365    }
3366   
3367    /*
3368      A hook to use clp to quickly explore some part of the tree.
3369    */
3370    if (fastNodeDepth_ == 1000 &&/*!parentModel_*/(specialOptions_&2048) == 0) {
3371        fastNodeDepth_ = -1;
3372        CbcObject * obj =
3373            new CbcFollowOn(this);
3374        obj->setPriority(1);
3375        addObjects(1, &obj);
3376        delete obj;
3377    }
3378    int saveNumberSolves = numberSolves_;
3379    int saveNumberIterations = numberIterations_;
3380    if ((fastNodeDepth_ >= 0||(moreSpecialOptions_&33554432)!=0)
3381        &&/*!parentModel_*/(specialOptions_&2048) == 0) {
3382        // add in a general depth object doClp
3383        int type = (fastNodeDepth_ <= 100) ? fastNodeDepth_ : -(fastNodeDepth_ - 100);
3384        if ((moreSpecialOptions_&33554432)!=0)
3385          type=12;
3386        else
3387          fastNodeDepth_ += 1000000;     // mark as done
3388        CbcObject * obj =
3389            new CbcGeneralDepth(this, type);
3390        addObjects(1, &obj);
3391        delete obj;
3392        // fake number of objects
3393        numberObjects_--;
3394        if (parallelMode() < -1) {
3395            // But make sure position is correct
3396            OsiObject * obj2 = object_[numberObjects_];
3397            obj = dynamic_cast<CbcObject *> (obj2);
3398            assert (obj);
3399            obj->setPosition(numberObjects_);
3400        }
3401    }
3402#ifdef COIN_HAS_CLP
3403#ifdef NO_CRUNCH
3404    if (true) {
3405        OsiClpSolverInterface * clpSolver
3406        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3407        if (clpSolver && !parentModel_) {
3408            ClpSimplex * clpSimplex = clpSolver->getModelPtr();
3409            clpSimplex->setSpecialOptions(clpSimplex->specialOptions() | 131072);
3410            //clpSimplex->startPermanentArrays();
3411            clpSimplex->setPersistenceFlag(2);
3412        }
3413    }
3414#endif
3415#endif
3416// Save copy of solver
3417    OsiSolverInterface * saveSolver = NULL;
3418    if (!parentModel_ && (specialOptions_&(512 + 32768)) != 0)
3419        saveSolver = solver_->clone();
3420    double checkCutoffForRestart = 1.0e100;
3421    saveModel(saveSolver, &checkCutoffForRestart, &feasible);
3422    if ((specialOptions_&262144) != 0 && !parentModel_) {
3423        // Save stuff and return!
3424        storedRowCuts_->saveStuff(bestObjective_, bestSolution_,
3425                                  solver_->getColLower(),
3426                                  solver_->getColUpper());
3427        delete [] lowerBefore;
3428        delete [] upperBefore;
3429        delete saveSolver;
3430        if (flipObjective)
3431          flipModel();
3432        return;
3433    }
3434    /*
3435      We've taken the continuous relaxation as far as we can. Time to branch.
3436      The first order of business is to actually create a node. chooseBranch
3437      currently uses strong branching to evaluate branch object candidates,
3438      unless forced back to simple branching. If chooseBranch concludes that a
3439      branching candidate is monotone (anyAction == -1) or infeasible (anyAction
3440      == -2) when forced to integer values, it returns here immediately.
3441
3442      Monotone variables trigger a call to resolve(). If the problem remains
3443      feasible, try again to choose a branching variable. At the end of the loop,
3444      resolved == true indicates that some variables were fixed.
3445
3446      Loss of feasibility will result in the deletion of newNode.
3447    */
3448
3449    bool resolved = false ;
3450    CbcNode *newNode = NULL ;
3451    numberFixedAtRoot_ = 0;
3452    numberFixedNow_ = 0;
3453    if (!parentModel_&&(moreSpecialOptions2_&2)!=0) {
3454#ifdef COIN_HAS_CLP
3455      OsiClpSolverInterface * clpSolver
3456        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3457      if (clpSolver) {
3458        if (getCutoff()>1.0e20) {
3459          printf("Zapping costs\n");
3460          int numberColumns=solver_->getNumCols();
3461          double * zeroCost = new double [numberColumns];
3462          // could make small random
3463          memset(zeroCost,0,numberColumns*sizeof(double));
3464          solver_->setObjective(zeroCost);
3465          double objValue = solver_->getObjValue();
3466          solver_->setDblParam(OsiObjOffset,-objValue);
3467          clpSolver->getModelPtr()->setObjectiveValue(objValue);
3468          delete [] zeroCost;
3469        } else {
3470          moreSpecialOptions2_ &= ~2;
3471        }
3472      } else {
3473#endif
3474          moreSpecialOptions2_ &= ~2;
3475#ifdef COIN_HAS_CLP
3476      }
3477#endif
3478    }
3479    int numberIterationsAtContinuous = numberIterations_;
3480    //solverCharacteristics_->setSolver(solver_);
3481    if (feasible) {
3482#define HOTSTART -1
3483#if HOTSTART<0
3484        if (bestSolution_ && !parentModel_ && !hotstartSolution_ &&
3485                (moreSpecialOptions_&1024) != 0) {
3486            // Set priorities so only branch on ones we need to
3487            // use djs and see if only few branches needed
3488#ifndef NDEBUG
3489            double integerTolerance = getIntegerTolerance() ;
3490#endif
3491            bool possible = true;
3492            const double * saveLower = continuousSolver_->getColLower();
3493            const double * saveUpper = continuousSolver_->getColUpper();
3494            for (int i = 0; i < numberObjects_; i++) {
3495                const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object_[i]);
3496                if (thisOne) {
3497                    int iColumn = thisOne->columnNumber();
3498                    if (saveUpper[iColumn] > saveLower[iColumn] + 1.5) {
3499                        possible = false;
3500                        break;
3501                    }
3502                } else {
3503                    possible = false;
3504                    break;
3505                }
3506            }
3507            if (possible) {
3508                OsiSolverInterface * solver = continuousSolver_->clone();
3509                int numberColumns = solver->getNumCols();
3510                for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
3511                    double value = bestSolution_[iColumn] ;
3512                    value = CoinMax(value, saveLower[iColumn]) ;
3513                    value = CoinMin(value, saveUpper[iColumn]) ;
3514                    value = floor(value + 0.5);
3515                    if (solver->isInteger(iColumn)) {
3516                        solver->setColLower(iColumn, value);
3517                        solver->setColUpper(iColumn, value);
3518                    }
3519                }
3520                solver->setHintParam(OsiDoDualInResolve, false, OsiHintTry);
3521                // objlim and all slack
3522                double direction = solver->getObjSense();
3523                solver->setDblParam(OsiDualObjectiveLimit, 1.0e50*direction);
3524                CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis *> (solver->getEmptyWarmStart());
3525                solver->setWarmStart(basis);
3526                delete basis;
3527                bool changed = true;
3528                hotstartPriorities_ = new int [numberColumns];
3529                for (int iColumn = 0; iColumn < numberColumns; iColumn++)
3530                    hotstartPriorities_[iColumn] = 1;
3531                while (changed) {
3532                    changed = false;
3533                    solver->resolve();
3534                    if (!solver->isProvenOptimal()) {
3535                        possible = false;
3536                        break;
3537                    }
3538                    const double * dj = solver->getReducedCost();
3539                    const double * colLower = solver->getColLower();
3540                    const double * colUpper = solver->getColUpper();
3541                    const double * solution = solver->getColSolution();
3542                    int nAtLbNatural = 0;
3543                    int nAtUbNatural = 0;
3544                    int nZeroDj = 0;
3545                    int nForced = 0;
3546                    for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
3547                        double value = solution[iColumn] ;
3548                        value = CoinMax(value, saveLower[iColumn]) ;
3549                        value = CoinMin(value, saveUpper[iColumn]) ;
3550                        if (solver->isInteger(iColumn)) {
3551                            assert(fabs(value - solution[iColumn]) <= integerTolerance) ;
3552                            if (hotstartPriorities_[iColumn] == 1) {
3553                                if (dj[iColumn] < -1.0e-6) {
3554                                    // negative dj
3555                                    if (saveUpper[iColumn] == colUpper[iColumn]) {
3556                                        nAtUbNatural++;
3557                                        hotstartPriorities_[iColumn] = 2;
3558                                        solver->setColLower(iColumn, saveLower[iColumn]);
3559                                        solver->setColUpper(iColumn, saveUpper[iColumn]);
3560                                    } else {
3561                                        nForced++;
3562                                    }
3563                                } else if (dj[iColumn] > 1.0e-6) {
3564                                    // positive dj
3565                                    if (saveLower[iColumn] == colLower[iColumn]) {
3566                                        nAtLbNatural++;
3567                                        hotstartPriorities_[iColumn] = 2;
3568                                        solver->setColLower(iColumn, saveLower[iColumn]);
3569                                        solver->setColUpper(iColumn, saveUpper[iColumn]);
3570                                    } else {
3571                                        nForced++;
3572                                    }
3573                                } else {
3574                                    // zero dj
3575                                    nZeroDj++;
3576                                }
3577                            }
3578                        }
3579                    }
3580#ifdef CLP_INVESTIGATE
3581                    printf("%d forced, %d naturally at lower, %d at upper - %d zero dj\n",
3582                           nForced, nAtLbNatural, nAtUbNatural, nZeroDj);
3583#endif
3584                    if (nAtLbNatural || nAtUbNatural) {
3585                        changed = true;
3586                    } else {
3587                        if (nForced + nZeroDj > 5000 ||
3588                                (nForced + nZeroDj)*2 > numberIntegers_)
3589                            possible = false;
3590                    }
3591                }
3592                delete solver;
3593            }
3594            if (possible) {
3595                setHotstartSolution(bestSolution_);
3596                if (!saveCompare) {
3597                    // create depth first comparison
3598                    saveCompare = nodeCompare_;
3599                    // depth first
3600                    nodeCompare_ = new CbcCompareDepth();
3601                    tree_->setComparison(*nodeCompare_) ;
3602                }
3603            } else {
3604                delete [] hotstartPriorities_;
3605                hotstartPriorities_ = NULL;
3606            }
3607        }
3608#endif
3609#if HOTSTART>0
3610        if (hotstartSolution_ && !hotstartPriorities_) {
3611            // Set up hot start
3612            OsiSolverInterface * solver = solver_->clone();
3613            double direction = solver_->getObjSense() ;
3614            int numberColumns = solver->getNumCols();
3615            double * saveLower = CoinCopyOfArray(solver->getColLower(), numberColumns);
3616            double * saveUpper = CoinCopyOfArray(solver->getColUpper(), numberColumns);
3617            // move solution
3618            solver->setColSolution(hotstartSolution_);
3619            // point to useful information
3620            const double * saveSolution = testSolution_;
3621            testSolution_ = solver->getColSolution();
3622            OsiBranchingInformation usefulInfo = usefulInformation();
3623            testSolution_ = saveSolution;
3624            /*
3625            Run through the objects and use feasibleRegion() to set variable bounds
3626            so as to fix the variables specified in the objects at their value in this
3627            solution. Since the object list contains (at least) one object for every
3628            integer variable, this has the effect of fixing all integer variables.
3629            */
3630            for (int i = 0; i < numberObjects_; i++)
3631                object_[i]->feasibleRegion(solver, &usefulInfo);
3632            solver->resolve();
3633            assert (solver->isProvenOptimal());
3634            double gap = CoinMax((solver->getObjValue() - solver_->getObjValue()) * direction, 0.0) ;
3635            const double * dj = solver->getReducedCost();
3636            const double * colLower = solver->getColLower();
3637            const double * colUpper = solver->getColUpper();
3638            const double * solution = solver->getColSolution();
3639            int nAtLbNatural = 0;
3640            int nAtUbNatural = 0;
3641            int nAtLbNaturalZero = 0;
3642            int nAtUbNaturalZero = 0;
3643            int nAtLbFixed = 0;
3644            int nAtUbFixed = 0;
3645            int nAtOther = 0;
3646            int nAtOtherNatural = 0;
3647            int nNotNeeded = 0;
3648            delete [] hotstartSolution_;
3649            hotstartSolution_ = new double [numberColumns];
3650            delete [] hotstartPriorities_;
3651            hotstartPriorities_ = new int [numberColumns];
3652            int * order = (int *) saveUpper;
3653            int nFix = 0;
3654            double bestRatio = COIN_DBL_MAX;
3655            for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
3656                double value = solution[iColumn] ;
3657                value = CoinMax(value, saveLower[iColumn]) ;
3658                value = CoinMin(value, saveUpper[iColumn]) ;
3659                double sortValue = COIN_DBL_MAX;
3660                if (solver->isInteger(iColumn)) {
3661                    assert(fabs(value - solution[iColumn]) <= 1.0e-5) ;
3662                    double value2 = floor(value + 0.5);
3663                    if (dj[iColumn] < -1.0e-6) {
3664                        // negative dj
3665                        //assert (value2==colUpper[iColumn]);
3666                        if (saveUpper[iColumn] == colUpper[iColumn]) {
3667                            nAtUbNatural++;
3668                            sortValue = 0.0;
3669                            double value = -dj[iColumn];
3670                            if (value > gap)
3671                                nFix++;
3672                            else if (gap < value*bestRatio)
3673                                bestRatio = gap / value;
3674                            if (saveLower[iColumn] != colLower[iColumn]) {
3675                                nNotNeeded++;
3676                                sortValue = 1.0e20;
3677                            }
3678                        } else if (saveLower[iColumn] == colUpper[iColumn]) {
3679                            nAtLbFixed++;
3680                            sortValue = dj[iColumn];
3681                        } else {
3682                            nAtOther++;
3683                            sortValue = 0.0;
3684                            if (saveLower[iColumn] != colLower[iColumn] &&
3685                                    saveUpper[iColumn] != colUpper[iColumn]) {
3686                                nNotNeeded++;
3687                                sortValue = 1.0e20;
3688                            }
3689                        }
3690                    } else if (dj[iColumn] > 1.0e-6) {
3691                        // positive dj
3692                        //assert (value2==colLower[iColumn]);
3693                        if (saveLower[iColumn] == colLower[iColumn]) {
3694                            nAtLbNatural++;
3695                            sortValue = 0.0;
3696                            double value = dj[iColumn];
3697                            if (value > gap)
3698                                nFix++;
3699                            else if (gap < value*bestRatio)
3700                                bestRatio = gap / value;
3701                            if (saveUpper[iColumn] != colUpper[iColumn]) {
3702                                nNotNeeded++;
3703                                sortValue = 1.0e20;
3704                            }
3705                        } else if (saveUpper[iColumn] == colLower[iColumn]) {
3706                            nAtUbFixed++;
3707                            sortValue = -dj[iColumn];
3708                        } else {
3709                            nAtOther++;
3710                            sortValue = 0.0;
3711                            if (saveLower[iColumn] != colLower[iColumn] &&
3712                                    saveUpper[iColumn] != colUpper[iColumn]) {
3713                                nNotNeeded++;
3714                                sortValue = 1.0e20;
3715                            }
3716                        }
3717                    } else {
3718                        // zero dj
3719                        if (value2 == saveUpper[iColumn]) {
3720                            nAtUbNaturalZero++;
3721                            sortValue = 0.0;
3722                            if (saveLower[iColumn] != colLower[iColumn]) {
3723                                nNotNeeded++;
3724                                sortValue = 1.0e20;
3725                            }
3726                        } else if (value2 == saveLower[iColumn]) {
3727                            nAtLbNaturalZero++;
3728                            sortValue = 0.0;
3729                        } else {
3730                            nAtOtherNatural++;
3731                            sortValue = 0.0;
3732                            if (saveLower[iColumn] != colLower[iColumn] &&
3733                                    saveUpper[iColumn] != colUpper[iColumn]) {
3734                                nNotNeeded++;
3735                                sortValue = 1.0e20;
3736                            }
3737                        }
3738                    }
3739#if HOTSTART==3
3740                    sortValue = -fabs(dj[iColumn]);
3741#endif
3742                }
3743                hotstartSolution_[iColumn] = value ;
3744                saveLower[iColumn] = sortValue;
3745                order[iColumn] = iColumn;
3746            }
3747            COIN_DETAIL_PRINT(printf("** can fix %d columns - best ratio for others is %g on gap of %g\n",
3748                                     nFix, bestRatio, gap));
3749            int nNeg = 0;
3750            CoinSort_2(saveLower, saveLower + numberColumns, order);
3751            for (int i = 0; i < numberColumns; i++) {
3752                if (saveLower[i] < 0.0) {
3753                    nNeg++;
3754#if HOTSTART==2||HOTSTART==3
3755                    // swap sign ?
3756                    saveLower[i] = -saveLower[i];
3757#endif
3758                }
3759            }
3760            CoinSort_2(saveLower, saveLower + nNeg, order);
3761            for (int i = 0; i < numberColumns; i++) {
3762#if HOTSTART==1
3763                hotstartPriorities_[order[i]] = 100;
3764#else
3765                hotstartPriorities_[order[i]] = -(i + 1);
3766#endif
3767            }
3768            COIN_DETAIL_PRINT(printf("nAtLbNat %d,nAtUbNat %d,nAtLbNatZero %d,nAtUbNatZero %d,nAtLbFixed %d,nAtUbFixed %d,nAtOther %d,nAtOtherNat %d, useless %d %d\n",
3769                   nAtLbNatural,
3770                   nAtUbNatural,
3771                   nAtLbNaturalZero,
3772                   nAtUbNaturalZero,
3773                   nAtLbFixed,
3774                   nAtUbFixed,
3775                   nAtOther,
3776                                     nAtOtherNatural, nNotNeeded, nNeg));
3777            delete [] saveLower;
3778            delete [] saveUpper;
3779            if (!saveCompare) {
3780                // create depth first comparison
3781                saveCompare = nodeCompare_;
3782                // depth first
3783                nodeCompare_ = new CbcCompareDepth();
3784                tree_->setComparison(*nodeCompare_) ;
3785            }
3786        }
3787#endif
3788        newNode = new CbcNode ;
3789        // Set objective value (not so obvious if NLP etc)
3790        setObjectiveValue(newNode, NULL);
3791        anyAction = -1 ;
3792        // To make depth available we may need a fake node
3793        CbcNode fakeNode;
3794        if (!currentNode_) {
3795            // Not true if sub trees assert (!numberNodes_);
3796            currentNode_ = &fakeNode;
3797        }
3798        phase_ = 3;
3799        // only allow 1000 passes
3800        int numberPassesLeft = 1000;
3801        // This is first crude step
3802        if (numberAnalyzeIterations_) {
3803            delete [] analyzeResults_;
3804            analyzeResults_ = new double [4*numberIntegers_];
3805            numberFixedAtRoot_ = newNode->analyze(this, analyzeResults_);
3806            if (numberFixedAtRoot_ > 0) {
3807              COIN_DETAIL_PRINT(printf("%d fixed by analysis\n", numberFixedAtRoot_));
3808                setPointers(solver_);
3809                numberFixedNow_ = numberFixedAtRoot_;
3810            } else if (numberFixedAtRoot_ < 0) {
3811              COIN_DETAIL_PRINT(printf("analysis found to be infeasible\n"));
3812                anyAction = -2;
3813                delete newNode ;
3814                newNode = NULL ;
3815                feasible = false ;
3816            }
3817        }
3818        OsiSolverBranch * branches = NULL;
3819        anyAction = chooseBranch(newNode, numberPassesLeft, NULL, cuts, resolved,
3820                                 NULL, NULL, NULL, branches);
3821        if (anyAction == -2 || newNode->objectiveValue() >= cutoff) {
3822            if (anyAction != -2) {
3823                // zap parent nodeInfo
3824#ifdef COIN_DEVELOP
3825                printf("zapping CbcNodeInfo %x\n", reinterpret_cast<int>(newNode->nodeInfo()->parent()));
3826#endif
3827                if (newNode->nodeInfo())
3828                    newNode->nodeInfo()->nullParent();
3829            }
3830            delete newNode ;
3831            newNode = NULL ;
3832            feasible = false ;
3833        }
3834    }
3835    if (newNode && probingInfo_) {
3836        int number01 = probingInfo_->numberIntegers();
3837        //const fixEntry * entry = probingInfo_->fixEntries();
3838        const int * toZero = probingInfo_->toZero();
3839        //const int * toOne = probingInfo_->toOne();
3840        //const int * integerVariable = probingInfo_->integerVariable();
3841        if (toZero[number01]) {
3842            CglTreeProbingInfo info(*probingInfo_);
3843            if ((moreSpecialOptions_&1048576)!=0&&!parentModel_) {
3844              /*
3845                Marginal idea. Further exploration probably good. Build some extra
3846                cliques from probing info. Not quite worth the effort?
3847              */
3848              CglProbing generator1;
3849              generator1.setUsingObjective(false);
3850              generator1.setMaxPass(1);
3851              generator1.setMaxPassRoot(1);
3852              generator1.setMaxLook(100);
3853              generator1.setRowCuts(3);
3854              generator1.setMaxElements(300);
3855              generator1.setMaxProbeRoot(solver_->getNumCols());
3856              CoinThreadRandom randomGenerator;
3857              //CglTreeProbingInfo info(solver_);
3858              info.level = 0;
3859              info.formulation_rows = solver_->getNumRows();
3860              info.inTree = false;
3861              info.randomNumberGenerator=&randomGenerator;
3862              info.pass=4;
3863              generator1.setMode(8);
3864              OsiCuts cs;
3865              generator1.generateCutsAndModify(*solver_,cs,&info);
3866              // very clunky
3867              OsiSolverInterface * temp = generator1.cliqueModel(solver_,2);
3868              CglPreProcess dummy;
3869              OsiSolverInterface * newSolver=dummy.cliqueIt(*temp,0.0001);
3870              delete temp;
3871              OsiSolverInterface * fake = NULL;
3872              if (newSolver) {
3873#if 0
3874                int numberCliques = generator1.numberCliques();
3875                cliqueEntry * entry = generator1.cliqueEntry();
3876                cliqueType * type = new cliqueType [numberCliques];
3877                int * start = new int [numberCliques+1];
3878                start[numberCliques]=2*numberCliques;
3879                int n=0;
3880                for (int i=0;i<numberCliques;i++) {
3881                  start[i]=2*i;
3882                  setOneFixesInCliqueEntry(entry[2*i],true);
3883                  setOneFixesInCliqueEntry(entry[2*i+1],true);
3884                  type[i]=0;
3885                }
3886                fake = info.analyze(*solver_, 1,numberCliques,start,
3887                                    entry,type);
3888                delete [] type;
3889                delete [] entry;
3890#else
3891                fake = info.analyze(*newSolver, 1,-1);
3892#endif
3893                delete newSolver;
3894              } else {
3895                fake = info.analyze(*solver_, 1);
3896              }
3897              if (fake) {
3898                //fake->writeMps("fake");
3899                CglFakeClique cliqueGen(fake);
3900                cliqueGen.setStarCliqueReport(false);
3901                cliqueGen.setRowCliqueReport(false);
3902                cliqueGen.setMinViolation(0.1);
3903                addCutGenerator(&cliqueGen, 1, "Fake cliques", true, false, false, -200);
3904                generator_[numberCutGenerators_-1]->setTiming(true);
3905                for (int i = 0; i < numberCutGenerators_; i++) {
3906                  CglKnapsackCover * cutGen =
3907                  dynamic_cast<CglKnapsackCover *>(generator_[i]->generator());
3908                  if (cutGen) {
3909                    cutGen->createCliques(*fake,2,200,false);
3910                  }
3911                }
3912              }
3913            }
3914            if (probingInfo_->packDown()) {
3915#ifdef CLP_INVESTIGATE
3916                printf("%d implications on %d 0-1\n", toZero[number01], number01);
3917#endif
3918                // Create a cut generator that remembers implications discovered at root.
3919                CglImplication implication(probingInfo_);
3920                addCutGenerator(&implication, 1, "ImplicationCuts", true, false, false, -200);
3921                generator_[numberCutGenerators_-1]->setGlobalCuts(true);
3922            } else {
3923                delete probingInfo_;
3924                probingInfo_ = NULL;
3925            }
3926        } else {
3927            delete probingInfo_;
3928
3929            probingInfo_ = NULL;
3930        }
3931    }
3932    /*
3933      At this point, the root subproblem is infeasible or fathomed by bound
3934      (newNode == NULL), or we're live with an objective value that satisfies the
3935      current objective cutoff.
3936    */
3937    assert (!newNode || newNode->objectiveValue() <= cutoff) ;
3938    // Save address of root node as we don't want to delete it
3939    /*
3940      The common case is that the lp relaxation is feasible but doesn't satisfy
3941      integrality (i.e., newNode->branchingObject(), indicating we've been able to
3942      select a branching variable). Remove any cuts that have gone slack due to
3943      forcing monotone variables. Then tack on an CbcFullNodeInfo object and full
3944      basis (via createInfo()) and stash the new cuts in the nodeInfo (via
3945      addCuts()). If, by some miracle, we have an integral solution at the root
3946      (newNode->branchingObject() is NULL), takeOffCuts() will ensure that the solver holds
3947      a valid solution for use by setBestSolution().
3948    */
3949    CoinWarmStartBasis *lastws = NULL ;
3950    if (feasible && newNode->branchingObject()) {
3951        if (resolved) {
3952            takeOffCuts(cuts, false, NULL) ;
3953#     ifdef CHECK_CUT_COUNTS
3954            { printf("Number of rows after chooseBranch fix (root)"
3955                         "(active only) %d\n",
3956                         numberRowsAtContinuous_ + numberNewCuts_ + numberOldActiveCuts_) ;
3957                const CoinWarmStartBasis* debugws =
3958                    dynamic_cast <const CoinWarmStartBasis*>(solver_->getWarmStart()) ;
3959                debugws->print() ;
3960                delete debugws ;
3961            }
3962#     endif
3963        }
3964        //newNode->createInfo(this,NULL,NULL,NULL,NULL,0,0) ;
3965        //newNode->nodeInfo()->addCuts(cuts,
3966        //                       newNode->numberBranches(),whichGenerator_) ;
3967        if (lastws) delete lastws ;
3968        lastws = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
3969    }
3970    /*
3971      Continuous data to be used later
3972    */
3973    continuousObjective_ = solver_->getObjValue() * solver_->getObjSense();
3974    continuousInfeasibilities_ = 0 ;
3975    if (newNode) {
3976        continuousObjective_ = newNode->objectiveValue() ;
3977        delete [] continuousSolution_;
3978        continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
3979                                              numberColumns);
3980        continuousInfeasibilities_ = newNode->numberUnsatisfied() ;
3981    }
3982    /*
3983      Bound may have changed so reset in objects
3984    */
3985    {
3986        int i ;
3987        for (i = 0; i < numberObjects_; i++)
3988            object_[i]->resetBounds(solver_) ;
3989    }
3990    /*
3991      Feasible? Then we should have either a live node prepped for future
3992      expansion (indicated by variable() >= 0), or (miracle of miracles) an
3993      integral solution at the root node.
3994
3995      initializeInfo sets the reference counts in the nodeInfo object.  Since
3996      this node is still live, push it onto the heap that holds the live set.
3997    */
3998    if (newNode) {
3999        if (newNode->branchingObject()) {
4000            newNode->initializeInfo() ;
4001            tree_->push(newNode) ;
4002            // save pointer to root node - so can pick up bounds
4003            if (!topOfTree_)
4004              topOfTree_ = dynamic_cast<CbcFullNodeInfo *>(newNode->nodeInfo()) ;
4005            if (statistics_) {
4006                if (numberNodes2_ == maximumStatistics_) {
4007                    maximumStatistics_ = 2 * maximumStatistics_;
4008                    CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_];
4009                    memset(temp, 0, maximumStatistics_*sizeof(CbcStatistics *));
4010                    memcpy(temp, statistics_, numberNodes2_*sizeof(CbcStatistics *));
4011                    delete [] statistics_;
4012                    statistics_ = temp;
4013                }
4014                assert (!statistics_[numberNodes2_]);
4015                statistics_[numberNodes2_] = new CbcStatistics(newNode, this);
4016            }
4017            numberNodes2_++;
4018#     ifdef CHECK_NODE
4019            printf("Node %x on tree\n", newNode) ;
4020#     endif
4021        } else {
4022            // continuous is integer
4023            double objectiveValue = newNode->objectiveValue();
4024            setBestSolution(CBC_SOLUTION, objectiveValue,
4025                            solver_->getColSolution()) ;
4026            if (eventHandler) {
4027              // we are stopping anyway so no need to test return code
4028              eventHandler->event(CbcEventHandler::solution);
4029            }
4030            delete newNode ;
4031            newNode = NULL ;
4032        }
4033    }
4034
4035    if (printFrequency_ <= 0) {
4036        printFrequency_ = 1000 ;
4037        if (getNumCols() > 2000)
4038            printFrequency_ = 100 ;
4039    }
4040    /*
4041      It is possible that strong branching fixes one variable and then the code goes round
4042      again and again.  This can take too long.  So we need to warn user - just once.
4043    */
4044    numberLongStrong_ = 0;
4045    CbcNode * createdNode = NULL;
4046#ifdef CBC_THREAD
4047    if ((specialOptions_&2048) != 0)
4048        numberThreads_ = 0;
4049    if (numberThreads_ ) {
4050        nodeCompare_->sayThreaded(); // need to use addresses
4051        master_ = new CbcBaseModel(*this,
4052                                   (parallelMode() < -1) ? 1 : 0);
4053        masterThread_ = master_->masterThread();
4054    }
4055#endif
4056#ifdef COIN_HAS_CLP
4057    {
4058        OsiClpSolverInterface * clpSolver
4059        = dynamic_cast<OsiClpSolverInterface *> (solver_);
4060        if (clpSolver && !parentModel_) {
4061            clpSolver->computeLargestAway();
4062        }
4063    }
4064#endif
4065    /*
4066      At last, the actual branch-and-cut search loop, which will iterate until
4067      the live set is empty or we hit some limit (integrality gap, time, node
4068      count, etc.). The overall flow is to rebuild a subproblem, reoptimise using
4069      solveWithCuts(), choose a branching pattern with chooseBranch(), and finally
4070      add the node to the live set.
4071
4072      The first action is to winnow the live set to remove nodes which are worse
4073      than the current objective cutoff.
4074    */
4075    if (solver_->getRowCutDebuggerAlways()) {
4076        OsiRowCutDebugger * debuggerX = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
4077        const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger() ;
4078        if (!debugger) {
4079            // infeasible!!
4080            printf("before search\n");
4081            const double * lower = solver_->getColLower();
4082            const double * upper = solver_->getColUpper();
4083            const double * solution = debuggerX->optimalSolution();
4084            int numberColumns = solver_->getNumCols();
4085            for (int i = 0; i < numberColumns; i++) {
4086                if (solver_->isInteger(i)) {
4087                    if (solution[i] < lower[i] - 1.0e-6 || solution[i] > upper[i] + 1.0e-6)
4088                        printf("**** ");
4089                    printf("%d %g <= %g <= %g\n",
4090                           i, lower[i], solution[i], upper[i]);
4091                }
4092            }
4093            //abort();
4094        }
4095    }
4096    {
4097        // may be able to change cutoff now
4098        double cutoff = getCutoff();
4099        double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
4100        if (cutoff > bestObjective_ - increment) {
4101            cutoff = bestObjective_ - increment ;
4102            setCutoff(cutoff) ;
4103        }
4104    }
4105#ifdef CBC_THREAD
4106    bool goneParallel = false;
4107#endif
4108#define MAX_DEL_NODE 1
4109    CbcNode * delNode[MAX_DEL_NODE+1];
4110    int nDeleteNode = 0;
4111    // For Printing etc when parallel
4112    int lastEvery1000 = 0;
4113    int lastPrintEvery = 0;
4114    int numberConsecutiveInfeasible = 0;
4115#define PERTURB_IN_FATHOM
4116#ifdef PERTURB_IN_FATHOM
4117    // allow in fathom
4118    if ((moreSpecialOptions_& 262144) != 0)
4119      specialOptions_ |= 131072;
4120#endif
4121    while (true) {
4122        lockThread();
4123#ifdef COIN_HAS_CLP
4124        // See if we want dantzig row choice
4125        goToDantzig(100, savePivotMethod);
4126#endif
4127        if (tree_->empty()) {
4128#ifdef CBC_THREAD
4129            if (parallelMode() > 0 && master_) {
4130                int anyLeft = master_->waitForThreadsInTree(0);
4131                if (!anyLeft) {
4132                    master_->stopThreads(-1);
4133                    break;
4134                }
4135            } else {
4136                break;
4137            }
4138#else
4139            break;
4140#endif
4141        } else {
4142            unlockThread();
4143        }
4144        // If done 100 nodes see if worth trying reduction
4145        if (numberNodes_ == 50 || numberNodes_ == 100) {
4146#ifdef COIN_HAS_CLP
4147            OsiClpSolverInterface * clpSolver
4148            = dynamic_cast<OsiClpSolverInterface *> (solver_);
4149            if (clpSolver && ((specialOptions_&131072) == 0) && true) {
4150                ClpSimplex * simplex = clpSolver->getModelPtr();
4151                int perturbation = simplex->perturbation();
4152#ifdef CLP_INVESTIGATE
4153                printf("Testing its n,s %d %d solves n,s %d %d - pert %d\n",
4154                       numberIterations_, saveNumberIterations,
4155                       numberSolves_, saveNumberSolves, perturbation);
4156#endif
4157                if (perturbation == 50 && (numberIterations_ - saveNumberIterations) <
4158                        8*(numberSolves_ - saveNumberSolves)) {
4159                    // switch off perturbation
4160                    simplex->setPerturbation(100);
4161#ifdef CLP_INVESTIGATE
4162                    printf("Perturbation switched off\n");
4163#endif
4164                }
4165            }
4166#endif
4167            /*
4168              Decide if we want to do a restart.
4169            */
4170            if (saveSolver && (specialOptions_&(512 + 32768)) != 0) {
4171                bool tryNewSearch = solverCharacteristics_->reducedCostsAccurate() &&
4172                                    (getCutoff() < 1.0e20 && getCutoff() < checkCutoffForRestart);
4173                int numberColumns = getNumCols();
4174                if (tryNewSearch) {
4175                    checkCutoffForRestart = getCutoff() ;
4176#ifdef CLP_INVESTIGATE
4177                    printf("after %d nodes, cutoff %g - looking\n",
4178                           numberNodes_, getCutoff());
4179#endif
4180                    saveSolver->resolve();
4181                    double direction = saveSolver->getObjSense() ;
4182                    double gap = checkCutoffForRestart - saveSolver->getObjValue() * direction ;
4183                    double tolerance;
4184                    saveSolver->getDblParam(OsiDualTolerance, tolerance) ;
4185                    if (gap <= 0.0)
4186                        gap = tolerance;
4187                    gap += 100.0 * tolerance;
4188                    double integerTolerance = getDblParam(CbcIntegerTolerance) ;
4189
4190                    const double *lower = saveSolver->getColLower() ;
4191                    const double *upper = saveSolver->getColUpper() ;
4192                    const double *solution = saveSolver->getColSolution() ;
4193                    const double *reducedCost = saveSolver->getReducedCost() ;
4194
4195                    int numberFixed = 0 ;
4196                    int numberFixed2 = 0;
4197#ifdef COIN_DEVELOP
4198                    printf("gap %g\n", gap);
4199#endif
4200                    for (int i = 0 ; i < numberIntegers_ ; i++) {
4201                        int iColumn = integerVariable_[i] ;
4202                        double djValue = direction * reducedCost[iColumn] ;
4203                        if (upper[iColumn] - lower[iColumn] > integerTolerance) {
4204                            if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) {
4205                                //printf("%d to lb on dj of %g - bounds %g %g\n",
4206                                //     iColumn,djValue,lower[iColumn],upper[iColumn]);
4207                                saveSolver->setColUpper(iColumn, lower[iColumn]) ;
4208                                numberFixed++ ;
4209                            } else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) {
4210                                //printf("%d to ub on dj of %g - bounds %g %g\n",
4211                                //     iColumn,djValue,lower[iColumn],upper[iColumn]);
4212                                saveSolver->setColLower(iColumn, upper[iColumn]) ;
4213                                numberFixed++ ;
4214                            }
4215                        } else {
4216                            //printf("%d has dj of %g - already fixed to %g\n",
4217                            //     iColumn,djValue,lower[iColumn]);
4218                            numberFixed2++;
4219                        }
4220                    }
4221#ifdef COIN_DEVELOP
4222                    if ((specialOptions_&1) != 0) {
4223                        const OsiRowCutDebugger *debugger = saveSolver->getRowCutDebugger() ;
4224                        if (debugger) {
4225                            printf("Contains optimal\n") ;
4226                            OsiSolverInterface * temp = saveSolver->clone();
4227                            const double * solution = debugger->optimalSolution();
4228                            const double *lower = temp->getColLower() ;
4229                            const double *upper = temp->getColUpper() ;
4230                            int n = temp->getNumCols();
4231                            for (int i = 0; i < n; i++) {
4232                                if (temp->isInteger(i)) {
4233                                    double value = floor(solution[i] + 0.5);
4234                                    assert (value >= lower[i] && value <= upper[i]);
4235                                    temp->setColLower(i, value);
4236                                    temp->setColUpper(i, value);
4237                                }
4238                            }
4239                            temp->writeMps("reduced_fix");
4240                            delete temp;
4241                            saveSolver->writeMps("reduced");
4242                        } else {
4243                            abort();
4244                        }
4245                    }
4246                    printf("Restart could fix %d integers (%d already fixed)\n",
4247                           numberFixed + numberFixed2, numberFixed2);
4248#endif
4249                    numberFixed += numberFixed2;
4250                    if (numberFixed*10 < numberColumns && numberFixed*4 < numberIntegers_)
4251                        tryNewSearch = false;
4252                }
4253                if (tryNewSearch) {
4254                    // back to solver without cuts?
4255                    OsiSolverInterface * solver2 = saveSolver->clone();
4256                    const double *lower = saveSolver->getColLower() ;
4257                    const double *upper = saveSolver->getColUpper() ;
4258                    for (int i = 0 ; i < numberIntegers_ ; i++) {
4259                        int iColumn = integerVariable_[i] ;
4260                        solver2->setColLower(iColumn, lower[iColumn]);
4261                        solver2->setColUpper(iColumn, upper[iColumn]);
4262                    }
4263                    // swap
4264                    delete saveSolver;
4265                    saveSolver = solver2;
4266                    double * newSolution = new double[numberColumns];
4267                    double objectiveValue = checkCutoffForRestart;
4268                    // Save the best solution so far.
4269                    CbcSerendipity heuristic(*this);
4270                    if (bestSolution_)
4271                        heuristic.setInputSolution(bestSolution_, bestObjective_);
4272                    // Magic number
4273                    heuristic.setFractionSmall(0.8);
4274                    // `pumpTune' to stand-alone solver for explanations.
4275                    heuristic.setFeasibilityPumpOptions(1008013);
4276                    // Use numberNodes to say how many are original rows
4277                    heuristic.setNumberNodes(continuousSolver_->getNumRows());
4278#ifdef COIN_DEVELOP
4279                    if (continuousSolver_->getNumRows() <
4280                            solver_->getNumRows())
4281                        printf("%d rows added ZZZZZ\n",
4282                               solver_->getNumRows() - continuousSolver_->getNumRows());
4283#endif
4284                    int returnCode = heuristic.smallBranchAndBound(saveSolver,
4285                                     -1, newSolution,
4286                                     objectiveValue,
4287                                     checkCutoffForRestart, "Reduce");
4288                    if (returnCode < 0) {
4289#ifdef COIN_DEVELOP
4290                        printf("Restart - not small enough to do search after fixing\n");
4291#endif
4292                        delete [] newSolution;
4293                    } else {
4294                        // 1 for sol'n, 2 for finished, 3 for both
4295                        if ((returnCode&1) != 0) {
4296                            // increment number of solutions so other heuristics can test
4297                            numberSolutions_++;
4298                            numberHeuristicSolutions_++;
4299                            lastHeuristic_ = NULL;
4300                            setBestSolution(CBC_ROUNDING, objectiveValue, newSolution) ;
4301                        }
4302                        delete [] newSolution;
4303#ifdef CBC_THREAD
4304                        if (master_) {
4305                            lockThread();
4306                            if (parallelMode() > 0) {
4307                                while (master_->waitForThreadsInTree(0)) {
4308                                    lockThread();
4309                                    double dummyBest;
4310                                    tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4311                                    //unlockThread();
4312                                }
4313                            } else {
4314                                double dummyBest;
4315                                tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4316                            }
4317                            master_->waitForThreadsInTree(2);
4318                            delete master_;
4319                            master_ = NULL;
4320                            masterThread_ = NULL;
4321                        }
4322#endif
4323                        if (tree_->size()) {
4324                            double dummyBest;
4325                            tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4326                        }
4327                        break;
4328                    }
4329                }
4330                delete saveSolver;
4331                saveSolver = NULL;
4332            }
4333        }
4334        /*
4335          Check for abort on limits: node count, solution count, time, integrality gap.
4336        */
4337        if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
4338                numberSolutions_ < intParam_[CbcMaxNumSol] &&
4339                !maximumSecondsReached() &&
4340                !stoppedOnGap_ && !eventHappened_ && (maximumNumberIterations_ < 0 ||
4341                                                      numberIterations_ < maximumNumberIterations_))) {
4342            // out of loop
4343            break;
4344        }
4345#ifdef BONMIN
4346        assert(!solverCharacteristics_->solutionAddsCuts() || solverCharacteristics_->mipFeasible());
4347#endif
4348// Sets percentage of time when we try diving. Diving requires a bit of heap reorganisation, because
4349// we need to replace the comparison function to dive, and that requires reordering to retain the
4350// heap property.
4351#define DIVE_WHEN 1000
4352#define DIVE_STOP 2000
4353        int kNode = numberNodes_ % 4000;
4354        if (numberNodes_<100000 && kNode>DIVE_WHEN && kNode <= DIVE_STOP) {
4355            if (!parallelMode()) {
4356                if (kNode == DIVE_WHEN + 1 || numberConsecutiveInfeasible > 1) {
4357                    CbcCompareDefault * compare = dynamic_cast<CbcCompareDefault *>
4358                                                  (nodeCompare_);
4359                    // Don't interfere if user has replaced the compare function.
4360                    if (compare) {
4361                        //printf("Redoing tree\n");
4362                        compare->startDive(this);
4363                        numberConsecutiveInfeasible = 0;
4364                    }
4365                }
4366            }
4367        }
4368        // replace current cutoff?
4369        if (cutoff > getCutoff()) {
4370            double newCutoff = getCutoff();
4371            if (analyzeResults_) {
4372                // see if we could fix any (more)
4373                int n = 0;
4374                double * newLower = analyzeResults_;
4375                double * objLower = newLower + numberIntegers_;
4376                double * newUpper = objLower + numberIntegers_;
4377                double * objUpper = newUpper + numberIntegers_;
4378                for (int i = 0; i < numberIntegers_; i++) {
4379                    if (objLower[i] > newCutoff) {
4380                        n++;
4381                        if (objUpper[i] > newCutoff) {
4382                            newCutoff = -COIN_DBL_MAX;
4383                            break;
4384                        }
4385                    } else if (objUpper[i] > newCutoff) {
4386                        n++;
4387                    }
4388                }
4389                if (newCutoff == -COIN_DBL_MAX) {
4390                  COIN_DETAIL_PRINT(printf("Root analysis says finished\n"));
4391                } else if (n > numberFixedNow_) {
4392                  COIN_DETAIL_PRINT(printf("%d more fixed by analysis - now %d\n", n - numberFixedNow_, n));
4393                    numberFixedNow_ = n;
4394                }
4395            }
4396            if (eventHandler) {
4397                if (!eventHandler->event(CbcEventHandler::solution)) {
4398                    eventHappened_ = true; // exit
4399                }
4400                newCutoff = getCutoff();
4401            }
4402            lockThread();
4403            /*
4404              Clean the tree to reflect the new solution, then see if the
4405              node comparison predicate wants to make any changes. If so,
4406              call setComparison for the side effect of rebuilding the heap.
4407            */
4408            tree_->cleanTree(this,newCutoff,bestPossibleObjective_) ;
4409            if (nodeCompare_->newSolution(this) ||
4410                nodeCompare_->newSolution(this,continuousObjective_,
4411                                          continuousInfeasibilities_)) {
4412              tree_->setComparison(*nodeCompare_) ;
4413            }
4414            if (tree_->empty()) {
4415                continue;
4416            }
4417            unlockThread();
4418        }
4419        cutoff = getCutoff() ;
4420        /*
4421            Periodic activities: Opportunities to
4422            + tweak the nodeCompare criteria,
4423            + check if we've closed the integrality gap enough to quit,
4424            + print a summary line to let the user know we're working
4425        */
4426        if (numberNodes_ >= lastEvery1000) {
4427            lockThread();
4428#ifdef COIN_HAS_CLP
4429            // See if we want dantzig row choice
4430            goToDantzig(1000, savePivotMethod);
4431#endif
4432            lastEvery1000 = numberNodes_ + 1000;
4433            bool redoTree = nodeCompare_->every1000Nodes(this, numberNodes_) ;
4434#ifdef CHECK_CUT_SIZE
4435            verifyCutSize (tree_, *this);
4436#endif
4437            // redo tree if requested
4438            if (redoTree)
4439                tree_->setComparison(*nodeCompare_) ;
4440            unlockThread();
4441        }
4442        // Had hotstart before, now switched off
4443        if (saveCompare && !hotstartSolution_) {
4444            // hotstart switched off
4445            delete nodeCompare_; // off depth first
4446            nodeCompare_ = saveCompare;
4447            saveCompare = NULL;
4448            // redo tree
4449            lockThread();
4450            tree_->setComparison(*nodeCompare_) ;
4451            unlockThread();
4452        }
4453        if (numberNodes_ >= lastPrintEvery) {
4454            lastPrintEvery = numberNodes_ + printFrequency_;
4455            lockThread();
4456            int nNodes = tree_->size() ;
4457
4458            //MODIF PIERRE
4459            bestPossibleObjective_ = tree_->getBestPossibleObjective();
4460            unlockThread();
4461#ifdef CLP_INVESTIGATE
4462            if (getCutoff() < 1.0e20) {
4463                if (fabs(getCutoff() - (bestObjective_ - getCutoffIncrement())) > 1.0e-6 &&
4464                        !parentModel_)
4465                    printf("model cutoff in status %g, best %g, increment %g\n",
4466                           getCutoff(), bestObjective_, getCutoffIncrement());
4467                assert (getCutoff() < bestObjective_ - getCutoffIncrement() +
4468                        1.0e-6 + 1.0e-10*fabs(bestObjective_));
4469            }
4470#endif
4471            if (!intParam_[CbcPrinting]) {
4472                messageHandler()->message(CBC_STATUS, messages())
4473                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4474                << getCurrentSeconds()
4475                << CoinMessageEol ;
4476            } else if (intParam_[CbcPrinting] == 1) {
4477                messageHandler()->message(CBC_STATUS2, messages())
4478                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4479                << tree_->lastDepth() << tree_->lastUnsatisfied() 
4480                << tree_->lastObjective() << numberIterations_
4481                << getCurrentSeconds()
4482                << CoinMessageEol ;
4483            } else if (!numberExtraIterations_) {
4484                messageHandler()->message(CBC_STATUS2, messages())
4485                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4486                << tree_->lastDepth() << tree_->lastUnsatisfied() << numberIterations_
4487                << getCurrentSeconds()
4488                << CoinMessageEol ;
4489            } else {
4490                messageHandler()->message(CBC_STATUS3, messages())
4491                << numberNodes_ << numberExtraNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4492                << tree_->lastDepth() << tree_->lastUnsatisfied() << numberIterations_ << numberExtraIterations_
4493                << getCurrentSeconds()
4494                << CoinMessageEol ;
4495            }
4496            if (eventHandler && !eventHandler->event(CbcEventHandler::treeStatus)) {
4497                eventHappened_ = true; // exit
4498            }
4499        }
4500        // See if can stop on gap
4501        if(canStopOnGap()) {
4502            stoppedOnGap_ = true ;
4503        }
4504
4505#ifdef CHECK_NODE_FULL
4506        verifyTreeNodes(tree_, *this) ;
4507#   endif
4508#   ifdef CHECK_CUT_COUNTS
4509        verifyCutCounts(tree_, *this) ;
4510#   endif
4511        /*
4512          Now we come to the meat of the loop. To create the active subproblem, we'll
4513          pop the most promising node in the live set, rebuild the subproblem it
4514          represents, and then execute the current arm of the branch to create the
4515          active subproblem.
4516        */
4517        CbcNode * node = NULL;
4518#ifdef CBC_THREAD
4519        if (!parallelMode() || parallelMode() == -1) {
4520#endif
4521            node = tree_->bestNode(cutoff) ;
4522            // Possible one on tree worse than cutoff
4523            // Weird comparison function can leave ineligible nodes on tree
4524            if (!node || node->objectiveValue() > cutoff)
4525                continue;
4526            // Do main work of solving node here
4527            doOneNode(this, node, createdNode);
4528#ifdef JJF_ZERO
4529            if (node) {
4530              if (createdNode) {
4531                printf("Node %d depth %d, created %d depth %d\n",
4532                       node->nodeNumber(), node->depth(),
4533                       createdNode->nodeNumber(), createdNode->depth());
4534              } else {
4535                printf("Node %d depth %d,  no created node\n",
4536                       node->nodeNumber(), node->depth());
4537              }
4538            } else if (createdNode) {
4539              printf("Node exhausted, created %d depth %d\n",
4540                     createdNode->nodeNumber(), createdNode->depth());
4541            } else {
4542              printf("Node exhausted,  no created node\n");
4543              numberConsecutiveInfeasible = 2;
4544            }
4545#endif
4546            //if (createdNode)
4547            //numberConsecutiveInfeasible=0;
4548            //else
4549            //numberConsecutiveInfeasible++;
4550#ifdef CBC_THREAD
4551        } else if (parallelMode() > 0) {
4552            //lockThread();
4553            //node = tree_->bestNode(cutoff) ;
4554            // Possible one on tree worse than cutoff
4555            if (true || !node || node->objectiveValue() > cutoff) {
4556                assert (master_);
4557                if (master_) {
4558                    int anyLeft = master_->waitForThreadsInTree(1);
4559                    // may need to go round again
4560                    if (anyLeft) {
4561                        continue;
4562                    } else {
4563                        master_->stopThreads(-1);
4564                    }
4565                }
4566            }
4567            //unlockThread();
4568        } else {
4569            // Deterministic parallel
4570          if ((tree_->size() < CoinMax(numberThreads_, 8)||
4571               hotstartSolution_) && !goneParallel) {
4572                node = tree_->bestNode(cutoff) ;
4573                // Possible one on tree worse than cutoff
4574                if (!node || node->objectiveValue() > cutoff)
4575                    continue;
4576                // Do main work of solving node here
4577                doOneNode(this, node, createdNode);
4578                assert (createdNode);
4579                if (!createdNode->active()) {
4580                    delete createdNode;
4581                    createdNode = NULL;
4582                } else {
4583                    // Say one more pointing to this
4584                    node->nodeInfo()->increment() ;
4585                    tree_->push(createdNode) ;
4586                }
4587                if (node->active()) {
4588                    assert (node->nodeInfo());
4589                    if (node->nodeInfo()->numberBranchesLeft()) {
4590                        tree_->push(node) ;
4591                    } else {
4592                        node->setActive(false);
4593                    }
4594                } else {
4595                    if (node->nodeInfo()) {
4596                        if (!node->nodeInfo()->numberBranchesLeft())
4597                            node->nodeInfo()->allBranchesGone(); // can clean up
4598                        // So will delete underlying stuff
4599                        node->setActive(true);
4600                    }
4601                    delNode[nDeleteNode++] = node;
4602                    node = NULL;
4603                }
4604                if (nDeleteNode >= MAX_DEL_NODE) {
4605                    for (int i = 0; i < nDeleteNode; i++) {
4606                        //printf("trying to del %d %x\n",i,delNode[i]);
4607                        delete delNode[i];
4608                        //printf("done to del %d %x\n",i,delNode[i]);
4609                    }
4610                    nDeleteNode = 0;
4611                }
4612            } else {
4613                // Split and solve
4614                master_->deterministicParallel();
4615                goneParallel = true;
4616            }
4617        }
4618#endif
4619    }
4620    if (nDeleteNode) {
4621        for (int i = 0; i < nDeleteNode; i++) {
4622            delete delNode[i];
4623        }
4624        nDeleteNode = 0;
4625    }
4626#ifdef CBC_THREAD
4627    if (master_) {
4628        master_->stopThreads(-1);
4629        master_->waitForThreadsInTree(2);
4630        delete master_;
4631        master_ = NULL;
4632        masterThread_ = NULL;
4633        // adjust time to allow for children on some systems
4634        //dblParam_[CbcStartSeconds] -= CoinCpuTimeJustChildren();
4635    }
4636#endif
4637    /*
4638      End of the non-abort actions. The next block of code is executed if we've
4639      aborted because we hit one of the limits. Clean up by deleting the live set
4640      and break out of the node processing loop. Note that on an abort, node may
4641      have been pushed back onto the tree for further processing, in which case
4642      it'll be deleted in cleanTree. We need to check.
4643    */
4644    if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
4645          numberSolutions_ < intParam_[CbcMaxNumSol] &&
4646          !maximumSecondsReached() &&
4647          !stoppedOnGap_ &&
4648          !eventHappened_ &&
4649          (maximumNumberIterations_ < 0 || numberIterations_ < maximumNumberIterations_))
4650         ) {
4651        if (tree_->size()) {
4652            double dummyBest;
4653            tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4654        }
4655        delete nextRowCut_;
4656        /* order is important here:
4657         * maximumSecondsReached() should be checked before eventHappened_ and
4658         * isNodeLimitReached() should be checked after eventHappened_
4659         * reason is, that at timelimit, eventHappened_ is set to true to make Cbc stop fast
4660         *   and if Ctrl+C is hit, then the nodelimit is set to -1 to make Cbc stop
4661         */
4662        if (stoppedOnGap_) {
4663            messageHandler()->message(CBC_GAP, messages())
4664            << bestObjective_ - bestPossibleObjective_
4665            << dblParam_[CbcAllowableGap]
4666            << dblParam_[CbcAllowableFractionGap]*100.0
4667            << CoinMessageEol ;
4668            secondaryStatus_ = 2;
4669            status_ = 0 ;
4670        } else if (maximumSecondsReached()) {
4671            handler_->message(CBC_MAXTIME, messages_) << CoinMessageEol ;
4672            secondaryStatus_ = 4;
4673            status_ = 1 ;
4674        } else if (numberSolutions_ >= intParam_[CbcMaxNumSol]) {
4675            handler_->message(CBC_MAXSOLS, messages_) << CoinMessageEol ;
4676            secondaryStatus_ = 6;
4677            status_ = 1 ;
4678        } else if (isNodeLimitReached()) {
4679            handler_->message(CBC_MAXNODES, messages_) << CoinMessageEol ;
4680            secondaryStatus_ = 3;
4681            status_ = 1 ;
4682        } else if (maximumNumberIterations_ >= 0 && numberIterations_ >= maximumNumberIterations_) {
4683            handler_->message(CBC_MAXITERS, messages_) << CoinMessageEol ;
4684            secondaryStatus_ = 8;
4685            status_ = 1 ;
4686        } else {
4687            handler_->message(CBC_EVENT, messages_) << CoinMessageEol ;
4688            secondaryStatus_ = 5;
4689            status_ = 5 ;
4690        }
4691    }
4692    /*
4693      That's it, we've exhausted the search tree, or broken out of the loop because
4694      we hit some limit on evaluation.
4695
4696      We may have got an intelligent tree so give it one more chance
4697    */
4698    // Tell solver we are not in Branch and Cut
4699    solver_->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo, NULL) ;
4700    tree_->endSearch();
4701    //  If we did any sub trees - did we give up on any?
4702    if ( numberStoppedSubTrees_)
4703        status_ = 1;
4704    numberNodes_ += numberExtraNodes_;
4705    numberIterations_ += numberExtraIterations_;
4706    if (eventHandler) {
4707        eventHandler->event(CbcEventHandler::endSearch);
4708    }
4709    if (!status_) {
4710        // Set best possible unless stopped on gap
4711        if (secondaryStatus_ != 2)
4712            bestPossibleObjective_ = bestObjective_;
4713        handler_->message(CBC_END_GOOD, messages_)
4714        << bestObjective_ << numberIterations_ << numberNodes_ << getCurrentSeconds()
4715        << CoinMessageEol ;
4716    } else {
4717        handler_->message(CBC_END, messages_)
4718        << bestObjective_ << bestPossibleObjective_
4719        << numberIterations_ << numberNodes_ << getCurrentSeconds()
4720        << CoinMessageEol ;
4721    }
4722    if (numberStrongIterations_)
4723        handler_->message(CBC_STRONG_STATS, messages_)
4724        << strongInfo_[0] << numberStrongIterations_ << strongInfo_[2]
4725        << strongInfo_[1] << CoinMessageEol ;
4726    if (!numberExtraNodes_)
4727        handler_->message(CBC_OTHER_STATS, messages_)
4728        << maximumDepthActual_
4729        << numberDJFixed_ << CoinMessageEol ;
4730    else
4731        handler_->message(CBC_OTHER_STATS2, messages_)
4732        << maximumDepthActual_
4733        << numberDJFixed_ << numberExtraNodes_ << numberExtraIterations_
4734        << CoinMessageEol ;
4735    if (doStatistics == 100) {
4736        for (int i = 0; i < numberObjects_; i++) {
4737            CbcSimpleIntegerDynamicPseudoCost * obj =
4738                dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object_[i]) ;
4739            if (obj)
4740                obj->print();
4741        }
4742    }
4743    if (statistics_) {
4744        // report in some way
4745        int * lookup = new int[numberObjects_];
4746        int i;
4747        for (i = 0; i < numberObjects_; i++)
4748            lookup[i] = -1;
4749        bool goodIds = false; //true;
4750        for (i = 0; i < numberObjects_; i++) {
4751            int iColumn = object_[i]->columnNumber();
4752            if (iColumn >= 0 && iColumn < numberColumns) {
4753                if (lookup[i] == -1) {
4754                    lookup[i] = iColumn;
4755                } else {
4756                    goodIds = false;
4757                    break;
4758                }
4759            } else {
4760                goodIds = false;
4761                break;
4762            }
4763        }
4764        if (!goodIds) {
4765            delete [] lookup;
4766            lookup = NULL;
4767        }
4768        if (doStatistics >= 3) {
4769            printf("  node parent depth column   value                    obj      inf\n");
4770            for ( i = 0; i < numberNodes2_; i++) {
4771                statistics_[i]->print(lookup);
4772            }
4773        }
4774        if (doStatistics > 1) {
4775            // Find last solution
4776            int k;
4777            for (k = numberNodes2_ - 1; k >= 0; k--) {
4778                if (statistics_[k]->endingObjective() != COIN_DBL_MAX &&
4779                        !statistics_[k]->endingInfeasibility())
4780                    break;
4781            }
4782            if (k >= 0) {
4783                int depth = statistics_[k]->depth();
4784                int * which = new int[depth+1];
4785                for (i = depth; i >= 0; i--) {
4786                    which[i] = k;
4787                    k = statistics_[k]->parentNode();
4788                }
4789                printf("  node parent depth column   value                    obj      inf\n");
4790                for (i = 0; i <= depth; i++) {
4791                    statistics_[which[i]]->print(lookup);
4792                }
4793                delete [] which;
4794            }
4795        }
4796        // now summary
4797        int maxDepth = 0;
4798        double averageSolutionDepth = 0.0;
4799        int numberSolutions = 0;
4800        double averageCutoffDepth = 0.0;
4801        double averageSolvedDepth = 0.0;
4802        int numberCutoff = 0;
4803        int numberDown = 0;
4804        int numberFirstDown = 0;
4805        double averageInfDown = 0.0;
4806        double averageObjDown = 0.0;
4807        int numberCutoffDown = 0;
4808        int numberUp = 0;
4809        int numberFirstUp = 0;
4810        double averageInfUp = 0.0;
4811        double averageObjUp = 0.0;
4812        int numberCutoffUp = 0;
4813        double averageNumberIterations1 = 0.0;
4814        double averageValue = 0.0;
4815        for ( i = 0; i < numberNodes2_; i++) {
4816            int depth =  statistics_[i]->depth();
4817            int way =  statistics_[i]->way();
4818            double value = statistics_[i]->value();
4819            double startingObjective =  statistics_[i]->startingObjective();
4820            int startingInfeasibility = statistics_[i]->startingInfeasibility();
4821            double endingObjective = statistics_[i]->endingObjective();
4822            int endingInfeasibility = statistics_[i]->endingInfeasibility();
4823            maxDepth = CoinMax(depth, maxDepth);
4824            // Only for completed
4825            averageNumberIterations1 += statistics_[i]->numberIterations();
4826            averageValue += value;
4827            if (endingObjective != COIN_DBL_MAX && !endingInfeasibility) {
4828                numberSolutions++;
4829                averageSolutionDepth += depth;
4830            }
4831            if (endingObjective == COIN_DBL_MAX) {
4832                numberCutoff++;
4833                averageCutoffDepth += depth;
4834                if (way < 0) {
4835                    numberDown++;
4836                    numberCutoffDown++;
4837                    if (way == -1)
4838                        numberFirstDown++;
4839                } else {
4840                    numberUp++;
4841                    numberCutoffUp++;
4842                    if (way == 1)
4843                        numberFirstUp++;
4844                }
4845            } else {
4846                averageSolvedDepth += depth;
4847                if (way < 0) {
4848                    numberDown++;
4849                    averageInfDown += startingInfeasibility - endingInfeasibility;
4850                    averageObjDown += endingObjective - startingObjective;
4851                    if (way == -1)
4852                        numberFirstDown++;
4853                } else {
4854                    numberUp++;
4855                    averageInfUp += startingInfeasibility - endingInfeasibility;
4856                    averageObjUp += endingObjective - startingObjective;
4857                    if (way == 1)
4858                        numberFirstUp++;
4859                }
4860            }
4861        }
4862        // Now print
4863        if (numberSolutions)
4864            averageSolutionDepth /= static_cast<double> (numberSolutions);
4865        int numberSolved = numberNodes2_ - numberCutoff;
4866        double averageNumberIterations2 = numberIterations_ - averageNumberIterations1
4867                                          - numberIterationsAtContinuous;
4868        if (numberCutoff) {
4869            averageCutoffDepth /= static_cast<double> (numberCutoff);
4870            averageNumberIterations2 /= static_cast<double> (numberCutoff);
4871        }
4872        if (numberNodes2_)
4873            averageValue /= static_cast<double> (numberNodes2_);
4874        if (numberSolved) {
4875            averageNumberIterations1 /= static_cast<double> (numberSolved);
4876            averageSolvedDepth /= static_cast<double> (numberSolved);
4877        }
4878        printf("%d solution(s) were found (by branching) at an average depth of %g\n",
4879               numberSolutions, averageSolutionDepth);
4880        printf("average value of variable being branched on was %g\n",
4881               averageValue);
4882        printf("%d nodes were cutoff at an average depth of %g with iteration count of %g\n",
4883               numberCutoff, averageCutoffDepth, averageNumberIterations2);
4884        printf("%d nodes were solved at an average depth of %g with iteration count of %g\n",
4885               numberSolved, averageSolvedDepth, averageNumberIterations1);
4886        if (numberDown) {
4887            averageInfDown /= static_cast<double> (numberDown);
4888            averageObjDown /= static_cast<double> (numberDown);
4889        }
4890        printf("Down %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
4891               numberDown, numberFirstDown, numberDown - numberFirstDown, numberCutoffDown,
4892               averageInfDown, averageObjDown);
4893        if (numberUp) {
4894            averageInfUp /= static_cast<double> (numberUp);
4895            averageObjUp /= static_cast<double> (numberUp);
4896        }
4897        printf("Up %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
4898               numberUp, numberFirstUp, numberUp - numberFirstUp, numberCutoffUp,
4899               averageInfUp, averageObjUp);
4900        for ( i = 0; i < numberNodes2_; i++)
4901            delete statistics_[i];
4902        delete [] statistics_;
4903        statistics_ = NULL;
4904        maximumStatistics_ = 0;
4905        delete [] lookup;
4906    }
4907    /*
4908      If we think we have a solution, restore and confirm it with a call to
4909      setBestSolution().  We need to reset the cutoff value so as not to fathom
4910      the solution on bounds.  Note that calling setBestSolution( ..., true)
4911      leaves the continuousSolver_ bounds vectors fixed at the solution value.
4912
4913      Running resolve() here is a failsafe --- setBestSolution has already
4914      reoptimised using the continuousSolver_. If for some reason we fail to
4915      prove optimality, run the problem again after instructing the solver to
4916      tell us more.
4917
4918      If all looks good, replace solver_ with continuousSolver_, so that the
4919      outside world will be able to obtain information about the solution using
4920      public methods.
4921
4922      Don't replace if we are trying to save cuts
4923    */
4924    if (bestSolution_ && (solverCharacteristics_->solverType() < 2 || solverCharacteristics_->solverType() == 4) && 
4925        ((specialOptions_&8388608)==0||(specialOptions_&2048)!=0)) {
4926        setCutoff(1.0e50) ; // As best solution should be worse than cutoff
4927        // change cutoff as constraint if wanted
4928        if (cutoffRowNumber_>=0) {
4929          if (solver_->getNumRows()>cutoffRowNumber_)
4930            solver_->setRowUpper(cutoffRowNumber_,1.0e50);
4931        }
4932        // also in continuousSolver_
4933        if (continuousSolver_) {
4934          // Solvers know about direction
4935          double direction = solver_->getObjSense();
4936          continuousSolver_->setDblParam(OsiDualObjectiveLimit, 1.0e50*direction);
4937        }
4938        phase_ = 5;
4939        double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
4940        if ((specialOptions_&4) == 0)
4941            bestObjective_ += 100.0 * increment + 1.0e-3; // only set if we are going to solve
4942        setBestSolution(CBC_END_SOLUTION, bestObjective_, bestSolution_, 1) ;
4943        continuousSolver_->resolve() ;
4944        if (!continuousSolver_->isProvenOptimal()) {
4945            continuousSolver_->messageHandler()->setLogLevel(2) ;
4946            continuousSolver_->initialSolve() ;
4947        }
4948        delete solver_ ;
4949        // above deletes solverCharacteristics_
4950        solverCharacteristics_ = NULL;
4951        solver_ = continuousSolver_ ;
4952        setPointers(solver_);
4953        continuousSolver_ = NULL ;
4954    }
4955    /*
4956      Clean up dangling objects. continuousSolver_ may already be toast.
4957    */
4958    delete lastws ;
4959    if (saveObjects) {
4960        for (int i = 0; i < numberObjects_; i++)
4961            delete saveObjects[i];
4962        delete [] saveObjects;
4963    }
4964    numberStrong_ = saveNumberStrong;
4965    numberBeforeTrust_ = saveNumberBeforeTrust;
4966    delete [] whichGenerator_ ;
4967    whichGenerator_ = NULL;
4968    delete [] lowerBefore ;
4969    delete [] upperBefore ;
4970    delete [] walkback_ ;
4971    walkback_ = NULL ;
4972    delete [] lastNodeInfo_ ;
4973    lastNodeInfo_ = NULL;
4974    delete [] lastNumberCuts_ ;
4975    lastNumberCuts_ = NULL;
4976    delete [] lastCut_;
4977    lastCut_ = NULL;
4978    delete [] addedCuts_ ;
4979    addedCuts_ = NULL ;
4980    //delete persistentInfo;
4981    // Get rid of characteristics
4982    solverCharacteristics_ = NULL;
4983    if (continuousSolver_) {
4984        delete continuousSolver_ ;
4985        continuousSolver_ = NULL ;
4986    }
4987    /*
4988      Destroy global cuts by replacing with an empty OsiCuts object.
4989    */
4990    globalCuts_ = CbcRowCuts() ;
4991    delete globalConflictCuts_;
4992    globalConflictCuts_=NULL;
4993    if (!bestSolution_ && (specialOptions_&8388608)==0) {
4994        // make sure lp solver is infeasible
4995        int numberColumns = solver_->getNumCols();
4996        const double * columnLower = solver_->getColLower();
4997        int iColumn;
4998        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4999            if (solver_->isInteger(iColumn))
5000                solver_->setColUpper(iColumn, columnLower[iColumn]);
5001        }
5002        solver_->initialSolve();
5003    }
5004#ifdef COIN_HAS_CLP
5005    {
5006        OsiClpSolverInterface * clpSolver
5007        = dynamic_cast<OsiClpSolverInterface *> (solver_);
5008        if (clpSolver) {
5009            // Possible restore of pivot method
5010            if (savePivotMethod) {
5011                // model may have changed
5012                savePivotMethod->setModel(NULL);
5013                clpSolver->getModelPtr()->setDualRowPivotAlgorithm(*savePivotMethod);
5014                delete savePivotMethod;
5015            }
5016            clpSolver->setLargestAway(-1.0);
5017        }
5018    }
5019#endif
5020    if ((fastNodeDepth_ >= 1000000 || (moreSpecialOptions_&33554432)!=0)
5021         && !parentModel_) {
5022      // delete object off end
5023      delete object_[numberObjects_];
5024      if ((moreSpecialOptions_&33554432)==0)
5025        fastNodeDepth_ -= 1000000;
5026    }
5027    delete saveSolver;
5028    // Undo preprocessing performed during BaB.
5029    if (strategy_ && strategy_->preProcessState() > 0) {
5030        // undo preprocessing
5031        CglPreProcess * process = strategy_->process();
5032        assert (process);
5033        int n = originalSolver->getNumCols();
5034        if (bestSolution_) {
5035            delete [] bestSolution_;
5036            bestSolution_ = new double [n];
5037            process->postProcess(*solver_);
5038        }
5039        strategy_->deletePreProcess();
5040        // Solution now back in originalSolver
5041        delete solver_;
5042        solver_ = originalSolver;
5043        if (bestSolution_) {
5044            bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
5045            memcpy(bestSolution_, solver_->getColSolution(), n*sizeof(double));
5046        }
5047        // put back original objects if there were any
5048        if (originalObject) {
5049            int iColumn;
5050            assert (ownObjects_);
5051            for (iColumn = 0; iColumn < numberObjects_; iColumn++)
5052                delete object_[iColumn];
5053            delete [] object_;
5054            numberObjects_ = numberOriginalObjects;
5055            object_ = originalObject;
5056            delete [] integerVariable_;
5057            numberIntegers_ = 0;
5058            for (iColumn = 0; iColumn < n; iColumn++) {
5059                if (solver_->isInteger(iColumn))
5060                    numberIntegers_++;
5061            }
5062            integerVariable_ = new int[numberIntegers_];
5063            numberIntegers_ = 0;
5064            for (iColumn = 0; iColumn < n; iColumn++) {
5065                if (solver_->isInteger(iColumn))
5066                    integerVariable_[numberIntegers_++] = iColumn;
5067            }
5068        }
5069    }
5070    if (flipObjective)
5071      flipModel();
5072#ifdef COIN_HAS_CLP
5073    {
5074        OsiClpSolverInterface * clpSolver
5075        = dynamic_cast<OsiClpSolverInterface *> (solver_);
5076        if (clpSolver)
5077            clpSolver->setFakeObjective(reinterpret_cast<double *> (NULL));
5078    }
5079#endif
5080    moreSpecialOptions_ = saveMoreSpecialOptions;
5081    return ;
5082}
5083
5084
5085// Solve the initial LP relaxation
5086void
5087CbcModel::initialSolve()
5088{
5089    assert (solver_);
5090    // Double check optimization directions line up
5091    dblParam_[CbcOptimizationDirection] = solver_->getObjSense();
5092    // Check if bounds are all integral (as may get messed up later)
5093    checkModel();
5094    if (!solverCharacteristics_) {
5095        OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
5096        if (solverCharacteristics) {
5097            solverCharacteristics_ = solverCharacteristics;
5098        } else {
5099            // replace in solver
5100            OsiBabSolver defaultC;
5101            solver_->setAuxiliaryInfo(&defaultC);
5102            solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
5103        }
5104    }
5105    solverCharacteristics_->setSolver(solver_);
5106    solver_->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
5107    solver_->initialSolve();
5108    solver_->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo, NULL) ;
5109    if (!solver_->isProvenOptimal())
5110        solver_->resolve();
5111    // But set up so Jon Lee will be happy
5112    status_ = -1;
5113    secondaryStatus_ = -1;
5114    originalContinuousObjective_ = solver_->getObjValue() * solver_->getObjSense();
5115    bestPossibleObjective_ = originalContinuousObjective_;
5116    if (solver_->isProvenDualInfeasible())
5117      originalContinuousObjective_ = -COIN_DBL_MAX;
5118    delete [] continuousSolution_;
5119    continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
5120                                          solver_->getNumCols());
5121    setPointers(solver_);
5122    solverCharacteristics_ = NULL;
5123}
5124
5125/*! \brief Get an empty basis object
5126
5127  Return an empty CoinWarmStartBasis object with the requested capacity,
5128  appropriate for the current solver. The object is cloned from the object
5129  cached as emptyWarmStart_. If there is no cached object, the routine
5130  queries the solver for a warm start object, empties it, and caches the
5131  result.
5132*/
5133
5134CoinWarmStartBasis *CbcModel::getEmptyBasis (int ns, int na) const
5135
5136{
5137    CoinWarmStartBasis *emptyBasis ;
5138    /*
5139      Acquire an empty basis object, if we don't yet have one.
5140    */
5141    if (emptyWarmStart_ == 0) {
5142        if (solver_ == 0) {
5143            throw CoinError("Cannot construct basis without solver!",
5144                            "getEmptyBasis", "CbcModel") ;
5145        }
5146        emptyBasis =
5147            dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart()) ;
5148        if (emptyBasis == 0) {
5149            throw CoinError(
5150                "Solver does not appear to use a basis-oriented warm start.",
5151                "getEmptyBasis", "CbcModel") ;
5152        }
5153        emptyBasis->setSize(0, 0) ;
5154        emptyWarmStart_ = dynamic_cast<CoinWarmStart *>(emptyBasis) ;
5155    }
5156    /*
5157      Clone the empty basis object, resize it as requested, and return.
5158    */
5159    emptyBasis = dynamic_cast<CoinWarmStartBasis *>(emptyWarmStart_->clone()) ;
5160    assert(emptyBasis) ;
5161    if (ns != 0 || na != 0) emptyBasis->setSize(ns, na) ;
5162
5163    return (emptyBasis) ;
5164}
5165
5166
5167/** Default Constructor
5168
5169  Creates an empty model without an associated solver.
5170*/
5171CbcModel::CbcModel()
5172
5173        :
5174        solver_(NULL),
5175        ownership_(0x80000000),
5176        continuousSolver_(NULL),
5177        referenceSolver_(NULL),
5178        defaultHandler_(true),
5179        emptyWarmStart_(NULL),
5180        bestObjective_(COIN_DBL_MAX),
5181        bestPossibleObjective_(COIN_DBL_MAX),
5182        sumChangeObjective1_(0.0),
5183        sumChangeObjective2_(0.0),
5184        bestSolution_(NULL),
5185        savedSolutions_(NULL),
5186        currentSolution_(NULL),
5187        testSolution_(NULL),
5188        globalConflictCuts_(NULL),
5189        minimumDrop_(1.0e-4),
5190        numberSolutions_(0),
5191        numberSavedSolutions_(0),
5192        maximumSavedSolutions_(0),
5193        stateOfSearch_(0),
5194        whenCuts_(-1),
5195        hotstartSolution_(NULL),
5196        hotstartPriorities_(NULL),
5197        numberHeuristicSolutions_(0),
5198        numberNodes_(0),
5199        numberNodes2_(0),
5200        numberIterations_(0),
5201        numberSolves_(0),
5202        status_(-1),
5203        secondaryStatus_(-1),
5204        numberIntegers_(0),
5205        numberRowsAtContinuous_(0),
5206        cutoffRowNumber_(-1),
5207        maximumNumberCuts_(0),
5208        phase_(0),
5209        currentNumberCuts_(0),
5210        maximumDepth_(0),
5211        walkback_(NULL),
5212        lastNodeInfo_(NULL),
5213        lastCut_(NULL),
5214        lastDepth_(0),
5215        lastNumberCuts2_(0),
5216        maximumCuts_(0),
5217        lastNumberCuts_(NULL),
5218        addedCuts_(NULL),
5219        nextRowCut_(NULL),
5220        currentNode_(NULL),
5221        integerVariable_(NULL),
5222        integerInfo_(NULL),
5223        continuousSolution_(NULL),
5224        usedInSolution_(NULL),
5225        specialOptions_(0),
5226        moreSpecialOptions_(0),
5227        moreSpecialOptions2_(0),
5228        topOfTree_(NULL),
5229        subTreeModel_(NULL),
5230        heuristicModel_(NULL),
5231        numberStoppedSubTrees_(0),
5232        presolve_(0),
5233        numberStrong_(5),
5234        numberBeforeTrust_(10),
5235        numberPenalties_(20),
5236        stopNumberIterations_(-1),
5237        penaltyScaleFactor_(3.0),
5238        numberAnalyzeIterations_(0),
5239        analyzeResults_(NULL),
5240        numberInfeasibleNodes_(0),
5241        problemType_(0),
5242        printFrequency_(0),
5243        numberCutGenerators_(0),
5244        generator_(NULL),
5245        virginGenerator_(NULL),
5246