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

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

fix stupid bug in orbital branching

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