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

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

allow doRootThread when serial

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