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

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

fix multiple solver conflict abort

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