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

Last change on this file since 2199 was 2199, checked in by forrest, 4 years ago

fix resizeWhichGenerator bug hopefully

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