source: stable/2.8/Cbc/src/CbcModel.cpp @ 2082

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

allow restarts with deterministic

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