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

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

more changes for AddIntegers?

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