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

Last change on this file since 2043 was 2043, checked in by forrest, 5 years ago

minor changes for SOS etc

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 740.0 KB
Line 
1/* $Id: CbcModel.cpp 2043 2014-07-01 13:03:51Z 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    {
2008      // check
2009      int numberOdd = 0;
2010      for (int i = 0; i < numberObjects_; i++) {
2011        CbcSimpleInteger * obj =
2012          dynamic_cast <CbcSimpleInteger *>(object_[i]) ;
2013        if (!obj)
2014          numberOdd++;
2015      }
2016      if (numberOdd)
2017        moreSpecialOptions_ |= 1073741824;
2018    }
2019    // If NLP then we assume already solved outside branchAndbound
2020    if (!solverCharacteristics_->solverType() || solverCharacteristics_->solverType() == 4) {
2021        feasible = resolve(NULL, 0) != 0 ;
2022    } else {
2023        // pick up given status
2024        feasible = (solver_->isProvenOptimal() &&
2025                    !solver_->isDualObjectiveLimitReached()) ;
2026    }
2027    if (problemFeasibility_->feasible(this, 0) < 0) {
2028        feasible = false; // pretend infeasible
2029    }
2030    numberSavedSolutions_ = 0;
2031    int saveNumberStrong = numberStrong_;
2032    int saveNumberBeforeTrust = numberBeforeTrust_;
2033    /*
2034      If the linear relaxation of the root is infeasible, bail out now. Otherwise,
2035      continue with processing the root node.
2036    */
2037    if (!feasible) {
2038        status_ = 0 ;
2039        if (!solver_->isProvenDualInfeasible()) {
2040            handler_->message(CBC_INFEAS, messages_) << CoinMessageEol ;
2041            secondaryStatus_ = 1;
2042        } else {
2043            handler_->message(CBC_UNBOUNDED, messages_) << CoinMessageEol ;
2044            secondaryStatus_ = 7;
2045        }
2046        originalContinuousObjective_ = COIN_DBL_MAX;
2047        solverCharacteristics_ = NULL;
2048        if (flipObjective)
2049          flipModel();
2050        return ;
2051    } else if (!numberObjects_) { 
2052        // nothing to do
2053        // Undo preprocessing performed during BaB.
2054        if (strategy_ && strategy_->preProcessState() > 0) {
2055          // undo preprocessing
2056          CglPreProcess * process = strategy_->process();
2057          assert (process);
2058          int n = originalSolver->getNumCols();
2059          if (bestSolution_) {
2060            delete [] bestSolution_;
2061            bestSolution_ = new double [n];
2062            process->postProcess(*solver_);
2063          }
2064          strategy_->deletePreProcess();
2065          // Solution now back in originalSolver
2066          delete solver_;
2067          solver_ = originalSolver;
2068          if (bestSolution_) {
2069            bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
2070            memcpy(bestSolution_, solver_->getColSolution(), n*sizeof(double));
2071          }
2072          // put back original objects if there were any
2073          if (originalObject) {
2074            int iColumn;
2075            assert (ownObjects_);
2076            for (iColumn = 0; iColumn < numberObjects_; iColumn++)
2077              delete object_[iColumn];
2078            delete [] object_;
2079            numberObjects_ = numberOriginalObjects;
2080            object_ = originalObject;
2081            delete [] integerVariable_;
2082            numberIntegers_ = 0;
2083            for (iColumn = 0; iColumn < n; iColumn++) {
2084              if (solver_->isInteger(iColumn))
2085                numberIntegers_++;
2086            }
2087            integerVariable_ = new int[numberIntegers_];
2088            numberIntegers_ = 0;
2089            for (iColumn = 0; iColumn < n; iColumn++) {
2090            if (solver_->isInteger(iColumn))
2091              integerVariable_[numberIntegers_++] = iColumn;
2092            }
2093          }
2094        }
2095        if (flipObjective)
2096          flipModel();
2097        solverCharacteristics_ = NULL;
2098        bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
2099        int numberColumns = solver_->getNumCols();
2100        delete [] bestSolution_;
2101        bestSolution_ = new double[numberColumns];
2102        CoinCopyN(solver_->getColSolution(), numberColumns, bestSolution_);
2103        return ;
2104    }
2105    /*
2106      See if we're using the Osi side of the branching hierarchy. If so, either
2107      convert existing CbcObjects to OsiObjects, or generate them fresh. In the
2108      first case, CbcModel owns the objects on the object_ list. In the second
2109      case, the solver holds the objects and object_ simply points to the
2110      solver's list.
2111
2112      080417 The conversion code here (the block protected by `if (obj)') cannot
2113      possibly be correct. On the Osi side, descent is OsiObject -> OsiObject2 ->
2114      all other Osi object classes. On the Cbc side, it's OsiObject -> CbcObject
2115      -> all other Cbc object classes. It's structurally impossible for any Osi
2116      object to descend from CbcObject. The only thing I can see is that this is
2117      really dead code, and object detection is now handled from the Osi side.
2118    */
2119    // Convert to Osi if wanted
2120    //OsiBranchingInformation * persistentInfo = NULL;
2121    if (branchingMethod_ && branchingMethod_->chooseMethod()) {
2122        //persistentInfo = new OsiBranchingInformation(solver_);
2123        if (numberOriginalObjects) {
2124            for (int iObject = 0 ; iObject < numberObjects_ ; iObject++) {
2125                CbcObject * obj =
2126                    dynamic_cast <CbcObject *>(object_[iObject]) ;
2127                if (obj) {
2128                    CbcSimpleInteger * obj2 =
2129                        dynamic_cast <CbcSimpleInteger *>(obj) ;
2130                    if (obj2) {
2131                        // back to Osi land
2132                        object_[iObject] = obj2->osiObject();
2133                        delete obj;
2134                    } else {
2135                        OsiSimpleInteger * obj3 =
2136                            dynamic_cast <OsiSimpleInteger *>(obj) ;
2137                        if (!obj3) {
2138                            OsiSOS * obj4 =
2139                                dynamic_cast <OsiSOS *>(obj) ;
2140                            if (!obj4) {
2141                                CbcSOS * obj5 =
2142                                    dynamic_cast <CbcSOS *>(obj) ;
2143                                if (obj5) {
2144                                    // back to Osi land
2145                                    object_[iObject] = obj5->osiObject(solver_);
2146                                } else {
2147                                    printf("Code up CbcObject type in Osi land\n");
2148                                    abort();
2149                                }
2150                            }
2151                        }
2152                    }
2153                }
2154            }
2155            // and add to solver
2156            //if (!solver_->numberObjects()) {
2157            solver_->addObjects(numberObjects_, object_);
2158            //} else {
2159            //if (solver_->numberObjects()!=numberOriginalObjects) {
2160            //printf("should have trapped that solver has objects before\n");
2161            //abort();
2162            //}
2163            //}
2164        } else {
2165            /*
2166              As of 080104, findIntegersAndSOS is misleading --- the default OSI
2167              implementation finds only integers.
2168            */
2169            // do from solver
2170            deleteObjects(false);
2171            solver_->findIntegersAndSOS(false);
2172            numberObjects_ = solver_->numberObjects();
2173            object_ = solver_->objects();
2174            ownObjects_ = false;
2175        }
2176        branchingMethod_->chooseMethod()->setSolver(solver_);
2177    }
2178    // take off heuristics if have to (some do not work with SOS, for example)
2179    // object should know what's safe.
2180    {
2181        int numberOdd = 0;
2182        int numberSOS = 0;
2183        for (int i = 0; i < numberObjects_; i++) {
2184            if (!object_[i]->canDoHeuristics())
2185                numberOdd++;
2186            CbcSOS * obj =
2187                dynamic_cast <CbcSOS *>(object_[i]) ;
2188            if (obj)
2189                numberSOS++;
2190        }
2191        if (numberOdd) {
2192          if (numberHeuristics_ && (specialOptions_&1024)==0 ) {
2193                int k = 0;
2194                for (int i = 0; i < numberHeuristics_; i++) {
2195                    if (!heuristic_[i]->canDealWithOdd())
2196                        delete heuristic_[i];
2197                    else
2198                        heuristic_[k++] = heuristic_[i];
2199                }
2200                if (!k) {
2201                    delete [] heuristic_;
2202                    heuristic_ = NULL;
2203                }
2204                numberHeuristics_ = k;
2205                handler_->message(CBC_HEURISTICS_OFF, messages_) << numberOdd << CoinMessageEol ;
2206            }
2207            // If odd switch off AddIntegers
2208            specialOptions_ &= ~65536;
2209        } else if (numberSOS) {
2210            specialOptions_ |= 128; // say can do SOS in dynamic mode
2211            // switch off fast nodes for now
2212            fastNodeDepth_ = -1;
2213            moreSpecialOptions_ &= ~33554432; // no diving
2214        }
2215        if (numberThreads_ > 0) {
2216            // switch off fast nodes for now
2217            fastNodeDepth_ = -1;
2218        }
2219    }
2220    // Save objective (just so user can access it)
2221    originalContinuousObjective_ = solver_->getObjValue()* solver_->getObjSense();
2222    bestPossibleObjective_ = originalContinuousObjective_;
2223    sumChangeObjective1_ = 0.0;
2224    sumChangeObjective2_ = 0.0;
2225    /*
2226      OsiRowCutDebugger knows an optimal answer for a subset of MIP problems.
2227      Assuming it recognises the problem, when called upon it will check a cut to
2228      see if it cuts off the optimal answer.
2229    */
2230    // If debugger exists set specialOptions_ bit
2231    if (solver_->getRowCutDebuggerAlways()) {
2232        specialOptions_ |= 1;
2233    }
2234
2235# ifdef CBC_DEBUG
2236    if ((specialOptions_&1) == 0)
2237        solver_->activateRowCutDebugger(problemName.c_str()) ;
2238    if (solver_->getRowCutDebuggerAlways())
2239        specialOptions_ |= 1;
2240# endif
2241
2242    /*
2243      Begin setup to process a feasible root node.
2244    */
2245    bestObjective_ = CoinMin(bestObjective_, 1.0e50) ;
2246    if (!bestSolution_) {
2247        numberSolutions_ = 0 ;
2248        numberHeuristicSolutions_ = 0 ;
2249    }
2250    stateOfSearch_ = 0;
2251    // Everything is minimization
2252    {
2253        // needed to sync cutoffs
2254        double value ;
2255        solver_->getDblParam(OsiDualObjectiveLimit, value) ;
2256        dblParam_[CbcCurrentCutoff] = value * solver_->getObjSense();
2257    }
2258    double cutoff = getCutoff() ;
2259    double direction = solver_->getObjSense() ;
2260    dblParam_[CbcOptimizationDirection] = direction;
2261    if (cutoff < 1.0e20 && direction < 0.0)
2262        messageHandler()->message(CBC_CUTOFF_WARNING1,
2263                                  messages())
2264        << cutoff << -cutoff << CoinMessageEol ;
2265    if (cutoff > bestObjective_)
2266        cutoff = bestObjective_ ;
2267    setCutoff(cutoff) ;
2268    /*
2269      We probably already have a current solution, but just in case ...
2270    */
2271    int numberColumns = getNumCols() ;
2272    if (!currentSolution_)
2273        currentSolution_ = new double[numberColumns] ;
2274    testSolution_ = currentSolution_;
2275    /*
2276      Create a copy of the solver, thus capturing the original (root node)
2277      constraint system (aka the continuous system).
2278    */
2279    continuousSolver_ = solver_->clone() ;
2280
2281    // add cutoff as constraint if wanted
2282    if (cutoffRowNumber_==-2) {
2283      if (!parentModel_) {
2284        int numberColumns=solver_->getNumCols();
2285        double * obj = CoinCopyOfArray(solver_->getObjCoefficients(),numberColumns);
2286        int * indices = new int [numberColumns];
2287        int n=0;
2288        for (int i=0;i<numberColumns;i++) {
2289          if (obj[i]) {
2290            indices[n]=i;
2291            obj[n++]=obj[i];
2292          }
2293        }
2294        if (n) {
2295          double cutoff=getCutoff();
2296          // relax a little bit
2297          cutoff += 1.0e-4;
2298          double offset;
2299          solver_->getDblParam(OsiObjOffset, offset);
2300          cutoffRowNumber_ = solver_->getNumRows();
2301          solver_->addRow(n,indices,obj,-COIN_DBL_MAX,CoinMin(cutoff,1.0e25)+offset);
2302        } else {
2303          // no objective!
2304          cutoffRowNumber_ = -1;
2305        }
2306        delete [] indices;
2307        delete [] obj;
2308      } else {
2309        // switch off
2310        cutoffRowNumber_ = -1;
2311      }
2312    }
2313    numberRowsAtContinuous_ = getNumRows() ;
2314    solver_->saveBaseModel();
2315    /*
2316      Check the objective to see if we can deduce a nontrivial increment. If
2317      it's better than the current value for CbcCutoffIncrement, it'll be
2318      installed.
2319    */
2320    if (solverCharacteristics_->reducedCostsAccurate())
2321        analyzeObjective() ;
2322    {
2323        // may be able to change cutoff now
2324        double cutoff = getCutoff();
2325        double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
2326        if (cutoff > bestObjective_ - increment) {
2327            cutoff = bestObjective_ - increment ;
2328            setCutoff(cutoff) ;
2329        }
2330    }
2331#ifdef COIN_HAS_CLP
2332    // Possible save of pivot method
2333    ClpDualRowPivot * savePivotMethod = NULL;
2334    {
2335        // pass tolerance and increment to solver
2336        OsiClpSolverInterface * clpSolver
2337        = dynamic_cast<OsiClpSolverInterface *> (solver_);
2338        if (clpSolver)
2339            clpSolver->setStuff(getIntegerTolerance(), getCutoffIncrement());
2340#ifdef CLP_RESOLVE
2341        if ((moreSpecialOptions_&1048576)!=0&&!parentModel_&&clpSolver) {
2342          resolveClp(clpSolver,0);
2343        }
2344#endif
2345    }
2346#endif
2347    /*
2348      Set up for cut generation. addedCuts_ holds the cuts which are relevant for
2349      the active subproblem. whichGenerator will be used to record the generator
2350      that produced a given cut.
2351    */
2352    maximumWhich_ = 1000 ;
2353    delete [] whichGenerator_;
2354    whichGenerator_ = new int[maximumWhich_] ;
2355    memset(whichGenerator_, 0, maximumWhich_*sizeof(int));
2356    maximumNumberCuts_ = 0 ;
2357    currentNumberCuts_ = 0 ;
2358    delete [] addedCuts_ ;
2359    addedCuts_ = NULL ;
2360    OsiObject ** saveObjects = NULL;
2361    maximumRows_ = numberRowsAtContinuous_;
2362    currentDepth_ = 0;
2363    workingBasis_.resize(maximumRows_, numberColumns);
2364    /*
2365      Set up an empty heap and associated data structures to hold the live set
2366      (problems which require further exploration).
2367    */
2368    CbcCompareDefault * compareActual
2369    = dynamic_cast<CbcCompareDefault *> (nodeCompare_);
2370    if (compareActual) {
2371        compareActual->setBestPossible(direction*solver_->getObjValue());
2372        compareActual->setCutoff(getCutoff());
2373#ifdef JJF_ZERO
2374        if (false && !numberThreads_ && !parentModel_) {
2375            printf("CbcTreeArray ? threads ? parentArray\n");
2376            // Setup new style tree
2377            delete tree_;
2378            tree_ = new CbcTreeArray();
2379        }
2380#endif
2381    }
2382    tree_->setComparison(*nodeCompare_) ;
2383    /*
2384      Used to record the path from a node to the root of the search tree, so that
2385      we can then traverse from the root to the node when restoring a subproblem.
2386    */
2387    maximumDepth_ = 10 ;
2388    delete [] walkback_ ;
2389    walkback_ = new CbcNodeInfo * [maximumDepth_] ;
2390    lastDepth_ = 0;
2391    delete [] lastNodeInfo_ ;
2392    lastNodeInfo_ = new CbcNodeInfo * [maximumDepth_] ;
2393    delete [] lastNumberCuts_ ;
2394    lastNumberCuts_ = new int [maximumDepth_] ;
2395    maximumCuts_ = 100;
2396    lastNumberCuts2_ = 0;
2397    delete [] lastCut_;
2398    lastCut_ = new const OsiRowCut * [maximumCuts_];
2399    /*
2400      Used to generate bound edits for CbcPartialNodeInfo.
2401    */
2402    double * lowerBefore = new double [numberColumns] ;
2403    double * upperBefore = new double [numberColumns] ;
2404    /*
2405    Set up to run heuristics and generate cuts at the root node. The heavy
2406    lifting is hidden inside the calls to doHeuristicsAtRoot and solveWithCuts.
2407
2408    To start, tell cut generators they can be a bit more aggressive at the
2409    root node.
2410
2411    QUESTION: phase_ = 0 is documented as `initial solve', phase = 1 as `solve
2412        with cuts at root'. Is phase_ = 1 the correct indication when
2413        doHeurisiticsAtRoot is called to run heuristics outside of the main
2414        cut / heurisitc / reoptimise loop in solveWithCuts?
2415
2416      Generate cuts at the root node and reoptimise. solveWithCuts does the heavy
2417      lifting. It will iterate a generate/reoptimise loop (including reduced cost
2418      fixing) until no cuts are generated, the change in objective falls off,  or
2419      the limit on the number of rounds of cut generation is exceeded.
2420
2421      At the end of all this, any cuts will be recorded in cuts and also
2422      installed in the solver's constraint system. We'll have reoptimised, and
2423      removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
2424      adjusted accordingly).
2425
2426      Tell cut generators they can be a bit more aggressive at root node
2427
2428      TODO: Why don't we make a copy of the solution after solveWithCuts?
2429      TODO: If numberUnsatisfied == 0, don't we have a solution?
2430    */
2431    phase_ = 1;
2432    int iCutGenerator;
2433    for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
2434        // If parallel switch off global cuts
2435        if (numberThreads_) {
2436            generator_[iCutGenerator]->setGlobalCuts(false);
2437            generator_[iCutGenerator]->setGlobalCutsAtRoot(false);
2438        }
2439        CglCutGenerator * generator = generator_[iCutGenerator]->generator();
2440        generator->setAggressiveness(generator->getAggressiveness() + 100);
2441        if (!generator->canDoGlobalCuts())
2442          generator->setGlobalCuts(false);
2443    }
2444    OsiCuts cuts ;
2445    int anyAction = -1 ;
2446    numberOldActiveCuts_ = 0 ;
2447    numberNewCuts_ = 0 ;
2448    // Array to mark solution
2449    delete [] usedInSolution_;
2450    usedInSolution_ = new int[numberColumns];
2451    CoinZeroN(usedInSolution_, numberColumns);
2452    /*
2453      For printing totals and for CbcNode (numberNodes_)
2454    */
2455    numberIterations_ = 0 ;
2456    numberNodes_ = 0 ;
2457    numberNodes2_ = 0 ;
2458    maximumStatistics_ = 0;
2459    maximumDepthActual_ = 0;
2460    numberDJFixed_ = 0.0;
2461    if (!parentModel_) {
2462        if ((specialOptions_&262144) != 0) {
2463            // create empty stored cuts
2464            //storedRowCuts_ = new CglStored(solver_->getNumCols());
2465        } else if ((specialOptions_&524288) != 0 && storedRowCuts_) {
2466            // tighten and set best solution
2467            // A) tight bounds on integer variables
2468            /*
2469                storedRowCuts_ are coming in from outside, probably for nonlinear.
2470              John was unsure about origin.
2471            */
2472            const double * lower = solver_->getColLower();
2473            const double * upper = solver_->getColUpper();
2474            const double * tightLower = storedRowCuts_->tightLower();
2475            const double * tightUpper = storedRowCuts_->tightUpper();
2476            int nTightened = 0;
2477            for (int i = 0; i < numberIntegers_; i++) {
2478                int iColumn = integerVariable_[i];
2479                if (tightLower[iColumn] > lower[iColumn]) {
2480                    nTightened++;
2481                    solver_->setColLower(iColumn, tightLower[iColumn]);
2482                }
2483                if (tightUpper[iColumn] < upper[iColumn]) {
2484                    nTightened++;
2485                    solver_->setColUpper(iColumn, tightUpper[iColumn]);
2486                }
2487            }
2488            if (nTightened)
2489              COIN_DETAIL_PRINT(printf("%d tightened by alternate cuts\n", nTightened));
2490            if (storedRowCuts_->bestObjective() < bestObjective_) {
2491                // B) best solution
2492                double objValue = storedRowCuts_->bestObjective();
2493                setBestSolution(CBC_SOLUTION, objValue,
2494                                storedRowCuts_->bestSolution()) ;
2495                // Do heuristics
2496                // Allow RINS
2497                for (int i = 0; i < numberHeuristics_; i++) {
2498                    CbcHeuristicRINS * rins
2499                    = dynamic_cast<CbcHeuristicRINS *> (heuristic_[i]);
2500                    if (rins) {
2501                        rins->setLastNode(-100);
2502                    }
2503                }
2504            }
2505        }
2506    }
2507#ifdef SWITCH_VARIABLES
2508    // see if any switching variables
2509    if (numberIntegers_<solver_->getNumCols())
2510      findSwitching();
2511#endif
2512    /*
2513      Run heuristics at the root. This is the only opportunity to run FPump; it
2514      will be removed from the heuristics list by doHeuristicsAtRoot.
2515    */
2516    // See if multiple runs wanted
2517    CbcModel ** rootModels=NULL;
2518    if (!parentModel_&&multipleRootTries_%100) {
2519      double rootTimeCpu=CoinCpuTime();
2520      double startTimeRoot=CoinGetTimeOfDay();
2521      int numberRootThreads=1;
2522      /* undocumented fine tuning
2523         aabbcc where cc is number of tries
2524         bb if nonzero is number of threads
2525         aa if nonzero just do heuristics
2526      */
2527      int numberModels = multipleRootTries_%100;
2528#ifdef CBC_THREAD
2529      numberRootThreads = (multipleRootTries_/100)%100;
2530      if (!numberRootThreads)
2531        numberRootThreads=numberModels;
2532#endif
2533      int otherOptions = (multipleRootTries_/10000)%100;
2534      rootModels = new CbcModel * [numberModels];
2535      unsigned int newSeed = randomSeed_;
2536      if (newSeed==0) {
2537        double time = fabs(CoinGetTimeOfDay());
2538        while (time>=COIN_INT_MAX)
2539          time *= 0.5;
2540        newSeed = static_cast<unsigned int>(time);
2541      } else if (newSeed<0) {
2542        newSeed = 123456789;
2543#ifdef COIN_HAS_CLP
2544        OsiClpSolverInterface * clpSolver
2545          = dynamic_cast<OsiClpSolverInterface *> (solver_);
2546        if (clpSolver) {
2547          newSeed += clpSolver->getModelPtr()->randomNumberGenerator()->getSeed();
2548        }
2549#endif
2550      }
2551      CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis *> (solver_->getEmptyWarmStart());
2552      int numberIterations=0;
2553      for (int i=0;i<numberModels;i++) { 
2554        rootModels[i]=new CbcModel(*this);
2555        rootModels[i]->setNumberThreads(0);
2556        rootModels[i]->setMaximumNodes(otherOptions ? -1 : 0);
2557        rootModels[i]->setRandomSeed(newSeed+10000000*i);
2558        rootModels[i]->randomNumberGenerator()->setSeed(newSeed+50000000*i);
2559        rootModels[i]->setMultipleRootTries(0);
2560        // use seed
2561        rootModels[i]->setSpecialOptions(specialOptions_ |(4194304|8388608));
2562        rootModels[i]->setMoreSpecialOptions(moreSpecialOptions_ & 
2563                                             (~134217728));
2564        rootModels[i]->solver_->setWarmStart(basis);
2565#ifdef COIN_HAS_CLP
2566        OsiClpSolverInterface * clpSolver
2567          = dynamic_cast<OsiClpSolverInterface *> (rootModels[i]->solver_);
2568        if (clpSolver) {
2569          ClpSimplex * simplex = clpSolver->getModelPtr();
2570          if (defaultHandler_)
2571            simplex->setDefaultMessageHandler();
2572          simplex->setRandomSeed(newSeed+20000000*i);
2573          simplex->allSlackBasis();
2574          int logLevel=simplex->logLevel();
2575          if (logLevel==1)
2576            simplex->setLogLevel(0);
2577          if (i!=0) {
2578            double random = simplex->randomNumberGenerator()->randomDouble();
2579            int bias = static_cast<int>(random*(numberIterations/4));
2580            simplex->setMaximumIterations(numberIterations/2+bias);
2581            simplex->primal();
2582            simplex->setMaximumIterations(COIN_INT_MAX);
2583            simplex->dual();
2584          } else {
2585            simplex->primal();
2586            numberIterations=simplex->numberIterations();
2587          }
2588          simplex->setLogLevel(logLevel);
2589          clpSolver->setWarmStart(NULL);
2590        }
2591#endif
2592        for (int j=0;j<numberHeuristics_;j++)
2593          rootModels[i]->heuristic_[j]->setSeed(rootModels[i]->heuristic_[j]->getSeed()+100000000*i);
2594        for (int j=0;j<numberCutGenerators_;j++)
2595          rootModels[i]->generator_[j]->generator()->refreshSolver(rootModels[i]->solver_);
2596      }
2597      delete basis;
2598#ifdef CBC_THREAD
2599      if (numberRootThreads==1) {
2600#endif
2601        for (int iModel=0;iModel<numberModels;iModel++) {
2602          doRootCbcThread(rootModels[iModel]);
2603          // see if solved at root node
2604          if (rootModels[iModel]->getMaximumNodes()) {
2605            feasible=false;
2606            break;
2607          }
2608        }
2609#ifdef CBC_THREAD
2610      } else {
2611        Coin_pthread_t * threadId = new Coin_pthread_t [numberRootThreads];
2612        for (int kModel=0;kModel<numberModels;kModel+=numberRootThreads) {
2613          bool finished=false;
2614          for (int iModel=kModel;iModel<CoinMin(numberModels,kModel+numberRootThreads);iModel++) {
2615            pthread_create(&(threadId[iModel-kModel].thr), NULL, 
2616                           doRootCbcThread,
2617                           rootModels[iModel]);
2618          }
2619          // wait
2620          for (int iModel=kModel;iModel<CoinMin(numberModels,kModel+numberRootThreads);iModel++) {
2621            pthread_join(threadId[iModel-kModel].thr, NULL);
2622          }
2623          // see if solved at root node
2624          for (int iModel=kModel;iModel<CoinMin(numberModels,kModel+numberRootThreads);iModel++) {
2625            if (rootModels[iModel]->getMaximumNodes())
2626              finished=true;
2627          }
2628          if (finished) {
2629            feasible=false;
2630            break;
2631          }
2632        }
2633        delete [] threadId;
2634      }
2635#endif
2636      // sort solutions
2637      int * which = new int [numberModels];
2638      double * value = new double [numberModels];
2639      int numberSolutions=0;
2640      for (int iModel=0;iModel<numberModels;iModel++) {
2641        if (rootModels[iModel]->bestSolution()) {
2642          which[numberSolutions]=iModel;
2643          value[numberSolutions++]=
2644            -rootModels[iModel]->getMinimizationObjValue();
2645        }
2646      }
2647      char general[100];
2648      rootTimeCpu=CoinCpuTime()-rootTimeCpu;
2649      if (numberRootThreads==1)
2650        sprintf(general,"Multiple root solvers took a total of %.2f seconds\n",
2651                rootTimeCpu);
2652      else
2653        sprintf(general,"Multiple root solvers took a total of %.2f seconds (%.2f elapsed)\n",
2654                rootTimeCpu,CoinGetTimeOfDay()-startTimeRoot);
2655      messageHandler()->message(CBC_GENERAL,
2656                                messages())
2657        << general << CoinMessageEol ;
2658      CoinSort_2(value,value+numberSolutions,which);
2659      // to get name
2660      CbcHeuristicRINS dummyHeuristic;
2661      dummyHeuristic.setHeuristicName("Multiple root solvers");
2662      lastHeuristic_=&dummyHeuristic;
2663      for (int i=0;i<numberSolutions;i++) {
2664        double objValue = - value[i];
2665        if (objValue<getCutoff()) {
2666          int iModel=which[i];
2667          setBestSolution(CBC_ROUNDING,objValue,
2668                          rootModels[iModel]->bestSolution());
2669        }
2670      }
2671      lastHeuristic_=NULL;
2672      delete [] which;
2673      delete [] value;
2674    }
2675    // Do heuristics
2676    if (numberObjects_&&!rootModels)
2677        doHeuristicsAtRoot();
2678    if (solverCharacteristics_->solutionAddsCuts()) {
2679      // With some heuristics solver needs a resolve here
2680      solver_->resolve();
2681      if(!isProvenOptimal()){
2682        solver_->initialSolve();
2683      }
2684    }
2685    /*
2686      Grepping through the code, it would appear that this is a command line
2687      debugging hook.  There's no obvious place in the code where this is set to
2688      a negative value.
2689
2690      User hook, says John.
2691    */
2692    if ( intParam_[CbcMaxNumNode] < 0
2693      ||numberSolutions_>=getMaximumSolutions())
2694      eventHappened_ = true; // stop as fast as possible
2695    stoppedOnGap_ = false ;
2696    // See if can stop on gap
2697    bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
2698    if(canStopOnGap()) {
2699        if (bestPossibleObjective_ < getCutoff())
2700            stoppedOnGap_ = true ;
2701        feasible = false;
2702        //eventHappened_=true; // stop as fast as possible
2703    }
2704    /*
2705      Set up for statistics collection, if requested. Standard values are
2706      documented in CbcModel.hpp. The magic number 100 will trigger a dump of
2707      CbcSimpleIntegerDynamicPseudoCost objects (no others). Looks like another
2708      command line debugging hook.
2709    */
2710    statistics_ = NULL;
2711    // Do on switch
2712    if (doStatistics > 0 && doStatistics <= 100) {
2713        maximumStatistics_ = 10000;
2714        statistics_ = new CbcStatistics * [maximumStatistics_];
2715        memset(statistics_, 0, maximumStatistics_*sizeof(CbcStatistics *));
2716    }
2717    // See if we can add integers
2718    if (noObjects && numberIntegers_ < solver_->getNumCols() && (specialOptions_&65536) != 0 && !parentModel_)
2719        AddIntegers();
2720    /*
2721      Do an initial round of cut generation for the root node. Depending on the
2722      type of underlying solver, we may want to do this even if the initial query
2723      to the objects indicates they're satisfied.
2724
2725      solveWithCuts does the heavy lifting. It will iterate a generate/reoptimise
2726      loop (including reduced cost fixing) until no cuts are generated, the
2727      change in objective falls off,  or the limit on the number of rounds of cut
2728      generation is exceeded.
2729
2730      At the end of all this, any cuts will be recorded in cuts and also
2731      installed in the solver's constraint system. We'll have reoptimised, and
2732      removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
2733      adjusted accordingly).
2734    */
2735    int iObject ;
2736    int preferredWay ;
2737    int numberUnsatisfied = 0 ;
2738    delete [] currentSolution_;
2739    currentSolution_ = new double [numberColumns];
2740    testSolution_ = currentSolution_;
2741    memcpy(currentSolution_, solver_->getColSolution(),
2742           numberColumns*sizeof(double)) ;
2743    // point to useful information
2744    OsiBranchingInformation usefulInfo = usefulInformation();
2745
2746    for (iObject = 0 ; iObject < numberObjects_ ; iObject++) {
2747        double infeasibility =
2748            object_[iObject]->infeasibility(&usefulInfo, preferredWay) ;
2749        if (infeasibility ) numberUnsatisfied++ ;
2750    }
2751    // replace solverType
2752    double * tightBounds = NULL;
2753    if (solverCharacteristics_->tryCuts())  {
2754
2755        if (numberUnsatisfied)   {
2756            // User event
2757            if (!eventHappened_ && feasible) {
2758                if (rootModels) {
2759                  // for fixings
2760                  int numberColumns=solver_->getNumCols();
2761                  tightBounds = new double [2*numberColumns];
2762                  {
2763                    const double * lower = solver_->getColLower();
2764                    const double * upper = solver_->getColUpper();
2765                    for (int i=0;i<numberColumns;i++) {
2766                      tightBounds[2*i+0]=lower[i];
2767                      tightBounds[2*i+1]=upper[i];
2768                    }
2769                  }
2770                  int numberModels = multipleRootTries_%100;
2771                  const OsiSolverInterface ** solvers = new 
2772                    const OsiSolverInterface * [numberModels];
2773                  int numberRows=continuousSolver_->getNumRows();
2774                  int maxCuts=0;
2775                  for (int i=0;i<numberModels;i++) {
2776                    solvers[i]=rootModels[i]->solver();
2777                    const double * lower = solvers[i]->getColLower();
2778                    const double * upper = solvers[i]->getColUpper();
2779                    for (int j=0;j<numberColumns;j++) {
2780                      tightBounds[2*j+0]=CoinMax(lower[j],tightBounds[2*j+0]);
2781                      tightBounds[2*j+1]=CoinMin(upper[j],tightBounds[2*j+1]);
2782                    }
2783                    int numberRows2=solvers[i]->getNumRows();
2784                    assert (numberRows2>=numberRows);
2785                    maxCuts += numberRows2-numberRows;
2786                    // accumulate statistics
2787                    for (int j=0;j<numberCutGenerators_;j++) {
2788                      generator_[j]->addStatistics(rootModels[i]->cutGenerator(j));
2789                    }
2790                  }
2791                  for (int j=0;j<numberCutGenerators_;j++) {
2792                    generator_[j]->scaleBackStatistics(numberModels);
2793                  }
2794                  //CbcRowCuts rowCut(maxCuts);
2795                  const OsiRowCutDebugger *debugger = NULL;
2796                  if ((specialOptions_&1) != 0) 
2797                    debugger = solver_->getRowCutDebugger() ;
2798                  for (int iModel=0;iModel<numberModels;iModel++) {
2799                    int numberRows2=solvers[iModel]->getNumRows();
2800                    const CoinPackedMatrix * rowCopy = solvers[iModel]->getMatrixByRow();
2801                    const int * rowLength = rowCopy->getVectorLengths(); 
2802                    const double * elements = rowCopy->getElements();
2803                    const int * column = rowCopy->getIndices();
2804                    const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
2805                    const double * rowLower = solvers[iModel]->getRowLower();
2806                    const double * rowUpper = solvers[iModel]->getRowUpper();
2807                    for (int iRow=numberRows;iRow<numberRows2;iRow++) {
2808                      OsiRowCut rc;
2809                      rc.setLb(rowLower[iRow]);
2810                      rc.setUb(rowUpper[iRow]);
2811                      CoinBigIndex start = rowStart[iRow];
2812                      rc.setRow(rowLength[iRow],column+start,elements+start,false);
2813                      if (debugger)
2814                        CoinAssert (!debugger->invalidCut(rc));
2815                      globalCuts_.addCutIfNotDuplicate(rc);
2816                    }
2817                    //int cutsAdded=globalCuts_.numberCuts()-numberCuts;
2818                    //numberCuts += cutsAdded;
2819                    //printf("Model %d gave %d cuts (out of %d possible)\n",
2820                    //     iModel,cutsAdded,numberRows2-numberRows);
2821                  }
2822                  // normally replace global cuts
2823                  //if (!globalCuts_.())
2824                  //globalCuts_=rowCutrowCut.addCuts(globalCuts_);
2825                  //rowCut.addCuts(globalCuts_);
2826                  int nTightened=0;
2827                  bool feasible=true;
2828                  {
2829                    double tolerance=1.0e-5;
2830                    const double * lower = solver_->getColLower();
2831                    const double * upper = solver_->getColUpper();
2832                    for (int i=0;i<numberColumns;i++) {
2833                      if (tightBounds[2*i+0]>tightBounds[2*i+1]) {
2834                        feasible=false;
2835                        printf("Bad bounds on %d %g,%g was %g,%g\n",
2836                               i,tightBounds[2*i+0],tightBounds[2*i+1],
2837                               lower[i],upper[i]);
2838                      }
2839                      //int k=0;
2840                      double oldLower=lower[i];
2841                      double oldUpper=upper[i];
2842                      if (tightBounds[2*i+0]>oldLower+tolerance) {
2843                        nTightened++;
2844                        //k++;
2845                        solver_->setColLower(i,tightBounds[2*i+0]);
2846                      }
2847                      if (tightBounds[2*i+1]<oldUpper-tolerance) {
2848                        nTightened++;
2849                        //k++;
2850                        solver_->setColUpper(i,tightBounds[2*i+1]);
2851                      }
2852                      //if (k)
2853                      //printf("new bounds on %d %g,%g was %g,%g\n",
2854                      //       i,tightBounds[2*i+0],tightBounds[2*i+1],
2855                      //       oldLower,oldUpper);
2856                    }
2857                    if (!feasible)
2858                      abort(); // deal with later
2859                  }
2860                  delete [] tightBounds;
2861                  tightBounds=NULL;
2862                  char printBuffer[200];
2863                  sprintf(printBuffer,"%d solvers added %d different cuts out of pool of %d",
2864                          numberModels,globalCuts_.sizeRowCuts(),maxCuts);
2865                  messageHandler()->message(CBC_GENERAL, messages())
2866                    << printBuffer << CoinMessageEol ;
2867                  if (nTightened) {
2868                    sprintf(printBuffer,"%d bounds were tightened",
2869                          nTightened);
2870                    messageHandler()->message(CBC_GENERAL, messages())
2871                      << printBuffer << CoinMessageEol ;
2872                  }
2873                  delete [] solvers;
2874                }
2875                if (!parentModel_&&(moreSpecialOptions_&67108864) != 0) {
2876                  // load cuts from file
2877                  FILE * fp = fopen("global.cuts","rb");
2878                  if (fp) {
2879                    size_t nRead;
2880                    int numberColumns=solver_->getNumCols();
2881                    int numCols;
2882                    nRead = fread(&numCols, sizeof(int), 1, fp);
2883                    if (nRead != 1)
2884                      throw("Error in fread");
2885                    if (numberColumns!=numCols) {
2886                      printf("Mismatch on columns %d %d\n",numberColumns,numCols);
2887                      fclose(fp);
2888                    } else {
2889                      // If rootModel just do some
2890                      double threshold=-1.0;
2891                      if (!multipleRootTries_)
2892                        threshold=0.5;
2893                      int initialCuts=0;
2894                      int initialGlobal = globalCuts_.sizeRowCuts();
2895                      double * elements = new double [numberColumns+2];
2896                      int * indices = new int [numberColumns];
2897                      int numberEntries=1;
2898                      while (numberEntries>0) {
2899                        nRead = fread(&numberEntries, sizeof(int), 1, fp);
2900                        if (nRead != 1)
2901                          throw("Error in fread");
2902                        double randomNumber=randomNumberGenerator_.randomDouble();
2903                        if (numberEntries>0) {
2904                          initialCuts++;
2905                          nRead = fread(elements, sizeof(double), numberEntries+2, fp);
2906                          if (nRead != static_cast<size_t>(numberEntries+2))
2907                            throw("Error in fread");
2908                          nRead = fread(indices, sizeof(int), numberEntries, fp);
2909                          if (nRead != static_cast<size_t>(numberEntries))
2910                            throw("Error in fread");
2911                          if (randomNumber>threshold) {
2912                            OsiRowCut rc;
2913                            rc.setLb(elements[numberEntries]);
2914                            rc.setUb(elements[numberEntries+1]);
2915                            rc.setRow(numberEntries,indices,elements,
2916                                      false);
2917                            rc.setGloballyValidAsInteger(2);
2918                            globalCuts_.addCutIfNotDuplicate(rc) ;
2919                          } 
2920                        }
2921                      }
2922                      fclose(fp);
2923                      // fixes
2924                      int nTightened=0;
2925                      fp = fopen("global.fix","rb");
2926                      if (fp) {
2927                        nRead = fread(indices, sizeof(int), 2, fp);
2928                        if (nRead != 2)
2929                          throw("Error in fread");
2930                        if (numberColumns!=indices[0]) {
2931                          printf("Mismatch on columns %d %d\n",numberColumns,
2932                                 indices[0]);
2933                        } else {
2934                          indices[0]=1;
2935                          while (indices[0]>=0) {
2936                            nRead = fread(indices, sizeof(int), 2, fp);
2937                            if (nRead != 2)
2938                              throw("Error in fread");
2939                            int iColumn=indices[0];
2940                            if (iColumn>=0) {
2941                              nTightened++;
2942                              nRead = fread(elements, sizeof(double), 4, fp);
2943                              if (nRead != 4)
2944                                throw("Error in fread");
2945                              solver_->setColLower(iColumn,elements[0]);
2946                              solver_->setColUpper(iColumn,elements[1]);
2947                            } 
2948                          }
2949                        }
2950                      }
2951                      if (fp)
2952                        fclose(fp);
2953                      char printBuffer[200];
2954                      sprintf(printBuffer,"%d cuts read in of which %d were unique, %d bounds tightened",
2955                             initialCuts,
2956                             globalCuts_.sizeRowCuts()-initialGlobal,nTightened); 
2957                      messageHandler()->message(CBC_GENERAL, messages())
2958                        << printBuffer << CoinMessageEol ;
2959                      delete [] elements;
2960                      delete [] indices;
2961                    }
2962                  }
2963                }
2964                feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
2965                                         NULL);
2966                if (multipleRootTries_&&
2967                    (moreSpecialOptions_&134217728)!=0) {
2968                  FILE * fp=NULL;
2969                  size_t nRead;
2970                  int numberColumns=solver_->getNumCols();
2971                  int initialCuts=0;
2972                  if ((moreSpecialOptions_&134217728)!=0) {
2973                    // append so go down to end
2974                    fp = fopen("global.cuts","r+b");
2975                    if (fp) {
2976                      int numCols;
2977                      nRead = fread(&numCols, sizeof(int), 1, fp);
2978                      if (nRead != 1)
2979                        throw("Error in fread");
2980                      if (numberColumns!=numCols) {
2981                        printf("Mismatch on columns %d %d\n",numberColumns,numCols);
2982                        fclose(fp);
2983                        fp=NULL;
2984                      }
2985                    }
2986                  }
2987                  double * elements = new double [numberColumns+2];
2988                  int * indices = new int [numberColumns];
2989                  if (fp) {
2990                    int numberEntries=1;
2991                    while (numberEntries>0) {
2992                      fpos_t position;
2993                      fgetpos(fp, &position);
2994                      nRead = fread(&numberEntries, sizeof(int), 1, fp);
2995                      if (nRead != 1)
2996                        throw("Error in fread");
2997                      if (numberEntries>0) {
2998                        initialCuts++;
2999                        nRead = fread(elements, sizeof(double), numberEntries+2, fp);
3000                        if (nRead != static_cast<size_t>(numberEntries+2))
3001                          throw("Error in fread");
3002                        nRead = fread(indices, sizeof(int), numberEntries, fp);
3003                        if (nRead != static_cast<size_t>(numberEntries))
3004                          throw("Error in fread");
3005                      } else {
3006                        // end
3007                        fsetpos(fp, &position);
3008                      }
3009                    }
3010                  } else {
3011                    fp = fopen("global.cuts","wb");
3012                    size_t nWrite;
3013                    nWrite=fwrite(&numberColumns,sizeof(int),1,fp);
3014                    if (nWrite != 1)
3015                      throw("Error in fwrite");
3016                  }
3017                  size_t nWrite;
3018                  // now append binding cuts
3019                  int numberC=continuousSolver_->getNumRows();
3020                  int numberRows=solver_->getNumRows();
3021                  printf("Saving %d cuts (up from %d)\n",
3022                         initialCuts+numberRows-numberC,initialCuts);
3023                  const double * rowLower = solver_->getRowLower();
3024                  const double * rowUpper = solver_->getRowUpper();
3025                  // Row copy
3026                  CoinPackedMatrix matrixByRow(*solver_->getMatrixByRow());
3027                  const double * elementByRow = matrixByRow.getElements();
3028                  const int * column = matrixByRow.getIndices();
3029                  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
3030                  const int * rowLength = matrixByRow.getVectorLengths();
3031                  for (int iRow=numberC;iRow<numberRows;iRow++) {
3032                    int n=rowLength[iRow];
3033                    assert (n);
3034                    CoinBigIndex start=rowStart[iRow];
3035                    memcpy(elements,elementByRow+start,n*sizeof(double));
3036                    memcpy(indices,column+start,n*sizeof(int));
3037                    elements[n]=rowLower[iRow];
3038                    elements[n+1]=rowUpper[iRow];
3039                    nWrite=fwrite(&n,sizeof(int),1,fp);
3040                    if (nWrite != 1)
3041                      throw("Error in fwrite");
3042                    nWrite=fwrite(elements,sizeof(double),n+2,fp);
3043                    if (nWrite != static_cast<size_t>(n+2))
3044                      throw("Error in fwrite");
3045                    nWrite=fwrite(indices,sizeof(int),n,fp);
3046                    if (nWrite != static_cast<size_t>(n))
3047                      throw("Error in fwrite");
3048                  }
3049                  // eof marker
3050                  int eofMarker=-1;
3051                  nWrite=fwrite(&eofMarker,sizeof(int),1,fp);
3052                  if (nWrite != 1)
3053                    throw("Error in fwrite");
3054                  fclose(fp);
3055                  // do tighter bounds (? later extra to original columns)
3056                  int nTightened=0;
3057                  const double * lower = solver_->getColLower();
3058                  const double * upper = solver_->getColUpper();
3059                  const double * originalLower = continuousSolver_->getColLower();
3060                  const double * originalUpper = continuousSolver_->getColUpper();
3061                  double tolerance=1.0e-5;
3062                  for (int i=0;i<numberColumns;i++) {
3063                    if (lower[i]>originalLower[i]+tolerance) {
3064                      nTightened++;
3065                    }
3066                    if (upper[i]<originalUpper[i]-tolerance) {
3067                      nTightened++;
3068                    }
3069                  }
3070                  if (nTightened) {
3071                    fp = fopen("global.fix","wb");
3072                    size_t nWrite;
3073                    indices[0]=numberColumns;
3074                    if (originalColumns_)
3075                      indices[1]=COIN_INT_MAX;
3076                    else
3077                      indices[1]=-1;
3078                    nWrite=fwrite(indices,sizeof(int),2,fp);
3079                    if (nWrite != 2)
3080                      throw("Error in fwrite");
3081                    for (int i=0;i<numberColumns;i++) {
3082                      int nTightened=0;
3083                      if (lower[i]>originalLower[i]+tolerance) {
3084                        nTightened++;
3085                      }
3086                      if (upper[i]<originalUpper[i]-tolerance) {
3087                        nTightened++;
3088                      }
3089                      if (nTightened) {
3090                        indices[0]=i;
3091                        if (originalColumns_)
3092                          indices[1]=originalColumns_[i];
3093                        elements[0]=lower[i];
3094                        elements[1]=upper[i];
3095                        elements[2]=originalLower[i];
3096                        elements[3]=originalUpper[i];
3097                        nWrite=fwrite(indices,sizeof(int),2,fp);
3098                        if (nWrite != 2)
3099                          throw("Error in fwrite");
3100                        nWrite=fwrite(elements,sizeof(double),4,fp);
3101                        if (nWrite != 4)
3102                          throw("Error in fwrite");
3103                      }
3104                    }
3105                    // eof marker
3106                    indices[0]=-1;
3107                    nWrite=fwrite(indices,sizeof(int),2,fp);
3108                    if (nWrite != 2)
3109                      throw("Error in fwrite");
3110                    fclose(fp);
3111                  }
3112                  delete [] elements;
3113                  delete [] indices;
3114                }
3115                if ((specialOptions_&524288) != 0 && !parentModel_
3116                        && storedRowCuts_) {
3117                    if (feasible) {
3118                        /* pick up stuff and try again
3119                        add cuts, maybe keep around
3120                        do best solution and if so new heuristics
3121                        obviously tighten bounds
3122                        */
3123                        // A and B probably done on entry
3124                        // A) tight bounds on integer variables
3125                        const double * lower = solver_->getColLower();
3126                        const double * upper = solver_->getColUpper();
3127                        const double * tightLower = storedRowCuts_->tightLower();
3128                        const double * tightUpper = storedRowCuts_->tightUpper();
3129                        int nTightened = 0;
3130                        for (int i = 0; i < numberIntegers_; i++) {
3131                            int iColumn = integerVariable_[i];
3132                            if (tightLower[iColumn] > lower[iColumn]) {
3133                                nTightened++;
3134                                solver_->setColLower(iColumn, tightLower[iColumn]);
3135                            }
3136                            if (tightUpper[iColumn] < upper[iColumn]) {
3137                                nTightened++;
3138                                solver_->setColUpper(iColumn, tightUpper[iColumn]);
3139                            }
3140                        }
3141                        if (nTightened)
3142                          COIN_DETAIL_PRINT(printf("%d tightened by alternate cuts\n", nTightened));
3143                        if (storedRowCuts_->bestObjective() < bestObjective_) {
3144                            // B) best solution
3145                            double objValue = storedRowCuts_->bestObjective();
3146                            setBestSolution(CBC_SOLUTION, objValue,
3147                                            storedRowCuts_->bestSolution()) ;
3148                            // Do heuristics
3149                            // Allow RINS
3150                            for (int i = 0; i < numberHeuristics_; i++) {
3151                                CbcHeuristicRINS * rins
3152                                = dynamic_cast<CbcHeuristicRINS *> (heuristic_[i]);
3153                                if (rins) {
3154                                    rins->setLastNode(-100);
3155                                }
3156                            }
3157                            doHeuristicsAtRoot();
3158                        }
3159#ifdef JJF_ZERO
3160                        int nCuts = storedRowCuts_->sizeRowCuts();
3161                        // add to global list
3162                        for (int i = 0; i < nCuts; i++) {
3163                            OsiRowCut newCut(*storedRowCuts_->rowCutPointer(i));
3164                            newCut.setGloballyValidAsInteger(2);
3165                            newCut.mutableRow().setTestForDuplicateIndex(false);
3166                            globalCuts_.insert(newCut) ;
3167                        }
3168#else
3169                        addCutGenerator(storedRowCuts_, -99, "Stored from previous run",
3170                                        true, false, false, -200);
3171#endif
3172                        // Set cuts as active
3173                        delete [] addedCuts_ ;
3174                        maximumNumberCuts_ = cuts.sizeRowCuts();
3175                        if (maximumNumberCuts_) {
3176                            addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
3177                        } else {
3178                            addedCuts_ = NULL;
3179                        }
3180                        for (int i = 0; i < maximumNumberCuts_; i++)
3181                            addedCuts_[i] = new CbcCountRowCut(*cuts.rowCutPtr(i),
3182                                                               NULL, -1, -1, 2);
3183                        COIN_DETAIL_PRINT(printf("size %d\n", cuts.sizeRowCuts()));
3184                        cuts = OsiCuts();
3185                        currentNumberCuts_ = maximumNumberCuts_;
3186                        feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
3187                                                 NULL);
3188                        for (int i = 0; i < maximumNumberCuts_; i++)
3189                            delete addedCuts_[i];
3190                    }
3191                    delete storedRowCuts_;
3192                    storedRowCuts_ = NULL;
3193                }
3194            } else {
3195                feasible = false;
3196            }
3197        }       else if (solverCharacteristics_->solutionAddsCuts() ||
3198                   solverCharacteristics_->alwaysTryCutsAtRootNode()) {
3199            // may generate cuts and turn the solution
3200            //to an infeasible one
3201            feasible = solveWithCuts(cuts, 2,
3202                                     NULL);
3203        }
3204    }
3205    if (rootModels) {
3206      int numberModels = multipleRootTries_%100;
3207      for (int i=0;i<numberModels;i++)
3208        delete rootModels[i];
3209      delete [] rootModels;
3210    }
3211    // check extra info on feasibility
3212    if (!solverCharacteristics_->mipFeasible())
3213        feasible = false;
3214    // If max nodes==0 - don't do strong branching
3215    if (!getMaximumNodes()) {
3216      if (feasible)
3217        feasible=false;
3218      else
3219        setMaximumNodes(1); //allow to stop on success
3220    }
3221    topOfTree_=NULL;
3222#ifdef CLP_RESOLVE
3223    if ((moreSpecialOptions_&2097152)!=0&&!parentModel_&&feasible) {
3224      OsiClpSolverInterface * clpSolver
3225        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3226      if (clpSolver)
3227        resolveClp(clpSolver,0);
3228    }
3229#endif
3230    // make cut generators less aggressive
3231    for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
3232        CglCutGenerator * generator = generator_[iCutGenerator]->generator();
3233        generator->setAggressiveness(generator->getAggressiveness() - 100);
3234    }
3235    currentNumberCuts_ = numberNewCuts_ ;
3236    if (solverCharacteristics_->solutionAddsCuts()) {
3237      // With some heuristics solver needs a resolve here (don't know if this is bug in heuristics)
3238      solver_->resolve();
3239      if(!isProvenOptimal()){
3240        solver_->initialSolve();
3241      }
3242    }
3243    // See if can stop on gap
3244    bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
3245    if(canStopOnGap()) {
3246        if (bestPossibleObjective_ < getCutoff())
3247            stoppedOnGap_ = true ;
3248        feasible = false;
3249    }
3250    // User event
3251    if (eventHappened_)
3252        feasible = false;
3253#if defined(COIN_HAS_CLP)&&defined(COIN_HAS_CPX)
3254    /*
3255      This is the notion of using Cbc stuff to get going, then calling cplex to
3256      finish off.
3257    */
3258    if (feasible && (specialOptions_&16384) != 0 && fastNodeDepth_ == -2 && !parentModel_) {
3259        // Use Cplex to do search!
3260        double time1 = CoinCpuTime();
3261        OsiClpSolverInterface * clpSolver
3262        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3263        OsiCpxSolverInterface cpxSolver;
3264        double direction = clpSolver->getObjSense();
3265        cpxSolver.setObjSense(direction);
3266        // load up cplex
3267        const CoinPackedMatrix * matrix = continuousSolver_->getMatrixByCol();
3268        const double * rowLower = continuousSolver_->getRowLower();
3269        const double * rowUpper = continuousSolver_->getRowUpper();
3270        const double * columnLower = continuousSolver_->getColLower();
3271        const double * columnUpper = continuousSolver_->getColUpper();
3272        const double * objective = continuousSolver_->getObjCoefficients();
3273        cpxSolver.loadProblem(*matrix, columnLower, columnUpper,
3274                              objective, rowLower, rowUpper);
3275        double * setSol = new double [numberIntegers_];
3276        int * setVar = new int [numberIntegers_];
3277        // cplex doesn't know about objective offset
3278        double offset = clpSolver->getModelPtr()->objectiveOffset();
3279        for (int i = 0; i < numberIntegers_; i++) {
3280            int iColumn = integerVariable_[i];
3281            cpxSolver.setInteger(iColumn);
3282            if (bestSolution_) {
3283                setSol[i] = bestSolution_[iColumn];
3284                setVar[i] = iColumn;
3285            }
3286        }
3287        CPXENVptr env = cpxSolver.getEnvironmentPtr();
3288        CPXLPptr lpPtr = cpxSolver.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL);
3289        cpxSolver.switchToMIP();
3290        if (bestSolution_) {
3291            CPXcopymipstart(env, lpPtr, numberIntegers_, setVar, setSol);
3292        }
3293        if (clpSolver->getNumRows() > continuousSolver_->getNumRows() && false) {
3294            // add cuts
3295            const CoinPackedMatrix * matrix = clpSolver->getMatrixByRow();
3296            const double * rhs = clpSolver->getRightHandSide();
3297            const char * rowSense = clpSolver->getRowSense();
3298            const double * elementByRow = matrix->getElements();
3299            const int * column = matrix->getIndices();
3300            const CoinBigIndex * rowStart = matrix->getVectorStarts();
3301            const int * rowLength = matrix->getVectorLengths();
3302            int nStart = continuousSolver_->getNumRows();
3303            int nRows = clpSolver->getNumRows();
3304            int size = rowStart[nRows-1] + rowLength[nRows-1] -
3305                       rowStart[nStart];
3306            int nAdd = 0;
3307            double * rmatval = new double [size];
3308            int * rmatind = new int [size];
3309            int * rmatbeg = new int [nRows-nStart+1];
3310            size = 0;
3311            rmatbeg[0] = 0;
3312            for (int i = nStart; i < nRows; i++) {
3313                for (int k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
3314                    rmatind[size] = column[k];
3315                    rmatval[size++] = elementByRow[k];
3316                }
3317                nAdd++;
3318                rmatbeg[nAdd] = size;
3319            }
3320            CPXaddlazyconstraints(env, lpPtr, nAdd, size,
3321                                  rhs, rowSense, rmatbeg,
3322                                  rmatind, rmatval, NULL);
3323            CPXsetintparam( env, CPX_PARAM_REDUCE,
3324                            // CPX_PREREDUCE_NOPRIMALORDUAL (0)
3325                            CPX_PREREDUCE_PRIMALONLY);
3326        }
3327        if (getCutoff() < 1.0e50) {
3328            double useCutoff = getCutoff() + offset;
3329            if (bestObjective_ < 1.0e50)
3330                useCutoff = bestObjective_ + offset + 1.0e-7;
3331            cpxSolver.setDblParam(OsiDualObjectiveLimit, useCutoff*
3332                                  direction);
3333            if ( direction > 0.0 )
3334                CPXsetdblparam( env, CPX_PARAM_CUTUP, useCutoff ) ; // min
3335            else
3336                CPXsetdblparam( env, CPX_PARAM_CUTLO, useCutoff ) ; // max
3337        }
3338        CPXsetdblparam(env, CPX_PARAM_EPGAP, dblParam_[CbcAllowableFractionGap]);
3339        delete [] setSol;
3340        delete [] setVar;
3341        char printBuffer[200];
3342        if (offset) {
3343            sprintf(printBuffer, "Add %g to all Cplex messages for true objective",
3344                    -offset);
3345            messageHandler()->message(CBC_GENERAL, messages())
3346            << printBuffer << CoinMessageEol ;
3347            cpxSolver.setDblParam(OsiObjOffset, offset);
3348        }
3349        cpxSolver.branchAndBound();
3350        double timeTaken = CoinCpuTime() - time1;
3351        sprintf(printBuffer, "Cplex took %g seconds",
3352                timeTaken);
3353        messageHandler()->message(CBC_GENERAL, messages())
3354        << printBuffer << CoinMessageEol ;
3355        numberExtraNodes_ = CPXgetnodecnt(env, lpPtr);
3356        numberExtraIterations_ = CPXgetmipitcnt(env, lpPtr);
3357        double value = cpxSolver.getObjValue() * direction;
3358        if (cpxSolver.isProvenOptimal() && value <= getCutoff()) {
3359            feasible = true;
3360            clpSolver->setWarmStart(NULL);
3361            // try and do solution
3362            double * newSolution =
3363                CoinCopyOfArray(cpxSolver.getColSolution(),
3364                                getNumCols());
3365            setBestSolution(CBC_STRONGSOL, value, newSolution) ;
3366            delete [] newSolution;
3367        }
3368        feasible = false;
3369    }
3370#endif
3371    if (!parentModel_&&(moreSpecialOptions_&268435456) != 0) {
3372      // try idiotic idea
3373      CbcObject * obj = new CbcIdiotBranch(this);
3374      obj->setPriority(1); // temp
3375      addObjects(1, &obj);
3376      delete obj;
3377    }
3378   
3379    /*
3380      A hook to use clp to quickly explore some part of the tree.
3381    */
3382    if (fastNodeDepth_ == 1000 &&/*!parentModel_*/(specialOptions_&2048) == 0) {
3383        fastNodeDepth_ = -1;
3384        CbcObject * obj =
3385            new CbcFollowOn(this);
3386        obj->setPriority(1);
3387        addObjects(1, &obj);
3388        delete obj;
3389    }
3390    int saveNumberSolves = numberSolves_;
3391    int saveNumberIterations = numberIterations_;
3392    if ((fastNodeDepth_ >= 0||(moreSpecialOptions_&33554432)!=0)
3393        &&/*!parentModel_*/(specialOptions_&2048) == 0) {
3394        // add in a general depth object doClp
3395        int type = (fastNodeDepth_ <= 100) ? fastNodeDepth_ : -(fastNodeDepth_ - 100);
3396        if ((moreSpecialOptions_&33554432)!=0)
3397          type=12;
3398        else
3399          fastNodeDepth_ += 1000000;     // mark as done
3400        CbcObject * obj =
3401            new CbcGeneralDepth(this, type);
3402        addObjects(1, &obj);
3403        delete obj;
3404        // fake number of objects
3405        numberObjects_--;
3406        if (parallelMode() < -1) {
3407            // But make sure position is correct
3408            OsiObject * obj2 = object_[numberObjects_];
3409            obj = dynamic_cast<CbcObject *> (obj2);
3410            assert (obj);
3411            obj->setPosition(numberObjects_);
3412        }
3413    }
3414#ifdef COIN_HAS_CLP
3415#ifdef NO_CRUNCH
3416    if (true) {
3417        OsiClpSolverInterface * clpSolver
3418        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3419        if (clpSolver && !parentModel_) {
3420            ClpSimplex * clpSimplex = clpSolver->getModelPtr();
3421            clpSimplex->setSpecialOptions(clpSimplex->specialOptions() | 131072);
3422            //clpSimplex->startPermanentArrays();
3423            clpSimplex->setPersistenceFlag(2);
3424        }
3425    }
3426#endif
3427#endif
3428// Save copy of solver
3429    OsiSolverInterface * saveSolver = NULL;
3430    if (!parentModel_ && (specialOptions_&(512 + 32768)) != 0)
3431        saveSolver = solver_->clone();
3432    double checkCutoffForRestart = 1.0e100;
3433    saveModel(saveSolver, &checkCutoffForRestart, &feasible);
3434    if ((specialOptions_&262144) != 0 && !parentModel_) {
3435        // Save stuff and return!
3436        storedRowCuts_->saveStuff(bestObjective_, bestSolution_,
3437                                  solver_->getColLower(),
3438                                  solver_->getColUpper());
3439        delete [] lowerBefore;
3440        delete [] upperBefore;
3441        delete saveSolver;
3442        if (flipObjective)
3443          flipModel();
3444        return;
3445    }
3446    /*
3447      We've taken the continuous relaxation as far as we can. Time to branch.
3448      The first order of business is to actually create a node. chooseBranch
3449      currently uses strong branching to evaluate branch object candidates,
3450      unless forced back to simple branching. If chooseBranch concludes that a
3451      branching candidate is monotone (anyAction == -1) or infeasible (anyAction
3452      == -2) when forced to integer values, it returns here immediately.
3453
3454      Monotone variables trigger a call to resolve(). If the problem remains
3455      feasible, try again to choose a branching variable. At the end of the loop,
3456      resolved == true indicates that some variables were fixed.
3457
3458      Loss of feasibility will result in the deletion of newNode.
3459    */
3460
3461    bool resolved = false ;
3462    CbcNode *newNode = NULL ;
3463    numberFixedAtRoot_ = 0;
3464    numberFixedNow_ = 0;
3465    if (!parentModel_&&(moreSpecialOptions2_&2)!=0) {
3466#ifdef COIN_HAS_CLP
3467      OsiClpSolverInterface * clpSolver
3468        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3469      if (clpSolver) {
3470        if (getCutoff()>1.0e20) {
3471          printf("Zapping costs\n");
3472          int numberColumns=solver_->getNumCols();
3473          double * zeroCost = new double [numberColumns];
3474          // could make small random
3475          memset(zeroCost,0,numberColumns*sizeof(double));
3476          solver_->setObjective(zeroCost);
3477          double objValue = solver_->getObjValue();
3478          solver_->setDblParam(OsiObjOffset,-objValue);
3479          clpSolver->getModelPtr()->setObjectiveValue(objValue);
3480          delete [] zeroCost;
3481        } else {
3482          moreSpecialOptions2_ &= ~2;
3483        }
3484      } else {
3485#endif
3486          moreSpecialOptions2_ &= ~2;
3487#ifdef COIN_HAS_CLP
3488      }
3489#endif
3490    }
3491    int numberIterationsAtContinuous = numberIterations_;
3492    //solverCharacteristics_->setSolver(solver_);
3493    if (feasible) {
3494#define HOTSTART -1
3495#if HOTSTART<0
3496        if (bestSolution_ && !parentModel_ && !hotstartSolution_ &&
3497                (moreSpecialOptions_&1024) != 0) {
3498            // Set priorities so only branch on ones we need to
3499            // use djs and see if only few branches needed
3500#ifndef NDEBUG
3501            double integerTolerance = getIntegerTolerance() ;
3502#endif
3503            bool possible = true;
3504            const double * saveLower = continuousSolver_->getColLower();
3505            const double * saveUpper = continuousSolver_->getColUpper();
3506            for (int i = 0; i < numberObjects_; i++) {
3507                const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object_[i]);
3508                if (thisOne) {
3509                    int iColumn = thisOne->columnNumber();
3510                    if (saveUpper[iColumn] > saveLower[iColumn] + 1.5) {
3511                        possible = false;
3512                        break;
3513                    }
3514                } else {
3515                    possible = false;
3516                    break;
3517                }
3518            }
3519            if (possible) {
3520                OsiSolverInterface * solver = continuousSolver_->clone();
3521                int numberColumns = solver->getNumCols();
3522                for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
3523                    double value = bestSolution_[iColumn] ;
3524                    value = CoinMax(value, saveLower[iColumn]) ;
3525                    value = CoinMin(value, saveUpper[iColumn]) ;
3526                    value = floor(value + 0.5);
3527                    if (solver->isInteger(iColumn)) {
3528                        solver->setColLower(iColumn, value);
3529                        solver->setColUpper(iColumn, value);
3530                    }
3531                }
3532                solver->setHintParam(OsiDoDualInResolve, false, OsiHintTry);
3533                // objlim and all slack
3534                double direction = solver->getObjSense();
3535                solver->setDblParam(OsiDualObjectiveLimit, 1.0e50*direction);
3536                CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis *> (solver->getEmptyWarmStart());
3537                solver->setWarmStart(basis);
3538                delete basis;
3539                bool changed = true;
3540                hotstartPriorities_ = new int [numberColumns];
3541                for (int iColumn = 0; iColumn < numberColumns; iColumn++)
3542                    hotstartPriorities_[iColumn] = 1;
3543                while (changed) {
3544                    changed = false;
3545                    solver->resolve();
3546                    if (!solver->isProvenOptimal()) {
3547                        possible = false;
3548                        break;
3549                    }
3550                    const double * dj = solver->getReducedCost();
3551                    const double * colLower = solver->getColLower();
3552                    const double * colUpper = solver->getColUpper();
3553                    const double * solution = solver->getColSolution();
3554                    int nAtLbNatural = 0;
3555                    int nAtUbNatural = 0;
3556                    int nZeroDj = 0;
3557                    int nForced = 0;
3558                    for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
3559                        double value = solution[iColumn] ;
3560                        value = CoinMax(value, saveLower[iColumn]) ;
3561                        value = CoinMin(value, saveUpper[iColumn]) ;
3562                        if (solver->isInteger(iColumn)) {
3563                            assert(fabs(value - solution[iColumn]) <= integerTolerance) ;
3564                            if (hotstartPriorities_[iColumn] == 1) {
3565                                if (dj[iColumn] < -1.0e-6) {
3566                                    // negative dj
3567                                    if (saveUpper[iColumn] == colUpper[iColumn]) {
3568                                        nAtUbNatural++;
3569                                        hotstartPriorities_[iColumn] = 2;
3570                                        solver->setColLower(iColumn, saveLower[iColumn]);
3571                                        solver->setColUpper(iColumn, saveUpper[iColumn]);
3572                                    } else {
3573                                        nForced++;
3574                                    }
3575                                } else if (dj[iColumn] > 1.0e-6) {
3576                                    // positive dj
3577                                    if (saveLower[iColumn] == colLower[iColumn]) {
3578                                        nAtLbNatural++;
3579                                        hotstartPriorities_[iColumn] = 2;
3580                                        solver->setColLower(iColumn, saveLower[iColumn]);
3581                                        solver->setColUpper(iColumn, saveUpper[iColumn]);
3582                                    } else {
3583                                        nForced++;
3584                                    }
3585                                } else {
3586                                    // zero dj
3587                                    nZeroDj++;
3588                                }
3589                            }
3590                        }
3591                    }
3592#ifdef CLP_INVESTIGATE
3593                    printf("%d forced, %d naturally at lower, %d at upper - %d zero dj\n",
3594                           nForced, nAtLbNatural, nAtUbNatural, nZeroDj);
3595#endif
3596                    if (nAtLbNatural || nAtUbNatural) {
3597                        changed = true;
3598                    } else {
3599                        if (nForced + nZeroDj > 5000 ||
3600                                (nForced + nZeroDj)*2 > numberIntegers_)
3601                            possible = false;
3602                    }
3603                }
3604                delete solver;
3605            }
3606            if (possible) {
3607                setHotstartSolution(bestSolution_);
3608                if (!saveCompare) {
3609                    // create depth first comparison
3610                    saveCompare = nodeCompare_;
3611                    // depth first
3612                    nodeCompare_ = new CbcCompareDepth();
3613                    tree_->setComparison(*nodeCompare_) ;
3614                }
3615            } else {
3616                delete [] hotstartPriorities_;
3617                hotstartPriorities_ = NULL;
3618            }
3619        }
3620#endif
3621#if HOTSTART>0
3622        if (hotstartSolution_ && !hotstartPriorities_) {
3623            // Set up hot start
3624            OsiSolverInterface * solver = solver_->clone();
3625            double direction = solver_->getObjSense() ;
3626            int numberColumns = solver->getNumCols();
3627            double * saveLower = CoinCopyOfArray(solver->getColLower(), numberColumns);
3628            double * saveUpper = CoinCopyOfArray(solver->getColUpper(), numberColumns);
3629            // move solution
3630            solver->setColSolution(hotstartSolution_);
3631            // point to useful information
3632            const double * saveSolution = testSolution_;
3633            testSolution_ = solver->getColSolution();
3634            OsiBranchingInformation usefulInfo = usefulInformation();
3635            testSolution_ = saveSolution;
3636            /*
3637            Run through the objects and use feasibleRegion() to set variable bounds
3638            so as to fix the variables specified in the objects at their value in this
3639            solution. Since the object list contains (at least) one object for every
3640            integer variable, this has the effect of fixing all integer variables.
3641            */
3642            for (int i = 0; i < numberObjects_; i++)
3643                object_[i]->feasibleRegion(solver, &usefulInfo);
3644            solver->resolve();
3645            assert (solver->isProvenOptimal());
3646            double gap = CoinMax((solver->getObjValue() - solver_->getObjValue()) * direction, 0.0) ;
3647            const double * dj = solver->getReducedCost();
3648            const double * colLower = solver->getColLower();
3649            const double * colUpper = solver->getColUpper();
3650            const double * solution = solver->getColSolution();
3651            int nAtLbNatural = 0;
3652            int nAtUbNatural = 0;
3653            int nAtLbNaturalZero = 0;
3654            int nAtUbNaturalZero = 0;
3655            int nAtLbFixed = 0;
3656            int nAtUbFixed = 0;
3657            int nAtOther = 0;
3658            int nAtOtherNatural = 0;
3659            int nNotNeeded = 0;
3660            delete [] hotstartSolution_;
3661            hotstartSolution_ = new double [numberColumns];
3662            delete [] hotstartPriorities_;
3663            hotstartPriorities_ = new int [numberColumns];
3664            int * order = (int *) saveUpper;
3665            int nFix = 0;
3666            double bestRatio = COIN_DBL_MAX;
3667            for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
3668                double value = solution[iColumn] ;
3669                value = CoinMax(value, saveLower[iColumn]) ;
3670                value = CoinMin(value, saveUpper[iColumn]) ;
3671                double sortValue = COIN_DBL_MAX;
3672                if (solver->isInteger(iColumn)) {
3673                    assert(fabs(value - solution[iColumn]) <= 1.0e-5) ;
3674                    double value2 = floor(value + 0.5);
3675                    if (dj[iColumn] < -1.0e-6) {
3676                        // negative dj
3677                        //assert (value2==colUpper[iColumn]);
3678                        if (saveUpper[iColumn] == colUpper[iColumn]) {
3679                            nAtUbNatural++;
3680                            sortValue = 0.0;
3681                            double value = -dj[iColumn];
3682                            if (value > gap)
3683                                nFix++;
3684                            else if (gap < value*bestRatio)
3685                                bestRatio = gap / value;
3686                            if (saveLower[iColumn] != colLower[iColumn]) {
3687                                nNotNeeded++;
3688                                sortValue = 1.0e20;
3689                            }
3690                        } else if (saveLower[iColumn] == colUpper[iColumn]) {
3691                            nAtLbFixed++;
3692                            sortValue = dj[iColumn];
3693                        } else {
3694                            nAtOther++;
3695                            sortValue = 0.0;
3696                            if (saveLower[iColumn] != colLower[iColumn] &&
3697                                    saveUpper[iColumn] != colUpper[iColumn]) {
3698                                nNotNeeded++;
3699                                sortValue = 1.0e20;
3700                            }
3701                        }
3702                    } else if (dj[iColumn] > 1.0e-6) {
3703                        // positive dj
3704                        //assert (value2==colLower[iColumn]);
3705                        if (saveLower[iColumn] == colLower[iColumn]) {
3706                            nAtLbNatural++;
3707                            sortValue = 0.0;
3708                            double value = dj[iColumn];
3709                            if (value > gap)
3710                                nFix++;
3711                            else if (gap < value*bestRatio)
3712                                bestRatio = gap / value;
3713                            if (saveUpper[iColumn] != colUpper[iColumn]) {
3714                                nNotNeeded++;
3715                                sortValue = 1.0e20;
3716                            }
3717                        } else if (saveUpper[iColumn] == colLower[iColumn]) {
3718                            nAtUbFixed++;
3719                            sortValue = -dj[iColumn];
3720                        } else {
3721                            nAtOther++;
3722                            sortValue = 0.0;
3723                            if (saveLower[iColumn] != colLower[iColumn] &&
3724                                    saveUpper[iColumn] != colUpper[iColumn]) {
3725                                nNotNeeded++;
3726                                sortValue = 1.0e20;
3727                            }
3728                        }
3729                    } else {
3730                        // zero dj
3731                        if (value2 == saveUpper[iColumn]) {
3732                            nAtUbNaturalZero++;
3733                            sortValue = 0.0;
3734                            if (saveLower[iColumn] != colLower[iColumn]) {
3735                                nNotNeeded++;
3736                                sortValue = 1.0e20;
3737                            }
3738                        } else if (value2 == saveLower[iColumn]) {
3739                            nAtLbNaturalZero++;
3740                            sortValue = 0.0;
3741                        } else {
3742                            nAtOtherNatural++;
3743                            sortValue = 0.0;
3744                            if (saveLower[iColumn] != colLower[iColumn] &&
3745                                    saveUpper[iColumn] != colUpper[iColumn]) {
3746                                nNotNeeded++;
3747                                sortValue = 1.0e20;
3748                            }
3749                        }
3750                    }
3751#if HOTSTART==3
3752                    sortValue = -fabs(dj[iColumn]);
3753#endif
3754                }
3755                hotstartSolution_[iColumn] = value ;
3756                saveLower[iColumn] = sortValue;
3757                order[iColumn] = iColumn;
3758            }
3759            COIN_DETAIL_PRINT(printf("** can fix %d columns - best ratio for others is %g on gap of %g\n",
3760                                     nFix, bestRatio, gap));
3761            int nNeg = 0;
3762            CoinSort_2(saveLower, saveLower + numberColumns, order);
3763            for (int i = 0; i < numberColumns; i++) {
3764                if (saveLower[i] < 0.0) {
3765                    nNeg++;
3766#if HOTSTART==2||HOTSTART==3
3767                    // swap sign ?
3768                    saveLower[i] = -saveLower[i];
3769#endif
3770                }
3771            }
3772            CoinSort_2(saveLower, saveLower + nNeg, order);
3773            for (int i = 0; i < numberColumns; i++) {
3774#if HOTSTART==1
3775                hotstartPriorities_[order[i]] = 100;
3776#else
3777                hotstartPriorities_[order[i]] = -(i + 1);
3778#endif
3779            }
3780            COIN_DETAIL_PRINT(printf("nAtLbNat %d,nAtUbNat %d,nAtLbNatZero %d,nAtUbNatZero %d,nAtLbFixed %d,nAtUbFixed %d,nAtOther %d,nAtOtherNat %d, useless %d %d\n",
3781                   nAtLbNatural,
3782                   nAtUbNatural,
3783                   nAtLbNaturalZero,
3784                   nAtUbNaturalZero,
3785                   nAtLbFixed,
3786                   nAtUbFixed,
3787                   nAtOther,
3788                                     nAtOtherNatural, nNotNeeded, nNeg));
3789            delete [] saveLower;
3790            delete [] saveUpper;
3791            if (!saveCompare) {
3792                // create depth first comparison
3793                saveCompare = nodeCompare_;
3794                // depth first
3795                nodeCompare_ = new CbcCompareDepth();
3796                tree_->setComparison(*nodeCompare_) ;
3797            }
3798        }
3799#endif
3800        newNode = new CbcNode ;
3801        // Set objective value (not so obvious if NLP etc)
3802        setObjectiveValue(newNode, NULL);
3803        anyAction = -1 ;
3804        // To make depth available we may need a fake node
3805        CbcNode fakeNode;
3806        if (!currentNode_) {
3807            // Not true if sub trees assert (!numberNodes_);
3808            currentNode_ = &fakeNode;
3809        }
3810        phase_ = 3;
3811        // only allow 1000 passes
3812        int numberPassesLeft = 1000;
3813        // This is first crude step
3814        if (numberAnalyzeIterations_) {
3815            delete [] analyzeResults_;
3816            analyzeResults_ = new double [4*numberIntegers_];
3817            numberFixedAtRoot_ = newNode->analyze(this, analyzeResults_);
3818            if (numberFixedAtRoot_ > 0) {
3819              COIN_DETAIL_PRINT(printf("%d fixed by analysis\n", numberFixedAtRoot_));
3820                setPointers(solver_);
3821                numberFixedNow_ = numberFixedAtRoot_;
3822            } else if (numberFixedAtRoot_ < 0) {
3823              COIN_DETAIL_PRINT(printf("analysis found to be infeasible\n"));
3824                anyAction = -2;
3825                delete newNode ;
3826                newNode = NULL ;
3827                feasible = false ;
3828            }
3829        }
3830        OsiSolverBranch * branches = NULL;
3831        anyAction = chooseBranch(newNode, numberPassesLeft, NULL, cuts, resolved,
3832                                 NULL, NULL, NULL, branches);
3833        if (anyAction == -2 || newNode->objectiveValue() >= cutoff) {
3834            if (anyAction != -2) {
3835                // zap parent nodeInfo
3836#ifdef COIN_DEVELOP
3837                printf("zapping CbcNodeInfo %x\n", reinterpret_cast<int>(newNode->nodeInfo()->parent()));
3838#endif
3839                if (newNode->nodeInfo())
3840                    newNode->nodeInfo()->nullParent();
3841            }
3842            delete newNode ;
3843            newNode = NULL ;
3844            feasible = false ;
3845        }
3846    }
3847    if (newNode && probingInfo_) {
3848        int number01 = probingInfo_->numberIntegers();
3849        //const fixEntry * entry = probingInfo_->fixEntries();
3850        const int * toZero = probingInfo_->toZero();
3851        //const int * toOne = probingInfo_->toOne();
3852        //const int * integerVariable = probingInfo_->integerVariable();
3853        if (toZero[number01]) {
3854            CglTreeProbingInfo info(*probingInfo_);
3855            if ((moreSpecialOptions_&1048576)!=0&&!parentModel_) {
3856              /*
3857                Marginal idea. Further exploration probably good. Build some extra
3858                cliques from probing info. Not quite worth the effort?
3859              */
3860              CglProbing generator1;
3861              generator1.setUsingObjective(false);
3862              generator1.setMaxPass(1);
3863              generator1.setMaxPassRoot(1);
3864              generator1.setMaxLook(100);
3865              generator1.setRowCuts(3);
3866              generator1.setMaxElements(300);
3867              generator1.setMaxProbeRoot(solver_->getNumCols());
3868              CoinThreadRandom randomGenerator;
3869              //CglTreeProbingInfo info(solver_);
3870              info.level = 0;
3871              info.formulation_rows = solver_->getNumRows();
3872              info.inTree = false;
3873              info.randomNumberGenerator=&randomGenerator;
3874              info.pass=4;
3875              generator1.setMode(8);
3876              OsiCuts cs;
3877              generator1.generateCutsAndModify(*solver_,cs,&info);
3878              // very clunky
3879              OsiSolverInterface * temp = generator1.cliqueModel(solver_,2);
3880              CglPreProcess dummy;
3881              OsiSolverInterface * newSolver=dummy.cliqueIt(*temp,0.0001);
3882              delete temp;
3883              OsiSolverInterface * fake = NULL;
3884              if (newSolver) {
3885#if 0
3886                int numberCliques = generator1.numberCliques();
3887                cliqueEntry * entry = generator1.cliqueEntry();
3888                cliqueType * type = new cliqueType [numberCliques];
3889                int * start = new int [numberCliques+1];
3890                start[numberCliques]=2*numberCliques;
3891                int n=0;
3892                for (int i=0;i<numberCliques;i++) {
3893                  start[i]=2*i;
3894                  setOneFixesInCliqueEntry(entry[2*i],true);
3895                  setOneFixesInCliqueEntry(entry[2*i+1],true);
3896                  type[i]=0;
3897                }
3898                fake = info.analyze(*solver_, 1,numberCliques,start,
3899                                    entry,type);
3900                delete [] type;
3901                delete [] entry;
3902#else
3903                fake = info.analyze(*newSolver, 1,-1);
3904#endif
3905                delete newSolver;
3906              } else {
3907                fake = info.analyze(*solver_, 1);
3908              }
3909              if (fake) {
3910                //fake->writeMps("fake");
3911                CglFakeClique cliqueGen(fake);
3912                cliqueGen.setStarCliqueReport(false);
3913                cliqueGen.setRowCliqueReport(false);
3914                cliqueGen.setMinViolation(0.1);
3915                addCutGenerator(&cliqueGen, 1, "Fake cliques", true, false, false, -200);
3916                generator_[numberCutGenerators_-1]->setTiming(true);
3917                for (int i = 0; i < numberCutGenerators_; i++) {
3918                  CglKnapsackCover * cutGen =
3919                  dynamic_cast<CglKnapsackCover *>(generator_[i]->generator());
3920                  if (cutGen) {
3921                    cutGen->createCliques(*fake,2,200,false);
3922                  }
3923                }
3924              }
3925            }
3926            if (probingInfo_->packDown()) {
3927#ifdef CLP_INVESTIGATE
3928                printf("%d implications on %d 0-1\n", toZero[number01], number01);
3929#endif
3930                // Create a cut generator that remembers implications discovered at root.
3931                CglImplication implication(probingInfo_);
3932                addCutGenerator(&implication, 1, "ImplicationCuts", true, false, false, -200);
3933                generator_[numberCutGenerators_-1]->setGlobalCuts(true);
3934            } else {
3935                delete probingInfo_;
3936                probingInfo_ = NULL;
3937            }
3938        } else {
3939            delete probingInfo_;
3940
3941            probingInfo_ = NULL;
3942        }
3943    }
3944    /*
3945      At this point, the root subproblem is infeasible or fathomed by bound
3946      (newNode == NULL), or we're live with an objective value that satisfies the
3947      current objective cutoff.
3948    */
3949    assert (!newNode || newNode->objectiveValue() <= cutoff) ;
3950    // Save address of root node as we don't want to delete it
3951    /*
3952      The common case is that the lp relaxation is feasible but doesn't satisfy
3953      integrality (i.e., newNode->branchingObject(), indicating we've been able to
3954      select a branching variable). Remove any cuts that have gone slack due to
3955      forcing monotone variables. Then tack on an CbcFullNodeInfo object and full
3956      basis (via createInfo()) and stash the new cuts in the nodeInfo (via
3957      addCuts()). If, by some miracle, we have an integral solution at the root
3958      (newNode->branchingObject() is NULL), takeOffCuts() will ensure that the solver holds
3959      a valid solution for use by setBestSolution().
3960    */
3961    CoinWarmStartBasis *lastws = NULL ;
3962    if (feasible && newNode->branchingObject()) {
3963        if (resolved) {
3964            takeOffCuts(cuts, false, NULL) ;
3965#     ifdef CHECK_CUT_COUNTS
3966            { printf("Number of rows after chooseBranch fix (root)"
3967                         "(active only) %d\n",
3968                         numberRowsAtContinuous_ + numberNewCuts_ + numberOldActiveCuts_) ;
3969                const CoinWarmStartBasis* debugws =
3970                    dynamic_cast <const CoinWarmStartBasis*>(solver_->getWarmStart()) ;
3971                debugws->print() ;
3972                delete debugws ;
3973            }
3974#     endif
3975        }
3976        //newNode->createInfo(this,NULL,NULL,NULL,NULL,0,0) ;
3977        //newNode->nodeInfo()->addCuts(cuts,
3978        //                       newNode->numberBranches(),whichGenerator_) ;
3979        if (lastws) delete lastws ;
3980        lastws = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
3981    }
3982    /*
3983      Continuous data to be used later
3984    */
3985    continuousObjective_ = solver_->getObjValue() * solver_->getObjSense();
3986    continuousInfeasibilities_ = 0 ;
3987    if (newNode) {
3988        continuousObjective_ = newNode->objectiveValue() ;
3989        delete [] continuousSolution_;
3990        continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
3991                                              numberColumns);
3992        continuousInfeasibilities_ = newNode->numberUnsatisfied() ;
3993    }
3994    /*
3995      Bound may have changed so reset in objects
3996    */
3997    {
3998        int i ;
3999        for (i = 0; i < numberObjects_; i++)
4000            object_[i]->resetBounds(solver_) ;
4001    }
4002    /*
4003      Feasible? Then we should have either a live node prepped for future
4004      expansion (indicated by variable() >= 0), or (miracle of miracles) an
4005      integral solution at the root node.
4006
4007      initializeInfo sets the reference counts in the nodeInfo object.  Since
4008      this node is still live, push it onto the heap that holds the live set.
4009    */
4010    if (newNode) {
4011        if (newNode->branchingObject()) {
4012            newNode->initializeInfo() ;
4013            tree_->push(newNode) ;
4014            // save pointer to root node - so can pick up bounds
4015            if (!topOfTree_)
4016              topOfTree_ = dynamic_cast<CbcFullNodeInfo *>(newNode->nodeInfo()) ;
4017            if (statistics_) {
4018                if (numberNodes2_ == maximumStatistics_) {
4019                    maximumStatistics_ = 2 * maximumStatistics_;
4020                    CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_];
4021                    memset(temp, 0, maximumStatistics_*sizeof(CbcStatistics *));
4022                    memcpy(temp, statistics_, numberNodes2_*sizeof(CbcStatistics *));
4023                    delete [] statistics_;
4024                    statistics_ = temp;
4025                }
4026                assert (!statistics_[numberNodes2_]);
4027                statistics_[numberNodes2_] = new CbcStatistics(newNode, this);
4028            }
4029            numberNodes2_++;
4030#     ifdef CHECK_NODE
4031            printf("Node %x on tree\n", newNode) ;
4032#     endif
4033        } else {
4034            // continuous is integer
4035            double objectiveValue = newNode->objectiveValue();
4036            setBestSolution(CBC_SOLUTION, objectiveValue,
4037                            solver_->getColSolution()) ;
4038            if (eventHandler) {
4039              // we are stopping anyway so no need to test return code
4040              eventHandler->event(CbcEventHandler::solution);
4041            }
4042            delete newNode ;
4043            newNode = NULL ;
4044        }
4045    }
4046
4047    if (printFrequency_ <= 0) {
4048        printFrequency_ = 1000 ;
4049        if (getNumCols() > 2000)
4050            printFrequency_ = 100 ;
4051    }
4052    /*
4053      It is possible that strong branching fixes one variable and then the code goes round
4054      again and again.  This can take too long.  So we need to warn user - just once.
4055    */
4056    numberLongStrong_ = 0;
4057    CbcNode * createdNode = NULL;
4058#ifdef CBC_THREAD
4059    if ((specialOptions_&2048) != 0)
4060        numberThreads_ = 0;
4061    if (numberThreads_ ) {
4062        nodeCompare_->sayThreaded(); // need to use addresses
4063        master_ = new CbcBaseModel(*this,
4064                                   (parallelMode() < -1) ? 1 : 0);
4065        masterThread_ = master_->masterThread();
4066    }
4067#endif
4068#ifdef COIN_HAS_CLP
4069    {
4070        OsiClpSolverInterface * clpSolver
4071        = dynamic_cast<OsiClpSolverInterface *> (solver_);
4072        if (clpSolver && !parentModel_) {
4073            clpSolver->computeLargestAway();
4074        }
4075    }
4076#endif
4077    /*
4078      At last, the actual branch-and-cut search loop, which will iterate until
4079      the live set is empty or we hit some limit (integrality gap, time, node
4080      count, etc.). The overall flow is to rebuild a subproblem, reoptimise using
4081      solveWithCuts(), choose a branching pattern with chooseBranch(), and finally
4082      add the node to the live set.
4083
4084      The first action is to winnow the live set to remove nodes which are worse
4085      than the current objective cutoff.
4086    */
4087    if (solver_->getRowCutDebuggerAlways()) {
4088        OsiRowCutDebugger * debuggerX = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
4089        const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger() ;
4090        if (!debugger) {
4091            // infeasible!!
4092            printf("before search\n");
4093            const double * lower = solver_->getColLower();
4094            const double * upper = solver_->getColUpper();
4095            const double * solution = debuggerX->optimalSolution();
4096            int numberColumns = solver_->getNumCols();
4097            for (int i = 0; i < numberColumns; i++) {
4098                if (solver_->isInteger(i)) {
4099                    if (solution[i] < lower[i] - 1.0e-6 || solution[i] > upper[i] + 1.0e-6)
4100                        printf("**** ");
4101                    printf("%d %g <= %g <= %g\n",
4102                           i, lower[i], solution[i], upper[i]);
4103                }
4104            }
4105            //abort();
4106        }
4107    }
4108    {
4109        // may be able to change cutoff now
4110        double cutoff = getCutoff();
4111        double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
4112        if (cutoff > bestObjective_ - increment) {
4113            cutoff = bestObjective_ - increment ;
4114            setCutoff(cutoff) ;
4115        }
4116    }
4117#ifdef CBC_THREAD
4118    bool goneParallel = false;
4119#endif
4120#define MAX_DEL_NODE 1
4121    CbcNode * delNode[MAX_DEL_NODE+1];
4122    int nDeleteNode = 0;
4123    // For Printing etc when parallel
4124    int lastEvery1000 = 0;
4125    int lastPrintEvery = 0;
4126    int numberConsecutiveInfeasible = 0;
4127#define PERTURB_IN_FATHOM
4128#ifdef PERTURB_IN_FATHOM
4129    // allow in fathom
4130    if ((moreSpecialOptions_& 262144) != 0)
4131      specialOptions_ |= 131072;
4132#endif
4133    while (true) {
4134        lockThread();
4135#ifdef COIN_HAS_CLP
4136        // See if we want dantzig row choice
4137        goToDantzig(100, savePivotMethod);
4138#endif
4139        if (tree_->empty()) {
4140#ifdef CBC_THREAD
4141            if (parallelMode() > 0 && master_) {
4142                int anyLeft = master_->waitForThreadsInTree(0);
4143                if (!anyLeft) {
4144                    master_->stopThreads(-1);
4145                    break;
4146                }
4147            } else {
4148                break;
4149            }
4150#else
4151            break;
4152#endif
4153        } else {
4154            unlockThread();
4155        }
4156        // If done 100 nodes see if worth trying reduction
4157        if (numberNodes_ == 50 || numberNodes_ == 100) {
4158#ifdef COIN_HAS_CLP
4159            OsiClpSolverInterface * clpSolver
4160            = dynamic_cast<OsiClpSolverInterface *> (solver_);
4161            if (clpSolver && ((specialOptions_&131072) == 0) && true) {
4162                ClpSimplex * simplex = clpSolver->getModelPtr();
4163                int perturbation = simplex->perturbation();
4164#ifdef CLP_INVESTIGATE
4165                printf("Testing its n,s %d %d solves n,s %d %d - pert %d\n",
4166                       numberIterations_, saveNumberIterations,
4167                       numberSolves_, saveNumberSolves, perturbation);
4168#endif
4169                if (perturbation == 50 && (numberIterations_ - saveNumberIterations) <
4170                        8*(numberSolves_ - saveNumberSolves)) {
4171                    // switch off perturbation
4172                    simplex->setPerturbation(100);
4173#ifdef CLP_INVESTIGATE
4174                    printf("Perturbation switched off\n");
4175#endif
4176                }
4177            }
4178#endif
4179            /*
4180              Decide if we want to do a restart.
4181            */
4182            if (saveSolver && (specialOptions_&(512 + 32768)) != 0) {
4183                bool tryNewSearch = solverCharacteristics_->reducedCostsAccurate() &&
4184                                    (getCutoff() < 1.0e20 && getCutoff() < checkCutoffForRestart);
4185                int numberColumns = getNumCols();
4186                if (tryNewSearch) {
4187                    checkCutoffForRestart = getCutoff() ;
4188#ifdef CLP_INVESTIGATE
4189                    printf("after %d nodes, cutoff %g - looking\n",
4190                           numberNodes_, getCutoff());
4191#endif
4192                    saveSolver->resolve();
4193                    double direction = saveSolver->getObjSense() ;
4194                    double gap = checkCutoffForRestart - saveSolver->getObjValue() * direction ;
4195                    double tolerance;
4196                    saveSolver->getDblParam(OsiDualTolerance, tolerance) ;
4197                    if (gap <= 0.0)
4198                        gap = tolerance;
4199                    gap += 100.0 * tolerance;
4200                    double integerTolerance = getDblParam(CbcIntegerTolerance) ;
4201
4202                    const double *lower = saveSolver->getColLower() ;
4203                    const double *upper = saveSolver->getColUpper() ;
4204                    const double *solution = saveSolver->getColSolution() ;
4205                    const double *reducedCost = saveSolver->getReducedCost() ;
4206
4207                    int numberFixed = 0 ;
4208                    int numberFixed2 = 0;
4209#ifdef COIN_DEVELOP
4210                    printf("gap %g\n", gap);
4211#endif
4212                    for (int i = 0 ; i < numberIntegers_ ; i++) {
4213                        int iColumn = integerVariable_[i] ;
4214                        double djValue = direction * reducedCost[iColumn] ;
4215                        if (upper[iColumn] - lower[iColumn] > integerTolerance) {
4216                            if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) {
4217                                //printf("%d to lb on dj of %g - bounds %g %g\n",
4218                                //     iColumn,djValue,lower[iColumn],upper[iColumn]);
4219                                saveSolver->setColUpper(iColumn, lower[iColumn]) ;
4220                                numberFixed++ ;
4221                            } else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) {
4222                                //printf("%d to ub on dj of %g - bounds %g %g\n",
4223                                //     iColumn,djValue,lower[iColumn],upper[iColumn]);
4224                                saveSolver->setColLower(iColumn, upper[iColumn]) ;
4225                                numberFixed++ ;
4226                            }
4227                        } else {
4228                            //printf("%d has dj of %g - already fixed to %g\n",
4229                            //     iColumn,djValue,lower[iColumn]);
4230                            numberFixed2++;
4231                        }
4232                    }
4233#ifdef COIN_DEVELOP
4234                    if ((specialOptions_&1) != 0) {
4235                        const OsiRowCutDebugger *debugger = saveSolver->getRowCutDebugger() ;
4236                        if (debugger) {
4237                            printf("Contains optimal\n") ;
4238                            OsiSolverInterface * temp = saveSolver->clone();
4239                            const double * solution = debugger->optimalSolution();
4240                            const double *lower = temp->getColLower() ;
4241                            const double *upper = temp->getColUpper() ;
4242                            int n = temp->getNumCols();
4243                            for (int i = 0; i < n; i++) {
4244                                if (temp->isInteger(i)) {
4245                                    double value = floor(solution[i] + 0.5);
4246                                    assert (value >= lower[i] && value <= upper[i]);
4247                                    temp->setColLower(i, value);
4248                                    temp->setColUpper(i, value);
4249                                }
4250                            }
4251                            temp->writeMps("reduced_fix");
4252                            delete temp;
4253                            saveSolver->writeMps("reduced");
4254                        } else {
4255                            abort();
4256                        }
4257                    }
4258                    printf("Restart could fix %d integers (%d already fixed)\n",
4259                           numberFixed + numberFixed2, numberFixed2);
4260#endif
4261                    numberFixed += numberFixed2;
4262                    if (numberFixed*10 < numberColumns && numberFixed*4 < numberIntegers_)
4263                        tryNewSearch = false;
4264                }
4265                if (tryNewSearch) {
4266                    // back to solver without cuts?
4267                    OsiSolverInterface * solver2 = saveSolver->clone();
4268                    const double *lower = saveSolver->getColLower() ;
4269                    const double *upper = saveSolver->getColUpper() ;
4270                    for (int i = 0 ; i < numberIntegers_ ; i++) {
4271                        int iColumn = integerVariable_[i] ;
4272                        solver2->setColLower(iColumn, lower[iColumn]);
4273                        solver2->setColUpper(iColumn, upper[iColumn]);
4274                    }
4275                    // swap
4276                    delete saveSolver;
4277                    saveSolver = solver2;
4278                    double * newSolution = new double[numberColumns];
4279                    double objectiveValue = checkCutoffForRestart;
4280                    // Save the best solution so far.
4281                    CbcSerendipity heuristic(*this);
4282                    if (bestSolution_)
4283                        heuristic.setInputSolution(bestSolution_, bestObjective_);
4284                    // Magic number
4285                    heuristic.setFractionSmall(0.8);
4286                    // `pumpTune' to stand-alone solver for explanations.
4287                    heuristic.setFeasibilityPumpOptions(1008013);
4288                    // Use numberNodes to say how many are original rows
4289                    heuristic.setNumberNodes(continuousSolver_->getNumRows());
4290#ifdef COIN_DEVELOP
4291                    if (continuousSolver_->getNumRows() <
4292                            solver_->getNumRows())
4293                        printf("%d rows added ZZZZZ\n",
4294                               solver_->getNumRows() - continuousSolver_->getNumRows());
4295#endif
4296                    int returnCode = heuristic.smallBranchAndBound(saveSolver,
4297                                     -1, newSolution,
4298                                     objectiveValue,
4299                                     checkCutoffForRestart, "Reduce");
4300                    if (returnCode < 0) {
4301#ifdef COIN_DEVELOP
4302                        printf("Restart - not small enough to do search after fixing\n");
4303#endif
4304                        delete [] newSolution;
4305                    } else {
4306                        // 1 for sol'n, 2 for finished, 3 for both
4307                        if ((returnCode&1) != 0) {
4308                            // increment number of solutions so other heuristics can test
4309                            numberSolutions_++;
4310                            numberHeuristicSolutions_++;
4311                            lastHeuristic_ = NULL;
4312                            setBestSolution(CBC_ROUNDING, objectiveValue, newSolution) ;
4313                        }
4314                        delete [] newSolution;
4315#ifdef CBC_THREAD
4316                        if (master_) {
4317                            lockThread();
4318                            if (parallelMode() > 0) {
4319                                while (master_->waitForThreadsInTree(0)) {
4320                                    lockThread();
4321                                    double dummyBest;
4322                                    tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4323                                    //unlockThread();
4324                                }
4325                            } else {
4326                                double dummyBest;
4327                                tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4328                            }
4329                            master_->waitForThreadsInTree(2);
4330                            delete master_;
4331                            master_ = NULL;
4332                            masterThread_ = NULL;
4333                        }
4334#endif
4335                        if (tree_->size()) {
4336                            double dummyBest;
4337                            tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4338                        }
4339                        break;
4340                    }
4341                }
4342                delete saveSolver;
4343                saveSolver = NULL;
4344            }
4345        }
4346        /*
4347          Check for abort on limits: node count, solution count, time, integrality gap.
4348        */
4349        if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
4350                numberSolutions_ < intParam_[CbcMaxNumSol] &&
4351                !maximumSecondsReached() &&
4352                !stoppedOnGap_ && !eventHappened_ && (maximumNumberIterations_ < 0 ||
4353                                                      numberIterations_ < maximumNumberIterations_))) {
4354            // out of loop
4355            break;
4356        }
4357#ifdef BONMIN
4358        assert(!solverCharacteristics_->solutionAddsCuts() || solverCharacteristics_->mipFeasible());
4359#endif
4360// Sets percentage of time when we try diving. Diving requires a bit of heap reorganisation, because
4361// we need to replace the comparison function to dive, and that requires reordering to retain the
4362// heap property.
4363#define DIVE_WHEN 1000
4364#define DIVE_STOP 2000
4365        int kNode = numberNodes_ % 4000;
4366        if (numberNodes_<100000 && kNode>DIVE_WHEN && kNode <= DIVE_STOP) {
4367            if (!parallelMode()) {
4368                if (kNode == DIVE_WHEN + 1 || numberConsecutiveInfeasible > 1) {
4369                    CbcCompareDefault * compare = dynamic_cast<CbcCompareDefault *>
4370                                                  (nodeCompare_);
4371                    // Don't interfere if user has replaced the compare function.
4372                    if (compare) {
4373                        //printf("Redoing tree\n");
4374                        compare->startDive(this);
4375                        numberConsecutiveInfeasible = 0;
4376                    }
4377                }
4378            }
4379        }
4380        // replace current cutoff?
4381        if (cutoff > getCutoff()) {
4382            double newCutoff = getCutoff();
4383            if (analyzeResults_) {
4384                // see if we could fix any (more)
4385                int n = 0;
4386                double * newLower = analyzeResults_;
4387                double * objLower = newLower + numberIntegers_;
4388                double * newUpper = objLower + numberIntegers_;
4389                double * objUpper = newUpper + numberIntegers_;
4390                for (int i = 0; i < numberIntegers_; i++) {
4391                    if (objLower[i] > newCutoff) {
4392                        n++;
4393                        if (objUpper[i] > newCutoff) {
4394                            newCutoff = -COIN_DBL_MAX;
4395                            break;
4396                        }
4397                    } else if (objUpper[i] > newCutoff) {
4398                        n++;
4399                    }
4400                }
4401                if (newCutoff == -COIN_DBL_MAX) {
4402                  COIN_DETAIL_PRINT(printf("Root analysis says finished\n"));
4403                } else if (n > numberFixedNow_) {
4404                  COIN_DETAIL_PRINT(printf("%d more fixed by analysis - now %d\n", n - numberFixedNow_, n));
4405                    numberFixedNow_ = n;
4406                }
4407            }
4408            if (eventHandler) {
4409                if (!eventHandler->event(CbcEventHandler::solution)) {
4410                    eventHappened_ = true; // exit
4411                }
4412                newCutoff = getCutoff();
4413            }
4414            lockThread();
4415            /*
4416              Clean the tree to reflect the new solution, then see if the
4417              node comparison predicate wants to make any changes. If so,
4418              call setComparison for the side effect of rebuilding the heap.
4419            */
4420            tree_->cleanTree(this,newCutoff,bestPossibleObjective_) ;
4421            if (nodeCompare_->newSolution(this) ||
4422                nodeCompare_->newSolution(this,continuousObjective_,
4423                                          continuousInfeasibilities_)) {
4424              tree_->setComparison(*nodeCompare_) ;
4425            }
4426            if (tree_->empty()) {
4427                continue;
4428            }
4429            unlockThread();
4430        }
4431        cutoff = getCutoff() ;
4432        /*
4433            Periodic activities: Opportunities to
4434            + tweak the nodeCompare criteria,
4435            + check if we've closed the integrality gap enough to quit,
4436            + print a summary line to let the user know we're working
4437        */
4438        if (numberNodes_ >= lastEvery1000) {
4439            lockThread();
4440#ifdef COIN_HAS_CLP
4441            // See if we want dantzig row choice
4442            goToDantzig(1000, savePivotMethod);
4443#endif
4444            lastEvery1000 = numberNodes_ + 1000;
4445            bool redoTree = nodeCompare_->every1000Nodes(this, numberNodes_) ;
4446#ifdef CHECK_CUT_SIZE
4447            verifyCutSize (tree_, *this);
4448#endif
4449            // redo tree if requested
4450            if (redoTree)
4451                tree_->setComparison(*nodeCompare_) ;
4452            unlockThread();
4453        }
4454        // Had hotstart before, now switched off
4455        if (saveCompare && !hotstartSolution_) {
4456            // hotstart switched off
4457            delete nodeCompare_; // off depth first
4458            nodeCompare_ = saveCompare;
4459            saveCompare = NULL;
4460            // redo tree
4461            lockThread();
4462            tree_->setComparison(*nodeCompare_) ;
4463            unlockThread();
4464        }
4465        if (numberNodes_ >= lastPrintEvery) {
4466            lastPrintEvery = numberNodes_ + printFrequency_;
4467            lockThread();
4468            int nNodes = tree_->size() ;
4469
4470            //MODIF PIERRE
4471            bestPossibleObjective_ = tree_->getBestPossibleObjective();
4472#ifdef CBC_THREAD
4473            if (parallelMode() > 0 && master_) {
4474              // need to adjust for ones not on tree
4475              int numberThreads = master_->numberThreads();
4476              for (int i=0;i<numberThreads;i++) {
4477                CbcThread * child = master_->child(i);
4478                if (child->node()) {
4479                  // adjust
4480                  double value = child->node()->objectiveValue();
4481                  bestPossibleObjective_ = CoinMin(bestPossibleObjective_, value);
4482                }
4483              }
4484            }
4485#endif
4486            unlockThread();
4487#ifdef CLP_INVESTIGATE
4488            if (getCutoff() < 1.0e20) {
4489                if (fabs(getCutoff() - (bestObjective_ - getCutoffIncrement())) > 1.0e-6 &&
4490                        !parentModel_)
4491                    printf("model cutoff in status %g, best %g, increment %g\n",
4492                           getCutoff(), bestObjective_, getCutoffIncrement());
4493                assert (getCutoff() < bestObjective_ - getCutoffIncrement() +
4494                        1.0e-6 + 1.0e-10*fabs(bestObjective_));
4495            }
4496#endif
4497            if (!intParam_[CbcPrinting]) {
4498                messageHandler()->message(CBC_STATUS, messages())
4499                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4500                << getCurrentSeconds()
4501                << CoinMessageEol ;
4502            } else if (intParam_[CbcPrinting] == 1) {
4503                messageHandler()->message(CBC_STATUS2, messages())
4504                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4505                << tree_->lastDepth() << tree_->lastUnsatisfied() 
4506                << tree_->lastObjective() << numberIterations_
4507                << getCurrentSeconds()
4508                << CoinMessageEol ;
4509            } else if (!numberExtraIterations_) {
4510                messageHandler()->message(CBC_STATUS2, messages())
4511                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4512                << tree_->lastDepth() << tree_->lastUnsatisfied() << numberIterations_
4513                << getCurrentSeconds()
4514                << CoinMessageEol ;
4515            } else {
4516                messageHandler()->message(CBC_STATUS3, messages())
4517                << numberNodes_ << numberExtraNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4518                << tree_->lastDepth() << tree_->lastUnsatisfied() << numberIterations_ << numberExtraIterations_
4519                << getCurrentSeconds()
4520                << CoinMessageEol ;
4521            }
4522            if (eventHandler && !eventHandler->event(CbcEventHandler::treeStatus)) {
4523                eventHappened_ = true; // exit
4524            }
4525        }
4526        // See if can stop on gap
4527        if(canStopOnGap()) {
4528            stoppedOnGap_ = true ;
4529        }
4530
4531#ifdef CHECK_NODE_FULL
4532        verifyTreeNodes(tree_, *this) ;
4533#   endif
4534#   ifdef CHECK_CUT_COUNTS
4535        verifyCutCounts(tree_, *this) ;
4536#   endif
4537        /*
4538          Now we come to the meat of the loop. To create the active subproblem, we'll
4539          pop the most promising node in the live set, rebuild the subproblem it
4540          represents, and then execute the current arm of the branch to create the
4541          active subproblem.
4542        */
4543        CbcNode * node = NULL;
4544#ifdef CBC_THREAD
4545        if (!parallelMode() || parallelMode() == -1) {
4546#endif
4547            node = tree_->bestNode(cutoff) ;
4548            // Possible one on tree worse than cutoff
4549            // Weird comparison function can leave ineligible nodes on tree
4550            if (!node || node->objectiveValue() > cutoff)
4551                continue;
4552            // Do main work of solving node here
4553            doOneNode(this, node, createdNode);
4554#ifdef JJF_ZERO
4555            if (node) {
4556              if (createdNode) {
4557                printf("Node %d depth %d, created %d depth %d\n",
4558                       node->nodeNumber(), node->depth(),
4559                       createdNode->nodeNumber(), createdNode->depth());
4560              } else {
4561                printf("Node %d depth %d,  no created node\n",
4562                       node->nodeNumber(), node->depth());
4563              }
4564            } else if (createdNode) {
4565              printf("Node exhausted, created %d depth %d\n",
4566                     createdNode->nodeNumber(), createdNode->depth());
4567            } else {
4568              printf("Node exhausted,  no created node\n");
4569              numberConsecutiveInfeasible = 2;
4570            }
4571#endif
4572            //if (createdNode)
4573            //numberConsecutiveInfeasible=0;
4574            //else
4575            //numberConsecutiveInfeasible++;
4576#ifdef CBC_THREAD
4577        } else if (parallelMode() > 0) {
4578            //lockThread();
4579            //node = tree_->bestNode(cutoff) ;
4580            // Possible one on tree worse than cutoff
4581            if (true || !node || node->objectiveValue() > cutoff) {
4582                assert (master_);
4583                if (master_) {
4584                    int anyLeft = master_->waitForThreadsInTree(1);
4585                    // may need to go round again
4586                    if (anyLeft) {
4587                        continue;
4588                    } else {
4589                        master_->stopThreads(-1);
4590                    }
4591                }
4592            }
4593            //unlockThread();
4594        } else {
4595            // Deterministic parallel
4596          if ((tree_->size() < CoinMax(numberThreads_, 8)||
4597               hotstartSolution_) && !goneParallel) {
4598                node = tree_->bestNode(cutoff) ;
4599                // Possible one on tree worse than cutoff
4600                if (!node || node->objectiveValue() > cutoff)
4601                    continue;
4602                // Do main work of solving node here
4603                doOneNode(this, node, createdNode);
4604                assert (createdNode);
4605                if (!createdNode->active()) {
4606                    delete createdNode;
4607                    createdNode = NULL;
4608                } else {
4609                    // Say one more pointing to this
4610                    node->nodeInfo()->increment() ;
4611                    tree_->push(createdNode) ;
4612                }
4613                if (node->active()) {
4614                    assert (node->nodeInfo());
4615                    if (node->nodeInfo()->numberBranchesLeft()) {
4616                        tree_->push(node) ;
4617                    } else {
4618                        node->setActive(false);
4619                    }
4620                } else {
4621                    if (node->nodeInfo()) {
4622                        if (!node->nodeInfo()->numberBranchesLeft())
4623                            node->nodeInfo()->allBranchesGone(); // can clean up
4624                        // So will delete underlying stuff
4625                        node->setActive(true);
4626                    }
4627                    delNode[nDeleteNode++] = node;
4628                    node = NULL;
4629                }
4630                if (nDeleteNode >= MAX_DEL_NODE) {
4631                    for (int i = 0; i < nDeleteNode; i++) {
4632                        //printf("trying to del %d %x\n",i,delNode[i]);
4633                        delete delNode[i];
4634                        //printf("done to del %d %x\n",i,delNode[i]);
4635                    }
4636                    nDeleteNode = 0;
4637                }
4638            } else {
4639                // Split and solve
4640                master_->deterministicParallel();
4641                goneParallel = true;
4642            }
4643        }
4644#endif
4645    }
4646    if (nDeleteNode) {
4647        for (int i = 0; i < nDeleteNode; i++) {
4648            delete delNode[i];
4649        }
4650        nDeleteNode = 0;
4651    }
4652#ifdef CBC_THREAD
4653    if (master_) {
4654        master_->stopThreads(-1);
4655        master_->waitForThreadsInTree(2);
4656        delete master_;
4657        master_ = NULL;
4658        masterThread_ = NULL;
4659        // adjust time to allow for children on some systems
4660        //dblParam_[CbcStartSeconds] -= CoinCpuTimeJustChildren();
4661    }
4662#endif
4663    /*
4664      End of the non-abort actions. The next block of code is executed if we've
4665      aborted because we hit one of the limits. Clean up by deleting the live set
4666      and break out of the node processing loop. Note that on an abort, node may
4667      have been pushed back onto the tree for further processing, in which case
4668      it'll be deleted in cleanTree. We need to check.
4669    */
4670    if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
4671          numberSolutions_ < intParam_[CbcMaxNumSol] &&
4672          !maximumSecondsReached() &&
4673          !stoppedOnGap_ &&
4674          !eventHappened_ &&
4675          (maximumNumberIterations_ < 0 || numberIterations_ < maximumNumberIterations_))
4676         ) {
4677        if (tree_->size()) {
4678            double dummyBest;
4679            tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4680        }
4681        delete nextRowCut_;
4682        /* order is important here:
4683         * maximumSecondsReached() should be checked before eventHappened_ and
4684         * isNodeLimitReached() should be checked after eventHappened_
4685         * reason is, that at timelimit, eventHappened_ is set to true to make Cbc stop fast
4686         *   and if Ctrl+C is hit, then the nodelimit is set to -1 to make Cbc stop
4687         */
4688        if (stoppedOnGap_) {
4689            messageHandler()->message(CBC_GAP, messages())
4690            << bestObjective_ - bestPossibleObjective_
4691            << dblParam_[CbcAllowableGap]
4692            << dblParam_[CbcAllowableFractionGap]*100.0
4693            << CoinMessageEol ;
4694            secondaryStatus_ = 2;
4695            status_ = 0 ;
4696        } else if (maximumSecondsReached()) {
4697            handler_->message(CBC_MAXTIME, messages_) << CoinMessageEol ;
4698            secondaryStatus_ = 4;
4699            status_ = 1 ;
4700        } else if (numberSolutions_ >= intParam_[CbcMaxNumSol]) {
4701            handler_->message(CBC_MAXSOLS, messages_) << CoinMessageEol ;
4702            secondaryStatus_ = 6;
4703            status_ = 1 ;
4704        } else if (isNodeLimitReached()) {
4705            handler_->message(CBC_MAXNODES, messages_) << CoinMessageEol ;
4706            secondaryStatus_ = 3;
4707            status_ = 1 ;
4708        } else if (maximumNumberIterations_ >= 0 && numberIterations_ >= maximumNumberIterations_) {
4709            handler_->message(CBC_MAXITERS, messages_) << CoinMessageEol ;
4710            secondaryStatus_ = 8;
4711            status_ = 1 ;
4712        } else {
4713            handler_->message(CBC_EVENT, messages_) << CoinMessageEol ;
4714            secondaryStatus_ = 5;
4715            status_ = 5 ;
4716        }
4717    }
4718    /*
4719      That's it, we've exhausted the search tree, or broken out of the loop because
4720      we hit some limit on evaluation.
4721
4722      We may have got an intelligent tree so give it one more chance
4723    */
4724    // Tell solver we are not in Branch and Cut
4725    solver_->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo, NULL) ;
4726    tree_->endSearch();
4727    //  If we did any sub trees - did we give up on any?
4728    if ( numberStoppedSubTrees_)
4729        status_ = 1;
4730    numberNodes_ += numberExtraNodes_;
4731    numberIterations_ += numberExtraIterations_;
4732    if (eventHandler) {
4733        eventHandler->event(CbcEventHandler::endSearch);
4734    }
4735    if (!status_) {
4736        // Set best possible unless stopped on gap
4737        if (secondaryStatus_ != 2)
4738            bestPossibleObjective_ = bestObjective_;
4739        handler_->message(CBC_END_GOOD, messages_)
4740        << bestObjective_ << numberIterations_ << numberNodes_ << getCurrentSeconds()
4741        << CoinMessageEol ;
4742    } else {
4743        handler_->message(CBC_END, messages_)
4744        << bestObjective_ << bestPossibleObjective_
4745        << numberIterations_ << numberNodes_ << getCurrentSeconds()
4746        << CoinMessageEol ;
4747    }
4748    if (numberStrongIterations_)
4749        handler_->message(CBC_STRONG_STATS, messages_)
4750        << strongInfo_[0] << numberStrongIterations_ << strongInfo_[2]
4751        << strongInfo_[1] << CoinMessageEol ;
4752    if (!numberExtraNodes_)
4753        handler_->message(CBC_OTHER_STATS, messages_)
4754        << maximumDepthActual_
4755        << numberDJFixed_ << CoinMessageEol ;
4756    else
4757        handler_->message(CBC_OTHER_STATS2, messages_)
4758        << maximumDepthActual_
4759        << numberDJFixed_ << numberExtraNodes_ << numberExtraIterations_
4760        << CoinMessageEol ;
4761    if (doStatistics == 100) {
4762        for (int i = 0; i < numberObjects_; i++) {
4763            CbcSimpleIntegerDynamicPseudoCost * obj =
4764                dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object_[i]) ;
4765            if (obj)
4766                obj->print();
4767        }
4768    }
4769    if (statistics_) {
4770        // report in some way
4771        int * lookup = new int[numberObjects_];
4772        int i;
4773        for (i = 0; i < numberObjects_; i++)
4774            lookup[i] = -1;
4775        bool goodIds = false; //true;
4776        for (i = 0; i < numberObjects_; i++) {
4777            int iColumn = object_[i]->columnNumber();
4778            if (iColumn >= 0 && iColumn < numberColumns) {
4779                if (lookup[i] == -1) {
4780                    lookup[i] = iColumn;
4781                } else {
4782                    goodIds = false;
4783                    break;
4784                }
4785            } else {
4786                goodIds = false;
4787                break;
4788            }
4789        }
4790        if (!goodIds) {
4791            delete [] lookup;
4792            lookup = NULL;
4793        }
4794        if (doStatistics >= 3) {
4795            printf("  node parent depth column   value                    obj      inf\n");
4796            for ( i = 0; i < numberNodes2_; i++) {
4797                statistics_[i]->print(lookup);
4798            }
4799        }
4800        if (doStatistics > 1) {
4801            // Find last solution
4802            int k;
4803            for (k = numberNodes2_ - 1; k >= 0; k--) {
4804                if (statistics_[k]->endingObjective() != COIN_DBL_MAX &&
4805                        !statistics_[k]->endingInfeasibility())
4806                    break;
4807            }
4808            if (k >= 0) {
4809                int depth = statistics_[k]->depth();
4810                int * which = new int[depth+1];
4811                for (i = depth; i >= 0; i--) {
4812                    which[i] = k;
4813                    k = statistics_[k]->parentNode();
4814                }
4815                printf("  node parent depth column   value                    obj      inf\n");
4816                for (i = 0; i <= depth; i++) {
4817                    statistics_[which[i]]->print(lookup);
4818                }
4819                delete [] which;
4820            }
4821        }
4822        // now summary
4823        int maxDepth = 0;
4824        double averageSolutionDepth = 0.0;
4825        int numberSolutions = 0;
4826        double averageCutoffDepth = 0.0;
4827        double averageSolvedDepth = 0.0;
4828        int numberCutoff = 0;
4829        int numberDown = 0;
4830        int numberFirstDown = 0;
4831        double averageInfDown = 0.0;
4832        double averageObjDown = 0.0;
4833        int numberCutoffDown = 0;
4834        int numberUp = 0;
4835        int numberFirstUp = 0;
4836        double averageInfUp = 0.0;
4837        double averageObjUp = 0.0;
4838        int numberCutoffUp = 0;
4839        double averageNumberIterations1 = 0.0;
4840        double averageValue = 0.0;
4841        for ( i = 0; i < numberNodes2_; i++) {
4842            int depth =  statistics_[i]->depth();
4843            int way =  statistics_[i]->way();
4844            double value = statistics_[i]->value();
4845            double startingObjective =  statistics_[i]->startingObjective();
4846            int startingInfeasibility = statistics_[i]->startingInfeasibility();
4847            double endingObjective = statistics_[i]->endingObjective();
4848            int endingInfeasibility = statistics_[i]->endingInfeasibility();
4849            maxDepth = CoinMax(depth, maxDepth);
4850            // Only for completed
4851            averageNumberIterations1 += statistics_[i]->numberIterations();
4852            averageValue += value;
4853            if (endingObjective != COIN_DBL_MAX && !endingInfeasibility) {
4854                numberSolutions++;
4855                averageSolutionDepth += depth;
4856            }
4857            if (endingObjective == COIN_DBL_MAX) {
4858                numberCutoff++;
4859                averageCutoffDepth += depth;
4860                if (way < 0) {
4861                    numberDown++;
4862                    numberCutoffDown++;
4863                    if (way == -1)
4864                        numberFirstDown++;
4865                } else {
4866                    numberUp++;
4867                    numberCutoffUp++;
4868                    if (way == 1)
4869                        numberFirstUp++;
4870                }
4871            } else {
4872                averageSolvedDepth += depth;
4873                if (way < 0) {
4874                    numberDown++;
4875                    averageInfDown += startingInfeasibility - endingInfeasibility;
4876                    averageObjDown += endingObjective - startingObjective;
4877                    if (way == -1)
4878                        numberFirstDown++;
4879                } else {
4880                    numberUp++;
4881                    averageInfUp += startingInfeasibility - endingInfeasibility;
4882                    averageObjUp += endingObjective - startingObjective;
4883                    if (way == 1)
4884                        numberFirstUp++;
4885                }
4886            }
4887        }
4888        // Now print
4889        if (numberSolutions)
4890            averageSolutionDepth /= static_cast<double> (numberSolutions);
4891        int numberSolved = numberNodes2_ - numberCutoff;
4892        double averageNumberIterations2 = numberIterations_ - averageNumberIterations1
4893                                          - numberIterationsAtContinuous;
4894        if (numberCutoff) {
4895            averageCutoffDepth /= static_cast<double> (numberCutoff);
4896            averageNumberIterations2 /= static_cast<double> (numberCutoff);
4897        }
4898        if (numberNodes2_)
4899            averageValue /= static_cast<double> (numberNodes2_);
4900        if (numberSolved) {
4901            averageNumberIterations1 /= static_cast<double> (numberSolved);
4902            averageSolvedDepth /= static_cast<double> (numberSolved);
4903        }
4904        printf("%d solution(s) were found (by branching) at an average depth of %g\n",
4905               numberSolutions, averageSolutionDepth);
4906        printf("average value of variable being branched on was %g\n",
4907               averageValue);
4908        printf("%d nodes were cutoff at an average depth of %g with iteration count of %g\n",
4909               numberCutoff, averageCutoffDepth, averageNumberIterations2);
4910        printf("%d nodes were solved at an average depth of %g with iteration count of %g\n",
4911               numberSolved, averageSolvedDepth, averageNumberIterations1);
4912        if (numberDown) {
4913            averageInfDown /= static_cast<double> (numberDown);
4914            averageObjDown /= static_cast<double> (numberDown);
4915        }
4916        printf("Down %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
4917               numberDown, numberFirstDown, numberDown - numberFirstDown, numberCutoffDown,
4918               averageInfDown, averageObjDown);
4919        if (numberUp) {
4920            averageInfUp /= static_cast<double> (numberUp);
4921            averageObjUp /= static_cast<double> (numberUp);
4922        }
4923        printf("Up %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
4924               numberUp, numberFirstUp, numberUp - numberFirstUp, numberCutoffUp,
4925               averageInfUp, averageObjUp);
4926        for ( i = 0; i < numberNodes2_; i++)
4927            delete statistics_[i];
4928        delete [] statistics_;
4929        statistics_ = NULL;
4930        maximumStatistics_ = 0;
4931        delete [] lookup;
4932    }
4933    /*
4934      If we think we have a solution, restore and confirm it with a call to
4935      setBestSolution().  We need to reset the cutoff value so as not to fathom
4936      the solution on bounds.  Note that calling setBestSolution( ..., true)
4937      leaves the continuousSolver_ bounds vectors fixed at the solution value.
4938
4939      Running resolve() here is a failsafe --- setBestSolution has already
4940      reoptimised using the continuousSolver_. If for some reason we fail to
4941      prove optimality, run the problem again after instructing the solver to
4942      tell us more.
4943
4944      If all looks good, replace solver_ with continuousSolver_, so that the
4945      outside world will be able to obtain information about the solution using
4946      public methods.
4947
4948      Don't replace if we are trying to save cuts
4949    */
4950    if (bestSolution_ && (solverCharacteristics_->solverType() < 2 || solverCharacteristics_->solverType() == 4) && 
4951        ((specialOptions_&8388608)==0||(specialOptions_&2048)!=0)) {
4952        setCutoff(1.0e50) ; // As best solution should be worse than cutoff
4953        // change cutoff as constraint if wanted
4954        if (cutoffRowNumber_>=0) {
4955          if (solver_->getNumRows()>cutoffRowNumber_)
4956            solver_->setRowUpper(cutoffRowNumber_,1.0e50);
4957        }
4958        // also in continuousSolver_
4959        if (continuousSolver_) {
4960          // Solvers know about direction
4961          double direction = solver_->getObjSense();
4962          continuousSolver_->setDblParam(OsiDualObjectiveLimit, 1.0e50*direction);
4963        }
4964        phase_ = 5;
4965        double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
4966        if ((specialOptions_&4) == 0)
4967            bestObjective_ += 100.0 * increment + 1.0e-3; // only set if we are going to solve
4968        setBestSolution(CBC_END_SOLUTION, bestObjective_, bestSolution_, 1) ;
4969        continuousSolver_->resolve() ;
4970        if (!continuousSolver_->isProvenOptimal()) {
4971            continuousSolver_->messageHandler()->setLogLevel(2) ;
4972            continuousSolver_->initialSolve() ;
4973        }
4974        delete solver_ ;
4975        // above deletes solverCharacteristics_
4976        solverCharacteristics_ = NULL;
4977        solver_ = continuousSolver_ ;
4978        setPointers(solver_);
4979        continuousSolver_ = NULL ;
4980    }
4981    /*
4982      Clean up dangling objects. continuousSolver_ may already be toast.
4983    */
4984    delete lastws ;
4985    if (saveObjects) {
4986        for (int i = 0; i < numberObjects_; i++)
4987            delete saveObjects[i];
4988        delete [] saveObjects;
4989    }
4990    numberStrong_ = saveNumberStrong;
4991    numberBeforeTrust_ = saveNumberBeforeTrust;
4992    delete [] whichGenerator_ ;
4993    whichGenerator_ = NULL;
4994    delete [] lowerBefore ;
4995    delete [] upperBefore ;
4996    delete [] walkback_ ;
4997    walkback_ = NULL ;
4998    delete [] lastNodeInfo_ ;
4999    lastNodeInfo_ = NULL;
5000    delete [] lastNumberCuts_ ;
5001    lastNumberCuts_ = NULL;
5002    delete [] lastCut_;
5003    lastCut_ = NULL;
5004    delete [] addedCuts_ ;
5005    addedCuts_ = NULL ;
5006    //delete persistentInfo;
5007    // Get rid of characteristics
5008    solverCharacteristics_ = NULL;
5009    if (continuousSolver_) {
5010        delete continuousSolver_ ;
5011        continuousSolver_ = NULL ;
5012    }
5013    /*
5014      Destroy global cuts by replacing with an empty OsiCuts object.
5015    */
5016    globalCuts_ = CbcRowCuts() ;
5017    delete globalConflictCuts_;
5018    globalConflictCuts_=NULL;
5019    if (!bestSolution_ && (specialOptions_&8388608)==0) {
5020        // make sure lp solver is infeasible
5021        int numberColumns = solver_->getNumCols();
5022        const double * columnLower = solver_->getColLower();
5023        int iColumn;
5024        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
5025            if (solver_->isInteger(iColumn))
5026                solver_->setColUpper(iColumn, columnLower[iColumn]);
5027        }
5028        solver_->initialSolve();
5029    }
5030#ifdef COIN_HAS_CLP
5031    {
5032        OsiClpSolverInterface * clpSolver
5033        = dynamic_cast<OsiClpSolverInterface *> (solver_);
5034        if (clpSolver) {
5035            // Possible restore of pivot method
5036            if (savePivotMethod) {
5037                // model may have changed
5038                savePivotMethod->setModel(NULL);
5039                clpSolver->getModelPtr()->setDualRowPivotAlgorithm(*savePivotMethod);
5040                delete savePivotMethod;
5041            }
5042            clpSolver->setLargestAway(-1.0);
5043        }
5044    }
5045#endif
5046    if ((fastNodeDepth_ >= 1000000 || (moreSpecialOptions_&33554432)!=0)
5047         && !parentModel_) {
5048      // delete object off end
5049      delete object_[numberObjects_];
5050      if ((moreSpecialOptions_&33554432)==0)
5051        fastNodeDepth_ -= 1000000;
5052    }
5053    delete saveSolver;
5054    // Undo preprocessing performed during BaB.
5055    if (strategy_ && strategy_->preProcessState() > 0) {
5056        // undo preprocessing
5057        CglPreProcess * process = strategy_->process();
5058        assert (process);
5059        int n = originalSolver->getNumCols();
5060        if (bestSolution_) {
5061            delete [] bestSolution_;
5062            bestSolution_ = new double [n];
5063            process->postProcess(*solver_);
5064        }
5065        strategy_->deletePreProcess();
5066        // Solution now back in originalSolver
5067        delete solver_;
5068        solver_ = originalSolver;
5069        if (bestSolution_) {
5070            bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
5071            memcpy(bestSolution_, solver_->getColSolution(), n*sizeof(double));
5072        }
5073        // put back original objects if there were any
5074        if (originalObject) {
5075            int iColumn;
5076            assert (ownObjects_);
5077            for (iColumn = 0; iColumn < numberObjects_; iColumn++)
5078                delete object_[iColumn];
5079            delete [] object_;
5080            numberObjects_ = numberOriginalObjects;
5081            object_ = originalObject;
5082            delete [] integerVariable_;
5083            numberIntegers_ = 0;
5084            for (iColumn = 0; iColumn < n; iColumn++) {
5085                if (solver_->isInteger(iColumn))
5086                    numberIntegers_++;
5087            }
5088            integerVariable_ = new int[numberIntegers_];
5089            numberIntegers_ = 0;
5090            for (iColumn = 0; iColumn < n; iColumn++) {
5091                if (solver_->isInteger(iColumn))
5092                    integerVariable_[numberIntegers_++] = iColumn;
5093            }
5094        }
5095    }
5096    if (flipObjective)
5097      flipModel();
5098#ifdef COIN_HAS_CLP
5099    {
5100        OsiClpSolverInterface * clpSolver
5101        = dynamic_cast<OsiClpSolverInterface *> (solver_);
5102        if (clpSolver)
5103            clpSolver->setFakeObjective(reinterpret_cast<double *> (NULL));
5104    }
5105#endif
5106    moreSpecialOptions_ = saveMoreSpecialOptions;
5107    return ;
5108}
5109
5110
5111// Solve the initial LP relaxation
5112void
5113CbcModel::initialSolve()
5114{
5115    assert (solver_);
5116    // Double check optimization directions line up
5117    dblParam_[CbcOptimizationDirection] = solver_->getObjSense();
5118    // Check if bounds are all integral (as may get messed up later)
5119    checkModel();
5120    if (!solverCharacteristics_) {
5121        OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
5122        if (solverCharacteristics) {
5123            solverCharacteristics_ = solverCharacteristics;
5124        } else {
5125            // replace in solver
5126            OsiBabSolver defaultC;
5127            solver_->setAuxiliaryInfo(&defaultC);
5128            solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
5129        }
5130    }
5131    solverCharacteristics_->setSolver(solver_);
5132    solver_->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
5133    solver_->initialSolve();
5134    solver_->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo, NULL) ;
5135    if (!solver_->isProvenOptimal())
5136        solver_->resolve();
5137    // But set up so Jon Lee will be happy
5138    status_ = -1;
5139    secondaryStatus_ = -1;
5140    originalContinuousObjective_ = solver_->getObjValue() * solver_->getObjSense();
5141    bestPossibleObjective_ = originalContinuousObjective_;
5142    if (solver_->isProvenDualInfeasible())
5143      originalContinuousObjective_ = -COIN_DBL_MAX;
5144    delete [] continuousSolution_;
5145    continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
5146                                          solver_->getNumCols());
5147    setPointers(solver_);
5148    solverCharacteristics_ = NULL;
5149}
5150
5151/*! \brief Get an empty basis object
5152
5153  Return an empty CoinWarmStartBasis object with the requested capacity,
5154  appropriate for the current solver. The object is cloned from the object
5155  cached as emptyWarmStart_. If there is no cached object, the routine
5156  queries the solver for a warm start object, empties it, and caches the
5157  result.
5158*/
5159
5160CoinWarmStartBasis *CbcModel::getEmptyBasis (int ns, int na) const
5161
5162{
5163    CoinWarmStartBasis *emptyBasis ;
5164    /*
5165      Acquire an empty basis object, if we don't yet have one.
5166    */
5167    if (emptyWarmStart_ == 0) {
5168        if (solver_ == 0) {
5169            throw CoinError("Cannot construct basis without solver!",
5170                            "getEmptyBasis", "CbcModel") ;
5171        }
5172        emptyBasis =
5173            dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart()) ;
5174        if (emptyBasis == 0) {
5175            throw CoinError(
5176                "Solver does not appear to use a basis-oriented warm start.",
5177                "getEmptyBasis", "CbcModel") ;
5178        }
5179        emptyBasis->setSize(0, 0) ;
5180        emptyWarmStart_ = dynamic_cast<CoinWarmStart *>(emptyBasis) ;
5181    }
5182    /*
5183      Clone the empty basis object, resize it as requested, and return.
5184    */
5185    emptyBasis = dynamic_cast<CoinWarmStartBasis *>(emptyWarmStart_->clone()) ;
5186    assert(emptyBasis) ;
5187    if (ns != 0 || na != 0) emptyBasis->setSize(ns, na) ;
5188
5189    return (emptyBasis) ;
5190}
5191
5192
5193/** Default Constructor
5194
5195  Creates an empty model without an associated solver.
5196*/
5197CbcModel::CbcModel()
5198
5199        :
5200        solver_(NULL),
5201        ownership_(0x80000000),
5202        continuousSolver_(NULL),
5203        referenceSolver_(NULL),
5204        defaultHandler_(true),
5205        emptyWarmStart_(NULL),
5206        bestObjective_(COIN_DBL_MAX),
5207        bestPossibleObjective_(COIN_DBL_MAX),
5208        sumChangeObjective1_(0.0),
5209        sumChangeObjective2_(0.0),
5210        bestSolution_(NULL),
5211        savedSolutions_(NULL),
5212        currentSolution_(NULL),
5213        testSolution_(NULL),
5214        globalConflictCuts_(NULL),
5215        minimumDrop_(1.0e-4),
5216        numberSolutions_(0),
5217        numberSavedSolutions_(0),
5218        maximumSavedSolutions_(0),
5219        stateOfSearch_(0),
5220        whenCuts_(-1),
5221        hotstartSolution_(NULL),
5222        hotstartPriorities_(NULL),
5223        numberHeuristicSolutions_(0),
5224        numberNodes_(0),
5225        numberNodes2_(0),
5226        numberIterations_(0),
5227        numberSolves_(0),
5228        status_(-1),
5229        secondaryStatus_(-1),
5230        numberIntegers_(0),
5231        numberRowsAtContinuous_(0),
5232        cutoffRowNumber_(-1),
5233        maximumNumberCuts_(0),
5234        phase_(0),
5235        currentNumberCuts_(0),
5236        maximumDepth_(0),
5237        walkback_(NULL),
5238        lastNodeInfo_(NULL),
5239        lastCut_(NULL),
5240        lastDepth_(0),
5241        lastNumberCuts2_(0