source: stable/2.8/Cbc/src/CbcModel.cpp @ 1942

Last change on this file since 1942 was 1926, checked in by stefan, 6 years ago

sync with trunk rev 1925

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