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

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

move bad static definition

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