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

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

First try at orbital branching

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 744.5 KB
Line 
1/* $Id: CbcModel.cpp 2048 2014-07-16 09:29:16Z 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) {
2290        delete symmetryInfo_;
2291        symmetryInfo_=NULL;
2292        moreSpecialOptions2_ &= ~128;
2293      }
2294    }
2295#endif
2296
2297    // add cutoff as constraint if wanted
2298    if (cutoffRowNumber_==-2) {
2299      if (!parentModel_) {
2300        int numberColumns=solver_->getNumCols();
2301        double * obj = CoinCopyOfArray(solver_->getObjCoefficients(),numberColumns);
2302        int * indices = new int [numberColumns];
2303        int n=0;
2304        for (int i=0;i<numberColumns;i++) {
2305          if (obj[i]) {
2306            indices[n]=i;
2307            obj[n++]=obj[i];
2308          }
2309        }
2310        if (n) {
2311          double cutoff=getCutoff();
2312          // relax a little bit
2313          cutoff += 1.0e-4;
2314          double offset;
2315          solver_->getDblParam(OsiObjOffset, offset);
2316          cutoffRowNumber_ = solver_->getNumRows();
2317          solver_->addRow(n,indices,obj,-COIN_DBL_MAX,CoinMin(cutoff,1.0e25)+offset);
2318        } else {
2319          // no objective!
2320          cutoffRowNumber_ = -1;
2321        }
2322        delete [] indices;
2323        delete [] obj;
2324      } else {
2325        // switch off
2326        cutoffRowNumber_ = -1;
2327      }
2328    }
2329    numberRowsAtContinuous_ = getNumRows() ;
2330    solver_->saveBaseModel();
2331    /*
2332      Check the objective to see if we can deduce a nontrivial increment. If
2333      it's better than the current value for CbcCutoffIncrement, it'll be
2334      installed.
2335    */
2336    if (solverCharacteristics_->reducedCostsAccurate())
2337        analyzeObjective() ;
2338    {
2339        // may be able to change cutoff now
2340        double cutoff = getCutoff();
2341        double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
2342        if (cutoff > bestObjective_ - increment) {
2343            cutoff = bestObjective_ - increment ;
2344            setCutoff(cutoff) ;
2345        }
2346    }
2347#ifdef COIN_HAS_CLP
2348    // Possible save of pivot method
2349    ClpDualRowPivot * savePivotMethod = NULL;
2350    {
2351        // pass tolerance and increment to solver
2352        OsiClpSolverInterface * clpSolver
2353        = dynamic_cast<OsiClpSolverInterface *> (solver_);
2354        if (clpSolver)
2355            clpSolver->setStuff(getIntegerTolerance(), getCutoffIncrement());
2356#ifdef CLP_RESOLVE
2357        if ((moreSpecialOptions_&1048576)!=0&&!parentModel_&&clpSolver) {
2358          resolveClp(clpSolver,0);
2359        }
2360#endif
2361    }
2362#endif
2363    /*
2364      Set up for cut generation. addedCuts_ holds the cuts which are relevant for
2365      the active subproblem. whichGenerator will be used to record the generator
2366      that produced a given cut.
2367    */
2368    maximumWhich_ = 1000 ;
2369    delete [] whichGenerator_;
2370    whichGenerator_ = new int[maximumWhich_] ;
2371    memset(whichGenerator_, 0, maximumWhich_*sizeof(int));
2372    maximumNumberCuts_ = 0 ;
2373    currentNumberCuts_ = 0 ;
2374    delete [] addedCuts_ ;
2375    addedCuts_ = NULL ;
2376    OsiObject ** saveObjects = NULL;
2377    maximumRows_ = numberRowsAtContinuous_;
2378    currentDepth_ = 0;
2379    workingBasis_.resize(maximumRows_, numberColumns);
2380    /*
2381      Set up an empty heap and associated data structures to hold the live set
2382      (problems which require further exploration).
2383    */
2384    CbcCompareDefault * compareActual
2385    = dynamic_cast<CbcCompareDefault *> (nodeCompare_);
2386    if (compareActual) {
2387        compareActual->setBestPossible(direction*solver_->getObjValue());
2388        compareActual->setCutoff(getCutoff());
2389#ifdef JJF_ZERO
2390        if (false && !numberThreads_ && !parentModel_) {
2391            printf("CbcTreeArray ? threads ? parentArray\n");
2392            // Setup new style tree
2393            delete tree_;
2394            tree_ = new CbcTreeArray();
2395        }
2396#endif
2397    }
2398    tree_->setComparison(*nodeCompare_) ;
2399    /*
2400      Used to record the path from a node to the root of the search tree, so that
2401      we can then traverse from the root to the node when restoring a subproblem.
2402    */
2403    maximumDepth_ = 10 ;
2404    delete [] walkback_ ;
2405    walkback_ = new CbcNodeInfo * [maximumDepth_] ;
2406    lastDepth_ = 0;
2407    delete [] lastNodeInfo_ ;
2408    lastNodeInfo_ = new CbcNodeInfo * [maximumDepth_] ;
2409    delete [] lastNumberCuts_ ;
2410    lastNumberCuts_ = new int [maximumDepth_] ;
2411    maximumCuts_ = 100;
2412    lastNumberCuts2_ = 0;
2413    delete [] lastCut_;
2414    lastCut_ = new const OsiRowCut * [maximumCuts_];
2415    /*
2416      Used to generate bound edits for CbcPartialNodeInfo.
2417    */
2418    double * lowerBefore = new double [numberColumns] ;
2419    double * upperBefore = new double [numberColumns] ;
2420    /*
2421    Set up to run heuristics and generate cuts at the root node. The heavy
2422    lifting is hidden inside the calls to doHeuristicsAtRoot and solveWithCuts.
2423
2424    To start, tell cut generators they can be a bit more aggressive at the
2425    root node.
2426
2427    QUESTION: phase_ = 0 is documented as `initial solve', phase = 1 as `solve
2428        with cuts at root'. Is phase_ = 1 the correct indication when
2429        doHeurisiticsAtRoot is called to run heuristics outside of the main
2430        cut / heurisitc / reoptimise loop in solveWithCuts?
2431
2432      Generate cuts at the root node and reoptimise. solveWithCuts does the heavy
2433      lifting. It will iterate a generate/reoptimise loop (including reduced cost
2434      fixing) until no cuts are generated, the change in objective falls off,  or
2435      the limit on the number of rounds of cut generation is exceeded.
2436
2437      At the end of all this, any cuts will be recorded in cuts and also
2438      installed in the solver's constraint system. We'll have reoptimised, and
2439      removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
2440      adjusted accordingly).
2441
2442      Tell cut generators they can be a bit more aggressive at root node
2443
2444      TODO: Why don't we make a copy of the solution after solveWithCuts?
2445      TODO: If numberUnsatisfied == 0, don't we have a solution?
2446    */
2447    phase_ = 1;
2448    int iCutGenerator;
2449    for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
2450        // If parallel switch off global cuts
2451        if (numberThreads_) {
2452            generator_[iCutGenerator]->setGlobalCuts(false);
2453            generator_[iCutGenerator]->setGlobalCutsAtRoot(false);
2454        }
2455        CglCutGenerator * generator = generator_[iCutGenerator]->generator();
2456        generator->setAggressiveness(generator->getAggressiveness() + 100);
2457        if (!generator->canDoGlobalCuts())
2458          generator->setGlobalCuts(false);
2459    }
2460    OsiCuts cuts ;
2461    int anyAction = -1 ;
2462    numberOldActiveCuts_ = 0 ;
2463    numberNewCuts_ = 0 ;
2464    // Array to mark solution
2465    delete [] usedInSolution_;
2466    usedInSolution_ = new int[numberColumns];
2467    CoinZeroN(usedInSolution_, numberColumns);
2468    /*
2469      For printing totals and for CbcNode (numberNodes_)
2470    */
2471    numberIterations_ = 0 ;
2472    numberNodes_ = 0 ;
2473    numberNodes2_ = 0 ;
2474    maximumStatistics_ = 0;
2475    maximumDepthActual_ = 0;
2476    numberDJFixed_ = 0.0;
2477    if (!parentModel_) {
2478        if ((specialOptions_&262144) != 0) {
2479            // create empty stored cuts
2480            //storedRowCuts_ = new CglStored(solver_->getNumCols());
2481        } else if ((specialOptions_&524288) != 0 && storedRowCuts_) {
2482            // tighten and set best solution
2483            // A) tight bounds on integer variables
2484            /*
2485                storedRowCuts_ are coming in from outside, probably for nonlinear.
2486              John was unsure about origin.
2487            */
2488            const double * lower = solver_->getColLower();
2489            const double * upper = solver_->getColUpper();
2490            const double * tightLower = storedRowCuts_->tightLower();
2491            const double * tightUpper = storedRowCuts_->tightUpper();
2492            int nTightened = 0;
2493            for (int i = 0; i < numberIntegers_; i++) {
2494                int iColumn = integerVariable_[i];
2495                if (tightLower[iColumn] > lower[iColumn]) {
2496                    nTightened++;
2497                    solver_->setColLower(iColumn, tightLower[iColumn]);
2498                }
2499                if (tightUpper[iColumn] < upper[iColumn]) {
2500                    nTightened++;
2501                    solver_->setColUpper(iColumn, tightUpper[iColumn]);
2502                }
2503            }
2504            if (nTightened)
2505              COIN_DETAIL_PRINT(printf("%d tightened by alternate cuts\n", nTightened));
2506            if (storedRowCuts_->bestObjective() < bestObjective_) {
2507                // B) best solution
2508                double objValue = storedRowCuts_->bestObjective();
2509                setBestSolution(CBC_SOLUTION, objValue,
2510                                storedRowCuts_->bestSolution()) ;
2511                // Do heuristics
2512                // Allow RINS
2513                for (int i = 0; i < numberHeuristics_; i++) {
2514                    CbcHeuristicRINS * rins
2515                    = dynamic_cast<CbcHeuristicRINS *> (heuristic_[i]);
2516                    if (rins) {
2517                        rins->setLastNode(-100);
2518                    }
2519                }
2520            }
2521        }
2522    }
2523#ifdef SWITCH_VARIABLES
2524    // see if any switching variables
2525    if (numberIntegers_<solver_->getNumCols())
2526      findSwitching();
2527#endif
2528    /*
2529      Run heuristics at the root. This is the only opportunity to run FPump; it
2530      will be removed from the heuristics list by doHeuristicsAtRoot.
2531    */
2532    // See if multiple runs wanted
2533    CbcModel ** rootModels=NULL;
2534    if (!parentModel_&&multipleRootTries_%100) {
2535      double rootTimeCpu=CoinCpuTime();
2536      double startTimeRoot=CoinGetTimeOfDay();
2537      int numberRootThreads=1;
2538      /* undocumented fine tuning
2539         aabbcc where cc is number of tries
2540         bb if nonzero is number of threads
2541         aa if nonzero just do heuristics
2542      */
2543      int numberModels = multipleRootTries_%100;
2544#ifdef CBC_THREAD
2545      numberRootThreads = (multipleRootTries_/100)%100;
2546      if (!numberRootThreads)
2547        numberRootThreads=numberModels;
2548#endif
2549      int otherOptions = (multipleRootTries_/10000)%100;
2550      rootModels = new CbcModel * [numberModels];
2551      unsigned int newSeed = randomSeed_;
2552      if (newSeed==0) {
2553        double time = fabs(CoinGetTimeOfDay());
2554        while (time>=COIN_INT_MAX)
2555          time *= 0.5;
2556        newSeed = static_cast<unsigned int>(time);
2557      } else if (newSeed<0) {
2558        newSeed = 123456789;
2559#ifdef COIN_HAS_CLP
2560        OsiClpSolverInterface * clpSolver
2561          = dynamic_cast<OsiClpSolverInterface *> (solver_);
2562        if (clpSolver) {
2563          newSeed += clpSolver->getModelPtr()->randomNumberGenerator()->getSeed();
2564        }
2565#endif
2566      }
2567      CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis *> (solver_->getEmptyWarmStart());
2568      int numberIterations=0;
2569      for (int i=0;i<numberModels;i++) { 
2570        rootModels[i]=new CbcModel(*this);
2571        rootModels[i]->setNumberThreads(0);
2572        rootModels[i]->setMaximumNodes(otherOptions ? -1 : 0);
2573        rootModels[i]->setRandomSeed(newSeed+10000000*i);
2574        rootModels[i]->randomNumberGenerator()->setSeed(newSeed+50000000*i);
2575        rootModels[i]->setMultipleRootTries(0);
2576        // use seed
2577        rootModels[i]->setSpecialOptions(specialOptions_ |(4194304|8388608));
2578        rootModels[i]->setMoreSpecialOptions(moreSpecialOptions_ & 
2579                                             (~134217728));
2580        rootModels[i]->setMoreSpecialOptions2(moreSpecialOptions2_ & 
2581                                             (~128));
2582        rootModels[i]->solver_->setWarmStart(basis);
2583#ifdef COIN_HAS_CLP
2584        OsiClpSolverInterface * clpSolver
2585          = dynamic_cast<OsiClpSolverInterface *> (rootModels[i]->solver_);
2586        if (clpSolver) {
2587          ClpSimplex * simplex = clpSolver->getModelPtr();
2588          if (defaultHandler_)
2589            simplex->setDefaultMessageHandler();
2590          simplex->setRandomSeed(newSeed+20000000*i);
2591          simplex->allSlackBasis();
2592          int logLevel=simplex->logLevel();
2593          if (logLevel==1)
2594            simplex->setLogLevel(0);
2595          if (i!=0) {
2596            double random = simplex->randomNumberGenerator()->randomDouble();
2597            int bias = static_cast<int>(random*(numberIterations/4));
2598            simplex->setMaximumIterations(numberIterations/2+bias);
2599            simplex->primal();
2600            simplex->setMaximumIterations(COIN_INT_MAX);
2601            simplex->dual();
2602          } else {
2603            simplex->primal();
2604            numberIterations=simplex->numberIterations();
2605          }
2606          simplex->setLogLevel(logLevel);
2607          clpSolver->setWarmStart(NULL);
2608        }
2609#endif
2610        for (int j=0;j<numberHeuristics_;j++)
2611          rootModels[i]->heuristic_[j]->setSeed(rootModels[i]->heuristic_[j]->getSeed()+100000000*i);
2612        for (int j=0;j<numberCutGenerators_;j++)
2613          rootModels[i]->generator_[j]->generator()->refreshSolver(rootModels[i]->solver_);
2614      }
2615      delete basis;
2616#ifdef CBC_THREAD
2617      if (numberRootThreads==1) {
2618#endif
2619        for (int iModel=0;iModel<numberModels;iModel++) {
2620          doRootCbcThread(rootModels[iModel]);
2621          // see if solved at root node
2622          if (rootModels[iModel]->getMaximumNodes()) {
2623            feasible=false;
2624            break;
2625          }
2626        }
2627#ifdef CBC_THREAD
2628      } else {
2629        Coin_pthread_t * threadId = new Coin_pthread_t [numberRootThreads];
2630        for (int kModel=0;kModel<numberModels;kModel+=numberRootThreads) {
2631          bool finished=false;
2632          for (int iModel=kModel;iModel<CoinMin(numberModels,kModel+numberRootThreads);iModel++) {
2633            pthread_create(&(threadId[iModel-kModel].thr), NULL, 
2634                           doRootCbcThread,
2635                           rootModels[iModel]);
2636          }
2637          // wait
2638          for (int iModel=kModel;iModel<CoinMin(numberModels,kModel+numberRootThreads);iModel++) {
2639            pthread_join(threadId[iModel-kModel].thr, NULL);
2640          }
2641          // see if solved at root node
2642          for (int iModel=kModel;iModel<CoinMin(numberModels,kModel+numberRootThreads);iModel++) {
2643            if (rootModels[iModel]->getMaximumNodes())
2644              finished=true;
2645          }
2646          if (finished) {
2647            feasible=false;
2648            break;
2649          }
2650        }
2651        delete [] threadId;
2652      }
2653#endif
2654      // sort solutions
2655      int * which = new int [numberModels];
2656      double * value = new double [numberModels];
2657      int numberSolutions=0;
2658      for (int iModel=0;iModel<numberModels;iModel++) {
2659        if (rootModels[iModel]->bestSolution()) {
2660          which[numberSolutions]=iModel;
2661          value[numberSolutions++]=
2662            -rootModels[iModel]->getMinimizationObjValue();
2663        }
2664      }
2665      char general[100];
2666      rootTimeCpu=CoinCpuTime()-rootTimeCpu;
2667      if (numberRootThreads==1)
2668        sprintf(general,"Multiple root solvers took a total of %.2f seconds\n",
2669                rootTimeCpu);
2670      else
2671        sprintf(general,"Multiple root solvers took a total of %.2f seconds (%.2f elapsed)\n",
2672                rootTimeCpu,CoinGetTimeOfDay()-startTimeRoot);
2673      messageHandler()->message(CBC_GENERAL,
2674                                messages())
2675        << general << CoinMessageEol ;
2676      CoinSort_2(value,value+numberSolutions,which);
2677      // to get name
2678      CbcHeuristicRINS dummyHeuristic;
2679      dummyHeuristic.setHeuristicName("Multiple root solvers");
2680      lastHeuristic_=&dummyHeuristic;
2681      for (int i=0;i<numberSolutions;i++) {
2682        double objValue = - value[i];
2683        if (objValue<getCutoff()) {
2684          int iModel=which[i];
2685          setBestSolution(CBC_ROUNDING,objValue,
2686                          rootModels[iModel]->bestSolution());
2687        }
2688      }
2689      lastHeuristic_=NULL;
2690      delete [] which;
2691      delete [] value;
2692    }
2693    // Do heuristics
2694    if (numberObjects_&&!rootModels)
2695        doHeuristicsAtRoot();
2696    if (solverCharacteristics_->solutionAddsCuts()) {
2697      // With some heuristics solver needs a resolve here
2698      solver_->resolve();
2699      if(!isProvenOptimal()){
2700        solver_->initialSolve();
2701      }
2702    }
2703    /*
2704      Grepping through the code, it would appear that this is a command line
2705      debugging hook.  There's no obvious place in the code where this is set to
2706      a negative value.
2707
2708      User hook, says John.
2709    */
2710    if ( intParam_[CbcMaxNumNode] < 0
2711      ||numberSolutions_>=getMaximumSolutions())
2712      eventHappened_ = true; // stop as fast as possible
2713    stoppedOnGap_ = false ;
2714    // See if can stop on gap
2715    bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
2716    if(canStopOnGap()) {
2717        if (bestPossibleObjective_ < getCutoff())
2718            stoppedOnGap_ = true ;
2719        feasible = false;
2720        //eventHappened_=true; // stop as fast as possible
2721    }
2722    /*
2723      Set up for statistics collection, if requested. Standard values are
2724      documented in CbcModel.hpp. The magic number 100 will trigger a dump of
2725      CbcSimpleIntegerDynamicPseudoCost objects (no others). Looks like another
2726      command line debugging hook.
2727    */
2728    statistics_ = NULL;
2729    // Do on switch
2730    if (doStatistics > 0 && doStatistics <= 100) {
2731        maximumStatistics_ = 10000;
2732        statistics_ = new CbcStatistics * [maximumStatistics_];
2733        memset(statistics_, 0, maximumStatistics_*sizeof(CbcStatistics *));
2734    }
2735    // See if we can add integers
2736    if (noObjects && numberIntegers_ < solver_->getNumCols() && (specialOptions_&65536) != 0 && !parentModel_)
2737        AddIntegers();
2738    /*
2739      Do an initial round of cut generation for the root node. Depending on the
2740      type of underlying solver, we may want to do this even if the initial query
2741      to the objects indicates they're satisfied.
2742
2743      solveWithCuts does the heavy lifting. It will iterate a generate/reoptimise
2744      loop (including reduced cost fixing) until no cuts are generated, the
2745      change in objective falls off,  or the limit on the number of rounds of cut
2746      generation is exceeded.
2747
2748      At the end of all this, any cuts will be recorded in cuts and also
2749      installed in the solver's constraint system. We'll have reoptimised, and
2750      removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
2751      adjusted accordingly).
2752    */
2753    int iObject ;
2754    int preferredWay ;
2755    int numberUnsatisfied = 0 ;
2756    delete [] currentSolution_;
2757    currentSolution_ = new double [numberColumns];
2758    testSolution_ = currentSolution_;
2759    memcpy(currentSolution_, solver_->getColSolution(),
2760           numberColumns*sizeof(double)) ;
2761    // point to useful information
2762    OsiBranchingInformation usefulInfo = usefulInformation();
2763
2764    for (iObject = 0 ; iObject < numberObjects_ ; iObject++) {
2765        double infeasibility =
2766            object_[iObject]->infeasibility(&usefulInfo, preferredWay) ;
2767        if (infeasibility ) numberUnsatisfied++ ;
2768    }
2769    // replace solverType
2770    double * tightBounds = NULL;
2771    if (solverCharacteristics_->tryCuts())  {
2772
2773        if (numberUnsatisfied)   {
2774            // User event
2775            if (!eventHappened_ && feasible) {
2776                if (rootModels) {
2777                  // for fixings
2778                  int numberColumns=solver_->getNumCols();
2779                  tightBounds = new double [2*numberColumns];
2780                  {
2781                    const double * lower = solver_->getColLower();
2782                    const double * upper = solver_->getColUpper();
2783                    for (int i=0;i<numberColumns;i++) {
2784                      tightBounds[2*i+0]=lower[i];
2785                      tightBounds[2*i+1]=upper[i];
2786                    }
2787                  }
2788                  int numberModels = multipleRootTries_%100;
2789                  const OsiSolverInterface ** solvers = new 
2790                    const OsiSolverInterface * [numberModels];
2791                  int numberRows=continuousSolver_->getNumRows();
2792                  int maxCuts=0;
2793                  for (int i=0;i<numberModels;i++) {
2794                    solvers[i]=rootModels[i]->solver();
2795                    const double * lower = solvers[i]->getColLower();
2796                    const double * upper = solvers[i]->getColUpper();
2797                    for (int j=0;j<numberColumns;j++) {
2798                      tightBounds[2*j+0]=CoinMax(lower[j],tightBounds[2*j+0]);
2799                      tightBounds[2*j+1]=CoinMin(upper[j],tightBounds[2*j+1]);
2800                    }
2801                    int numberRows2=solvers[i]->getNumRows();
2802                    assert (numberRows2>=numberRows);
2803                    maxCuts += numberRows2-numberRows;
2804                    // accumulate statistics
2805                    for (int j=0;j<numberCutGenerators_;j++) {
2806                      generator_[j]->addStatistics(rootModels[i]->cutGenerator(j));
2807                    }
2808                  }
2809                  for (int j=0;j<numberCutGenerators_;j++) {
2810                    generator_[j]->scaleBackStatistics(numberModels);
2811                  }
2812                  //CbcRowCuts rowCut(maxCuts);
2813                  const OsiRowCutDebugger *debugger = NULL;
2814                  if ((specialOptions_&1) != 0) 
2815                    debugger = solver_->getRowCutDebugger() ;
2816                  for (int iModel=0;iModel<numberModels;iModel++) {
2817                    int numberRows2=solvers[iModel]->getNumRows();
2818                    const CoinPackedMatrix * rowCopy = solvers[iModel]->getMatrixByRow();
2819                    const int * rowLength = rowCopy->getVectorLengths(); 
2820                    const double * elements = rowCopy->getElements();
2821                    const int * column = rowCopy->getIndices();
2822                    const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
2823                    const double * rowLower = solvers[iModel]->getRowLower();
2824                    const double * rowUpper = solvers[iModel]->getRowUpper();
2825                    for (int iRow=numberRows;iRow<numberRows2;iRow++) {
2826                      OsiRowCut rc;
2827                      rc.setLb(rowLower[iRow]);
2828                      rc.setUb(rowUpper[iRow]);
2829                      CoinBigIndex start = rowStart[iRow];
2830                      rc.setRow(rowLength[iRow],column+start,elements+start,false);
2831                      if (debugger)
2832                        CoinAssert (!debugger->invalidCut(rc));
2833                      globalCuts_.addCutIfNotDuplicate(rc);
2834                    }
2835                    //int cutsAdded=globalCuts_.numberCuts()-numberCuts;
2836                    //numberCuts += cutsAdded;
2837                    //printf("Model %d gave %d cuts (out of %d possible)\n",
2838                    //     iModel,cutsAdded,numberRows2-numberRows);
2839                  }
2840                  // normally replace global cuts
2841                  //if (!globalCuts_.())
2842                  //globalCuts_=rowCutrowCut.addCuts(globalCuts_);
2843                  //rowCut.addCuts(globalCuts_);
2844                  int nTightened=0;
2845                  bool feasible=true;
2846                  {
2847                    double tolerance=1.0e-5;
2848                    const double * lower = solver_->getColLower();
2849                    const double * upper = solver_->getColUpper();
2850                    for (int i=0;i<numberColumns;i++) {
2851                      if (tightBounds[2*i+0]>tightBounds[2*i+1]) {
2852                        feasible=false;
2853                        printf("Bad bounds on %d %g,%g was %g,%g\n",
2854                               i,tightBounds[2*i+0],tightBounds[2*i+1],
2855                               lower[i],upper[i]);
2856                      }
2857                      //int k=0;
2858                      double oldLower=lower[i];
2859                      double oldUpper=upper[i];
2860                      if (tightBounds[2*i+0]>oldLower+tolerance) {
2861                        nTightened++;
2862                        //k++;
2863                        solver_->setColLower(i,tightBounds[2*i+0]);
2864                      }
2865                      if (tightBounds[2*i+1]<oldUpper-tolerance) {
2866                        nTightened++;
2867                        //k++;
2868                        solver_->setColUpper(i,tightBounds[2*i+1]);
2869                      }
2870                      //if (k)
2871                      //printf("new bounds on %d %g,%g was %g,%g\n",
2872                      //       i,tightBounds[2*i+0],tightBounds[2*i+1],
2873                      //       oldLower,oldUpper);
2874                    }
2875                    if (!feasible)
2876                      abort(); // deal with later
2877                  }
2878                  delete [] tightBounds;
2879                  tightBounds=NULL;
2880                  char printBuffer[200];
2881                  sprintf(printBuffer,"%d solvers added %d different cuts out of pool of %d",
2882                          numberModels,globalCuts_.sizeRowCuts(),maxCuts);
2883                  messageHandler()->message(CBC_GENERAL, messages())
2884                    << printBuffer << CoinMessageEol ;
2885                  if (nTightened) {
2886                    sprintf(printBuffer,"%d bounds were tightened",
2887                          nTightened);
2888                    messageHandler()->message(CBC_GENERAL, messages())
2889                      << printBuffer << CoinMessageEol ;
2890                  }
2891                  delete [] solvers;
2892                }
2893                if (!parentModel_&&(moreSpecialOptions_&67108864) != 0) {
2894                  // load cuts from file
2895                  FILE * fp = fopen("global.cuts","rb");
2896                  if (fp) {
2897                    size_t nRead;
2898                    int numberColumns=solver_->getNumCols();
2899                    int numCols;
2900                    nRead = fread(&numCols, sizeof(int), 1, fp);
2901                    if (nRead != 1)
2902                      throw("Error in fread");
2903                    if (numberColumns!=numCols) {
2904                      printf("Mismatch on columns %d %d\n",numberColumns,numCols);
2905                      fclose(fp);
2906                    } else {
2907                      // If rootModel just do some
2908                      double threshold=-1.0;
2909                      if (!multipleRootTries_)
2910                        threshold=0.5;
2911                      int initialCuts=0;
2912                      int initialGlobal = globalCuts_.sizeRowCuts();
2913                      double * elements = new double [numberColumns+2];
2914                      int * indices = new int [numberColumns];
2915                      int numberEntries=1;
2916                      while (numberEntries>0) {
2917                        nRead = fread(&numberEntries, sizeof(int), 1, fp);
2918                        if (nRead != 1)
2919                          throw("Error in fread");
2920                        double randomNumber=randomNumberGenerator_.randomDouble();
2921                        if (numberEntries>0) {
2922                          initialCuts++;
2923                          nRead = fread(elements, sizeof(double), numberEntries+2, fp);
2924                          if (nRead != static_cast<size_t>(numberEntries+2))
2925                            throw("Error in fread");
2926                          nRead = fread(indices, sizeof(int), numberEntries, fp);
2927                          if (nRead != static_cast<size_t>(numberEntries))
2928                            throw("Error in fread");
2929                          if (randomNumber>threshold) {
2930                            OsiRowCut rc;
2931                            rc.setLb(elements[numberEntries]);
2932                            rc.setUb(elements[numberEntries+1]);
2933                            rc.setRow(numberEntries,indices,elements,
2934                                      false);
2935                            rc.setGloballyValidAsInteger(2);
2936                            globalCuts_.addCutIfNotDuplicate(rc) ;
2937                          } 
2938                        }
2939                      }
2940                      fclose(fp);
2941                      // fixes
2942                      int nTightened=0;
2943                      fp = fopen("global.fix","rb");
2944                      if (fp) {
2945                        nRead = fread(indices, sizeof(int), 2, fp);
2946                        if (nRead != 2)
2947                          throw("Error in fread");
2948                        if (numberColumns!=indices[0]) {
2949                          printf("Mismatch on columns %d %d\n",numberColumns,
2950                                 indices[0]);
2951                        } else {
2952                          indices[0]=1;
2953                          while (indices[0]>=0) {
2954                            nRead = fread(indices, sizeof(int), 2, fp);
2955                            if (nRead != 2)
2956                              throw("Error in fread");
2957                            int iColumn=indices[0];
2958                            if (iColumn>=0) {
2959                              nTightened++;
2960                              nRead = fread(elements, sizeof(double), 4, fp);
2961                              if (nRead != 4)
2962                                throw("Error in fread");
2963                              solver_->setColLower(iColumn,elements[0]);
2964                              solver_->setColUpper(iColumn,elements[1]);
2965                            } 
2966                          }
2967                        }
2968                      }
2969                      if (fp)
2970                        fclose(fp);
2971                      char printBuffer[200];
2972                      sprintf(printBuffer,"%d cuts read in of which %d were unique, %d bounds tightened",
2973                             initialCuts,
2974                             globalCuts_.sizeRowCuts()-initialGlobal,nTightened); 
2975                      messageHandler()->message(CBC_GENERAL, messages())
2976                        << printBuffer << CoinMessageEol ;
2977                      delete [] elements;
2978                      delete [] indices;
2979                    }
2980                  }
2981                }
2982                feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
2983                                         NULL);
2984                if (multipleRootTries_&&
2985                    (moreSpecialOptions_&134217728)!=0) {
2986                  FILE * fp=NULL;
2987                  size_t nRead;
2988                  int numberColumns=solver_->getNumCols();
2989                  int initialCuts=0;
2990                  if ((moreSpecialOptions_&134217728)!=0) {
2991                    // append so go down to end
2992                    fp = fopen("global.cuts","r+b");
2993                    if (fp) {
2994                      int numCols;
2995                      nRead = fread(&numCols, sizeof(int), 1, fp);
2996                      if (nRead != 1)
2997                        throw("Error in fread");
2998                      if (numberColumns!=numCols) {
2999                        printf("Mismatch on columns %d %d\n",numberColumns,numCols);
3000                        fclose(fp);
3001                        fp=NULL;
3002                      }
3003                    }
3004                  }
3005                  double * elements = new double [numberColumns+2];
3006                  int * indices = new int [numberColumns];
3007                  if (fp) {
3008                    int numberEntries=1;
3009                    while (numberEntries>0) {
3010                      fpos_t position;
3011                      fgetpos(fp, &position);
3012                      nRead = fread(&numberEntries, sizeof(int), 1, fp);
3013                      if (nRead != 1)
3014                        throw("Error in fread");
3015                      if (numberEntries>0) {
3016                        initialCuts++;
3017                        nRead = fread(elements, sizeof(double), numberEntries+2, fp);
3018                        if (nRead != static_cast<size_t>(numberEntries+2))
3019                          throw("Error in fread");
3020                        nRead = fread(indices, sizeof(int), numberEntries, fp);
3021                        if (nRead != static_cast<size_t>(numberEntries))
3022                          throw("Error in fread");
3023                      } else {
3024                        // end
3025                        fsetpos(fp, &position);
3026                      }
3027                    }
3028                  } else {
3029                    fp = fopen("global.cuts","wb");
3030                    size_t nWrite;
3031                    nWrite=fwrite(&numberColumns,sizeof(int),1,fp);
3032                    if (nWrite != 1)
3033                      throw("Error in fwrite");
3034                  }
3035                  size_t nWrite;
3036                  // now append binding cuts
3037                  int numberC=continuousSolver_->getNumRows();
3038                  int numberRows=solver_->getNumRows();
3039                  printf("Saving %d cuts (up from %d)\n",
3040                         initialCuts+numberRows-numberC,initialCuts);
3041                  const double * rowLower = solver_->getRowLower();
3042                  const double * rowUpper = solver_->getRowUpper();
3043                  // Row copy
3044                  CoinPackedMatrix matrixByRow(*solver_->getMatrixByRow());
3045                  const double * elementByRow = matrixByRow.getElements();
3046                  const int * column = matrixByRow.getIndices();
3047                  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
3048                  const int * rowLength = matrixByRow.getVectorLengths();
3049                  for (int iRow=numberC;iRow<numberRows;iRow++) {
3050                    int n=rowLength[iRow];
3051                    assert (n);
3052                    CoinBigIndex start=rowStart[iRow];
3053                    memcpy(elements,elementByRow+start,n*sizeof(double));
3054                    memcpy(indices,column+start,n*sizeof(int));
3055                    elements[n]=rowLower[iRow];
3056                    elements[n+1]=rowUpper[iRow];
3057                    nWrite=fwrite(&n,sizeof(int),1,fp);
3058                    if (nWrite != 1)
3059                      throw("Error in fwrite");
3060                    nWrite=fwrite(elements,sizeof(double),n+2,fp);
3061                    if (nWrite != static_cast<size_t>(n+2))
3062                      throw("Error in fwrite");
3063                    nWrite=fwrite(indices,sizeof(int),n,fp);
3064                    if (nWrite != static_cast<size_t>(n))
3065                      throw("Error in fwrite");
3066                  }
3067                  // eof marker
3068                  int eofMarker=-1;
3069                  nWrite=fwrite(&eofMarker,sizeof(int),1,fp);
3070                  if (nWrite != 1)
3071                    throw("Error in fwrite");
3072                  fclose(fp);
3073                  // do tighter bounds (? later extra to original columns)
3074                  int nTightened=0;
3075                  const double * lower = solver_->getColLower();
3076                  const double * upper = solver_->getColUpper();
3077                  const double * originalLower = continuousSolver_->getColLower();
3078                  const double * originalUpper = continuousSolver_->getColUpper();
3079                  double tolerance=1.0e-5;
3080                  for (int i=0;i<numberColumns;i++) {
3081                    if (lower[i]>originalLower[i]+tolerance) {
3082                      nTightened++;
3083                    }
3084                    if (upper[i]<originalUpper[i]-tolerance) {
3085                      nTightened++;
3086                    }
3087                  }
3088                  if (nTightened) {
3089                    fp = fopen("global.fix","wb");
3090                    size_t nWrite;
3091                    indices[0]=numberColumns;
3092                    if (originalColumns_)
3093                      indices[1]=COIN_INT_MAX;
3094                    else
3095                      indices[1]=-1;
3096                    nWrite=fwrite(indices,sizeof(int),2,fp);
3097                    if (nWrite != 2)
3098                      throw("Error in fwrite");
3099                    for (int i=0;i<numberColumns;i++) {
3100                      int nTightened=0;
3101                      if (lower[i]>originalLower[i]+tolerance) {
3102                        nTightened++;
3103                      }
3104                      if (upper[i]<originalUpper[i]-tolerance) {
3105                        nTightened++;
3106                      }
3107                      if (nTightened) {
3108                        indices[0]=i;
3109                        if (originalColumns_)
3110                          indices[1]=originalColumns_[i];
3111                        elements[0]=lower[i];
3112                        elements[1]=upper[i];
3113                        elements[2]=originalLower[i];
3114                        elements[3]=originalUpper[i];
3115                        nWrite=fwrite(indices,sizeof(int),2,fp);
3116                        if (nWrite != 2)
3117                          throw("Error in fwrite");
3118                        nWrite=fwrite(elements,sizeof(double),4,fp);
3119                        if (nWrite != 4)
3120                          throw("Error in fwrite");
3121                      }
3122                    }
3123                    // eof marker
3124                    indices[0]=-1;
3125                    nWrite=fwrite(indices,sizeof(int),2,fp);
3126                    if (nWrite != 2)
3127                      throw("Error in fwrite");
3128                    fclose(fp);
3129                  }
3130                  delete [] elements;
3131                  delete [] indices;
3132                }
3133                if ((specialOptions_&524288) != 0 && !parentModel_
3134                        && storedRowCuts_) {
3135                    if (feasible) {
3136                        /* pick up stuff and try again
3137                        add cuts, maybe keep around
3138                        do best solution and if so new heuristics
3139                        obviously tighten bounds
3140                        */
3141                        // A and B probably done on entry
3142                        // A) tight bounds on integer variables
3143                        const double * lower = solver_->getColLower();
3144                        const double * upper = solver_->getColUpper();
3145                        const double * tightLower = storedRowCuts_->tightLower();
3146                        const double * tightUpper = storedRowCuts_->tightUpper();
3147                        int nTightened = 0;
3148                        for (int i = 0; i < numberIntegers_; i++) {
3149                            int iColumn = integerVariable_[i];
3150                            if (tightLower[iColumn] > lower[iColumn]) {
3151                                nTightened++;
3152                                solver_->setColLower(iColumn, tightLower[iColumn]);
3153                            }
3154                            if (tightUpper[iColumn] < upper[iColumn]) {
3155                                nTightened++;
3156                                solver_->setColUpper(iColumn, tightUpper[iColumn]);
3157                            }
3158                        }
3159                        if (nTightened)
3160                          COIN_DETAIL_PRINT(printf("%d tightened by alternate cuts\n", nTightened));
3161                        if (storedRowCuts_->bestObjective() < bestObjective_) {
3162                            // B) best solution
3163                            double objValue = storedRowCuts_->bestObjective();
3164                            setBestSolution(CBC_SOLUTION, objValue,
3165                                            storedRowCuts_->bestSolution()) ;
3166                            // Do heuristics
3167                            // Allow RINS
3168                            for (int i = 0; i < numberHeuristics_; i++) {
3169                                CbcHeuristicRINS * rins
3170                                = dynamic_cast<CbcHeuristicRINS *> (heuristic_[i]);
3171                                if (rins) {
3172                                    rins->setLastNode(-100);
3173                                }
3174                            }
3175                            doHeuristicsAtRoot();
3176                        }
3177#ifdef JJF_ZERO
3178                        int nCuts = storedRowCuts_->sizeRowCuts();
3179                        // add to global list
3180                        for (int i = 0; i < nCuts; i++) {
3181                            OsiRowCut newCut(*storedRowCuts_->rowCutPointer(i));
3182                            newCut.setGloballyValidAsInteger(2);
3183                            newCut.mutableRow().setTestForDuplicateIndex(false);
3184                            globalCuts_.insert(newCut) ;
3185                        }
3186#else
3187                        addCutGenerator(storedRowCuts_, -99, "Stored from previous run",
3188                                        true, false, false, -200);
3189#endif
3190                        // Set cuts as active
3191                        delete [] addedCuts_ ;
3192                        maximumNumberCuts_ = cuts.sizeRowCuts();
3193                        if (maximumNumberCuts_) {
3194                            addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
3195                        } else {
3196                            addedCuts_ = NULL;
3197                        }
3198                        for (int i = 0; i < maximumNumberCuts_; i++)
3199                            addedCuts_[i] = new CbcCountRowCut(*cuts.rowCutPtr(i),
3200                                                               NULL, -1, -1, 2);
3201                        COIN_DETAIL_PRINT(printf("size %d\n", cuts.sizeRowCuts()));
3202                        cuts = OsiCuts();
3203                        currentNumberCuts_ = maximumNumberCuts_;
3204                        feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
3205                                                 NULL);
3206                        for (int i = 0; i < maximumNumberCuts_; i++)
3207                            delete addedCuts_[i];
3208                    }
3209                    delete storedRowCuts_;
3210                    storedRowCuts_ = NULL;
3211                }
3212            } else {
3213                feasible = false;
3214            }
3215        }       else if (solverCharacteristics_->solutionAddsCuts() ||
3216                   solverCharacteristics_->alwaysTryCutsAtRootNode()) {
3217            // may generate cuts and turn the solution
3218            //to an infeasible one
3219            feasible = solveWithCuts(cuts, 2,
3220                                     NULL);
3221        }
3222    }
3223    if (rootModels) {
3224      int numberModels = multipleRootTries_%100;
3225      for (int i=0;i<numberModels;i++)
3226        delete rootModels[i];
3227      delete [] rootModels;
3228    }
3229    // check extra info on feasibility
3230    if (!solverCharacteristics_->mipFeasible())
3231        feasible = false;
3232    // If max nodes==0 - don't do strong branching
3233    if (!getMaximumNodes()) {
3234      if (feasible)
3235        feasible=false;
3236      else
3237        setMaximumNodes(1); //allow to stop on success
3238    }
3239    topOfTree_=NULL;
3240#ifdef CLP_RESOLVE
3241    if ((moreSpecialOptions_&2097152)!=0&&!parentModel_&&feasible) {
3242      OsiClpSolverInterface * clpSolver
3243        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3244      if (clpSolver)
3245        resolveClp(clpSolver,0);
3246    }
3247#endif
3248    // make cut generators less aggressive
3249    for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
3250        CglCutGenerator * generator = generator_[iCutGenerator]->generator();
3251        generator->setAggressiveness(generator->getAggressiveness() - 100);
3252    }
3253    currentNumberCuts_ = numberNewCuts_ ;
3254    if (solverCharacteristics_->solutionAddsCuts()) {
3255      // With some heuristics solver needs a resolve here (don't know if this is bug in heuristics)
3256      solver_->resolve();
3257      if(!isProvenOptimal()){
3258        solver_->initialSolve();
3259      }
3260    }
3261    // See if can stop on gap
3262    bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
3263    if(canStopOnGap()) {
3264        if (bestPossibleObjective_ < getCutoff())
3265            stoppedOnGap_ = true ;
3266        feasible = false;
3267    }
3268    // User event
3269    if (eventHappened_)
3270        feasible = false;
3271#if defined(COIN_HAS_CLP)&&defined(COIN_HAS_CPX)
3272    /*
3273      This is the notion of using Cbc stuff to get going, then calling cplex to
3274      finish off.
3275    */
3276    if (feasible && (specialOptions_&16384) != 0 && fastNodeDepth_ == -2 && !parentModel_) {
3277        // Use Cplex to do search!
3278        double time1 = CoinCpuTime();
3279        OsiClpSolverInterface * clpSolver
3280        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3281        OsiCpxSolverInterface cpxSolver;
3282        double direction = clpSolver->getObjSense();
3283        cpxSolver.setObjSense(direction);
3284        // load up cplex
3285        const CoinPackedMatrix * matrix = continuousSolver_->getMatrixByCol();
3286        const double * rowLower = continuousSolver_->getRowLower();
3287        const double * rowUpper = continuousSolver_->getRowUpper();
3288        const double * columnLower = continuousSolver_->getColLower();
3289        const double * columnUpper = continuousSolver_->getColUpper();
3290        const double * objective = continuousSolver_->getObjCoefficients();
3291        cpxSolver.loadProblem(*matrix, columnLower, columnUpper,
3292                              objective, rowLower, rowUpper);
3293        double * setSol = new double [numberIntegers_];
3294        int * setVar = new int [numberIntegers_];
3295        // cplex doesn't know about objective offset
3296        double offset = clpSolver->getModelPtr()->objectiveOffset();
3297        for (int i = 0; i < numberIntegers_; i++) {
3298            int iColumn = integerVariable_[i];
3299            cpxSolver.setInteger(iColumn);
3300            if (bestSolution_) {
3301                setSol[i] = bestSolution_[iColumn];
3302                setVar[i] = iColumn;
3303            }
3304        }
3305        CPXENVptr env = cpxSolver.getEnvironmentPtr();
3306        CPXLPptr lpPtr = cpxSolver.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL);
3307        cpxSolver.switchToMIP();
3308        if (bestSolution_) {
3309            CPXcopymipstart(env, lpPtr, numberIntegers_, setVar, setSol);
3310        }
3311        if (clpSolver->getNumRows() > continuousSolver_->getNumRows() && false) {
3312            // add cuts
3313            const CoinPackedMatrix * matrix = clpSolver->getMatrixByRow();
3314            const double * rhs = clpSolver->getRightHandSide();
3315            const char * rowSense = clpSolver->getRowSense();
3316            const double * elementByRow = matrix->getElements();
3317            const int * column = matrix->getIndices();
3318            const CoinBigIndex * rowStart = matrix->getVectorStarts();
3319            const int * rowLength = matrix->getVectorLengths();
3320            int nStart = continuousSolver_->getNumRows();
3321            int nRows = clpSolver->getNumRows();
3322            int size = rowStart[nRows-1] + rowLength[nRows-1] -
3323                       rowStart[nStart];
3324            int nAdd = 0;
3325            double * rmatval = new double [size];
3326            int * rmatind = new int [size];
3327            int * rmatbeg = new int [nRows-nStart+1];
3328            size = 0;
3329            rmatbeg[0] = 0;
3330            for (int i = nStart; i < nRows; i++) {
3331                for (int k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
3332                    rmatind[size] = column[k];
3333                    rmatval[size++] = elementByRow[k];
3334                }
3335                nAdd++;
3336                rmatbeg[nAdd] = size;
3337            }
3338            CPXaddlazyconstraints(env, lpPtr, nAdd, size,
3339                                  rhs, rowSense, rmatbeg,
3340                                  rmatind, rmatval, NULL);
3341            CPXsetintparam( env, CPX_PARAM_REDUCE,
3342                            // CPX_PREREDUCE_NOPRIMALORDUAL (0)
3343                            CPX_PREREDUCE_PRIMALONLY);
3344        }
3345        if (getCutoff() < 1.0e50) {
3346            double useCutoff = getCutoff() + offset;
3347            if (bestObjective_ < 1.0e50)
3348                useCutoff = bestObjective_ + offset + 1.0e-7;
3349            cpxSolver.setDblParam(OsiDualObjectiveLimit, useCutoff*
3350                                  direction);
3351            if ( direction > 0.0 )
3352                CPXsetdblparam( env, CPX_PARAM_CUTUP, useCutoff ) ; // min
3353            else
3354                CPXsetdblparam( env, CPX_PARAM_CUTLO, useCutoff ) ; // max
3355        }
3356        CPXsetdblparam(env, CPX_PARAM_EPGAP, dblParam_[CbcAllowableFractionGap]);
3357        delete [] setSol;
3358        delete [] setVar;
3359        char printBuffer[200];
3360        if (offset) {
3361            sprintf(printBuffer, "Add %g to all Cplex messages for true objective",
3362                    -offset);
3363            messageHandler()->message(CBC_GENERAL, messages())
3364            << printBuffer << CoinMessageEol ;
3365            cpxSolver.setDblParam(OsiObjOffset, offset);
3366        }
3367        cpxSolver.branchAndBound();
3368        double timeTaken = CoinCpuTime() - time1;
3369        sprintf(printBuffer, "Cplex took %g seconds",
3370                timeTaken);
3371        messageHandler()->message(CBC_GENERAL, messages())
3372        << printBuffer << CoinMessageEol ;
3373        numberExtraNodes_ = CPXgetnodecnt(env, lpPtr);
3374        numberExtraIterations_ = CPXgetmipitcnt(env, lpPtr);
3375        double value = cpxSolver.getObjValue() * direction;
3376        if (cpxSolver.isProvenOptimal() && value <= getCutoff()) {
3377            feasible = true;
3378            clpSolver->setWarmStart(NULL);
3379            // try and do solution
3380            double * newSolution =
3381                CoinCopyOfArray(cpxSolver.getColSolution(),
3382                                getNumCols());
3383            setBestSolution(CBC_STRONGSOL, value, newSolution) ;
3384            delete [] newSolution;
3385        }
3386        feasible = false;
3387    }
3388#endif
3389    if (!parentModel_&&(moreSpecialOptions_&268435456) != 0) {
3390      // try idiotic idea
3391      CbcObject * obj = new CbcIdiotBranch(this);
3392      obj->setPriority(1); // temp
3393      addObjects(1, &obj);
3394      delete obj;
3395    }
3396   
3397    /*
3398      A hook to use clp to quickly explore some part of the tree.
3399    */
3400    if (fastNodeDepth_ == 1000 &&/*!parentModel_*/(specialOptions_&2048) == 0) {
3401        fastNodeDepth_ = -1;
3402        CbcObject * obj =
3403            new CbcFollowOn(this);
3404        obj->setPriority(1);
3405        addObjects(1, &obj);
3406        delete obj;
3407    }
3408    int saveNumberSolves = numberSolves_;
3409    int saveNumberIterations = numberIterations_;
3410    if ((fastNodeDepth_ >= 0||(moreSpecialOptions_&33554432)!=0)
3411        &&/*!parentModel_*/(specialOptions_&2048) == 0) {
3412        // add in a general depth object doClp
3413        int type = (fastNodeDepth_ <= 100) ? fastNodeDepth_ : -(fastNodeDepth_ - 100);
3414        if ((moreSpecialOptions_&33554432)!=0)
3415          type=12;
3416        else
3417          fastNodeDepth_ += 1000000;     // mark as done
3418        CbcObject * obj =
3419            new CbcGeneralDepth(this, type);
3420        addObjects(1, &obj);
3421        delete obj;
3422        // fake number of objects
3423        numberObjects_--;
3424        if (parallelMode() < -1) {
3425            // But make sure position is correct
3426            OsiObject * obj2 = object_[numberObjects_];
3427            obj = dynamic_cast<CbcObject *> (obj2);
3428            assert (obj);
3429            obj->setPosition(numberObjects_);
3430        }
3431    }
3432#ifdef COIN_HAS_CLP
3433#ifdef NO_CRUNCH
3434    if (true) {
3435        OsiClpSolverInterface * clpSolver
3436        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3437        if (clpSolver && !parentModel_) {
3438            ClpSimplex * clpSimplex = clpSolver->getModelPtr();
3439            clpSimplex->setSpecialOptions(clpSimplex->specialOptions() | 131072);
3440            //clpSimplex->startPermanentArrays();
3441            clpSimplex->setPersistenceFlag(2);
3442        }
3443    }
3444#endif
3445#endif
3446// Save copy of solver
3447    OsiSolverInterface * saveSolver = NULL;
3448    if (!parentModel_ && (specialOptions_&(512 + 32768)) != 0)
3449        saveSolver = solver_->clone();
3450    double checkCutoffForRestart = 1.0e100;
3451    saveModel(saveSolver, &checkCutoffForRestart, &feasible);
3452    if ((specialOptions_&262144) != 0 && !parentModel_) {
3453        // Save stuff and return!
3454        storedRowCuts_->saveStuff(bestObjective_, bestSolution_,
3455                                  solver_->getColLower(),
3456                                  solver_->getColUpper());
3457        delete [] lowerBefore;
3458        delete [] upperBefore;
3459        delete saveSolver;
3460        if (flipObjective)
3461          flipModel();
3462        return;
3463    }
3464    /*
3465      We've taken the continuous relaxation as far as we can. Time to branch.
3466      The first order of business is to actually create a node. chooseBranch
3467      currently uses strong branching to evaluate branch object candidates,
3468      unless forced back to simple branching. If chooseBranch concludes that a
3469      branching candidate is monotone (anyAction == -1) or infeasible (anyAction
3470      == -2) when forced to integer values, it returns here immediately.
3471
3472      Monotone variables trigger a call to resolve(). If the problem remains
3473      feasible, try again to choose a branching variable. At the end of the loop,
3474      resolved == true indicates that some variables were fixed.
3475
3476      Loss of feasibility will result in the deletion of newNode.
3477    */
3478
3479    bool resolved = false ;
3480    CbcNode *newNode = NULL ;
3481    numberFixedAtRoot_ = 0;
3482    numberFixedNow_ = 0;
3483    if (!parentModel_&&(moreSpecialOptions2_&2)!=0) {
3484#ifdef COIN_HAS_CLP
3485      OsiClpSolverInterface * clpSolver
3486        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3487      if (clpSolver) {
3488        if (getCutoff()>1.0e20) {
3489          printf("Zapping costs\n");
3490          int numberColumns=solver_->getNumCols();
3491          double * zeroCost = new double [numberColumns];
3492          // could make small random
3493          memset(zeroCost,0,numberColumns*sizeof(double));
3494          solver_->setObjective(zeroCost);
3495          double objValue = solver_->getObjValue();
3496          solver_->setDblParam(OsiObjOffset,-objValue);
3497          clpSolver->getModelPtr()->setObjectiveValue(objValue);
3498          delete [] zeroCost;
3499        } else {
3500          moreSpecialOptions2_ &= ~2;
3501        }
3502      } else {
3503#endif
3504          moreSpecialOptions2_ &= ~2;
3505#ifdef COIN_HAS_CLP
3506      }
3507#endif
3508    }
3509    int numberIterationsAtContinuous = numberIterations_;
3510    //solverCharacteristics_->setSolver(solver_);
3511    if (feasible) {
3512#define HOTSTART -1
3513#if HOTSTART<0
3514        if (bestSolution_ && !parentModel_ && !hotstartSolution_ &&
3515                (moreSpecialOptions_&1024) != 0) {
3516            // Set priorities so only branch on ones we need to
3517            // use djs and see if only few branches needed
3518#ifndef NDEBUG
3519            double integerTolerance = getIntegerTolerance() ;
3520#endif
3521            bool possible = true;
3522            const double * saveLower = continuousSolver_->getColLower();
3523            const double * saveUpper = continuousSolver_->getColUpper();
3524            for (int i = 0; i < numberObjects_; i++) {
3525                const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object_[i]);
3526                if (thisOne) {
3527                    int iColumn = thisOne->columnNumber();
3528                    if (saveUpper[iColumn] > saveLower[iColumn] + 1.5) {
3529                        possible = false;
3530                        break;
3531                    }
3532                } else {
3533                    possible = false;
3534                    break;
3535                }
3536            }
3537            if (possible) {
3538                OsiSolverInterface * solver = continuousSolver_->clone();
3539                int numberColumns = solver->getNumCols();
3540                for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
3541                    double value = bestSolution_[iColumn] ;
3542                    value = CoinMax(value, saveLower[iColumn]) ;
3543                    value = CoinMin(value, saveUpper[iColumn]) ;
3544                    value = floor(value + 0.5);
3545                    if (solver->isInteger(iColumn)) {
3546                        solver->setColLower(iColumn, value);
3547                        solver->setColUpper(iColumn, value);
3548                    }
3549                }
3550                solver->setHintParam(OsiDoDualInResolve, false, OsiHintTry);
3551                // objlim and all slack
3552                double direction = solver->getObjSense();
3553                solver->setDblParam(OsiDualObjectiveLimit, 1.0e50*direction);
3554                CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis *> (solver->getEmptyWarmStart());
3555                solver->setWarmStart(basis);
3556                delete basis;
3557                bool changed = true;
3558                hotstartPriorities_ = new int [numberColumns];
3559                for (int iColumn = 0; iColumn < numberColumns; iColumn++)
3560                    hotstartPriorities_[iColumn] = 1;
3561                while (changed) {
3562                    changed = false;
3563                    solver->resolve();
3564                    if (!solver->isProvenOptimal()) {
3565                        possible = false;
3566                        break;
3567                    }
3568                    const double * dj = solver->getReducedCost();
3569                    const double * colLower = solver->getColLower();
3570                    const double * colUpper = solver->getColUpper();
3571                    const double * solution = solver->getColSolution();
3572                    int nAtLbNatural = 0;
3573                    int nAtUbNatural = 0;
3574                    int nZeroDj = 0;
3575                    int nForced = 0;
3576                    for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
3577                        double value = solution[iColumn] ;
3578                        value = CoinMax(value, saveLower[iColumn]) ;
3579                        value = CoinMin(value, saveUpper[iColumn]) ;
3580                        if (solver->isInteger(iColumn)) {
3581                            assert(fabs(value - solution[iColumn]) <= integerTolerance) ;
3582                            if (hotstartPriorities_[iColumn] == 1) {
3583                                if (dj[iColumn] < -1.0e-6) {
3584                                    // negative dj
3585                                    if (saveUpper[iColumn] == colUpper[iColumn]) {
3586                                        nAtUbNatural++;
3587                                        hotstartPriorities_[iColumn] = 2;
3588                                        solver->setColLower(iColumn, saveLower[iColumn]);
3589                                        solver->setColUpper(iColumn, saveUpper[iColumn]);
3590                                    } else {
3591                                        nForced++;
3592                                    }
3593                                } else if (dj[iColumn] > 1.0e-6) {
3594                                    // positive dj
3595                                    if (saveLower[iColumn] == colLower[iColumn]) {
3596                                        nAtLbNatural++;
3597                                        hotstartPriorities_[iColumn] = 2;
3598                                        solver->setColLower(iColumn, saveLower[iColumn]);
3599                                        solver->setColUpper(iColumn, saveUpper[iColumn]);
3600                                    } else {
3601                                        nForced++;
3602                                    }
3603                                } else {
3604                                    // zero dj
3605                                    nZeroDj++;
3606                                }
3607                            }
3608                        }
3609                    }
3610#ifdef CLP_INVESTIGATE
3611                    printf("%d forced, %d naturally at lower, %d at upper - %d zero dj\n",
3612                           nForced, nAtLbNatural, nAtUbNatural, nZeroDj);
3613#endif
3614                    if (nAtLbNatural || nAtUbNatural) {
3615                        changed = true;
3616                    } else {
3617                        if (nForced + nZeroDj > 5000 ||
3618                                (nForced + nZeroDj)*2 > numberIntegers_)
3619                            possible = false;
3620                    }
3621                }
3622                delete solver;
3623            }
3624            if (possible) {
3625                setHotstartSolution(bestSolution_);
3626                if (!saveCompare) {
3627                    // create depth first comparison
3628                    saveCompare = nodeCompare_;
3629                    // depth first
3630                    nodeCompare_ = new CbcCompareDepth();
3631                    tree_->setComparison(*nodeCompare_) ;
3632                }
3633            } else {
3634                delete [] hotstartPriorities_;
3635                hotstartPriorities_ = NULL;
3636            }
3637        }
3638#endif
3639#if HOTSTART>0
3640        if (hotstartSolution_ && !hotstartPriorities_) {
3641            // Set up hot start
3642            OsiSolverInterface * solver = solver_->clone();
3643            double direction = solver_->getObjSense() ;
3644            int numberColumns = solver->getNumCols();
3645            double * saveLower = CoinCopyOfArray(solver->getColLower(), numberColumns);
3646            double * saveUpper = CoinCopyOfArray(solver->getColUpper(), numberColumns);
3647            // move solution
3648            solver->setColSolution(hotstartSolution_);
3649            // point to useful information
3650            const double * saveSolution = testSolution_;
3651            testSolution_ = solver->getColSolution();
3652            OsiBranchingInformation usefulInfo = usefulInformation();
3653            testSolution_ = saveSolution;
3654            /*
3655            Run through the objects and use feasibleRegion() to set variable bounds
3656            so as to fix the variables specified in the objects at their value in this
3657            solution. Since the object list contains (at least) one object for every
3658            integer variable, this has the effect of fixing all integer variables.
3659            */
3660            for (int i = 0; i < numberObjects_; i++)
3661                object_[i]->feasibleRegion(solver, &usefulInfo);
3662            solver->resolve();
3663            assert (solver->isProvenOptimal());
3664            double gap = CoinMax((solver->getObjValue() - solver_->getObjValue()) * direction, 0.0) ;
3665            const double * dj = solver->getReducedCost();
3666            const double * colLower = solver->getColLower();
3667            const double * colUpper = solver->getColUpper();
3668            const double * solution = solver->getColSolution();
3669            int nAtLbNatural = 0;
3670            int nAtUbNatural = 0;
3671            int nAtLbNaturalZero = 0;
3672            int nAtUbNaturalZero = 0;
3673            int nAtLbFixed = 0;
3674            int nAtUbFixed = 0;
3675            int nAtOther = 0;
3676            int nAtOtherNatural = 0;
3677            int nNotNeeded = 0;
3678            delete [] hotstartSolution_;
3679            hotstartSolution_ = new double [numberColumns];
3680            delete [] hotstartPriorities_;
3681            hotstartPriorities_ = new int [numberColumns];
3682            int * order = (int *) saveUpper;
3683            int nFix = 0;
3684            double bestRatio = COIN_DBL_MAX;
3685            for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
3686                double value = solution[iColumn] ;
3687                value = CoinMax(value, saveLower[iColumn]) ;
3688                value = CoinMin(value, saveUpper[iColumn]) ;
3689                double sortValue = COIN_DBL_MAX;
3690                if (solver->isInteger(iColumn)) {
3691                    assert(fabs(value - solution[iColumn]) <= 1.0e-5) ;
3692                    double value2 = floor(value + 0.5);
3693                    if (dj[iColumn] < -1.0e-6) {
3694                        // negative dj
3695                        //assert (value2==colUpper[iColumn]);
3696                        if (saveUpper[iColumn] == colUpper[iColumn]) {
3697                            nAtUbNatural++;
3698                            sortValue = 0.0;
3699                            double value = -dj[iColumn];
3700                            if (value > gap)
3701                                nFix++;
3702                            else if (gap < value*bestRatio)
3703                                bestRatio = gap / value;
3704                            if (saveLower[iColumn] != colLower[iColumn]) {
3705                                nNotNeeded++;
3706                                sortValue = 1.0e20;
3707                            }
3708                        } else if (saveLower[iColumn] == colUpper[iColumn]) {
3709                            nAtLbFixed++;
3710                            sortValue = dj[iColumn];
3711                        } else {
3712                            nAtOther++;
3713                            sortValue = 0.0;
3714                            if (saveLower[iColumn] != colLower[iColumn] &&
3715                                    saveUpper[iColumn] != colUpper[iColumn]) {
3716                                nNotNeeded++;
3717                                sortValue = 1.0e20;
3718                            }
3719                        }
3720                    } else if (dj[iColumn] > 1.0e-6) {
3721                        // positive dj
3722                        //assert (value2==colLower[iColumn]);
3723                        if (saveLower[iColumn] == colLower[iColumn]) {
3724                            nAtLbNatural++;
3725                            sortValue = 0.0;
3726                            double value = dj[iColumn];
3727                            if (value > gap)
3728                                nFix++;
3729                            else if (gap < value*bestRatio)
3730                                bestRatio = gap / value;
3731                            if (saveUpper[iColumn] != colUpper[iColumn]) {
3732                                nNotNeeded++;
3733                                sortValue = 1.0e20;
3734                            }
3735                        } else if (saveUpper[iColumn] == colLower[iColumn]) {
3736                            nAtUbFixed++;
3737                            sortValue = -dj[iColumn];
3738                        } else {
3739                            nAtOther++;
3740                            sortValue = 0.0;
3741                            if (saveLower[iColumn] != colLower[iColumn] &&
3742                                    saveUpper[iColumn] != colUpper[iColumn]) {
3743                                nNotNeeded++;
3744                                sortValue = 1.0e20;
3745                            }
3746                        }
3747                    } else {
3748                        // zero dj
3749                        if (value2 == saveUpper[iColumn]) {
3750                            nAtUbNaturalZero++;
3751                            sortValue = 0.0;
3752                            if (saveLower[iColumn] != colLower[iColumn]) {
3753                                nNotNeeded++;
3754                                sortValue = 1.0e20;
3755                            }
3756                        } else if (value2 == saveLower[iColumn]) {
3757                            nAtLbNaturalZero++;
3758                            sortValue = 0.0;
3759                        } else {
3760                            nAtOtherNatural++;
3761                            sortValue = 0.0;
3762                            if (saveLower[iColumn] != colLower[iColumn] &&
3763                                    saveUpper[iColumn] != colUpper[iColumn]) {
3764                                nNotNeeded++;
3765                                sortValue = 1.0e20;
3766                            }
3767                        }
3768                    }
3769#if HOTSTART==3
3770                    sortValue = -fabs(dj[iColumn]);
3771#endif
3772                }
3773                hotstartSolution_[iColumn] = value ;
3774                saveLower[iColumn] = sortValue;
3775                order[iColumn] = iColumn;
3776            }
3777            COIN_DETAIL_PRINT(printf("** can fix %d columns - best ratio for others is %g on gap of %g\n",
3778                                     nFix, bestRatio, gap));
3779            int nNeg = 0;
3780            CoinSort_2(saveLower, saveLower + numberColumns, order);
3781            for (int i = 0; i < numberColumns; i++) {
3782                if (saveLower[i] < 0.0) {
3783                    nNeg++;
3784#if HOTSTART==2||HOTSTART==3
3785                    // swap sign ?
3786                    saveLower[i] = -saveLower[i];
3787#endif
3788                }
3789            }
3790            CoinSort_2(saveLower, saveLower + nNeg, order);
3791            for (int i = 0; i < numberColumns; i++) {
3792#if HOTSTART==1
3793                hotstartPriorities_[order[i]] = 100;
3794#else
3795                hotstartPriorities_[order[i]] = -(i + 1);
3796#endif
3797            }
3798            COIN_DETAIL_PRINT(printf("nAtLbNat %d,nAtUbNat %d,nAtLbNatZero %d,nAtUbNatZero %d,nAtLbFixed %d,nAtUbFixed %d,nAtOther %d,nAtOtherNat %d, useless %d %d\n",
3799                   nAtLbNatural,
3800                   nAtUbNatural,
3801                   nAtLbNaturalZero,
3802                   nAtUbNaturalZero,
3803                   nAtLbFixed,
3804                   nAtUbFixed,
3805                   nAtOther,
3806                                     nAtOtherNatural, nNotNeeded, nNeg));
3807            delete [] saveLower;
3808            delete [] saveUpper;
3809            if (!saveCompare) {
3810                // create depth first comparison
3811                saveCompare = nodeCompare_;
3812                // depth first
3813                nodeCompare_ = new CbcCompareDepth();
3814                tree_->setComparison(*nodeCompare_) ;
3815            }
3816        }
3817#endif
3818        newNode = new CbcNode ;
3819        // Set objective value (not so obvious if NLP etc)
3820        setObjectiveValue(newNode, NULL);
3821        anyAction = -1 ;
3822        // To make depth available we may need a fake node
3823        CbcNode fakeNode;
3824        if (!currentNode_) {
3825            // Not true if sub trees assert (!numberNodes_);
3826            currentNode_ = &fakeNode;
3827        }
3828        phase_ = 3;
3829        // only allow 1000 passes
3830        int numberPassesLeft = 1000;
3831        // This is first crude step
3832        if (numberAnalyzeIterations_) {
3833            delete [] analyzeResults_;
3834            analyzeResults_ = new double [4*numberIntegers_];
3835            numberFixedAtRoot_ = newNode->analyze(this, analyzeResults_);
3836            if (numberFixedAtRoot_ > 0) {
3837              COIN_DETAIL_PRINT(printf("%d fixed by analysis\n", numberFixedAtRoot_));
3838                setPointers(solver_);
3839                numberFixedNow_ = numberFixedAtRoot_;
3840            } else if (numberFixedAtRoot_ < 0) {
3841              COIN_DETAIL_PRINT(printf("analysis found to be infeasible\n"));
3842                anyAction = -2;
3843                delete newNode ;
3844                newNode = NULL ;
3845                feasible = false ;
3846            }
3847        }
3848        OsiSolverBranch * branches = NULL;
3849        anyAction = chooseBranch(newNode, numberPassesLeft, NULL, cuts, resolved,
3850                                 NULL, NULL, NULL, branches);
3851        if (anyAction == -2 || newNode->objectiveValue() >= cutoff) {
3852            if (anyAction != -2) {
3853                // zap parent nodeInfo
3854#ifdef COIN_DEVELOP
3855                printf("zapping CbcNodeInfo %x\n", reinterpret_cast<int>(newNode->nodeInfo()->parent()));
3856#endif
3857                if (newNode->nodeInfo())
3858                    newNode->nodeInfo()->nullParent();
3859            }
3860            delete newNode ;
3861            newNode = NULL ;
3862            feasible = false ;
3863        }
3864    }
3865    if (newNode && probingInfo_) {
3866        int number01 = probingInfo_->numberIntegers();
3867        //const fixEntry * entry = probingInfo_->fixEntries();
3868        const int * toZero = probingInfo_->toZero();
3869        //const int * toOne = probingInfo_->toOne();
3870        //const int * integerVariable = probingInfo_->integerVariable();
3871        if (toZero[number01]) {
3872            CglTreeProbingInfo info(*probingInfo_);
3873            if ((moreSpecialOptions2_&64)!=0&&!parentModel_) {
3874              /*
3875                Marginal idea. Further exploration probably good. Build some extra
3876                cliques from probing info. Not quite worth the effort?
3877              */
3878              CglProbing generator1;
3879              generator1.setUsingObjective(false);
3880              generator1.setMaxPass(1);
3881              generator1.setMaxPassRoot(1);
3882              generator1.setMaxLook(100);
3883              generator1.setRowCuts(3);
3884              generator1.setMaxElements(300);
3885              generator1.setMaxProbeRoot(solver_->getNumCols());
3886              CoinThreadRandom randomGenerator;
3887              //CglTreeProbingInfo info(solver_);
3888              info.level = 0;
3889              info.formulation_rows = solver_->getNumRows();
3890              info.inTree = false;
3891              info.randomNumberGenerator=&randomGenerator;
3892              info.pass=4;
3893              generator1.setMode(8);
3894              OsiCuts cs;
3895              generator1.generateCutsAndModify(*solver_,cs,&info);
3896              // very clunky
3897              OsiSolverInterface * temp = generator1.cliqueModel(solver_,2);
3898              CglPreProcess dummy;
3899              OsiSolverInterface * newSolver=dummy.cliqueIt(*temp,0.0001);
3900              delete temp;
3901              OsiSolverInterface * fake = NULL;
3902              if (newSolver) {
3903#if 0
3904                int numberCliques = generator1.numberCliques();
3905                cliqueEntry * entry = generator1.cliqueEntry();
3906                cliqueType * type = new cliqueType [numberCliques];
3907                int * start = new int [numberCliques+1];
3908                start[numberCliques]=2*numberCliques;
3909                int n=0;
3910                for (int i=0;i<numberCliques;i++) {
3911                  start[i]=2*i;
3912                  setOneFixesInCliqueEntry(entry[2*i],true);
3913                  setOneFixesInCliqueEntry(entry[2*i+1],true);
3914                  type[i]=0;
3915                }
3916                fake = info.analyze(*solver_, 1,numberCliques,start,
3917                                    entry,type);
3918                delete [] type;
3919                delete [] entry;
3920#else
3921                fake = info.analyze(*newSolver, 1,-1);
3922#endif
3923                delete newSolver;
3924              } else {
3925                fake = info.analyze(*solver_, 1);
3926              }
3927              if (fake) {
3928                //fake->writeMps("fake");
3929                CglFakeClique cliqueGen(fake);
3930                cliqueGen.setStarCliqueReport(false);
3931                cliqueGen.setRowCliqueReport(false);
3932                cliqueGen.setMinViolation(0.1);
3933                addCutGenerator(&cliqueGen, 1, "Fake cliques", true, false, false, -200);
3934                generator_[numberCutGenerators_-1]->setTiming(true);
3935                for (int i = 0; i < numberCutGenerators_; i++) {
3936                  CglKnapsackCover * cutGen =
3937                  dynamic_cast<CglKnapsackCover *>(generator_[i]->generator());
3938                  if (cutGen) {
3939                    cutGen->createCliques(*fake,2,200,false);
3940                  }
3941                }
3942              }
3943            }
3944            if (probingInfo_->packDown()) {
3945#ifdef CLP_INVESTIGATE
3946                printf("%d implications on %d 0-1\n", toZero[number01], number01);
3947#endif
3948                // Create a cut generator that remembers implications discovered at root.
3949                CglImplication implication(probingInfo_);
3950                addCutGenerator(&implication, 1, "ImplicationCuts", true, false, false, -200);
3951                generator_[numberCutGenerators_-1]->setGlobalCuts(true);
3952            } else {
3953                delete probingInfo_;
3954                probingInfo_ = NULL;
3955            }
3956        } else {
3957            delete probingInfo_;
3958
3959            probingInfo_ = NULL;
3960        }
3961    }
3962    /*
3963      At this point, the root subproblem is infeasible or fathomed by bound
3964      (newNode == NULL), or we're live with an objective value that satisfies the
3965      current objective cutoff.
3966    */
3967    assert (!newNode || newNode->objectiveValue() <= cutoff) ;
3968    // Save address of root node as we don't want to delete it
3969    /*
3970      The common case is that the lp relaxation is feasible but doesn't satisfy
3971      integrality (i.e., newNode->branchingObject(), indicating we've been able to
3972      select a branching variable). Remove any cuts that have gone slack due to
3973      forcing monotone variables. Then tack on an CbcFullNodeInfo object and full
3974      basis (via createInfo()) and stash the new cuts in the nodeInfo (via
3975      addCuts()). If, by some miracle, we have an integral solution at the root
3976      (newNode->branchingObject() is NULL), takeOffCuts() will ensure that the solver holds
3977      a valid solution for use by setBestSolution().
3978    */
3979    CoinWarmStartBasis *lastws = NULL ;
3980    if (feasible && newNode->branchingObject()) {
3981        if (resolved) {
3982            takeOffCuts(cuts, false, NULL) ;
3983#     ifdef CHECK_CUT_COUNTS
3984            { printf("Number of rows after chooseBranch fix (root)"
3985                         "(active only) %d\n",
3986                         numberRowsAtContinuous_ + numberNewCuts_ + numberOldActiveCuts_) ;
3987                const CoinWarmStartBasis* debugws =
3988                    dynamic_cast <const CoinWarmStartBasis*>(solver_->getWarmStart()) ;
3989                debugws->print() ;
3990                delete debugws ;
3991            }
3992#     endif
3993        }
3994        //newNode->createInfo(this,NULL,NULL,NULL,NULL,0,0) ;
3995        //newNode->nodeInfo()->addCuts(cuts,
3996        //                       newNode->numberBranches(),whichGenerator_) ;
3997        if (lastws) delete lastws ;
3998        lastws = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
3999    }
4000    /*
4001      Continuous data to be used later
4002    */
4003    continuousObjective_ = solver_->getObjValue() * solver_->getObjSense();
4004    continuousInfeasibilities_ = 0 ;
4005    if (newNode) {
4006        continuousObjective_ = newNode->objectiveValue() ;
4007        delete [] continuousSolution_;
4008        continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
4009                                              numberColumns);
4010        continuousInfeasibilities_ = newNode->numberUnsatisfied() ;
4011    }
4012    /*
4013      Bound may have changed so reset in objects
4014    */
4015    {
4016        int i ;
4017        for (i = 0; i < numberObjects_; i++)
4018            object_[i]->resetBounds(solver_) ;
4019    }
4020    /*
4021      Feasible? Then we should have either a live node prepped for future
4022      expansion (indicated by variable() >= 0), or (miracle of miracles) an
4023      integral solution at the root node.
4024
4025      initializeInfo sets the reference counts in the nodeInfo object.  Since
4026      this node is still live, push it onto the heap that holds the live set.
4027    */
4028    if (newNode) {
4029        if (newNode->branchingObject()) {
4030            newNode->initializeInfo() ;
4031            tree_->push(newNode) ;
4032            // save pointer to root node - so can pick up bounds
4033            if (!topOfTree_)
4034              topOfTree_ = dynamic_cast<CbcFullNodeInfo *>(newNode->nodeInfo()) ;
4035            if (statistics_) {
4036                if (numberNodes2_ == maximumStatistics_) {
4037                    maximumStatistics_ = 2 * maximumStatistics_;
4038                    CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_];
4039                    memset(temp, 0, maximumStatistics_*sizeof(CbcStatistics *));
4040                    memcpy(temp, statistics_, numberNodes2_*sizeof(CbcStatistics *));
4041                    delete [] statistics_;
4042                    statistics_ = temp;
4043                }
4044                assert (!statistics_[numberNodes2_]);
4045                statistics_[numberNodes2_] = new CbcStatistics(newNode, this);
4046            }
4047            numberNodes2_++;
4048#     ifdef CHECK_NODE
4049            printf("Node %x on tree\n", newNode) ;
4050#     endif
4051        } else {
4052            // continuous is integer
4053            double objectiveValue = newNode->objectiveValue();
4054            setBestSolution(CBC_SOLUTION, objectiveValue,
4055                            solver_->getColSolution()) ;
4056            if (eventHandler) {
4057              // we are stopping anyway so no need to test return code
4058              eventHandler->event(CbcEventHandler::solution);
4059            }
4060            delete newNode ;
4061            newNode = NULL ;
4062        }
4063    }
4064
4065    if (printFrequency_ <= 0) {
4066        printFrequency_ = 1000 ;
4067        if (getNumCols() > 2000)
4068            printFrequency_ = 100 ;
4069    }
4070    /*
4071      It is possible that strong branching fixes one variable and then the code goes round
4072      again and again.  This can take too long.  So we need to warn user - just once.
4073    */
4074    numberLongStrong_ = 0;
4075    CbcNode * createdNode = NULL;
4076#ifdef CBC_THREAD
4077    if ((specialOptions_&2048) != 0)
4078        numberThreads_ = 0;
4079    if (numberThreads_ ) {
4080        nodeCompare_->sayThreaded(); // need to use addresses
4081        master_ = new CbcBaseModel(*this,
4082                                   (parallelMode() < -1) ? 1 : 0);
4083        masterThread_ = master_->masterThread();
4084    }
4085#endif
4086#ifdef COIN_HAS_CLP
4087    {
4088        OsiClpSolverInterface * clpSolver
4089        = dynamic_cast<OsiClpSolverInterface *> (solver_);
4090        if (clpSolver && !parentModel_) {
4091            clpSolver->computeLargestAway();
4092        }
4093    }
4094#endif
4095    /*
4096      At last, the actual branch-and-cut search loop, which will iterate until
4097      the live set is empty or we hit some limit (integrality gap, time, node
4098      count, etc.). The overall flow is to rebuild a subproblem, reoptimise using
4099      solveWithCuts(), choose a branching pattern with chooseBranch(), and finally
4100      add the node to the live set.
4101
4102      The first action is to winnow the live set to remove nodes which are worse
4103      than the current objective cutoff.
4104    */
4105    if (solver_->getRowCutDebuggerAlways()) {
4106        OsiRowCutDebugger * debuggerX = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
4107        const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger() ;
4108        if (!debugger) {
4109            // infeasible!!
4110            printf("before search\n");
4111            const double * lower = solver_->getColLower();
4112            const double * upper = solver_->getColUpper();
4113            const double * solution = debuggerX->optimalSolution();
4114            int numberColumns = solver_->getNumCols();
4115            for (int i = 0; i < numberColumns; i++) {
4116                if (solver_->isInteger(i)) {
4117                    if (solution[i] < lower[i] - 1.0e-6 || solution[i] > upper[i] + 1.0e-6)
4118                        printf("**** ");
4119                    printf("%d %g <= %g <= %g\n",
4120                           i, lower[i], solution[i], upper[i]);
4121                }
4122            }
4123            //abort();
4124        }
4125    }
4126    {
4127        // may be able to change cutoff now
4128        double cutoff = getCutoff();
4129        double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
4130        if (cutoff > bestObjective_ - increment) {
4131            cutoff = bestObjective_ - increment ;
4132            setCutoff(cutoff) ;
4133        }
4134    }
4135#ifdef CBC_THREAD
4136    bool goneParallel = false;
4137#endif
4138#define MAX_DEL_NODE 1
4139    CbcNode * delNode[MAX_DEL_NODE+1];
4140    int nDeleteNode = 0;
4141    // For Printing etc when parallel
4142    int lastEvery1000 = 0;
4143    int lastPrintEvery = 0;
4144    int numberConsecutiveInfeasible = 0;
4145#define PERTURB_IN_FATHOM
4146#ifdef PERTURB_IN_FATHOM
4147    // allow in fathom
4148    if ((moreSpecialOptions_& 262144) != 0)
4149      specialOptions_ |= 131072;
4150#endif
4151    while (true) {
4152        lockThread();
4153#ifdef COIN_HAS_CLP
4154        // See if we want dantzig row choice
4155        goToDantzig(100, savePivotMethod);
4156#endif
4157        if (tree_->empty()) {
4158#ifdef CBC_THREAD
4159            if (parallelMode() > 0 && master_) {
4160                int anyLeft = master_->waitForThreadsInTree(0);
4161                if (!anyLeft) {
4162                    master_->stopThreads(-1);
4163                    break;
4164                }
4165            } else {
4166                break;
4167            }
4168#else
4169            break;
4170#endif
4171        } else {
4172            unlockThread();
4173        }
4174        // If done 100 nodes see if worth trying reduction
4175        if (numberNodes_ == 50 || numberNodes_ == 100) {
4176#ifdef COIN_HAS_CLP
4177            OsiClpSolverInterface * clpSolver
4178            = dynamic_cast<OsiClpSolverInterface *> (solver_);
4179            if (clpSolver && ((specialOptions_&131072) == 0) && true) {
4180                ClpSimplex * simplex = clpSolver->getModelPtr();
4181                int perturbation = simplex->perturbation();
4182#ifdef CLP_INVESTIGATE
4183                printf("Testing its n,s %d %d solves n,s %d %d - pert %d\n",
4184                       numberIterations_, saveNumberIterations,
4185                       numberSolves_, saveNumberSolves, perturbation);
4186#endif
4187                if (perturbation == 50 && (numberIterations_ - saveNumberIterations) <
4188                        8*(numberSolves_ - saveNumberSolves)) {
4189                    // switch off perturbation
4190                    simplex->setPerturbation(100);
4191#ifdef CLP_INVESTIGATE
4192                    printf("Perturbation switched off\n");
4193#endif
4194                }
4195            }
4196#endif
4197            /*
4198              Decide if we want to do a restart.
4199            */
4200            if (saveSolver && (specialOptions_&(512 + 32768)) != 0) {
4201                bool tryNewSearch = solverCharacteristics_->reducedCostsAccurate() &&
4202                                    (getCutoff() < 1.0e20 && getCutoff() < checkCutoffForRestart);
4203                int numberColumns = getNumCols();
4204                if (tryNewSearch) {
4205                    checkCutoffForRestart = getCutoff() ;
4206#ifdef CLP_INVESTIGATE
4207                    printf("after %d nodes, cutoff %g - looking\n",
4208                           numberNodes_, getCutoff());
4209#endif
4210                    saveSolver->resolve();
4211                    double direction = saveSolver->getObjSense() ;
4212                    double gap = checkCutoffForRestart - saveSolver->getObjValue() * direction ;
4213                    double tolerance;
4214                    saveSolver->getDblParam(OsiDualTolerance, tolerance) ;
4215                    if (gap <= 0.0)
4216                        gap = tolerance;
4217                    gap += 100.0 * tolerance;
4218                    double integerTolerance = getDblParam(CbcIntegerTolerance) ;
4219
4220                    const double *lower = saveSolver->getColLower() ;
4221                    const double *upper = saveSolver->getColUpper() ;
4222                    const double *solution = saveSolver->getColSolution() ;
4223                    const double *reducedCost = saveSolver->getReducedCost() ;
4224
4225                    int numberFixed = 0 ;
4226                    int numberFixed2 = 0;
4227#ifdef COIN_DEVELOP
4228                    printf("gap %g\n", gap);
4229#endif
4230                    for (int i = 0 ; i < numberIntegers_ ; i++) {
4231                        int iColumn = integerVariable_[i] ;
4232                        double djValue = direction * reducedCost[iColumn] ;
4233                        if (upper[iColumn] - lower[iColumn] > integerTolerance) {
4234                            if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) {
4235                                //printf("%d to lb on dj of %g - bounds %g %g\n",
4236                                //     iColumn,djValue,lower[iColumn],upper[iColumn]);
4237                                saveSolver->setColUpper(iColumn, lower[iColumn]) ;
4238                                numberFixed++ ;
4239                            } else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) {
4240                                //printf("%d to ub on dj of %g - bounds %g %g\n",
4241                                //     iColumn,djValue,lower[iColumn],upper[iColumn]);
4242                                saveSolver->setColLower(iColumn, upper[iColumn]) ;
4243                                numberFixed++ ;
4244                            }
4245                        } else {
4246                            //printf("%d has dj of %g - already fixed to %g\n",
4247                            //     iColumn,djValue,lower[iColumn]);
4248                            numberFixed2++;
4249                        }
4250                    }
4251#ifdef COIN_DEVELOP
4252                    if ((specialOptions_&1) != 0) {
4253                        const OsiRowCutDebugger *debugger = saveSolver->getRowCutDebugger() ;
4254                        if (debugger) {
4255                            printf("Contains optimal\n") ;
4256                            OsiSolverInterface * temp = saveSolver->clone();
4257                            const double * solution = debugger->optimalSolution();
4258                            const double *lower = temp->getColLower() ;
4259                            const double *upper = temp->getColUpper() ;
4260                            int n = temp->getNumCols();
4261                            for (int i = 0; i < n; i++) {
4262                                if (temp->isInteger(i)) {
4263                                    double value = floor(solution[i] + 0.5);
4264                                    assert (value >= lower[i] && value <= upper[i]);
4265                                    temp->setColLower(i, value);
4266                                    temp->setColUpper(i, value);
4267                                }
4268                            }
4269                            temp->writeMps("reduced_fix");
4270                            delete temp;
4271                            saveSolver->writeMps("reduced");
4272                        } else {
4273                            abort();
4274                        }
4275                    }
4276                    printf("Restart could fix %d integers (%d already fixed)\n",
4277                           numberFixed + numberFixed2, numberFixed2);
4278#endif
4279                    numberFixed += numberFixed2;
4280                    if (numberFixed*10 < numberColumns && numberFixed*4 < numberIntegers_)
4281                        tryNewSearch = false;
4282                }
4283                if (tryNewSearch) {
4284                    // back to solver without cuts?
4285                    OsiSolverInterface * solver2 = saveSolver->clone();
4286                    const double *lower = saveSolver->getColLower() ;
4287                    const double *upper = saveSolver->getColUpper() ;
4288                    for (int i = 0 ; i < numberIntegers_ ; i++) {
4289                        int iColumn = integerVariable_[i] ;
4290                        solver2->setColLower(iColumn, lower[iColumn]);
4291                        solver2->setColUpper(iColumn, upper[iColumn]);
4292                    }
4293                    // swap
4294                    delete saveSolver;
4295                    saveSolver = solver2;
4296                    double * newSolution = new double[numberColumns];
4297                    double objectiveValue = checkCutoffForRestart;
4298                    // Save the best solution so far.
4299                    CbcSerendipity heuristic(*this);
4300                    if (bestSolution_)
4301                        heuristic.setInputSolution(bestSolution_, bestObjective_);
4302                    // Magic number
4303                    heuristic.setFractionSmall(0.8);
4304                    // `pumpTune' to stand-alone solver for explanations.
4305                    heuristic.setFeasibilityPumpOptions(1008013);
4306                    // Use numberNodes to say how many are original rows
4307                    heuristic.setNumberNodes(continuousSolver_->getNumRows());
4308#ifdef COIN_DEVELOP
4309                    if (continuousSolver_->getNumRows() <
4310                            solver_->getNumRows())
4311                        printf("%d rows added ZZZZZ\n",
4312                               solver_->getNumRows() - continuousSolver_->getNumRows());
4313#endif
4314                    int returnCode = heuristic.smallBranchAndBound(saveSolver,
4315                                     -1, newSolution,
4316                                     objectiveValue,
4317                                     checkCutoffForRestart, "Reduce");
4318                    if (returnCode < 0) {
4319#ifdef COIN_DEVELOP
4320                        printf("Restart - not small enough to do search after fixing\n");
4321#endif
4322                        delete [] newSolution;
4323                    } else {
4324                        // 1 for sol'n, 2 for finished, 3 for both
4325                        if ((returnCode&1) != 0) {
4326                            // increment number of solutions so other heuristics can test
4327                            numberSolutions_++;
4328                            numberHeuristicSolutions_++;
4329                            lastHeuristic_ = NULL;
4330                            setBestSolution(CBC_ROUNDING, objectiveValue, newSolution) ;
4331                        }
4332                        delete [] newSolution;
4333#ifdef CBC_THREAD
4334                        if (master_) {
4335                            lockThread();
4336                            if (parallelMode() > 0) {
4337                                while (master_->waitForThreadsInTree(0)) {
4338                                    lockThread();
4339                                    double dummyBest;
4340                                    tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4341                                    //unlockThread();
4342                                }
4343                            } else {
4344                                double dummyBest;
4345                                tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4346                            }
4347                            master_->waitForThreadsInTree(2);
4348                            delete master_;
4349                            master_ = NULL;
4350                            masterThread_ = NULL;
4351                        }
4352#endif
4353                        if (tree_->size()) {
4354                            double dummyBest;
4355                            tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4356                        }
4357                        break;
4358                    }
4359                }
4360                delete saveSolver;
4361                saveSolver = NULL;
4362            }
4363        }
4364        /*
4365          Check for abort on limits: node count, solution count, time, integrality gap.
4366        */
4367        if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
4368                numberSolutions_ < intParam_[CbcMaxNumSol] &&
4369                !maximumSecondsReached() &&
4370                !stoppedOnGap_ && !eventHappened_ && (maximumNumberIterations_ < 0 ||
4371                                                      numberIterations_ < maximumNumberIterations_))) {
4372            // out of loop
4373            break;
4374        }
4375#ifdef BONMIN
4376        assert(!solverCharacteristics_->solutionAddsCuts() || solverCharacteristics_->mipFeasible());
4377#endif
4378// Sets percentage of time when we try diving. Diving requires a bit of heap reorganisation, because
4379// we need to replace the comparison function to dive, and that requires reordering to retain the
4380// heap property.
4381#define DIVE_WHEN 1000
4382#define DIVE_STOP 2000
4383        int kNode = numberNodes_ % 4000;
4384        if (numberNodes_<100000 && kNode>DIVE_WHEN && kNode <= DIVE_STOP) {
4385            if (!parallelMode()) {
4386                if (kNode == DIVE_WHEN + 1 || numberConsecutiveInfeasible > 1) {
4387                    CbcCompareDefault * compare = dynamic_cast<CbcCompareDefault *>
4388                                                  (nodeCompare_);
4389                    // Don't interfere if user has replaced the compare function.
4390                    if (compare) {
4391                        //printf("Redoing tree\n");
4392                        compare->startDive(this);
4393                        numberConsecutiveInfeasible = 0;
4394                    }
4395                }
4396            }
4397        }
4398        // replace current cutoff?
4399        if (cutoff > getCutoff()) {
4400            double newCutoff = getCutoff();
4401            if (analyzeResults_) {
4402                // see if we could fix any (more)
4403                int n = 0;
4404                double * newLower = analyzeResults_;
4405                double * objLower = newLower + numberIntegers_;
4406                double * newUpper = objLower + numberIntegers_;
4407                double * objUpper = newUpper + numberIntegers_;
4408                for (int i = 0; i < numberIntegers_; i++) {
4409                    if (objLower[i] > newCutoff) {
4410                        n++;
4411                        if (objUpper[i] > newCutoff) {
4412                            newCutoff = -COIN_DBL_MAX;
4413                            break;
4414                        }
4415                    } else if (objUpper[i] > newCutoff) {
4416                        n++;
4417                    }
4418                }
4419                if (newCutoff == -COIN_DBL_MAX) {
4420                  COIN_DETAIL_PRINT(printf("Root analysis says finished\n"));
4421                } else if (n > numberFixedNow_) {
4422                  COIN_DETAIL_PRINT(printf("%d more fixed by analysis - now %d\n", n - numberFixedNow_, n));
4423                    numberFixedNow_ = n;
4424                }
4425            }
4426            if (eventHandler) {
4427                if (!eventHandler->event(CbcEventHandler::solution)) {
4428                    eventHappened_ = true; // exit
4429                }
4430                newCutoff = getCutoff();
4431            }
4432            lockThread();
4433            /*
4434              Clean the tree to reflect the new solution, then see if the
4435              node comparison predicate wants to make any changes. If so,
4436              call setComparison for the side effect of rebuilding the heap.
4437            */
4438            tree_->cleanTree(this,newCutoff,bestPossibleObjective_) ;
4439            if (nodeCompare_->newSolution(this) ||
4440                nodeCompare_->newSolution(this,continuousObjective_,
4441                                          continuousInfeasibilities_)) {
4442              tree_->setComparison(*nodeCompare_) ;
4443            }
4444            if (tree_->empty()) {
4445                continue;
4446            }
4447            unlockThread();
4448        }
4449        cutoff = getCutoff() ;
4450        /*
4451            Periodic activities: Opportunities to
4452            + tweak the nodeCompare criteria,
4453            + check if we've closed the integrality gap enough to quit,
4454            + print a summary line to let the user know we're working
4455        */
4456        if (numberNodes_ >= lastEvery1000) {
4457            lockThread();
4458#ifdef COIN_HAS_CLP
4459            // See if we want dantzig row choice
4460            goToDantzig(1000, savePivotMethod);
4461#endif
4462            lastEvery1000 = numberNodes_ + 1000;
4463            bool redoTree = nodeCompare_->every1000Nodes(this, numberNodes_) ;
4464#ifdef CHECK_CUT_SIZE
4465            verifyCutSize (tree_, *this);
4466#endif
4467            // redo tree if requested
4468            if (redoTree)
4469                tree_->setComparison(*nodeCompare_) ;
4470            unlockThread();
4471        }
4472        // Had hotstart before, now switched off
4473        if (saveCompare && !hotstartSolution_) {
4474            // hotstart switched off
4475            delete nodeCompare_; // off depth first
4476            nodeCompare_ = saveCompare;
4477            saveCompare = NULL;
4478            // redo tree
4479            lockThread();
4480            tree_->setComparison(*nodeCompare_) ;
4481            unlockThread();
4482        }
4483        if (numberNodes_ >= lastPrintEvery) {
4484            lastPrintEvery = numberNodes_ + printFrequency_;
4485            lockThread();
4486            int nNodes = tree_->size() ;
4487
4488            //MODIF PIERRE
4489            bestPossibleObjective_ = tree_->getBestPossibleObjective();
4490#ifdef CBC_THREAD
4491            if (parallelMode() > 0 && master_) {
4492              // need to adjust for ones not on tree
4493              int numberThreads = master_->numberThreads();
4494              for (int i=0;i<numberThreads;i++) {
4495                CbcThread * child = master_->child(i);
4496                if (child->node()) {
4497                  // adjust
4498                  double value = child->node()->objectiveValue();
4499                  bestPossibleObjective_ = CoinMin(bestPossibleObjective_, value);
4500                }
4501              }
4502            }
4503#endif
4504            unlockThread();
4505#ifdef CLP_INVESTIGATE
4506            if (getCutoff() < 1.0e20) {
4507                if (fabs(getCutoff() - (bestObjective_ - getCutoffIncrement())) > 1.0e-6 &&
4508                        !parentModel_)
4509                    printf("model cutoff in status %g, best %g, increment %g\n",
4510                           getCutoff(), bestObjective_, getCutoffIncrement());
4511                assert (getCutoff() < bestObjective_ - getCutoffIncrement() +
4512                        1.0e-6 + 1.0e-10*fabs(bestObjective_));
4513            }
4514#endif
4515            if (!intParam_[CbcPrinting]) {
4516                messageHandler()->message(CBC_STATUS, messages())
4517                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4518                << getCurrentSeconds()
4519                << CoinMessageEol ;
4520            } else if (intParam_[CbcPrinting] == 1) {
4521                messageHandler()->message(CBC_STATUS2, messages())
4522                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4523                << tree_->lastDepth() << tree_->lastUnsatisfied() 
4524                << tree_->lastObjective() << numberIterations_
4525                << getCurrentSeconds()
4526                << CoinMessageEol ;
4527            } else if (!numberExtraIterations_) {
4528                messageHandler()->message(CBC_STATUS2, messages())
4529                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4530                << tree_->lastDepth() << tree_->lastUnsatisfied() << numberIterations_
4531                << getCurrentSeconds()
4532                << CoinMessageEol ;
4533            } else {
4534                messageHandler()->message(CBC_STATUS3, messages())
4535                << numberNodes_ << numberExtraNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4536                << tree_->lastDepth() << tree_->lastUnsatisfied() << numberIterations_ << numberExtraIterations_
4537                << getCurrentSeconds()
4538                << CoinMessageEol ;
4539            }
4540#ifdef COIN_HAS_NTY
4541            if (symmetryInfo_) 
4542              symmetryInfo_->statsOrbits(this,1);
4543#endif
4544            if (eventHandler && !eventHandler->event(CbcEventHandler::treeStatus)) {
4545                eventHappened_ = true; // exit
4546            }
4547        }
4548        // See if can stop on gap
4549        if(canStopOnGap()) {
4550            stoppedOnGap_ = true ;
4551        }
4552
4553#ifdef CHECK_NODE_FULL
4554        verifyTreeNodes(tree_, *this) ;
4555#   endif
4556#   ifdef CHECK_CUT_COUNTS
4557        verifyCutCounts(tree_, *this) ;
4558#   endif
4559        /*
4560          Now we come to the meat of the loop. To create the active subproblem, we'll
4561          pop the most promising node in the live set, rebuild the subproblem it
4562          represents, and then execute the current arm of the branch to create the
4563          active subproblem.
4564        */
4565        CbcNode * node = NULL;
4566#ifdef CBC_THREAD
4567        if (!parallelMode() || parallelMode() == -1) {
4568#endif
4569            node = tree_->bestNode(cutoff) ;
4570            // Possible one on tree worse than cutoff
4571            // Weird comparison function can leave ineligible nodes on tree
4572            if (!node || node->objectiveValue() > cutoff)
4573                continue;
4574            // Do main work of solving node here
4575            doOneNode(this, node, createdNode);
4576#ifdef JJF_ZERO
4577            if (node) {
4578              if (createdNode) {
4579                printf("Node %d depth %d, created %d depth %d\n",
4580                       node->nodeNumber(), node->depth(),
4581                       createdNode->nodeNumber(), createdNode->depth());
4582              } else {
4583                printf("Node %d depth %d,  no created node\n",
4584                       node->nodeNumber(), node->depth());
4585              }
4586            } else if (createdNode) {
4587              printf("Node exhausted, created %d depth %d\n",
4588                     createdNode->nodeNumber(), createdNode->depth());
4589            } else {
4590              printf("Node exhausted,  no created node\n");
4591              numberConsecutiveInfeasible = 2;
4592            }
4593#endif
4594            //if (createdNode)
4595            //numberConsecutiveInfeasible=0;
4596            //else
4597            //numberConsecutiveInfeasible++;
4598#ifdef CBC_THREAD
4599        } else if (parallelMode() > 0) {
4600            //lockThread();
4601            //node = tree_->bestNode(cutoff) ;
4602            // Possible one on tree worse than cutoff
4603            if (true || !node || node->objectiveValue() > cutoff) {
4604                assert (master_);
4605                if (master_) {
4606                    int anyLeft = master_->waitForThreadsInTree(1);
4607                    // may need to go round again
4608                    if (anyLeft) {
4609                        continue;
4610                    } else {
4611                        master_->stopThreads(-1);
4612                    }
4613                }
4614            }
4615            //unlockThread();
4616        } else {
4617            // Deterministic parallel
4618          if ((tree_->size() < CoinMax(numberThreads_, 8)||
4619               hotstartSolution_) && !goneParallel) {
4620                node = tree_->bestNode(cutoff) ;
4621                // Possible one on tree worse than cutoff
4622                if (!node || node->objectiveValue() > cutoff)
4623                    continue;
4624                // Do main work of solving node here
4625                doOneNode(this, node, createdNode);
4626                assert (createdNode);
4627                if (!createdNode->active()) {
4628                    delete createdNode;
4629                    createdNode = NULL;
4630                } else {
4631                    // Say one more pointing to this
4632                    node->nodeInfo()->increment() ;
4633                    tree_->push(createdNode) ;
4634                }
4635                if (node->active()) {
4636                    assert (node->nodeInfo());
4637                    if (node->nodeInfo()->numberBranchesLeft()) {
4638                        tree_->push(node) ;
4639                    } else {
4640                        node->setActive(false);
4641                    }
4642                } else {
4643                    if (node->nodeInfo()) {
4644                        if (!node->nodeInfo()->numberBranchesLeft())
4645                            node->nodeInfo()->allBranchesGone(); // can clean up
4646                        // So will delete underlying stuff
4647                        node->setActive(true);
4648                    }
4649                    delNode[nDeleteNode++] = node;
4650                    node = NULL;
4651                }
4652                if (nDeleteNode >= MAX_DEL_NODE) {
4653                    for (int i = 0; i < nDeleteNode; i++) {
4654                        //printf("trying to del %d %x\n",i,delNode[i]);
4655                        delete delNode[i];
4656                        //printf("done to del %d %x\n",i,delNode[i]);
4657                    }
4658                    nDeleteNode = 0;
4659                }
4660            } else {
4661                // Split and solve
4662                master_->deterministicParallel();
4663                goneParallel = true;
4664            }
4665        }
4666#endif
4667    }
4668    if (nDeleteNode) {
4669        for (int i = 0; i < nDeleteNode; i++) {
4670            delete delNode[i];
4671        }
4672        nDeleteNode = 0;
4673    }
4674#ifdef CBC_THREAD
4675    if (master_) {
4676        master_->stopThreads(-1);
4677        master_->waitForThreadsInTree(2);
4678        delete master_;
4679        master_ = NULL;
4680        masterThread_ = NULL;
4681        // adjust time to allow for children on some systems
4682        //dblParam_[CbcStartSeconds] -= CoinCpuTimeJustChildren();
4683    }
4684#endif
4685    /*
4686      End of the non-abort actions. The next block of code is executed if we've
4687      aborted because we hit one of the limits. Clean up by deleting the live set
4688      and break out of the node processing loop. Note that on an abort, node may
4689      have been pushed back onto the tree for further processing, in which case
4690      it'll be deleted in cleanTree. We need to check.
4691    */
4692    if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
4693          numberSolutions_ < intParam_[CbcMaxNumSol] &&
4694          !maximumSecondsReached() &&
4695          !stoppedOnGap_ &&
4696          !eventHappened_ &&
4697          (maximumNumberIterations_ < 0 || numberIterations_ < maximumNumberIterations_))
4698         ) {
4699        if (tree_->size()) {
4700            double dummyBest;
4701            tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4702        }
4703        delete nextRowCut_;
4704        /* order is important here:
4705         * maximumSecondsReached() should be checked before eventHappened_ and
4706         * isNodeLimitReached() should be checked after eventHappened_
4707         * reason is, that at timelimit, eventHappened_ is set to true to make Cbc stop fast
4708         *   and if Ctrl+C is hit, then the nodelimit is set to -1 to make Cbc stop
4709         */
4710        if (stoppedOnGap_) {
4711            messageHandler()->message(CBC_GAP, messages())
4712            << bestObjective_ - bestPossibleObjective_
4713            << dblParam_[CbcAllowableGap]
4714            << dblParam_[CbcAllowableFractionGap]*100.0
4715            << CoinMessageEol ;
4716            secondaryStatus_ = 2;
4717            status_ = 0 ;
4718        } else if (maximumSecondsReached()) {
4719            handler_->message(CBC_MAXTIME, messages_) << CoinMessageEol ;
4720            secondaryStatus_ = 4;
4721            status_ = 1 ;
4722        } else if (numberSolutions_ >= intParam_[CbcMaxNumSol]) {
4723            handler_->message(CBC_MAXSOLS, messages_) << CoinMessageEol ;
4724            secondaryStatus_ = 6;
4725            status_ = 1 ;
4726        } else if (isNodeLimitReached()) {
4727            handler_->message(CBC_MAXNODES, messages_) << CoinMessageEol ;
4728            secondaryStatus_ = 3;
4729            status_ = 1 ;
4730        } else if (maximumNumberIterations_ >= 0 && numberIterations_ >= maximumNumberIterations_) {
4731            handler_->message(CBC_MAXITERS, messages_) << CoinMessageEol ;
4732            secondaryStatus_ = 8;
4733            status_ = 1 ;
4734        } else {
4735            handler_->message(CBC_EVENT, messages_) << CoinMessageEol ;
4736            secondaryStatus_ = 5;
4737            status_ = 5 ;
4738        }
4739    }
4740    /*
4741      That's it, we've exhausted the search tree, or broken out of the loop because
4742      we hit some limit on evaluation.
4743
4744      We may have got an intelligent tree so give it one more chance
4745    */
4746    // Tell solver we are not in Branch and Cut
4747    solver_->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo, NULL) ;
4748    tree_->endSearch();
4749    //  If we did any sub trees - did we give up on any?
4750    if ( numberStoppedSubTrees_)
4751        status_ = 1;
4752    numberNodes_ += numberExtraNodes_;
4753    numberIterations_ += numberExtraIterations_;
4754    if (eventHandler) {
4755        eventHandler->event(CbcEventHandler::endSearch);
4756    }
4757    if (!status_) {
4758        // Set best possible unless stopped on gap
4759        if (secondaryStatus_ != 2)
4760            bestPossibleObjective_ = bestObjective_;
4761        handler_->message(CBC_END_GOOD, messages_)
4762        << bestObjective_ << numberIterations_ << numberNodes_ << getCurrentSeconds()
4763        << CoinMessageEol ;
4764    } else {
4765        handler_->message(CBC_END, messages_)
4766        << bestObjective_ << bestPossibleObjective_
4767        << numberIterations_ << numberNodes_ << getCurrentSeconds()
4768        << CoinMessageEol ;
4769    }
4770    if (numberStrongIterations_)
4771        handler_->message(CBC_STRONG_STATS, messages_)
4772        << strongInfo_[0] << numberStrongIterations_ << strongInfo_[2]
4773        << strongInfo_[1] << CoinMessageEol ;
4774    if (!numberExtraNodes_)
4775        handler_->message(CBC_OTHER_STATS, messages_)
4776        << maximumDepthActual_
4777        << numberDJFixed_ << CoinMessageEol ;
4778    else
4779        handler_->message(CBC_OTHER_STATS2, messages_)
4780        << maximumDepthActual_
4781        << numberDJFixed_ << numberExtraNodes_ << numberExtraIterations_
4782        << CoinMessageEol ;
4783#ifdef COIN_HAS_NTY
4784    if (symmetryInfo_) 
4785      symmetryInfo_->statsOrbits(this,1);
4786#endif
4787    if (doStatistics == 100) {
4788        for (int i = 0; i < numberObjects_; i++) {
4789            CbcSimpleIntegerDynamicPseudoCost * obj =
4790                dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object_[i]) ;
4791            if (obj)
4792                obj->print();
4793        }
4794    }
4795    if (statistics_) {
4796        // report in some way
4797        int * lookup = new int[numberObjects_];
4798        int i;
4799        for (i = 0; i < numberObjects_; i++)
4800            lookup[i] = -1;
4801        bool goodIds = false; //true;
4802        for (i = 0; i < numberObjects_; i++) {
4803            int iColumn = object_[i]->columnNumber();
4804            if (iColumn >= 0 && iColumn < numberColumns) {
4805                if (lookup[i] == -1) {
4806                    lookup[i] = iColumn;
4807                } else {
4808                    goodIds = false;
4809                    break;
4810                }
4811            } else {
4812                goodIds = false;
4813                break;
4814            }
4815        }
4816        if (!goodIds) {
4817            delete [] lookup;
4818            lookup = NULL;
4819        }
4820        if (doStatistics >= 3) {
4821            printf("  node parent depth column   value                    obj      inf\n");
4822            for ( i = 0; i < numberNodes2_; i++) {
4823                statistics_[i]->print(lookup);
4824            }
4825        }
4826        if (doStatistics > 1) {
4827            // Find last solution
4828            int k;
4829            for (k = numberNodes2_ - 1; k >= 0; k--) {
4830                if (statistics_[k]->endingObjective() != COIN_DBL_MAX &&
4831                        !statistics_[k]->endingInfeasibility())
4832                    break;
4833            }
4834            if (k >= 0) {
4835                int depth = statistics_[k]->depth();
4836                int * which = new int[depth+1];
4837                for (i = depth; i >= 0; i--) {
4838                    which[i] = k;
4839                    k = statistics_[k]->parentNode();
4840                }
4841                printf("  node parent depth column   value                    obj      inf\n");
4842                for (i = 0; i <= depth; i++) {
4843                    statistics_[which[i]]->print(lookup);
4844                }
4845                delete [] which;
4846            }
4847        }
4848        // now summary
4849        int maxDepth = 0;
4850        double averageSolutionDepth = 0.0;
4851        int numberSolutions = 0;
4852        double averageCutoffDepth = 0.0;
4853        double averageSolvedDepth = 0.0;
4854        int numberCutoff = 0;
4855        int numberDown = 0;
4856        int numberFirstDown = 0;
4857        double averageInfDown = 0.0;
4858        double averageObjDown = 0.0;
4859        int numberCutoffDown = 0;
4860        int numberUp = 0;
4861        int numberFirstUp = 0;
4862        double averageInfUp = 0.0;
4863        double averageObjUp = 0.0;
4864        int numberCutoffUp = 0;
4865        double averageNumberIterations1 = 0.0;
4866        double averageValue = 0.0;
4867        for ( i = 0; i < numberNodes2_; i++) {
4868            int depth =  statistics_[i]->depth();
4869            int way =  statistics_[i]->way();
4870            double value = statistics_[i]->value();
4871            double startingObjective =  statistics_[i]->startingObjective();
4872            int startingInfeasibility = statistics_[i]->startingInfeasibility();
4873            double endingObjective = statistics_[i]->endingObjective();
4874            int endingInfeasibility = statistics_[i]->endingInfeasibility();
4875            maxDepth = CoinMax(depth, maxDepth);
4876            // Only for completed
4877            averageNumberIterations1 += statistics_[i]->numberIterations();
4878            averageValue += value;
4879            if (endingObjective != COIN_DBL_MAX && !endingInfeasibility) {
4880                numberSolutions++;
4881                averageSolutionDepth += depth;
4882            }
4883            if (endingObjective == COIN_DBL_MAX) {
4884                numberCutoff++;
4885                averageCutoffDepth += depth;
4886                if (way < 0) {
4887                    numberDown++;
4888                    numberCutoffDown++;
4889                    if (way == -1)
4890                        numberFirstDown++;
4891                } else {
4892                    numberUp++;
4893                    numberCutoffUp++;
4894                    if (way == 1)
4895                        numberFirstUp++;
4896                }
4897            } else {
4898                averageSolvedDepth += depth;
4899                if (way < 0) {
4900                    numberDown++;
4901                    averageInfDown += startingInfeasibility - endingInfeasibility;
4902                    averageObjDown += endingObjective - startingObjective;
4903                    if (way == -1)
4904                        numberFirstDown++;
4905                } else {
4906                    numberUp++;
4907                    averageInfUp += startingInfeasibility - endingInfeasibility;
4908                    averageObjUp += endingObjective - startingObjective;
4909                    if (way == 1)
4910                        numberFirstUp++;
4911                }
4912            }
4913        }
4914        // Now print
4915        if (numberSolutions)
4916            averageSolutionDepth /= static_cast<double> (numberSolutions);
4917        int numberSolved = numberNodes2_ - numberCutoff;
4918        double averageNumberIterations2 = numberIterations_ - averageNumberIterations1
4919                                          - numberIterationsAtContinuous;
4920        if (numberCutoff) {
4921            averageCutoffDepth /= static_cast<double> (numberCutoff);
4922            averageNumberIterations2 /= static_cast<double> (numberCutoff);
4923        }
4924        if (numberNodes2_)
4925            averageValue /= static_cast<double> (numberNodes2_);
4926        if (numberSolved) {
4927            averageNumberIterations1 /= static_cast<double> (numberSolved);
4928            averageSolvedDepth /= static_cast<double> (numberSolved);
4929        }
4930        printf("%d solution(s) were found (by branching) at an average depth of %g\n",
4931               numberSolutions, averageSolutionDepth);
4932        printf("average value of variable being branched on was %g\n",
4933               averageValue);
4934        printf("%d nodes were cutoff at an average depth of %g with iteration count of %g\n",
4935               numberCutoff, averageCutoffDepth, averageNumberIterations2);
4936        printf("%d nodes were solved at an average depth of %g with iteration count of %g\n",
4937               numberSolved, averageSolvedDepth, averageNumberIterations1);
4938        if (numberDown) {
4939            averageInfDown /= static_cast<double> (numberDown);
4940            averageObjDown /= static_cast<double> (numberDown);
4941        }
4942        printf("Down %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
4943               numberDown, numberFirstDown, numberDown - numberFirstDown, numberCutoffDown,
4944               averageInfDown, averageObjDown);
4945        if (numberUp) {
4946            averageInfUp /= static_cast<double> (numberUp);
4947            averageObjUp /= static_cast<double> (numberUp);
4948        }
4949        printf("Up %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
4950               numberUp, numberFirstUp, numberUp - numberFirstUp, numberCutoffUp,
4951               averageInfUp, averageObjUp);
4952        for ( i = 0; i < numberNodes2_; i++)
4953            delete statistics_[i];
4954        delete [] statistics_;
4955        statistics_ = NULL;
4956        maximumStatistics_ = 0;
4957        delete [] lookup;
4958    }
4959    /*
4960      If we think we have a solution, restore and confirm it with a call to
4961      setBestSolution().  We need to reset the cutoff value so as not to fathom
4962      the solution on bounds.  Note that calling setBestSolution( ..., true)
4963      leaves the continuousSolver_ bounds vectors fixed at the solution value.
4964
4965      Running resolve() here is a failsafe --- setBestSolution has already
4966      reoptimised using the continuousSolver_. If for some reason we fail to
4967      prove optimality, run the problem again after instructing the solver to
4968      tell us more.
4969
4970      If all looks good, replace solver_ with continuousSolver_, so that the
4971      outside world will be able to obtain information about the solution using
4972      public methods.
4973
4974      Don't replace if we are trying to save cuts
4975    */
4976    if (bestSolution_ && (solverCharacteristics_->solverType() < 2 || solverCharacteristics_->solverType() == 4) && 
4977        ((specialOptions_&8388608)==0||(specialOptions_&2048)!=0)) {
4978        setCutoff(1.0e50) ; // As best solution should be worse than cutoff
4979        // change cutoff as constraint if wanted
4980        if (cutoffRowNumber_>=0) {
4981          if (solver_->getNumRows()>cutoffRowNumber_)
4982            solver_->setRowUpper(cutoffRowNumber_,1.0e50);
4983        }
4984        // also in continuousSolver_
4985        if (continuousSolver_) {
4986          // Solvers know about direction
4987          double direction = solver_->getObjSense();
4988          continuousSolver_->setDblParam(OsiDualObjectiveLimit, 1.0e50*direction);
4989        }
4990        phase_ = 5;
4991        double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
4992        if ((specialOptions_&4) == 0)
4993            bestObjective_ += 100.0 * increment + 1.0e-3; // only set if we are going to solve
4994        setBestSolution(CBC_END_SOLUTION, bestObjective_, bestSolution_, 1) ;
4995        continuousSolver_->resolve() ;
4996        if (!continuousSolver_->isProvenOptimal()) {
4997            continuousSolver_->messageHandler()->setLogLevel(2) ;
4998            continuousSolver_->initialSolve() ;
4999        }
5000        delete solver_ ;
5001        // above deletes solverCharacteristics_
5002        solverCharacteristics_ = NULL;
5003        solver_ = continuousSolver_ ;
5004        setPointers(solver_);
5005        continuousSolver_ = NULL ;
5006    }
5007    /*
5008      Clean up dangling objects. continuousSolver_ may already be toast.
5009    */
5010    delete lastws ;
5011    if (saveObjects) {
5012        for (int i = 0; i < numberObjects_; i++)
5013            delete saveObjects[i];
5014        delete [] saveObjects;
5015    }
5016    numberStrong_ = saveNumberStrong;
5017    numberBeforeTrust_ = saveNumberBeforeTrust;
5018    delete [] whichGenerator_ ;
5019    whichGenerator_ = NULL;
5020    delete [] lowerBefore ;
5021    delete [] upperBefore ;
5022    delete [] walkback_ ;
5023    walkback_ = NULL ;
5024    delete [] lastNodeInfo_ ;
5025    lastNodeInfo_ = NULL;
5026    delete [] lastNumberCuts_ ;
5027    lastNumberCuts_ = NULL;
5028    delete [] lastCut_;
5029    lastCut_ = NULL;
5030    delete [] addedCuts_ ;
5031    addedCuts_ = NULL ;
5032    //delete persistentInfo;
5033    // Get rid of characteristics
5034    solverCharacteristics_ = NULL;
5035    if (continuousSolver_) {
5036        delete continuousSolver_ ;
5037        continuousSolver_ = NULL ;
5038    }
5039    /*
5040      Destroy global cuts by replacing with an empty OsiCuts object.
5041    */
5042    globalCuts_ = CbcRowCuts() ;
5043    delete globalConflictCuts_;
5044    globalConflictCuts_=NULL;
5045    if (!bestSolution_ && (specialOptions_&8388608)==0) {
5046        // make sure lp solver is infeasible
5047        int numberColumns = solver_->getNumCols();
5048        const double * columnLower = solver_->getColLower();
5049        int iColumn;
5050        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
5051            if (solver_->isInteger(iColumn))
5052                solver_->setColUpper(iColumn, columnLower[iColumn]);
5053        }
5054        solver_->initialSolve();
5055    }
5056#ifdef COIN_HAS_CLP
5057    {
5058        OsiClpSolverInterface * clpSolver
5059        = dynamic_cast<OsiClpSolverInterface *> (solver_);
5060        if (clpSolver) {
5061            // Possible restore of pivot method
5062            if (savePivotMethod) {
5063                // model may have changed
5064                savePivotMethod->setModel(NULL);
5065                clpSolver->getModelPtr()->setDualRowPivotAlgorithm(*savePivotMethod);
5066                delete savePivotMethod;
5067            }
5068            clpSolver->setLargestAway(-1.0);
5069        }
5070    }
5071#endif
5072    if ((fastNodeDepth_ >= 1000000 || (moreSpecialOptions_&33554432)!=0)
5073         && !parentModel_) {
5074      // delete object off end
5075      delete object_[numberObjects_];
5076      if ((moreSpecialOptions_&33554432)==0)
5077        fastNodeDepth_ -= 1000000;
5078    }
5079    delete saveSolver;
5080    // Undo preprocessing performed during BaB.
5081    if (strategy_ && strategy_->preProcessState() > 0) {
5082        // undo preprocessing
5083        CglPreProcess * process = strategy_->process();
5084        assert (process);
5085        int n = originalSolver->getNumCols();
5086        if (bestSolution_) {
5087            delete [] bestSolution_;
5088            bestSolution_ = new double [n];
5089            process->postProcess(*solver_);
5090        }
5091        strategy_->deletePreProcess();
5092        // Solution now back in originalSolver
5093        delete solver_;
5094        solver_ = originalSolver;
5095        if (bestSolution_) {
5096            bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
5097            memcpy(bestSolution_, solver_->getColSolution(), n*sizeof(double));
5098        }
5099        // put back original objects if there were any
5100        if (originalObject) {
5101            int iColumn;
5102            assert (ownObjects_);
5103            for (iColumn = 0; iColumn < numberObjects_; iColumn++)
5104                delete object_[iColumn];
5105            delete [] object_;
5106            numberObjects_ = numberOriginalObjects;
5107            object_ = originalObject;
5108            delete [] integerVariable_;
5109            numberIntegers_ = 0;
5110            for (iColumn = 0; iColumn < n; iColumn++) {
5111                if (solver_->isInteger(iColumn))
5112                    numberIntegers_++;
5113            }
5114            integerVariable_ = new int[numberIntegers_];
5115            numberIntegers_ = 0;
5116            for (iColumn = 0; iColumn < n; iColumn++) {
5117                if (solver_->isInteger(iColumn))
5118                    integerVariable_[numberIntegers_++] = iColumn;
5119            }
5120        }
5121    }
5122    if (flipObjective)
5123      flipModel();
5124#ifdef COIN_HAS_CLP
5125    {
5126        OsiClpSolverInterface * clpSolver
5127        = dynamic_cast<OsiClpSolverInterface *> (solver_);
5128        if (clpSolver)
5129            clpSolver->setFakeObjective(reinterpret_cast<double *> (NULL));
5130    }
5131#endif
5132    moreSpecialOptions_ = saveMoreSpecialOptions;
5133    return ;
5134}
5135
5136
5137// Solve the initial LP relaxation
5138void
5139CbcModel::initialSolve()
5140{
5141    assert (solver_);
5142    // Double check optimization directions line up
5143    dblParam_[CbcOptimizationDirection] = solver_->getObjSense();
5144    // Check if bounds are all integral (as may get messed up later)
5145    checkModel();
5146    if (!solverCharacteristics_) {
5147        OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
5148        if (solverCharacteristics) {
5149            solverCharacteristics_ = solverCharacteristics;
5150        } else {
5151            // replace in solver
5152            OsiBabSolver defaultC;
5153            solver_->setAuxiliaryInfo(&defaultC);
5154            solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
5155        }
5156    }
5157    solverCharacteristics_->setSolver(solver_);
5158    solver_->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
5159    solver_->initialSolve();
5160    solver_->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo, NULL) ;
5161    if (!solver_->isProvenOptimal())
5162        solver_->resolve();
5163    // But set up so Jon Lee will be happy
5164    status_ = -1;
5165    secondaryStatus_ = -1;
5166    originalContinuousObjective_ = solver_->getObjValue() * solver_->getObjSense();
5167    bestPossibleObjective_ = originalContinuousObjective_;
5168    if (solver_->isProvenDualInfeasible())
5169      originalContinuousObjective_ = -COIN_DBL_MAX;
5170    delete [] continuousSolution_;
5171    continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
5172                                          solver_->getNumCols());
5173    setPointers(solver_);
5174    solverCharacteristics_ = NULL;
5175}
5176
5177/*! \brief Get an empty basis object
5178
5179  Return an empty CoinWarmStartBasis object with the requested capacity,
5180  appropriate for the current solver. The object is cloned from the object
5181  cached as emptyWarmStart_. If there is no cached object, the routine
5182  queries the solver for a warm start object, empties it, and caches the
5183  result.
5184*/
5185
5186CoinWarmStartBasis *CbcModel::getEmptyBasis (int ns, int na) const
5187
5188{
5189    CoinWarmStartBasis *emptyBasis ;
5190    /*
5191      Acquire an empty basis object, if we don't yet have one.
5192    */
5193    if (emptyWarmStart_ == 0) {
5194        if (solver_ == 0) {
5195            throw CoinError("Cannot construct basis without solver!",
5196                            "getEmptyBasis", "CbcModel") ;
5197        }
5198        emptyBasis =
5199            dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart()) ;
5200        if (emptyBasis == 0) {
5201            throw CoinError(
5202                "Solver does not appear to use a basis-oriented warm start.",
5203                "getEmptyBasis", "CbcModel") ;
5204        }
5205        emptyBasis->setSize(0, 0) ;
5206        emptyWarmStart_ = dynamic_cast<CoinWarmStart *>(emptyBasis) ;
5207    }
5208    /*
5209      Clone the empty basis object, resize it as requested, and return.
5210    */
5211    emptyBasis = dynamic_cast<CoinWarmStartBasis *>(emptyWarmStart_->clone()) ;
5212    assert(emptyBasis) ;
5213    if (ns != 0 || na != 0) emptyBasis->setSize(ns, na) ;
5214
5215    return (emptyBasis) ;
5216}
5217
5218
5219/** Default Constructor
5220
5221  Creates an empty model without an associated solver.
5222*/
5223CbcModel::CbcModel()
5224
5225        :
5226        solver_(NULL),
5227        ownership_(0x80000000),
5228        continuousSolver_(NULL),
5229        referenceSolver_(NULL),
5230        defaultHandler_(true),
5231        emptyWarmStart_(NULL),
5232        bestObjective_(COIN_DBL_MAX),
5233        bestPossibleObjective_(COIN_DBL_MAX),
5234        sumChangeObjective1_(0.0),
5235        sumChangeObjective2_(0.0),
5236        bestSolution_(NULL),
5237        savedSolutions_(NULL),
5238        currentSolution_(NULL),
5239        testSolution_(NULL),
5240        globalConflictCuts_(NULL