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

Last change on this file since 2371 was 2371, checked in by forrest, 3 years ago

finger trouble - out new multiple

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