source: stable/2.6/Cbc/src/CbcModel.cpp @ 1627

Last change on this file since 1627 was 1627, checked in by tkr, 8 years ago

Merging r1623 from trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 648.9 KB
RevLine 
[1271]1/* $Id: CbcModel.cpp 1627 2011-03-26 15:13:57Z tkr $ */
[2]2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4#if defined(_MSC_VER)
5// Turn off compiler warning about long names
6#  pragma warning(disable:4786)
7#endif
[311]8
[325]9#include "CbcConfig.h"
[311]10
[2]11#include <string>
12//#define CBC_DEBUG 1
13//#define CHECK_CUT_COUNTS
[1412]14//#define CHECK_NODE
[2]15//#define CHECK_NODE_FULL
[264]16//#define NODE_LOG
[640]17//#define GLOBAL_CUTS_JUST_POINTERS
[1271]18#ifdef CGL_DEBUG_GOMORY
19extern int gomory_try;
[687]20#endif
[2]21#include <cassert>
22#include <cmath>
23#include <cfloat>
[277]24
[311]25#ifdef COIN_HAS_CLP
[208]26// include Presolve from Clp
27#include "ClpPresolve.hpp"
28#include "OsiClpSolverInterface.hpp"
[931]29#include "ClpNode.hpp"
30#include "ClpDualRowDantzig.hpp"
31#include "ClpSimplexPrimal.hpp"
[277]32#endif
33
34#include "CbcEventHandler.hpp"
[2]35
36#include "OsiSolverInterface.hpp"
[264]37#include "OsiAuxInfo.hpp"
[222]38#include "OsiSolverBranch.hpp"
[640]39#include "OsiChooseVariable.hpp"
[2]40#include "CoinWarmStartBasis.hpp"
41#include "CoinPackedMatrix.hpp"
[140]42#include "CoinHelperFunctions.hpp"
[2]43#include "CbcBranchActual.hpp"
[135]44#include "CbcBranchDynamic.hpp"
[2]45#include "CbcHeuristic.hpp"
[687]46#include "CbcHeuristicFPump.hpp"
[1315]47#include "CbcHeuristicRINS.hpp"
[1271]48#include "CbcHeuristicDive.hpp"
[2]49#include "CbcModel.hpp"
[186]50#include "CbcTreeLocal.hpp"
[150]51#include "CbcStatistics.hpp"
[97]52#include "CbcStrategy.hpp"
[2]53#include "CbcMessage.hpp"
54#include "OsiRowCut.hpp"
55#include "OsiColCut.hpp"
56#include "OsiRowCutDebugger.hpp"
57#include "OsiCuts.hpp"
58#include "CbcCountRowCut.hpp"
59#include "CbcCutGenerator.hpp"
[165]60#include "CbcFeasibilityBase.hpp"
[687]61#include "CbcFathom.hpp"
[2]62// include Probing
63#include "CglProbing.hpp"
[1271]64#include "CglGomory.hpp"
65#include "CglTwomir.hpp"
[222]66// include preprocessing
67#include "CglPreProcess.hpp"
[640]68#include "CglDuplicateRow.hpp"
69#include "CglStored.hpp"
[833]70#include "CglClique.hpp"
[2]71
72#include "CoinTime.hpp"
[640]73#include "CoinMpsIO.hpp"
[2]74
75#include "CbcCompareActual.hpp"
76#include "CbcTree.hpp"
[1409]77// This may be dummy
78#include "CbcThread.hpp"
[2]79/* Various functions local to CbcModel.cpp */
80
81namespace {
82
83//-------------------------------------------------------------------
[1286]84// Returns the greatest common denominator of two
85// positive integers, a and b, found using Euclid's algorithm
[2]86//-------------------------------------------------------------------
[1286]87static int gcd(int a, int b)
[2]88{
[1286]89    int remainder = -1;
90    // make sure a<=b (will always remain so)
91    if (a > b) {
92        // Swap a and b
93        int temp = a;
94        a = b;
95        b = temp;
[2]96    }
[1286]97    // if zero then gcd is nonzero (zero may occur in rhs of packed)
98    if (!a) {
99        if (b) {
100            return b;
101        } else {
102            printf("**** gcd given two zeros!!\n");
103            abort();
104        }
105    }
106    while (remainder) {
107        remainder = b % a;
108        b = a;
109        a = remainder;
110    }
111    return b;
[2]112}
113
[931]114
[2]115
116#ifdef CHECK_NODE_FULL
117
118/*
119  Routine to verify that tree linkage is correct. The invariant that is tested
120  is
121
122  reference count = (number of actual references) + (number of branches left)
123
124  The routine builds a set of paired arrays, info and count, by traversing the
125  tree. Each CbcNodeInfo is recorded in info, and the number of times it is
126  referenced (via the parent field) is recorded in count. Then a final check is
127  made to see if the numberPointingToThis_ field agrees.
128*/
129
130void verifyTreeNodes (const CbcTree * branchingTree, const CbcModel &model)
131
[1286]132{
133    if (model.getNodeCount() == 661) return;
134    printf("*** CHECKING tree after %d nodes\n", model.getNodeCount()) ;
135
136    int j ;
137    int nNodes = branchingTree->size() ;
[2]138# define MAXINFO 1000
[1286]139    int *count = new int [MAXINFO] ;
140    CbcNodeInfo **info = new CbcNodeInfo*[MAXINFO] ;
141    int nInfo = 0 ;
142    /*
143      Collect all CbcNodeInfo objects in info, by starting from each live node and
144      traversing back to the root. Nodes in the live set should have unexplored
145      branches remaining.
[2]146
[1286]147      TODO: The `while (nodeInfo)' loop could be made to break on reaching a
148        common ancester (nodeInfo is found in info[k]). Alternatively, the
149        check could change to signal an error if nodeInfo is not found above a
150        common ancestor.
151    */
152    for (j = 0 ; j < nNodes ; j++) {
153        CbcNode *node = branchingTree->nodePointer(j) ;
154        if (!node)
155            continue;
156        CbcNodeInfo *nodeInfo = node->nodeInfo() ;
157        int change = node->nodeInfo()->numberBranchesLeft() ;
158        assert(change) ;
159        while (nodeInfo) {
160            int k ;
161            for (k = 0 ; k < nInfo ; k++) {
162                if (nodeInfo == info[k]) break ;
163            }
164            if (k == nInfo) {
165                assert(nInfo < MAXINFO) ;
166                nInfo++ ;
167                info[k] = nodeInfo ;
168                count[k] = 0 ;
169            }
170            nodeInfo = nodeInfo->parent() ;
171        }
172    }
173    /*
174      Walk the info array. For each nodeInfo, look up its parent in info and
175      increment the corresponding count.
176    */
177    for (j = 0 ; j < nInfo ; j++) {
178        CbcNodeInfo *nodeInfo = info[j] ;
179        nodeInfo = nodeInfo->parent() ;
180        if (nodeInfo) {
181            int k ;
182            for (k = 0 ; k < nInfo ; k++) {
183                if (nodeInfo == info[k]) break ;
184            }
185            assert (k < nInfo) ;
186            count[k]++ ;
187        }
188    }
189    /*
190      Walk the info array one more time and check that the invariant holds. The
191      number of references (numberPointingToThis()) should equal the sum of the
192      number of actual references (held in count[]) plus the number of potential
193      references (unexplored branches, numberBranchesLeft()).
194    */
195    for (j = 0; j < nInfo; j++) {
196        CbcNodeInfo * nodeInfo = info[j] ;
197        if (nodeInfo) {
198            int k ;
199            for (k = 0; k < nInfo; k++)
200                if (nodeInfo == info[k])
201                    break ;
202            printf("Nodeinfo %x - %d left, %d count\n",
203                   nodeInfo,
204                   nodeInfo->numberBranchesLeft(),
205                   nodeInfo->numberPointingToThis()) ;
206            assert(nodeInfo->numberPointingToThis() ==
207                   count[k] + nodeInfo->numberBranchesLeft()) ;
208        }
209    }
[2]210
[1286]211    delete [] count ;
212    delete [] info ;
[2]213
[1286]214    return ;
215}
216
[36]217#endif  /* CHECK_NODE_FULL */
[2]218
[931]219
[2]220
221#ifdef CHECK_CUT_COUNTS
222
223/*
224  Routine to verify that cut reference counts are correct.
225*/
226void verifyCutCounts (const CbcTree * branchingTree, CbcModel &model)
227
[1286]228{
229    printf("*** CHECKING cuts after %d nodes\n", model.getNodeCount()) ;
[2]230
[1286]231    int j ;
232    int nNodes = branchingTree->size() ;
[2]233
[1286]234    /*
235      cut.tempNumber_ exists for the purpose of doing this verification. Clear it
236      in all cuts. We traverse the tree by starting from each live node and working
237      back to the root. At each CbcNodeInfo, check for cuts.
238    */
239    for (j = 0 ; j < nNodes ; j++) {
240        CbcNode *node = branchingTree->nodePointer(j) ;
241        CbcNodeInfo * nodeInfo = node->nodeInfo() ;
242        assert (node->nodeInfo()->numberBranchesLeft()) ;
243        while (nodeInfo) {
244            int k ;
245            for (k = 0 ; k < nodeInfo->numberCuts() ; k++) {
246                CbcCountRowCut *cut = nodeInfo->cuts()[k] ;
247                if (cut) cut->tempNumber_ = 0;
248            }
249            nodeInfo = nodeInfo->parent() ;
250        }
251    }
252    /*
253      Walk the live set again, this time collecting the list of cuts in use at each
254      node. addCuts1 will collect the cuts in model.addedCuts_. Take into account
255      that when we recreate the basis for a node, we compress out the slack cuts.
256    */
257    for (j = 0 ; j < nNodes ; j++) {
258        CoinWarmStartBasis *debugws = model.getEmptyBasis() ;
259        CbcNode *node = branchingTree->nodePointer(j) ;
260        CbcNodeInfo *nodeInfo = node->nodeInfo();
261        int change = node->nodeInfo()->numberBranchesLeft() ;
262        printf("Node %d %x (info %x) var %d way %d obj %g", j, node,
263               node->nodeInfo(), node->columnNumber(), node->way(),
264               node->objectiveValue()) ;
[2]265
[1286]266        model.addCuts1(node, debugws) ;
[2]267
[1286]268        int i ;
269        int numberRowsAtContinuous = model.numberRowsAtContinuous() ;
270        CbcCountRowCut **addedCuts = model.addedCuts() ;
271        for (i = 0 ; i < model.currentNumberCuts() ; i++) {
272            CoinWarmStartBasis::Status status =
273                debugws->getArtifStatus(i + numberRowsAtContinuous) ;
274            if (status != CoinWarmStartBasis::basic && addedCuts[i]) {
275                addedCuts[i]->tempNumber_ += change ;
276            }
277        }
[2]278
[1286]279        while (nodeInfo) {
280            nodeInfo = nodeInfo->parent() ;
281            if (nodeInfo) printf(" -> %x", nodeInfo);
282        }
283        printf("\n") ;
284        delete debugws ;
285    }
286    /*
287      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.
[2]288
[1286]289      TODO: Rewrite to check and print mismatch only when tempNumber_ == 0?
290    */
291    for (j = 0 ; j < nNodes ; j++) {
292        CbcNode *node = branchingTree->nodePointer(j) ;
293        CbcNodeInfo *nodeInfo = node->nodeInfo();
294        while (nodeInfo) {
295            int k ;
296            for (k = 0 ; k < nodeInfo->numberCuts() ; k++) {
297                CbcCountRowCut *cut = nodeInfo->cuts()[k] ;
298                if (cut && cut->tempNumber_ >= 0) {
299                    if (cut->tempNumber_ != cut->numberPointingToThis())
300                        printf("mismatch %x %d %x %d %d\n", nodeInfo, k,
301                               cut, cut->tempNumber_, cut->numberPointingToThis()) ;
302                    else
303                        printf("   match %x %d %x %d %d\n", nodeInfo, k,
304                               cut, cut->tempNumber_, cut->numberPointingToThis()) ;
305                    cut->tempNumber_ = -1 ;
306                }
307            }
308            nodeInfo = nodeInfo->parent() ;
309        }
310    }
[2]311
[1286]312    return ;
313}
[2]314
315#endif /* CHECK_CUT_COUNTS */
316
[931]317
[640]318#ifdef CHECK_CUT_SIZE
319
320/*
321  Routine to verify that cut reference counts are correct.
322*/
323void verifyCutSize (const CbcTree * branchingTree, CbcModel &model)
[1286]324{
[640]325
[1286]326    int j ;
327    int nNodes = branchingTree->size() ;
328    int totalCuts = 0;
[640]329
[1286]330    /*
331      cut.tempNumber_ exists for the purpose of doing this verification. Clear it
332      in all cuts. We traverse the tree by starting from each live node and working
333      back to the root. At each CbcNodeInfo, check for cuts.
334    */
335    for (j = 0 ; j < nNodes ; j++) {
336        CbcNode *node = branchingTree->nodePointer(j) ;
337        CbcNodeInfo * nodeInfo = node->nodeInfo() ;
338        assert (node->nodeInfo()->numberBranchesLeft()) ;
339        while (nodeInfo) {
340            totalCuts += nodeInfo->numberCuts();
341            nodeInfo = nodeInfo->parent() ;
342        }
[640]343    }
[1286]344    printf("*** CHECKING cuts (size) after %d nodes - %d cuts\n", model.getNodeCount(), totalCuts) ;
345    return ;
[2]346}
347
[640]348#endif /* CHECK_CUT_SIZE */
349
350}
351
[1286]352/* End unnamed namespace for CbcModel.cpp */
[2]353
[931]354
[1286]355void
[2]356CbcModel::analyzeObjective ()
357/*
358  Try to find a minimum change in the objective function. The first scan
359  checks that there are no continuous variables with non-zero coefficients,
360  and grabs the largest objective coefficient associated with an unfixed
361  integer variable. The second scan attempts to scale up the objective
362  coefficients to a point where they are sufficiently close to integer that
363  we can pretend they are integer, and calculate a gcd over the coefficients
364  of interest. This will be the minimum increment for the scaled coefficients.
365  The final action is to scale the increment back for the original coefficients
366  and install it, if it's better than the existing value.
367
368  John's note: We could do better than this.
369
370  John's second note - apologies for changing s to z
371*/
[1286]372{
373    const double *objective = getObjCoefficients() ;
374    const double *lower = getColLower() ;
375    const double *upper = getColUpper() ;
376    /*
377      Scan continuous and integer variables to see if continuous
378      are cover or network with integral rhs.
379    */
380    double continuousMultiplier = 1.0;
381    double * coeffMultiplier = NULL;
382    double largestObj = 0.0;
383    double smallestObj = COIN_DBL_MAX;
384    {
385        const double *rowLower = getRowLower() ;
386        const double *rowUpper = getRowUpper() ;
387        int numberRows = solver_->getNumRows() ;
388        double * rhs = new double [numberRows];
389        memset(rhs, 0, numberRows*sizeof(double));
390        int iColumn;
391        int numberColumns = solver_->getNumCols() ;
392        // Column copy of matrix
393        bool allPlusOnes = true;
394        bool allOnes = true;
395        int problemType = -1;
396        const double * element = solver_->getMatrixByCol()->getElements();
397        const int * row = solver_->getMatrixByCol()->getIndices();
398        const CoinBigIndex * columnStart = solver_->getMatrixByCol()->getVectorStarts();
399        const int * columnLength = solver_->getMatrixByCol()->getVectorLengths();
400        int numberInteger = 0;
401        int numberIntegerObj = 0;
402        int numberGeneralIntegerObj = 0;
403        int numberIntegerWeight = 0;
404        int numberContinuousObj = 0;
405        double cost = COIN_DBL_MAX;
406        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
407            if (upper[iColumn] == lower[iColumn]) {
408                CoinBigIndex start = columnStart[iColumn];
409                CoinBigIndex end = start + columnLength[iColumn];
410                for (CoinBigIndex j = start; j < end; j++) {
411                    int iRow = row[j];
412                    rhs[iRow] += lower[iColumn] * element[j];
413                }
414            } else {
415                double objValue = objective[iColumn];
416                if (solver_->isInteger(iColumn))
417                    numberInteger++;
418                if (objValue) {
419                    if (!solver_->isInteger(iColumn)) {
420                        numberContinuousObj++;
421                    } else {
422                        largestObj = CoinMax(largestObj, fabs(objValue));
423                        smallestObj = CoinMin(smallestObj, fabs(objValue));
424                        numberIntegerObj++;
425                        if (cost == COIN_DBL_MAX)
426                            cost = objValue;
427                        else if (cost != objValue)
428                            cost = -COIN_DBL_MAX;
429                        int gap = static_cast<int> (upper[iColumn] - lower[iColumn]);
430                        if (gap > 1) {
431                            numberGeneralIntegerObj++;
432                            numberIntegerWeight += gap;
433                        }
434                    }
435                }
436            }
437        }
438        int iType = 0;
439        if (!numberContinuousObj && numberIntegerObj <= 5 && numberIntegerWeight <= 100 &&
440                numberIntegerObj*3 < numberObjects_ && !parentModel_ && solver_->getNumRows() > 100)
441            iType = 3 + 4;
442        else if (!numberContinuousObj && numberIntegerObj <= 100 &&
443                 numberIntegerObj*5 < numberObjects_ && numberIntegerWeight <= 100 &&
444                 !parentModel_ &&
445                 solver_->getNumRows() > 100 && cost != -COIN_DBL_MAX)
446            iType = 2 + 4;
447        else if (!numberContinuousObj && numberIntegerObj <= 100 &&
448                 numberIntegerObj*5 < numberObjects_ &&
449                 !parentModel_ &&
450                 solver_->getNumRows() > 100 && cost != -COIN_DBL_MAX)
451            iType = 8;
452        int iTest = getMaximumNodes();
453        if (iTest >= 987654320 && iTest < 987654330 && numberObjects_ && !parentModel_) {
454            iType = iTest - 987654320;
455            printf("Testing %d integer variables out of %d objects (%d integer) have cost of %g - %d continuous\n",
456                   numberIntegerObj, numberObjects_, numberInteger, cost, numberContinuousObj);
457            if (iType == 9)
458                exit(77);
459            if (numberContinuousObj)
460                iType = 0;
461        }
462
463        //if (!numberContinuousObj&&(numberIntegerObj<=5||cost!=-COIN_DBL_MAX)&&
464        //numberIntegerObj*3<numberObjects_&&!parentModel_&&solver_->getNumRows()>100) {
465        if (iType) {
466            /*
467            A) put high priority on (if none)
468            B) create artificial objective (if clp)
469            */
470            int iPriority = -1;
471            for (int i = 0; i < numberObjects_; i++) {
472                int k = object_[i]->priority();
473                if (iPriority == -1)
474                    iPriority = k;
475                else if (iPriority != k)
476                    iPriority = -2;
477            }
478            bool branchOnSatisfied = ((iType & 1) != 0);
479            bool createFake = ((iType & 2) != 0);
480            bool randomCost = ((iType & 4) != 0);
481            if (iPriority >= 0) {
482                char general[200];
483                if (cost == -COIN_DBL_MAX) {
484                    sprintf(general, "%d integer variables out of %d objects (%d integer) have costs - high priority",
485                            numberIntegerObj, numberObjects_, numberInteger);
486                } else if (cost == COIN_DBL_MAX) {
487                    sprintf(general, "No integer variables out of %d objects (%d integer) have costs",
488                            numberObjects_, numberInteger);
489                    branchOnSatisfied = false;
490                } else {
491                    sprintf(general, "%d integer variables out of %d objects (%d integer) have cost of %g - high priority",
492                            numberIntegerObj, numberObjects_, numberInteger, cost);
493                }
494                messageHandler()->message(CBC_GENERAL,
495                                          messages())
496                << general << CoinMessageEol ;
497                sprintf(general, "branch on satisfied %c create fake objective %c random cost %c",
498                        branchOnSatisfied ? 'Y' : 'N',
499                        createFake ? 'Y' : 'N',
500                        randomCost ? 'Y' : 'N');
501                messageHandler()->message(CBC_GENERAL,
502                                          messages())
503                << general << CoinMessageEol ;
504                // switch off clp type branching
505                fastNodeDepth_ = -1;
506                int highPriority = (branchOnSatisfied) ? -999 : 100;
507                for (int i = 0; i < numberObjects_; i++) {
508                    CbcSimpleInteger * thisOne = dynamic_cast <CbcSimpleInteger *> (object_[i]);
509                    object_[i]->setPriority(1000);
510                    if (thisOne) {
511                        int iColumn = thisOne->columnNumber();
512                        if (objective[iColumn])
513                            thisOne->setPriority(highPriority);
514                    }
515                }
516            }
[1088]517#ifdef COIN_HAS_CLP
[1286]518            OsiClpSolverInterface * clpSolver
519            = dynamic_cast<OsiClpSolverInterface *> (solver_);
520            if (clpSolver && createFake) {
521                // Create artificial objective to be used when all else fixed
522                int numberColumns = clpSolver->getNumCols();
523                double * fakeObj = new double [numberColumns];
524                // Column copy
525                const CoinPackedMatrix  * matrixByCol = clpSolver->getMatrixByCol();
526                //const double * element = matrixByCol.getElements();
527                //const int * row = matrixByCol.getIndices();
528                //const CoinBigIndex * columnStart = matrixByCol.getVectorStarts();
529                const int * columnLength = matrixByCol->getVectorLengths();
530                const double * solution = clpSolver->getColSolution();
[1393]531#ifdef JJF_ZERO
[1286]532                int nAtBound = 0;
533                for (int i = 0; i < numberColumns; i++) {
534                    double lowerValue = lower[i];
535                    double upperValue = upper[i];
536                    if (clpSolver->isInteger(i)) {
537                        double lowerValue = lower[i];
538                        double upperValue = upper[i];
539                        double value = solution[i];
540                        if (value < lowerValue + 1.0e-6 ||
541                                value > upperValue - 1.0e-6)
542                            nAtBound++;
543                    }
544                }
[1271]545#endif
[1409]546                /*
547                  Generate a random objective function for problems where the given objective
548                  function is not terribly useful. (Nearly feasible, single integer variable,
549                  that sort of thing.
550                */
[1286]551                CoinDrand48(true, 1234567);
552                for (int i = 0; i < numberColumns; i++) {
553                    double lowerValue = lower[i];
554                    double upperValue = upper[i];
555                    double value = (randomCost) ? ceil((CoinDrand48() + 0.5) * 1000)
556                                   : i + 1 + columnLength[i] * 1000;
557                    value *= 0.001;
558                    //value += columnLength[i];
559                    if (lowerValue > -1.0e5 || upperValue < 1.0e5) {
560                        if (fabs(lowerValue) > fabs(upperValue))
561                            value = - value;
562                        if (clpSolver->isInteger(i)) {
563                            double solValue = solution[i];
564                            // Better to add in 0.5 or 1.0??
565                            if (solValue < lowerValue + 1.0e-6)
566                                value = fabs(value) + 0.5; //fabs(value*1.5);
567                            else if (solValue > upperValue - 1.0e-6)
568                                value = -fabs(value) - 0.5; //-fabs(value*1.5);
569                        }
570                    } else {
571                        value = 0.0;
572                    }
573                    fakeObj[i] = value;
574                }
575                // pass to solver
576                clpSolver->setFakeObjective(fakeObj);
577                delete [] fakeObj;
578            }
[1088]579#endif
[1286]580        } else if (largestObj < smallestObj*5.0 && !parentModel_ &&
581                   !numberContinuousObj &&
582                   !numberGeneralIntegerObj &&
583                   numberIntegerObj*2 < numberColumns) {
584            // up priorities on costed
585            int iPriority = -1;
586            for (int i = 0; i < numberObjects_; i++) {
587                int k = object_[i]->priority();
588                if (iPriority == -1)
589                    iPriority = k;
590                else if (iPriority != k)
591                    iPriority = -2;
592            }
593            if (iPriority >= 100) {
[1271]594#ifdef CLP_INVESTIGATE
[1286]595                printf("Setting variables with obj to high priority\n");
[1271]596#endif
[1286]597                for (int i = 0; i < numberObjects_; i++) {
598                    CbcSimpleInteger * obj =
599                        dynamic_cast <CbcSimpleInteger *>(object_[i]) ;
600                    if (obj) {
601                        int iColumn = obj->columnNumber();
602                        if (objective[iColumn])
603                            object_[i]->setPriority(iPriority - 1);
604                    }
605                }
606            }
607        }
608        int iRow;
609        for (iRow = 0; iRow < numberRows; iRow++) {
610            if (rowLower[iRow] > -1.0e20 &&
611                    fabs(rowLower[iRow] - rhs[iRow] - floor(rowLower[iRow] - rhs[iRow] + 0.5)) > 1.0e-10) {
612                continuousMultiplier = 0.0;
613                break;
614            }
615            if (rowUpper[iRow] < 1.0e20 &&
616                    fabs(rowUpper[iRow] - rhs[iRow] - floor(rowUpper[iRow] - rhs[iRow] + 0.5)) > 1.0e-10) {
617                continuousMultiplier = 0.0;
618                break;
619            }
620            // set rhs to limiting value
621            if (rowLower[iRow] != rowUpper[iRow]) {
622                if (rowLower[iRow] > -1.0e20) {
623                    if (rowUpper[iRow] < 1.0e20) {
624                        // no good
625                        continuousMultiplier = 0.0;
626                        break;
627                    } else {
628                        rhs[iRow] = rowLower[iRow] - rhs[iRow];
629                        if (problemType < 0)
630                            problemType = 3; // set cover
631                        else if (problemType != 3)
632                            problemType = 4;
633                    }
634                } else {
635                    rhs[iRow] = rowUpper[iRow] - rhs[iRow];
636                    if (problemType < 0)
637                        problemType = 1; // set partitioning <=
638                    else if (problemType != 1)
639                        problemType = 4;
640                }
641            } else {
642                rhs[iRow] = rowUpper[iRow] - rhs[iRow];
643                if (problemType < 0)
644                    problemType = 3; // set partitioning ==
645                else if (problemType != 2)
646                    problemType = 2;
647            }
648            if (fabs(rhs[iRow] - 1.0) > 1.0e-12)
649                problemType = 4;
650        }
651        if (continuousMultiplier) {
652            // 1 network, 2 cover, 4 negative cover
653            int possible = 7;
654            bool unitRhs = true;
655            // See which rows could be set cover
656            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
657                if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
658                    CoinBigIndex start = columnStart[iColumn];
659                    CoinBigIndex end = start + columnLength[iColumn];
660                    for (CoinBigIndex j = start; j < end; j++) {
661                        double value = element[j];
662                        if (value == 1.0) {
663                        } else if (value == -1.0) {
664                            rhs[row[j]] = -0.5;
665                            allPlusOnes = false;
666                        } else {
667                            rhs[row[j]] = -COIN_DBL_MAX;
668                            allOnes = false;
669                        }
670                    }
671                }
672            }
673            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
674                if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
675                    if (!isInteger(iColumn)) {
676                        CoinBigIndex start = columnStart[iColumn];
677                        CoinBigIndex end = start + columnLength[iColumn];
678                        double rhsValue = 0.0;
679                        // 1 all ones, -1 all -1s, 2 all +- 1, 3 no good
680                        int type = 0;
681                        for (CoinBigIndex j = start; j < end; j++) {
682                            double value = element[j];
683                            if (fabs(value) != 1.0) {
684                                type = 3;
685                                break;
686                            } else if (value == 1.0) {
687                                if (!type)
688                                    type = 1;
689                                else if (type != 1)
690                                    type = 2;
691                            } else {
692                                if (!type)
693                                    type = -1;
694                                else if (type != -1)
695                                    type = 2;
696                            }
697                            int iRow = row[j];
698                            if (rhs[iRow] == -COIN_DBL_MAX) {
699                                type = 3;
700                                break;
701                            } else if (rhs[iRow] == -0.5) {
702                                // different values
703                                unitRhs = false;
704                            } else if (rhsValue) {
705                                if (rhsValue != rhs[iRow])
706                                    unitRhs = false;
707                            } else {
708                                rhsValue = rhs[iRow];
709                            }
710                        }
711                        // if no elements OK
712                        if (type == 3) {
713                            // no good
714                            possible = 0;
715                            break;
716                        } else if (type == 2) {
717                            if (end - start > 2) {
718                                // no good
719                                possible = 0;
720                                break;
721                            } else {
722                                // only network
723                                possible &= 1;
724                                if (!possible)
725                                    break;
726                            }
727                        } else if (type == 1) {
728                            // only cover
729                            possible &= 2;
730                            if (!possible)
731                                break;
732                        } else if (type == -1) {
733                            // only negative cover
734                            possible &= 4;
735                            if (!possible)
736                                break;
737                        }
738                    }
739                }
740            }
741            if ((possible == 2 || possible == 4) && !unitRhs) {
[931]742#if COIN_DEVELOP>1
[1286]743                printf("XXXXXX Continuous all +1 but different rhs\n");
[810]744#endif
[1286]745                possible = 0;
746            }
747            // may be all integer
748            if (possible != 7) {
749                if (!possible)
750                    continuousMultiplier = 0.0;
751                else if (possible == 1)
752                    continuousMultiplier = 1.0;
753                else
754                    continuousMultiplier = 0.0; // 0.5 was incorrect;
[931]755#if COIN_DEVELOP>1
[1286]756                if (continuousMultiplier)
757                    printf("XXXXXX multiplier of %g\n", continuousMultiplier);
[810]758#endif
[1286]759                if (continuousMultiplier == 0.5) {
760                    coeffMultiplier = new double [numberColumns];
761                    bool allOne = true;
762                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
763                        coeffMultiplier[iColumn] = 1.0;
764                        if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
765                            if (!isInteger(iColumn)) {
766                                CoinBigIndex start = columnStart[iColumn];
767                                int iRow = row[start];
768                                double value = rhs[iRow];
769                                assert (value >= 0.0);
770                                if (value != 0.0 && value != 1.0)
771                                    allOne = false;
772                                coeffMultiplier[iColumn] = 0.5 * value;
773                            }
774                        }
775                    }
776                    if (allOne) {
777                        // back to old way
778                        delete [] coeffMultiplier;
779                        coeffMultiplier = NULL;
780                    }
781                }
782            } else {
783                // all integer
784                problemType_ = problemType;
[931]785#if COIN_DEVELOP>1
[1286]786                printf("Problem type is %d\n", problemType_);
[814]787#endif
[1286]788            }
789        }
790
791        // But try again
792        if (continuousMultiplier < 1.0) {
793            memset(rhs, 0, numberRows*sizeof(double));
794            int * count = new int [numberRows];
795            memset(count, 0, numberRows*sizeof(int));
796            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
797                CoinBigIndex start = columnStart[iColumn];
798                CoinBigIndex end = start + columnLength[iColumn];
799                if (upper[iColumn] == lower[iColumn]) {
800                    for (CoinBigIndex j = start; j < end; j++) {
801                        int iRow = row[j];
802                        rhs[iRow] += lower[iColumn] * element[j];
803                    }
804                } else if (solver_->isInteger(iColumn)) {
805                    for (CoinBigIndex j = start; j < end; j++) {
806                        int iRow = row[j];
807                        if (fabs(element[j] - floor(element[j] + 0.5)) > 1.0e-10)
808                            rhs[iRow]  = COIN_DBL_MAX;
809                    }
810                } else {
811                    for (CoinBigIndex j = start; j < end; j++) {
812                        int iRow = row[j];
813                        count[iRow]++;
814                        if (fabs(element[j]) != 1.0)
815                            rhs[iRow]  = COIN_DBL_MAX;
816                    }
817                }
818            }
819            // now look at continuous
820            bool allGood = true;
821            double direction = solver_->getObjSense() ;
822            int numberObj = 0;
823            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
824                if (upper[iColumn] > lower[iColumn]) {
825                    double objValue = objective[iColumn] * direction;
826                    if (objValue && !solver_->isInteger(iColumn)) {
827                        numberObj++;
828                        CoinBigIndex start = columnStart[iColumn];
829                        CoinBigIndex end = start + columnLength[iColumn];
830                        if (objValue > 0.0) {
831                            // wants to be as low as possible
832                            if (lower[iColumn] < -1.0e10 || fabs(lower[iColumn] - floor(lower[iColumn] + 0.5)) > 1.0e-10) {
833                                allGood = false;
834                                break;
835                            } else if (upper[iColumn] < 1.0e10 && fabs(upper[iColumn] - floor(upper[iColumn] + 0.5)) > 1.0e-10) {
836                                allGood = false;
837                                break;
838                            }
839                            bool singletonRow = true;
840                            bool equality = false;
841                            for (CoinBigIndex j = start; j < end; j++) {
842                                int iRow = row[j];
843                                if (count[iRow] > 1)
844                                    singletonRow = false;
845                                else if (rowLower[iRow] == rowUpper[iRow])
846                                    equality = true;
847                                double rhsValue = rhs[iRow];
848                                double lowerValue = rowLower[iRow];
849                                double upperValue = rowUpper[iRow];
850                                if (rhsValue < 1.0e20) {
851                                    if (lowerValue > -1.0e20)
852                                        lowerValue -= rhsValue;
853                                    if (upperValue < 1.0e20)
854                                        upperValue -= rhsValue;
855                                }
856                                if (fabs(rhsValue) > 1.0e20 || fabs(rhsValue - floor(rhsValue + 0.5)) > 1.0e-10
857                                        || fabs(element[j]) != 1.0) {
858                                    // no good
859                                    allGood = false;
860                                    break;
861                                }
862                                if (element[j] > 0.0) {
863                                    if (lowerValue > -1.0e20 && fabs(lowerValue - floor(lowerValue + 0.5)) > 1.0e-10) {
864                                        // no good
865                                        allGood = false;
866                                        break;
867                                    }
868                                } else {
869                                    if (upperValue < 1.0e20 && fabs(upperValue - floor(upperValue + 0.5)) > 1.0e-10) {
870                                        // no good
871                                        allGood = false;
872                                        break;
873                                    }
874                                }
875                            }
876                            if (!singletonRow && end > start + 1 && !equality)
877                                allGood = false;
878                            if (!allGood)
879                                break;
880                        } else {
881                            // wants to be as high as possible
882                            if (upper[iColumn] > 1.0e10 || fabs(upper[iColumn] - floor(upper[iColumn] + 0.5)) > 1.0e-10) {
883                                allGood = false;
884                                break;
885                            } else if (lower[iColumn] > -1.0e10 && fabs(lower[iColumn] - floor(lower[iColumn] + 0.5)) > 1.0e-10) {
886                                allGood = false;
887                                break;
888                            }
889                            bool singletonRow = true;
890                            bool equality = false;
891                            for (CoinBigIndex j = start; j < end; j++) {
892                                int iRow = row[j];
893                                if (count[iRow] > 1)
894                                    singletonRow = false;
895                                else if (rowLower[iRow] == rowUpper[iRow])
896                                    equality = true;
897                                double rhsValue = rhs[iRow];
898                                double lowerValue = rowLower[iRow];
899                                double upperValue = rowUpper[iRow];
900                                if (rhsValue < 1.0e20) {
901                                    if (lowerValue > -1.0e20)
902                                        lowerValue -= rhsValue;
903                                    if (upperValue < 1.0e20)
904                                        upperValue -= rhsValue;
905                                }
906                                if (fabs(rhsValue) > 1.0e20 || fabs(rhsValue - floor(rhsValue + 0.5)) > 1.0e-10
907                                        || fabs(element[j]) != 1.0) {
908                                    // no good
909                                    allGood = false;
910                                    break;
911                                }
912                                if (element[j] < 0.0) {
913                                    if (lowerValue > -1.0e20 && fabs(lowerValue - floor(lowerValue + 0.5)) > 1.0e-10) {
914                                        // no good
915                                        allGood = false;
916                                        break;
917                                    }
918                                } else {
919                                    if (upperValue < 1.0e20 && fabs(upperValue - floor(upperValue + 0.5)) > 1.0e-10) {
920                                        // no good
921                                        allGood = false;
922                                        break;
923                                    }
924                                }
925                            }
926                            if (!singletonRow && end > start + 1 && !equality)
927                                allGood = false;
928                            if (!allGood)
929                                break;
930                        }
931                    }
932                }
933            }
934            delete [] count;
935            if (allGood) {
[931]936#if COIN_DEVELOP>1
[1286]937                if (numberObj)
938                    printf("YYYY analysis says all continuous with costs will be integer\n");
[837]939#endif
[1286]940                continuousMultiplier = 1.0;
941            }
942        }
943        delete [] rhs;
[833]944    }
[1286]945    /*
946      Take a first scan to see if there are unfixed continuous variables in the
947      objective.  If so, the minimum objective change could be arbitrarily small.
948      Also pick off the maximum coefficient of an unfixed integer variable.
[2]949
[1286]950      If the objective is found to contain only integer variables, set the
951      fathoming discipline to strict.
952    */
953    double maximumCost = 0.0 ;
954    //double trueIncrement=0.0;
955    int iColumn ;
956    int numberColumns = getNumCols() ;
957    double scaleFactor = 1.0; // due to rhs etc
[1409]958    /*
959      Original model did not have integer bounds.
960    */
[1286]961    if ((specialOptions_&65536) == 0) {
962        /* be on safe side (later look carefully as may be able to
963           to get 0.5 say if bounds are multiples of 0.5 */
964        for (iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
965            if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
966                double value;
967                value = fabs(lower[iColumn]);
968                if (floor(value + 0.5) != value) {
969                    scaleFactor = CoinMin(scaleFactor, 0.5);
970                    if (floor(2.0*value + 0.5) != 2.0*value) {
971                        scaleFactor = CoinMin(scaleFactor, 0.25);
972                        if (floor(4.0*value + 0.5) != 4.0*value) {
973                            scaleFactor = 0.0;
974                        }
975                    }
976                }
977                value = fabs(upper[iColumn]);
978                if (floor(value + 0.5) != value) {
979                    scaleFactor = CoinMin(scaleFactor, 0.5);
980                    if (floor(2.0*value + 0.5) != 2.0*value) {
981                        scaleFactor = CoinMin(scaleFactor, 0.25);
982                        if (floor(4.0*value + 0.5) != 4.0*value) {
983                            scaleFactor = 0.0;
984                        }
985                    }
986                }
987            }
988        }
989    }
990    bool possibleMultiple = continuousMultiplier != 0.0 && scaleFactor != 0.0 ;
991    if (possibleMultiple) {
992        for (iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
993            if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
994                maximumCost = CoinMax(maximumCost, fabs(objective[iColumn])) ;
995            }
996        }
997    }
998    setIntParam(CbcModel::CbcFathomDiscipline, possibleMultiple) ;
999    /*
1000      If a nontrivial increment is possible, try and figure it out. We're looking
1001      for gcd(c<j>) for all c<j> that are coefficients of unfixed integer
1002      variables. Since the c<j> might not be integers, try and inflate them
1003      sufficiently that they look like integers (and we'll deflate the gcd
1004      later).
[2]1005
[1286]1006      2520.0 is used as it is a nice multiple of 2,3,5,7
1007    */
1008    if (possibleMultiple && maximumCost) {
1009        int increment = 0 ;
1010        double multiplier = 2520.0 ;
1011        while (10.0*multiplier*maximumCost < 1.0e8)
1012            multiplier *= 10.0 ;
1013        int bigIntegers = 0; // Count of large costs which are integer
1014        for (iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
1015            if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
1016                double objValue = fabs(objective[iColumn]);
1017                if (!isInteger(iColumn)) {
1018                    if (!coeffMultiplier)
1019                        objValue *= continuousMultiplier;
1020                    else
1021                        objValue *= coeffMultiplier[iColumn];
1022                }
1023                if (objValue) {
1024                    double value = objValue * multiplier ;
1025                    if (value < 2.1e9) {
1026                        int nearest = static_cast<int> (floor(value + 0.5)) ;
1027                        if (fabs(value - floor(value + 0.5)) > 1.0e-8) {
1028                            increment = 0 ;
1029                            break ;
1030                        } else if (!increment) {
1031                            increment = nearest ;
1032                        } else {
1033                            increment = gcd(increment, nearest) ;
1034                        }
1035                    } else {
1036                        // large value - may still be multiple of 1.0
1037                        if (fabs(objValue - floor(objValue + 0.5)) > 1.0e-8) {
1038                            increment = 0;
1039                            break;
1040                        } else {
1041                            bigIntegers++;
1042                        }
1043                    }
1044                }
1045            }
1046        }
1047        delete [] coeffMultiplier;
1048        /*
1049          If the increment beats the current value for objective change, install it.
1050        */
1051        if (increment) {
1052            double value = increment ;
1053            double cutoff = getDblParam(CbcModel::CbcCutoffIncrement) ;
1054            if (bigIntegers) {
1055                // allow for 1.0
1056                increment = gcd(increment, static_cast<int> (multiplier));
1057                value = increment;
1058            }
1059            value /= multiplier ;
1060            value *= scaleFactor;
1061            //trueIncrement=CoinMax(cutoff,value);;
1062            if (value*0.999 > cutoff) {
1063                messageHandler()->message(CBC_INTEGERINCREMENT,
1064                                          messages())
1065                << value << CoinMessageEol ;
1066                setDblParam(CbcModel::CbcCutoffIncrement, value*0.999) ;
1067            }
1068        }
[474]1069    }
[2]1070
[1286]1071    return ;
[2]1072}
1073
[1332]1074/*
1075saveModel called (carved out of) BranchandBound
1076*/
[1330]1077void CbcModel::saveModel(OsiSolverInterface * saveSolver, double * checkCutoffForRestart, bool * feasible)
1078{
1079    if (saveSolver && (specialOptions_&32768) != 0) {
1080        // See if worth trying reduction
1081        *checkCutoffForRestart = getCutoff();
1082        bool tryNewSearch = solverCharacteristics_->reducedCostsAccurate() &&
1083                            (*checkCutoffForRestart < 1.0e20);
1084        int numberColumns = getNumCols();
1085        if (tryNewSearch) {
1086#ifdef CLP_INVESTIGATE
1087            printf("after %d nodes, cutoff %g - looking\n",
1088                   numberNodes_, getCutoff());
1089#endif
1090            saveSolver->resolve();
1091            double direction = saveSolver->getObjSense() ;
1092            double gap = *checkCutoffForRestart - saveSolver->getObjValue() * direction ;
1093            double tolerance;
1094            saveSolver->getDblParam(OsiDualTolerance, tolerance) ;
1095            if (gap <= 0.0)
1096                gap = tolerance;
1097            gap += 100.0 * tolerance;
1098            double integerTolerance = getDblParam(CbcIntegerTolerance) ;
1099
1100            const double *lower = saveSolver->getColLower() ;
1101            const double *upper = saveSolver->getColUpper() ;
1102            const double *solution = saveSolver->getColSolution() ;
1103            const double *reducedCost = saveSolver->getReducedCost() ;
1104
1105            int numberFixed = 0 ;
1106            int numberFixed2 = 0;
1107            for (int i = 0 ; i < numberIntegers_ ; i++) {
1108                int iColumn = integerVariable_[i] ;
1109                double djValue = direction * reducedCost[iColumn] ;
1110                if (upper[iColumn] - lower[iColumn] > integerTolerance) {
1111                    if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) {
1112                        saveSolver->setColUpper(iColumn, lower[iColumn]) ;
1113                        numberFixed++ ;
1114                    } else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) {
1115                        saveSolver->setColLower(iColumn, upper[iColumn]) ;
1116                        numberFixed++ ;
1117                    }
1118                } else {
1119                    numberFixed2++;
1120                }
1121            }
1122#ifdef COIN_DEVELOP
[1409]1123            /*
1124              We're debugging. (specialOptions 1)
1125            */
[1330]1126            if ((specialOptions_&1) != 0) {
1127                const OsiRowCutDebugger *debugger = saveSolver->getRowCutDebugger() ;
1128                if (debugger) {
1129                    printf("Contains optimal\n") ;
1130                    OsiSolverInterface * temp = saveSolver->clone();
1131                    const double * solution = debugger->optimalSolution();
1132                    const double *lower = temp->getColLower() ;
1133                    const double *upper = temp->getColUpper() ;
1134                    int n = temp->getNumCols();
1135                    for (int i = 0; i < n; i++) {
1136                        if (temp->isInteger(i)) {
1137                            double value = floor(solution[i] + 0.5);
1138                            assert (value >= lower[i] && value <= upper[i]);
1139                            temp->setColLower(i, value);
1140                            temp->setColUpper(i, value);
1141                        }
1142                    }
1143                    temp->writeMps("reduced_fix");
1144                    delete temp;
1145                    saveSolver->writeMps("reduced");
1146                } else {
1147                    abort();
1148                }
1149            }
1150            printf("Restart could fix %d integers (%d already fixed)\n",
1151                   numberFixed + numberFixed2, numberFixed2);
1152#endif
1153            numberFixed += numberFixed2;
1154            if (numberFixed*20 < numberColumns)
1155                tryNewSearch = false;
1156        }
1157        if (tryNewSearch) {
1158            // back to solver without cuts?
1159            OsiSolverInterface * solver2 = continuousSolver_->clone();
1160            const double *lower = saveSolver->getColLower() ;
1161            const double *upper = saveSolver->getColUpper() ;
1162            for (int i = 0 ; i < numberIntegers_ ; i++) {
1163                int iColumn = integerVariable_[i] ;
1164                solver2->setColLower(iColumn, lower[iColumn]);
1165                solver2->setColUpper(iColumn, upper[iColumn]);
1166            }
1167            // swap
1168            delete saveSolver;
1169            saveSolver = solver2;
1170            double * newSolution = new double[numberColumns];
1171            double objectiveValue = *checkCutoffForRestart;
1172            CbcSerendipity heuristic(*this);
1173            if (bestSolution_)
1174                heuristic.setInputSolution(bestSolution_, bestObjective_);
1175            heuristic.setFractionSmall(0.9);
1176            heuristic.setFeasibilityPumpOptions(1008013);
1177            // Use numberNodes to say how many are original rows
1178            heuristic.setNumberNodes(continuousSolver_->getNumRows());
1179#ifdef COIN_DEVELOP
1180            if (continuousSolver_->getNumRows() <
1181                    saveSolver->getNumRows())
1182                printf("%d rows added ZZZZZ\n",
1183                       solver_->getNumRows() - continuousSolver_->getNumRows());
1184#endif
1185            int returnCode = heuristic.smallBranchAndBound(saveSolver,
1186                             -1, newSolution,
1187                             objectiveValue,
1188                             *checkCutoffForRestart, "Reduce");
1189            if (returnCode < 0) {
1190#ifdef COIN_DEVELOP
1191                printf("Restart - not small enough to do search after fixing\n");
1192#endif
1193                delete [] newSolution;
1194            } else {
1195                if ((returnCode&1) != 0) {
1196                    // increment number of solutions so other heuristics can test
1197                    numberSolutions_++;
1198                    numberHeuristicSolutions_++;
1199                    lastHeuristic_ = NULL;
1200                    setBestSolution(CBC_ROUNDING, objectiveValue, newSolution) ;
1201                }
1202                delete [] newSolution;
1203                *feasible = false; // stop search
1204            }
[1412]1205#if 0 // probably not needed def CBC_THREAD
1206            if (master_) {
1207                lockThread();
1208                if (parallelMode() > 0) {
1209                    while (master_->waitForThreadsInTree(0)) {
1210                        lockThread();
1211                        double dummyBest;
1212                        tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
1213                        //unlockThread();
1214                    }
1215                }
1216                master_->waitForThreadsInTree(2);
1217                delete master_;
1218                master_ = NULL;
1219                masterThread_ = NULL;
1220            }
1221#endif
[1330]1222        }
1223    }
1224}
1225/*
1226Adds integers, called from BranchandBound()
1227*/
1228void CbcModel::AddIntegers()
1229{
1230    int numberColumns = continuousSolver_->getNumCols();
1231    int numberRows = continuousSolver_->getNumRows();
1232    int * del = new int [CoinMax(numberColumns, numberRows)];
1233    int * original = new int [numberColumns];
1234    char * possibleRow = new char [numberRows];
1235    {
1236        const CoinPackedMatrix * rowCopy = continuousSolver_->getMatrixByRow();
1237        const int * column = rowCopy->getIndices();
1238        const int * rowLength = rowCopy->getVectorLengths();
1239        const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
1240        const double * rowLower = continuousSolver_->getRowLower();
1241        const double * rowUpper = continuousSolver_->getRowUpper();
1242        const double * element = rowCopy->getElements();
1243        for (int i = 0; i < numberRows; i++) {
1244            int nLeft = 0;
1245            bool possible = false;
1246            if (rowLower[i] < -1.0e20) {
1247                double value = rowUpper[i];
1248                if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1249                    possible = true;
1250            } else if (rowUpper[i] > 1.0e20) {
1251                double value = rowLower[i];
1252                if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1253                    possible = true;
1254            } else {
1255                double value = rowUpper[i];
1256                if (rowLower[i] == rowUpper[i] &&
1257                        fabs(value - floor(value + 0.5)) < 1.0e-8)
1258                    possible = true;
1259            }
1260            for (CoinBigIndex j = rowStart[i];
1261                    j < rowStart[i] + rowLength[i]; j++) {
1262                int iColumn = column[j];
1263                if (continuousSolver_->isInteger(iColumn)) {
1264                    if (fabs(element[j]) != 1.0)
1265                        possible = false;
1266                } else {
1267                    nLeft++;
1268                }
1269            }
1270            if (possible || !nLeft)
1271                possibleRow[i] = 1;
1272            else
1273                possibleRow[i] = 0;
1274        }
1275    }
1276    int nDel = 0;
1277    for (int i = 0; i < numberColumns; i++) {
1278        original[i] = i;
1279        if (continuousSolver_->isInteger(i))
1280            del[nDel++] = i;
1281    }
1282    int nExtra = 0;
1283    OsiSolverInterface * copy1 = continuousSolver_->clone();
1284    int nPass = 0;
1285    while (nDel && nPass < 10) {
1286        nPass++;
1287        OsiSolverInterface * copy2 = copy1->clone();
1288        int nLeft = 0;
1289        for (int i = 0; i < nDel; i++)
1290            original[del[i]] = -1;
1291        for (int i = 0; i < numberColumns; i++) {
1292            int kOrig = original[i];
1293            if (kOrig >= 0)
1294                original[nLeft++] = kOrig;
1295        }
1296        assert (nLeft == numberColumns - nDel);
1297        copy2->deleteCols(nDel, del);
1298        numberColumns = copy2->getNumCols();
1299        const CoinPackedMatrix * rowCopy = copy2->getMatrixByRow();
1300        numberRows = rowCopy->getNumRows();
1301        const int * column = rowCopy->getIndices();
1302        const int * rowLength = rowCopy->getVectorLengths();
1303        const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
1304        const double * rowLower = copy2->getRowLower();
1305        const double * rowUpper = copy2->getRowUpper();
1306        const double * element = rowCopy->getElements();
1307        const CoinPackedMatrix * columnCopy = copy2->getMatrixByCol();
1308        const int * columnLength = columnCopy->getVectorLengths();
1309        nDel = 0;
1310        // Could do gcd stuff on ones with costs
1311        for (int i = 0; i < numberRows; i++) {
1312            if (!rowLength[i]) {
1313                del[nDel++] = i;
1314                possibleRow[i] = 1;
1315            } else if (possibleRow[i]) {
1316                if (rowLength[i] == 1) {
1317                    int k = rowStart[i];
1318                    int iColumn = column[k];
1319                    if (!copy2->isInteger(iColumn)) {
1320                        double mult = 1.0 / fabs(element[k]);
1321                        if (rowLower[i] < -1.0e20) {
1322                            double value = rowUpper[i] * mult;
1323                            if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
1324                                del[nDel++] = i;
1325                                if (columnLength[iColumn] == 1) {
1326                                    copy2->setInteger(iColumn);
1327                                    int kOrig = original[iColumn];
1328                                    setOptionalInteger(kOrig);
1329                                }
1330                            }
1331                        } else if (rowUpper[i] > 1.0e20) {
1332                            double value = rowLower[i] * mult;
1333                            if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
1334                                del[nDel++] = i;
1335                                if (columnLength[iColumn] == 1) {
1336                                    copy2->setInteger(iColumn);
1337                                    int kOrig = original[iColumn];
1338                                    setOptionalInteger(kOrig);
1339                                }
1340                            }
1341                        } else {
1342                            double value = rowUpper[i] * mult;
1343                            if (rowLower[i] == rowUpper[i] &&
1344                                    fabs(value - floor(value + 0.5)) < 1.0e-8) {
1345                                del[nDel++] = i;
1346                                copy2->setInteger(iColumn);
1347                                int kOrig = original[iColumn];
1348                                setOptionalInteger(kOrig);
1349                            }
1350                        }
1351                    }
1352                } else {
1353                    // only if all singletons
1354                    bool possible = false;
1355                    if (rowLower[i] < -1.0e20) {
1356                        double value = rowUpper[i];
1357                        if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1358                            possible = true;
1359                    } else if (rowUpper[i] > 1.0e20) {
1360                        double value = rowLower[i];
1361                        if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1362                            possible = true;
1363                    } else {
1364                        double value = rowUpper[i];
1365                        if (rowLower[i] == rowUpper[i] &&
1366                                fabs(value - floor(value + 0.5)) < 1.0e-8)
1367                            possible = true;
1368                    }
1369                    if (possible) {
1370                        for (CoinBigIndex j = rowStart[i];
1371                                j < rowStart[i] + rowLength[i]; j++) {
1372                            int iColumn = column[j];
1373                            if (columnLength[iColumn] != 1 || fabs(element[j]) != 1.0) {
1374                                possible = false;
1375                                break;
1376                            }
1377                        }
1378                        if (possible) {
1379                            for (CoinBigIndex j = rowStart[i];
1380                                    j < rowStart[i] + rowLength[i]; j++) {
1381                                int iColumn = column[j];
1382                                if (!copy2->isInteger(iColumn)) {
1383                                    copy2->setInteger(iColumn);
1384                                    int kOrig = original[iColumn];
1385                                    setOptionalInteger(kOrig);
1386                                }
1387                            }
1388                            del[nDel++] = i;
1389                        }
1390                    }
1391                }
1392            }
1393        }
1394        if (nDel) {
1395            copy2->deleteRows(nDel, del);
1396        }
1397        if (nDel != numberRows) {
1398            nDel = 0;
1399            for (int i = 0; i < numberColumns; i++) {
1400                if (copy2->isInteger(i)) {
1401                    del[nDel++] = i;
1402                    nExtra++;
1403                }
1404            }
1405        } else {
1406            nDel = 0;
1407        }
1408        delete copy1;
1409        copy1 = copy2->clone();
1410        delete copy2;
1411    }
1412    // See if what's left is a network
1413    bool couldBeNetwork = false;
1414    if (copy1->getNumRows() && copy1->getNumCols()) {
1415#ifdef COIN_HAS_CLP
1416        OsiClpSolverInterface * clpSolver
1417        = dynamic_cast<OsiClpSolverInterface *> (copy1);
1418        if (false && clpSolver) {
1419            numberRows = clpSolver->getNumRows();
1420            char * rotate = new char[numberRows];
1421            int n = clpSolver->getModelPtr()->findNetwork(rotate, 1.0);
1422            delete [] rotate;
1423#ifdef CLP_INVESTIGATE
1424            printf("INTA network %d rows out of %d\n", n, numberRows);
1425#endif
1426            if (CoinAbs(n) == numberRows) {
1427                couldBeNetwork = true;
1428                for (int i = 0; i < numberRows; i++) {
1429                    if (!possibleRow[i]) {
1430                        couldBeNetwork = false;
1431#ifdef CLP_INVESTIGATE
1432                        printf("but row %d is bad\n", i);
1433#endif
1434                        break;
1435                    }
1436                }
1437            }
1438        } else
1439#endif
1440        {
1441            numberColumns = copy1->getNumCols();
1442            numberRows = copy1->getNumRows();
1443            const double * rowLower = copy1->getRowLower();
1444            const double * rowUpper = copy1->getRowUpper();
1445            couldBeNetwork = true;
1446            for (int i = 0; i < numberRows; i++) {
1447                if (rowLower[i] > -1.0e20 && fabs(rowLower[i] - floor(rowLower[i] + 0.5)) > 1.0e-12) {
1448                    couldBeNetwork = false;
1449                    break;
1450                }
1451                if (rowUpper[i] < 1.0e20 && fabs(rowUpper[i] - floor(rowUpper[i] + 0.5)) > 1.0e-12) {
1452                    couldBeNetwork = false;
1453                    break;
1454                }
1455            }
1456            if (couldBeNetwork) {
1457                const CoinPackedMatrix  * matrixByCol = copy1->getMatrixByCol();
1458                const double * element = matrixByCol->getElements();
1459                //const int * row = matrixByCol->getIndices();
1460                const CoinBigIndex * columnStart = matrixByCol->getVectorStarts();
1461                const int * columnLength = matrixByCol->getVectorLengths();
1462                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1463                    CoinBigIndex start = columnStart[iColumn];
1464                    CoinBigIndex end = start + columnLength[iColumn];
1465                    if (end > start + 2) {
1466                        couldBeNetwork = false;
1467                        break;
1468                    }
1469                    int type = 0;
1470                    for (CoinBigIndex j = start; j < end; j++) {
1471                        double value = element[j];
1472                        if (fabs(value) != 1.0) {
1473                            couldBeNetwork = false;
1474                            break;
1475                        } else if (value == 1.0) {
1476                            if ((type&1) == 0)
1477                                type |= 1;
1478                            else
1479                                type = 7;
1480                        } else if (value == -1.0) {
1481                            if ((type&2) == 0)
1482                                type |= 2;
1483                            else
1484                                type = 7;
1485                        }
1486                    }
1487                    if (type > 3) {
1488                        couldBeNetwork = false;
1489                        break;
1490                    }
1491                }
1492            }
1493        }
1494    }
1495    if (couldBeNetwork) {
1496        for (int i = 0; i < numberColumns; i++)
1497            setOptionalInteger(original[i]);
1498    }
1499    if (nExtra || couldBeNetwork) {
1500        numberColumns = copy1->getNumCols();
1501        numberRows = copy1->getNumRows();
1502        if (!numberColumns || !numberRows) {
1503            int numberColumns = solver_->getNumCols();
1504            for (int i = 0; i < numberColumns; i++)
1505                assert(solver_->isInteger(i));
1506        }
1507#ifdef CLP_INVESTIGATE
1508        if (couldBeNetwork || nExtra)
1509            printf("INTA %d extra integers, %d left%s\n", nExtra,
1510                   numberColumns,
1511                   couldBeNetwork ? ", all network" : "");
1512#endif
1513        findIntegers(true, 2);
1514        convertToDynamic();
1515    }
1516#ifdef CLP_INVESTIGATE
1517    if (!couldBeNetwork && copy1->getNumCols() &&
1518            copy1->getNumRows()) {
1519        printf("INTA %d rows and %d columns remain\n",
1520               copy1->getNumRows(), copy1->getNumCols());
1521        if (copy1->getNumCols() < 200) {
1522            copy1->writeMps("moreint");
1523            printf("INTA Written remainder to moreint.mps.gz %d rows %d cols\n",
1524                   copy1->getNumRows(), copy1->getNumCols());
1525        }
1526    }
1527#endif
1528    delete copy1;
1529    delete [] del;
1530    delete [] original;
1531    delete [] possibleRow;
1532    // double check increment
1533    analyzeObjective();
1534}
[2]1535/**
1536  \todo
1537  Normally, it looks like we enter here from command dispatch in the main
1538  routine, after calling the solver for an initial solution
1539  (CbcModel::initialSolve, which simply calls the solver's initialSolve
1540  routine.) The first thing we do is call resolve. Presumably there are
1541  circumstances where this is nontrivial? There's also a call from
1542  CbcModel::originalModel (tied up with integer presolve), which should be
1543  checked.
1544
1545*/
1546
1547/*
1548  The overall flow can be divided into three stages:
1549    * Prep: Check that the lp relaxation remains feasible at the root. If so,
1550      do all the setup for B&C.
1551    * Process the root node: Generate cuts, apply heuristics, and in general do
1552      the best we can to resolve the problem without B&C.
1553    * Do B&C search until we hit a limit or exhaust the search tree.
[1286]1554
[2]1555  Keep in mind that in general there is no node in the search tree that
1556  corresponds to the active subproblem. The active subproblem is represented
1557  by the current state of the model,  of the solver, and of the constraint
1558  system held by the solver.
1559*/
[1132]1560#ifdef COIN_HAS_CPX
1561#include "OsiCpxSolverInterface.hpp"
[1271]1562#include "cplex.h"
[1132]1563#endif
[1286]1564void CbcModel::branchAndBound(int doStatistics)
[2]1565
1566{
[1286]1567    /*
[1627]1568      Capture a time stamp before we start (unless set).
[1286]1569    */
[1627]1570    if (!dblParam_[CbcStartSeconds]) {
1571      if (!useElapsedTime())
1572        dblParam_[CbcStartSeconds] = CoinCpuTime();
1573      else
1574        dblParam_[CbcStartSeconds] = CoinGetTimeOfDay();
1575    }
[1286]1576    dblParam_[CbcSmallestChange] = COIN_DBL_MAX;
1577    dblParam_[CbcSumChange] = 0.0;
1578    dblParam_[CbcLargestChange] = 0.0;
1579    intParam_[CbcNumberBranches] = 0;
[1315]1580    // Double check optimization directions line up
1581    dblParam_[CbcOptimizationDirection] = solver_->getObjSense();
[1286]1582    strongInfo_[0] = 0;
1583    strongInfo_[1] = 0;
1584    strongInfo_[2] = 0;
1585    strongInfo_[3] = 0;
1586    strongInfo_[4] = 0;
1587    strongInfo_[5] = 0;
1588    strongInfo_[6] = 0;
1589    numberStrongIterations_ = 0;
1590    currentNode_ = NULL;
1591    // See if should do cuts old way
1592    if (parallelMode() < 0) {
1593        specialOptions_ |= 4096 + 8192;
1594    } else if (parallelMode() > 0) {
1595        specialOptions_ |= 4096;
1596    }
1597    int saveMoreSpecialOptions = moreSpecialOptions_;
1598    if (dynamic_cast<CbcTreeLocal *> (tree_))
1599        specialOptions_ |= 4096 + 8192;
[838]1600#ifdef COIN_HAS_CLP
[1286]1601    {
1602        OsiClpSolverInterface * clpSolver
1603        = dynamic_cast<OsiClpSolverInterface *> (solver_);
1604        if (clpSolver) {
1605            // pass in disaster handler
1606            CbcDisasterHandler handler(this);
1607            clpSolver->passInDisasterHandler(&handler);
1608            // Initialise solvers seed
1609            clpSolver->getModelPtr()->setRandomSeed(1234567);
[1393]1610#ifdef JJF_ZERO
[1286]1611            // reduce factorization frequency
1612            int frequency = clpSolver->getModelPtr()->factorizationFrequency();
1613            clpSolver->getModelPtr()->setFactorizationFrequency(CoinMin(frequency, 120));
[1016]1614#endif
[1286]1615        }
1616    }
[838]1617#endif
[1286]1618    // original solver (only set if pre-processing)
1619    OsiSolverInterface * originalSolver = NULL;
1620    int numberOriginalObjects = numberObjects_;
1621    OsiObject ** originalObject = NULL;
1622    // Save whether there were any objects
1623    bool noObjects = (numberObjects_ == 0);
1624    // Set up strategies
[1409]1625    /*
1626      See if the user has supplied a strategy object and deal with it if present.
1627      The call to setupOther will set numberStrong_ and numberBeforeTrust_, and
1628      perform integer preprocessing, if requested.
[1364]1629
[1409]1630      We need to hang on to a pointer to solver_. setupOther will assign a
1631      preprocessed solver to model, but will instruct assignSolver not to trash the
1632      existing one.
1633    */
[1357]1634    if (strategy_) {
1635        // May do preprocessing
1636        originalSolver = solver_;
1637        strategy_->setupOther(*this);
1638        if (strategy_->preProcessState()) {
1639            // pre-processing done
1640            if (strategy_->preProcessState() < 0) {
1641                // infeasible
1642                handler_->message(CBC_INFEAS, messages_) << CoinMessageEol ;
1643                status_ = 0 ;
1644                secondaryStatus_ = 1;
1645                originalContinuousObjective_ = COIN_DBL_MAX;
1646                return ;
1647            } else if (numberObjects_ && object_) {
1648                numberOriginalObjects = numberObjects_;
1649                // redo sequence
1650                numberIntegers_ = 0;
1651                int numberColumns = getNumCols();
1652                int nOrig = originalSolver->getNumCols();
1653                CglPreProcess * process = strategy_->process();
1654                assert (process);
1655                const int * originalColumns = process->originalColumns();
1656                // allow for cliques etc
1657                nOrig = CoinMax(nOrig, originalColumns[numberColumns-1] + 1);
1658                // try and redo debugger
1659                OsiRowCutDebugger * debugger = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
1660                if (debugger)
1661                    debugger->redoSolution(numberColumns, originalColumns);
[1409]1662                // User-provided solution might have been best. Synchronise.
[1357]1663                if (bestSolution_) {
1664                    // need to redo - in case no better found in BAB
1665                    // just get integer part right
1666                    for (int i = 0; i < numberColumns; i++) {
1667                        int jColumn = originalColumns[i];
1668                        bestSolution_[i] = bestSolution_[jColumn];
1669                    }
1670                }
1671                originalObject = object_;
1672                // object number or -1
1673                int * temp = new int[nOrig];
1674                int iColumn;
1675                for (iColumn = 0; iColumn < nOrig; iColumn++)
1676                    temp[iColumn] = -1;
1677                int iObject;
1678                int nNonInt = 0;
1679                for (iObject = 0; iObject < numberOriginalObjects; iObject++) {
1680                    iColumn = originalObject[iObject]->columnNumber();
1681                    if (iColumn < 0) {
1682                        nNonInt++;
1683                    } else {
1684                        temp[iColumn] = iObject;
1685                    }
1686                }
1687                int numberNewIntegers = 0;
1688                int numberOldIntegers = 0;
1689                int numberOldOther = 0;
1690                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1691                    int jColumn = originalColumns[iColumn];
1692                    if (temp[jColumn] >= 0) {
1693                        int iObject = temp[jColumn];
1694                        CbcSimpleInteger * obj =
1695                            dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ;
1696                        if (obj)
1697                            numberOldIntegers++;
1698                        else
1699                            numberOldOther++;
1700                    } else if (isInteger(iColumn)) {
1701                        numberNewIntegers++;
1702                    }
1703                }
1704                /*
1705                  Allocate an array to hold the indices of the integer variables.
1706                  Make a large enough array for all objects
1707                */
1708                numberObjects_ = numberNewIntegers + numberOldIntegers + numberOldOther + nNonInt;
1709                object_ = new OsiObject * [numberObjects_];
1710                delete [] integerVariable_;
1711                integerVariable_ = new int [numberNewIntegers+numberOldIntegers];
1712                /*
1713                  Walk the variables again, filling in the indices and creating objects for
1714                  the integer variables. Initially, the objects hold the index and upper &
1715                  lower bounds.
1716                */
1717                numberIntegers_ = 0;
1718                int n = originalColumns[numberColumns-1] + 1;
1719                int * backward = new int[n];
1720                int i;
1721                for ( i = 0; i < n; i++)
1722                    backward[i] = -1;
1723                for (i = 0; i < numberColumns; i++)
1724                    backward[originalColumns[i]] = i;
1725                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1726                    int jColumn = originalColumns[iColumn];
1727                    if (temp[jColumn] >= 0) {
1728                        int iObject = temp[jColumn];
1729                        CbcSimpleInteger * obj =
1730                            dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ;
1731                        if (obj) {
1732                            object_[numberIntegers_] = originalObject[iObject]->clone();
1733                            // redo ids etc
1734                            //object_[numberIntegers_]->resetSequenceEtc(numberColumns,originalColumns);
1735                            object_[numberIntegers_]->resetSequenceEtc(numberColumns, backward);
1736                            integerVariable_[numberIntegers_++] = iColumn;
1737                        }
1738                    } else if (isInteger(iColumn)) {
1739                        object_[numberIntegers_] =
1740                            new CbcSimpleInteger(this, iColumn);
1741                        integerVariable_[numberIntegers_++] = iColumn;
1742                    }
1743                }
1744                delete [] backward;
1745                numberObjects_ = numberIntegers_;
1746                // Now append other column stuff
1747                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1748                    int jColumn = originalColumns[iColumn];
1749                    if (temp[jColumn] >= 0) {
1750                        int iObject = temp[jColumn];
1751                        CbcSimpleInteger * obj =
1752                            dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ;
1753                        if (!obj) {
1754                            object_[numberObjects_] = originalObject[iObject]->clone();
1755                            // redo ids etc
1756                            CbcObject * obj =
1757                                dynamic_cast <CbcObject *>(object_[numberObjects_]) ;
1758                            assert (obj);
1759                            obj->redoSequenceEtc(this, numberColumns, originalColumns);
1760                            numberObjects_++;
1761                        }
1762                    }
1763                }
1764                // now append non column stuff
1765                for (iObject = 0; iObject < numberOriginalObjects; iObject++) {
1766                    iColumn = originalObject[iObject]->columnNumber();
1767                    if (iColumn < 0) {
1768                        // already has column numbers changed
1769                        object_[numberObjects_] = originalObject[iObject]->clone();
[1393]1770#ifdef JJF_ZERO
[1357]1771                        // redo ids etc
1772                        CbcObject * obj =
1773                            dynamic_cast <CbcObject *>(object_[numberObjects_]) ;
1774                        assert (obj);
1775                        obj->redoSequenceEtc(this, numberColumns, originalColumns);
1776#endif
1777                        numberObjects_++;
1778                    }
1779                }
1780                delete [] temp;
1781                if (!numberObjects_)
1782                    handler_->message(CBC_NOINT, messages_) << CoinMessageEol ;
1783            } else {
1784                int numberColumns = getNumCols();
1785                CglPreProcess * process = strategy_->process();
1786                assert (process);
1787                const int * originalColumns = process->originalColumns();
1788                // try and redo debugger
1789                OsiRowCutDebugger * debugger = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
1790                if (debugger)
1791                    debugger->redoSolution(numberColumns, originalColumns);
1792            }
1793        } else {
1794            //no preprocessing
1795            originalSolver = NULL;
1796        }
1797        strategy_->setupCutGenerators(*this);
1798        strategy_->setupHeuristics(*this);
1799        // Set strategy print level to models
1800        strategy_->setupPrinting(*this, handler_->logLevel());
1801    }
[1286]1802    eventHappened_ = false;
1803    CbcEventHandler *eventHandler = getEventHandler() ;
1804    if (eventHandler)
1805        eventHandler->setModel(this);
[1016]1806#define CLIQUE_ANALYSIS
[833]1807#ifdef CLIQUE_ANALYSIS
[1286]1808    // set up for probing
[1387]1809    // If we're doing clever stuff with cliques, additional info here.
[1286]1810    if (!parentModel_)
1811        probingInfo_ = new CglTreeProbingInfo(solver_);
1812    else
1813        probingInfo_ = NULL;
[833]1814#else
[1286]1815    probingInfo_ = NULL;
[833]1816#endif
[640]1817
[1286]1818    // Try for dominated columns
1819    if ((specialOptions_&64) != 0) {
1820        CglDuplicateRow dupcuts(solver_);
1821        dupcuts.setMode(2);
1822        CglStored * storedCuts = dupcuts.outDuplicates(solver_);
1823        if (storedCuts) {
1824            printf("adding dup cuts\n");
1825            addCutGenerator(storedCuts, 1, "StoredCuts from dominated",
1826                            true, false, false, -200);
1827        }
[1271]1828    }
[1286]1829    if (!nodeCompare_)
1830        nodeCompare_ = new CbcCompareDefault();;
1831    // See if hot start wanted
1832    CbcCompareBase * saveCompare = NULL;
[1387]1833    // User supplied hotstart. Adapt for preprocessing.
[1286]1834    if (hotstartSolution_) {
1835        if (strategy_ && strategy_->preProcessState() > 0) {
1836            CglPreProcess * process = strategy_->process();
1837            assert (process);
1838            int n = solver_->getNumCols();
1839            const int * originalColumns = process->originalColumns();
1840            // columns should be in order ... but
1841            double * tempS = new double[n];
1842            for (int i = 0; i < n; i++) {
1843                int iColumn = originalColumns[i];
1844                tempS[i] = hotstartSolution_[iColumn];
1845            }
1846            delete [] hotstartSolution_;
1847            hotstartSolution_ = tempS;
1848            if (hotstartPriorities_) {
1849                int * tempP = new int [n];
1850                for (int i = 0; i < n; i++) {
1851                    int iColumn = originalColumns[i];
1852                    tempP[i] = hotstartPriorities_[iColumn];
1853                }
1854                delete [] hotstartPriorities_;
1855                hotstartPriorities_ = tempP;
1856            }
[231]1857        }
[1286]1858        saveCompare = nodeCompare_;
1859        // depth first
1860        nodeCompare_ = new CbcCompareDepth();
[231]1861    }
[1286]1862    if (!problemFeasibility_)
1863        problemFeasibility_ = new CbcFeasibilityBase();
[2]1864# ifdef CBC_DEBUG
[1286]1865    std::string problemName ;
1866    solver_->getStrParam(OsiProbName, problemName) ;
1867    printf("Problem name - %s\n", problemName.c_str()) ;
1868    solver_->setHintParam(OsiDoReducePrint, false, OsiHintDo, 0) ;
[2]1869# endif
[1286]1870    /*
1871      Assume we're done, and see if we're proven wrong.
1872    */
1873    status_ = 0 ;
1874    secondaryStatus_ = 0;
1875    phase_ = 0;
1876    /*
1877      Scan the variables, noting the integer variables. Create an
1878      CbcSimpleInteger object for each integer variable.
1879    */
1880    findIntegers(false) ;
1881    // Say not dynamic pseudo costs
1882    ownership_ &= ~0x40000000;
1883    // If dynamic pseudo costs then do
1884    if (numberBeforeTrust_)
1885        convertToDynamic();
[1387]1886    // Set up char array to say if integer (speed)
[1286]1887    delete [] integerInfo_;
1888    {
1889        int n = solver_->getNumCols();
1890        integerInfo_ = new char [n];
1891        for (int i = 0; i < n; i++) {
1892            if (solver_->isInteger(i))
1893                integerInfo_[i] = 1;
1894            else
1895                integerInfo_[i] = 0;
1896        }
[197]1897    }
[1286]1898    if (preferredWay_) {
1899        // set all unset ones
1900        for (int iObject = 0 ; iObject < numberObjects_ ; iObject++) {
1901            CbcObject * obj =
1902                dynamic_cast <CbcObject *>(object_[iObject]) ;
1903            if (obj && !obj->preferredWay())
1904                obj->setPreferredWay(preferredWay_);
1905        }
[640]1906    }
[1286]1907    /*
1908      Ensure that objects on the lists of OsiObjects, heuristics, and cut
1909      generators attached to this model all refer to this model.
1910    */
1911    synchronizeModel() ;
1912    if (!solverCharacteristics_) {
1913        OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
1914        if (solverCharacteristics) {
1915            solverCharacteristics_ = solverCharacteristics;
1916        } else {
1917            // replace in solver
1918            OsiBabSolver defaultC;
1919            solver_->setAuxiliaryInfo(&defaultC);
1920            solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
1921        }
[395]1922    }
[862]1923
[1286]1924    solverCharacteristics_->setSolver(solver_);
1925    // Set so we can tell we are in initial phase in resolve
1926    continuousObjective_ = -COIN_DBL_MAX ;
1927    /*
1928      Solve the relaxation.
[2]1929
[1286]1930      Apparently there are circumstances where this will be non-trivial --- i.e.,
1931      we've done something since initialSolve that's trashed the solution to the
1932      continuous relaxation.
1933    */
1934    /* Tell solver we are in Branch and Cut
1935       Could use last parameter for subtle differences */
1936    solver_->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
[1053]1937#ifdef COIN_HAS_CLP
[1286]1938    {
1939        OsiClpSolverInterface * clpSolver
1940        = dynamic_cast<OsiClpSolverInterface *> (solver_);
1941        if (clpSolver) {
1942            ClpSimplex * clpSimplex = clpSolver->getModelPtr();
1943            if ((specialOptions_&32) == 0) {
1944                // take off names
1945                clpSimplex->dropNames();
1946            }
1947            // no crunch if mostly continuous
1948            if ((clpSolver->specialOptions()&(1 + 8)) != (1 + 8)) {
1949                int numberColumns = solver_->getNumCols();
1950                if (numberColumns > 1000 && numberIntegers_*4 < numberColumns)
1951                    clpSolver->setSpecialOptions(clpSolver->specialOptions()&(~1));
1952            }
1953            //#define NO_CRUNCH
[1053]1954#ifdef NO_CRUNCH
[1286]1955            printf("TEMP switching off crunch\n");
1956            int iOpt = clpSolver->specialOptions();
1957            iOpt &= ~1;
1958            iOpt |= 65536;
1959            clpSolver->setSpecialOptions(iOpt);
[1053]1960#endif
[1286]1961        }
[1053]1962    }
1963#endif
[1286]1964    bool feasible;
1965    // If NLP then we assume already solved outside branchAndbound
1966    if (!solverCharacteristics_->solverType() || solverCharacteristics_->solverType() == 4) {
1967        feasible = resolve(NULL, 0) != 0 ;
[480]1968    } else {
[1286]1969        // pick up given status
1970        feasible = (solver_->isProvenOptimal() &&
1971                    !solver_->isDualObjectiveLimitReached()) ;
[480]1972    }
[1286]1973    if (problemFeasibility_->feasible(this, 0) < 0) {
1974        feasible = false; // pretend infeasible
[640]1975    }
[1286]1976    numberSavedSolutions_ = 0;
1977    int saveNumberStrong = numberStrong_;
1978    int saveNumberBeforeTrust = numberBeforeTrust_;
1979    /*
1980      If the linear relaxation of the root is infeasible, bail out now. Otherwise,
1981      continue with processing the root node.
1982    */
1983    if (!feasible) {
1984        status_ = 0 ;
1985        if (!solver_->isProvenDualInfeasible()) {
1986            handler_->message(CBC_INFEAS, messages_) << CoinMessageEol ;
1987            secondaryStatus_ = 1;
1988        } else {
1989            handler_->message(CBC_UNBOUNDED, messages_) << CoinMessageEol ;
1990            secondaryStatus_ = 7;
1991        }
1992        originalContinuousObjective_ = COIN_DBL_MAX;
1993        solverCharacteristics_ = NULL;
1994        return ;
1995    } else if (!numberObjects_) {
1996        // nothing to do
1997        solverCharacteristics_ = NULL;
1998        bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
1999        int numberColumns = solver_->getNumCols();
2000        delete [] bestSolution_;
2001        bestSolution_ = new double[numberColumns];
2002        CoinCopyN(solver_->getColSolution(), numberColumns, bestSolution_);
2003        return ;
[640]2004    }
[1409]2005    /*
2006      See if we're using the Osi side of the branching hierarchy. If so, either
2007      convert existing CbcObjects to OsiObjects, or generate them fresh. In the
2008      first case, CbcModel owns the objects on the object_ list. In the second
2009      case, the solver holds the objects and object_ simply points to the
2010      solver's list.
[1364]2011
[1409]2012      080417 The conversion code here (the block protected by `if (obj)') cannot
2013      possibly be correct. On the Osi side, descent is OsiObject -> OsiObject2 ->
2014      all other Osi object classes. On the Cbc side, it's OsiObject -> CbcObject
2015      -> all other Cbc object classes. It's structurally impossible for any Osi
2016      object to descend from CbcObject. The only thing I can see is that this is
2017      really dead code, and object detection is now handled from the Osi side.
2018    */
[1286]2019    // Convert to Osi if wanted
2020    bool useOsiBranching = false;
2021    //OsiBranchingInformation * persistentInfo = NULL;
2022    if (branchingMethod_ && branchingMethod_->chooseMethod()) {
2023        useOsiBranching = true;
2024        //persistentInfo = new OsiBranchingInformation(solver_);
2025        if (numberOriginalObjects) {
2026            for (int iObject = 0 ; iObject < numberObjects_ ; iObject++) {
2027                CbcObject * obj =
2028                    dynamic_cast <CbcObject *>(object_[iObject]) ;
2029                if (obj) {
2030                    CbcSimpleInteger * obj2 =
2031                        dynamic_cast <CbcSimpleInteger *>(obj) ;
2032                    if (obj2) {
2033                        // back to Osi land
2034                        object_[iObject] = obj2->osiObject();
2035                        delete obj;
2036                    } else {
2037                        OsiSimpleInteger * obj3 =
2038                            dynamic_cast <OsiSimpleInteger *>(obj) ;
2039                        if (!obj3) {
2040                            OsiSOS * obj4 =
2041                                dynamic_cast <OsiSOS *>(obj) ;
2042                            if (!obj4) {
2043                                CbcSOS * obj5 =
2044                                    dynamic_cast <CbcSOS *>(obj) ;
2045                                if (obj5) {
2046                                    // back to Osi land
2047                                    object_[iObject] = obj5->osiObject(solver_);
2048                                } else {
2049                                    printf("Code up CbcObject type in Osi land\n");
2050                                    abort();
2051                                }
2052                            }
2053                        }
2054                    }
2055                }
2056            }
2057            // and add to solver
2058            //if (!solver_->numberObjects()) {
2059            solver_->addObjects(numberObjects_, object_);
2060            //} else {
2061            //if (solver_->numberObjects()!=numberOriginalObjects) {
2062            //printf("should have trapped that solver has objects before\n");
2063            //abort();
2064            //}
2065            //}
2066        } else {
[1409]2067            /*
2068              As of 080104, findIntegersAndSOS is misleading --- the default OSI
2069              implementation finds only integers.
2070            */
[1286]2071            // do from solver
2072            deleteObjects(false);
2073            solver_->findIntegersAndSOS(false);
2074            numberObjects_ = solver_->numberObjects();
2075            object_ = solver_->objects();
2076            ownObjects_ = false;
2077        }
2078        branchingMethod_->chooseMethod()->setSolver(solver_);
[640]2079    }
[1387]2080    // take off heuristics if have to (some do not work with SOS, for example)
2081    // object should know what's safe.
[1286]2082    {
2083        int numberOdd = 0;
2084        int numberSOS = 0;
2085        for (int i = 0; i < numberObjects_; i++) {
2086            if (!object_[i]->canDoHeuristics())
2087                numberOdd++;
2088            CbcSOS * obj =
2089                dynamic_cast <CbcSOS *>(object_[i]) ;
2090            if (obj)
2091                numberSOS++;
2092        }
2093        if (numberOdd) {
2094            if (numberHeuristics_) {
2095                int k = 0;
2096                for (int i = 0; i < numberHeuristics_; i++) {
2097                    if (!heuristic_[i]->canDealWithOdd())
2098                        delete heuristic_[i];
2099                    else
2100                        heuristic_[k++] = heuristic_[i];
2101                }
2102                if (!k) {
2103                    delete [] heuristic_;
2104                    heuristic_ = NULL;
2105                }
2106                numberHeuristics_ = k;
2107                handler_->message(CBC_HEURISTICS_OFF, messages_) << numberOdd << CoinMessageEol ;
2108            }
2109        } else if (numberSOS) {
2110            specialOptions_ |= 128; // say can do SOS in dynamic mode
2111            // switch off fast nodes for now
2112            fastNodeDepth_ = -1;
2113        }
2114        if (numberThreads_ > 0) {
2115            // switch off fast nodes for now
2116            fastNodeDepth_ = -1;
2117        }
[931]2118    }
[1286]2119    // Save objective (just so user can access it)
2120    originalContinuousObjective_ = solver_->getObjValue();
2121    bestPossibleObjective_ = originalContinuousObjective_;
2122    sumChangeObjective1_ = 0.0;
2123    sumChangeObjective2_ = 0.0;
2124    /*
2125      OsiRowCutDebugger knows an optimal answer for a subset of MIP problems.
2126      Assuming it recognises the problem, when called upon it will check a cut to
2127      see if it cuts off the optimal answer.
2128    */
2129    // If debugger exists set specialOptions_ bit
2130    if (solver_->getRowCutDebuggerAlways()) {
2131        specialOptions_ |= 1;
2132    }
[56]2133
2134# ifdef CBC_DEBUG
[1286]2135    if ((specialOptions_&1) == 0)
2136        solver_->activateRowCutDebugger(problemName.c_str()) ;
2137    if (solver_->getRowCutDebuggerAlways())
2138        specialOptions_ |= 1;
[2]2139# endif
2140
[1286]2141    /*
2142      Begin setup to process a feasible root node.
2143    */
2144    bestObjective_ = CoinMin(bestObjective_, 1.0e50) ;
2145    if (!bestSolution_) {
2146        numberSolutions_ = 0 ;
2147        numberHeuristicSolutions_ = 0 ;
2148    }
2149    stateOfSearch_ = 0;
2150    // Everything is minimization
2151    {
2152        // needed to sync cutoffs
2153        double value ;
2154        solver_->getDblParam(OsiDualObjectiveLimit, value) ;
2155        dblParam_[CbcCurrentCutoff] = value * solver_->getObjSense();
2156    }
2157    double cutoff = getCutoff() ;
2158    double direction = solver_->getObjSense() ;
2159    dblParam_[CbcOptimizationDirection] = direction;
2160    if (cutoff < 1.0e20 && direction < 0.0)
2161        messageHandler()->message(CBC_CUTOFF_WARNING1,
2162                                  messages())
2163        << cutoff << -cutoff << CoinMessageEol ;
2164    if (cutoff > bestObjective_)
2165        cutoff = bestObjective_ ;
2166    setCutoff(cutoff) ;
2167    /*
2168      We probably already have a current solution, but just in case ...
2169    */
2170    int numberColumns = getNumCols() ;
2171    if (!currentSolution_)
2172        currentSolution_ = new double[numberColumns] ;
2173    testSolution_ = currentSolution_;
2174    /*
2175      Create a copy of the solver, thus capturing the original (root node)
2176      constraint system (aka the continuous system).
2177    */
2178    continuousSolver_ = solver_->clone() ;
[152]2179
[1286]2180    numberRowsAtContinuous_ = getNumRows() ;
2181    solver_->saveBaseModel();
2182    /*
2183      Check the objective to see if we can deduce a nontrivial increment. If
2184      it's better than the current value for CbcCutoffIncrement, it'll be
2185      installed.
2186    */
2187    if (solverCharacteristics_->reducedCostsAccurate())
2188        analyzeObjective() ;
2189    {
2190        // may be able to change cutoff now
2191        double cutoff = getCutoff();
2192        double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
2193        if (cutoff > bestObjective_ - increment) {
2194            cutoff = bestObjective_ - increment ;
2195            setCutoff(cutoff) ;
2196        }
[1271]2197    }
[931]2198#ifdef COIN_HAS_CLP
[1286]2199    // Possible save of pivot method
2200    ClpDualRowPivot * savePivotMethod = NULL;
2201    {
2202        // pass tolerance and increment to solver
2203        OsiClpSolverInterface * clpSolver
2204        = dynamic_cast<OsiClpSolverInterface *> (solver_);
2205        if (clpSolver)
2206            clpSolver->setStuff(getIntegerTolerance(), getCutoffIncrement());
2207    }
[931]2208#endif
[1286]2209    /*
2210      Set up for cut generation. addedCuts_ holds the cuts which are relevant for
2211      the active subproblem. whichGenerator will be used to record the generator
2212      that produced a given cut.
2213    */
2214    maximumWhich_ = 1000 ;
2215    delete [] whichGenerator_;
2216    whichGenerator_ = new int[maximumWhich_] ;
2217    memset(whichGenerator_, 0, maximumWhich_*sizeof(int));
2218    maximumNumberCuts_ = 0 ;
2219    currentNumberCuts_ = 0 ;
2220    delete [] addedCuts_ ;
2221    addedCuts_ = NULL ;
2222    OsiObject ** saveObjects = NULL;
2223    maximumRows_ = numberRowsAtContinuous_;
2224    currentDepth_ = 0;
2225    workingBasis_.resize(maximumRows_, numberColumns);
2226    /*
2227      Set up an empty heap and associated data structures to hold the live set
2228      (problems which require further exploration).
2229    */
2230    CbcCompareDefault * compareActual
2231    = dynamic_cast<CbcCompareDefault *> (nodeCompare_);
2232    if (compareActual) {
2233        compareActual->setBestPossible(direction*solver_->getObjValue());
2234        compareActual->setCutoff(getCutoff());
[1393]2235#ifdef JJF_ZERO
[1357]2236        if (false && !numberThreads_ && !parentModel_) {
[1286]2237            printf("CbcTreeArray ? threads ? parentArray\n");
2238            // Setup new style tree
2239            delete tree_;
2240            tree_ = new CbcTreeArray();
2241        }
[1341]2242#endif
[1132]2243    }
[1286]2244    tree_->setComparison(*nodeCompare_) ;
2245    /*
2246      Used to record the path from a node to the root of the search tree, so that
2247      we can then traverse from the root to the node when restoring a subproblem.
2248    */
2249    maximumDepth_ = 10 ;
2250    delete [] walkback_ ;
2251    walkback_ = new CbcNodeInfo * [maximumDepth_] ;
2252    lastDepth_ = 0;
2253    delete [] lastNodeInfo_ ;
2254    lastNodeInfo_ = new CbcNodeInfo * [maximumDepth_] ;
2255    delete [] lastNumberCuts_ ;
2256    lastNumberCuts_ = new int [maximumDepth_] ;
2257    maximumCuts_ = 100;
2258    lastNumberCuts2_ = 0;
2259    delete [] lastCut_;
2260    lastCut_ = new const OsiRowCut * [maximumCuts_];
2261    /*
2262      Used to generate bound edits for CbcPartialNodeInfo.
2263    */
2264    double * lowerBefore = new double [numberColumns] ;
2265    double * upperBefore = new double [numberColumns] ;
2266    /*
[1409]2267    Set up to run heuristics and generate cuts at the root node. The heavy
2268    lifting is hidden inside the calls to doHeuristicsAtRoot and solveWithCuts.
[2]2269
[1409]2270    To start, tell cut generators they can be a bit more aggressive at the
2271    root node.
[1364]2272
[1409]2273    QUESTION: phase_ = 0 is documented as `initial solve', phase = 1 as `solve
2274        with cuts at root'. Is phase_ = 1 the correct indication when
2275        doHeurisiticsAtRoot is called to run heuristics outside of the main
2276        cut / heurisitc / reoptimise loop in solveWithCuts?
[1364]2277
[1286]2278      Generate cuts at the root node and reoptimise. solveWithCuts does the heavy
2279      lifting. It will iterate a generate/reoptimise loop (including reduced cost
2280      fixing) until no cuts are generated, the change in objective falls off,  or
2281      the limit on the number of rounds of cut generation is exceeded.
[2]2282
[1286]2283      At the end of all this, any cuts will be recorded in cuts and also
2284      installed in the solver's constraint system. We'll have reoptimised, and
2285      removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
2286      adjusted accordingly).
[2]2287
[1286]2288      Tell cut generators they can be a bit more aggressive at root node
2289
2290      TODO: Why don't we make a copy of the solution after solveWithCuts?
2291      TODO: If numberUnsatisfied == 0, don't we have a solution?
2292    */
2293    phase_ = 1;
2294    int iCutGenerator;
2295    for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
[1409]2296        // If parallel switch off global cuts
2297        if (numberThreads_) {
2298            generator_[iCutGenerator]->setGlobalCuts(false);
2299            generator_[iCutGenerator]->setGlobalCutsAtRoot(false);
2300        }
[1286]2301        CglCutGenerator * generator = generator_[iCutGenerator]->generator();
2302        generator->setAggressiveness(generator->getAggressiveness() + 100);
[1271]2303    }
[1286]2304    OsiCuts cuts ;
2305    int anyAction = -1 ;
2306    numberOldActiveCuts_ = 0 ;
2307    numberNewCuts_ = 0 ;
2308    // Array to mark solution
2309    delete [] usedInSolution_;
2310    usedInSolution_ = new int[numberColumns];
2311    CoinZeroN(usedInSolution_, numberColumns);
2312    /*
2313      For printing totals and for CbcNode (numberNodes_)
2314    */
2315    numberIterations_ = 0 ;
2316    numberSolves_ = 0 ;
2317    numberNodes_ = 0 ;
2318    numberNodes2_ = 0 ;
2319    maximumStatistics_ = 0;
2320    maximumDepthActual_ = 0;
2321    numberDJFixed_ = 0.0;
[1315]2322    if (!parentModel_) {
2323        if ((specialOptions_&262144) != 0) {
2324            // create empty stored cuts
2325            //storedRowCuts_ = new CglStored(solver_->getNumCols());
2326        } else if ((specialOptions_&524288) != 0 && storedRowCuts_) {
2327            // tighten and set best solution
2328            // A) tight bounds on integer variables
[1409]2329            /*
2330                storedRowCuts_ are coming in from outside, probably for nonlinear.
2331              John was unsure about origin.
2332            */
[1315]2333            const double * lower = solver_->getColLower();
2334            const double * upper = solver_->getColUpper();
2335            const double * tightLower = storedRowCuts_->tightLower();
2336            const double * tightUpper = storedRowCuts_->tightUpper();
2337            int nTightened = 0;
2338            for (int i = 0; i < numberIntegers_; i++) {
2339                int iColumn = integerVariable_[i];
2340                if (tightLower[iColumn] > lower[iColumn]) {
2341                    nTightened++;
2342                    solver_->setColLower(iColumn, tightLower[iColumn]);
2343                }
2344                if (tightUpper[iColumn] < upper[iColumn]) {
2345                    nTightened++;
2346                    solver_->setColUpper(iColumn, tightUpper[iColumn]);
2347                }
2348            }
2349            if (nTightened)
2350                printf("%d tightened by alternate cuts\n", nTightened);
2351            if (storedRowCuts_->bestObjective() < bestObjective_) {
2352                // B) best solution
2353                double objValue = storedRowCuts_->bestObjective();
2354                setBestSolution(CBC_SOLUTION, objValue,
2355                                storedRowCuts_->bestSolution()) ;
2356                // Do heuristics
2357                // Allow RINS
2358                for (int i = 0; i < numberHeuristics_; i++) {
2359                    CbcHeuristicRINS * rins
2360                    = dynamic_cast<CbcHeuristicRINS *> (heuristic_[i]);
2361                    if (rins) {
2362                        rins->setLastNode(-100);
2363                    }
2364                }
2365            }
2366        }
2367    }
[1409]2368    /*
2369      Run heuristics at the root. This is the only opportunity to run FPump; it
2370      will be removed from the heuristics list by doHeuristicsAtRoot.
2371    */
[1286]2372    // Do heuristics
2373    doHeuristicsAtRoot();
[1409]2374    /*
2375      Grepping through the code, it would appear that this is a command line
2376      debugging hook.  There's no obvious place in the code where this is set to
2377      a negative value.
[1387]2378
[1409]2379      User hook, says John.
2380    */
[1286]2381    if ( intParam_[CbcMaxNumNode] < 0)
2382        eventHappened_ = true; // stop as fast as possible
2383    stoppedOnGap_ = false ;
2384    // See if can stop on gap
2385    bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
2386    double testGap = CoinMax(dblParam_[CbcAllowableGap],
2387                             CoinMax(fabs(bestObjective_), fabs(bestPossibleObjective_))
2388                             * dblParam_[CbcAllowableFractionGap]);
2389    if (bestObjective_ - bestPossibleObjective_ < testGap && getCutoffIncrement() >= 0.0) {
2390        if (bestPossibleObjective_ < getCutoff())
2391            stoppedOnGap_ = true ;
2392        feasible = false;
2393        //eventHappened_=true; // stop as fast as possible
[1271]2394    }
[1409]2395    /*
2396      Set up for statistics collection, if requested. Standard values are
2397      documented in CbcModel.hpp. The magic number 100 will trigger a dump of
2398      CbcSimpleIntegerDynamicPseudoCost objects (no others). Looks like another
2399      command line debugging hook.
2400    */
[1286]2401    statistics_ = NULL;
2402    // Do on switch
2403    if (doStatistics > 0 && doStatistics <= 100) {
2404        maximumStatistics_ = 10000;
2405        statistics_ = new CbcStatistics * [maximumStatistics_];
2406        memset(statistics_, 0, maximumStatistics_*sizeof(CbcStatistics *));
[1271]2407    }
[1286]2408    // See if we can add integers
[1330]2409    if (noObjects && numberIntegers_ < solver_->getNumCols() && (specialOptions_&65536) != 0 && !parentModel_)
[1357]2410        AddIntegers();
[1409]2411    /*
2412      Do an initial round of cut generation for the root node. Depending on the
2413      type of underlying solver, we may want to do this even if the initial query
2414      to the objects indicates they're satisfied.
[1330]2415
[1409]2416      solveWithCuts does the heavy lifting. It will iterate a generate/reoptimise
2417      loop (including reduced cost fixing) until no cuts are generated, the
2418      change in objective falls off,  or the limit on the number of rounds of cut
2419      generation is exceeded.
[1364]2420
[1409]2421      At the end of all this, any cuts will be recorded in cuts and also
2422      installed in the solver's constraint system. We'll have reoptimised, and
2423      removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
2424      adjusted accordingly).
2425    */
[1330]2426    int iObject ;
2427    int preferredWay ;
2428    int numberUnsatisfied = 0 ;
2429    delete [] currentSolution_;
2430    currentSolution_ = new double [numberColumns];
2431    testSolution_ = currentSolution_;
2432    memcpy(currentSolution_, solver_->getColSolution(),
2433           numberColumns*sizeof(double)) ;
2434    // point to useful information
2435    OsiBranchingInformation usefulInfo = usefulInformation();
2436
2437    for (iObject = 0 ; iObject < numberObjects_ ; iObject++) {
2438        double infeasibility =
2439            object_[iObject]->infeasibility(&usefulInfo, preferredWay) ;
2440        if (infeasibility ) numberUnsatisfied++ ;
2441    }
2442    // replace solverType
2443    if (solverCharacteristics_->tryCuts())  {
2444
2445        if (numberUnsatisfied)   {
2446            // User event
2447            if (!eventHappened_ && feasible) {
2448                feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
2449                                         NULL);
2450                if ((specialOptions_&524288) != 0 && !parentModel_
2451                        && storedRowCuts_) {
2452                    if (feasible) {
2453                        /* pick up stuff and try again
2454                        add cuts, maybe keep around
2455                        do best solution and if so new heuristics
2456                        obviously tighten bounds
2457                        */
2458                        // A and B probably done on entry
2459                        // A) tight bounds on integer variables
2460                        const double * lower = solver_->getColLower();
2461                        const double * upper = solver_->getColUpper();
2462                        const double * tightLower = storedRowCuts_->tightLower();
2463                        const double * tightUpper = storedRowCuts_->tightUpper();
2464                        int nTightened = 0;
2465                        for (int i = 0; i < numberIntegers_; i++) {
2466                            int iColumn = integerVariable_[i];
2467                            if (tightLower[iColumn] > lower[iColumn]) {
2468                                nTightened++;
2469                                solver_->setColLower(iColumn, tightLower[iColumn]);
[1286]2470                            }
[1330]2471                            if (tightUpper[iColumn] < upper[iColumn]) {
2472                                nTightened++;
2473                                solver_->setColUpper(iColumn, tightUpper[iColumn]);
2474                            }
[1286]2475                        }
[1330]2476                        if (nTightened)
2477                            printf("%d tightened by alternate cuts\n", nTightened);
2478                        if (storedRowCuts_->bestObjective() < bestObjective_) {
2479                            // B) best solution
2480                            double objValue = storedRowCuts_->bestObjective();
2481                            setBestSolution(CBC_SOLUTION, objValue,
2482                                            storedRowCuts_->bestSolution()) ;
2483                            // Do heuristics
2484                            // Allow RINS
2485                            for (int i = 0; i < numberHeuristics_; i++) {
2486                                CbcHeuristicRINS * rins
2487                                = dynamic_cast<CbcHeuristicRINS *> (heuristic_[i]);
2488                                if (rins) {
2489                                    rins->setLastNode(-100);
[1286]2490                                }
2491                            }
[1330]2492                            doHeuristicsAtRoot();
[1286]2493                        }
[1393]2494#ifdef JJF_ZERO
[1330]2495                        int nCuts = storedRowCuts_->sizeRowCuts();
2496                        // add to global list
2497                        for (int i = 0; i < nCuts; i++) {
2498                            OsiRowCut newCut(*storedRowCuts_->rowCutPointer(i));
2499                            newCut.setGloballyValidAsInteger(2);
2500                            newCut.mutableRow().setTestForDuplicateIndex(false);
2501                            globalCuts_.insert(newCut) ;
[1286]2502                        }
[1330]2503#else
2504                        addCutGenerator(storedRowCuts_, -99, "Stored from previous run",
2505                                        true, false, false, -200);
[1271]2506#endif
[1330]2507                        // Set cuts as active
2508                        delete [] addedCuts_ ;
2509                        maximumNumberCuts_ = cuts.sizeRowCuts();
2510                        if (maximumNumberCuts_) {
2511                            addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
2512                        } else {
2513                            addedCuts_ = NULL;
[1286]2514                        }
[1330]2515                        for (int i = 0; i < maximumNumberCuts_; i++)
2516                            addedCuts_[i] = new CbcCountRowCut(*cuts.rowCutPtr(i),
2517                                                               NULL, -1, -1, 2);
2518                        printf("size %d\n", cuts.sizeRowCuts());
2519                        cuts = OsiCuts();
2520                        currentNumberCuts_ = maximumNumberCuts_;
2521                        feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
2522                                                 NULL);
2523                        for (int i = 0; i < maximumNumberCuts_; i++)
2524                            delete addedCuts_[i];
[1286]2525                    }
[1330]2526                    delete storedRowCuts_;
2527                    storedRowCuts_ = NULL;
[1286]2528                }
[1330]2529            } else {
2530                feasible = false;
[1286]2531            }
[1330]2532        }       else if (solverCharacteristics_->solutionAddsCuts() ||
2533                   solverCharacteristics_->alwaysTryCutsAtRootNode()) {
2534            // may generate cuts and turn the solution
2535            //to an infeasible one
2536            feasible = solveWithCuts(cuts, 1,
2537                                     NULL);
[1286]2538        }
[1271]2539    }
[1330]2540    // check extra info on feasibility
2541    if (!solverCharacteristics_->mipFeasible())
2542        feasible = false;
[2]2543
[1286]2544    // make cut generators less aggressive
2545    for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
2546        CglCutGenerator * generator = generator_[iCutGenerator]->generator();
2547        generator->setAggressiveness(generator->getAggressiveness() - 100);
[1132]2548    }
[1286]2549    currentNumberCuts_ = numberNewCuts_ ;
2550    // See if can stop on gap
2551    bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
2552    testGap = CoinMax(dblParam_[CbcAllowableGap],
2553                      CoinMax(fabs(bestObjective_), fabs(bestPossibleObjective_))
2554                      * dblParam_[CbcAllowableFractionGap]);
2555    if (bestObjective_ - bestPossibleObjective_ < testGap && getCutoffIncrement() >= 0.0) {
2556        if (bestPossibleObjective_ < getCutoff())
2557            stoppedOnGap_ = true ;
2558        feasible = false;
[1132]2559    }
[1286]2560    // User event
2561    if (eventHappened_)
2562        feasible = false;
2563#if defined(COIN_HAS_CLP)&&defined(COIN_HAS_CPX)
[1409]2564    /*
2565      This is the notion of using Cbc stuff to get going, then calling cplex to
2566      finish off.
2567    */
[1286]2568    if (feasible && (specialOptions_&16384) != 0 && fastNodeDepth_ == -2 && !parentModel_) {
2569        // Use Cplex to do search!
2570        double time1 = CoinCpuTime();
2571        OsiClpSolverInterface * clpSolver
2572        = dynamic_cast<OsiClpSolverInterface *> (solver_);
2573        OsiCpxSolverInterface cpxSolver;
2574        double direction = clpSolver->getObjSense();
2575        cpxSolver.setObjSense(direction);
2576        // load up cplex
2577        const CoinPackedMatrix * matrix = continuousSolver_->getMatrixByCol();
2578        const double * rowLower = continuousSolver_->getRowLower();
2579        const double * rowUpper = continuousSolver_->getRowUpper();
2580        const double * columnLower = continuousSolver_->getColLower();
2581        const double * columnUpper = continuousSolver_->getColUpper();
2582        const double * objective = continuousSolver_->getObjCoefficients();
2583        cpxSolver.loadProblem(*matrix, columnLower, columnUpper,
2584                              objective, rowLower, rowUpper);
2585        double * setSol = new double [numberIntegers_];
2586        int * setVar = new int [numberIntegers_];
2587        // cplex doesn't know about objective offset
2588        double offset = clpSolver->getModelPtr()->objectiveOffset();
2589        for (int i = 0; i < numberIntegers_; i++) {
2590            int iColumn = integerVariable_[i];
2591            cpxSolver.setInteger(iColumn);
2592            if (bestSolution_) {
2593                setSol[i] = bestSolution_[iColumn];
2594                setVar[i] = iColumn;
2595            }
2596        }
2597        CPXENVptr env = cpxSolver.getEnvironmentPtr();
2598        CPXLPptr lpPtr = cpxSolver.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL);
2599        cpxSolver.switchToMIP();
2600        if (bestSolution_) {
2601            CPXcopymipstart(env, lpPtr, numberIntegers_, setVar, setSol);
2602        }
2603        if (clpSolver->getNumRows() > continuousSolver_->getNumRows() && false) {
2604            // add cuts
2605            const CoinPackedMatrix * matrix = clpSolver->getMatrixByRow();
2606            const double * rhs = clpSolver->getRightHandSide();
2607            const char * rowSense = clpSolver->getRowSense();
2608            const double * elementByRow = matrix->getElements();
2609            const int * column = matrix->getIndices();
2610            const CoinBigIndex * rowStart = matrix->getVectorStarts();
2611            const int * rowLength = matrix->getVectorLengths();
2612            int nStart = continuousSolver_->getNumRows();
2613            int nRows = clpSolver->getNumRows();
2614            int size = rowStart[nRows-1] + rowLength[nRows-1] -
2615                       rowStart[nStart];
2616            int nAdd = 0;
2617            double * rmatval = new double [size];
2618            int * rmatind = new int [size];
2619            int * rmatbeg = new int [nRows-nStart+1];
2620            size = 0;
2621            rmatbeg[0] = 0;
2622            for (int i = nStart; i < nRows; i++) {
2623                for (int k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2624                    rmatind[size] = column[k];
2625                    rmatval[size++] = elementByRow[k];
2626                }
2627                nAdd++;
2628                rmatbeg[nAdd] = size;
2629            }
2630            CPXaddlazyconstraints(env, lpPtr, nAdd, size,
2631                                  rhs, rowSense, rmatbeg,
2632                                  rmatind, rmatval, NULL);
2633            CPXsetintparam( env, CPX_PARAM_REDUCE,
2634                            // CPX_PREREDUCE_NOPRIMALORDUAL (0)
2635                            CPX_PREREDUCE_PRIMALONLY);
2636        }
2637        if (getCutoff() < 1.0e50) {
2638            double useCutoff = getCutoff() + offset;
2639            if (bestObjective_ < 1.0e50)
2640                useCutoff = bestObjective_ + offset + 1.0e-7;
2641            cpxSolver.setDblParam(OsiDualObjectiveLimit, useCutoff*
2642                                  direction);
2643            if ( direction > 0.0 )
2644                CPXsetdblparam( env, CPX_PARAM_CUTUP, useCutoff ) ; // min
2645            else
2646                CPXsetdblparam( env, CPX_PARAM_CUTLO, useCutoff ) ; // max
2647        }
2648        CPXsetdblparam(env, CPX_PARAM_EPGAP, dblParam_[CbcAllowableFractionGap]);
2649        delete [] setSol;
2650        delete [] setVar;
2651        char printBuffer[200];
2652        if (offset) {
2653            sprintf(printBuffer, "Add %g to all Cplex messages for true objective",
2654                    -offset);
2655            messageHandler()->message(CBC_GENERAL, messages())
2656            << printBuffer << CoinMessageEol ;
2657            cpxSolver.setDblParam(OsiObjOffset, offset);
2658        }
2659        cpxSolver.branchAndBound();
2660        double timeTaken = CoinCpuTime() - time1;
2661        sprintf(printBuffer, "Cplex took %g seconds",
2662                timeTaken);
2663        messageHandler()->message(CBC_GENERAL, messages())
2664        << printBuffer << CoinMessageEol ;
2665        numberExtraNodes_ = CPXgetnodecnt(env, lpPtr);
2666        numberExtraIterations_ = CPXgetmipitcnt(env, lpPtr);
2667        double value = cpxSolver.getObjValue() * direction;
2668        if (cpxSolver.isProvenOptimal() && value <= getCutoff()) {
2669            feasible = true;
2670            clpSolver->setWarmStart(NULL);
2671            // try and do solution
2672            double * newSolution =
2673                CoinCopyOfArray(cpxSolver.getColSolution(),
2674                                getNumCols());
2675            setBestSolution(CBC_STRONGSOL, value, newSolution) ;
2676            delete [] newSolution;
2677        }
2678        feasible = false;
[1132]2679    }
[1286]2680#endif
[1409]2681    /*
2682      A hook to use clp to quickly explore some part of the tree.
2683    */
[1286]2684    if (fastNodeDepth_ == 1000 &&/*!parentModel_*/(specialOptions_&2048) == 0) {
2685        fastNodeDepth_ = -1;
2686        CbcObject * obj =
2687            new CbcFollowOn(this);
2688        obj->setPriority(1);
2689        addObjects(1, &obj);
[1132]2690    }
[1286]2691    int saveNumberSolves = numberSolves_;
2692    int saveNumberIterations = numberIterations_;
2693    if (fastNodeDepth_ >= 0 &&/*!parentModel_*/(specialOptions_&2048) == 0) {
2694        // add in a general depth object doClp
2695        int type = (fastNodeDepth_ <= 100) ? fastNodeDepth_ : -(fastNodeDepth_ - 100);
2696        CbcObject * obj =
2697            new CbcGeneralDepth(this, type);
2698        addObjects(1, &obj);
2699        // mark as done
2700        fastNodeDepth_ += 1000000;
2701        delete obj;
2702        // fake number of objects
2703        numberObjects_--;
2704        if (parallelMode() < -1) {
2705            // But make sure position is correct
2706            OsiObject * obj2 = object_[numberObjects_];
2707            obj = dynamic_cast<CbcObject *> (obj2);
2708            assert (obj);
2709            obj->setPosition(numberObjects_);
2710        }
[1132]2711    }
[1286]2712#ifdef COIN_HAS_CLP
2713#ifdef NO_CRUNCH
2714    if (true) {
2715        OsiClpSolverInterface * clpSolver
2716        = dynamic_cast<OsiClpSolverInterface *> (solver_);
2717        if (clpSolver && !parentModel_) {
2718            ClpSimplex * clpSimplex = clpSolver->getModelPtr();
2719            clpSimplex->setSpecialOptions(clpSimplex->specialOptions() | 131072);
2720            //clpSimplex->startPermanentArrays();
2721            clpSimplex->setPersistenceFlag(2);
2722        }
[1121]2723    }
[931]2724#endif
2725#endif
[1286]2726// Save copy of solver
2727    OsiSolverInterface * saveSolver = NULL;
2728    if (!parentModel_ && (specialOptions_&(512 + 32768)) != 0)
2729        saveSolver = solver_->clone();
2730    double checkCutoffForRestart = 1.0e100;
[1330]2731    saveModel(saveSolver, &checkCutoffForRestart, &feasible);
[1315]2732    if ((specialOptions_&262144) != 0 && !parentModel_) {
2733        // Save stuff and return!
2734        storedRowCuts_->saveStuff(bestObjective_, bestSolution_,
2735                                  solver_->getColLower(),
2736                                  solver_->getColUpper());
2737        delete [] lowerBefore;
2738        delete [] upperBefore;
2739        delete saveSolver;
2740        return;
2741    }
[1286]2742    /*
2743      We've taken the continuous relaxation as far as we can. Time to branch.
2744      The first order of business is to actually create a node. chooseBranch
2745      currently uses strong branching to evaluate branch object candidates,
2746      unless forced back to simple branching. If chooseBranch concludes that a
2747      branching candidate is monotone (anyAction == -1) or infeasible (anyAction
2748      == -2) when forced to integer values, it returns here immediately.
[2]2749
[1286]2750      Monotone variables trigger a call to resolve(). If the problem remains
2751      feasible, try again to choose a branching variable. At the end of the loop,
2752      resolved == true indicates that some variables were fixed.
[2]2753
[1286]2754      Loss of feasibility will result in the deletion of newNode.
2755    */
[2]2756
[1286]2757    bool resolved = false ;
2758    CbcNode *newNode = NULL ;
2759    numberFixedAtRoot_ = 0;
2760    numberFixedNow_ = 0;
2761    int numberIterationsAtContinuous = numberIterations_;
2762    //solverCharacteristics_->setSolver(solver_);
2763    if (feasible) {
2764        //#define HOTSTART -1
[1271]2765#if HOTSTART<0
[1286]2766        if (bestSolution_ && !parentModel_ && !hotstartSolution_ &&
2767                (moreSpecialOptions_&1024) != 0) {
2768            // Set priorities so only branch on ones we need to
2769            // use djs and see if only few branches needed
[1271]2770#ifndef NDEBUG
[1286]2771            double integerTolerance = getIntegerTolerance() ;
[1271]2772#endif
[1286]2773            bool possible = true;
2774            const double * saveLower = continuousSolver_->getColLower();
2775            const double * saveUpper = continuousSolver_->getColUpper();
2776            for (int i = 0; i < numberObjects_; i++) {
2777                const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object_[i]);
2778                if (thisOne) {
2779                    int iColumn = thisOne->columnNumber();
2780                    if (saveUpper[iColumn] > saveLower[iColumn] + 1.5) {
2781                        possible = false;
2782                        break;
2783                    }
2784                } else {
2785                    possible = false;
2786                    break;
2787                }
2788            }
2789            if (possible) {
2790                OsiSolverInterface * solver = continuousSolver_->clone();
2791                int numberColumns = solver->getNumCols();
2792                for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
2793                    double value = bestSolution_[iColumn] ;
2794                    value = CoinMax(value, saveLower[iColumn]) ;
2795                    value = CoinMin(value, saveUpper[iColumn]) ;
2796                    value = floor(value + 0.5);
2797                    if (solver->isInteger(iColumn)) {
2798                        solver->setColLower(iColumn, value);
2799                        solver->setColUpper(iColumn, value);
2800                    }
2801                }
2802                solver->setHintParam(OsiDoDualInResolve, false, OsiHintTry);
2803                // objlim and all slack
2804                double direction = solver->getObjSense();
2805                solver->setDblParam(OsiDualObjectiveLimit, 1.0e50*direction);
2806                CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis *> (solver->getEmptyWarmStart());
2807                solver->setWarmStart(basis);
2808                delete basis;
2809                bool changed = true;
2810                hotstartPriorities_ = new int [numberColumns];
2811                for (int iColumn = 0; iColumn < numberColumns; iColumn++)
2812                    hotstartPriorities_[iColumn] = 1;
2813                while (changed) {
2814                    changed = false;
2815                    solver->resolve();
2816                    if (!solver->isProvenOptimal()) {
2817                        possible = false;
2818                        break;
2819                    }
2820                    const double * dj = solver->getReducedCost();
2821                    const double * colLower = solver->getColLower();
2822                    const double * colUpper = solver->getColUpper();
2823                    const double * solution = solver->getColSolution();
2824                    int nAtLbNatural = 0;
2825                    int nAtUbNatural = 0;
2826                    int nZeroDj = 0;
2827                    int nForced = 0;
2828                    for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
2829                        double value = solution[iColumn] ;
2830                        value = CoinMax(value, saveLower[iColumn]) ;
2831                        value = CoinMin(value, saveUpper[iColumn]) ;
2832                        if (solver->isInteger(iColumn)) {
2833                            assert(fabs(value - solution[iColumn]) <= integerTolerance) ;
2834                            if (hotstartPriorities_[iColumn] == 1) {
2835                                if (dj[iColumn] < -1.0e-6) {
2836                                    // negative dj
2837                                    if (saveUpper[iColumn] == colUpper[iColumn]) {
2838                                        nAtUbNatural++;
2839                                        hotstartPriorities_[iColumn] = 2;
2840                                        solver->setColLower(iColumn, saveLower[iColumn]);
2841                                        solver->setColUpper(iColumn, saveUpper[iColumn]);
2842                                    } else {
2843                                        nForced++;
2844                                    }
2845                                } else if (dj[iColumn] > 1.0e-6) {
2846                                    // positive dj
2847                                    if (saveLower[iColumn] == colLower[iColumn]) {
2848                                        nAtLbNatural++;
2849                                        hotstartPriorities_[iColumn] = 2;
2850                                        solver->setColLower(iColumn, saveLower[iColumn]);
2851                                        solver->setColUpper(iColumn, saveUpper[iColumn]);
2852                                    } else {
2853                                        nForced++;
2854                                    }
2855                                } else {
2856                                    // zero dj
2857                                    nZeroDj++;
2858                                }
2859                            }
2860                        }
2861                    }
[1271]2862#ifdef CLP_INVESTIGATE
[1286]2863                    printf("%d forced, %d naturally at lower, %d at upper - %d zero dj\n",
2864                           nForced, nAtLbNatural, nAtUbNatural, nZeroDj);
[1271]2865#endif
[1286]2866                    if (nAtLbNatural || nAtUbNatural) {
2867                        changed = true;
2868                    } else {
2869                        if (nForced + nZeroDj > 50 ||
2870                                (nForced + nZeroDj)*4 > numberIntegers_)
2871                            possible = false;
2872                    }
2873                }
2874                delete solver;
2875            }
2876            if (possible) {
2877                setHotstartSolution(bestSolution_);
2878                if (!saveCompare) {
2879                    // create depth first comparison
2880                    saveCompare = nodeCompare_;
2881                    // depth first
2882                    nodeCompare_ = new CbcCompareDepth();
2883                    tree_->setComparison(*nodeCompare_) ;
2884                }
2885            } else {
2886                delete [] hotstartPriorities_;
2887                hotstartPriorities_ = NULL;
2888            }
2889        }
[1271]2890#endif
2891#if HOTSTART>0
[1286]2892        if (hotstartSolution_ && !hotstartPriorities_) {
2893            // Set up hot start
2894            OsiSolverInterface * solver = solver_->clone();
2895            double direction = solver_->getObjSense() ;
2896            int numberColumns = solver->getNumCols();
2897            double * saveLower = CoinCopyOfArray(solver->getColLower(), numberColumns);
2898            double * saveUpper = CoinCopyOfArray(solver->getColUpper(), numberColumns);
2899            // move solution
2900            solver->setColSolution(hotstartSolution_);
2901            // point to useful information
2902            const double * saveSolution = testSolution_;
2903            testSolution_ = solver->getColSolution();
2904            OsiBranchingInformation usefulInfo = usefulInformation();
2905            testSolution_ = saveSolution;
2906            /*
2907            Run through the objects and use feasibleRegion() to set variable bounds
2908            so as to fix the variables specified in the objects at their value in this
2909            solution. Since the object list contains (at least) one object for every
2910            integer variable, this has the effect of fixing all integer variables.
2911            */
2912            for (int i = 0; i < numberObjects_; i++)
2913                object_[i]->feasibleRegion(solver, &usefulInfo);
2914            solver->resolve();
2915            assert (solver->isProvenOptimal());
2916            double gap = CoinMax((solver->getObjValue() - solver_->getObjValue()) * direction, 0.0) ;
2917            const double * dj = solver->getReducedCost();
2918            const double * colLower = solver->getColLower();
2919            const double * colUpper = solver->getColUpper();
2920            const double * solution = solver->getColSolution();
2921            int nAtLbNatural = 0;
2922            int nAtUbNatural = 0;
2923            int nAtLbNaturalZero = 0;
2924            int nAtUbNaturalZero = 0;
2925            int nAtLbFixed = 0;
2926            int nAtUbFixed = 0;
2927            int nAtOther = 0;
2928            int nAtOtherNatural = 0;
2929            int nNotNeeded = 0;
2930            delete [] hotstartSolution_;
2931            hotstartSolution_ = new double [numberColumns];
2932            delete [] hotstartPriorities_;
2933            hotstartPriorities_ = new int [numberColumns];
2934            int * order = (int *) saveUpper;
2935            int nFix = 0;
2936            double bestRatio = COIN_DBL_MAX;
2937            for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
2938                double value = solution[iColumn] ;
2939                value = CoinMax(value, saveLower[iColumn]) ;
2940                value = CoinMin(value, saveUpper[iColumn]) ;
2941                double sortValue = COIN_DBL_MAX;
2942                if (solver->isInteger(iColumn)) {
2943                    assert(fabs(value - solution[iColumn]) <= 1.0e-5) ;
2944                    double value2 = floor(value + 0.5);
2945                    if (dj[iColumn] < -1.0e-6) {
2946                        // negative dj
2947                        //assert (value2==colUpper[iColumn]);
2948                        if (saveUpper[iColumn] == colUpper[iColumn]) {
2949                            nAtUbNatural++;
2950                            sortValue = 0.0;
2951                            double value = -dj[iColumn];
2952                            if (value > gap)
2953                                nFix++;
2954                            else if (gap < value*bestRatio)
2955                                bestRatio = gap / value;
2956                            if (saveLower[iColumn] != colLower[iColumn]) {
2957                                nNotNeeded++;
2958                                sortValue = 1.0e20;
2959                            }
2960                        } else if (saveLower[iColumn] == colUpper[iColumn]) {
2961                            nAtLbFixed++;
2962                            sortValue = dj[iColumn];
2963                        } else {
2964                            nAtOther++;
2965                            sortValue = 0.0;
2966                            if (saveLower[iColumn] != colLower[iColumn] &&
2967                                    saveUpper[iColumn] != colUpper[iColumn]) {
2968                                nNotNeeded++;
2969                                sortValue = 1.0e20;
2970                            }
2971                        }
2972                    } else if (dj[iColumn] > 1.0e-6) {
2973                        // positive dj
2974                        //assert (value2==colLower[iColumn]);
2975                        if (saveLower[iColumn] == colLower[iColumn]) {
2976                            nAtLbNatural++;
2977                            sortValue = 0.0;
2978                            double value = dj[iColumn];
2979                            if (value > gap)
2980                                nFix++;
2981                            else if (gap < value*bestRatio)
2982                                bestRatio = gap / value;
2983                            if (saveUpper[iColumn] != colUpper[iColumn]) {
2984                                nNotNeeded++;
2985                                sortValue = 1.0e20;
2986                            }
2987                        } else if (saveUpper[iColumn] == colLower[iColumn]) {
2988                            nAtUbFixed++;
2989                            sortValue = -dj[iColumn];
2990                        } else {
2991                            nAtOther++;
2992                            sortValue = 0.0;
2993                            if (saveLower[iColumn] != colLower[iColumn] &&
2994                                    saveUpper[iColumn] != colUpper[iColumn]) {
2995                                nNotNeeded++;
2996                                sortValue = 1.0e20;
2997                            }
2998                        }
2999                    } else {
3000                        // zero dj
3001                        if (value2 == saveUpper[iColumn]) {
3002                            nAtUbNaturalZero++;
3003                            sortValue = 0.0;
3004                            if (saveLower[iColumn] != colLower[iColumn]) {
3005                                nNotNeeded++;
3006                                sortValue = 1.0e20;
3007                            }
3008                        } else if (value2 == saveLower[iColumn]) {
3009                            nAtLbNaturalZero++;
3010                            sortValue = 0.0;
3011                        } else {
3012                            nAtOtherNatural++;
3013                            sortValue = 0.0;
3014                            if (saveLower[iColumn] != colLower[iColumn] &&
3015                                    saveUpper[iColumn] != colUpper[iColumn]) {
3016                                nNotNeeded++;
3017                                sortValue = 1.0e20;
3018                            }
3019                        }
3020                    }
[1016]3021#if HOTSTART==3
[1286]3022                    sortValue = -fabs(dj[iColumn]);
[1016]3023#endif
[1286]3024                }
3025                hotstartSolution_[iColumn] = value ;
3026                saveLower[iColumn] = sortValue;
3027                order[iColumn] = iColumn;
3028            }
3029            printf("** can fix %d columns - best ratio for others is %g on gap of %g\n",
3030                   nFix, bestRatio, gap);
3031            int nNeg = 0;
3032            CoinSort_2(saveLower, saveLower + numberColumns, order);
3033            for (int i = 0; i < numberColumns; i++) {
3034                if (saveLower[i] < 0.0) {
3035                    nNeg++;
[1016]3036#if HOTSTART==2||HOTSTART==3
[1286]3037                    // swap sign ?
3038                    saveLower[i] = -saveLower[i];
[1016]3039#endif
[1286]3040                }
3041            }
3042            CoinSort_2(saveLower, saveLower + nNeg, order);
3043            for (int i = 0; i < numberColumns; i++) {
[1016]3044#if HOTSTART==1
[1286]3045                hotstartPriorities_[order[i]] = 100;
[1016]3046#else
[1286]3047                hotstartPriorities_[order[i]] = -(i + 1);
[1016]3048#endif
[1286]3049            }
3050            printf("nAtLbNat %d,nAtUbNat %d,nAtLbNatZero %d,nAtUbNatZero %d,nAtLbFixed %d,nAtUbFixed %d,nAtOther %d,nAtOtherNat %d, useless %d %d\n",
3051                   nAtLbNatural,
3052                   nAtUbNatural,
3053                   nAtLbNaturalZero,
3054                   nAtUbNaturalZero,
3055                   nAtLbFixed,
3056                   nAtUbFixed,
3057                   nAtOther,
3058                   nAtOtherNatural, nNotNeeded, nNeg);
3059            delete [] saveLower;
3060            delete [] saveUpper;
3061            if (!saveCompare) {
3062                // create depth first comparison
3063                saveCompare = nodeCompare_;
3064                // depth first
3065                nodeCompare_ = new CbcCompareDepth();
3066                tree_->setComparison(*nodeCompare_) ;
3067            }
3068        }
[1016]3069#endif
[1286]3070        newNode = new CbcNode ;
3071        // Set objective value (not so obvious if NLP etc)
3072        setObjectiveValue(newNode, NULL);
3073        anyAction = -1 ;
3074        // To make depth available we may need a fake node
3075        CbcNode fakeNode;
3076        if (!currentNode_) {
3077            // Not true if sub trees assert (!numberNodes_);
3078            currentNode_ = &fakeNode;
3079        }
3080        phase_ = 3;
3081        // only allow 1000 passes
3082        int numberPassesLeft = 1000;
3083        // This is first crude step
3084        if (numberAnalyzeIterations_) {
3085            delete [] analyzeResults_;
3086            analyzeResults_ = new double [4*numberIntegers_];
3087            numberFixedAtRoot_ = newNode->analyze(this, analyzeResults_);
3088            if (numberFixedAtRoot_ > 0) {
3089                printf("%d fixed by analysis\n", numberFixedAtRoot_);
3090                setPointers(solver_);
3091                numberFixedNow_ = numberFixedAtRoot_;
3092            } else if (numberFixedAtRoot_ < 0) {
3093                printf("analysis found to be infeasible\n");
3094                anyAction = -2;
3095                delete newNode ;
3096                newNode = NULL ;
3097                feasible = false ;
3098            }
3099        }
3100        OsiSolverBranch * branches = NULL;
3101        anyAction = chooseBranch(newNode, numberPassesLeft, NULL, cuts, resolved,
3102                                 NULL, NULL, NULL, branches);
3103        if (anyAction == -2 || newNode->objectiveValue() >= cutoff) {
3104            if (anyAction != -2) {
3105                // zap parent nodeInfo
[640]3106#ifdef COIN_DEVELOP
[1286]3107                printf("zapping CbcNodeInfo %x\n", reinterpret_cast<int>(newNode->nodeInfo()->parent()));
[640]3108#endif
[1286]3109                if (newNode->nodeInfo())
3110                    newNode->nodeInfo()->nullParent();
3111            }
3112            delete newNode ;
3113            newNode = NULL ;
3114            feasible = false ;
3115        }
[640]3116    }
[1315]3117    if (newNode && probingInfo_) {
3118        int number01 = probingInfo_->numberIntegers();
3119        //const fixEntry * entry = probingInfo_->fixEntries();
3120        const int * toZero = probingInfo_->toZero();
3121        //const int * toOne = probingInfo_->toOne();
3122        //const int * integerVariable = probingInfo_->integerVariable();
3123        if (toZero[number01]) {
3124            CglTreeProbingInfo info(*probingInfo_);
[1393]3125#ifdef JJF_ZERO
[1409]3126            /*
3127              Marginal idea. Further exploration probably good. Build some extra
3128              cliques from probing info. Not quite worth the effort?
3129            */
[1315]3130            OsiSolverInterface * fake = info.analyze(*solver_, 1);
3131            if (fake) {
3132                fake->writeMps("fake");
3133                CglFakeClique cliqueGen(fake);
3134                cliqueGen.setStarCliqueReport(false);
3135                cliqueGen.setRowCliqueReport(false);
3136                cliqueGen.setMinViolation(0.1);
3137                addCutGenerator(&cliqueGen, 1, "Fake cliques", true, false, false, -200);
3138                generator_[numberCutGenerators_-1]->setTiming(true);
3139            }
3140#endif
3141            if (probingInfo_->packDown()) {
3142#ifdef CLP_INVESTIGATE
3143                printf("%d implications on %d 0-1\n", toZero[number01], number01);
3144#endif
[1409]3145                // Create a cut generator that remembers implications discovered at root.
[1315]3146                CglImplication implication(probingInfo_);
3147                addCutGenerator(&implication, 1, "ImplicationCuts", true, false, false, -200);
3148            } else {
3149                delete probingInfo_;
3150                probingInfo_ = NULL;
3151            }
3152        } else {
3153            delete probingInfo_;
3154
3155            probingInfo_ = NULL;
3156        }
3157    }
[1286]3158    /*
3159      At this point, the root subproblem is infeasible or fathomed by bound
3160      (newNode == NULL), or we're live with an objective value that satisfies the
3161      current objective cutoff.
3162    */
3163    assert (!newNode || newNode->objectiveValue() <= cutoff) ;
3164    // Save address of root node as we don't want to delete it
3165    // initialize for print out
3166    int lastDepth = 0;
3167    int lastUnsatisfied = 0;
3168    if (newNode)
3169        lastUnsatisfied = newNode->numberUnsatisfied();
3170    /*
3171      The common case is that the lp relaxation is feasible but doesn't satisfy
3172      integrality (i.e., newNode->branchingObject(), indicating we've been able to
3173      select a branching variable). Remove any cuts that have gone slack due to
3174      forcing monotone variables. Then tack on an CbcFullNodeInfo object and full
3175      basis (via createInfo()) and stash the new cuts in the nodeInfo (via
3176      addCuts()). If, by some miracle, we have an integral solution at the root
3177      (newNode->branchingObject() is NULL), takeOffCuts() will ensure that the solver holds
3178      a valid solution for use by setBestSolution().
3179    */
3180    CoinWarmStartBasis *lastws = NULL ;
3181    if (feasible && newNode->branchingObject()) {
3182        if (resolved) {
3183            takeOffCuts(cuts, false, NULL) ;
[2]3184#     ifdef CHECK_CUT_COUNTS
[1286]3185            { printf("Number of rows after chooseBranch fix (root)"
3186                         "(active only) %d\n",
3187                         numberRowsAtContinuous_ + numberNewCuts_ + numberOldActiveCuts_) ;
3188                const CoinWarmStartBasis* debugws =
3189                    dynamic_cast <const CoinWarmStartBasis*>(solver_->getWarmStart()) ;
3190                debugws->print() ;
3191                delete debugws ;
3192            }
[2]3193#     endif
[1286]3194        }
3195        //newNode->createInfo(this,NULL,NULL,NULL,NULL,0,0) ;
3196        //newNode->nodeInfo()->addCuts(cuts,
3197        //                       newNode->numberBranches(),whichGenerator_) ;
3198        if (lastws) delete lastws ;
3199        lastws = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
[2]3200    }
[1286]3201    /*
3202      Continuous data to be used later
3203    */
3204    continuousObjective_ = solver_->getObjValue() * solver_->getObjSense();
3205    continuousInfeasibilities_ = 0 ;
3206    if (newNode) {
3207        continuousObjective_ = newNode->objectiveValue() ;
3208        delete [] continuousSolution_;
3209        continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
3210                                              numberColumns);
3211        continuousInfeasibilities_ = newNode->numberUnsatisfied() ;
3212    }
3213    /*
3214      Bound may have changed so reset in objects
3215    */
3216    {
3217        int i ;
3218        for (i = 0; i < numberObjects_; i++)
3219            object_[i]->resetBounds(solver_) ;
3220    }
3221    /*
3222      Feasible? Then we should have either a live node prepped for future
3223      expansion (indicated by variable() >= 0), or (miracle of miracles) an
3224      integral solution at the root node.
[2]3225
[1286]3226      initializeInfo sets the reference counts in the nodeInfo object.  Since
3227      this node is still live, push it onto the heap that holds the live set.
3228    */
3229    double bestValue = 0.0 ;
3230    if (newNode) {
3231        bestValue = newNode->objectiveValue();
3232        if (newNode->branchingObject()) {
3233            newNode->initializeInfo() ;
3234            tree_->push(newNode) ;
3235            if (statistics_) {
3236                if (numberNodes2_ == maximumStatistics_) {
3237                    maximumStatistics_ = 2 * maximumStatistics_;
3238                    CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_];
3239                    memset(temp, 0, maximumStatistics_*sizeof(CbcStatistics *));
3240                    memcpy(temp, statistics_, numberNodes2_*sizeof(CbcStatistics *));
3241                    delete [] statistics_;
3242                    statistics_ = temp;
3243                }
3244                assert (!statistics_[numberNodes2_]);
3245                statistics_[numberNodes2_] = new CbcStatistics(newNode, this);
3246            }
3247            numberNodes2_++;
[2]3248#     ifdef CHECK_NODE
[1286]3249            printf("Node %x on tree\n", newNode) ;
[2]3250#     endif
[1286]3251        } else {
3252            // continuous is integer
3253            double objectiveValue = newNode->objectiveValue();
3254            setBestSolution(CBC_SOLUTION, objectiveValue,
3255                            solver_->getColSolution()) ;
3256            delete newNode ;
3257            newNode = NULL ;
3258        }
[2]3259    }
3260
[1286]3261    if (printFrequency_ <= 0) {
3262        printFrequency_ = 1000 ;
3263        if (getNumCols() > 2000)
3264            printFrequency_ = 100 ;
3265    }
3266    /*
3267      It is possible that strong branching fixes one variable and then the code goes round
3268      again and again.  This can take too long.  So we need to warn user - just once.
3269    */
3270    numberLongStrong_ = 0;
3271    CbcNode * createdNode = NULL;
[1122]3272#ifdef CBC_THREAD
[1409]3273    if ((specialOptions_&2048) != 0)
3274        numberThreads_ = 0;
3275    if (numberThreads_ ) {
[1286]3276        nodeCompare_->sayThreaded(); // need to use addresses
[1409]3277        master_ = new CbcBaseModel(*this,
3278                                   (parallelMode() < -1) ? 1 : 0);
3279        masterThread_ = master_->masterThread();
[687]3280    }
3281#endif
[931]3282#ifdef COIN_HAS_CLP
[1286]3283    {
3284        OsiClpSolverInterface * clpSolver
3285        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3286        if (clpSolver && !parentModel_) {
3287            clpSolver->computeLargestAway();
3288        }
[931]3289    }
3290#endif
[1286]3291    /*
3292      At last, the actual branch-and-cut search loop, which will iterate until
3293      the live set is empty or we hit some limit (integrality gap, time, node
3294      count, etc.). The overall flow is to rebuild a subproblem, reoptimise using
3295      solveWithCuts(), choose a branching pattern with chooseBranch(), and finally
3296      add the node to the live set.
[2]3297
[1286]3298      The first action is to winnow the live set to remove nodes which are worse
3299      than the current objective cutoff.
3300    */
3301    if (solver_->getRowCutDebuggerAlways()) {
3302        OsiRowCutDebugger * debuggerX = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
3303        const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger() ;
3304        if (!debugger) {
3305            // infeasible!!
3306            printf("before search\n");
3307            const double * lower = solver_->getColLower();
3308            const double * upper = solver_->getColUpper();
3309            const double * solution = debuggerX->optimalSolution();
3310            int numberColumns = solver_->getNumCols();
3311            for (int i = 0; i < numberColumns; i++) {
3312                if (solver_->isInteger(i)) {
3313                    if (solution[i] < lower[i] - 1.0e-6 || solution[i] > upper[i] + 1.0e-6)
3314                        printf("**** ");
3315                    printf("%d %g <= %g <= %g\n",
3316                           i, lower[i], solution[i], upper[i]);
3317                }
3318            }
3319            //abort();
3320        }
[789]3321    }
[1286]3322    {
3323        // may be able to change cutoff now
3324        double cutoff = getCutoff();
3325        double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
3326        if (cutoff > bestObjective_ - increment) {
3327            cutoff = bestObjective_ - increment ;
3328            setCutoff(cutoff) ;
3329        }
[1271]3330    }
[1121]3331#ifdef CBC_THREAD
[1286]3332    bool goneParallel = false;
[1121]3333#endif
[838]3334#define MAX_DEL_NODE 1
[1286]3335    CbcNode * delNode[MAX_DEL_NODE+1];
3336    int nDeleteNode = 0;
3337    // For Printing etc when parallel
3338    int lastEvery1000 = 0;
3339    int lastPrintEvery = 0;
[1315]3340    int numberConsecutiveInfeasible = 0;
[1286]3341    while (true) {
[1409]3342        lockThread();
[931]3343#ifdef COIN_HAS_CLP
[1409]3344        // See if we want dantzig row choice
3345        goToDantzig(100, savePivotMethod);
[1357]3346#endif
[1286]3347        if (tree_->empty()) {
[1122]3348#ifdef CBC_THREAD
[1412]3349            if (parallelMode() > 0 && master_) {
[1409]3350                int anyLeft = master_->waitForThreadsInTree(0);
[1412]3351                if (!anyLeft) {
3352                    master_->stopThreads(-1);
[1409]3353                    break;
[1412]3354                }
[1409]3355            } else {
3356                break;
[1286]3357            }
[1409]3358#else
3359            break;
[1122]3360#endif
[1409]3361        } else {
[1286]3362            unlockThread();
3363        }
3364        // If done 100 nodes see if worth trying reduction
3365        if (numberNodes_ == 50 || numberNodes_ == 100) {
[1271]3366#ifdef COIN_HAS_CLP
[1286]3367            OsiClpSolverInterface * clpSolver
3368            = dynamic_cast<OsiClpSolverInterface *> (solver_);
3369            if (clpSolver && ((specialOptions_&131072) == 0) && true) {
3370                ClpSimplex * simplex = clpSolver->getModelPtr();
3371                int perturbation = simplex->perturbation();
[1132]3372#ifdef CLP_INVESTIGATE
[1286]3373                printf("Testing its n,s %d %d solves n,s %d %d - pert %d\n",
3374                       numberIterations_, saveNumberIterations,
3375                       numberSolves_, saveNumberSolves, perturbation);
[1132]3376#endif
[1286]3377                if (perturbation == 50 && (numberIterations_ - saveNumberIterations) <
3378                        8*(numberSolves_ - saveNumberSolves)) {
3379                    // switch off perturbation
3380                    simplex->setPerturbation(100);
[1271]3381#ifdef PERTURB_IN_FATHOM
[1286]3382                    // but allow in fathom
3383                    specialOptions_ |= 131072;
[1271]3384#endif
3385#ifdef CLP_INVESTIGATE
[1286]3386                    printf("Perturbation switched off\n");
[1271]3387#endif
[1286]3388                }
3389            }