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

Last change on this file since 1886 was 1886, checked in by stefan, 7 years ago

fix compiler (gcc 4.6.2) warnings in optimized mode, mainly about unused variables

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