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

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

changes for cuts and extra variables

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