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

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

changes for rays

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 773.2 KB
Line 
1/* $Id: CbcModel.cpp 2370 2018-03-31 16:45: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      char general[200];
2758      /*
2759        1 - bit set - first come first served
2760        2 - bit set - pass on cutoffs (not deterministic if 1 set)
2761       */
2762#define MULTIPLE_TYPE 3
2763#if MULTIPLE_TYPE == 1 || MULTIPLE_TYPE == 3
2764#define FIRST_COME
2765#endif
2766#if MULTIPLE_TYPE == 2 || MULTIPLE_TYPE == 3
2767#define SET_CUTOFF
2768#endif
2769#ifndef MAX_HEURISTIC_TIMES
2770#define MAX_HEURISTIC_TIMES 16
2771#endif
2772#ifdef SET_CUTOFF
2773      double bestCutoff = getCutoff();
2774#endif
2775#ifdef CBC_THREAD
2776      if (numberRootThreads==1) {
2777#endif
2778        for (int iModel=0;iModel<numberModels;iModel++) {
2779          rootBundle bundle;
2780          bundle.useCutoff=COIN_DBL_MAX;
2781          bundle.model=rootModels[iModel];
2782          bundle.switches=0;
2783#ifdef SET_CUTOFF
2784          bundle.useCutoff=bestCutoff;
2785          rootModels[iModel]->setCutoff(bestCutoff);
2786#endif
2787          doRootCbcThread(&bundle);
2788          sprintf(general,"Ending multiple root solver %d out of %d",
2789                  iModel+1,numberModels);
2790          messageHandler()->message(CBC_GENERAL,messages())
2791            << general << CoinMessageEol ;
2792#ifdef SET_CUTOFF
2793          if (rootModels[iModel]->getCutoff()<bestCutoff)
2794            bestCutoff = rootModels[iModel]->getCutoff();
2795#endif
2796          // see if solved at root node
2797          if (rootModels[iModel]->getMaximumNodes()) {
2798            feasible=false;
2799            break;
2800          }
2801        }
2802#ifdef CBC_THREAD
2803      } else {
2804        // redone so can use numberThreads_ threads more efficiently
2805        // but then if changing cutoff - not deterministic
2806        // BUT changing cutoff has a bug (what bug - debug)
2807        Coin_pthread_t * threadId = new Coin_pthread_t [numberRootThreads];
2808        rootBundle * bundle = new rootBundle [numberRootThreads];
2809        int numberDone=0;
2810        int kModel=numberRootThreads;
2811        for (int iModel=0;iModel<kModel;iModel++) {
2812          threadId[iModel].status=1000+iModel;
2813          bundle[iModel].useCutoff=COIN_DBL_MAX;
2814          if (numberDone>=MAX_HEURISTIC_TIMES) {
2815#ifdef SET_CUTOFF
2816            bundle[iModel].useCutoff=bestCutoff;
2817            rootModels[iModel]->setCutoff(bestCutoff);
2818#endif
2819            for (int i=0;i<rootModels[iModel]->numberHeuristics_;i++)
2820              delete rootModels[iModel]->heuristic_[i];
2821            delete [] rootModels[iModel]->heuristic_;
2822            rootModels[iModel]->heuristic_=NULL;
2823            numberHeuristics_=0;
2824          }
2825          bundle[iModel].model=rootModels[iModel];
2826          bundle[iModel].switches=0;
2827          pthread_create(&(threadId[iModel].thr), NULL, 
2828                         doRootCbcThread,
2829                         bundle+iModel);
2830        }
2831        bool finished=false;
2832        while (numberDone<numberModels) {
2833#ifdef FIRST_COME
2834          usleep(100000);
2835#else
2836          for (int iModel=0;iModel<numberRootThreads;iModel++) {
2837            pthread_join(threadId[iModel].thr, NULL);
2838#ifdef SET_CUTOFF
2839            int jModel = static_cast<int>(threadId[iModel].status-1000);
2840            if (rootModels[jModel]->getCutoff()<bestCutoff)
2841              bestCutoff = rootModels[jModel]->getCutoff();
2842#endif
2843          }
2844#endif
2845          for (int iModel=0;iModel<numberRootThreads;iModel++) {
2846            if (threadId[iModel].status>0) {
2847              int jModel = static_cast<int>(threadId[iModel].status-1000);
2848              if (rootModels[jModel]->secondaryStatus()==123456789) {
2849#ifdef FIRST_COME
2850                pthread_join(threadId[iModel].thr, NULL);
2851#ifdef SET_CUTOFF
2852                if (rootModels[jModel]->getCutoff()<bestCutoff)
2853                  bestCutoff = rootModels[jModel]->getCutoff();
2854#endif
2855#endif
2856                sprintf(general,"Ending multiple root solver %d (%d) out of %d",
2857                        numberDone+1,jModel+1,numberModels);
2858                messageHandler()->message(CBC_GENERAL,messages())
2859                  << general << CoinMessageEol ;
2860                if (rootModels[jModel]->getMaximumNodes()) {
2861                  finished=true;
2862                  numberDone += numberModels-kModel;
2863                  kModel=numberModels;
2864                }
2865                numberDone++;
2866                if (kModel<numberModels) {
2867                  bundle[iModel].useCutoff=COIN_DBL_MAX;
2868                  bundle[iModel].model=rootModels[kModel];
2869                  bundle[iModel].switches=0;
2870#ifdef SET_CUTOFF
2871                  bundle[iModel].useCutoff=bestCutoff;
2872#endif
2873                  if (numberDone>=MAX_HEURISTIC_TIMES) {
2874                    for (int i=0;i<rootModels[kModel]->numberHeuristics_;i++)
2875                      delete rootModels[kModel]->heuristic_[i];
2876                    delete [] rootModels[kModel]->heuristic_;
2877                    rootModels[kModel]->heuristic_=NULL;
2878                    rootModels[kModel]->numberHeuristics_=0;
2879                  }
2880                  threadId[iModel].status=1000+kModel;
2881                  pthread_create(&(threadId[iModel].thr), NULL, 
2882                                 doRootCbcThread,
2883                                 bundle+iModel);
2884                  kModel++;
2885                } else {
2886                  threadId[iModel].status=-1;
2887                }
2888              }
2889            }
2890          }
2891        }
2892        // see if solved at root node
2893        if (finished) 
2894          feasible=false;
2895        delete [] threadId;
2896        delete [] bundle;
2897      }
2898#endif
2899      // sort solutions
2900      int * which = new int [numberModels];
2901      double * value = new double [numberModels];
2902      int numberSolutions=0;
2903      for (int iModel=0;iModel<numberModels;iModel++) {
2904        if (rootModels[iModel]->bestSolution()) {
2905          which[numberSolutions]=iModel;
2906          value[numberSolutions++]=
2907            -rootModels[iModel]->getMinimizationObjValue();
2908        }
2909      }
2910      rootTimeCpu=CoinCpuTime()-rootTimeCpu;
2911      if (numberRootThreads==1)
2912        sprintf(general,"Multiple root solvers took a total of %.2f seconds\n",
2913                rootTimeCpu);
2914      else
2915        sprintf(general,"Multiple root solvers took a total of %.2f seconds (%.2f elapsed)\n",
2916                rootTimeCpu,CoinGetTimeOfDay()-startTimeRoot);
2917      messageHandler()->message(CBC_GENERAL,
2918                                messages())
2919        << general << CoinMessageEol ;
2920      CoinSort_2(value,value+numberSolutions,which);
2921      // to get name
2922      CbcHeuristicRINS dummyHeuristic;
2923      dummyHeuristic.setHeuristicName("Multiple root solvers");
2924      lastHeuristic_=&dummyHeuristic;
2925      for (int i=0;i<numberSolutions;i++) {
2926        double objValue = - value[i];
2927        if (objValue<getCutoff()) {
2928          int iModel=which[i];
2929          setBestSolution(CBC_ROUNDING,objValue,
2930                          rootModels[iModel]->bestSolution());
2931        }
2932      }
2933      lastHeuristic_=NULL;
2934      delete [] which;
2935      delete [] value;
2936    }
2937    // Do heuristics
2938    if (numberObjects_&&!rootModels)
2939        doHeuristicsAtRoot();
2940    if (solverCharacteristics_->solutionAddsCuts()) {
2941      // With some heuristics solver needs a resolve here
2942      solver_->resolve();
2943      if(!isProvenOptimal()){
2944        solver_->initialSolve();
2945      }
2946    }
2947    /*
2948      Grepping through the code, it would appear that this is a command line
2949      debugging hook.  There's no obvious place in the code where this is set to
2950      a negative value.
2951
2952      User hook, says John.
2953    */
2954    if ( intParam_[CbcMaxNumNode] < 0
2955      ||numberSolutions_>=getMaximumSolutions())
2956      eventHappened_ = true; // stop as fast as possible
2957    stoppedOnGap_ = false ;
2958    // See if can stop on gap
2959    bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
2960    if(canStopOnGap()) {
2961        if (bestPossibleObjective_ < getCutoff())
2962            stoppedOnGap_ = true ;
2963        feasible = false;
2964        //eventHappened_=true; // stop as fast as possible
2965    }
2966    /*
2967      Set up for statistics collection, if requested. Standard values are
2968      documented in CbcModel.hpp. The magic number 100 will trigger a dump of
2969      CbcSimpleIntegerDynamicPseudoCost objects (no others). Looks like another
2970      command line debugging hook.
2971    */
2972    statistics_ = NULL;
2973    // Do on switch
2974    if (doStatistics > 0 && doStatistics <= 100) {
2975        maximumStatistics_ = 10000;
2976        statistics_ = new CbcStatistics * [maximumStatistics_];
2977        memset(statistics_, 0, maximumStatistics_*sizeof(CbcStatistics *));
2978    }
2979    // See if we can add integers
2980    if (noObjects && numberIntegers_ < solver_->getNumCols() && (specialOptions_&65536) != 0 && !parentModel_ && false) {
2981      int numberIntegers1=0;
2982      int numberColumns = solver_->getNumCols();
2983      for (int i=0;i<numberColumns;i++) {
2984        if (solver_->isInteger(i))
2985          numberIntegers1++;
2986      }
2987      AddIntegers();
2988      // make sure in sync
2989      int numberIntegers2=0;
2990      for (int i=0;i<numberColumns;i++) {
2991        if (solver_->isInteger(i))
2992          numberIntegers2++;
2993      }
2994      if (numberIntegers1<numberIntegers2) {
2995        findIntegers(true,2);
2996        convertToDynamic();
2997      }
2998    }
2999
3000    /*
3001      Do an initial round of cut generation for the root node. Depending on the
3002      type of underlying solver, we may want to do this even if the initial query
3003      to the objects indicates they're satisfied.
3004
3005      solveWithCuts does the heavy lifting. It will iterate a generate/reoptimise
3006      loop (including reduced cost fixing) until no cuts are generated, the
3007      change in objective falls off,  or the limit on the number of rounds of cut
3008      generation is exceeded.
3009
3010      At the end of all this, any cuts will be recorded in cuts and also
3011      installed in the solver's constraint system. We'll have reoptimised, and
3012      removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
3013      adjusted accordingly).
3014    */
3015    int iObject ;
3016    int numberUnsatisfied = 0 ;
3017    delete [] currentSolution_;
3018    currentSolution_ = new double [numberColumns];
3019    testSolution_ = currentSolution_;
3020    memcpy(currentSolution_, solver_->getColSolution(),
3021           numberColumns*sizeof(double)) ;
3022    // point to useful information
3023    OsiBranchingInformation usefulInfo = usefulInformation();
3024
3025    for (iObject = 0 ; iObject < numberObjects_ ; iObject++) {
3026        double infeasibility =
3027            object_[iObject]->checkInfeasibility(&usefulInfo) ;
3028        if (infeasibility ) numberUnsatisfied++ ;
3029    }
3030    // replace solverType
3031    double * tightBounds = NULL;
3032    if (solverCharacteristics_->tryCuts())  {
3033
3034        if (numberUnsatisfied)   {
3035            // User event
3036            if (!eventHappened_ && feasible) {
3037                if (rootModels) {
3038                  // for fixings
3039                  int numberColumns=solver_->getNumCols();
3040                  tightBounds = new double [2*numberColumns];
3041                  {
3042                    const double * lower = solver_->getColLower();
3043                    const double * upper = solver_->getColUpper();
3044                    for (int i=0;i<numberColumns;i++) {
3045                      tightBounds[2*i+0]=lower[i];
3046                      tightBounds[2*i+1]=upper[i];
3047                    }
3048                  }
3049                  int numberModels = multipleRootTries_%100;
3050                  const OsiSolverInterface ** solvers = new 
3051                    const OsiSolverInterface * [numberModels];
3052                  int numberRows=continuousSolver_->getNumRows();
3053                  int maxCuts=0;
3054                  for (int i=0;i<numberModels;i++) {
3055                    solvers[i]=rootModels[i]->solver();
3056                    const double * lower = solvers[i]->getColLower();
3057                    const double * upper = solvers[i]->getColUpper();
3058                    for (int j=0;j<numberColumns;j++) {
3059                      tightBounds[2*j+0]=CoinMax(lower[j],tightBounds[2*j+0]);
3060                      tightBounds[2*j+1]=CoinMin(upper[j],tightBounds[2*j+1]);
3061                    }
3062                    int numberRows2=solvers[i]->getNumRows();
3063                    assert (numberRows2>=numberRows);
3064                    maxCuts += numberRows2-numberRows;
3065                    // accumulate statistics
3066                    for (int j=0;j<numberCutGenerators_;j++) {
3067                      generator_[j]->addStatistics(rootModels[i]->cutGenerator(j));
3068                    }
3069                  }
3070                  for (int j=0;j<numberCutGenerators_;j++) {
3071                    generator_[j]->scaleBackStatistics(numberModels);
3072                  }
3073                  //CbcRowCuts rowCut(maxCuts);
3074                  const OsiRowCutDebugger *debugger = NULL;
3075                  if ((specialOptions_&1) != 0) 
3076                    debugger = solver_->getRowCutDebugger() ;
3077                  for (int iModel=0;iModel<numberModels;iModel++) {
3078                    int numberRows2=solvers[iModel]->getNumRows();
3079                    const CoinPackedMatrix * rowCopy = solvers[iModel]->getMatrixByRow();
3080                    const int * rowLength = rowCopy->getVectorLengths(); 
3081                    const double * elements = rowCopy->getElements();
3082                    const int * column = rowCopy->getIndices();
3083                    const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3084                    const double * rowLower = solvers[iModel]->getRowLower();
3085                    const double * rowUpper = solvers[iModel]->getRowUpper();
3086                    for (int iRow=numberRows;iRow<numberRows2;iRow++) {
3087                      OsiRowCut rc;
3088                      rc.setLb(rowLower[iRow]);
3089                      rc.setUb(rowUpper[iRow]);
3090                      CoinBigIndex start = rowStart[iRow];
3091                      rc.setRow(rowLength[iRow],column+start,elements+start,false);
3092                      if (debugger)
3093                        CoinAssert (!debugger->invalidCut(rc));
3094                      globalCuts_.addCutIfNotDuplicate(rc);
3095                    }
3096                    //int cutsAdded=globalCuts_.numberCuts()-numberCuts;
3097                    //numberCuts += cutsAdded;
3098                    //printf("Model %d gave %d cuts (out of %d possible)\n",
3099                    //     iModel,cutsAdded,numberRows2-numberRows);
3100                  }
3101                  // normally replace global cuts
3102                  //if (!globalCuts_.())
3103                  //globalCuts_=rowCutrowCut.addCuts(globalCuts_);
3104                  //rowCut.addCuts(globalCuts_);
3105                  int nTightened=0;
3106                  assert(feasible);
3107                  {
3108                    double tolerance=1.0e-5;
3109                    const double * lower = solver_->getColLower();
3110                    const double * upper = solver_->getColUpper();
3111                    for (int i=0;i<numberColumns;i++) {
3112                      if (tightBounds[2*i+0]>tightBounds[2*i+1]+1.0e-9) {
3113                        feasible=false;
3114                        char general[200];
3115                        sprintf(general,"Solvers give infeasible bounds on %d %g,%g was %g,%g - search finished\n",
3116                               i,tightBounds[2*i+0],tightBounds[2*i+1],lower[i],upper[i]);
3117                        messageHandler()->message(CBC_GENERAL,messages())
3118                          << general << CoinMessageEol ;
3119                        break;
3120                      }
3121                      double oldLower=lower[i];
3122                      double oldUpper=upper[i];
3123                      if (tightBounds[2*i+0]>oldLower+tolerance) {
3124                        nTightened++;
3125                        solver_->setColLower(i,tightBounds[2*i+0]);
3126                      }
3127                      if (tightBounds[2*i+1]<oldUpper-tolerance) {
3128                        nTightened++;
3129                        solver_->setColUpper(i,tightBounds[2*i+1]);
3130                      }
3131                    }
3132                  }
3133                  delete [] tightBounds;
3134                  tightBounds=NULL;
3135                  char printBuffer[200];
3136                  sprintf(printBuffer,"%d solvers added %d different cuts out of pool of %d",
3137                          numberModels,globalCuts_.sizeRowCuts(),maxCuts);
3138                  messageHandler()->message(CBC_GENERAL, messages())
3139                    << printBuffer << CoinMessageEol ;
3140                  if (nTightened) {
3141                    sprintf(printBuffer,"%d bounds were tightened",
3142                          nTightened);
3143                    messageHandler()->message(CBC_GENERAL, messages())
3144                      << printBuffer << CoinMessageEol ;
3145                  }
3146                  delete [] solvers;
3147                }
3148                if (!parentModel_&&(moreSpecialOptions_&67108864) != 0) {
3149                  // load cuts from file
3150                  FILE * fp = fopen("global.cuts","rb");
3151                  if (fp) {
3152                    size_t nRead;
3153                    int numberColumns=solver_->getNumCols();
3154                    int numCols;
3155                    nRead = fread(&numCols, sizeof(int), 1, fp);
3156                    if (nRead != 1)
3157                      throw("Error in fread");
3158                    if (numberColumns!=numCols) {
3159                      printf("Mismatch on columns %d %d\n",numberColumns,numCols);
3160                      fclose(fp);
3161                    } else {
3162                      // If rootModel just do some
3163                      double threshold=-1.0;
3164                      if (!multipleRootTries_)
3165                        threshold=0.5;
3166                      int initialCuts=0;
3167                      int initialGlobal = globalCuts_.sizeRowCuts();
3168                      double * elements = new double [numberColumns+2];
3169                      int * indices = new int [numberColumns];
3170                      int numberEntries=1;
3171                      while (numberEntries>0) {
3172                        nRead = fread(&numberEntries, sizeof(int), 1, fp);
3173                        if (nRead != 1)
3174                          throw("Error in fread");
3175                        double randomNumber=randomNumberGenerator_.randomDouble();
3176                        if (numberEntries>0) {
3177                          initialCuts++;
3178                          nRead = fread(elements, sizeof(double), numberEntries+2, fp);
3179                          if (nRead != static_cast<size_t>(numberEntries+2))
3180                            throw("Error in fread");
3181                          nRead = fread(indices, sizeof(int), numberEntries, fp);
3182                          if (nRead != static_cast<size_t>(numberEntries))
3183                            throw("Error in fread");
3184                          if (randomNumber>threshold) {
3185                            OsiRowCut rc;
3186                            rc.setLb(elements[numberEntries]);
3187                            rc.setUb(elements[numberEntries+1]);
3188                            rc.setRow(numberEntries,indices,elements,
3189                                      false);
3190                            rc.setGloballyValidAsInteger(2);
3191                            globalCuts_.addCutIfNotDuplicate(rc) ;
3192                          } 
3193                        }
3194                      }
3195                      fclose(fp);
3196                      // fixes
3197                      int nTightened=0;
3198                      fp = fopen("global.fix","rb");
3199                      if (fp) {
3200                        nRead = fread(indices, sizeof(int), 2, fp);
3201                        if (nRead != 2)
3202                          throw("Error in fread");
3203                        if (numberColumns!=indices[0]) {
3204                          printf("Mismatch on columns %d %d\n",numberColumns,
3205                                 indices[0]);
3206                        } else {
3207                          indices[0]=1;
3208                          while (indices[0]>=0) {
3209                            nRead = fread(indices, sizeof(int), 2, fp);
3210                            if (nRead != 2)
3211                              throw("Error in fread");
3212                            int iColumn=indices[0];
3213                            if (iColumn>=0) {
3214                              nTightened++;
3215                              nRead = fread(elements, sizeof(double), 4, fp);
3216                              if (nRead != 4)
3217                                throw("Error in fread");
3218                              solver_->setColLower(iColumn,elements[0]);
3219                              solver_->setColUpper(iColumn,elements[1]);
3220                            } 
3221                          }
3222                        }
3223                      }
3224                      if (fp)
3225                        fclose(fp);
3226                      char printBuffer[200];
3227                      sprintf(printBuffer,"%d cuts read in of which %d were unique, %d bounds tightened",
3228                             initialCuts,
3229                             globalCuts_.sizeRowCuts()-initialGlobal,nTightened); 
3230                      messageHandler()->message(CBC_GENERAL, messages())
3231                        << printBuffer << CoinMessageEol ;
3232                      delete [] elements;
3233                      delete [] indices;
3234                    }
3235                  }
3236                }
3237                if (feasible)
3238                  feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
3239                                         NULL);
3240                if (multipleRootTries_&&
3241                    (moreSpecialOptions_&134217728)!=0) {
3242                  FILE * fp=NULL;
3243                  size_t nRead;
3244                  int numberColumns=solver_->getNumCols();
3245                  int initialCuts=0;
3246                  if ((moreSpecialOptions_&134217728)!=0) {
3247                    // append so go down to end
3248                    fp = fopen("global.cuts","r+b");
3249                    if (fp) {
3250                      int numCols;
3251                      nRead = fread(&numCols, sizeof(int), 1, fp);
3252                      if (nRead != 1)
3253                        throw("Error in fread");
3254                      if (numberColumns!=numCols) {
3255                        printf("Mismatch on columns %d %d\n",numberColumns,numCols);
3256                        fclose(fp);
3257                        fp=NULL;
3258                      }
3259                    }
3260                  }
3261                  double * elements = new double [numberColumns+2];
3262                  int * indices = new int [numberColumns];
3263                  if (fp) {
3264                    int numberEntries=1;
3265                    while (numberEntries>0) {
3266                      fpos_t position;
3267                      fgetpos(fp, &position);
3268                      nRead = fread(&numberEntries, sizeof(int), 1, fp);
3269                      if (nRead != 1)
3270                        throw("Error in fread");
3271                      if (numberEntries>0) {
3272                        initialCuts++;
3273                        nRead = fread(elements, sizeof(double), numberEntries+2, fp);
3274                        if (nRead != static_cast<size_t>(numberEntries+2))
3275                          throw("Error in fread");
3276                        nRead = fread(indices, sizeof(int), numberEntries, fp);
3277                        if (nRead != static_cast<size_t>(numberEntries))
3278                          throw("Error in fread");
3279                      } else {
3280                        // end
3281                        fsetpos(fp, &position);
3282                      }
3283                    }
3284                  } else {
3285                    fp = fopen("global.cuts","wb");
3286                    size_t nWrite;
3287                    nWrite=fwrite(&numberColumns,sizeof(int),1,fp);
3288                    if (nWrite != 1)
3289                      throw("Error in fwrite");
3290                  }
3291                  size_t nWrite;
3292                  // now append binding cuts
3293                  int numberC=continuousSolver_->getNumRows();
3294                  int numberRows=solver_->getNumRows();
3295                  printf("Saving %d cuts (up from %d)\n",
3296                         initialCuts+numberRows-numberC,initialCuts);
3297                  const double * rowLower = solver_->getRowLower();
3298                  const double * rowUpper = solver_->getRowUpper();
3299                  // Row copy
3300                  CoinPackedMatrix matrixByRow(*solver_->getMatrixByRow());
3301                  const double * elementByRow = matrixByRow.getElements();
3302                  const int * column = matrixByRow.getIndices();
3303                  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
3304                  const int * rowLength = matrixByRow.getVectorLengths();
3305                  for (int iRow=numberC;iRow<numberRows;iRow++) {
3306                    int n=rowLength[iRow];
3307                    assert (n);
3308                    CoinBigIndex start=rowStart[iRow];
3309                    memcpy(elements,elementByRow+start,n*sizeof(double));
3310                    memcpy(indices,column+start,n*sizeof(int));
3311                    elements[n]=rowLower[iRow];
3312                    elements[n+1]=rowUpper[iRow];
3313                    nWrite=fwrite(&n,sizeof(int),1,fp);
3314                    if (nWrite != 1)
3315                      throw("Error in fwrite");
3316                    nWrite=fwrite(elements,sizeof(double),n+2,fp);
3317                    if (nWrite != static_cast<size_t>(n+2))
3318                      throw("Error in fwrite");
3319                    nWrite=fwrite(indices,sizeof(int),n,fp);
3320                    if (nWrite != static_cast<size_t>(n))
3321                      throw("Error in fwrite");
3322                  }
3323                  // eof marker
3324                  int eofMarker=-1;
3325                  nWrite=fwrite(&eofMarker,sizeof(int),1,fp);
3326                  if (nWrite != 1)
3327                    throw("Error in fwrite");
3328                  fclose(fp);
3329                  // do tighter bounds (? later extra to original columns)
3330                  int nTightened=0;
3331                  const double * lower = solver_->getColLower();
3332                  const double * upper = solver_->getColUpper();
3333                  const double * originalLower = continuousSolver_->getColLower();
3334                  const double * originalUpper = continuousSolver_->getColUpper();
3335                  double tolerance=1.0e-5;
3336                  for (int i=0;i<numberColumns;i++) {
3337                    if (lower[i]>originalLower[i]+tolerance) {
3338                      nTightened++;
3339                    }
3340                    if (upper[i]<originalUpper[i]-tolerance) {
3341                      nTightened++;
3342                    }
3343                  }
3344                  if (nTightened) {
3345                    fp = fopen("global.fix","wb");
3346                    size_t nWrite;
3347                    indices[0]=numberColumns;
3348                    if (originalColumns_)
3349                      indices[1]=COIN_INT_MAX;
3350                    else
3351                      indices[1]=-1;
3352                    nWrite=fwrite(indices,sizeof(int),2,fp);
3353                    if (nWrite != 2)
3354                      throw("Error in fwrite");
3355                    for (int i=0;i<numberColumns;i++) {
3356                      int nTightened=0;
3357                      if (lower[i]>originalLower[i]+tolerance) {
3358                        nTightened++;
3359                      }
3360                      if (upper[i]<originalUpper[i]-tolerance) {
3361                        nTightened++;
3362                      }
3363                      if (nTightened) {
3364                        indices[0]=i;
3365                        if (originalColumns_)
3366                          indices[1]=originalColumns_[i];
3367                        elements[0]=lower[i];
3368                        elements[1]=upper[i];
3369                        elements[2]=originalLower[i];
3370                        elements[3]=originalUpper[i];
3371                        nWrite=fwrite(indices,sizeof(int),2,fp);
3372                        if (nWrite != 2)
3373                          throw("Error in fwrite");
3374                        nWrite=fwrite(elements,sizeof(double),4,fp);
3375                        if (nWrite != 4)
3376                          throw("Error in fwrite");
3377                      }
3378                    }
3379                    // eof marker
3380                    indices[0]=-1;
3381                    nWrite=fwrite(indices,sizeof(int),2,fp);
3382                    if (nWrite != 2)
3383                      throw("Error in fwrite");
3384                    fclose(fp);
3385                  }
3386                  delete [] elements;
3387                  delete [] indices;
3388                }
3389                if ((specialOptions_&524288) != 0 && !parentModel_
3390                        && storedRowCuts_) {
3391                    if (feasible) {
3392                        /* pick up stuff and try again
3393                        add cuts, maybe keep around
3394                        do best solution and if so new heuristics
3395                        obviously tighten bounds
3396                        */
3397                        // A and B probably done on entry
3398                        // A) tight bounds on integer variables
3399                        const double * lower = solver_->getColLower();
3400                        const double * upper = solver_->getColUpper();
3401                        const double * tightLower = storedRowCuts_->tightLower();
3402                        const double * tightUpper = storedRowCuts_->tightUpper();
3403                        int nTightened = 0;
3404                        for (int i = 0; i < numberIntegers_; i++) {
3405                            int iColumn = integerVariable_[i];
3406                            if (tightLower[iColumn] > lower[iColumn]) {
3407                                nTightened++;
3408                                solver_->setColLower(iColumn, tightLower[iColumn]);
3409                            }
3410                            if (tightUpper[iColumn] < upper[iColumn]) {
3411                                nTightened++;
3412                                solver_->setColUpper(iColumn, tightUpper[iColumn]);
3413                            }
3414                        }
3415                        if (nTightened)
3416                          COIN_DETAIL_PRINT(printf("%d tightened by alternate cuts\n", nTightened));
3417                        if (storedRowCuts_->bestObjective() < bestObjective_) {
3418                            // B) best solution
3419                            double objValue = storedRowCuts_->bestObjective();
3420                            setBestSolution(CBC_SOLUTION, objValue,
3421                                            storedRowCuts_->bestSolution()) ;
3422                            // Do heuristics
3423                            // Allow RINS
3424                            for (int i = 0; i < numberHeuristics_; i++) {
3425                                CbcHeuristicRINS * rins
3426                                = dynamic_cast<CbcHeuristicRINS *> (heuristic_[i]);
3427                                if (rins) {
3428                                    rins->setLastNode(-100);
3429                                }
3430                            }
3431                            doHeuristicsAtRoot();
3432                        }
3433#ifdef JJF_ZERO
3434                        int nCuts = storedRowCuts_->sizeRowCuts();
3435                        // add to global list
3436                        for (int i = 0; i < nCuts; i++) {
3437                            OsiRowCut newCut(*storedRowCuts_->rowCutPointer(i));
3438                            newCut.setGloballyValidAsInteger(2);
3439                            newCut.mutableRow().setTestForDuplicateIndex(false);
3440                            globalCuts_.insert(newCut) ;
3441                        }
3442#else
3443                        addCutGenerator(storedRowCuts_, -99, "Stored from previous run",
3444                                        true, false, false, -200);
3445#endif
3446                        // Set cuts as active
3447                        delete [] addedCuts_ ;
3448                        maximumNumberCuts_ = cuts.sizeRowCuts();
3449                        if (maximumNumberCuts_) {
3450                            addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
3451                        } else {
3452                            addedCuts_ = NULL;
3453                        }
3454                        for (int i = 0; i < maximumNumberCuts_; i++)
3455                            addedCuts_[i] = new CbcCountRowCut(*cuts.rowCutPtr(i),
3456                                                               NULL, -1, -1, 2);
3457                        COIN_DETAIL_PRINT(printf("size %d\n", cuts.sizeRowCuts()));
3458                        cuts = OsiCuts();
3459                        currentNumberCuts_ = maximumNumberCuts_;
3460                        feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
3461                                                 NULL);
3462                        for (int i = 0; i < maximumNumberCuts_; i++)
3463                            delete addedCuts_[i];
3464                    }
3465                    delete storedRowCuts_;
3466                    storedRowCuts_ = NULL;
3467                }
3468            } else {
3469                feasible = false;
3470            }
3471        }       else if (solverCharacteristics_->solutionAddsCuts() ||
3472                   solverCharacteristics_->alwaysTryCutsAtRootNode()) {
3473            // may generate cuts and turn the solution
3474            //to an infeasible one
3475            feasible = solveWithCuts(cuts, 2,
3476                                     NULL);
3477        }
3478    }
3479    if (rootModels) {
3480      int numberModels = multipleRootTries_%100;
3481      for (int i=0;i<numberModels;i++)
3482        delete rootModels[i];
3483      delete [] rootModels;
3484    }
3485    // check extra info on feasibility
3486    if (!solverCharacteristics_->mipFeasible())
3487        feasible = false;
3488    // If max nodes==0 - don't do strong branching
3489    if (!getMaximumNodes()) {
3490      if (feasible)
3491        feasible=false;
3492      else
3493        setMaximumNodes(1); //allow to stop on success
3494    }
3495    topOfTree_=NULL;
3496#ifdef CLP_RESOLVE
3497    if ((moreSpecialOptions_&2097152)!=0&&!parentModel_&&feasible) {
3498      OsiClpSolverInterface * clpSolver
3499        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3500      if (clpSolver)
3501        resolveClp(clpSolver,0);
3502    }
3503#endif
3504    // make cut generators less aggressive
3505    for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
3506        CglCutGenerator * generator = generator_[iCutGenerator]->generator();
3507        generator->setAggressiveness(generator->getAggressiveness() - 100);
3508    }
3509    currentNumberCuts_ = numberNewCuts_ ;
3510    if (solverCharacteristics_->solutionAddsCuts()) {
3511      // With some heuristics solver needs a resolve here (don't know if this is bug in heuristics)
3512      solver_->resolve();
3513      if(!isProvenOptimal()){
3514        solver_->initialSolve();
3515      }
3516    }
3517    // See if can stop on gap
3518    bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
3519    if(canStopOnGap()) {
3520        if (bestPossibleObjective_ < getCutoff())
3521            stoppedOnGap_ = true ;
3522        feasible = false;
3523    }
3524    // User event
3525    if (eventHappened_)
3526        feasible = false;
3527#if defined(COIN_HAS_CLP)&&defined(COIN_HAS_CPX)
3528    /*
3529      This is the notion of using Cbc stuff to get going, then calling cplex to
3530      finish off.
3531    */
3532    if (feasible && (specialOptions_&16384) != 0 && fastNodeDepth_ == -2 && !parentModel_) {
3533        // Use Cplex to do search!
3534        double time1 = CoinCpuTime();
3535        OsiClpSolverInterface * clpSolver
3536        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3537        OsiCpxSolverInterface cpxSolver;
3538        double direction = clpSolver->getObjSense();
3539        cpxSolver.setObjSense(direction);
3540        // load up cplex
3541        const CoinPackedMatrix * matrix = continuousSolver_->getMatrixByCol();
3542        const double * rowLower = continuousSolver_->getRowLower();
3543        const double * rowUpper = continuousSolver_->getRowUpper();
3544        const double * columnLower = continuousSolver_->getColLower();
3545        const double * columnUpper = continuousSolver_->getColUpper();
3546        const double * objective = continuousSolver_->getObjCoefficients();
3547        cpxSolver.loadProblem(*matrix, columnLower, columnUpper,
3548                              objective, rowLower, rowUpper);
3549        double * setSol = new double [numberIntegers_];
3550        int * setVar = new int [numberIntegers_];
3551        // cplex doesn't know about objective offset
3552        double offset = clpSolver->getModelPtr()->objectiveOffset();
3553        for (int i = 0; i < numberIntegers_; i++) {
3554            int iColumn = integerVariable_[i];
3555            cpxSolver.setInteger(iColumn);
3556            if (bestSolution_) {
3557                setSol[i] = bestSolution_[iColumn];
3558                setVar[i] = iColumn;
3559            }
3560        }
3561        CPXENVptr env = cpxSolver.getEnvironmentPtr();
3562        CPXLPptr lpPtr = cpxSolver.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL);
3563        cpxSolver.switchToMIP();
3564        if (bestSolution_) {
3565#if 0
3566            CPXcopymipstart(env, lpPtr, numberIntegers_, setVar, setSol);
3567#else
3568            int zero = 0;
3569            CPXaddmipstarts(env, lpPtr, 1, numberIntegers_, &zero, setVar, setSol, NULL, NULL);
3570#endif
3571        }
3572        if (clpSolver->getNumRows() > continuousSolver_->getNumRows() && false) {
3573            // add cuts
3574            const CoinPackedMatrix * matrix = clpSolver->getMatrixByRow();
3575            const double * rhs = clpSolver->getRightHandSide();
3576            const char * rowSense = clpSolver->getRowSense();
3577            const double * elementByRow = matrix->getElements();
3578            const int * column = matrix->getIndices();
3579            const CoinBigIndex * rowStart = matrix->getVectorStarts();
3580            const int * rowLength = matrix->getVectorLengths();
3581            int nStart = continuousSolver_->getNumRows();
3582            int nRows = clpSolver->getNumRows();
3583            int size = rowStart[nRows-1] + rowLength[nRows-1] -
3584                       rowStart[nStart];
3585            int nAdd = 0;
3586            double * rmatval = new double [size];
3587            int * rmatind = new int [size];
3588            int * rmatbeg = new int [nRows-nStart+1];
3589            size = 0;
3590            rmatbeg[0] = 0;
3591            for (int i = nStart; i < nRows; i++) {
3592                for (int k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
3593                    rmatind[size] = column[k];
3594                    rmatval[size++] = elementByRow[k];
3595                }
3596                nAdd++;
3597                rmatbeg[nAdd] = size;
3598            }
3599            CPXaddlazyconstraints(env, lpPtr, nAdd, size,
3600                                  rhs, rowSense, rmatbeg,
3601                                  rmatind, rmatval, NULL);
3602            CPXsetintparam( env, CPX_PARAM_REDUCE,
3603                            // CPX_PREREDUCE_NOPRIMALORDUAL (0)
3604                            CPX_PREREDUCE_PRIMALONLY);
3605        }
3606        if (getCutoff() < 1.0e50) {
3607            double useCutoff = getCutoff() + offset;
3608            if (bestObjective_ < 1.0e50)
3609                useCutoff = bestObjective_ + offset + 1.0e-7;
3610            cpxSolver.setDblParam(OsiDualObjectiveLimit, useCutoff*
3611                                  direction);
3612            if ( direction > 0.0 )
3613                CPXsetdblparam( env, CPX_PARAM_CUTUP, useCutoff ) ; // min
3614            else
3615                CPXsetdblparam( env, CPX_PARAM_CUTLO, useCutoff ) ; // max
3616        }
3617        CPXsetdblparam(env, CPX_PARAM_EPGAP, dblParam_[CbcAllowableFractionGap]);
3618        delete [] setSol;
3619        delete [] setVar;
3620        char printBuffer[200];
3621        if (offset) {
3622            sprintf(printBuffer, "Add %g to all Cplex messages for true objective",
3623                    -offset);
3624            messageHandler()->message(CBC_GENERAL, messages())
3625            << printBuffer << CoinMessageEol ;
3626            cpxSolver.setDblParam(OsiObjOffset, offset);
3627        }
3628        cpxSolver.branchAndBound();
3629        double timeTaken = CoinCpuTime() - time1;
3630        sprintf(printBuffer, "Cplex took %g seconds",
3631                timeTaken);
3632        messageHandler()->message(CBC_GENERAL, messages())
3633        << printBuffer << CoinMessageEol ;
3634        numberExtraNodes_ = CPXgetnodecnt(env, lpPtr);
3635        numberExtraIterations_ = CPXgetmipitcnt(env, lpPtr);
3636        double value = cpxSolver.getObjValue() * direction;
3637        if (cpxSolver.isProvenOptimal() && value <= getCutoff()) {
3638            feasible = true;
3639            clpSolver->setWarmStart(NULL);
3640            // try and do solution
3641            double * newSolution =
3642                CoinCopyOfArray(cpxSolver.getColSolution(),
3643                                getNumCols());
3644            setBestSolution(CBC_STRONGSOL, value, newSolution) ;
3645            delete [] newSolution;
3646        }
3647        feasible = false;
3648    }
3649#endif
3650    if (!parentModel_&&(moreSpecialOptions_&268435456) != 0) {
3651      // try idiotic idea
3652      CbcObject * obj = new CbcIdiotBranch(this);
3653      obj->setPriority(1); // temp
3654      addObjects(1, &obj);
3655      delete obj;
3656    }
3657   
3658    /*
3659      A hook to use clp to quickly explore some part of the tree.
3660    */
3661    if (fastNodeDepth_ == 1000 &&/*!parentModel_*/(specialOptions_&2048) == 0) {
3662        fastNodeDepth_ = -1;
3663        CbcObject * obj =
3664            new CbcFollowOn(this);
3665        obj->setPriority(1);
3666        addObjects(1, &obj);
3667        delete obj;
3668    }
3669    int saveNumberSolves = numberSolves_;
3670    int saveNumberIterations = numberIterations_;
3671    if ((fastNodeDepth_ >= 0||(moreSpecialOptions_&33554432)!=0)
3672        &&/*!parentModel_*/(specialOptions_&2048) == 0) {
3673        // add in a general depth object doClp
3674        int type = (fastNodeDepth_ <= 100) ? fastNodeDepth_ : -(fastNodeDepth_ - 100);
3675        if ((moreSpecialOptions_&33554432)!=0)
3676          type=12;
3677        else
3678          fastNodeDepth_ += 1000000;     // mark as done
3679        CbcObject * obj =
3680            new CbcGeneralDepth(this, type);
3681        addObjects(1, &obj);
3682        delete obj;
3683        // fake number of objects
3684        numberObjects_--;
3685        if (parallelMode() < -1) {
3686            // But make sure position is correct
3687            OsiObject * obj2 = object_[numberObjects_];
3688            obj = dynamic_cast<CbcObject *> (obj2);
3689            assert (obj);
3690            obj->setPosition(numberObjects_);
3691        }
3692    }
3693#ifdef COIN_HAS_CLP
3694#ifdef NO_CRUNCH
3695    if (true) {
3696        OsiClpSolverInterface * clpSolver
3697        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3698        if (clpSolver && !parentModel_) {
3699            ClpSimplex * clpSimplex = clpSolver->getModelPtr();
3700            clpSimplex->setSpecialOptions(clpSimplex->specialOptions() | 131072);
3701            //clpSimplex->startPermanentArrays();
3702            clpSimplex->setPersistenceFlag(2);
3703        }
3704    }
3705#endif
3706#endif
3707// Save copy of solver
3708    OsiSolverInterface * saveSolver = NULL;
3709    if (!parentModel_ && (specialOptions_&(512 + 32768)) != 0)
3710        saveSolver = solver_->clone();
3711    double checkCutoffForRestart = 1.0e100;
3712    saveModel(saveSolver, &checkCutoffForRestart, &feasible);
3713    if ((specialOptions_&262144) != 0 && !parentModel_) {
3714        // Save stuff and return!
3715        storedRowCuts_->saveStuff(bestObjective_, bestSolution_,
3716                                  solver_->getColLower(),
3717                                  solver_->getColUpper());
3718        delete [] lowerBefore;
3719        delete [] upperBefore;
3720        delete saveSolver;
3721        if (flipObjective)
3722          flipModel();
3723        return;
3724    }
3725    /*
3726      We've taken the continuous relaxation as far as we can. Time to branch.
3727      The first order of business is to actually create a node. chooseBranch
3728      currently uses strong branching to evaluate branch object candidates,
3729      unless forced back to simple branching. If chooseBranch concludes that a
3730      branching candidate is monotone (anyAction == -1) or infeasible (anyAction
3731      == -2) when forced to integer values, it returns here immediately.
3732
3733      Monotone variables trigger a call to resolve(). If the problem remains
3734      feasible, try again to choose a branching variable. At the end of the loop,
3735      resolved == true indicates that some variables were fixed.
3736
3737      Loss of feasibility will result in the deletion of newNode.
3738    */
3739
3740    bool resolved = false ;
3741    CbcNode *newNode = NULL ;
3742    numberFixedAtRoot_ = 0;
3743    numberFixedNow_ = 0;
3744    if (!parentModel_&&(moreSpecialOptions2_&2)!=0) {
3745#ifdef COIN_HAS_CLP
3746      OsiClpSolverInterface * clpSolver
3747        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3748      if (clpSolver) {
3749        if (getCutoff()>1.0e20) {
3750          printf("Zapping costs\n");
3751          int numberColumns=solver_->getNumCols();
3752          double * zeroCost = new double [numberColumns];
3753          // could make small random
3754          memset(zeroCost,0,numberColumns*sizeof(double));
3755          solver_->setObjective(zeroCost);
3756          double objValue = solver_->getObjValue();
3757          solver_->setDblParam(OsiObjOffset,-objValue);
3758          clpSolver->getModelPtr()->setObjectiveValue(objValue);
3759          delete [] zeroCost;
3760        } else {
3761          moreSpecialOptions2_ &= ~2;
3762        }
3763      } else {
3764#endif
3765          moreSpecialOptions2_ &= ~2;
3766#ifdef COIN_HAS_CLP
3767      }
3768#endif
3769    }
3770    int numberIterationsAtContinuous = numberIterations_;
3771    //solverCharacteristics_->setSolver(solver_);
3772    if (feasible) {
3773      // mark all cuts as globally valid
3774      int numberCuts=cuts.sizeRowCuts();
3775      resizeWhichGenerator(0,numberCuts);
3776      for (int i=0;i<numberCuts;i++) {
3777        cuts.rowCutPtr(i)->setGloballyValid();
3778        whichGenerator_[i]=20000+(whichGenerator_[i]%10000);
3779      }
3780#define HOTSTART -1
3781#if HOTSTART<0
3782        if (bestSolution_ && !parentModel_ && !hotstartSolution_ &&
3783                (moreSpecialOptions_&1024) != 0) {
3784            // Set priorities so only branch on ones we need to
3785            // use djs and see if only few branches needed
3786#ifndef NDEBUG
3787            double integerTolerance = getIntegerTolerance() ;
3788#endif
3789            bool possible = true;
3790            const double * saveLower = continuousSolver_->getColLower();
3791            const double * saveUpper = continuousSolver_->getColUpper();
3792            for (int i = 0; i < numberObjects_; i++) {
3793                const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object_[i]);
3794                if (thisOne) {
3795                    int iColumn = thisOne->columnNumber();
3796                    if (saveUpper[iColumn] > saveLower[iColumn] + 1.5) {
3797                        possible = false;
3798                        break;
3799                    }
3800                } else {
3801                    possible = false;
3802                    break;
3803                }
3804            }
3805            if (possible) {
3806                OsiSolverInterface * solver = continuousSolver_->clone();
3807                int numberColumns = solver->getNumCols();
3808                for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
3809                    double value = bestSolution_[iColumn] ;
3810                    value = CoinMax(value, saveLower[iColumn]) ;
3811                    value = CoinMin(value, saveUpper[iColumn]) ;
3812                    value = floor(value + 0.5);
3813                    if (solver->isInteger(iColumn)) {
3814                        solver->setColLower(iColumn, value);
3815                        solver->setColUpper(iColumn, value);
3816                    }
3817                }
3818                solver->setHintParam(OsiDoDualInResolve, false, OsiHintTry);
3819                // objlim and all slack
3820                double direction = solver->getObjSense();
3821                solver->setDblParam(OsiDualObjectiveLimit, 1.0e50*direction);
3822                CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis *> (solver->getEmptyWarmStart());
3823                solver->setWarmStart(basis);
3824                delete basis;
3825                bool changed = true;
3826                hotstartPriorities_ = new int [numberColumns];
3827                for (int iColumn = 0; iColumn < numberColumns; iColumn++)
3828                    hotstartPriorities_[iColumn] = 1;
3829                while (changed) {
3830                    changed = false;
3831                    solver->resolve();
3832                    if (!solver->isProvenOptimal()) {
3833                        possible = false;
3834                        break;
3835                    }
3836                    const double * dj = solver->getReducedCost();
3837                    const double * colLower = solver->getColLower();
3838                    const double * colUpper = solver->getColUpper();
3839                    const double * solution = solver->getColSolution();
3840                    int nAtLbNatural = 0;
3841                    int nAtUbNatural = 0;
3842                    int nZeroDj = 0;
3843                    int nForced = 0;
3844                    for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
3845                        double value = solution[iColumn] ;
3846                        value = CoinMax(value, saveLower[iColumn]) ;
3847                        value = CoinMin(value, saveUpper[iColumn]) ;
3848                        if (solver->isInteger(iColumn)) {
3849                            assert(fabs(value - solution[iColumn]) <= integerTolerance) ;
3850                            if (hotstartPriorities_[iColumn] == 1) {
3851                                if (dj[iColumn] < -1.0e-6) {
3852                                    // negative dj
3853                                    if (saveUpper[iColumn] == colUpper[iColumn]) {
3854                                        nAtUbNatural++;
3855                                        hotstartPriorities_[iColumn] = 2;
3856                                        solver->setColLower(iColumn, saveLower[iColumn]);
3857                                        solver->setColUpper(iColumn, saveUpper[iColumn]);
3858                                    } else {
3859                                        nForced++;
3860                                    }
3861                                } else if (dj[iColumn] > 1.0e-6) {
3862                                    // positive dj
3863                                    if (saveLower[iColumn] == colLower[iColumn]) {
3864                                        nAtLbNatural++;
3865                                        hotstartPriorities_[iColumn] = 2;
3866                                        solver->setColLower(iColumn, saveLower[iColumn]);
3867                                        solver->setColUpper(iColumn, saveUpper[iColumn]);
3868                                    } else {
3869                                        nForced++;
3870                                    }
3871                                } else {
3872                                    // zero dj
3873                                    nZeroDj++;
3874                                }
3875                            }
3876                        }
3877                    }
3878#if CBC_USEFUL_PRINTING>1
3879                    printf("%d forced, %d naturally at lower, %d at upper - %d zero dj\n",
3880                           nForced, nAtLbNatural, nAtUbNatural, nZeroDj);
3881#endif
3882                    if (nAtLbNatural || nAtUbNatural) {
3883                        changed = true;
3884                    } else {
3885                        if (nForced + nZeroDj > 5000 ||
3886                                (nForced + nZeroDj)*2 > numberIntegers_)
3887                            possible = false;
3888                    }
3889                }
3890                delete solver;
3891            }
3892            if (possible) {
3893                setHotstartSolution(bestSolution_);
3894                if (!saveCompare) {
3895                    // create depth first comparison
3896                    saveCompare = nodeCompare_;
3897                    // depth first
3898                    nodeCompare_ = new CbcCompareDepth();
3899                    tree_->setComparison(*nodeCompare_) ;
3900                }
3901            } else {
3902                delete [] hotstartPriorities_;
3903                hotstartPriorities_ = NULL;
3904            }
3905        }
3906#endif
3907#if HOTSTART>0
3908        if (hotstartSolution_ && !hotstartPriorities_) {
3909            // Set up hot start
3910            OsiSolverInterface * solver = solver_->clone();
3911            double direction = solver_->getObjSense() ;
3912            int numberColumns = solver->getNumCols();
3913            double * saveLower = CoinCopyOfArray(solver->getColLower(), numberColumns);
3914            double * saveUpper = CoinCopyOfArray(solver->getColUpper(), numberColumns);
3915            // move solution
3916            solver->setColSolution(hotstartSolution_);
3917            // point to useful information
3918            const double * saveSolution = testSolution_;
3919            testSolution_ = solver->getColSolution();
3920            OsiBranchingInformation usefulInfo = usefulInformation();
3921            testSolution_ = saveSolution;
3922            /*
3923            Run through the objects and use feasibleRegion() to set variable bounds
3924            so as to fix the variables specified in the objects at their value in this
3925            solution. Since the object list contains (at least) one object for every
3926            integer variable, this has the effect of fixing all integer variables.
3927            */
3928            for (int i = 0; i < numberObjects_; i++)
3929                object_[i]->feasibleRegion(solver, &usefulInfo);
3930            solver->resolve();
3931            assert (solver->isProvenOptimal());
3932            double gap = CoinMax((solver->getObjValue() - solver_->getObjValue()) * direction, 0.0) ;
3933            const double * dj = solver->getReducedCost();
3934            const double * colLower = solver->getColLower();
3935            const double * colUpper = solver->getColUpper();
3936            const double * solution = solver->getColSolution();
3937            int nAtLbNatural = 0;
3938            int nAtUbNatural = 0;
3939            int nAtLbNaturalZero = 0;
3940            int nAtUbNaturalZero = 0;
3941            int nAtLbFixed = 0;
3942            int nAtUbFixed = 0;
3943            int nAtOther = 0;
3944            int nAtOtherNatural = 0;
3945            int nNotNeeded = 0;
3946            delete [] hotstartSolution_;
3947            hotstartSolution_ = new double [numberColumns];
3948            delete [] hotstartPriorities_;
3949            hotstartPriorities_ = new int [numberColumns];
3950            int * order = (int *) saveUpper;
3951            int nFix = 0;
3952            double bestRatio = COIN_DBL_MAX;
3953            for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
3954                double value = solution[iColumn] ;
3955                value = CoinMax(value, saveLower[iColumn]) ;
3956                value = CoinMin(value, saveUpper[iColumn]) ;
3957                double sortValue = COIN_DBL_MAX;
3958                if (solver->isInteger(iColumn)) {
3959                    assert(fabs(value - solution[iColumn]) <= 1.0e-5) ;
3960                    double value2 = floor(value + 0.5);
3961                    if (dj[iColumn] < -1.0e-6) {
3962                        // negative dj
3963                        //assert (value2==colUpper[iColumn]);
3964                        if (saveUpper[iColumn] == colUpper[iColumn]) {
3965                            nAtUbNatural++;
3966                            sortValue = 0.0;
3967                            double value = -dj[iColumn];
3968                            if (value > gap)
3969                                nFix++;
3970                            else if (gap < value*bestRatio)
3971                                bestRatio = gap / value;
3972                            if (saveLower[iColumn] != colLower[iColumn]) {
3973                                nNotNeeded++;
3974                                sortValue = 1.0e20;
3975                            }
3976                        } else if (saveLower[iColumn] == colUpper[iColumn]) {
3977                            nAtLbFixed++;
3978                            sortValue = dj[iColumn];
3979                        } else {
3980                            nAtOther++;
3981                            sortValue = 0.0;
3982                            if (saveLower[iColumn] != colLower[iColumn] &&
3983                                    saveUpper[iColumn] != colUpper[iColumn]) {
3984                                nNotNeeded++;
3985                                sortValue = 1.0e20;
3986                            }
3987                        }
3988                    } else if (dj[iColumn] > 1.0e-6) {
3989                        // positive dj
3990                        //assert (value2==colLower[iColumn]);
3991                        if (saveLower[iColumn] == colLower[iColumn]) {
3992                            nAtLbNatural++;
3993                            sortValue = 0.0;
3994                            double value = dj[iColumn];
3995                            if (value > gap)
3996                                nFix++;
3997                            else if (gap < value*bestRatio)
3998                                bestRatio = gap / value;
3999                            if (saveUpper[iColumn] != colUpper[iColumn]) {
4000                                nNotNeeded++;
4001                                sortValue = 1.0e20;
4002                            }
4003                        } else if (saveUpper[iColumn] == colLower[iColumn]) {
4004                            nAtUbFixed++;
4005                            sortValue = -dj[iColumn];
4006                        } else {
4007                            nAtOther++;
4008                            sortValue = 0.0;
4009                            if (saveLower[iColumn] != colLower[iColumn] &&
4010                                    saveUpper[iColumn] != colUpper[iColumn]) {
4011                                nNotNeeded++;
4012                                sortValue = 1.0e20;
4013                            }
4014                        }
4015                    } else {
4016                        // zero dj
4017                        if (value2 == saveUpper[iColumn]) {
4018                            nAtUbNaturalZero++;
4019                            sortValue = 0.0;
4020                            if (saveLower[iColumn] != colLower[iColumn]) {
4021                                nNotNeeded++;
4022                                sortValue = 1.0e20;
4023                            }
4024                        } else if (value2 == saveLower[iColumn]) {
4025                            nAtLbNaturalZero++;
4026                            sortValue = 0.0;
4027                        } else {
4028                            nAtOtherNatural++;
4029                            sortValue = 0.0;
4030                            if (saveLower[iColumn] != colLower[iColumn] &&
4031                                    saveUpper[iColumn] != colUpper[iColumn]) {
4032                                nNotNeeded++;
4033                                sortValue = 1.0e20;
4034                            }
4035                        }
4036                    }
4037#if HOTSTART==3
4038                    sortValue = -fabs(dj[iColumn]);
4039#endif
4040                }
4041                hotstartSolution_[iColumn] = value ;
4042                saveLower[iColumn] = sortValue;
4043                order[iColumn] = iColumn;
4044            }
4045            COIN_DETAIL_PRINT(printf("** can fix %d columns - best ratio for others is %g on gap of %g\n",
4046                                     nFix, bestRatio, gap));
4047            int nNeg = 0;
4048            CoinSort_2(saveLower, saveLower + numberColumns, order);
4049            for (int i = 0; i < numberColumns; i++) {
4050                if (saveLower[i] < 0.0) {
4051                    nNeg++;
4052#if HOTSTART==2||HOTSTART==3
4053                    // swap sign ?
4054                    saveLower[i] = -saveLower[i];
4055#endif
4056                }
4057            }
4058            CoinSort_2(saveLower, saveLower + nNeg, order);
4059            for (int i = 0; i < numberColumns; i++) {
4060#if HOTSTART==1
4061                hotstartPriorities_[order[i]] = 100;
4062#else
4063                hotstartPriorities_[order[i]] = -(i + 1);
4064#endif
4065            }
4066            COIN_DETAIL_PRINT(printf("nAtLbNat %d,nAtUbNat %d,nAtLbNatZero %d,nAtUbNatZero %d,nAtLbFixed %d,nAtUbFixed %d,nAtOther %d,nAtOtherNat %d, useless %d %d\n",
4067                   nAtLbNatural,
4068                   nAtUbNatural,
4069                   nAtLbNaturalZero,
4070                   nAtUbNaturalZero,
4071                   nAtLbFixed,
4072                   nAtUbFixed,
4073                   nAtOther,
4074                                     nAtOtherNatural, nNotNeeded, nNeg));
4075            delete [] saveLower;
4076            delete [] saveUpper;
4077            if (!saveCompare) {
4078                // create depth first comparison
4079                saveCompare = nodeCompare_;
4080                // depth first
4081                nodeCompare_ = new CbcCompareDepth();
4082                tree_->setComparison(*nodeCompare_) ;
4083            }
4084        }
4085#endif
4086        newNode = new CbcNode ;
4087        // Set objective value (not so obvious if NLP etc)
4088        setObjectiveValue(newNode, NULL);
4089        anyAction = -1 ;
4090        // To make depth available we may need a fake node
4091        CbcNode fakeNode;
4092        if (!currentNode_) {
4093            // Not true if sub trees assert (!numberNodes_);
4094            currentNode_ = &fakeNode;
4095        }
4096        phase_ = 3;
4097        // only allow 1000 passes
4098        int numberPassesLeft = 1000;
4099        // This is first crude step
4100        if (numberAnalyzeIterations_ && !parentModel_) {
4101            delete [] analyzeResults_;
4102            //int numberColumns = solver_->getNumCols();
4103            analyzeResults_ = new double [5*numberIntegers_];
4104            numberFixedAtRoot_ = newNode->analyze(this, analyzeResults_);
4105            if (numberFixedAtRoot_ > 0) {
4106              COIN_DETAIL_PRINT(printf("%d fixed by analysis\n", numberFixedAtRoot_));
4107                setPointers(solver_);
4108                numberFixedNow_ = numberFixedAtRoot_;
4109            } else if (numberFixedAtRoot_ < 0) {
4110              COIN_DETAIL_PRINT(printf("analysis found to be infeasible\n"));
4111                anyAction = -2;
4112                delete newNode ;
4113                newNode = NULL ;
4114                feasible = false ;
4115            }
4116        }
4117        OsiSolverBranch * branches = NULL;
4118        if (feasible)
4119          anyAction = chooseBranch(newNode, numberPassesLeft, NULL, cuts, resolved,
4120                                 NULL, NULL, NULL, branches);
4121        if (anyAction == -2 || newNode->objectiveValue() >= cutoff) {
4122            if (anyAction != -2) {
4123                // zap parent nodeInfo
4124#ifdef COIN_DEVELOP
4125                printf("zapping CbcNodeInfo %x\n", newNode->nodeInfo()->parent());
4126#endif
4127                if (newNode->nodeInfo())
4128                    newNode->nodeInfo()->nullParent();
4129            }
4130            delete newNode ;
4131            newNode = NULL ;
4132            feasible = false ;
4133        }
4134    }
4135    if (newNode && probingInfo_) {
4136        int number01 = probingInfo_->numberIntegers();
4137        //const fixEntry * entry = probingInfo_->fixEntries();
4138        const int * toZero = probingInfo_->toZero();
4139        //const int * toOne = probingInfo_->toOne();
4140        //const int * integerVariable = probingInfo_->integerVariable();
4141        if (toZero[number01]) {
4142            CglTreeProbingInfo info(*probingInfo_);
4143            if ((moreSpecialOptions2_&64)!=0&&!parentModel_) {
4144              /*
4145                Marginal idea. Further exploration probably good. Build some extra
4146                cliques from probing info. Not quite worth the effort?
4147              */
4148              CglProbing generator1;
4149              generator1.setUsingObjective(false);
4150              generator1.setMaxPass(1);
4151              generator1.setMaxPassRoot(1);
4152              generator1.setMaxLook(100);
4153              generator1.setRowCuts(3);
4154              generator1.setMaxElements(300);
4155              generator1.setMaxProbeRoot(solver_->getNumCols());
4156              CoinThreadRandom randomGenerator;
4157              //CglTreeProbingInfo info(solver_);
4158              info.level = 0;
4159              info.formulation_rows = solver_->getNumRows();
4160              info.inTree = false;
4161              if (parentModel_) {
4162                info.parentSolver=parentModel_->continuousSolver();
4163                // indicate if doing full search
4164                info.hasParent = ((specialOptions_&67108864)==0) ? 1 : 2; 
4165              } else {
4166                info.hasParent=0;
4167                info.parentSolver=NULL;
4168              }
4169              info.originalColumns=originalColumns();
4170              info.randomNumberGenerator=&randomGenerator;
4171              info.pass=4;
4172              generator1.setMode(8);
4173              OsiCuts cs;
4174              generator1.generateCutsAndModify(*solver_,cs,&info);
4175              // very clunky
4176              OsiSolverInterface * temp = generator1.cliqueModel(solver_,2);
4177              CglPreProcess dummy;
4178              OsiSolverInterface * newSolver=dummy.cliqueIt(*temp,0.0001);
4179              delete temp;
4180              OsiSolverInterface * fake = NULL;
4181              if (newSolver) {
4182#if 0
4183                int numberCliques = generator1.numberCliques();
4184                cliqueEntry * entry = generator1.cliqueEntry();
4185                cliqueType * type = new cliqueType [numberCliques];
4186                int * start = new int [numberCliques+1];
4187                start[numberCliques]=2*numberCliques;
4188                int n=0;
4189                for (int i=0;i<numberCliques;i++) {
4190                  start[i]=2*i;
4191                  setOneFixesInCliqueEntry(entry[2*i],true);
4192                  setOneFixesInCliqueEntry(entry[2*i+1],true);
4193                  type[i]=0;
4194                }
4195                fake = info.analyze(*solver_, 1,numberCliques,start,
4196                                    entry,type);
4197                delete [] type;
4198                delete [] entry;
4199#else
4200                fake = info.analyze(*newSolver, 1,-1);
4201#endif
4202                delete newSolver;
4203              } else {
4204                fake = info.analyze(*solver_, 1);
4205              }
4206              if (fake) {
4207                //fake->writeMps("fake");
4208                CglFakeClique cliqueGen(fake);
4209                cliqueGen.setStarCliqueReport(false);
4210                cliqueGen.setRowCliqueReport(false);
4211                cliqueGen.setMinViolation(0.1);
4212                addCutGenerator(&cliqueGen, 1, "Fake cliques", true, false, false, -200);
4213                generator_[numberCutGenerators_-1]->setTiming(true);
4214                for (int i = 0; i < numberCutGenerators_; i++) {
4215                  CglKnapsackCover * cutGen =
4216                  dynamic_cast<CglKnapsackCover *>(generator_[i]->generator());
4217                  if (cutGen) {
4218                    cutGen->createCliques(*fake,2,200,false);
4219                  }
4220                }
4221              }
4222            }
4223            if (probingInfo_->packDown()) {
4224#if CBC_USEFUL_PRINTING>1
4225                printf("%d implications on %d 0-1\n", toZero[number01], number01);
4226#endif
4227                // Create a cut generator that remembers implications discovered at root.
4228                CglImplication implication(probingInfo_);
4229                addCutGenerator(&implication, 1, "ImplicationCuts", true, false, false, -200);
4230                generator_[numberCutGenerators_-1]->setGlobalCuts(true);
4231                generator_[numberCutGenerators_-1]->setTiming(true);
4232            } else {
4233                delete probingInfo_;
4234                probingInfo_ = NULL;
4235            }
4236        } else {
4237            delete probingInfo_;
4238
4239            probingInfo_ = NULL;
4240        }
4241    }
4242    /*
4243      At this point, the root subproblem is infeasible or fathomed by bound
4244      (newNode == NULL), or we're live with an objective value that satisfies the
4245      current objective cutoff.
4246    */
4247    assert (!newNode || newNode->objectiveValue() <= cutoff) ;
4248    // Save address of root node as we don't want to delete it
4249    /*
4250      The common case is that the lp relaxation is feasible but doesn't satisfy
4251      integrality (i.e., newNode->branchingObject(), indicating we've been able to
4252      select a branching variable). Remove any cuts that have gone slack due to
4253      forcing monotone variables. Then tack on an CbcFullNodeInfo object and full
4254      basis (via createInfo()) and stash the new cuts in the nodeInfo (via
4255      addCuts()). If, by some miracle, we have an integral solution at the root
4256      (newNode->branchingObject() is NULL), takeOffCuts() will ensure that the solver holds
4257      a valid solution for use by setBestSolution().
4258    */
4259    CoinWarmStartBasis *lastws = NULL ;
4260    if (feasible && newNode->branchingObject()) {
4261        if (resolved) {
4262            takeOffCuts(cuts, false, NULL) ;
4263#     ifdef CHECK_CUT_COUNTS
4264            { printf("Number of rows after chooseBranch fix (root)"
4265                         "(active only) %d\n",
4266                         numberRowsAtContinuous_ + numberNewCuts_ + numberOldActiveCuts_) ;
4267                const CoinWarmStartBasis* debugws =
4268                    dynamic_cast <const CoinWarmStartBasis*>(solver_->getWarmStart()) ;
4269                debugws->print() ;
4270                delete debugws ;
4271            }
4272#     endif
4273        }
4274        //newNode->createInfo(this,NULL,NULL,NULL,NULL,0,0) ;
4275        //newNode->nodeInfo()->addCuts(cuts,
4276        //                       newNode->numberBranches(),whichGenerator_) ;
4277        if (lastws) delete lastws ;
4278        lastws = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
4279    }
4280    /*
4281      Continuous data to be used later
4282    */
4283    continuousObjective_ = solver_->getObjValue() * solver_->getObjSense();
4284    continuousInfeasibilities_ = 0 ;
4285    if (newNode) {
4286        continuousObjective_ = newNode->objectiveValue() ;
4287        delete [] continuousSolution_;
4288        continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
4289                                              numberColumns);
4290        continuousInfeasibilities_ = newNode->numberUnsatisfied() ;
4291    }
4292    /*
4293      Bound may have changed so reset in objects
4294    */
4295    {
4296        int i ;
4297        for (i = 0; i < numberObjects_; i++)
4298            object_[i]->resetBounds(solver_) ;
4299    }
4300    /*
4301      Feasible? Then we should have either a live node prepped for future
4302      expansion (indicated by variable() >= 0), or (miracle of miracles) an
4303      integral solution at the root node.
4304
4305      initializeInfo sets the reference counts in the nodeInfo object.  Since
4306      this node is still live, push it onto the heap that holds the live set.
4307    */
4308    if (newNode) {
4309        if (newNode->branchingObject()) {
4310            newNode->initializeInfo() ;
4311            tree_->push(newNode) ;
4312            // save pointer to root node - so can pick up bounds
4313            if (!topOfTree_)
4314              topOfTree_ = dynamic_cast<CbcFullNodeInfo *>(newNode->nodeInfo()) ;
4315            if (statistics_) {
4316                if (numberNodes2_ == maximumStatistics_) {
4317                    maximumStatistics_ = 2 * maximumStatistics_;
4318                    CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_];
4319                    memset(temp, 0, maximumStatistics_*sizeof(CbcStatistics *));
4320                    memcpy(temp, statistics_, numberNodes2_*sizeof(CbcStatistics *));
4321                    delete [] statistics_;
4322                    statistics_ = temp;
4323                }
4324                assert (!statistics_[numberNodes2_]);
4325                statistics_[numberNodes2_] = new CbcStatistics(newNode, this);
4326            }
4327            numberNodes2_++;
4328#     ifdef CHECK_NODE
4329            printf("Node %x on tree\n", newNode) ;
4330#     endif
4331        } else {
4332            // continuous is integer
4333            double objectiveValue = newNode->objectiveValue();
4334            setBestSolution(CBC_SOLUTION, objectiveValue,
4335                            solver_->getColSolution()) ;
4336            if (eventHandler) {
4337              // we are stopping anyway so no need to test return code
4338              eventHandler->event(CbcEventHandler::solution);
4339            }
4340            delete newNode ;
4341            newNode = NULL ;
4342        }
4343    }
4344
4345    if (printFrequency_ <= 0) {
4346        printFrequency_ = 1000 ;
4347        if (getNumCols() > 2000)
4348            printFrequency_ = 100 ;
4349    }
4350    /*
4351      It is possible that strong branching fixes one variable and then the code goes round
4352      again and again.  This can take too long.  So we need to warn user - just once.
4353    */
4354    numberLongStrong_ = 0;
4355    CbcNode * createdNode = NULL;
4356#ifdef CBC_THREAD
4357    if ((specialOptions_&2048) != 0)
4358        numberThreads_ = 0;
4359    if (numberThreads_ ) {
4360        nodeCompare_->sayThreaded(); // need to use addresses
4361        master_ = new CbcBaseModel(*this,
4362                                   (parallelMode() < -1) ? 1 : 0);
4363        masterThread_ = master_->masterThread();
4364    }
4365#endif
4366#ifdef COIN_HAS_CLP
4367    {
4368        OsiClpSolverInterface * clpSolver
4369        = dynamic_cast<OsiClpSolverInterface *> (solver_);
4370        if (clpSolver && !parentModel_) {
4371            clpSolver->computeLargestAway();
4372        }
4373    }
4374#endif
4375    /*
4376      At last, the actual branch-and-cut search loop, which will iterate until
4377      the live set is empty or we hit some limit (integrality gap, time, node
4378      count, etc.). The overall flow is to rebuild a subproblem, reoptimise using
4379      solveWithCuts(), choose a branching pattern with chooseBranch(), and finally
4380      add the node to the live set.
4381
4382      The first action is to winnow the live set to remove nodes which are worse
4383      than the current objective cutoff.
4384    */
4385    if (solver_->getRowCutDebuggerAlways()) {
4386        OsiRowCutDebugger * debuggerX = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
4387        const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger() ;
4388        if (!debugger) {
4389            // infeasible!!
4390            printf("before search\n");
4391            const double * lower = solver_->getColLower();
4392            const double * upper = solver_->getColUpper();
4393            const double * solution = debuggerX->optimalSolution();
4394            int numberColumns = solver_->getNumCols();
4395            for (int i = 0; i < numberColumns; i++) {
4396                if (solver_->isInteger(i)) {
4397                    if (solution[i] < lower[i] - 1.0e-6 || solution[i] > upper[i] + 1.0e-6)
4398                        printf("**** ");
4399                    printf("%d %g <= %g <= %g\n",
4400                           i, lower[i], solution[i], upper[i]);
4401                }
4402            }
4403            //abort();
4404        }
4405    }
4406    {
4407        // may be able to change cutoff now
4408        double cutoff = getCutoff();
4409        double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
4410        if (cutoff > bestObjective_ - increment) {
4411            cutoff = bestObjective_ - increment ;
4412            setCutoff(cutoff) ;
4413        }
4414    }
4415#ifdef CBC_THREAD
4416    bool goneParallel = false;
4417#endif
4418#define MAX_DEL_NODE 1
4419    CbcNode * delNode[MAX_DEL_NODE+1];
4420    int nDeleteNode = 0;
4421    // For Printing etc when parallel
4422    int lastEvery1000 = 0;
4423    int lastPrintEvery = 0;
4424    int numberConsecutiveInfeasible = 0;
4425#define PERTURB_IN_FATHOM
4426#ifdef PERTURB_IN_FATHOM
4427    // allow in fathom
4428    if ((moreSpecialOptions_& 262144) != 0)
4429      specialOptions_ |= 131072;
4430#endif
4431    while (true) {
4432        lockThread();
4433#ifdef COIN_HAS_CLP
4434        // See if we want dantzig row choice
4435        goToDantzig(100, savePivotMethod);
4436#endif
4437        //#define REPORT_DYNAMIC 2
4438#if REPORT_DYNAMIC
4439        if (numberNodes_&&!parentModel_&&(tree_->empty()||(numberNodes_%10000)==0)) {
4440          // Get average up and down costs
4441          double averageUp = 0.0;
4442          double averageDown = 0.0;
4443          int numberUp = 0;
4444          int numberDown = 0;
4445          int minTimesDown = COIN_INT_MAX;
4446          int maxTimesDown = 0;
4447          int neverBranchedDown = 0;
4448          int infeasibleTimesDown = 0;
4449          int minTimesUp = COIN_INT_MAX;
4450          int maxTimesUp = 0;
4451          int infeasibleTimesUp = 0;
4452          int neverBranchedUp = 0;
4453          int neverBranched = 0;
4454          int i;
4455          int numberInts=0;
4456          bool endOfSearch = tree_->empty();
4457          int  numberUp2 = 0;
4458          int numberDown2 = 0;
4459          for ( i = 0; i < numberObjects_; i++) {
4460            OsiObject * object = object_[i];
4461            CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
4462              dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
4463            if (dynamicObject) {
4464              numberInts++;
4465              if (dynamicObject->numberTimesUp()||
4466                  dynamicObject->numberTimesDown()) {
4467                int  nUp = 0;
4468                int nDown = 0;
4469                double up = 0.0;
4470                double down = 0.0;
4471                if (dynamicObject->numberTimesUp()) {
4472                  numberUp++;
4473                  nUp = dynamicObject->numberTimesUp();
4474                  minTimesUp = CoinMin(minTimesUp,nUp);
4475                  maxTimesUp = CoinMax(maxTimesUp,nUp);
4476                  up = dynamicObject->upDynamicPseudoCost();
4477                  averageUp += up;
4478                  numberUp2 += nUp;
4479                  infeasibleTimesUp += dynamicObject->numberTimesUpInfeasible();
4480                } else {
4481                  neverBranchedUp++;
4482                }
4483                if (dynamicObject->numberTimesDown()) {
4484                  numberDown++;
4485                  nDown = dynamicObject->numberTimesDown();
4486                  minTimesDown = CoinMin(minTimesDown,nDown);
4487                  maxTimesDown = CoinMax(maxTimesDown,nDown);
4488                  down = dynamicObject->downDynamicPseudoCost();
4489                  averageDown += down;
4490                  numberDown2 += dynamicObject->numberTimesDown();
4491                  infeasibleTimesDown += dynamicObject->numberTimesDownInfeasible();
4492                } else {
4493                  neverBranchedDown++;
4494                }
4495#if REPORT_DYNAMIC > 1
4496#if REPORT_DYNAMIC == 2
4497                if (endOfSearch&&numberIntegers_<400) {
4498#elif REPORT_DYNAMIC == 3
4499                if (endOfSearch) {
4500#else
4501                  {
4502#endif
4503                  dynamicObject->print(0,0.0);
4504                }
4505#endif
4506              } else {
4507                neverBranched++;
4508#if REPORT_DYNAMIC > 2
4509#if REPORT_DYNAMIC == 3
4510                if (endOfSearch&&numberIntegers_<400) {
4511#elif REPORT_DYNAMIC == 4
4512                if (endOfSearch) {
4513#else
4514                  {
4515#endif
4516                  printf("col %d - never branched on\n",dynamicObject->columnNumber());
4517                }
4518#endif
4519              }
4520            }
4521          }
4522          if (numberUp)
4523            averageUp /= static_cast<double> (numberUp);
4524          else
4525            averageUp = 0.0;
4526          if (numberDown)
4527            averageDown /= static_cast<double> (numberDown);
4528          else
4529            averageDown = 0.0;
4530          printf("Report for %d variables (%d never branched on) after %d nodes - total solves down %d up %d\n",
4531                 numberInts,neverBranched,numberNodes_,numberDown2,numberUp2);
4532          if ((neverBranchedDown||neverBranchedUp)&&endOfSearch)
4533            printf("odd %d never branched down and %d never branched up\n",
4534                   neverBranchedDown,neverBranchedUp);
4535          printf("down average %g times (%d infeasible) average increase %g min/max times (%d,%d)\n",
4536                 static_cast<double>(numberDown2)/numberDown,infeasibleTimesDown,averageDown,
4537                 minTimesDown,maxTimesDown);
4538          printf("up average %g times (%d infeasible) average increase %g min/max times (%d,%d)\n",
4539                 static_cast<double>(numberUp2)/numberUp,infeasibleTimesUp,averageUp,
4540                 minTimesUp,maxTimesUp);
4541        }
4542#endif
4543        if (tree_->empty()) {
4544#ifdef CBC_THREAD
4545            if (parallelMode() > 0 && master_) {
4546                int anyLeft = master_->waitForThreadsInTree(0);
4547                if (!anyLeft) {
4548                    master_->stopThreads(-1);
4549                    break;
4550                }
4551            } else {
4552                break;
4553            }
4554#else
4555            break;
4556#endif
4557        } else {
4558            unlockThread();
4559        }
4560        // If done 50/100 nodes see if worth trying reduction
4561        if (numberNodes_ >= nextCheckRestart) {
4562          if (nextCheckRestart<100)
4563            nextCheckRestart=100;
4564          else
4565            nextCheckRestart=COIN_INT_MAX;
4566#ifdef COIN_HAS_CLP
4567            OsiClpSolverInterface * clpSolver
4568            = dynamic_cast<OsiClpSolverInterface *> (solver_);
4569            if (clpSolver && ((specialOptions_&131072) == 0) && true) {
4570                ClpSimplex * simplex = clpSolver->getModelPtr();
4571                int perturbation = simplex->perturbation();
4572#if CBC_USEFUL_PRINTING>1
4573                printf("Testing its n,s %d %d solves n,s %d %d - pert %d\n",
4574                       numberIterations_, saveNumberIterations,
4575                       numberSolves_, saveNumberSolves, perturbation);
4576#endif
4577                if (perturbation == 50 && (numberIterations_ - saveNumberIterations) <
4578                        8*(numberSolves_ - saveNumberSolves)) {
4579                    // switch off perturbation
4580                    simplex->setPerturbation(100);
4581#if CBC_USEFUL_PRINTING>1
4582                    printf("Perturbation switched off\n");
4583#endif
4584                }
4585            }
4586#endif
4587            /*
4588              Decide if we want to do a restart.
4589            */
4590            if (saveSolver && (specialOptions_&(512 + 32768)) != 0) {
4591                bool tryNewSearch = solverCharacteristics_->reducedCostsAccurate() &&
4592                                    (getCutoff() < 1.0e20 && getCutoff() < checkCutoffForRestart);
4593                int numberColumns = getNumCols();
4594                if (tryNewSearch) {
4595                    // adding increment back allows current best - tiny bit weaker
4596                    checkCutoffForRestart = getCutoff() + getCutoffIncrement() ;
4597#if CBC_USEFUL_PRINTING>1
4598                    printf("after %d nodes, cutoff %g - looking\n",
4599                           numberNodes_, getCutoff());
4600#endif
4601                    saveSolver->resolve();
4602                    double direction = saveSolver->getObjSense() ;
4603                    double gap = checkCutoffForRestart - saveSolver->getObjValue() * direction ;
4604                    double tolerance;
4605                    saveSolver->getDblParam(OsiDualTolerance, tolerance) ;
4606                    if (gap <= 0.0)
4607                        gap = tolerance;
4608                    gap += 100.0 * tolerance;
4609                    double integerTolerance = getDblParam(CbcIntegerTolerance) ;
4610
4611                    const double *lower = saveSolver->getColLower() ;
4612                    const double *upper = saveSolver->getColUpper() ;
4613                    const double *solution = saveSolver->getColSolution() ;
4614                    const double *reducedCost = saveSolver->getReducedCost() ;
4615
4616                    int numberFixed = 0 ;
4617                    int numberFixed2 = 0;
4618#ifdef COIN_DEVELOP
4619                    printf("gap %g\n", gap);
4620#endif
4621                    for (int i = 0 ; i < numberIntegers_ ; i++) {
4622                        int iColumn = integerVariable_[i] ;
4623                        double djValue = direction * reducedCost[iColumn] ;
4624                        if (upper[iColumn] - lower[iColumn] > integerTolerance) {
4625                            if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) {
4626                                //printf("%d to lb on dj of %g - bounds %g %g\n",
4627                                //     iColumn,djValue,lower[iColumn],upper[iColumn]);
4628                                saveSolver->setColUpper(iColumn, lower[iColumn]) ;
4629                                numberFixed++ ;
4630                            } else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) {
4631                                //printf("%d to ub on dj of %g - bounds %g %g\n",
4632                                //     iColumn,djValue,lower[iColumn],upper[iColumn]);
4633                                saveSolver->setColLower(iColumn, upper[iColumn]) ;
4634                                numberFixed++ ;
4635                            }
4636                        } else {
4637                            //printf("%d has dj of %g - already fixed to %g\n",
4638                            //     iColumn,djValue,lower[iColumn]);
4639                            numberFixed2++;
4640                        }
4641                    }
4642#ifdef COIN_DEVELOP
4643                    if ((specialOptions_&1) != 0) {
4644                        const OsiRowCutDebugger *debugger = saveSolver->getRowCutDebugger() ;
4645                        if (debugger) {
4646                            printf("Contains optimal\n") ;
4647                            OsiSolverInterface * temp = saveSolver->clone();
4648                            const double * solution = debugger->optimalSolution();
4649                            const double *lower = temp->getColLower() ;
4650                            const double *upper = temp->getColUpper() ;
4651                            int n = temp->getNumCols();
4652                            for (int i = 0; i < n; i++) {
4653                                if (temp->isInteger(i)) {
4654                                    double value = floor(solution[i] + 0.5);
4655                                    assert (value >= lower[i] && value <= upper[i]);
4656                                    temp->setColLower(i, value);
4657                                    temp->setColUpper(i, value);
4658                                }
4659                            }
4660                            temp->writeMps("reduced_fix");
4661                            delete temp;
4662                            saveSolver->writeMps("reduced");
4663                        } else {
4664                            abort();
4665                        }
4666                    }
4667                    printf("Restart could fix %d integers (%d already fixed)\n",
4668                           numberFixed + numberFixed2, numberFixed2);
4669#endif
4670                    numberFixed += numberFixed2;
4671                    if (numberFixed*10 < numberColumns && numberFixed*4 < numberIntegers_)
4672                        tryNewSearch = false;
4673                }
4674#ifdef CONFLICT_CUTS
4675                // temporary
4676                //if ((moreSpecialOptions_&4194304)!=0)
4677                //tryNewSearch=false;
4678#endif
4679                if (tryNewSearch) {
4680                    // back to solver without cuts?
4681                    OsiSolverInterface * solver2 = saveSolver->clone();
4682                    const double *lower = saveSolver->getColLower() ;
4683                    const double *upper = saveSolver->getColUpper() ;
4684                    for (int i = 0 ; i < numberIntegers_ ; i++) {
4685                        int iColumn = integerVariable_[i] ;
4686                        solver2->setColLower(iColumn, lower[iColumn]);
4687                        solver2->setColUpper(iColumn, upper[iColumn]);
4688                    }
4689                    // swap
4690                    delete saveSolver;
4691                    saveSolver = solver2;
4692                    double * newSolution = new double[numberColumns];
4693                    double objectiveValue = checkCutoffForRestart;
4694                    // Save the best solution so far.
4695                    CbcSerendipity heuristic(*this);
4696                    if (bestSolution_)
4697                        heuristic.setInputSolution(bestSolution_, bestObjective_);
4698                    // Magic number
4699                    heuristic.setFractionSmall(0.8);
4700                    // `pumpTune' to stand-alone solver for explanations.
4701                    heuristic.setFeasibilityPumpOptions(1008013);
4702                    // Use numberNodes to say how many are original rows
4703                    heuristic.setNumberNodes(continuousSolver_->getNumRows());
4704#ifdef COIN_DEVELOP
4705                    if (continuousSolver_->getNumRows() <
4706                            solver_->getNumRows())
4707                        printf("%d rows added ZZZZZ\n",
4708                               solver_->getNumRows() - continuousSolver_->getNumRows());
4709#endif
4710                    int returnCode = heuristic.smallBranchAndBound(saveSolver,
4711                                     -1, newSolution,
4712                                     objectiveValue,
4713                                     checkCutoffForRestart, "Reduce");
4714                    if (returnCode < 0) {
4715#ifdef COIN_DEVELOP
4716                        printf("Restart - not small enough to do search after fixing\n");
4717#endif
4718                        delete [] newSolution;
4719                    } else {
4720                        // 1 for sol'n, 2 for finished, 3 for both
4721                        if ((returnCode&1) != 0) {
4722                            // increment number of solutions so other heuristics can test
4723                            numberSolutions_++;
4724                            numberHeuristicSolutions_++;
4725                            lastHeuristic_ = NULL;
4726                            setBestSolution(CBC_ROUNDING, objectiveValue, newSolution) ;
4727                        }
4728                        delete [] newSolution;
4729#ifdef CBC_THREAD
4730                        if (master_) {
4731                            lockThread();
4732                            if (parallelMode() > 0) {
4733                                while (master_->waitForThreadsInTree(0)) {
4734                                    lockThread();
4735                                    double dummyBest;
4736                                    tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4737                                    //unlockThread();
4738                                }
4739                            } else {
4740                                double dummyBest;
4741                                tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4742                            }
4743                            master_->waitForThreadsInTree(2);
4744                            delete master_;
4745                            master_ = NULL;
4746                            masterThread_ = NULL;
4747                        }
4748#endif
4749                        if (tree_->size()) {
4750                            double dummyBest;
4751                            tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4752                        }
4753                        break;
4754                    }
4755                }
4756                delete saveSolver;
4757                saveSolver = NULL;
4758            }
4759        }
4760        /*
4761          Check for abort on limits: node count, solution count, time, integrality gap.
4762        */
4763        if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
4764                numberSolutions_ < intParam_[CbcMaxNumSol] &&
4765                !maximumSecondsReached() &&
4766                !stoppedOnGap_ && !eventHappened_ && (maximumNumberIterations_ < 0 ||
4767                                                      numberIterations_ < maximumNumberIterations_))) {
4768            // out of loop
4769            break;
4770        }
4771#ifdef BONMIN
4772        assert(!solverCharacteristics_->solutionAddsCuts() || solverCharacteristics_->mipFeasible());
4773#endif
4774// Sets percentage of time when we try diving. Diving requires a bit of heap reorganisation, because
4775// we need to replace the comparison function to dive, and that requires reordering to retain the
4776// heap property.
4777#define DIVE_WHEN 1000
4778#define DIVE_STOP 2000
4779        int kNode = numberNodes_ % 4000;
4780        if (numberNodes_<100000 && kNode>DIVE_WHEN && kNode <= DIVE_STOP) {
4781            if (!parallelMode()) {
4782                if (kNode == DIVE_WHEN + 1 || numberConsecutiveInfeasible > 1) {
4783                    CbcCompareDefault * compare = dynamic_cast<CbcCompareDefault *>
4784                                                  (nodeCompare_);
4785                    // Don't interfere if user has replaced the compare function.
4786                    if (compare) {
4787                        //printf("Redoing tree\n");
4788                        compare->startDive(this);
4789                        numberConsecutiveInfeasible = 0;
4790                    }
4791                }
4792            }
4793        }
4794        // replace current cutoff?
4795        if (cutoff > getCutoff()) {
4796            double newCutoff = getCutoff();
4797            if (analyzeResults_) {
4798                // see if we could fix any (more)
4799                int n = 0;
4800                double * newLower = analyzeResults_;
4801                double * objLower = newLower + numberIntegers_;
4802                double * newUpper = objLower + numberIntegers_;
4803                double * objUpper = newUpper + numberIntegers_;
4804                for (int i = 0; i < numberIntegers_; i++) {
4805                    if (objLower[i] > newCutoff) {
4806                        n++;
4807                        if (objUpper[i] > newCutoff) {
4808                            newCutoff = -COIN_DBL_MAX;
4809                            break;
4810                        }
4811                        // add as global cut
4812                        objLower[i]=-COIN_DBL_MAX;
4813                        OsiRowCut rc;
4814                        rc.setLb(newLower[i]);
4815                        rc.setUb(COIN_DBL_MAX);
4816                        double one=1.0;
4817                        rc.setRow(1,integerVariable_+i,&one,false);
4818                        rc.setGloballyValidAsInteger(2);
4819                        globalCuts_.addCutIfNotDuplicate(rc) ;
4820                    } else if (objUpper[i] > newCutoff) {
4821                        n++;
4822                        // add as global cut
4823                        objUpper[i]=-COIN_DBL_MAX;
4824                        OsiRowCut rc;
4825                        rc.setLb(-COIN_DBL_MAX);
4826                        rc.setUb(newUpper[i]);
4827                        double one=1.0;
4828                        rc.setRow(1,integerVariable_+i,&one,false);
4829                        rc.setGloballyValidAsInteger(2);
4830                        globalCuts_.addCutIfNotDuplicate(rc) ;
4831                    }
4832                }
4833                if (newCutoff == -COIN_DBL_MAX) {
4834                  COIN_DETAIL_PRINT(printf("Root analysis says finished\n"));
4835                } else if (n > numberFixedNow_) {
4836                  COIN_DETAIL_PRINT(printf("%d more fixed by analysis - now %d\n", n - numberFixedNow_, n));
4837                    numberFixedNow_ = n;
4838                }
4839            }
4840            if (eventHandler) {
4841                if (!eventHandler->event(CbcEventHandler::solution)) {
4842                    eventHappened_ = true; // exit
4843                }
4844                newCutoff = getCutoff();
4845            }
4846            lockThread();
4847            /*
4848              Clean the tree to reflect the new solution, then see if the
4849              node comparison predicate wants to make any changes. If so,
4850              call setComparison for the side effect of rebuilding the heap.
4851            */
4852            tree_->cleanTree(this,newCutoff,bestPossibleObjective_) ;
4853            if (nodeCompare_->newSolution(this) ||
4854                nodeCompare_->newSolution(this,continuousObjective_,
4855                                          continuousInfeasibilities_)) {
4856              tree_->setComparison(*nodeCompare_) ;
4857            }
4858            if (tree_->empty()) {
4859                continue;
4860            }
4861            unlockThread();
4862        }
4863        cutoff = getCutoff() ;
4864        /*
4865            Periodic activities: Opportunities to
4866            + tweak the nodeCompare criteria,
4867            + check if we've closed the integrality gap enough to quit,
4868            + print a summary line to let the user know we're working
4869        */
4870        if (numberNodes_ >= lastEvery1000) {
4871            lockThread();
4872#ifdef COIN_HAS_CLP
4873            // See if we want dantzig row choice
4874            goToDantzig(1000, savePivotMethod);
4875#endif
4876            lastEvery1000 = numberNodes_ + 1000;
4877            bool redoTree = nodeCompare_->every1000Nodes(this, numberNodes_) ;
4878#ifdef CHECK_CUT_SIZE
4879            verifyCutSize (tree_, *this);
4880#endif
4881            // redo tree if requested
4882            if (redoTree)
4883                tree_->setComparison(*nodeCompare_) ;
4884            unlockThread();
4885        }
4886        // Had hotstart before, now switched off
4887        if (saveCompare && !hotstartSolution_) {
4888            // hotstart switched off
4889            delete nodeCompare_; // off depth first
4890            nodeCompare_ = saveCompare;
4891            saveCompare = NULL;
4892            // redo tree
4893            lockThread();
4894            tree_->setComparison(*nodeCompare_) ;
4895            unlockThread();
4896        }
4897        if (numberNodes_ >= lastPrintEvery) {
4898            lastPrintEvery = numberNodes_ + printFrequency_;
4899            lockThread();
4900            int nNodes = tree_->size() ;
4901
4902            //MODIF PIERRE
4903            bestPossibleObjective_ = tree_->getBestPossibleObjective();
4904#ifdef CBC_THREAD
4905            if (parallelMode() > 0 && master_) {
4906              // need to adjust for ones not on tree
4907              int numberThreads = master_->numberThreads();
4908              for (int i=0;i<numberThreads;i++) {
4909                CbcThread * child = master_->child(i);
4910                if (child->node()) {
4911                  // adjust
4912                  double value = child->node()->objectiveValue();
4913                  bestPossibleObjective_ = CoinMin(bestPossibleObjective_, value);
4914                }
4915              }
4916            }
4917#endif
4918            unlockThread();
4919#if CBC_USEFUL_PRINTING>1
4920            if (getCutoff() < 1.0e20) {
4921                if (fabs(getCutoff() - (bestObjective_ - getCutoffIncrement())) > 1.0e-6 &&
4922                        !parentModel_)
4923                    printf("model cutoff in status %g, best %g, increment %g\n",
4924                           getCutoff(), bestObjective_, getCutoffIncrement());
4925                assert (getCutoff() < bestObjective_ - getCutoffIncrement() +
4926                        1.0e-6 + 1.0e-10*fabs(bestObjective_));
4927            }
4928#endif
4929            if (!intParam_[CbcPrinting]) {
4930                // Parallel may not have any nodes
4931                if (!nNodes) 
4932                  bestPossibleObjective_ = lastBestPossibleObjective;
4933                else
4934                  lastBestPossibleObjective = bestPossibleObjective_;
4935                messageHandler()->message(CBC_STATUS, messages())
4936                  << numberNodes_ << CoinMax(nNodes,1) << bestObjective_ << bestPossibleObjective_
4937                << getCurrentSeconds()
4938                << CoinMessageEol ;
4939            } else if (intParam_[CbcPrinting] == 1) {
4940                messageHandler()->message(CBC_STATUS2, messages())
4941                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4942                << tree_->lastDepth() << tree_->lastUnsatisfied() 
4943                << tree_->lastObjective() << numberIterations_
4944                << getCurrentSeconds()
4945                << CoinMessageEol ;
4946            } else if (!numberExtraIterations_) {
4947                messageHandler()->message(CBC_STATUS2, messages())
4948                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4949                << tree_->lastDepth() << tree_->lastUnsatisfied() << numberIterations_
4950                << getCurrentSeconds()
4951                << CoinMessageEol ;
4952            } else {
4953                messageHandler()->message(CBC_STATUS3, messages())
4954                << numberNodes_ << numberFathoms_ << numberExtraNodes_ << nNodes
4955                << bestObjective_ << bestPossibleObjective_
4956                << tree_->lastDepth() << tree_->lastUnsatisfied() << numberIterations_ << numberExtraIterations_
4957                << getCurrentSeconds()
4958                << CoinMessageEol ;
4959            }
4960#ifdef COIN_HAS_NTY
4961            if (symmetryInfo_) 
4962              symmetryInfo_->statsOrbits(this,1);
4963#endif
4964#if PRINT_CONFLICT==1
4965            if (numberConflictCuts>lastNumberConflictCuts) {
4966              double length = lengthConflictCuts/numberConflictCuts;
4967              printf("%d new conflict cuts - total %d - average length %g\n",
4968                     numberConflictCuts-lastNumberConflictCuts,
4969                     numberConflictCuts,length);
4970              lastNumberConflictCuts = numberConflictCuts;
4971            }
4972#endif
4973            if (eventHandler && !eventHandler->event(CbcEventHandler::treeStatus)) {
4974                eventHappened_ = true; // exit
4975            }
4976        }
4977        // See if can stop on gap
4978        if(canStopOnGap()) {
4979            stoppedOnGap_ = true ;
4980        }
4981
4982#ifdef CHECK_NODE_FULL
4983        verifyTreeNodes(tree_, *this) ;
4984#   endif
4985#   ifdef CHECK_CUT_COUNTS
4986        verifyCutCounts(tree_, *this) ;
4987#   endif
4988        /*
4989          Now we come to the meat of the loop. To create the active subproblem, we'll
4990          pop the most promising node in the live set, rebuild the subproblem it
4991          represents, and then execute the current arm of the branch to create the
4992          active subproblem.
4993        */
4994        CbcNode * node = NULL;
4995#ifdef CBC_THREAD
4996        if (!parallelMode() || parallelMode() == -1) {
4997#endif
4998            node = tree_->bestNode(cutoff) ;
4999            // Possible one on tree worse than cutoff
5000            // Weird comparison function can leave ineligible nodes on tree
5001            if (!node || node->objectiveValue() > cutoff)
5002                continue;
5003            // Do main work of solving node here
5004            doOneNode(this, node, createdNode);
5005#ifdef JJF_ZERO
5006            if (node) {
5007              if (createdNode) {
5008                printf("Node %d depth %d, created %d depth %d\n",
5009                       node->nodeNumber(), node->depth(),
5010                       createdNode->nodeNumber(), createdNode->depth());
5011              } else {
5012                printf("Node %d depth %d,  no created node\n",
5013                       node->nodeNumber(), node->depth());
5014              }
5015            } else if (createdNode) {
5016              printf("Node exhausted, created %d depth %d\n",
5017                     createdNode->nodeNumber(), createdNode->depth());
5018            } else {
5019              printf("Node exhausted,  no created node\n");
5020              numberConsecutiveInfeasible = 2;
5021            }
5022#endif
5023            //if (createdNode)
5024            //numberConsecutiveInfeasible=0;
5025            //else
5026            //numberConsecutiveInfeasible++;
5027#ifdef CBC_THREAD
5028        } else if (parallelMode() > 0) {
5029            //lockThread();
5030            //node = tree_->bestNode(cutoff) ;
5031            // Possible one on tree worse than cutoff
5032            if (true || !node || node->objectiveValue() > cutoff) {
5033                assert (master_);
5034                if (master_) {
5035                    int anyLeft = master_->waitForThreadsInTree(1);
5036                    // may need to go round again
5037                    if (anyLeft) {
5038                        continue;
5039                    } else {
5040                        master_->stopThreads(-1);
5041                    }
5042                }
5043            }
5044            //unlockThread();
5045        } else {
5046            // Deterministic parallel
5047          if ((tree_->size() < CoinMax(numberThreads_, 8)||
5048               hotstartSolution_) && !goneParallel) {
5049                node = tree_->bestNode(cutoff) ;
5050                // Possible one on tree worse than cutoff
5051                if (!node || node->objectiveValue() > cutoff)
5052                    continue;
5053                // Do main work of solving node here
5054                doOneNode(this, node, createdNode);
5055                assert (createdNode);
5056                if (!createdNode->active()) {
5057                    delete createdNode;
5058                    createdNode = NULL;
5059                } else {
5060                    // Say one more pointing to this
5061                    node->nodeInfo()->increment() ;
5062                    tree_->push(createdNode) ;
5063                }
5064                if (node->active()) {
5065                    assert (node->nodeInfo());
5066                    if (node->nodeInfo()->numberBranchesLeft()) {
5067                        tree_->push(node) ;
5068                    } else {
5069                        node->setActive(false);
5070                    }
5071                } else {
5072                    if (node->nodeInfo()) {
5073                        if (!node->nodeInfo()->numberBranchesLeft())
5074                            node->nodeInfo()->allBranchesGone(); // can clean up
5075                        // So will delete underlying stuff
5076                        node->setActive(true);
5077                    }
5078                    delNode[nDeleteNode++] = node;
5079                    node = NULL;
5080                }
5081                if (nDeleteNode >= MAX_DEL_NODE) {
5082                    for (int i = 0; i < nDeleteNode; i++) {
5083                        //printf("trying to del %d %x\n",i,delNode[i]);
5084                        delete delNode[i];
5085                        //printf("done to del %d %x\n",i,delNode[i]);
5086                    }
5087                    nDeleteNode = 0;
5088                }
5089            } else {
5090                // Split and solve
5091                master_->deterministicParallel();
5092                goneParallel = true;
5093            }
5094        }
5095#endif
5096    }
5097    if (nDeleteNode) {
5098        for (int i = 0; i < nDeleteNode; i++) {
5099            delete delNode[i];
5100        }
5101        nDeleteNode = 0;
5102    }
5103#ifdef CBC_THREAD
5104    if (master_) {
5105        master_->stopThreads(-1);
5106        master_->waitForThreadsInTree(2);
5107        // adjust time to allow for children on some systems
5108        //dblParam_[CbcStartSeconds] -= CoinCpuTimeJustChildren();
5109    }
5110#endif
5111    /*
5112      End of the non-abort actions. The next block of code is executed if we've
5113      aborted because we hit one of the limits. Clean up by deleting the live set
5114      and break out of the node processing loop. Note that on an abort, node may
5115      have been pushed back onto the tree for further processing, in which case
5116      it'll be deleted in cleanTree. We need to check.
5117    */
5118    if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
5119          numberSolutions_ < intParam_[CbcMaxNumSol] &&
5120          !maximumSecondsReached() &&
5121          !stoppedOnGap_ &&
5122          !eventHappened_ &&
5123          (maximumNumberIterations_ < 0 || numberIterations_ < maximumNumberIterations_))
5124         ) {
5125        if (tree_->size()) {
5126            double dummyBest;
5127            tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
5128#if 0 // Does not seem to be needed def CBC_THREAD
5129            if (parallelMode() > 0 && master_) {
5130              // see if any dangling nodes
5131              int numberThreads = master_->numberThreads();
5132              for (int i=0;i<numberThreads;i++) {
5133                CbcThread * child = master_->child(i);
5134                //if (child->createdNode())
5135                //printf("CHILD_NODE %p\n",child->createdNode());
5136                delete child->createdNode();
5137              }
5138            }
5139#endif
5140        }
5141        delete nextRowCut_;
5142        /* order is important here:
5143         * maximumSecondsReached() should be checked before eventHappened_ and
5144         * isNodeLimitReached() should be checked after eventHappened_
5145         * reason is, that at timelimit, eventHappened_ is set to true to make Cbc stop fast
5146         *   and if Ctrl+C is hit, then the nodelimit is set to -1 to make Cbc stop
5147         */
5148        if (stoppedOnGap_) {
5149            messageHandler()->message(CBC_GAP, messages())
5150            << bestObjective_ - bestPossibleObjective_
5151            << dblParam_[CbcAllowableGap]
5152            << dblParam_[CbcAllowableFractionGap]*100.0
5153            << CoinMessageEol ;
5154            secondaryStatus_ = 2;
5155            status_ = 0 ;
5156        } else if (maximumSecondsReached()) {
5157            handler_->message(CBC_MAXTIME, messages_) << CoinMessageEol ;
5158            secondaryStatus_ = 4;
5159            status_ = 1 ;
5160        } else if (numberSolutions_ >= intParam_[CbcMaxNumSol]) {
5161            handler_->message(CBC_MAXSOLS, messages_) << CoinMessageEol ;
5162            secondaryStatus_ = 6;
5163            status_ = 1 ;
5164        } else if (isNodeLimitReached()) {
5165            handler_->message(CBC_MAXNODES, messages_) << CoinMessageEol ;
5166            secondaryStatus_ = 3;
5167            status_ = 1 ;
5168        } else if (maximumNumberIterations_ >= 0 && numberIterations_ >= maximumNumberIterations_) {
5169            handler_->message(CBC_MAXITERS, messages_) << CoinMessageEol ;
5170            secondaryStatus_ = 8;
5171            status_ = 1 ;
5172        } else {
5173            handler_->message(CBC_EVENT, messages_) << CoinMessageEol ;
5174            secondaryStatus_ = 5;
5175            status_ = 5 ;
5176        }
5177    }
5178#ifdef CBC_THREAD
5179    if (master_) {
5180        delete master_;
5181        master_ = NULL;
5182        masterThread_ = NULL;
5183    }
5184#endif
5185    /*
5186      That's it, we've exhausted the search tree, or broken out of the loop because
5187      we hit some limit on evaluation.
5188
5189      We may have got an intelligent tree so give it one more chance
5190    */
5191    // Tell solver we are not in Branch and Cut
5192    solver_->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo, NULL) ;
5193    tree_->endSearch();
5194    //  If we did any sub trees - did we give up on any?
5195    if ( numberStoppedSubTrees_)
5196        status_ = 1;
5197    numberNodes_ += numberExtraNodes_;
5198    numberIterations_ += numberExtraIterations_;
5199    if (eventHandler) {
5200        eventHandler->event(CbcEventHandler::endSearch);
5201    }
5202    if (!status_) {
5203        // Set best possible unless stopped on gap
5204        if (secondaryStatus_ != 2)
5205            bestPossibleObjective_ = bestObjective_;
5206        handler_->message(CBC_END_GOOD, messages_)
5207        << bestObjective_ << numberIterations_ << numberNodes_ << getCurrentSeconds()
5208        << CoinMessageEol ;
5209    } else {
5210        handler_->message(CBC_END, messages_)
5211        << bestObjective_ << bestPossibleObjective_
5212        << numberIterations_ << numberNodes_ << getCurrentSeconds()
5213        << CoinMessageEol ;
5214    }
5215    if ((moreSpecialOptions_&4194304)!=0) {
5216      // Conflict cuts
5217      int numberCuts = globalCuts_.sizeRowCuts();
5218      int nConflict=0;
5219      double sizeConflict = 0.0;
5220      for (int i=0;i<numberCuts;i++) {
5221        OsiRowCut2 * cut = globalCuts_.cut(i);
5222        if (cut->whichRow()==1) {
5223          nConflict++;