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

Last change on this file since 1910 was 1910, checked in by tkr, 6 years ago

Merging r1909 from trunk

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