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

Last change on this file since 2003 was 2003, checked in by forrest, 6 years ago

change from stable 2002

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