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

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

more options, copy statistics structure analysis
start coding of "switch" variables i.e. badly scaled ints or hi/lo
changes to allow more influence on small branch and bound
changes to get correct printout with threads

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