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

Last change on this file since 2362 was 2362, checked in by forrest, 3 years ago

put back switch off mini bab whenparallel

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