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

Last change on this file since 2180 was 2165, checked in by forrest, 5 years ago

try and fix segfault with threads when exiting

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 753.5 KB
Line 
1/* $Id: CbcModel.cpp 2165 2015-03-17 08:58:35Z 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) var %d way %d obj %g", j, node,
271               node->nodeInfo(), node->columnNumber(), 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    // when to check for restart
1644    int nextCheckRestart=50;
1645    // Force minimization !!!!
1646    bool flipObjective = (solver_->getObjSense()<0.0);
1647    if (flipObjective)
1648      flipModel();
1649    dblParam_[CbcOptimizationDirection] = 1.0; // was solver_->getObjSense();
1650    strongInfo_[0] = 0;
1651    strongInfo_[1] = 0;
1652    strongInfo_[2] = 0;
1653    strongInfo_[3] = 0;
1654    strongInfo_[4] = 0;
1655    strongInfo_[5] = 0;
1656    strongInfo_[6] = 0;
1657    numberStrongIterations_ = 0;
1658    currentNode_ = NULL;
1659    // See if should do cuts old way
1660    if (parallelMode() < 0) {
1661        specialOptions_ |= 4096 + 8192;
1662    } else if (parallelMode() > 0) {
1663        specialOptions_ |= 4096;
1664    }
1665    int saveMoreSpecialOptions = moreSpecialOptions_;
1666    if (dynamic_cast<CbcTreeLocal *> (tree_))
1667        specialOptions_ |= 4096 + 8192;
1668#ifdef COIN_HAS_CLP
1669    {
1670        OsiClpSolverInterface * clpSolver
1671        = dynamic_cast<OsiClpSolverInterface *> (solver_);
1672        if (clpSolver) {
1673            // pass in disaster handler
1674            CbcDisasterHandler handler(this);
1675            clpSolver->passInDisasterHandler(&handler);
1676            // Initialise solvers seed (unless users says not)
1677            if ((specialOptions_&4194304)==0)
1678              clpSolver->getModelPtr()->setRandomSeed(1234567);
1679#ifdef JJF_ZERO
1680            // reduce factorization frequency
1681            int frequency = clpSolver->getModelPtr()->factorizationFrequency();
1682            clpSolver->getModelPtr()->setFactorizationFrequency(CoinMin(frequency, 120));
1683#endif
1684        }
1685    }
1686#endif
1687    // original solver (only set if pre-processing)
1688    OsiSolverInterface * originalSolver = NULL;
1689    int numberOriginalObjects = numberObjects_;
1690    OsiObject ** originalObject = NULL;
1691    // Save whether there were any objects
1692    bool noObjects = (numberObjects_ == 0);
1693    // Set up strategies
1694    /*
1695      See if the user has supplied a strategy object and deal with it if present.
1696      The call to setupOther will set numberStrong_ and numberBeforeTrust_, and
1697      perform integer preprocessing, if requested.
1698
1699      We need to hang on to a pointer to solver_. setupOther will assign a
1700      preprocessed solver to model, but will instruct assignSolver not to trash the
1701      existing one.
1702    */
1703    if (strategy_) {
1704        // May do preprocessing
1705        originalSolver = solver_;
1706        strategy_->setupOther(*this);
1707        if (strategy_->preProcessState()) {
1708            // pre-processing done
1709            if (strategy_->preProcessState() < 0) {
1710                // infeasible (or unbounded)
1711                status_ = 0 ;
1712                if (!solver_->isProvenDualInfeasible()) {
1713                  handler_->message(CBC_INFEAS, messages_) << CoinMessageEol ;
1714                  secondaryStatus_ = 1;
1715                } else {
1716                  handler_->message(CBC_UNBOUNDED, 
1717                                    messages_) << CoinMessageEol ;
1718                  secondaryStatus_ = 7;
1719                }
1720                originalContinuousObjective_ = COIN_DBL_MAX;
1721                if (flipObjective)
1722                  flipModel();
1723                return ;
1724            } else if (numberObjects_ && object_) {
1725                numberOriginalObjects = numberObjects_;
1726                // redo sequence
1727                numberIntegers_ = 0;
1728                int numberColumns = getNumCols();
1729                int nOrig = originalSolver->getNumCols();
1730                CglPreProcess * process = strategy_->process();
1731                assert (process);
1732                const int * originalColumns = process->originalColumns();
1733                // allow for cliques etc
1734                nOrig = CoinMax(nOrig, originalColumns[numberColumns-1] + 1);
1735                // try and redo debugger
1736                OsiRowCutDebugger * debugger = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
1737                if (debugger) {
1738                  if (numberColumns<=debugger->numberColumns())
1739                    debugger->redoSolution(numberColumns, originalColumns);
1740                  else
1741                    debugger=NULL; // no idea how to handle (SOS?)
1742                }
1743                // User-provided solution might have been best. Synchronise.
1744                if (bestSolution_) {
1745                    // need to redo - in case no better found in BAB
1746                    // just get integer part right
1747                    for (int i = 0; i < numberColumns; i++) {
1748                        int jColumn = originalColumns[i];
1749                        bestSolution_[i] = bestSolution_[jColumn];
1750                    }
1751                }
1752                originalObject = object_;
1753                // object number or -1
1754                int * temp = new int[nOrig];
1755                int iColumn;
1756                for (iColumn = 0; iColumn < nOrig; iColumn++)
1757                    temp[iColumn] = -1;
1758                int iObject;
1759                int nNonInt = 0;
1760                for (iObject = 0; iObject < numberOriginalObjects; iObject++) {
1761                    iColumn = originalObject[iObject]->columnNumber();
1762                    if (iColumn < 0) {
1763                        nNonInt++;
1764                    } else {
1765                        temp[iColumn] = iObject;
1766                    }
1767                }
1768                int numberNewIntegers = 0;
1769                int numberOldIntegers = 0;
1770                int numberOldOther = 0;
1771                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1772                    int jColumn = originalColumns[iColumn];
1773                    if (temp[jColumn] >= 0) {
1774                        int iObject = temp[jColumn];
1775                        CbcSimpleInteger * obj =
1776                            dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ;
1777                        if (obj)
1778                            numberOldIntegers++;
1779                        else
1780                            numberOldOther++;
1781                    } else if (isInteger(iColumn)) {
1782                        numberNewIntegers++;
1783                    }
1784                }
1785                /*
1786                  Allocate an array to hold the indices of the integer variables.
1787                  Make a large enough array for all objects
1788                */
1789                numberObjects_ = numberNewIntegers + numberOldIntegers + numberOldOther + nNonInt;
1790                object_ = new OsiObject * [numberObjects_];
1791                delete [] integerVariable_;
1792                integerVariable_ = new int [numberNewIntegers+numberOldIntegers];
1793                /*
1794                  Walk the variables again, filling in the indices and creating objects for
1795                  the integer variables. Initially, the objects hold the index and upper &
1796                  lower bounds.
1797                */
1798                numberIntegers_ = 0;
1799                int n = originalColumns[numberColumns-1] + 1;
1800                int * backward = new int[n];
1801                int i;
1802                for ( i = 0; i < n; i++)
1803                    backward[i] = -1;
1804                for (i = 0; i < numberColumns; i++)
1805                    backward[originalColumns[i]] = i;
1806                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1807                    int jColumn = originalColumns[iColumn];
1808                    if (temp[jColumn] >= 0) {
1809                        int iObject = temp[jColumn];
1810                        CbcSimpleInteger * obj =
1811                            dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ;
1812                        if (obj) {
1813                            object_[numberIntegers_] = originalObject[iObject]->clone();
1814                            // redo ids etc
1815                            //object_[numberIntegers_]->resetSequenceEtc(numberColumns,originalColumns);
1816                            object_[numberIntegers_]->resetSequenceEtc(numberColumns, backward);
1817                            integerVariable_[numberIntegers_++] = iColumn;
1818                        }
1819                    } else if (isInteger(iColumn)) {
1820                        object_[numberIntegers_] =
1821                            new CbcSimpleInteger(this, iColumn);
1822                        integerVariable_[numberIntegers_++] = iColumn;
1823                    }
1824                }
1825                delete [] backward;
1826                numberObjects_ = numberIntegers_;
1827                // Now append other column stuff
1828                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1829                    int jColumn = originalColumns[iColumn];
1830                    if (temp[jColumn] >= 0) {
1831                        int iObject = temp[jColumn];
1832                        CbcSimpleInteger * obj =
1833                            dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ;
1834                        if (!obj) {
1835                            object_[numberObjects_] = originalObject[iObject]->clone();
1836                            // redo ids etc
1837                            CbcObject * obj =
1838                                dynamic_cast <CbcObject *>(object_[numberObjects_]) ;
1839                            assert (obj);
1840                            obj->redoSequenceEtc(this, numberColumns, originalColumns);
1841                            numberObjects_++;
1842                        }
1843                    }
1844                }
1845                // now append non column stuff
1846                for (iObject = 0; iObject < numberOriginalObjects; iObject++) {
1847                    iColumn = originalObject[iObject]->columnNumber();
1848                    if (iColumn < 0) {
1849                        // already has column numbers changed
1850                        object_[numberObjects_] = originalObject[iObject]->clone();
1851#ifdef JJF_ZERO
1852                        // redo ids etc
1853                        CbcObject * obj =
1854                            dynamic_cast <CbcObject *>(object_[numberObjects_]) ;
1855                        assert (obj);
1856                        obj->redoSequenceEtc(this, numberColumns, originalColumns);
1857#endif
1858                        numberObjects_++;
1859                    }
1860                }
1861                delete [] temp;
1862                if (!numberObjects_)
1863                    handler_->message(CBC_NOINT, messages_) << CoinMessageEol ;
1864            } else {
1865                int numberColumns = getNumCols();
1866                CglPreProcess * process = strategy_->process();
1867                assert (process);
1868                const int * originalColumns = process->originalColumns();
1869                // try and redo debugger
1870                OsiRowCutDebugger * debugger = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
1871                if (debugger)
1872                    debugger->redoSolution(numberColumns, originalColumns);
1873            }
1874        } else {
1875            //no preprocessing
1876            originalSolver = NULL;
1877        }
1878        strategy_->setupCutGenerators(*this);
1879        strategy_->setupHeuristics(*this);
1880        // Set strategy print level to models
1881        strategy_->setupPrinting(*this, handler_->logLevel());
1882    }
1883    eventHappened_ = false;
1884    CbcEventHandler *eventHandler = getEventHandler() ;
1885    if (eventHandler)
1886        eventHandler->setModel(this);
1887#define CLIQUE_ANALYSIS
1888#ifdef CLIQUE_ANALYSIS
1889    // set up for probing
1890    // If we're doing clever stuff with cliques, additional info here.
1891    if (!parentModel_)
1892        probingInfo_ = new CglTreeProbingInfo(solver_);
1893    else
1894        probingInfo_ = NULL;
1895#else
1896    probingInfo_ = NULL;
1897#endif
1898
1899    // Try for dominated columns
1900    if ((specialOptions_&64) != 0) {
1901        CglDuplicateRow dupcuts(solver_);
1902        dupcuts.setMode(2);
1903        CglStored * storedCuts = dupcuts.outDuplicates(solver_);
1904        if (storedCuts) {
1905          COIN_DETAIL_PRINT(printf("adding dup cuts\n"));
1906            addCutGenerator(storedCuts, 1, "StoredCuts from dominated",
1907                            true, false, false, -200);
1908        }
1909    }
1910    if (!nodeCompare_)
1911        nodeCompare_ = new CbcCompareDefault();;
1912    // See if hot start wanted
1913    CbcCompareBase * saveCompare = NULL;
1914    // User supplied hotstart. Adapt for preprocessing.
1915    if (hotstartSolution_) {
1916        if (strategy_ && strategy_->preProcessState() > 0) {
1917            CglPreProcess * process = strategy_->process();
1918            assert (process);
1919            int n = solver_->getNumCols();
1920            const int * originalColumns = process->originalColumns();
1921            // columns should be in order ... but
1922            double * tempS = new double[n];
1923            for (int i = 0; i < n; i++) {
1924                int iColumn = originalColumns[i];
1925                tempS[i] = hotstartSolution_[iColumn];
1926            }
1927            delete [] hotstartSolution_;
1928            hotstartSolution_ = tempS;
1929            if (hotstartPriorities_) {
1930                int * tempP = new int [n];
1931                for (int i = 0; i < n; i++) {
1932                    int iColumn = originalColumns[i];
1933                    tempP[i] = hotstartPriorities_[iColumn];
1934                }
1935                delete [] hotstartPriorities_;
1936                hotstartPriorities_ = tempP;
1937            }
1938        }
1939        saveCompare = nodeCompare_;
1940        // depth first
1941        nodeCompare_ = new CbcCompareDepth();
1942    }
1943    if (!problemFeasibility_)
1944        problemFeasibility_ = new CbcFeasibilityBase();
1945# ifdef CBC_DEBUG
1946    std::string problemName ;
1947    solver_->getStrParam(OsiProbName, problemName) ;
1948    printf("Problem name - %s\n", problemName.c_str()) ;
1949    solver_->setHintParam(OsiDoReducePrint, false, OsiHintDo, 0) ;
1950# endif
1951    /*
1952      Assume we're done, and see if we're proven wrong.
1953    */
1954    status_ = 0 ;
1955    secondaryStatus_ = 0;
1956    phase_ = 0;
1957    /*
1958      Scan the variables, noting the integer variables. Create an
1959      CbcSimpleInteger object for each integer variable.
1960    */
1961    findIntegers(false) ;
1962    // Say not dynamic pseudo costs
1963    ownership_ &= ~0x40000000;
1964    // If dynamic pseudo costs then do
1965    if (numberBeforeTrust_)
1966        convertToDynamic();
1967    // Set up char array to say if integer (speed)
1968    delete [] integerInfo_;
1969    {
1970        int n = solver_->getNumCols();
1971        integerInfo_ = new char [n];
1972        for (int i = 0; i < n; i++) {
1973            if (solver_->isInteger(i))
1974                integerInfo_[i] = 1;
1975            else
1976                integerInfo_[i] = 0;
1977        }
1978    }
1979    if (preferredWay_) {
1980        // set all unset ones
1981        for (int iObject = 0 ; iObject < numberObjects_ ; iObject++) {
1982            CbcObject * obj =
1983                dynamic_cast <CbcObject *>(object_[iObject]) ;
1984            if (obj && !obj->preferredWay())
1985                obj->setPreferredWay(preferredWay_);
1986        }
1987    }
1988    /*
1989      Ensure that objects on the lists of OsiObjects, heuristics, and cut
1990      generators attached to this model all refer to this model.
1991    */
1992    synchronizeModel() ;
1993    if (!solverCharacteristics_) {
1994        OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
1995        if (solverCharacteristics) {
1996            solverCharacteristics_ = solverCharacteristics;
1997        } else {
1998            // replace in solver
1999            OsiBabSolver defaultC;
2000            solver_->setAuxiliaryInfo(&defaultC);
2001            solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
2002        }
2003    }
2004
2005    solverCharacteristics_->setSolver(solver_);
2006    // Set so we can tell we are in initial phase in resolve
2007    continuousObjective_ = -COIN_DBL_MAX ;
2008    /*
2009      Solve the relaxation.
2010
2011      Apparently there are circumstances where this will be non-trivial --- i.e.,
2012      we've done something since initialSolve that's trashed the solution to the
2013      continuous relaxation.
2014    */
2015    /* Tell solver we are in Branch and Cut
2016       Could use last parameter for subtle differences */
2017    solver_->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
2018#ifdef COIN_HAS_CLP
2019    {
2020        OsiClpSolverInterface * clpSolver
2021        = dynamic_cast<OsiClpSolverInterface *> (solver_);
2022        if (clpSolver) {
2023            ClpSimplex * clpSimplex = clpSolver->getModelPtr();
2024            if ((specialOptions_&32) == 0) {
2025                // take off names (unless going to be saving)
2026                if (numberAnalyzeIterations_>=0||(-numberAnalyzeIterations_&64)==0)
2027                  clpSimplex->dropNames();
2028            }
2029            // no crunch if mostly continuous
2030            if ((clpSolver->specialOptions()&(1 + 8)) != (1 + 8)) {
2031                int numberColumns = solver_->getNumCols();
2032                if (numberColumns > 1000 && numberIntegers_*4 < numberColumns)
2033                    clpSolver->setSpecialOptions(clpSolver->specialOptions()&(~1));
2034            }
2035            //#define NO_CRUNCH
2036#ifdef NO_CRUNCH
2037            printf("TEMP switching off crunch\n");
2038            int iOpt = clpSolver->specialOptions();
2039            iOpt &= ~1;
2040            iOpt |= 65536;
2041            clpSolver->setSpecialOptions(iOpt);
2042#endif
2043        }
2044    }
2045#endif
2046    bool feasible;
2047    numberSolves_ = 0 ;
2048    {
2049      // check
2050      int numberOdd = 0;
2051      for (int i = 0; i < numberObjects_; i++) {
2052        CbcSimpleInteger * obj =
2053          dynamic_cast <CbcSimpleInteger *>(object_[i]) ;
2054        if (!obj)
2055          numberOdd++;
2056      }
2057      if (numberOdd)
2058        moreSpecialOptions_ |= 1073741824;
2059    }
2060    // If NLP then we assume already solved outside branchAndbound
2061    if (!solverCharacteristics_->solverType() || solverCharacteristics_->solverType() == 4) {
2062        feasible = resolve(NULL, 0) != 0 ;
2063    } else {
2064        // pick up given status
2065        feasible = (solver_->isProvenOptimal() &&
2066                    !solver_->isDualObjectiveLimitReached()) ;
2067    }
2068    if (problemFeasibility_->feasible(this, 0) < 0) {
2069        feasible = false; // pretend infeasible
2070    }
2071    numberSavedSolutions_ = 0;
2072    int saveNumberStrong = numberStrong_;
2073    int saveNumberBeforeTrust = numberBeforeTrust_;
2074    /*
2075      If the linear relaxation of the root is infeasible, bail out now. Otherwise,
2076      continue with processing the root node.
2077    */
2078    if (!feasible) {
2079        status_ = 0 ;
2080        if (!solver_->isProvenDualInfeasible()) {
2081            handler_->message(CBC_INFEAS, messages_) << CoinMessageEol ;
2082            secondaryStatus_ = 1;
2083        } else {
2084            handler_->message(CBC_UNBOUNDED, messages_) << CoinMessageEol ;
2085            secondaryStatus_ = 7;
2086        }
2087        originalContinuousObjective_ = COIN_DBL_MAX;
2088        if (bestSolution_ && 
2089            ((specialOptions_&8388608)==0||(specialOptions_&2048)!=0)) {
2090          // best solution found by various heuristics - set solution
2091          char general[200];
2092          sprintf(general,"Solution of %g already found by heuristic",
2093                  bestObjective_);
2094          messageHandler()->message(CBC_GENERAL,
2095                                    messages())
2096            << general << CoinMessageEol ;
2097          setCutoff(1.0e50) ; // As best solution should be worse than cutoff
2098          // change cutoff as constraint if wanted
2099          if (cutoffRowNumber_>=0) {
2100            if (solver_->getNumRows()>cutoffRowNumber_)
2101              solver_->setRowUpper(cutoffRowNumber_,1.0e50);
2102          }
2103          // also in continuousSolver_
2104          if (continuousSolver_) {
2105            // Solvers know about direction
2106            double direction = solver_->getObjSense();
2107            continuousSolver_->setDblParam(OsiDualObjectiveLimit, 1.0e50*direction);
2108          } else {
2109            continuousSolver_ = solver_->clone();
2110          }
2111          phase_ = 5;
2112          double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
2113          if ((specialOptions_&4) == 0)
2114            bestObjective_ += 100.0 * increment + 1.0e-3; // only set if we are going to solve
2115          setBestSolution(CBC_END_SOLUTION, bestObjective_, bestSolution_, 1) ;
2116          continuousSolver_->resolve() ;
2117          if (!continuousSolver_->isProvenOptimal()) {
2118            continuousSolver_->messageHandler()->setLogLevel(2) ;
2119            continuousSolver_->initialSolve() ;
2120          }
2121          delete solver_ ;
2122          solver_ = continuousSolver_ ;
2123          setPointers(solver_);
2124          continuousSolver_ = NULL ;
2125        }
2126        solverCharacteristics_ = NULL;
2127        if (flipObjective)
2128          flipModel();
2129        return ;
2130    } else if (!numberObjects_) { 
2131        // nothing to do
2132        // Undo preprocessing performed during BaB.
2133        if (strategy_ && strategy_->preProcessState() > 0) {
2134          // undo preprocessing
2135          CglPreProcess * process = strategy_->process();
2136          assert (process);
2137          int n = originalSolver->getNumCols();
2138          if (bestSolution_) {
2139            delete [] bestSolution_;
2140            bestSolution_ = new double [n];
2141            process->postProcess(*solver_);
2142          }
2143          strategy_->deletePreProcess();
2144          // Solution now back in originalSolver
2145          delete solver_;
2146          solver_ = originalSolver;
2147          if (bestSolution_) {
2148            bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
2149            memcpy(bestSolution_, solver_->getColSolution(), n*sizeof(double));
2150          }
2151          // put back original objects if there were any
2152          if (originalObject) {
2153            int iColumn;
2154            assert (ownObjects_);
2155            for (iColumn = 0; iColumn < numberObjects_; iColumn++)
2156              delete object_[iColumn];
2157            delete [] object_;
2158            numberObjects_ = numberOriginalObjects;
2159            object_ = originalObject;
2160            delete [] integerVariable_;
2161            numberIntegers_ = 0;
2162            for (iColumn = 0; iColumn < n; iColumn++) {
2163              if (solver_->isInteger(iColumn))
2164                numberIntegers_++;
2165            }
2166            integerVariable_ = new int[numberIntegers_];
2167            numberIntegers_ = 0;
2168            for (iColumn = 0; iColumn < n; iColumn++) {
2169            if (solver_->isInteger(iColumn))
2170              integerVariable_[numberIntegers_++] = iColumn;
2171            }
2172          }
2173        }
2174        if (flipObjective)
2175          flipModel();
2176        solverCharacteristics_ = NULL;
2177        bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
2178        int numberColumns = solver_->getNumCols();
2179        delete [] bestSolution_;
2180        bestSolution_ = new double[numberColumns];
2181        CoinCopyN(solver_->getColSolution(), numberColumns, bestSolution_);
2182        return ;
2183    }
2184    /*
2185      See if we're using the Osi side of the branching hierarchy. If so, either
2186      convert existing CbcObjects to OsiObjects, or generate them fresh. In the
2187      first case, CbcModel owns the objects on the object_ list. In the second
2188      case, the solver holds the objects and object_ simply points to the
2189      solver's list.
2190
2191      080417 The conversion code here (the block protected by `if (obj)') cannot
2192      possibly be correct. On the Osi side, descent is OsiObject -> OsiObject2 ->
2193      all other Osi object classes. On the Cbc side, it's OsiObject -> CbcObject
2194      -> all other Cbc object classes. It's structurally impossible for any Osi
2195      object to descend from CbcObject. The only thing I can see is that this is
2196      really dead code, and object detection is now handled from the Osi side.
2197    */
2198    // Convert to Osi if wanted
2199    //OsiBranchingInformation * persistentInfo = NULL;
2200    if (branchingMethod_ && branchingMethod_->chooseMethod()) {
2201        //persistentInfo = new OsiBranchingInformation(solver_);
2202        if (numberOriginalObjects) {
2203            for (int iObject = 0 ; iObject < numberObjects_ ; iObject++) {
2204                CbcObject * obj =
2205                    dynamic_cast <CbcObject *>(object_[iObject]) ;
2206                if (obj) {
2207                    CbcSimpleInteger * obj2 =
2208                        dynamic_cast <CbcSimpleInteger *>(obj) ;
2209                    if (obj2) {
2210                        // back to Osi land
2211                        object_[iObject] = obj2->osiObject();
2212                        delete obj;
2213                    } else {
2214                        OsiSimpleInteger * obj3 =
2215                            dynamic_cast <OsiSimpleInteger *>(obj) ;
2216                        if (!obj3) {
2217                            OsiSOS * obj4 =
2218                                dynamic_cast <OsiSOS *>(obj) ;
2219                            if (!obj4) {
2220                                CbcSOS * obj5 =
2221                                    dynamic_cast <CbcSOS *>(obj) ;
2222                                if (obj5) {
2223                                    // back to Osi land
2224                                    object_[iObject] = obj5->osiObject(solver_);
2225                                } else {
2226                                    printf("Code up CbcObject type in Osi land\n");
2227                                    abort();
2228                                }
2229                            }
2230                        }
2231                    }
2232                }
2233            }
2234            // and add to solver
2235            //if (!solver_->numberObjects()) {
2236            solver_->addObjects(numberObjects_, object_);
2237            //} else {
2238            //if (solver_->numberObjects()!=numberOriginalObjects) {
2239            //printf("should have trapped that solver has objects before\n");
2240            //abort();
2241            //}
2242            //}
2243        } else {
2244            /*
2245              As of 080104, findIntegersAndSOS is misleading --- the default OSI
2246              implementation finds only integers.
2247            */
2248            // do from solver
2249            deleteObjects(false);
2250            solver_->findIntegersAndSOS(false);
2251            numberObjects_ = solver_->numberObjects();
2252            object_ = solver_->objects();
2253            ownObjects_ = false;
2254        }
2255        branchingMethod_->chooseMethod()->setSolver(solver_);
2256    }
2257    // take off heuristics if have to (some do not work with SOS, for example)
2258    // object should know what's safe.
2259    {
2260        int numberOdd = 0;
2261        int numberSOS = 0;
2262        for (int i = 0; i < numberObjects_; i++) {
2263            if (!object_[i]->canDoHeuristics())
2264                numberOdd++;
2265            CbcSOS * obj =
2266                dynamic_cast <CbcSOS *>(object_[i]) ;
2267            if (obj)
2268                numberSOS++;
2269        }
2270        if (numberOdd) {
2271          if (numberHeuristics_ && (specialOptions_&1024)==0 ) {
2272                int k = 0;
2273                for (int i = 0; i < numberHeuristics_; i++) {
2274                    if (!heuristic_[i]->canDealWithOdd())
2275                        delete heuristic_[i];
2276                    else
2277                        heuristic_[k++] = heuristic_[i];
2278                }
2279                if (!k) {
2280                    delete [] heuristic_;
2281                    heuristic_ = NULL;
2282                }
2283                numberHeuristics_ = k;
2284                handler_->message(CBC_HEURISTICS_OFF, messages_) << numberOdd << CoinMessageEol ;
2285            }
2286            // If odd switch off AddIntegers
2287            specialOptions_ &= ~65536;
2288            // switch off fast nodes for now
2289            fastNodeDepth_ = -1;
2290            moreSpecialOptions_ &= ~33554432; // no diving
2291        } else if (numberSOS) {
2292            specialOptions_ |= 128; // say can do SOS in dynamic mode
2293            // switch off fast nodes for now
2294            fastNodeDepth_ = -1;
2295            moreSpecialOptions_ &= ~33554432; // no diving
2296        }
2297        if (numberThreads_ > 0) {
2298            // switch off fast nodes for now
2299            //fastNodeDepth_ = -1;
2300        }
2301    }
2302    // Save objective (just so user can access it)
2303    originalContinuousObjective_ = solver_->getObjValue()* solver_->getObjSense();
2304    bestPossibleObjective_ = originalContinuousObjective_;
2305    sumChangeObjective1_ = 0.0;
2306    sumChangeObjective2_ = 0.0;
2307    /*
2308      OsiRowCutDebugger knows an optimal answer for a subset of MIP problems.
2309      Assuming it recognises the problem, when called upon it will check a cut to
2310      see if it cuts off the optimal answer.
2311    */
2312    // If debugger exists set specialOptions_ bit
2313    if (solver_->getRowCutDebuggerAlways()) {
2314        specialOptions_ |= 1;
2315    }
2316
2317# ifdef CBC_DEBUG
2318    if ((specialOptions_&1) == 0)
2319        solver_->activateRowCutDebugger(problemName.c_str()) ;
2320    if (solver_->getRowCutDebuggerAlways())
2321        specialOptions_ |= 1;
2322# endif
2323
2324    /*
2325      Begin setup to process a feasible root node.
2326    */
2327    bestObjective_ = CoinMin(bestObjective_, 1.0e50) ;
2328    if (!bestSolution_) {
2329        numberSolutions_ = 0 ;
2330        numberHeuristicSolutions_ = 0 ;
2331    }
2332    stateOfSearch_ = 0;
2333    // Everything is minimization
2334    {
2335        // needed to sync cutoffs
2336        double value ;
2337        solver_->getDblParam(OsiDualObjectiveLimit, value) ;
2338        dblParam_[CbcCurrentCutoff] = value * solver_->getObjSense();
2339    }
2340    double cutoff = getCutoff() ;
2341    double direction = solver_->getObjSense() ;
2342    dblParam_[CbcOptimizationDirection] = direction;
2343    if (cutoff < 1.0e20 && direction < 0.0)
2344        messageHandler()->message(CBC_CUTOFF_WARNING1,
2345                                  messages())
2346        << cutoff << -cutoff << CoinMessageEol ;
2347    if (cutoff > bestObjective_)
2348        cutoff = bestObjective_ ;
2349    setCutoff(cutoff) ;
2350    /*
2351      We probably already have a current solution, but just in case ...
2352    */
2353    int numberColumns = getNumCols() ;
2354    if (!currentSolution_)
2355        currentSolution_ = new double[numberColumns] ;
2356    testSolution_ = currentSolution_;
2357    /*
2358      Create a copy of the solver, thus capturing the original (root node)
2359      constraint system (aka the continuous system).
2360    */
2361    delete continuousSolver_;
2362    continuousSolver_ = solver_->clone() ;
2363#ifdef COIN_HAS_NTY
2364    // maybe allow on fix and restart later
2365    if ((moreSpecialOptions2_&(128|256))!=0&&!parentModel_) {
2366      symmetryInfo_ = new CbcSymmetry();
2367      symmetryInfo_->setupSymmetry(*continuousSolver_);
2368      int numberGenerators = symmetryInfo_->statsOrbits(this,0);
2369      if (!symmetryInfo_->numberUsefulOrbits()&&(moreSpecialOptions2_&(128|256))!=(128|256)) {
2370        delete symmetryInfo_;
2371        symmetryInfo_=NULL;
2372        moreSpecialOptions2_ &= ~(128|256);
2373      }
2374      if ((moreSpecialOptions2_&(128|256))==(128|256)) {
2375        //moreSpecialOptions2_ &= ~256;
2376      }
2377    }
2378#endif
2379
2380    // add cutoff as constraint if wanted
2381    if (cutoffRowNumber_==-2) {
2382      if (!parentModel_) {
2383        int numberColumns=solver_->getNumCols();
2384        double * obj = CoinCopyOfArray(solver_->getObjCoefficients(),numberColumns);
2385        int * indices = new int [numberColumns];
2386        int n=0;
2387        for (int i=0;i<numberColumns;i++) {
2388          if (obj[i]) {
2389            indices[n]=i;
2390            obj[n++]=obj[i];
2391          }
2392        }
2393        if (n) {
2394          double cutoff=getCutoff();
2395          // relax a little bit
2396          cutoff += 1.0e-4;
2397          double offset;
2398          solver_->getDblParam(OsiObjOffset, offset);
2399          cutoffRowNumber_ = solver_->getNumRows();
2400          solver_->addRow(n,indices,obj,-COIN_DBL_MAX,CoinMin(cutoff,1.0e25)+offset);
2401        } else {
2402          // no objective!
2403          cutoffRowNumber_ = -1;
2404        }
2405        delete [] indices;
2406        delete [] obj;
2407      } else {
2408        // switch off
2409        cutoffRowNumber_ = -1;
2410      }
2411    }
2412    numberRowsAtContinuous_ = getNumRows() ;
2413    solver_->saveBaseModel();
2414    /*
2415      Check the objective to see if we can deduce a nontrivial increment. If
2416      it's better than the current value for CbcCutoffIncrement, it'll be
2417      installed.
2418    */
2419    if (solverCharacteristics_->reducedCostsAccurate())
2420        analyzeObjective() ;
2421    {
2422        // may be able to change cutoff now
2423        double cutoff = getCutoff();
2424        double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
2425        if (cutoff > bestObjective_ - increment) {
2426            cutoff = bestObjective_ - increment ;
2427            setCutoff(cutoff) ;
2428        }
2429    }
2430#ifdef COIN_HAS_CLP
2431    // Possible save of pivot method
2432    ClpDualRowPivot * savePivotMethod = NULL;
2433    {
2434        // pass tolerance and increment to solver
2435        OsiClpSolverInterface * clpSolver
2436        = dynamic_cast<OsiClpSolverInterface *> (solver_);
2437        if (clpSolver)
2438            clpSolver->setStuff(getIntegerTolerance(), getCutoffIncrement());
2439#ifdef CLP_RESOLVE
2440        if ((moreSpecialOptions_&1048576)!=0&&!parentModel_&&clpSolver) {
2441          resolveClp(clpSolver,0);
2442        }
2443#endif
2444    }
2445#endif
2446    /*
2447      Set up for cut generation. addedCuts_ holds the cuts which are relevant for
2448      the active subproblem. whichGenerator will be used to record the generator
2449      that produced a given cut.
2450    */
2451    maximumWhich_ = 1000 ;
2452    delete [] whichGenerator_;
2453    whichGenerator_ = new int[maximumWhich_] ;
2454    memset(whichGenerator_, 0, maximumWhich_*sizeof(int));
2455    maximumNumberCuts_ = 0 ;
2456    currentNumberCuts_ = 0 ;
2457    delete [] addedCuts_ ;
2458    addedCuts_ = NULL ;
2459    OsiObject ** saveObjects = NULL;
2460    maximumRows_ = numberRowsAtContinuous_;
2461    currentDepth_ = 0;
2462    workingBasis_.resize(maximumRows_, numberColumns);
2463    /*
2464      Set up an empty heap and associated data structures to hold the live set
2465      (problems which require further exploration).
2466    */
2467    CbcCompareDefault * compareActual
2468    = dynamic_cast<CbcCompareDefault *> (nodeCompare_);
2469    if (compareActual) {
2470        compareActual->setBestPossible(direction*solver_->getObjValue());
2471        compareActual->setCutoff(getCutoff());
2472#ifdef JJF_ZERO
2473        if (false && !numberThreads_ && !parentModel_) {
2474            printf("CbcTreeArray ? threads ? parentArray\n");
2475            // Setup new style tree
2476            delete tree_;
2477            tree_ = new CbcTreeArray();
2478        }
2479#endif
2480    }
2481    tree_->setComparison(*nodeCompare_) ;
2482    /*
2483      Used to record the path from a node to the root of the search tree, so that
2484      we can then traverse from the root to the node when restoring a subproblem.
2485    */
2486    maximumDepth_ = 10 ;
2487    delete [] walkback_ ;
2488    walkback_ = new CbcNodeInfo * [maximumDepth_] ;
2489    lastDepth_ = 0;
2490    delete [] lastNodeInfo_ ;
2491    lastNodeInfo_ = new CbcNodeInfo * [maximumDepth_] ;
2492    delete [] lastNumberCuts_ ;
2493    lastNumberCuts_ = new int [maximumDepth_] ;
2494    maximumCuts_ = 100;
2495    lastNumberCuts2_ = 0;
2496    delete [] lastCut_;
2497    lastCut_ = new const OsiRowCut * [maximumCuts_];
2498    /*
2499      Used to generate bound edits for CbcPartialNodeInfo.
2500    */
2501    double * lowerBefore = new double [numberColumns] ;
2502    double * upperBefore = new double [numberColumns] ;
2503    /*
2504    Set up to run heuristics and generate cuts at the root node. The heavy
2505    lifting is hidden inside the calls to doHeuristicsAtRoot and solveWithCuts.
2506
2507    To start, tell cut generators they can be a bit more aggressive at the
2508    root node.
2509
2510    QUESTION: phase_ = 0 is documented as `initial solve', phase = 1 as `solve
2511        with cuts at root'. Is phase_ = 1 the correct indication when
2512        doHeurisiticsAtRoot is called to run heuristics outside of the main
2513        cut / heurisitc / reoptimise loop in solveWithCuts?
2514
2515      Generate cuts at the root node and reoptimise. solveWithCuts does the heavy
2516      lifting. It will iterate a generate/reoptimise loop (including reduced cost
2517      fixing) until no cuts are generated, the change in objective falls off,  or
2518      the limit on the number of rounds of cut generation is exceeded.
2519
2520      At the end of all this, any cuts will be recorded in cuts and also
2521      installed in the solver's constraint system. We'll have reoptimised, and
2522      removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
2523      adjusted accordingly).
2524
2525      Tell cut generators they can be a bit more aggressive at root node
2526
2527      TODO: Why don't we make a copy of the solution after solveWithCuts?
2528      TODO: If numberUnsatisfied == 0, don't we have a solution?
2529    */
2530    phase_ = 1;
2531    int iCutGenerator;
2532    for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
2533        // If parallel switch off global cuts
2534        if (numberThreads_) {
2535            generator_[iCutGenerator]->setGlobalCuts(false);
2536            generator_[iCutGenerator]->setGlobalCutsAtRoot(false);
2537        }
2538        CglCutGenerator * generator = generator_[iCutGenerator]->generator();
2539        generator->setAggressiveness(generator->getAggressiveness() + 100);
2540        if (!generator->canDoGlobalCuts())
2541          generator->setGlobalCuts(false);
2542    }
2543    OsiCuts cuts ;
2544    int anyAction = -1 ;
2545    numberOldActiveCuts_ = 0 ;
2546    numberNewCuts_ = 0 ;
2547    // Array to mark solution
2548    delete [] usedInSolution_;
2549    usedInSolution_ = new int[numberColumns];
2550    CoinZeroN(usedInSolution_, numberColumns);
2551    /*
2552      For printing totals and for CbcNode (numberNodes_)
2553    */
2554    numberIterations_ = 0 ;
2555    numberNodes_ = 0 ;
2556    numberNodes2_ = 0 ;
2557    maximumStatistics_ = 0;
2558    maximumDepthActual_ = 0;
2559    numberDJFixed_ = 0.0;
2560    if (!parentModel_) {
2561        if ((specialOptions_&262144) != 0) {
2562            // create empty stored cuts
2563            //storedRowCuts_ = new CglStored(solver_->getNumCols());
2564        } else if ((specialOptions_&524288) != 0 && storedRowCuts_) {
2565            // tighten and set best solution
2566            // A) tight bounds on integer variables
2567            /*
2568                storedRowCuts_ are coming in from outside, probably for nonlinear.
2569              John was unsure about origin.
2570            */
2571            const double * lower = solver_->getColLower();
2572            const double * upper = solver_->getColUpper();
2573            const double * tightLower = storedRowCuts_->tightLower();
2574            const double * tightUpper = storedRowCuts_->tightUpper();
2575            int nTightened = 0;
2576            for (int i = 0; i < numberIntegers_; i++) {
2577                int iColumn = integerVariable_[i];
2578                if (tightLower[iColumn] > lower[iColumn]) {
2579                    nTightened++;
2580                    solver_->setColLower(iColumn, tightLower[iColumn]);
2581                }
2582                if (tightUpper[iColumn] < upper[iColumn]) {
2583                    nTightened++;
2584                    solver_->setColUpper(iColumn, tightUpper[iColumn]);
2585                }
2586            }
2587            if (nTightened)
2588              COIN_DETAIL_PRINT(printf("%d tightened by alternate cuts\n", nTightened));
2589            if (storedRowCuts_->bestObjective() < bestObjective_) {
2590                // B) best solution
2591                double objValue = storedRowCuts_->bestObjective();
2592                setBestSolution(CBC_SOLUTION, objValue,
2593                                storedRowCuts_->bestSolution()) ;
2594                // Do heuristics
2595                // Allow RINS
2596                for (int i = 0; i < numberHeuristics_; i++) {
2597                    CbcHeuristicRINS * rins
2598                    = dynamic_cast<CbcHeuristicRINS *> (heuristic_[i]);
2599                    if (rins) {
2600                        rins->setLastNode(-100);
2601                    }
2602                }
2603            }
2604        }
2605    }
2606#ifdef SWITCH_VARIABLES
2607    // see if any switching variables
2608    if (numberIntegers_<solver_->getNumCols())
2609      findSwitching();
2610#endif
2611    /*
2612      Run heuristics at the root. This is the only opportunity to run FPump; it
2613      will be removed from the heuristics list by doHeuristicsAtRoot.
2614    */
2615    // See if multiple runs wanted
2616    CbcModel ** rootModels=NULL;
2617    if (!parentModel_&&multipleRootTries_%100) {
2618      double rootTimeCpu=CoinCpuTime();
2619      double startTimeRoot=CoinGetTimeOfDay();
2620      int numberRootThreads=1;
2621      /* undocumented fine tuning
2622         aabbcc where cc is number of tries
2623         bb if nonzero is number of threads
2624         aa if nonzero just do heuristics
2625      */
2626      int numberModels = multipleRootTries_%100;
2627#ifdef CBC_THREAD
2628      numberRootThreads = (multipleRootTries_/100)%100;
2629      if (!numberRootThreads)
2630        numberRootThreads=numberModels;
2631#endif
2632      int otherOptions = (multipleRootTries_/10000)%100;
2633      rootModels = new CbcModel * [numberModels];
2634      unsigned int newSeed = randomSeed_;
2635      if (newSeed==0) {
2636        double time = fabs(CoinGetTimeOfDay());
2637        while (time>=COIN_INT_MAX)
2638          time *= 0.5;
2639        newSeed = static_cast<unsigned int>(time);
2640      } else if (newSeed<0) {
2641        newSeed = 123456789;
2642#ifdef COIN_HAS_CLP
2643        OsiClpSolverInterface * clpSolver
2644          = dynamic_cast<OsiClpSolverInterface *> (solver_);
2645        if (clpSolver) {
2646          newSeed += clpSolver->getModelPtr()->randomNumberGenerator()->getSeed();
2647        }
2648#endif
2649      }
2650      CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis *> (solver_->getEmptyWarmStart());
2651      for (int i=0;i<numberModels;i++) { 
2652        rootModels[i]=new CbcModel(*this);
2653        rootModels[i]->setNumberThreads(0);
2654        rootModels[i]->setMaximumNodes(otherOptions ? -1 : 0);
2655        rootModels[i]->setRandomSeed(newSeed+10000000*i);
2656        rootModels[i]->randomNumberGenerator()->setSeed(newSeed+50000000*i);
2657        rootModels[i]->setMultipleRootTries(0);
2658        // use seed
2659        rootModels[i]->setSpecialOptions(specialOptions_ |(4194304|8388608));
2660        rootModels[i]->setMoreSpecialOptions(moreSpecialOptions_ & 
2661                                             (~(134217728|4194304)));
2662        rootModels[i]->setMoreSpecialOptions2(moreSpecialOptions2_ & 
2663                                              (~(128|256)));
2664        rootModels[i]->solver_->setWarmStart(basis);
2665#ifdef COIN_HAS_CLP
2666        OsiClpSolverInterface * clpSolver
2667          = dynamic_cast<OsiClpSolverInterface *> (rootModels[i]->solver_);
2668#define NEW_RANDOM_BASIS
2669#ifdef NEW_RANDOM_BASIS
2670        if (i==0)
2671          continue;
2672#endif
2673        if (clpSolver) {
2674          ClpSimplex * simplex = clpSolver->getModelPtr();
2675          if (defaultHandler_)
2676            simplex->setDefaultMessageHandler();
2677          simplex->setRandomSeed(newSeed+20000000*i);
2678          simplex->allSlackBasis();
2679          int logLevel=simplex->logLevel();
2680          if (logLevel==1)
2681            simplex->setLogLevel(0);
2682          if (i!=0) {
2683#ifdef NEW_RANDOM_BASIS
2684            int numberRows = simplex->numberRows();
2685            int throwOut=20;//2+numberRows/100;
2686            for (int iThrow=0;iThrow<throwOut;iThrow++) {
2687              double random = simplex->randomNumberGenerator()->randomDouble();
2688              int iStart=static_cast<int>(random*numberRows);
2689              for (int j=iStart;j<numberRows;j++) {
2690                if (simplex->getRowStatus(j)!=ClpSimplex::basic) {
2691                  simplex->setRowStatus(j,ClpSimplex::basic);
2692                  break;
2693                }
2694              }
2695            }
2696            clpSolver->setWarmStart(NULL);
2697#else
2698            double random = simplex->randomNumberGenerator()->randomDouble();
2699            int bias = static_cast<int>(random*(numberIterations/4));
2700            simplex->setMaximumIterations(numberIterations/2+bias);
2701            simplex->primal();
2702            simplex->setMaximumIterations(COIN_INT_MAX);
2703            simplex->dual();
2704#endif
2705          } else {
2706#ifndef NEW_RANDOM_BASIS
2707            simplex->primal();
2708#endif
2709#endif
2710          }
2711#ifdef NEW_RANDOM_BASIS
2712          simplex->setLogLevel(logLevel);
2713          clpSolver->setWarmStart(NULL);
2714#endif
2715        }
2716        for (int j=0;j<numberHeuristics_;j++)
2717          rootModels[i]->heuristic_[j]->setSeed(rootModels[i]->heuristic_[j]->getSeed()+100000000*i);
2718        for (int j=0;j<numberCutGenerators_;j++)
2719          rootModels[i]->generator_[j]->generator()->refreshSolver(rootModels[i]->solver_);
2720      }
2721      delete basis;
2722#ifdef CBC_THREAD
2723      if (numberRootThreads==1) {
2724#endif
2725        for (int iModel=0;iModel<numberModels;iModel++) {
2726          doRootCbcThread(rootModels[iModel]);
2727          // see if solved at root node
2728          if (rootModels[iModel]->getMaximumNodes()) {
2729            feasible=false;
2730            break;
2731          }
2732        }
2733#ifdef CBC_THREAD
2734      } else {
2735        Coin_pthread_t * threadId = new Coin_pthread_t [numberRootThreads];
2736        for (int kModel=0;kModel<numberModels;kModel+=numberRootThreads) {
2737          bool finished=false;
2738          for (int iModel=kModel;iModel<CoinMin(numberModels,kModel+numberRootThreads);iModel++) {
2739            pthread_create(&(threadId[iModel-kModel].thr), NULL, 
2740                           doRootCbcThread,
2741                           rootModels[iModel]);
2742          }
2743          // wait
2744          for (int iModel=kModel;iModel<CoinMin(numberModels,kModel+numberRootThreads);iModel++) {
2745            pthread_join(threadId[iModel-kModel].thr, NULL);
2746          }
2747          // see if solved at root node
2748          for (int iModel=kModel;iModel<CoinMin(numberModels,kModel+numberRootThreads);iModel++) {
2749            if (rootModels[iModel]->getMaximumNodes())
2750              finished=true;
2751          }
2752          if (finished) {
2753            feasible=false;
2754            break;
2755          }
2756        }
2757        delete [] threadId;
2758      }
2759#endif
2760      // sort solutions
2761      int * which = new int [numberModels];
2762      double * value = new double [numberModels];
2763      int numberSolutions=0;
2764      for (int iModel=0;iModel<numberModels;iModel++) {
2765        if (rootModels[iModel]->bestSolution()) {
2766          which[numberSolutions]=iModel;
2767          value[numberSolutions++]=
2768            -rootModels[iModel]->getMinimizationObjValue();
2769        }
2770      }
2771      char general[100];
2772      rootTimeCpu=CoinCpuTime()-rootTimeCpu;
2773      if (numberRootThreads==1)
2774        sprintf(general,"Multiple root solvers took a total of %.2f seconds\n",
2775                rootTimeCpu);
2776      else
2777        sprintf(general,"Multiple root solvers took a total of %.2f seconds (%.2f elapsed)\n",
2778                rootTimeCpu,CoinGetTimeOfDay()-startTimeRoot);
2779      messageHandler()->message(CBC_GENERAL,
2780                                messages())
2781        << general << CoinMessageEol ;
2782      CoinSort_2(value,value+numberSolutions,which);
2783      // to get name
2784      CbcHeuristicRINS dummyHeuristic;
2785      dummyHeuristic.setHeuristicName("Multiple root solvers");
2786      lastHeuristic_=&dummyHeuristic;
2787      for (int i=0;i<numberSolutions;i++) {
2788        double objValue = - value[i];
2789        if (objValue<getCutoff()) {
2790          int iModel=which[i];
2791          setBestSolution(CBC_ROUNDING,objValue,
2792                          rootModels[iModel]->bestSolution());
2793        }
2794      }
2795      lastHeuristic_=NULL;
2796      delete [] which;
2797      delete [] value;
2798    }
2799    // Do heuristics
2800    if (numberObjects_&&!rootModels)
2801        doHeuristicsAtRoot();
2802    if (solverCharacteristics_->solutionAddsCuts()) {
2803      // With some heuristics solver needs a resolve here
2804      solver_->resolve();
2805      if(!isProvenOptimal()){
2806        solver_->initialSolve();
2807      }
2808    }
2809    /*
2810      Grepping through the code, it would appear that this is a command line
2811      debugging hook.  There's no obvious place in the code where this is set to
2812      a negative value.
2813
2814      User hook, says John.
2815    */
2816    if ( intParam_[CbcMaxNumNode] < 0
2817      ||numberSolutions_>=getMaximumSolutions())
2818      eventHappened_ = true; // stop as fast as possible
2819    stoppedOnGap_ = false ;
2820    // See if can stop on gap
2821    bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
2822    if(canStopOnGap()) {
2823        if (bestPossibleObjective_ < getCutoff())
2824            stoppedOnGap_ = true ;
2825        feasible = false;
2826        //eventHappened_=true; // stop as fast as possible
2827    }
2828    /*
2829      Set up for statistics collection, if requested. Standard values are
2830      documented in CbcModel.hpp. The magic number 100 will trigger a dump of
2831      CbcSimpleIntegerDynamicPseudoCost objects (no others). Looks like another
2832      command line debugging hook.
2833    */
2834    statistics_ = NULL;
2835    // Do on switch
2836    if (doStatistics > 0 && doStatistics <= 100) {
2837        maximumStatistics_ = 10000;
2838        statistics_ = new CbcStatistics * [maximumStatistics_];
2839        memset(statistics_, 0, maximumStatistics_*sizeof(CbcStatistics *));
2840    }
2841    // See if we can add integers
2842    if (noObjects && numberIntegers_ < solver_->getNumCols() && (specialOptions_&65536) != 0 && !parentModel_)
2843        AddIntegers();
2844    /*
2845      Do an initial round of cut generation for the root node. Depending on the
2846      type of underlying solver, we may want to do this even if the initial query
2847      to the objects indicates they're satisfied.
2848
2849      solveWithCuts does the heavy lifting. It will iterate a generate/reoptimise
2850      loop (including reduced cost fixing) until no cuts are generated, the
2851      change in objective falls off,  or the limit on the number of rounds of cut
2852      generation is exceeded.
2853
2854      At the end of all this, any cuts will be recorded in cuts and also
2855      installed in the solver's constraint system. We'll have reoptimised, and
2856      removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
2857      adjusted accordingly).
2858    */
2859    int iObject ;
2860    int numberUnsatisfied = 0 ;
2861    delete [] currentSolution_;
2862    currentSolution_ = new double [numberColumns];
2863    testSolution_ = currentSolution_;
2864    memcpy(currentSolution_, solver_->getColSolution(),
2865           numberColumns*sizeof(double)) ;
2866    // point to useful information
2867    OsiBranchingInformation usefulInfo = usefulInformation();
2868
2869    for (iObject = 0 ; iObject < numberObjects_ ; iObject++) {
2870        double infeasibility =
2871            object_[iObject]->checkInfeasibility(&usefulInfo) ;
2872        if (infeasibility ) numberUnsatisfied++ ;
2873    }
2874    // replace solverType
2875    double * tightBounds = NULL;
2876    if (solverCharacteristics_->tryCuts())  {
2877
2878        if (numberUnsatisfied)   {
2879            // User event
2880            if (!eventHappened_ && feasible) {
2881                if (rootModels) {
2882                  // for fixings
2883                  int numberColumns=solver_->getNumCols();
2884                  tightBounds = new double [2*numberColumns];
2885                  {
2886                    const double * lower = solver_->getColLower();
2887                    const double * upper = solver_->getColUpper();
2888                    for (int i=0;i<numberColumns;i++) {
2889                      tightBounds[2*i+0]=lower[i];
2890                      tightBounds[2*i+1]=upper[i];
2891                    }
2892                  }
2893                  int numberModels = multipleRootTries_%100;
2894                  const OsiSolverInterface ** solvers = new 
2895                    const OsiSolverInterface * [numberModels];
2896                  int numberRows=continuousSolver_->getNumRows();
2897                  int maxCuts=0;
2898                  for (int i=0;i<numberModels;i++) {
2899                    solvers[i]=rootModels[i]->solver();
2900                    const double * lower = solvers[i]->getColLower();
2901                    const double * upper = solvers[i]->getColUpper();
2902                    for (int j=0;j<numberColumns;j++) {
2903                      tightBounds[2*j+0]=CoinMax(lower[j],tightBounds[2*j+0]);
2904                      tightBounds[2*j+1]=CoinMin(upper[j],tightBounds[2*j+1]);
2905                    }
2906                    int numberRows2=solvers[i]->getNumRows();
2907                    assert (numberRows2>=numberRows);
2908                    maxCuts += numberRows2-numberRows;
2909                    // accumulate statistics
2910                    for (int j=0;j<numberCutGenerators_;j++) {
2911                      generator_[j]->addStatistics(rootModels[i]->cutGenerator(j));
2912                    }
2913                  }
2914                  for (int j=0;j<numberCutGenerators_;j++) {
2915                    generator_[j]->scaleBackStatistics(numberModels);
2916                  }
2917                  //CbcRowCuts rowCut(maxCuts);
2918                  const OsiRowCutDebugger *debugger = NULL;
2919                  if ((specialOptions_&1) != 0) 
2920                    debugger = solver_->getRowCutDebugger() ;
2921                  for (int iModel=0;iModel<numberModels;iModel++) {
2922                    int numberRows2=solvers[iModel]->getNumRows();
2923                    const CoinPackedMatrix * rowCopy = solvers[iModel]->getMatrixByRow();
2924                    const int * rowLength = rowCopy->getVectorLengths(); 
2925                    const double * elements = rowCopy->getElements();
2926                    const int * column = rowCopy->getIndices();
2927                    const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
2928                    const double * rowLower = solvers[iModel]->getRowLower();
2929                    const double * rowUpper = solvers[iModel]->getRowUpper();
2930                    for (int iRow=numberRows;iRow<numberRows2;iRow++) {
2931                      OsiRowCut rc;
2932                      rc.setLb(rowLower[iRow]);
2933                      rc.setUb(rowUpper[iRow]);
2934                      CoinBigIndex start = rowStart[iRow];
2935                      rc.setRow(rowLength[iRow],column+start,elements+start,false);
2936                      if (debugger)
2937                        CoinAssert (!debugger->invalidCut(rc));
2938                      globalCuts_.addCutIfNotDuplicate(rc);
2939                    }
2940                    //int cutsAdded=globalCuts_.numberCuts()-numberCuts;
2941                    //numberCuts += cutsAdded;
2942                    //printf("Model %d gave %d cuts (out of %d possible)\n",
2943                    //     iModel,cutsAdded,numberRows2-numberRows);
2944                  }
2945                  // normally replace global cuts
2946                  //if (!globalCuts_.())
2947                  //globalCuts_=rowCutrowCut.addCuts(globalCuts_);
2948                  //rowCut.addCuts(globalCuts_);
2949                  int nTightened=0;
2950                  assert(feasible);
2951                  {
2952                    double tolerance=1.0e-5;
2953                    const double * lower = solver_->getColLower();
2954                    const double * upper = solver_->getColUpper();
2955                    for (int i=0;i<numberColumns;i++) {
2956                      if (tightBounds[2*i+0]>tightBounds[2*i+1]+1.0e-9) {
2957                        feasible=false;
2958                        char general[200];
2959                        sprintf(general,"Solvers give infeasible bounds on %d %g,%g was %g,%g - search finished\n",
2960                               i,tightBounds[2*i+0],tightBounds[2*i+1],lower[i],upper[i]);
2961                        messageHandler()->message(CBC_GENERAL,messages())
2962                          << general << CoinMessageEol ;
2963                        break;
2964                      }
2965                      double oldLower=lower[i];
2966                      double oldUpper=upper[i];
2967                      if (tightBounds[2*i+0]>oldLower+tolerance) {
2968                        nTightened++;
2969                        solver_->setColLower(i,tightBounds[2*i+0]);
2970                      }
2971                      if (tightBounds[2*i+1]<oldUpper-tolerance) {
2972                        nTightened++;
2973                        solver_->setColUpper(i,tightBounds[2*i+1]);
2974                      }
2975                    }
2976                  }
2977                  delete [] tightBounds;
2978                  tightBounds=NULL;
2979                  char printBuffer[200];
2980                  sprintf(printBuffer,"%d solvers added %d different cuts out of pool of %d",
2981                          numberModels,globalCuts_.sizeRowCuts(),maxCuts);
2982                  messageHandler()->message(CBC_GENERAL, messages())
2983                    << printBuffer << CoinMessageEol ;
2984                  if (nTightened) {
2985                    sprintf(printBuffer,"%d bounds were tightened",
2986                          nTightened);
2987                    messageHandler()->message(CBC_GENERAL, messages())
2988                      << printBuffer << CoinMessageEol ;
2989                  }
2990                  delete [] solvers;
2991                }
2992                if (!parentModel_&&(moreSpecialOptions_&67108864) != 0) {
2993                  // load cuts from file
2994                  FILE * fp = fopen("global.cuts","rb");
2995                  if (fp) {
2996                    size_t nRead;
2997                    int numberColumns=solver_->getNumCols();
2998                    int numCols;
2999                    nRead = fread(&numCols, sizeof(int), 1, fp);
3000                    if (nRead != 1)
3001                      throw("Error in fread");
3002                    if (numberColumns!=numCols) {
3003                      printf("Mismatch on columns %d %d\n",numberColumns,numCols);
3004                      fclose(fp);
3005                    } else {
3006                      // If rootModel just do some
3007                      double threshold=-1.0;
3008                      if (!multipleRootTries_)
3009                        threshold=0.5;
3010                      int initialCuts=0;
3011                      int initialGlobal = globalCuts_.sizeRowCuts();
3012                      double * elements = new double [numberColumns+2];
3013                      int * indices = new int [numberColumns];
3014                      int numberEntries=1;
3015                      while (numberEntries>0) {
3016                        nRead = fread(&numberEntries, sizeof(int), 1, fp);
3017                        if (nRead != 1)
3018                          throw("Error in fread");
3019                        double randomNumber=randomNumberGenerator_.randomDouble();
3020                        if (numberEntries>0) {
3021                          initialCuts++;
3022                          nRead = fread(elements, sizeof(double), numberEntries+2, fp);
3023                          if (nRead != static_cast<size_t>(numberEntries+2))
3024                            throw("Error in fread");
3025                          nRead = fread(indices, sizeof(int), numberEntries, fp);
3026                          if (nRead != static_cast<size_t>(numberEntries))
3027                            throw("Error in fread");
3028                          if (randomNumber>threshold) {
3029                            OsiRowCut rc;
3030                            rc.setLb(elements[numberEntries]);
3031                            rc.setUb(elements[numberEntries+1]);
3032                            rc.setRow(numberEntries,indices,elements,
3033                                      false);
3034                            rc.setGloballyValidAsInteger(2);
3035                            globalCuts_.addCutIfNotDuplicate(rc) ;
3036                          } 
3037                        }
3038                      }
3039                      fclose(fp);
3040                      // fixes
3041                      int nTightened=0;
3042                      fp = fopen("global.fix","rb");
3043                      if (fp) {
3044                        nRead = fread(indices, sizeof(int), 2, fp);
3045                        if (nRead != 2)
3046                          throw("Error in fread");
3047                        if (numberColumns!=indices[0]) {
3048                          printf("Mismatch on columns %d %d\n",numberColumns,
3049                                 indices[0]);
3050                        } else {
3051                          indices[0]=1;
3052                          while (indices[0]>=0) {
3053                            nRead = fread(indices, sizeof(int), 2, fp);
3054                            if (nRead != 2)
3055                              throw("Error in fread");
3056                            int iColumn=indices[0];
3057                            if (iColumn>=0) {
3058                              nTightened++;
3059                              nRead = fread(elements, sizeof(double), 4, fp);
3060                              if (nRead != 4)
3061                                throw("Error in fread");
3062                              solver_->setColLower(iColumn,elements[0]);
3063                              solver_->setColUpper(iColumn,elements[1]);
3064                            } 
3065                          }
3066                        }
3067                      }
3068                      if (fp)
3069                        fclose(fp);
3070                      char printBuffer[200];
3071                      sprintf(printBuffer,"%d cuts read in of which %d were unique, %d bounds tightened",
3072                             initialCuts,
3073                             globalCuts_.sizeRowCuts()-initialGlobal,nTightened); 
3074                      messageHandler()->message(CBC_GENERAL, messages())
3075                        << printBuffer << CoinMessageEol ;
3076                      delete [] elements;
3077                      delete [] indices;
3078                    }
3079                  }
3080                }
3081                if (feasible)
3082                  feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
3083                                         NULL);
3084                if (multipleRootTries_&&
3085                    (moreSpecialOptions_&134217728)!=0) {
3086                  FILE * fp=NULL;
3087                  size_t nRead;
3088                  int numberColumns=solver_->getNumCols();
3089                  int initialCuts=0;
3090                  if ((moreSpecialOptions_&134217728)!=0) {
3091                    // append so go down to end
3092                    fp = fopen("global.cuts","r+b");
3093                    if (fp) {
3094                      int numCols;
3095                      nRead = fread(&numCols, sizeof(int), 1, fp);
3096                      if (nRead != 1)
3097                        throw("Error in fread");
3098                      if (numberColumns!=numCols) {
3099                        printf("Mismatch on columns %d %d\n",numberColumns,numCols);
3100                        fclose(fp);
3101                        fp=NULL;
3102                      }
3103                    }
3104                  }
3105                  double * elements = new double [numberColumns+2];
3106                  int * indices = new int [numberColumns];
3107                  if (fp) {
3108                    int numberEntries=1;
3109                    while (numberEntries>0) {
3110                      fpos_t position;
3111                      fgetpos(fp, &position);
3112                      nRead = fread(&numberEntries, sizeof(int), 1, fp);
3113                      if (nRead != 1)
3114                        throw("Error in fread");
3115                      if (numberEntries>0) {
3116                        initialCuts++;
3117                        nRead = fread(elements, sizeof(double), numberEntries+2, fp);
3118                        if (nRead != static_cast<size_t>(numberEntries+2))
3119                          throw("Error in fread");
3120                        nRead = fread(indices, sizeof(int), numberEntries, fp);
3121                        if (nRead != static_cast<size_t>(numberEntries))
3122                          throw("Error in fread");
3123                      } else {
3124                        // end
3125                        fsetpos(fp, &position);
3126                      }
3127                    }
3128                  } else {
3129                    fp = fopen("global.cuts","wb");
3130                    size_t nWrite;
3131                    nWrite=fwrite(&numberColumns,sizeof(int),1,fp);
3132                    if (nWrite != 1)
3133                      throw("Error in fwrite");
3134                  }
3135                  size_t nWrite;
3136                  // now append binding cuts
3137                  int numberC=continuousSolver_->getNumRows();
3138                  int numberRows=solver_->getNumRows();
3139                  printf("Saving %d cuts (up from %d)\n",
3140                         initialCuts+numberRows-numberC,initialCuts);
3141                  const double * rowLower = solver_->getRowLower();
3142                  const double * rowUpper = solver_->getRowUpper();
3143                  // Row copy
3144                  CoinPackedMatrix matrixByRow(*solver_->getMatrixByRow());
3145                  const double * elementByRow = matrixByRow.getElements();
3146                  const int * column = matrixByRow.getIndices();
3147                  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
3148                  const int * rowLength = matrixByRow.getVectorLengths();
3149                  for (int iRow=numberC;iRow<numberRows;iRow++) {
3150                    int n=rowLength[iRow];
3151                    assert (n);
3152                    CoinBigIndex start=rowStart[iRow];
3153                    memcpy(elements,elementByRow+start,n*sizeof(double));
3154                    memcpy(indices,column+start,n*sizeof(int));
3155                    elements[n]=rowLower[iRow];
3156                    elements[n+1]=rowUpper[iRow];
3157                    nWrite=fwrite(&n,sizeof(int),1,fp);
3158                    if (nWrite != 1)
3159                      throw("Error in fwrite");
3160                    nWrite=fwrite(elements,sizeof(double),n+2,fp);
3161                    if (nWrite != static_cast<size_t>(n+2))
3162                      throw("Error in fwrite");
3163                    nWrite=fwrite(indices,sizeof(int),n,fp);
3164                    if (nWrite != static_cast<size_t>(n))
3165                      throw("Error in fwrite");
3166                  }
3167                  // eof marker
3168                  int eofMarker=-1;
3169                  nWrite=fwrite(&eofMarker,sizeof(int),1,fp);
3170                  if (nWrite != 1)
3171                    throw("Error in fwrite");
3172                  fclose(fp);
3173                  // do tighter bounds (? later extra to original columns)
3174                  int nTightened=0;
3175                  const double * lower = solver_->getColLower();
3176                  const double * upper = solver_->getColUpper();
3177                  const double * originalLower = continuousSolver_->getColLower();
3178                  const double * originalUpper = continuousSolver_->getColUpper();
3179                  double tolerance=1.0e-5;
3180                  for (int i=0;i<numberColumns;i++) {
3181                    if (lower[i]>originalLower[i]+tolerance) {
3182                      nTightened++;
3183                    }
3184                    if (upper[i]<originalUpper[i]-tolerance) {
3185                      nTightened++;
3186                    }
3187                  }
3188                  if (nTightened) {
3189                    fp = fopen("global.fix","wb");
3190                    size_t nWrite;
3191                    indices[0]=numberColumns;
3192                    if (originalColumns_)
3193                      indices[1]=COIN_INT_MAX;
3194                    else
3195                      indices[1]=-1;
3196                    nWrite=fwrite(indices,sizeof(int),2,fp);
3197                    if (nWrite != 2)
3198                      throw("Error in fwrite");
3199                    for (int i=0;i<numberColumns;i++) {
3200                      int nTightened=0;
3201                      if (lower[i]>originalLower[i]+tolerance) {
3202                        nTightened++;
3203                      }
3204                      if (upper[i]<originalUpper[i]-tolerance) {
3205                        nTightened++;
3206                      }
3207                      if (nTightened) {
3208                        indices[0]=i;
3209                        if (originalColumns_)
3210                          indices[1]=originalColumns_[i];
3211                        elements[0]=lower[i];
3212                        elements[1]=upper[i];
3213                        elements[2]=originalLower[i];
3214                        elements[3]=originalUpper[i];
3215                        nWrite=fwrite(indices,sizeof(int),2,fp);
3216                        if (nWrite != 2)
3217                          throw("Error in fwrite");
3218                        nWrite=fwrite(elements,sizeof(double),4,fp);
3219                        if (nWrite != 4)
3220                          throw("Error in fwrite");
3221                      }
3222                    }
3223                    // eof marker
3224                    indices[0]=-1;
3225                    nWrite=fwrite(indices,sizeof(int),2,fp);
3226                    if (nWrite != 2)
3227                      throw("Error in fwrite");
3228                    fclose(fp);
3229                  }
3230                  delete [] elements;
3231                  delete [] indices;
3232                }
3233                if ((specialOptions_&524288) != 0 && !parentModel_
3234                        && storedRowCuts_) {
3235                    if (feasible) {
3236                        /* pick up stuff and try again
3237                        add cuts, maybe keep around
3238                        do best solution and if so new heuristics
3239                        obviously tighten bounds
3240                        */
3241                        // A and B probably done on entry
3242                        // A) tight bounds on integer variables
3243                        const double * lower = solver_->getColLower();
3244                        const double * upper = solver_->getColUpper();
3245                        const double * tightLower = storedRowCuts_->tightLower();
3246                        const double * tightUpper = storedRowCuts_->tightUpper();
3247                        int nTightened = 0;
3248                        for (int i = 0; i < numberIntegers_; i++) {
3249                            int iColumn = integerVariable_[i];
3250                            if (tightLower[iColumn] > lower[iColumn]) {
3251                                nTightened++;
3252                                solver_->setColLower(iColumn, tightLower[iColumn]);
3253                            }
3254                            if (tightUpper[iColumn] < upper[iColumn]) {
3255                                nTightened++;
3256                                solver_->setColUpper(iColumn, tightUpper[iColumn]);
3257                            }
3258                        }
3259                        if (nTightened)
3260                          COIN_DETAIL_PRINT(printf("%d tightened by alternate cuts\n", nTightened));
3261                        if (storedRowCuts_->bestObjective() < bestObjective_) {
3262                            // B) best solution
3263                            double objValue = storedRowCuts_->bestObjective();
3264                            setBestSolution(CBC_SOLUTION, objValue,
3265                                            storedRowCuts_->bestSolution()) ;
3266                            // Do heuristics
3267                            // Allow RINS
3268                            for (int i = 0; i < numberHeuristics_; i++) {
3269                                CbcHeuristicRINS * rins
3270                                = dynamic_cast<CbcHeuristicRINS *> (heuristic_[i]);
3271                                if (rins) {
3272                                    rins->setLastNode(-100);
3273                                }
3274                            }
3275                            doHeuristicsAtRoot();
3276                        }
3277#ifdef JJF_ZERO
3278                        int nCuts = storedRowCuts_->sizeRowCuts();
3279                        // add to global list
3280                        for (int i = 0; i < nCuts; i++) {
3281                            OsiRowCut newCut(*storedRowCuts_->rowCutPointer(i));
3282                            newCut.setGloballyValidAsInteger(2);
3283                            newCut.mutableRow().setTestForDuplicateIndex(false);
3284                            globalCuts_.insert(newCut) ;
3285                        }
3286#else
3287                        addCutGenerator(storedRowCuts_, -99, "Stored from previous run",
3288                                        true, false, false, -200);
3289#endif
3290                        // Set cuts as active
3291                        delete [] addedCuts_ ;
3292                        maximumNumberCuts_ = cuts.sizeRowCuts();
3293                        if (maximumNumberCuts_) {
3294                            addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
3295                        } else {
3296                            addedCuts_ = NULL;
3297                        }
3298                        for (int i = 0; i < maximumNumberCuts_; i++)
3299                            addedCuts_[i] = new CbcCountRowCut(*cuts.rowCutPtr(i),
3300                                                               NULL, -1, -1, 2);
3301                        COIN_DETAIL_PRINT(printf("size %d\n", cuts.sizeRowCuts()));
3302                        cuts = OsiCuts();
3303                        currentNumberCuts_ = maximumNumberCuts_;
3304                        feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
3305                                                 NULL);
3306                        for (int i = 0; i < maximumNumberCuts_; i++)
3307                            delete addedCuts_[i];
3308                    }
3309                    delete storedRowCuts_;
3310                    storedRowCuts_ = NULL;
3311                }
3312            } else {
3313                feasible = false;
3314            }
3315        }       else if (solverCharacteristics_->solutionAddsCuts() ||
3316                   solverCharacteristics_->alwaysTryCutsAtRootNode()) {
3317            // may generate cuts and turn the solution
3318            //to an infeasible one
3319            feasible = solveWithCuts(cuts, 2,
3320                                     NULL);
3321        }
3322    }
3323    if (rootModels) {
3324      int numberModels = multipleRootTries_%100;
3325      for (int i=0;i<numberModels;i++)
3326        delete rootModels[i];
3327      delete [] rootModels;
3328    }
3329    // check extra info on feasibility
3330    if (!solverCharacteristics_->mipFeasible())
3331        feasible = false;
3332    // If max nodes==0 - don't do strong branching
3333    if (!getMaximumNodes()) {
3334      if (feasible)
3335        feasible=false;
3336      else
3337        setMaximumNodes(1); //allow to stop on success
3338    }
3339    topOfTree_=NULL;
3340#ifdef CLP_RESOLVE
3341    if ((moreSpecialOptions_&2097152)!=0&&!parentModel_&&feasible) {
3342      OsiClpSolverInterface * clpSolver
3343        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3344      if (clpSolver)
3345        resolveClp(clpSolver,0);
3346    }
3347#endif
3348    // make cut generators less aggressive
3349    for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
3350        CglCutGenerator * generator = generator_[iCutGenerator]->generator();
3351        generator->setAggressiveness(generator->getAggressiveness() - 100);
3352    }
3353    currentNumberCuts_ = numberNewCuts_ ;
3354    if (solverCharacteristics_->solutionAddsCuts()) {
3355      // With some heuristics solver needs a resolve here (don't know if this is bug in heuristics)
3356      solver_->resolve();
3357      if(!isProvenOptimal()){
3358        solver_->initialSolve();
3359      }
3360    }
3361    // See if can stop on gap
3362    bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
3363    if(canStopOnGap()) {
3364        if (bestPossibleObjective_ < getCutoff())
3365            stoppedOnGap_ = true ;
3366        feasible = false;
3367    }
3368    // User event
3369    if (eventHappened_)
3370        feasible = false;
3371#if defined(COIN_HAS_CLP)&&defined(COIN_HAS_CPX)
3372    /*
3373      This is the notion of using Cbc stuff to get going, then calling cplex to
3374      finish off.
3375    */
3376    if (feasible && (specialOptions_&16384) != 0 && fastNodeDepth_ == -2 && !parentModel_) {
3377        // Use Cplex to do search!
3378        double time1 = CoinCpuTime();
3379        OsiClpSolverInterface * clpSolver
3380        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3381        OsiCpxSolverInterface cpxSolver;
3382        double direction = clpSolver->getObjSense();
3383        cpxSolver.setObjSense(direction);
3384        // load up cplex
3385        const CoinPackedMatrix * matrix = continuousSolver_->getMatrixByCol();
3386        const double * rowLower = continuousSolver_->getRowLower();
3387        const double * rowUpper = continuousSolver_->getRowUpper();
3388        const double * columnLower = continuousSolver_->getColLower();
3389        const double * columnUpper = continuousSolver_->getColUpper();
3390        const double * objective = continuousSolver_->getObjCoefficients();
3391        cpxSolver.loadProblem(*matrix, columnLower, columnUpper,
3392                              objective, rowLower, rowUpper);
3393        double * setSol = new double [numberIntegers_];
3394        int * setVar = new int [numberIntegers_];
3395        // cplex doesn't know about objective offset
3396        double offset = clpSolver->getModelPtr()->objectiveOffset();
3397        for (int i = 0; i < numberIntegers_; i++) {
3398            int iColumn = integerVariable_[i];
3399            cpxSolver.setInteger(iColumn);
3400            if (bestSolution_) {
3401                setSol[i] = bestSolution_[iColumn];
3402                setVar[i] = iColumn;
3403            }
3404        }
3405        CPXENVptr env = cpxSolver.getEnvironmentPtr();
3406        CPXLPptr lpPtr = cpxSolver.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL);
3407        cpxSolver.switchToMIP();
3408        if (bestSolution_) {
3409            CPXcopymipstart(env, lpPtr, numberIntegers_, setVar, setSol);
3410        }
3411        if (clpSolver->getNumRows() > continuousSolver_->getNumRows() && false) {
3412            // add cuts
3413            const CoinPackedMatrix * matrix = clpSolver->getMatrixByRow();
3414            const double * rhs = clpSolver->getRightHandSide();
3415            const char * rowSense = clpSolver->getRowSense();
3416            const double * elementByRow = matrix->getElements();
3417            const int * column = matrix->getIndices();
3418            const CoinBigIndex * rowStart = matrix->getVectorStarts();
3419            const int * rowLength = matrix->getVectorLengths();
3420            int nStart = continuousSolver_->getNumRows();
3421            int nRows = clpSolver->getNumRows();
3422            int size = rowStart[nRows-1] + rowLength[nRows-1] -
3423                       rowStart[nStart];
3424            int nAdd = 0;
3425            double * rmatval = new double [size];
3426            int * rmatind = new int [size];
3427            int * rmatbeg = new int [nRows-nStart+1];
3428            size = 0;
3429            rmatbeg[0] = 0;
3430            for (int i = nStart; i < nRows; i++) {
3431                for (int k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
3432                    rmatind[size] = column[k];
3433                    rmatval[size++] = elementByRow[k];
3434                }
3435                nAdd++;
3436                rmatbeg[nAdd] = size;
3437            }
3438            CPXaddlazyconstraints(env, lpPtr, nAdd, size,
3439                                  rhs, rowSense, rmatbeg,
3440                                  rmatind, rmatval, NULL);
3441            CPXsetintparam( env, CPX_PARAM_REDUCE,
3442                            // CPX_PREREDUCE_NOPRIMALORDUAL (0)
3443                            CPX_PREREDUCE_PRIMALONLY);
3444        }
3445        if (getCutoff() < 1.0e50) {
3446            double useCutoff = getCutoff() + offset;
3447            if (bestObjective_ < 1.0e50)
3448                useCutoff = bestObjective_ + offset + 1.0e-7;
3449            cpxSolver.setDblParam(OsiDualObjectiveLimit, useCutoff*
3450                                  direction);
3451            if ( direction > 0.0 )
3452                CPXsetdblparam( env, CPX_PARAM_CUTUP, useCutoff ) ; // min
3453            else
3454                CPXsetdblparam( env, CPX_PARAM_CUTLO, useCutoff ) ; // max
3455        }
3456        CPXsetdblparam(env, CPX_PARAM_EPGAP, dblParam_[CbcAllowableFractionGap]);
3457        delete [] setSol;
3458        delete [] setVar;
3459        char printBuffer[200];
3460        if (offset) {
3461            sprintf(printBuffer, "Add %g to all Cplex messages for true objective",
3462                    -offset);
3463            messageHandler()->message(CBC_GENERAL, messages())
3464            << printBuffer << CoinMessageEol ;
3465            cpxSolver.setDblParam(OsiObjOffset, offset);
3466        }
3467        cpxSolver.branchAndBound();
3468        double timeTaken = CoinCpuTime() - time1;
3469        sprintf(printBuffer, "Cplex took %g seconds",
3470                timeTaken);
3471        messageHandler()->message(CBC_GENERAL, messages())
3472        << printBuffer << CoinMessageEol ;
3473        numberExtraNodes_ = CPXgetnodecnt(env, lpPtr);
3474        numberExtraIterations_ = CPXgetmipitcnt(env, lpPtr);
3475        double value = cpxSolver.getObjValue() * direction;
3476        if (cpxSolver.isProvenOptimal() && value <= getCutoff()) {
3477            feasible = true;
3478            clpSolver->setWarmStart(NULL);
3479            // try and do solution
3480            double * newSolution =
3481                CoinCopyOfArray(cpxSolver.getColSolution(),
3482                                getNumCols());
3483            setBestSolution(CBC_STRONGSOL, value, newSolution) ;
3484            delete [] newSolution;
3485        }
3486        feasible = false;
3487    }
3488#endif
3489    if (!parentModel_&&(moreSpecialOptions_&268435456) != 0) {
3490      // try idiotic idea
3491      CbcObject * obj = new CbcIdiotBranch(this);
3492      obj->setPriority(1); // temp
3493      addObjects(1, &obj);
3494      delete obj;
3495    }
3496   
3497    /*
3498      A hook to use clp to quickly explore some part of the tree.
3499    */
3500    if (fastNodeDepth_ == 1000 &&/*!parentModel_*/(specialOptions_&2048) == 0) {
3501        fastNodeDepth_ = -1;
3502        CbcObject * obj =
3503            new CbcFollowOn(this);
3504        obj->setPriority(1);
3505        addObjects(1, &obj);
3506        delete obj;
3507    }
3508    int saveNumberSolves = numberSolves_;
3509    int saveNumberIterations = numberIterations_;
3510    if ((fastNodeDepth_ >= 0||(moreSpecialOptions_&33554432)!=0)
3511        &&/*!parentModel_*/(specialOptions_&2048) == 0) {
3512        // add in a general depth object doClp
3513        int type = (fastNodeDepth_ <= 100) ? fastNodeDepth_ : -(fastNodeDepth_ - 100);
3514        if ((moreSpecialOptions_&33554432)!=0)
3515          type=12;
3516        else
3517          fastNodeDepth_ += 1000000;     // mark as done
3518        CbcObject * obj =
3519            new CbcGeneralDepth(this, type);
3520        addObjects(1, &obj);
3521        delete obj;
3522        // fake number of objects
3523        numberObjects_--;
3524        if (parallelMode() < -1) {
3525            // But make sure position is correct
3526            OsiObject * obj2 = object_[numberObjects_];
3527            obj = dynamic_cast<CbcObject *> (obj2);
3528            assert (obj);
3529            obj->setPosition(numberObjects_);
3530        }
3531    }
3532#ifdef COIN_HAS_CLP
3533#ifdef NO_CRUNCH
3534    if (true) {
3535        OsiClpSolverInterface * clpSolver
3536        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3537        if (clpSolver && !parentModel_) {
3538            ClpSimplex * clpSimplex = clpSolver->getModelPtr();
3539            clpSimplex->setSpecialOptions(clpSimplex->specialOptions() | 131072);
3540            //clpSimplex->startPermanentArrays();
3541            clpSimplex->setPersistenceFlag(2);
3542        }
3543    }
3544#endif
3545#endif
3546// Save copy of solver
3547    OsiSolverInterface * saveSolver = NULL;
3548    if (!parentModel_ && (specialOptions_&(512 + 32768)) != 0)
3549        saveSolver = solver_->clone();
3550    double checkCutoffForRestart = 1.0e100;
3551    saveModel(saveSolver, &checkCutoffForRestart, &feasible);
3552    if ((specialOptions_&262144) != 0 && !parentModel_) {
3553        // Save stuff and return!
3554        storedRowCuts_->saveStuff(bestObjective_, bestSolution_,
3555                                  solver_->getColLower(),
3556                                  solver_->getColUpper());
3557        delete [] lowerBefore;
3558        delete [] upperBefore;
3559        delete saveSolver;
3560        if (flipObjective)
3561          flipModel();
3562        return;
3563    }
3564    /*
3565      We've taken the continuous relaxation as far as we can. Time to branch.
3566      The first order of business is to actually create a node. chooseBranch
3567      currently uses strong branching to evaluate branch object candidates,
3568      unless forced back to simple branching. If chooseBranch concludes that a
3569      branching candidate is monotone (anyAction == -1) or infeasible (anyAction
3570      == -2) when forced to integer values, it returns here immediately.
3571
3572      Monotone variables trigger a call to resolve(). If the problem remains
3573      feasible, try again to choose a branching variable. At the end of the loop,
3574      resolved == true indicates that some variables were fixed.
3575
3576      Loss of feasibility will result in the deletion of newNode.
3577    */
3578
3579    bool resolved = false ;
3580    CbcNode *newNode = NULL ;
3581    numberFixedAtRoot_ = 0;
3582    numberFixedNow_ = 0;
3583    if (!parentModel_&&(moreSpecialOptions2_&2)!=0) {
3584#ifdef COIN_HAS_CLP
3585      OsiClpSolverInterface * clpSolver
3586        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3587      if (clpSolver) {
3588        if (getCutoff()>1.0e20) {
3589          printf("Zapping costs\n");
3590          int numberColumns=solver_->getNumCols();
3591          double * zeroCost = new double [numberColumns];
3592          // could make small random
3593          memset(zeroCost,0,numberColumns*sizeof(double));
3594          solver_->setObjective(zeroCost);
3595          double objValue = solver_->getObjValue();
3596          solver_->setDblParam(OsiObjOffset,-objValue);
3597          clpSolver->getModelPtr()->setObjectiveValue(objValue);
3598          delete [] zeroCost;
3599        } else {
3600          moreSpecialOptions2_ &= ~2;
3601        }
3602      } else {
3603#endif
3604          moreSpecialOptions2_ &= ~2;
3605#ifdef COIN_HAS_CLP
3606      }
3607#endif
3608    }
3609    int numberIterationsAtContinuous = numberIterations_;
3610    //solverCharacteristics_->setSolver(solver_);
3611    if (feasible) {
3612      // mark all cuts as globally valid
3613      int numberCuts=cuts.sizeRowCuts();
3614      for (int i=0;i<numberCuts;i++) {
3615        cuts.rowCutPtr(i)->setGloballyValid();
3616        whichGenerator_[i]=20000+(whichGenerator_[i]%10000);
3617      }
3618#define HOTSTART -1
3619#if HOTSTART<0
3620        if (bestSolution_ && !parentModel_ && !hotstartSolution_ &&
3621                (moreSpecialOptions_&1024) != 0) {
3622            // Set priorities so only branch on ones we need to
3623            // use djs and see if only few branches needed
3624#ifndef NDEBUG
3625            double integerTolerance = getIntegerTolerance() ;
3626#endif
3627            bool possible = true;
3628            const double * saveLower = continuousSolver_->getColLower();
3629            const double * saveUpper = continuousSolver_->getColUpper();
3630            for (int i = 0; i < numberObjects_; i++) {
3631                const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object_[i]);
3632                if (thisOne) {
3633                    int iColumn = thisOne->columnNumber();
3634                    if (saveUpper[iColumn] > saveLower[iColumn] + 1.5) {
3635                        possible = false;
3636                        break;
3637                    }
3638                } else {
3639                    possible = false;
3640                    break;
3641                }
3642            }
3643            if (possible) {
3644                OsiSolverInterface * solver = continuousSolver_->clone();
3645                int numberColumns = solver->getNumCols();
3646                for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
3647                    double value = bestSolution_[iColumn] ;
3648                    value = CoinMax(value, saveLower[iColumn]) ;
3649                    value = CoinMin(value, saveUpper[iColumn]) ;
3650                    value = floor(value + 0.5);
3651                    if (solver->isInteger(iColumn)) {
3652                        solver->setColLower(iColumn, value);
3653                        solver->setColUpper(iColumn, value);
3654                    }
3655                }
3656                solver->setHintParam(OsiDoDualInResolve, false, OsiHintTry);
3657                // objlim and all slack
3658                double direction = solver->getObjSense();
3659                solver->setDblParam(OsiDualObjectiveLimit, 1.0e50*direction);
3660                CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis *> (solver->getEmptyWarmStart());
3661                solver->setWarmStart(basis);
3662                delete basis;
3663                bool changed = true;
3664                hotstartPriorities_ = new int [numberColumns];
3665                for (int iColumn = 0; iColumn < numberColumns; iColumn++)
3666                    hotstartPriorities_[iColumn] = 1;
3667                while (changed) {
3668                    changed = false;
3669                    solver->resolve();
3670                    if (!solver->isProvenOptimal()) {
3671                        possible = false;
3672                        break;
3673                    }
3674                    const double * dj = solver->getReducedCost();
3675                    const double * colLower = solver->getColLower();
3676                    const double * colUpper = solver->getColUpper();
3677                    const double * solution = solver->getColSolution();
3678                    int nAtLbNatural = 0;
3679                    int nAtUbNatural = 0;
3680                    int nZeroDj = 0;
3681                    int nForced = 0;
3682                    for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
3683                        double value = solution[iColumn] ;
3684                        value = CoinMax(value, saveLower[iColumn]) ;
3685                        value = CoinMin(value, saveUpper[iColumn]) ;
3686                        if (solver->isInteger(iColumn)) {
3687                            assert(fabs(value - solution[iColumn]) <= integerTolerance) ;
3688                            if (hotstartPriorities_[iColumn] == 1) {
3689                                if (dj[iColumn] < -1.0e-6) {
3690                                    // negative dj
3691                                    if (saveUpper[iColumn] == colUpper[iColumn]) {
3692                                        nAtUbNatural++;
3693                                        hotstartPriorities_[iColumn] = 2;
3694                                        solver->setColLower(iColumn, saveLower[iColumn]);
3695                                        solver->setColUpper(iColumn, saveUpper[iColumn]);
3696                                    } else {
3697                                        nForced++;
3698                                    }
3699                                } else if (dj[iColumn] > 1.0e-6) {
3700                                    // positive dj
3701                                    if (saveLower[iColumn] == colLower[iColumn]) {
3702                                        nAtLbNatural++;
3703                                        hotstartPriorities_[iColumn] = 2;
3704                                        solver->setColLower(iColumn, saveLower[iColumn]);
3705                                        solver->setColUpper(iColumn, saveUpper[iColumn]);
3706                                    } else {
3707                                        nForced++;
3708                                    }
3709                                } else {
3710                                    // zero dj
3711                                    nZeroDj++;
3712                                }
3713                            }
3714                        }
3715                    }
3716#if CBC_USEFUL_PRINTING>1
3717                    printf("%d forced, %d naturally at lower, %d at upper - %d zero dj\n",
3718                           nForced, nAtLbNatural, nAtUbNatural, nZeroDj);
3719#endif
3720                    if (nAtLbNatural || nAtUbNatural) {
3721                        changed = true;
3722                    } else {
3723                        if (nForced + nZeroDj > 5000 ||
3724                                (nForced + nZeroDj)*2 > numberIntegers_)
3725                            possible = false;
3726                    }
3727                }
3728                delete solver;
3729            }
3730            if (possible) {
3731                setHotstartSolution(bestSolution_);
3732                if (!saveCompare) {
3733                    // create depth first comparison
3734                    saveCompare = nodeCompare_;
3735                    // depth first
3736                    nodeCompare_ = new CbcCompareDepth();
3737                    tree_->setComparison(*nodeCompare_) ;
3738                }
3739            } else {
3740                delete [] hotstartPriorities_;
3741                hotstartPriorities_ = NULL;
3742            }
3743        }
3744#endif
3745#if HOTSTART>0
3746        if (hotstartSolution_ && !hotstartPriorities_) {
3747            // Set up hot start
3748            OsiSolverInterface * solver = solver_->clone();
3749            double direction = solver_->getObjSense() ;
3750            int numberColumns = solver->getNumCols();
3751            double * saveLower = CoinCopyOfArray(solver->getColLower(), numberColumns);
3752            double * saveUpper = CoinCopyOfArray(solver->getColUpper(), numberColumns);
3753            // move solution
3754            solver->setColSolution(hotstartSolution_);
3755            // point to useful information
3756            const double * saveSolution = testSolution_;
3757            testSolution_ = solver->getColSolution();
3758            OsiBranchingInformation usefulInfo = usefulInformation();
3759            testSolution_ = saveSolution;
3760            /*
3761            Run through the objects and use feasibleRegion() to set variable bounds
3762            so as to fix the variables specified in the objects at their value in this
3763            solution. Since the object list contains (at least) one object for every
3764            integer variable, this has the effect of fixing all integer variables.
3765            */
3766            for (int i = 0; i < numberObjects_; i++)
3767                object_[i]->feasibleRegion(solver, &usefulInfo);
3768            solver->resolve();
3769            assert (solver->isProvenOptimal());
3770            double gap = CoinMax((solver->getObjValue() - solver_->getObjValue()) * direction, 0.0) ;
3771            const double * dj = solver->getReducedCost();
3772            const double * colLower = solver->getColLower();
3773            const double * colUpper = solver->getColUpper();
3774            const double * solution = solver->getColSolution();
3775            int nAtLbNatural = 0;
3776            int nAtUbNatural = 0;
3777            int nAtLbNaturalZero = 0;
3778            int nAtUbNaturalZero = 0;
3779            int nAtLbFixed = 0;
3780            int nAtUbFixed = 0;
3781            int nAtOther = 0;
3782            int nAtOtherNatural = 0;
3783            int nNotNeeded = 0;
3784            delete [] hotstartSolution_;
3785            hotstartSolution_ = new double [numberColumns];
3786            delete [] hotstartPriorities_;
3787            hotstartPriorities_ = new int [numberColumns];
3788            int * order = (int *) saveUpper;
3789            int nFix = 0;
3790            double bestRatio = COIN_DBL_MAX;
3791            for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
3792                double value = solution[iColumn] ;
3793                value = CoinMax(value, saveLower[iColumn]) ;
3794                value = CoinMin(value, saveUpper[iColumn]) ;
3795                double sortValue = COIN_DBL_MAX;
3796                if (solver->isInteger(iColumn)) {
3797                    assert(fabs(value - solution[iColumn]) <= 1.0e-5) ;
3798                    double value2 = floor(value + 0.5);
3799                    if (dj[iColumn] < -1.0e-6) {
3800                        // negative dj
3801                        //assert (value2==colUpper[iColumn]);
3802                        if (saveUpper[iColumn] == colUpper[iColumn]) {
3803                            nAtUbNatural++;
3804                            sortValue = 0.0;
3805                            double value = -dj[iColumn];
3806                            if (value > gap)
3807                                nFix++;
3808                            else if (gap < value*bestRatio)
3809                                bestRatio = gap / value;
3810                            if (saveLower[iColumn] != colLower[iColumn]) {
3811                                nNotNeeded++;
3812                                sortValue = 1.0e20;
3813                            }
3814                        } else if (saveLower[iColumn] == colUpper[iColumn]) {
3815                            nAtLbFixed++;
3816                            sortValue = dj[iColumn];
3817                        } else {
3818                            nAtOther++;
3819                            sortValue = 0.0;
3820                            if (saveLower[iColumn] != colLower[iColumn] &&
3821                                    saveUpper[iColumn] != colUpper[iColumn]) {
3822                                nNotNeeded++;
3823                                sortValue = 1.0e20;
3824                            }
3825                        }
3826                    } else if (dj[iColumn] > 1.0e-6) {
3827                        // positive dj
3828                        //assert (value2==colLower[iColumn]);
3829                        if (saveLower[iColumn] == colLower[iColumn]) {
3830                            nAtLbNatural++;
3831                            sortValue = 0.0;
3832                            double value = dj[iColumn];
3833                            if (value > gap)
3834                                nFix++;
3835                            else if (gap < value*bestRatio)
3836                                bestRatio = gap / value;
3837                            if (saveUpper[iColumn] != colUpper[iColumn]) {
3838                                nNotNeeded++;
3839                                sortValue = 1.0e20;
3840                            }
3841                        } else if (saveUpper[iColumn] == colLower[iColumn]) {
3842                            nAtUbFixed++;
3843                            sortValue = -dj[iColumn];
3844                        } else {
3845                            nAtOther++;
3846                            sortValue = 0.0;
3847                            if (saveLower[iColumn] != colLower[iColumn] &&
3848                                    saveUpper[iColumn] != colUpper[iColumn]) {
3849                                nNotNeeded++;
3850                                sortValue = 1.0e20;
3851                            }
3852                        }
3853                    } else {
3854                        // zero dj
3855                        if (value2 == saveUpper[iColumn]) {
3856                            nAtUbNaturalZero++;
3857                            sortValue = 0.0;
3858                            if (saveLower[iColumn] != colLower[iColumn]) {
3859                                nNotNeeded++;
3860                                sortValue = 1.0e20;
3861                            }
3862                        } else if (value2 == saveLower[iColumn]) {
3863                            nAtLbNaturalZero++;
3864                            sortValue = 0.0;
3865                        } else {
3866                            nAtOtherNatural++;
3867                            sortValue = 0.0;
3868                            if (saveLower[iColumn] != colLower[iColumn] &&
3869                                    saveUpper[iColumn] != colUpper[iColumn]) {
3870                                nNotNeeded++;
3871                                sortValue = 1.0e20;
3872                            }
3873                        }
3874                    }
3875#if HOTSTART==3
3876                    sortValue = -fabs(dj[iColumn]);
3877#endif
3878                }
3879                hotstartSolution_[iColumn] = value ;
3880                saveLower[iColumn] = sortValue;
3881                order[iColumn] = iColumn;
3882            }
3883            COIN_DETAIL_PRINT(printf("** can fix %d columns - best ratio for others is %g on gap of %g\n",
3884                                     nFix, bestRatio, gap));
3885            int nNeg = 0;
3886            CoinSort_2(saveLower, saveLower + numberColumns, order);
3887            for (int i = 0; i < numberColumns; i++) {
3888                if (saveLower[i] < 0.0) {
3889                    nNeg++;
3890#if HOTSTART==2||HOTSTART==3
3891                    // swap sign ?
3892                    saveLower[i] = -saveLower[i];
3893#endif
3894                }
3895            }
3896            CoinSort_2(saveLower, saveLower + nNeg, order);
3897            for (int i = 0; i < numberColumns; i++) {
3898#if HOTSTART==1
3899                hotstartPriorities_[order[i]] = 100;
3900#else
3901                hotstartPriorities_[order[i]] = -(i + 1);
3902#endif
3903            }
3904            COIN_DETAIL_PRINT(printf("nAtLbNat %d,nAtUbNat %d,nAtLbNatZero %d,nAtUbNatZero %d,nAtLbFixed %d,nAtUbFixed %d,nAtOther %d,nAtOtherNat %d, useless %d %d\n",
3905                   nAtLbNatural,
3906                   nAtUbNatural,
3907                   nAtLbNaturalZero,
3908                   nAtUbNaturalZero,
3909                   nAtLbFixed,
3910                   nAtUbFixed,
3911                   nAtOther,
3912                                     nAtOtherNatural, nNotNeeded, nNeg));
3913            delete [] saveLower;
3914            delete [] saveUpper;
3915            if (!saveCompare) {
3916                // create depth first comparison
3917                saveCompare = nodeCompare_;
3918                // depth first
3919                nodeCompare_ = new CbcCompareDepth();
3920                tree_->setComparison(*nodeCompare_) ;
3921            }
3922        }
3923#endif
3924        newNode = new CbcNode ;
3925        // Set objective value (not so obvious if NLP etc)
3926        setObjectiveValue(newNode, NULL);
3927        anyAction = -1 ;
3928        // To make depth available we may need a fake node
3929        CbcNode fakeNode;
3930        if (!currentNode_) {
3931            // Not true if sub trees assert (!numberNodes_);
3932            currentNode_ = &fakeNode;
3933        }
3934        phase_ = 3;
3935        // only allow 1000 passes
3936        int numberPassesLeft = 1000;
3937        // This is first crude step
3938        if (numberAnalyzeIterations_ && !parentModel_) {
3939            delete [] analyzeResults_;
3940            //int numberColumns = solver_->getNumCols();
3941            analyzeResults_ = new double [5*numberIntegers_];
3942            numberFixedAtRoot_ = newNode->analyze(this, analyzeResults_);
3943            if (numberFixedAtRoot_ > 0) {
3944              COIN_DETAIL_PRINT(printf("%d fixed by analysis\n", numberFixedAtRoot_));
3945                setPointers(solver_);
3946                numberFixedNow_ = numberFixedAtRoot_;
3947            } else if (numberFixedAtRoot_ < 0) {
3948              COIN_DETAIL_PRINT(printf("analysis found to be infeasible\n"));
3949                anyAction = -2;
3950                delete newNode ;
3951                newNode = NULL ;
3952                feasible = false ;
3953            }
3954        }
3955        OsiSolverBranch * branches = NULL;
3956        if (feasible)
3957          anyAction = chooseBranch(newNode, numberPassesLeft, NULL, cuts, resolved,
3958                                 NULL, NULL, NULL, branches);
3959        if (anyAction == -2 || newNode->objectiveValue() >= cutoff) {
3960            if (anyAction != -2) {
3961                // zap parent nodeInfo
3962#ifdef COIN_DEVELOP
3963                printf("zapping CbcNodeInfo %x\n", newNode->nodeInfo()->parent());
3964#endif
3965                if (newNode->nodeInfo())
3966                    newNode->nodeInfo()->nullParent();
3967            }
3968            delete newNode ;
3969            newNode = NULL ;
3970            feasible = false ;
3971        }
3972    }
3973    if (newNode && probingInfo_) {
3974        int number01 = probingInfo_->numberIntegers();
3975        //const fixEntry * entry = probingInfo_->fixEntries();
3976        const int * toZero = probingInfo_->toZero();
3977        //const int * toOne = probingInfo_->toOne();
3978        //const int * integerVariable = probingInfo_->integerVariable();
3979        if (toZero[number01]) {
3980            CglTreeProbingInfo info(*probingInfo_);
3981            if ((moreSpecialOptions2_&64)!=0&&!parentModel_) {
3982              /*
3983                Marginal idea. Further exploration probably good. Build some extra
3984                cliques from probing info. Not quite worth the effort?
3985              */
3986              CglProbing generator1;
3987              generator1.setUsingObjective(false);
3988              generator1.setMaxPass(1);
3989              generator1.setMaxPassRoot(1);
3990              generator1.setMaxLook(100);
3991              generator1.setRowCuts(3);
3992              generator1.setMaxElements(300);
3993              generator1.setMaxProbeRoot(solver_->getNumCols());
3994              CoinThreadRandom randomGenerator;
3995              //CglTreeProbingInfo info(solver_);
3996              info.level = 0;
3997              info.formulation_rows = solver_->getNumRows();
3998              info.inTree = false;
3999              info.randomNumberGenerator=&randomGenerator;
4000              info.pass=4;
4001              generator1.setMode(8);
4002              OsiCuts cs;
4003              generator1.generateCutsAndModify(*solver_,cs,&info);
4004              // very clunky
4005              OsiSolverInterface * temp = generator1.cliqueModel(solver_,2);
4006              CglPreProcess dummy;
4007              OsiSolverInterface * newSolver=dummy.cliqueIt(*temp,0.0001);
4008              delete temp;
4009              OsiSolverInterface * fake = NULL;
4010              if (newSolver) {
4011#if 0
4012                int numberCliques = generator1.numberCliques();
4013                cliqueEntry * entry = generator1.cliqueEntry();
4014                cliqueType * type = new cliqueType [numberCliques];
4015                int * start = new int [numberCliques+1];
4016                start[numberCliques]=2*numberCliques;
4017                int n=0;
4018                for (int i=0;i<numberCliques;i++) {
4019                  start[i]=2*i;
4020                  setOneFixesInCliqueEntry(entry[2*i],true);
4021                  setOneFixesInCliqueEntry(entry[2*i+1],true);
4022                  type[i]=0;
4023                }
4024                fake = info.analyze(*solver_, 1,numberCliques,start,
4025                                    entry,type);
4026                delete [] type;
4027                delete [] entry;
4028#else
4029                fake = info.analyze(*newSolver, 1,-1);
4030#endif
4031                delete newSolver;
4032              } else {
4033                fake = info.analyze(*solver_, 1);
4034              }
4035              if (fake) {
4036                //fake->writeMps("fake");
4037                CglFakeClique cliqueGen(fake);
4038                cliqueGen.setStarCliqueReport(false);
4039                cliqueGen.setRowCliqueReport(false);
4040                cliqueGen.setMinViolation(0.1);
4041                addCutGenerator(&cliqueGen, 1, "Fake cliques", true, false, false, -200);
4042                generator_[numberCutGenerators_-1]->setTiming(true);
4043                for (int i = 0; i < numberCutGenerators_; i++) {
4044                  CglKnapsackCover * cutGen =
4045                  dynamic_cast<CglKnapsackCover *>(generator_[i]->generator());
4046                  if (cutGen) {
4047                    cutGen->createCliques(*fake,2,200,false);
4048                  }
4049                }
4050              }
4051            }
4052            if (probingInfo_->packDown()) {
4053#if CBC_USEFUL_PRINTING>1
4054                printf("%d implications on %d 0-1\n", toZero[number01], number01);
4055#endif
4056                // Create a cut generator that remembers implications discovered at root.
4057                CglImplication implication(probingInfo_);
4058                addCutGenerator(&implication, 1, "ImplicationCuts", true, false, false, -200);
4059                generator_[numberCutGenerators_-1]->setGlobalCuts(true);
4060                generator_[numberCutGenerators_-1]->setTiming(true);
4061            } else {
4062                delete probingInfo_;
4063                probingInfo_ = NULL;
4064            }
4065        } else {
4066            delete probingInfo_;
4067
4068            probingInfo_ = NULL;
4069        }
4070    }
4071    /*
4072      At this point, the root subproblem is infeasible or fathomed by bound
4073      (newNode == NULL), or we're live with an objective value that satisfies the
4074      current objective cutoff.
4075    */
4076    assert (!newNode || newNode->objectiveValue() <= cutoff) ;
4077    // Save address of root node as we don't want to delete it
4078    /*
4079      The common case is that the lp relaxation is feasible but doesn't satisfy
4080      integrality (i.e., newNode->branchingObject(), indicating we've been able to
4081      select a branching variable). Remove any cuts that have gone slack due to
4082      forcing monotone variables. Then tack on an CbcFullNodeInfo object and full
4083      basis (via createInfo()) and stash the new cuts in the nodeInfo (via
4084      addCuts()). If, by some miracle, we have an integral solution at the root
4085      (newNode->branchingObject() is NULL), takeOffCuts() will ensure that the solver holds
4086      a valid solution for use by setBestSolution().
4087    */
4088    CoinWarmStartBasis *lastws = NULL ;
4089    if (feasible && newNode->branchingObject()) {
4090        if (resolved) {
4091            takeOffCuts(cuts, false, NULL) ;
4092#     ifdef CHECK_CUT_COUNTS
4093            { printf("Number of rows after chooseBranch fix (root)"
4094                         "(active only) %d\n",
4095                         numberRowsAtContinuous_ + numberNewCuts_ + numberOldActiveCuts_) ;
4096                const CoinWarmStartBasis* debugws =
4097                    dynamic_cast <const CoinWarmStartBasis*>(solver_->getWarmStart()) ;
4098                debugws->print() ;
4099                delete debugws ;
4100            }
4101#     endif
4102        }
4103        //newNode->createInfo(this,NULL,NULL,NULL,NULL,0,0) ;
4104        //newNode->nodeInfo()->addCuts(cuts,
4105        //                       newNode->numberBranches(),whichGenerator_) ;
4106        if (lastws) delete lastws ;
4107        lastws = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
4108    }
4109    /*
4110      Continuous data to be used later
4111    */
4112    continuousObjective_ = solver_->getObjValue() * solver_->getObjSense();
4113    continuousInfeasibilities_ = 0 ;
4114    if (newNode) {
4115        continuousObjective_ = newNode->objectiveValue() ;
4116        delete [] continuousSolution_;
4117        continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
4118                                              numberColumns);
4119        continuousInfeasibilities_ = newNode->numberUnsatisfied() ;
4120    }
4121    /*
4122      Bound may have changed so reset in objects
4123    */
4124    {
4125        int i ;
4126        for (i = 0; i < numberObjects_; i++)
4127            object_[i]->resetBounds(solver_) ;
4128    }
4129    /*
4130      Feasible? Then we should have either a live node prepped for future
4131      expansion (indicated by variable() >= 0), or (miracle of miracles) an
4132      integral solution at the root node.
4133
4134      initializeInfo sets the reference counts in the nodeInfo object.  Since
4135      this node is still live, push it onto the heap that holds the live set.
4136    */
4137    if (newNode) {
4138        if (newNode->branchingObject()) {
4139            newNode->initializeInfo() ;
4140            tree_->push(newNode) ;
4141            // save pointer to root node - so can pick up bounds
4142            if (!topOfTree_)
4143              topOfTree_ = dynamic_cast<CbcFullNodeInfo *>(newNode->nodeInfo()) ;
4144            if (statistics_) {
4145                if (numberNodes2_ == maximumStatistics_) {
4146                    maximumStatistics_ = 2 * maximumStatistics_;
4147                    CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_];
4148                    memset(temp, 0, maximumStatistics_*sizeof(CbcStatistics *));
4149                    memcpy(temp, statistics_, numberNodes2_*sizeof(CbcStatistics *));
4150                    delete [] statistics_;
4151                    statistics_ = temp;
4152                }
4153                assert (!statistics_[numberNodes2_]);
4154                statistics_[numberNodes2_] = new CbcStatistics(newNode, this);
4155            }
4156            numberNodes2_++;
4157#     ifdef CHECK_NODE
4158            printf("Node %x on tree\n", newNode) ;
4159#     endif
4160        } else {
4161            // continuous is integer
4162            double objectiveValue = newNode->objectiveValue();
4163            setBestSolution(CBC_SOLUTION, objectiveValue,
4164                            solver_->getColSolution()) ;
4165            if (eventHandler) {
4166              // we are stopping anyway so no need to test return code
4167              eventHandler->event(CbcEventHandler::solution);
4168            }
4169            delete newNode ;
4170            newNode = NULL ;
4171        }
4172    }
4173
4174    if (printFrequency_ <= 0) {
4175        printFrequency_ = 1000 ;
4176        if (getNumCols() > 2000)
4177            printFrequency_ = 100 ;
4178    }
4179    /*
4180      It is possible that strong branching fixes one variable and then the code goes round
4181      again and again.  This can take too long.  So we need to warn user - just once.
4182    */
4183    numberLongStrong_ = 0;
4184    CbcNode * createdNode = NULL;
4185#ifdef CBC_THREAD
4186    if ((specialOptions_&2048) != 0)
4187        numberThreads_ = 0;
4188    if (numberThreads_ ) {
4189        nodeCompare_->sayThreaded(); // need to use addresses
4190        master_ = new CbcBaseModel(*this,
4191                                   (parallelMode() < -1) ? 1 : 0);
4192        masterThread_ = master_->masterThread();
4193    }
4194#endif
4195#ifdef COIN_HAS_CLP
4196    {
4197        OsiClpSolverInterface * clpSolver
4198        = dynamic_cast<OsiClpSolverInterface *> (solver_);
4199        if (clpSolver && !parentModel_) {
4200            clpSolver->computeLargestAway();
4201        }
4202    }
4203#endif
4204    /*
4205      At last, the actual branch-and-cut search loop, which will iterate until
4206      the live set is empty or we hit some limit (integrality gap, time, node
4207      count, etc.). The overall flow is to rebuild a subproblem, reoptimise using
4208      solveWithCuts(), choose a branching pattern with chooseBranch(), and finally
4209      add the node to the live set.
4210
4211      The first action is to winnow the live set to remove nodes which are worse
4212      than the current objective cutoff.
4213    */
4214    if (solver_->getRowCutDebuggerAlways()) {
4215        OsiRowCutDebugger * debuggerX = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
4216        const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger() ;
4217        if (!debugger) {
4218            // infeasible!!
4219            printf("before search\n");
4220            const double * lower = solver_->getColLower();
4221            const double * upper = solver_->getColUpper();
4222            const double * solution = debuggerX->optimalSolution();
4223            int numberColumns = solver_->getNumCols();
4224            for (int i = 0; i < numberColumns; i++) {
4225                if (solver_->isInteger(i)) {
4226                    if (solution[i] < lower[i] - 1.0e-6 || solution[i] > upper[i] + 1.0e-6)
4227                        printf("**** ");
4228                    printf("%d %g <= %g <= %g\n",
4229                           i, lower[i], solution[i], upper[i]);
4230                }
4231            }
4232            //abort();
4233        }
4234    }
4235    {
4236        // may be able to change cutoff now
4237        double cutoff = getCutoff();
4238        double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
4239        if (cutoff > bestObjective_ - increment) {
4240            cutoff = bestObjective_ - increment ;
4241            setCutoff(cutoff) ;
4242        }
4243    }
4244#ifdef CBC_THREAD
4245    bool goneParallel = false;
4246#endif
4247#define MAX_DEL_NODE 1
4248    CbcNode * delNode[MAX_DEL_NODE+1];
4249    int nDeleteNode = 0;
4250    // For Printing etc when parallel
4251    int lastEvery1000 = 0;
4252    int lastPrintEvery = 0;
4253    int numberConsecutiveInfeasible = 0;
4254#define PERTURB_IN_FATHOM
4255#ifdef PERTURB_IN_FATHOM
4256    // allow in fathom
4257    if ((moreSpecialOptions_& 262144) != 0)
4258      specialOptions_ |= 131072;
4259#endif
4260    while (true) {
4261        lockThread();
4262#ifdef COIN_HAS_CLP
4263        // See if we want dantzig row choice
4264        goToDantzig(100, savePivotMethod);
4265#endif
4266        if (tree_->empty()) {
4267#ifdef CBC_THREAD
4268            if (parallelMode() > 0 && master_) {
4269                int anyLeft = master_->waitForThreadsInTree(0);
4270                if (!anyLeft) {
4271                    master_->stopThreads(-1);
4272                    break;
4273                }
4274            } else {
4275                break;
4276            }
4277#else
4278            break;
4279#endif
4280        } else {
4281            unlockThread();
4282        }
4283        // If done 50/100 nodes see if worth trying reduction
4284        if (numberNodes_ >= nextCheckRestart) {
4285          if (nextCheckRestart<100)
4286            nextCheckRestart=100;
4287          else
4288            nextCheckRestart=COIN_INT_MAX;
4289#ifdef COIN_HAS_CLP
4290            OsiClpSolverInterface * clpSolver
4291            = dynamic_cast<OsiClpSolverInterface *> (solver_);
4292            if (clpSolver && ((specialOptions_&131072) == 0) && true) {
4293                ClpSimplex * simplex = clpSolver->getModelPtr();
4294                int perturbation = simplex->perturbation();
4295#if CBC_USEFUL_PRINTING>1
4296                printf("Testing its n,s %d %d solves n,s %d %d - pert %d\n",
4297                       numberIterations_, saveNumberIterations,
4298                       numberSolves_, saveNumberSolves, perturbation);
4299#endif
4300                if (perturbation == 50 && (numberIterations_ - saveNumberIterations) <
4301                        8*(numberSolves_ - saveNumberSolves)) {
4302                    // switch off perturbation
4303                    simplex->setPerturbation(100);
4304#if CBC_USEFUL_PRINTING>1
4305                    printf("Perturbation switched off\n");
4306#endif
4307                }
4308            }
4309#endif
4310            /*
4311              Decide if we want to do a restart.
4312            */
4313            if (saveSolver && (specialOptions_&(512 + 32768)) != 0) {
4314                bool tryNewSearch = solverCharacteristics_->reducedCostsAccurate() &&
4315                                    (getCutoff() < 1.0e20 && getCutoff() < checkCutoffForRestart);
4316                int numberColumns = getNumCols();
4317                if (tryNewSearch) {
4318                    // adding increment back allows current best - tiny bit weaker
4319                    checkCutoffForRestart = getCutoff() + getCutoffIncrement() ;
4320#if CBC_USEFUL_PRINTING>1
4321                    printf("after %d nodes, cutoff %g - looking\n",
4322                           numberNodes_, getCutoff());
4323#endif
4324                    saveSolver->resolve();
4325                    double direction = saveSolver->getObjSense() ;
4326                    double gap = checkCutoffForRestart - saveSolver->getObjValue() * direction ;
4327                    double tolerance;
4328                    saveSolver->getDblParam(OsiDualTolerance, tolerance) ;
4329                    if (gap <= 0.0)
4330                        gap = tolerance;
4331                    gap += 100.0 * tolerance;
4332                    double integerTolerance = getDblParam(CbcIntegerTolerance) ;
4333
4334                    const double *lower = saveSolver->getColLower() ;
4335                    const double *upper = saveSolver->getColUpper() ;
4336                    const double *solution = saveSolver->getColSolution() ;
4337                    const double *reducedCost = saveSolver->getReducedCost() ;
4338
4339                    int numberFixed = 0 ;
4340                    int numberFixed2 = 0;
4341#ifdef COIN_DEVELOP
4342                    printf("gap %g\n", gap);
4343#endif
4344                    for (int i = 0 ; i < numberIntegers_ ; i++) {
4345                        int iColumn = integerVariable_[i] ;
4346                        double djValue = direction * reducedCost[iColumn] ;
4347                        if (upper[iColumn] - lower[iColumn] > integerTolerance) {
4348                            if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) {
4349                                //printf("%d to lb on dj of %g - bounds %g %g\n",
4350                                //     iColumn,djValue,lower[iColumn],upper[iColumn]);
4351                                saveSolver->setColUpper(iColumn, lower[iColumn]) ;
4352                                numberFixed++ ;
4353                            } else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) {
4354                                //printf("%d to ub on dj of %g - bounds %g %g\n",
4355                                //     iColumn,djValue,lower[iColumn],upper[iColumn]);
4356                                saveSolver->setColLower(iColumn, upper[iColumn]) ;
4357                                numberFixed++ ;
4358                            }
4359                        } else {
4360                            //printf("%d has dj of %g - already fixed to %g\n",
4361                            //     iColumn,djValue,lower[iColumn]);
4362                            numberFixed2++;
4363                        }
4364                    }
4365#ifdef COIN_DEVELOP
4366                    if ((specialOptions_&1) != 0) {
4367                        const OsiRowCutDebugger *debugger = saveSolver->getRowCutDebugger() ;
4368                        if (debugger) {
4369                            printf("Contains optimal\n") ;
4370                            OsiSolverInterface * temp = saveSolver->clone();
4371                            const double * solution = debugger->optimalSolution();
4372                            const double *lower = temp->getColLower() ;
4373                            const double *upper = temp->getColUpper() ;
4374                            int n = temp->getNumCols();
4375                            for (int i = 0; i < n; i++) {
4376                                if (temp->isInteger(i)) {
4377                                    double value = floor(solution[i] + 0.5);
4378                                    assert (value >= lower[i] && value <= upper[i]);
4379                                    temp->setColLower(i, value);
4380                                    temp->setColUpper(i, value);
4381                                }
4382                            }
4383                            temp->writeMps("reduced_fix");
4384                            delete temp;
4385                            saveSolver->writeMps("reduced");
4386                        } else {
4387                            abort();
4388                        }
4389                    }
4390                    printf("Restart could fix %d integers (%d already fixed)\n",
4391                           numberFixed + numberFixed2, numberFixed2);
4392#endif
4393                    numberFixed += numberFixed2;
4394                    if (numberFixed*10 < numberColumns && numberFixed*4 < numberIntegers_)
4395                        tryNewSearch = false;
4396                }
4397#ifdef CONFLICT_CUTS
4398                // temporary
4399                if ((moreSpecialOptions_&4194304)!=0)
4400                  tryNewSearch=false;
4401#endif
4402                if (tryNewSearch) {
4403                    // back to solver without cuts?
4404                    OsiSolverInterface * solver2 = saveSolver->clone();
4405                    const double *lower = saveSolver->getColLower() ;
4406                    const double *upper = saveSolver->getColUpper() ;
4407                    for (int i = 0 ; i < numberIntegers_ ; i++) {
4408                        int iColumn = integerVariable_[i] ;
4409                        solver2->setColLower(iColumn, lower[iColumn]);
4410                        solver2->setColUpper(iColumn, upper[iColumn]);
4411                    }
4412                    // swap
4413                    delete saveSolver;
4414                    saveSolver = solver2;
4415                    double * newSolution = new double[numberColumns];
4416                    double objectiveValue = checkCutoffForRestart;
4417                    // Save the best solution so far.
4418                    CbcSerendipity heuristic(*this);
4419                    if (bestSolution_)
4420                        heuristic.setInputSolution(bestSolution_, bestObjective_);
4421                    // Magic number
4422                    heuristic.setFractionSmall(0.8);
4423                    // `pumpTune' to stand-alone solver for explanations.
4424                    heuristic.setFeasibilityPumpOptions(1008013);
4425                    // Use numberNodes to say how many are original rows
4426                    heuristic.setNumberNodes(continuousSolver_->getNumRows());
4427#ifdef COIN_DEVELOP
4428                    if (continuousSolver_->getNumRows() <
4429                            solver_->getNumRows())
4430                        printf("%d rows added ZZZZZ\n",
4431                               solver_->getNumRows() - continuousSolver_->getNumRows());
4432#endif
4433                    int returnCode = heuristic.smallBranchAndBound(saveSolver,
4434                                     -1, newSolution,
4435                                     objectiveValue,
4436                                     checkCutoffForRestart, "Reduce");
4437                    if (returnCode < 0) {
4438#ifdef COIN_DEVELOP
4439                        printf("Restart - not small enough to do search after fixing\n");
4440#endif
4441                        delete [] newSolution;
4442                    } else {
4443                        // 1 for sol'n, 2 for finished, 3 for both
4444                        if ((returnCode&1) != 0) {
4445                            // increment number of solutions so other heuristics can test
4446                            numberSolutions_++;
4447                            numberHeuristicSolutions_++;
4448                            lastHeuristic_ = NULL;
4449                            setBestSolution(CBC_ROUNDING, objectiveValue, newSolution) ;
4450                        }
4451                        delete [] newSolution;
4452#ifdef CBC_THREAD
4453                        if (master_) {
4454                            lockThread();
4455                            if (parallelMode() > 0) {
4456                                while (master_->waitForThreadsInTree(0)) {
4457                                    lockThread();
4458                                    double dummyBest;
4459                                    tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4460                                    //unlockThread();
4461                                }
4462                            } else {
4463                                double dummyBest;
4464                                tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4465                            }
4466                            master_->waitForThreadsInTree(2);
4467                            delete master_;
4468                            master_ = NULL;
4469                            masterThread_ = NULL;
4470                        }
4471#endif
4472                        if (tree_->size()) {
4473                            double dummyBest;
4474                            tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4475                        }
4476                        break;
4477                    }
4478                }
4479                delete saveSolver;
4480                saveSolver = NULL;
4481            }
4482        }
4483        /*
4484          Check for abort on limits: node count, solution count, time, integrality gap.
4485        */
4486        if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
4487                numberSolutions_ < intParam_[CbcMaxNumSol] &&
4488                !maximumSecondsReached() &&
4489                !stoppedOnGap_ && !eventHappened_ && (maximumNumberIterations_ < 0 ||
4490                                                      numberIterations_ < maximumNumberIterations_))) {
4491            // out of loop
4492            break;
4493        }
4494#ifdef BONMIN
4495        assert(!solverCharacteristics_->solutionAddsCuts() || solverCharacteristics_->mipFeasible());
4496#endif
4497// Sets percentage of time when we try diving. Diving requires a bit of heap reorganisation, because
4498// we need to replace the comparison function to dive, and that requires reordering to retain the
4499// heap property.
4500#define DIVE_WHEN 1000
4501#define DIVE_STOP 2000
4502        int kNode = numberNodes_ % 4000;
4503        if (numberNodes_<100000 && kNode>DIVE_WHEN && kNode <= DIVE_STOP) {
4504            if (!parallelMode()) {
4505                if (kNode == DIVE_WHEN + 1 || numberConsecutiveInfeasible > 1) {
4506                    CbcCompareDefault * compare = dynamic_cast<CbcCompareDefault *>
4507                                                  (nodeCompare_);
4508                    // Don't interfere if user has replaced the compare function.
4509                    if (compare) {
4510                        //printf("Redoing tree\n");
4511                        compare->startDive(this);
4512                        numberConsecutiveInfeasible = 0;
4513                    }
4514                }
4515            }
4516        }
4517        // replace current cutoff?
4518        if (cutoff > getCutoff()) {
4519            double newCutoff = getCutoff();
4520            if (analyzeResults_) {
4521                // see if we could fix any (more)
4522                int n = 0;
4523                double * newLower = analyzeResults_;
4524                double * objLower = newLower + numberIntegers_;
4525                double * newUpper = objLower + numberIntegers_;
4526                double * objUpper = newUpper + numberIntegers_;
4527                for (int i = 0; i < numberIntegers_; i++) {
4528                    if (objLower[i] > newCutoff) {
4529                        n++;
4530                        if (objUpper[i] > newCutoff) {
4531                            newCutoff = -COIN_DBL_MAX;
4532                            break;
4533                        }
4534                        // add as global cut
4535                        objLower[i]=-COIN_DBL_MAX;
4536                        OsiRowCut rc;
4537                        rc.setLb(newLower[i]);
4538                        rc.setUb(COIN_DBL_MAX);
4539                        double one=1.0;
4540                        rc.setRow(1,integerVariable_+i,&one,false);
4541                        rc.setGloballyValidAsInteger(2);
4542                        globalCuts_.addCutIfNotDuplicate(rc) ;
4543                    } else if (objUpper[i] > newCutoff) {
4544                        n++;
4545                        // add as global cut
4546                        objUpper[i]=-COIN_DBL_MAX;
4547                        OsiRowCut rc;
4548                        rc.setLb(-COIN_DBL_MAX);
4549                        rc.setUb(newUpper[i]);
4550                        double one=1.0;
4551                        rc.setRow(1,integerVariable_+i,&one,false);
4552                        rc.setGloballyValidAsInteger(2);
4553                        globalCuts_.addCutIfNotDuplicate(rc) ;
4554                    }
4555                }
4556                if (newCutoff == -COIN_DBL_MAX) {
4557                  COIN_DETAIL_PRINT(printf("Root analysis says finished\n"));
4558                } else if (n > numberFixedNow_) {
4559                  COIN_DETAIL_PRINT(printf("%d more fixed by analysis - now %d\n", n - numberFixedNow_, n));
4560                    numberFixedNow_ = n;
4561                }
4562            }
4563            if (eventHandler) {
4564                if (!eventHandler->event(CbcEventHandler::solution)) {
4565                    eventHappened_ = true; // exit
4566                }
4567                newCutoff = getCutoff();
4568            }
4569            lockThread();
4570            /*
4571              Clean the tree to reflect the new solution, then see if the
4572              node comparison predicate wants to make any changes. If so,
4573              call setComparison for the side effect of rebuilding the heap.
4574            */
4575            tree_->cleanTree(this,newCutoff,bestPossibleObjective_) ;
4576            if (nodeCompare_->newSolution(this) ||
4577                nodeCompare_->newSolution(this,continuousObjective_,
4578                                          continuousInfeasibilities_)) {
4579              tree_->setComparison(*nodeCompare_) ;
4580            }
4581            if (tree_->empty()) {
4582                continue;
4583            }
4584            unlockThread();
4585        }
4586        cutoff = getCutoff() ;
4587        /*
4588            Periodic activities: Opportunities to
4589            + tweak the nodeCompare criteria,
4590            + check if we've closed the integrality gap enough to quit,
4591            + print a summary line to let the user know we're working
4592        */
4593        if (numberNodes_ >= lastEvery1000) {
4594            lockThread();
4595#ifdef COIN_HAS_CLP
4596            // See if we want dantzig row choice
4597            goToDantzig(1000, savePivotMethod);
4598#endif
4599            lastEvery1000 = numberNodes_ + 1000;
4600            bool redoTree = nodeCompare_->every1000Nodes(this, numberNodes_) ;
4601#ifdef CHECK_CUT_SIZE
4602            verifyCutSize (tree_, *this);
4603#endif
4604            // redo tree if requested
4605            if (redoTree)
4606                tree_->setComparison(*nodeCompare_) ;
4607            unlockThread();
4608        }
4609        // Had hotstart before, now switched off
4610        if (saveCompare && !hotstartSolution_) {
4611            // hotstart switched off
4612            delete nodeCompare_; // off depth first
4613            nodeCompare_ = saveCompare;
4614            saveCompare = NULL;
4615            // redo tree
4616            lockThread();
4617            tree_->setComparison(*nodeCompare_) ;
4618            unlockThread();
4619        }
4620        if (numberNodes_ >= lastPrintEvery) {
4621            lastPrintEvery = numberNodes_ + printFrequency_;
4622            lockThread();
4623            int nNodes = tree_->size() ;
4624
4625            //MODIF PIERRE
4626            bestPossibleObjective_ = tree_->getBestPossibleObjective();
4627#ifdef CBC_THREAD
4628            if (parallelMode() > 0 && master_) {
4629              // need to adjust for ones not on tree
4630              int numberThreads = master_->numberThreads();
4631              for (int i=0;i<numberThreads;i++) {
4632                CbcThread * child = master_->child(i);
4633                if (child->node()) {
4634                  // adjust
4635                  double value = child->node()->objectiveValue();
4636                  bestPossibleObjective_ = CoinMin(bestPossibleObjective_, value);
4637                }
4638              }
4639            }
4640#endif
4641            unlockThread();
4642#if CBC_USEFUL_PRINTING>1
4643            if (getCutoff() < 1.0e20) {
4644                if (fabs(getCutoff() - (bestObjective_ - getCutoffIncrement())) > 1.0e-6 &&
4645                        !parentModel_)
4646                    printf("model cutoff in status %g, best %g, increment %g\n",
4647                           getCutoff(), bestObjective_, getCutoffIncrement());
4648                assert (getCutoff() < bestObjective_ - getCutoffIncrement() +
4649                        1.0e-6 + 1.0e-10*fabs(bestObjective_));
4650            }
4651#endif
4652            if (!intParam_[CbcPrinting]) {
4653                messageHandler()->message(CBC_STATUS, messages())
4654                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4655                << getCurrentSeconds()
4656                << CoinMessageEol ;
4657            } else if (intParam_[CbcPrinting] == 1) {
4658                messageHandler()->message(CBC_STATUS2, messages())
4659                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4660                << tree_->lastDepth() << tree_->lastUnsatisfied() 
4661                << tree_->lastObjective() << numberIterations_
4662                << getCurrentSeconds()
4663                << CoinMessageEol ;
4664            } else if (!numberExtraIterations_) {
4665                messageHandler()->message(CBC_STATUS2, messages())
4666                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4667                << tree_->lastDepth() << tree_->lastUnsatisfied() << numberIterations_
4668                << getCurrentSeconds()
4669                << CoinMessageEol ;
4670            } else {
4671                messageHandler()->message(CBC_STATUS3, messages())
4672                << numberNodes_ << numberFathoms_ << numberExtraNodes_ << nNodes
4673                << bestObjective_ << bestPossibleObjective_
4674                << tree_->lastDepth() << tree_->lastUnsatisfied() << numberIterations_ << numberExtraIterations_
4675                << getCurrentSeconds()
4676                << CoinMessageEol ;
4677            }
4678#ifdef COIN_HAS_NTY
4679            if (symmetryInfo_) 
4680              symmetryInfo_->statsOrbits(this,1);
4681#endif
4682#if PRINT_CONFLICT==1
4683            if (numberConflictCuts>lastNumberConflictCuts) {
4684              double length = lengthConflictCuts/numberConflictCuts;
4685              printf("%d new conflict cuts - total %d - average length %g\n",
4686                     numberConflictCuts-lastNumberConflictCuts,
4687                     numberConflictCuts,length);
4688              lastNumberConflictCuts = numberConflictCuts;
4689            }
4690#endif
4691            if (eventHandler && !eventHandler->event(CbcEventHandler::treeStatus)) {
4692                eventHappened_ = true; // exit
4693            }
4694        }
4695        // See if can stop on gap
4696        if(canStopOnGap()) {
4697            stoppedOnGap_ = true ;
4698        }
4699
4700#ifdef CHECK_NODE_FULL
4701        verifyTreeNodes(tree_, *this) ;
4702#   endif
4703#   ifdef CHECK_CUT_COUNTS
4704        verifyCutCounts(tree_, *this) ;
4705#   endif
4706        /*
4707          Now we come to the meat of the loop. To create the active subproblem, we'll
4708          pop the most promising node in the live set, rebuild the subproblem it
4709          represents, and then execute the current arm of the branch to create the
4710          active subproblem.
4711        */
4712        CbcNode * node = NULL;
4713#ifdef CBC_THREAD
4714        if (!parallelMode() || parallelMode() == -1) {
4715#endif
4716            node = tree_->bestNode(cutoff) ;
4717            // Possible one on tree worse than cutoff
4718            // Weird comparison function can leave ineligible nodes on tree
4719            if (!node || node->objectiveValue() > cutoff)
4720                continue;
4721            // Do main work of solving node here
4722            doOneNode(this, node, createdNode);
4723#ifdef JJF_ZERO
4724            if (node) {
4725              if (createdNode) {
4726                printf("Node %d depth %d, created %d depth %d\n",
4727                       node->nodeNumber(), node->depth(),
4728                       createdNode->nodeNumber(), createdNode->depth());
4729              } else {
4730                printf("Node %d depth %d,  no created node\n",
4731                       node->nodeNumber(), node->depth());
4732              }
4733            } else if (createdNode) {
4734              printf("Node exhausted, created %d depth %d\n",
4735                     createdNode->nodeNumber(), createdNode->depth());
4736            } else {
4737              printf("Node exhausted,  no created node\n");
4738              numberConsecutiveInfeasible = 2;
4739            }
4740#endif
4741            //if (createdNode)
4742            //numberConsecutiveInfeasible=0;
4743            //else
4744            //numberConsecutiveInfeasible++;
4745#ifdef CBC_THREAD
4746        } else if (parallelMode() > 0) {
4747            //lockThread();
4748            //node = tree_->bestNode(cutoff) ;
4749            // Possible one on tree worse than cutoff
4750            if (true || !node || node->objectiveValue() > cutoff) {
4751                assert (master_);
4752                if (master_) {
4753                    int anyLeft = master_->waitForThreadsInTree(1);
4754                    // may need to go round again
4755                    if (anyLeft) {
4756                        continue;
4757                    } else {
4758                        master_->stopThreads(-1);
4759                    }
4760                }
4761            }
4762            //unlockThread();
4763        } else {
4764            // Deterministic parallel
4765          if ((tree_->size() < CoinMax(numberThreads_, 8)||
4766               hotstartSolution_) && !goneParallel) {
4767                node = tree_->bestNode(cutoff) ;
4768                // Possible one on tree worse than cutoff
4769                if (!node || node->objectiveValue() > cutoff)
4770                    continue;
4771                // Do main work of solving node here
4772                doOneNode(this, node, createdNode);
4773                assert (createdNode);
4774                if (!createdNode->active()) {
4775                    delete createdNode;
4776                    createdNode = NULL;
4777                } else {
4778                    // Say one more pointing to this
4779                    node->nodeInfo()->increment() ;
4780                    tree_->push(createdNode) ;
4781                }
4782                if (node->active()) {
4783                    assert (node->nodeInfo());
4784                    if (node->nodeInfo()->numberBranchesLeft()) {
4785                        tree_->push(node) ;
4786                    } else {
4787                        node->setActive(false);
4788                    }
4789                } else {
4790                    if (node->nodeInfo()) {
4791                        if (!node->nodeInfo()->numberBranchesLeft())
4792                            node->nodeInfo()->allBranchesGone(); // can clean up
4793                        // So will delete underlying stuff
4794                        node->setActive(true);
4795                    }
4796                    delNode[nDeleteNode++] = node;
4797                    node = NULL;
4798                }
4799                if (nDeleteNode >= MAX_DEL_NODE) {
4800                    for (int i = 0; i < nDeleteNode; i++) {
4801                        //printf("trying to del %d %x\n",i,delNode[i]);
4802                        delete delNode[i];
4803                        //printf("done to del %d %x\n",i,delNode[i]);
4804                    }
4805                    nDeleteNode = 0;
4806                }
4807            } else {
4808                // Split and solve
4809                master_->deterministicParallel();
4810                goneParallel = true;
4811            }
4812        }
4813#endif
4814    }
4815    if (nDeleteNode) {
4816        for (int i = 0; i < nDeleteNode; i++) {
4817            delete delNode[i];
4818        }
4819        nDeleteNode = 0;
4820    }
4821#ifdef CBC_THREAD
4822    if (master_) {
4823        master_->stopThreads(-1);
4824        master_->waitForThreadsInTree(2);
4825        // adjust time to allow for children on some systems
4826        //dblParam_[CbcStartSeconds] -= CoinCpuTimeJustChildren();
4827    }
4828#endif
4829    /*
4830      End of the non-abort actions. The next block of code is executed if we've
4831      aborted because we hit one of the limits. Clean up by deleting the live set
4832      and break out of the node processing loop. Note that on an abort, node may
4833      have been pushed back onto the tree for further processing, in which case
4834      it'll be deleted in cleanTree. We need to check.
4835    */
4836    if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
4837          numberSolutions_ < intParam_[CbcMaxNumSol] &&
4838          !maximumSecondsReached() &&
4839          !stoppedOnGap_ &&
4840          !eventHappened_ &&
4841          (maximumNumberIterations_ < 0 || numberIterations_ < maximumNumberIterations_))
4842         ) {
4843        if (tree_->size()) {
4844            double dummyBest;
4845            tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4846#if 0 // Does not seem to be needed def CBC_THREAD
4847            if (parallelMode() > 0 && master_) {
4848              // see if any dangling nodes
4849              int numberThreads = master_->numberThreads();
4850              for (int i=0;i<numberThreads;i++) {
4851                CbcThread * child = master_->child(i);
4852                //if (child->createdNode())
4853                //printf("CHILD_NODE %p\n",child->createdNode());
4854                delete child->createdNode();
4855              }
4856            }
4857#endif
4858        }
4859        delete nextRowCut_;
4860        /* order is important here:
4861         * maximumSecondsReached() should be checked before eventHappened_ and
4862         * isNodeLimitReached() should be checked after eventHappened_
4863         * reason is, that at timelimit, eventHappened_ is set to true to make Cbc stop fast
4864         *   and if Ctrl+C is hit, then the nodelimit is set to -1 to make Cbc stop
4865         */
4866        if (stoppedOnGap_) {
4867            messageHandler()->message(CBC_GAP, messages())
4868            << bestObjective_ - bestPossibleObjective_
4869            << dblParam_[CbcAllowableGap]
4870            << dblParam_[CbcAllowableFractionGap]*100.0
4871            << CoinMessageEol ;
4872            secondaryStatus_ = 2;
4873            status_ = 0 ;
4874        } else if (maximumSecondsReached()) {
4875            handler_->message(CBC_MAXTIME, messages_) << CoinMessageEol ;
4876            secondaryStatus_ = 4;
4877            status_ = 1 ;
4878        } else if (numberSolutions_ >= intParam_[CbcMaxNumSol]) {
4879            handler_->message(CBC_MAXSOLS, messages_) << CoinMessageEol ;
4880            secondaryStatus_ = 6;
4881            status_ = 1 ;
4882        } else if (isNodeLimitReached()) {
4883            handler_->message(CBC_MAXNODES, messages_) << CoinMessageEol ;
4884            secondaryStatus_ = 3;
4885            status_ = 1 ;
4886        } else if (maximumNumberIterations_ >= 0 && numberIterations_ >= maximumNumberIterations_) {
4887            handler_->message(CBC_MAXITERS, messages_) << CoinMessageEol ;
4888            secondaryStatus_ = 8;
4889            status_ = 1 ;
4890        } else {
4891            handler_->message(CBC_EVENT, messages_) << CoinMessageEol ;
4892            secondaryStatus_ = 5;
4893            status_ = 5 ;
4894        }
4895    }
4896#ifdef CBC_THREAD
4897    if (master_) {
4898        delete master_;
4899        master_ = NULL;
4900        masterThread_ = NULL;
4901    }
4902#endif
4903    /*
4904      That's it, we've exhausted the search tree, or broken out of the loop because
4905      we hit some limit on evaluation.
4906
4907      We may have got an intelligent tree so give it one more chance
4908    */
4909    // Tell solver we are not in Branch and Cut
4910    solver_->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo, NULL) ;
4911    tree_->endSearch();
4912    //  If we did any sub trees - did we give up on any?
4913    if ( numberStoppedSubTrees_)
4914        status_ = 1;
4915    numberNodes_ += numberExtraNodes_;
4916    numberIterations_ += numberExtraIterations_;
4917    if (eventHandler) {
4918        eventHandler->event(CbcEventHandler::endSearch);
4919    }
4920    if (!status_) {
4921        // Set best possible unless stopped on gap
4922        if (secondaryStatus_ != 2)
4923            bestPossibleObjective_ = bestObjective_;
4924        handler_->message(CBC_END_GOOD, messages_)
4925        << bestObjective_ << numberIterations_ << numberNodes_ << getCurrentSeconds()
4926        << CoinMessageEol ;
4927    } else {
4928        handler_->message(CBC_END, messages_)
4929        << bestObjective_ << bestPossibleObjective_
4930        << numberIterations_ << numberNodes_ << getCurrentSeconds()
4931        << CoinMessageEol ;
4932    }
4933    if ((moreSpecialOptions_&4194304)!=0) {
4934      // Conflict cuts
4935      int numberCuts = globalCuts_.sizeRowCuts();
4936      int nConflict=0;
4937      double sizeConflict = 0.0;
4938      for (int i=0;i<numberCuts;i++) {
4939        OsiRowCut2 * cut = globalCuts_.cut(i);
4940        if (cut->whichRow()==1) {
4941          nConflict++;
4942          sizeConflict += cut->row().getNumElements();
4943        }
4944      }
4945      if (nConflict) {
4946        sizeConflict /= nConflict;
4947        char general[200];
4948        sprintf(general, "%d conflict cuts generated - average length %g",
4949                nConflict,sizeConflict);
4950        messageHandler()->message(CBC_GENERAL,
4951                                  messages())
4952          << general << CoinMessageEol ;
4953      }
4954    }
4955    if (numberStrongIterations_)
4956        handler_->message(CBC_STRONG_STATS, messages_)
4957        << strongInfo_[0] << numberStrongIterations_ << strongInfo_[2]
4958        << strongInfo_[1] << CoinMessageEol ;
4959    if (!numberExtraNodes_)
4960        handler_->message(CBC_OTHER_STATS, messages_)
4961        << maximumDepthActual_
4962        << numberDJFixed_ << CoinMessageEol ;
4963    else
4964        handler_->message(CBC_OTHER_STATS2, messages_)
4965        << maximumDepthActual_
4966        << numberDJFixed_ << numberFathoms_ << numberExtraNodes_ << numberExtraIterations_
4967        << CoinMessageEol ;
4968#ifdef COIN_HAS_NTY
4969    if (symmetryInfo_) 
4970      symmetryInfo_->statsOrbits(this,1);
4971#endif
4972    if (doStatistics == 100) {
4973        for (int i = 0; i < numberObjects_; i++) {
4974            CbcSimpleIntegerDynamicPseudoCost * obj =
4975                dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object_[i]) ;
4976            if (obj)
4977                obj->print();
4978        }
4979    }
4980    if (statistics_) {
4981        // report in some way
4982        int * lookup = new int[numberObjects_];
4983        int i;
4984        for (i = 0; i < numberObjects_; i++)
4985            lookup[i] = -1;
4986        bool goodIds = false; //true;
4987        for (i = 0; i < numberObjects_; i++) {
4988            int iColumn = object_[i]->columnNumber();
4989            if (iColumn >= 0 && iColumn < numberColumns) {
4990                if (lookup[i] == -1) {
4991                    lookup[i] = iColumn;
4992                } else {
4993                    goodIds = false;
4994                    break;
4995                }
4996            } else {
4997                goodIds = false;
4998                break;
4999            }
5000        }
5001        if (!goodIds) {
5002            delete [] lookup;
5003            lookup = NULL;
5004        }
5005        if (doStatistics >= 3) {
5006            printf("  node parent depth column   value                    obj      inf\n");
5007            for ( i = 0; i < numberNodes2_; i++) {
5008                statistics_[i]->print(lookup);
5009            }
5010        }
5011        if (doStatistics > 1) {
5012            // Find last solution
5013            int k;
5014            for (k = numberNodes2_ - 1; k >= 0; k--) {
5015                if (statistics_[k]->endingObjective() != COIN_DBL_MAX &&
5016                        !statistics_[k]->endingInfeasibility())
5017                    break;
5018            }
5019            if (k >= 0) {
5020                int depth = statistics_[k]->depth();
5021                int * which = new int[depth+1];
5022                for (i = depth; i >= 0; i--) {
5023                    which[i] = k;
5024                    k = statistics_[k]->parentNode();
5025                }
5026                printf("  node parent depth column   value                    obj      inf\n");
5027                for (i = 0; i <= depth; i++) {
5028                    statistics_[which[i]]->print(lookup);
5029                }
5030                delete [] which;
5031            }
5032        }
5033        // now summary
5034        int maxDepth = 0;
5035        double averageSolutionDepth = 0.0;
5036        int numberSolutions = 0;
5037        double averageCutoffDepth = 0.0;
5038        double averageSolvedDepth = 0.0;
5039        int numberCutoff = 0;
5040        int numberDown = 0;
5041        int numberFirstDown = 0;
5042        double averageInfDown = 0.0;
5043        double averageObjDown = 0.0;
5044        int numberCutoffDown = 0;
5045        int numberUp = 0;
5046        int numberFirstUp = 0;
5047        double averageInfUp = 0.0;
5048        double averageObjUp = 0.0;
5049        int numberCutoffUp = 0;
5050        double averageNumberIterations1 = 0.0;
5051        double averageValue = 0.0;
5052        for ( i = 0; i < numberNodes2_; i++) {
5053            int depth =  statistics_[i]->depth();
5054            int way =  statistics_[i]->way();
5055            double value = statistics_[i]->value();
5056            double startingObjective =  statistics_[i]->startingObjective();
5057            int startingInfeasibility = statistics_[i]->startingInfeasibility();
5058            double endingObjective = statistics_[i]->endingObjective();
5059            int endingInfeasibility = statistics_[i]->endingInfeasibility();
5060            maxDepth = CoinMax(depth, maxDepth);
5061            // Only for completed
5062            averageNumberIterations1 += statistics_[i]->numberIterations();
5063            averageValue += value;
5064            if (endingObjective != COIN_DBL_MAX && !endingInfeasibility) {
5065                numberSolutions++;
5066                averageSolutionDepth += depth;
5067            }
5068            if (endingObjective == COIN_DBL_MAX) {
5069                numberCutoff++;
5070                averageCutoffDepth += depth;
5071                if (way < 0) {
5072                    numberDown++;
5073                    numberCutoffDown++;
5074                    if (way == -1)
5075                        numberFirstDown++;
5076                } else {
5077                    numberUp++;
5078                    numberCutoffUp++;
5079                    if (way == 1)
5080                        numberFirstUp++;
5081                }
5082            } else {
5083                averageSolvedDepth += depth;
5084                if (way < 0) {
5085                    numberDown++;
5086                    averageInfDown += startingInfeasibility - endingInfeasibility;
5087                    averageObjDown += endingObjective - startingObjective;
5088                    if (way == -1)
5089                        numberFirstDown++;
5090                } else {
5091                    numberUp++;
5092                    averageInfUp += startingInfeasibility - endingInfeasibility;
5093                    averageObjUp += endingObjective - startingObjective;
5094                    if (way == 1)
5095                        numberFirstUp++;
5096                }
5097            }
5098        }
5099        // Now print
5100        if (numberSolutions)
5101            averageSolutionDepth /= static_cast<double> (numberSolutions);
5102        int numberSolved = numberNodes2_ - numberCutoff;
5103        double averageNumberIterations2 = numberIterations_ - averageNumberIterations1
5104                                          - numberIterationsAtContinuous;
5105        if (numberCutoff) {
5106            averageCutoffDepth /= static_cast<double> (numberCutoff);
5107            averageNumberIterations2 /= static_cast<double> (numberCutoff);
5108        }
5109        if (numberNodes2_)
5110            averageValue /= static_cast<double> (numberNodes2_);
5111        if (numberSolved) {
5112            averageNumberIterations1 /= static_cast<double> (numberSolved);
5113            averageSolvedDepth /= static_cast<double> (numberSolved);
5114        }
5115        printf("%d solution(s) were found (by branching) at an average depth of %g\n",
5116               numberSolutions, averageSolutionDepth);
5117        printf("average value of variable being branched on was %g\n",
5118               averageValue);
5119        printf("%d nodes were cutoff at an average depth of %g with iteration count of %g\n",
5120               numberCutoff, averageCutoffDepth, averageNumberIterations2);
5121        printf("%d nodes were solved at an average depth of %g with iteration count of %g\n",
5122               numberSolved, averageSolvedDepth, averageNumberIterations1);
5123        if (numberDown) {
5124            averageInfDown /= static_cast<double> (numberDown);
5125            averageObjDown /= static_cast<double> (numberDown);
5126        }
5127        printf("Down %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
5128               numberDown, numberFirstDown, numberDown - numberFirstDown, numberCutoffDown,
5129               averageInfDown, averageObjDown);
5130        if (numberUp) {
5131            averageInfUp /= static_cast<double> (numberUp);
5132            averageObjUp /= static_cast<double> (numberUp);
5133        }
5134        printf("Up %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
5135               numberUp, numberFirstUp, numberUp - numberFirstUp, numberCutoffUp,
5136               averageInfUp, averageObjUp);
5137        for ( i = 0; i < numberNodes2_; i++)
5138            delete statistics_[i];
5139        delete [] statistics_;
5140        statistics_ = NULL;
5141        maximumStatistics_ = 0;
5142        delete [] lookup;
5143    }
5144    /*
5145      If we think we have a solution, restore and confirm it with a call to
5146      setBestSolution().  We need to reset the cutoff value so as not to fathom
5147      the solution on bounds.  Note that calling setBestSolution( ..., true)
5148      leaves the continuousSolver_ bounds vectors fixed at the solution value.
5149
5150      Running resolve() here is a failsafe --- setBestSolution has already
5151      reoptimised using the continuousSolver_. If for some reason we fail to
5152      prove optimality, run the problem again after instructing the solver to
5153      tell us more.
5154
5155      If all looks good, replace solver_ with continuousSolver_, so that the
5156      outside world will be able to obtain information about the solution using
5157      public methods.
5158
5159      Don't replace if we are trying to save cuts
5160    */
5161    if (bestSolution_ && (solverCharacteristics_->solverType() < 2 || solverCharacteristics_->solverType() == 4) && 
5162        ((specialOptions_&8388608)==0||(specialOptions_&2048)!=0)) {
5163        setCutoff(1.0e50) ; // As best solution should be worse than cutoff
5164        // change cutoff as constraint if wanted
5165        if (cutoffRowNumber_>=0) {
5166          if (solver_->getNumRows()>cutoffRowNumber_)
5167            solver_->setRowUpper(cutoffRowNumber_,1.0e50);
5168        }
5169        // also in continuousSolver_
5170        if (continuousSolver_) {
5171          // Solvers know about direction
5172          double direction = solver_->getObjSense();
5173          continuousSolver_->setDblParam(OsiDualObjectiveLimit, 1.0e50*direction);
5174        }
5175        phase_ = 5;
5176        double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
5177        if ((specialOptions_&4) == 0)
5178            bestObjective_ += 100.0 * increment + 1.0e-3; // only set if we are going to solve
5179        setBestSolution(CBC_END_SOLUTION, bestObjective_, bestSolution_, 1) ;
5180        continuousSolver_->resolve() ;
5181        if (!continuousSolver_->isProvenOptimal()) {
5182            continuousSolver_->messageHandler()->setLogLevel(2) ;
5183            continuousSolver_->initialSolve() ;
5184        }
5185        delete solver_ ;
5186        // above deletes solverCharacteristics_
5187        solverCharacteristics_ = NULL;
5188        solver_ = continuousSolver_ ;
5189        setPointers(solver_);
5190        continuousSolver_ = NULL ;
5191    }
5192    /*
5193      Clean up dangling objects. continuousSolver_ may already be toast.
5194    */
5195    delete lastws ;
5196    if (saveObjects) {
5197        for (int i = 0; i < numberObjects_; i++)
5198            delete saveObjects[i];
5199        delete [] saveObjects;
5200    }
5201    numberStrong_ = saveNumberStrong;
5202    numberBeforeTrust_ = saveNumberBeforeTrust;
5203    delete [] whichGenerator_ ;
5204    whichGenerator_ = NULL;
5205    delete [] lowerBefore