source: stable/2.9/Cbc/src/CbcModel.cpp @ 2143

Last change on this file since 2143 was 2106, checked in by forrest, 4 years ago

take out xxx ? 0 : 0

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