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

Last change on this file since 1844 was 1844, checked in by tkr, 6 years ago

Fixing forward declaration of static function

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