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

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

make message handlers threadsafer (also modify multiple root solver)

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