source: stable/2.8/Cbc/src/CbcModel.cpp @ 1883

Last change on this file since 1883 was 1883, checked in by stefan, 7 years ago

sync with trunk rev 1882

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