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

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

Fix SOS problem and also bad bestSolution_

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