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

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

changes for non simplex solvers

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