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

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

zerohalf thread safety plus restart with deterministic parallel

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