source: stable/2.8/Cbc/src/CbcModel.cpp @ 1856

Last change on this file since 1856 was 1845, checked in by stefan, 7 years ago

remove duplicate function declaration

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