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

Last change on this file since 1450 was 1450, checked in by forrest, 9 years ago

for BonMin?

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