source: stable/2.9/Cbc/src/CbcModel.cpp @ 2541

Last change on this file since 2541 was 2541, checked in by forrest, 8 months ago

set currentNode_ NULL

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