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

Last change on this file since 1880 was 1880, checked in by forrest, 7 years ago

make it easier to use slow exotic cuts
more on cutoff as constraint and multiple root solver fixes
general fixing of bugs found on MIQP etc

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 714.4 KB
Line 
1/* $Id: CbcModel.cpp 1880 2013-04-01 17:09:22Z forrest $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#if defined(_MSC_VER)
7// Turn off compiler warning about long names
8#  pragma warning(disable:4786)
9#endif
10
11#include "CbcConfig.h"
12
13#include <string>
14//#define CBC_DEBUG 1
15//#define CHECK_CUT_COUNTS
16//#define CHECK_NODE
17//#define CHECK_NODE_FULL
18//#define NODE_LOG
19//#define GLOBAL_CUTS_JUST_POINTERS
20#ifdef CGL_DEBUG_GOMORY
21extern int gomory_try;
22#endif
23#include <cassert>
24#include <cmath>
25#include <cfloat>
26#ifdef COIN_HAS_CLP
27// include Presolve from Clp
28#include "ClpPresolve.hpp"
29#include "OsiClpSolverInterface.hpp"
30#include "ClpNode.hpp"
31#include "ClpDualRowDantzig.hpp"
32#include "ClpSimplexPrimal.hpp"
33#endif
34
35#include "CbcEventHandler.hpp"
36
37#include "OsiSolverInterface.hpp"
38#include "OsiAuxInfo.hpp"
39#include "OsiSolverBranch.hpp"
40#include "OsiChooseVariable.hpp"
41#include "CoinWarmStartBasis.hpp"
42#include "CoinPackedMatrix.hpp"
43#include "CoinHelperFunctions.hpp"
44#include "CbcBranchActual.hpp"
45#include "CbcBranchDynamic.hpp"
46#include "CbcHeuristic.hpp"
47#include "CbcHeuristicFPump.hpp"
48#include "CbcHeuristicRINS.hpp"
49#include "CbcHeuristicDive.hpp"
50#include "CbcModel.hpp"
51#include "CbcTreeLocal.hpp"
52#include "CbcStatistics.hpp"
53#include "CbcStrategy.hpp"
54#include "CbcMessage.hpp"
55#include "OsiRowCut.hpp"
56#include "OsiColCut.hpp"
57#include "OsiRowCutDebugger.hpp"
58#include "OsiCuts.hpp"
59#include "CbcCountRowCut.hpp"
60#include "CbcCutGenerator.hpp"
61#include "CbcFeasibilityBase.hpp"
62#include "CbcFathom.hpp"
63#include "CbcFullNodeInfo.hpp"
64// include Probing
65#include "CglProbing.hpp"
66#include "CglGomory.hpp"
67#include "CglTwomir.hpp"
68// include preprocessing
69#include "CglPreProcess.hpp"
70#include "CglDuplicateRow.hpp"
71#include "CglStored.hpp"
72#include "CglClique.hpp"
73
74#include "CoinTime.hpp"
75#include "CoinMpsIO.hpp"
76
77#include "CbcCompareActual.hpp"
78#include "CbcTree.hpp"
79// This may be dummy
80#include "CbcThread.hpp"
81/* Various functions local to CbcModel.cpp */
82
83static void * doRootCbcThread(void * voidInfo);
84
85namespace {
86
87//-------------------------------------------------------------------
88// Returns the greatest common denominator of two
89// positive integers, a and b, found using Euclid's algorithm
90//-------------------------------------------------------------------
91static int gcd(int a, int b)
92{
93    int remainder = -1;
94    // make sure a<=b (will always remain so)
95    if (a > b) {
96        // Swap a and b
97        int temp = a;
98        a = b;
99        b = temp;
100    }
101    // if zero then gcd is nonzero (zero may occur in rhs of packed)
102    if (!a) {
103        if (b) {
104            return b;
105        } else {
106            printf("**** gcd given two zeros!!\n");
107            abort();
108        }
109    }
110    while (remainder) {
111        remainder = b % a;
112        b = a;
113        a = remainder;
114    }
115    return b;
116}
117
118
119
120#ifdef CHECK_NODE_FULL
121
122/*
123  Routine to verify that tree linkage is correct. The invariant that is tested
124  is
125
126  reference count = (number of actual references) + (number of branches left)
127
128  The routine builds a set of paired arrays, info and count, by traversing the
129  tree. Each CbcNodeInfo is recorded in info, and the number of times it is
130  referenced (via the parent field) is recorded in count. Then a final check is
131  made to see if the numberPointingToThis_ field agrees.
132*/
133
134void verifyTreeNodes (const CbcTree * branchingTree, const CbcModel &model)
135
136{
137    if (model.getNodeCount() == 661) return;
138    printf("*** CHECKING tree after %d nodes\n", model.getNodeCount()) ;
139
140    int j ;
141    int nNodes = branchingTree->size() ;
142# define MAXINFO 1000
143    int *count = new int [MAXINFO] ;
144    CbcNodeInfo **info = new CbcNodeInfo*[MAXINFO] ;
145    int nInfo = 0 ;
146    /*
147      Collect all CbcNodeInfo objects in info, by starting from each live node and
148      traversing back to the root. Nodes in the live set should have unexplored
149      branches remaining.
150
151      TODO: The `while (nodeInfo)' loop could be made to break on reaching a
152        common ancester (nodeInfo is found in info[k]). Alternatively, the
153        check could change to signal an error if nodeInfo is not found above a
154        common ancestor.
155    */
156    for (j = 0 ; j < nNodes ; j++) {
157        CbcNode *node = branchingTree->nodePointer(j) ;
158        if (!node)
159            continue;
160        CbcNodeInfo *nodeInfo = node->nodeInfo() ;
161        int change = node->nodeInfo()->numberBranchesLeft() ;
162        assert(change) ;
163        while (nodeInfo) {
164            int k ;
165            for (k = 0 ; k < nInfo ; k++) {
166                if (nodeInfo == info[k]) break ;
167            }
168            if (k == nInfo) {
169                assert(nInfo < MAXINFO) ;
170                nInfo++ ;
171                info[k] = nodeInfo ;
172                count[k] = 0 ;
173            }
174            nodeInfo = nodeInfo->parent() ;
175        }
176    }
177    /*
178      Walk the info array. For each nodeInfo, look up its parent in info and
179      increment the corresponding count.
180    */
181    for (j = 0 ; j < nInfo ; j++) {
182        CbcNodeInfo *nodeInfo = info[j] ;
183        nodeInfo = nodeInfo->parent() ;
184        if (nodeInfo) {
185            int k ;
186            for (k = 0 ; k < nInfo ; k++) {
187                if (nodeInfo == info[k]) break ;
188            }
189            assert (k < nInfo) ;
190            count[k]++ ;
191        }
192    }
193    /*
194      Walk the info array one more time and check that the invariant holds. The
195      number of references (numberPointingToThis()) should equal the sum of the
196      number of actual references (held in count[]) plus the number of potential
197      references (unexplored branches, numberBranchesLeft()).
198    */
199    for (j = 0; j < nInfo; j++) {
200        CbcNodeInfo * nodeInfo = info[j] ;
201        if (nodeInfo) {
202            int k ;
203            for (k = 0; k < nInfo; k++)
204                if (nodeInfo == info[k])
205                    break ;
206            printf("Nodeinfo %x - %d left, %d count\n",
207                   nodeInfo,
208                   nodeInfo->numberBranchesLeft(),
209                   nodeInfo->numberPointingToThis()) ;
210            assert(nodeInfo->numberPointingToThis() ==
211                   count[k] + nodeInfo->numberBranchesLeft()) ;
212        }
213    }
214
215    delete [] count ;
216    delete [] info ;
217
218    return ;
219}
220
221#endif  /* CHECK_NODE_FULL */
222
223
224
225#ifdef CHECK_CUT_COUNTS
226
227/*
228  Routine to verify that cut reference counts are correct.
229*/
230void verifyCutCounts (const CbcTree * branchingTree, CbcModel &model)
231
232{
233    printf("*** CHECKING cuts after %d nodes\n", model.getNodeCount()) ;
234
235    int j ;
236    int nNodes = branchingTree->size() ;
237
238    /*
239      cut.tempNumber_ exists for the purpose of doing this verification. Clear it
240      in all cuts. We traverse the tree by starting from each live node and working
241      back to the root. At each CbcNodeInfo, check for cuts.
242    */
243    for (j = 0 ; j < nNodes ; j++) {
244        CbcNode *node = branchingTree->nodePointer(j) ;
245        CbcNodeInfo * nodeInfo = node->nodeInfo() ;
246        assert (node->nodeInfo()->numberBranchesLeft()) ;
247        while (nodeInfo) {
248            int k ;
249            for (k = 0 ; k < nodeInfo->numberCuts() ; k++) {
250                CbcCountRowCut *cut = nodeInfo->cuts()[k] ;
251                if (cut) cut->tempNumber_ = 0;
252            }
253            nodeInfo = nodeInfo->parent() ;
254        }
255    }
256    /*
257      Walk the live set again, this time collecting the list of cuts in use at each
258      node. addCuts1 will collect the cuts in model.addedCuts_. Take into account
259      that when we recreate the basis for a node, we compress out the slack cuts.
260    */
261    for (j = 0 ; j < nNodes ; j++) {
262        CoinWarmStartBasis *debugws = model.getEmptyBasis() ;
263        CbcNode *node = branchingTree->nodePointer(j) ;
264        CbcNodeInfo *nodeInfo = node->nodeInfo();
265        int change = node->nodeInfo()->numberBranchesLeft() ;
266        printf("Node %d %x (info %x) var %d way %d obj %g", j, node,
267               node->nodeInfo(), node->columnNumber(), node->way(),
268               node->objectiveValue()) ;
269
270        model.addCuts1(node, debugws) ;
271
272        int i ;
273        int numberRowsAtContinuous = model.numberRowsAtContinuous() ;
274        CbcCountRowCut **addedCuts = model.addedCuts() ;
275        for (i = 0 ; i < model.currentNumberCuts() ; i++) {
276            CoinWarmStartBasis::Status status =
277                debugws->getArtifStatus(i + numberRowsAtContinuous) ;
278            if (status != CoinWarmStartBasis::basic && addedCuts[i]) {
279                addedCuts[i]->tempNumber_ += change ;
280            }
281        }
282
283        while (nodeInfo) {
284            nodeInfo = nodeInfo->parent() ;
285            if (nodeInfo) printf(" -> %x", nodeInfo);
286        }
287        printf("\n") ;
288        delete debugws ;
289    }
290    /*
291      The moment of truth: We've tallied up the references by direct scan of the  search tree. Check for agreement with the count in the cut.
292
293      TODO: Rewrite to check and print mismatch only when tempNumber_ == 0?
294    */
295    for (j = 0 ; j < nNodes ; j++) {
296        CbcNode *node = branchingTree->nodePointer(j) ;
297        CbcNodeInfo *nodeInfo = node->nodeInfo();
298        while (nodeInfo) {
299            int k ;
300            for (k = 0 ; k < nodeInfo->numberCuts() ; k++) {
301                CbcCountRowCut *cut = nodeInfo->cuts()[k] ;
302                if (cut && cut->tempNumber_ >= 0) {
303                    if (cut->tempNumber_ != cut->numberPointingToThis())
304                        printf("mismatch %x %d %x %d %d\n", nodeInfo, k,
305                               cut, cut->tempNumber_, cut->numberPointingToThis()) ;
306                    else
307                        printf("   match %x %d %x %d %d\n", nodeInfo, k,
308                               cut, cut->tempNumber_, cut->numberPointingToThis()) ;
309                    cut->tempNumber_ = -1 ;
310                }
311            }
312            nodeInfo = nodeInfo->parent() ;
313        }
314    }
315
316    return ;
317}
318
319#endif /* CHECK_CUT_COUNTS */
320
321
322#ifdef CHECK_CUT_SIZE
323
324/*
325  Routine to verify that cut reference counts are correct.
326*/
327void verifyCutSize (const CbcTree * branchingTree, CbcModel &model)
328{
329
330    int j ;
331    int nNodes = branchingTree->size() ;
332    int totalCuts = 0;
333
334    /*
335      cut.tempNumber_ exists for the purpose of doing this verification. Clear it
336      in all cuts. We traverse the tree by starting from each live node and working
337      back to the root. At each CbcNodeInfo, check for cuts.
338    */
339    for (j = 0 ; j < nNodes ; j++) {
340        CbcNode *node = branchingTree->nodePointer(j) ;
341        CbcNodeInfo * nodeInfo = node->nodeInfo() ;
342        assert (node->nodeInfo()->numberBranchesLeft()) ;
343        while (nodeInfo) {
344            totalCuts += nodeInfo->numberCuts();
345            nodeInfo = nodeInfo->parent() ;
346        }
347    }
348    printf("*** CHECKING cuts (size) after %d nodes - %d cuts\n", model.getNodeCount(), totalCuts) ;
349    return ;
350}
351
352#endif /* CHECK_CUT_SIZE */
353
354}
355
356/* End unnamed namespace for CbcModel.cpp */
357
358
359void
360CbcModel::analyzeObjective ()
361/*
362  Try to find a minimum change in the objective function. The first scan
363  checks that there are no continuous variables with non-zero coefficients,
364  and grabs the largest objective coefficient associated with an unfixed
365  integer variable. The second scan attempts to scale up the objective
366  coefficients to a point where they are sufficiently close to integer that
367  we can pretend they are integer, and calculate a gcd over the coefficients
368  of interest. This will be the minimum increment for the scaled coefficients.
369  The final action is to scale the increment back for the original coefficients
370  and install it, if it's better than the existing value.
371
372  John's note: We could do better than this.
373
374  John's second note - apologies for changing s to z
375*/
376{
377    const double *objective = getObjCoefficients() ;
378    const double *lower = getColLower() ;
379    const double *upper = getColUpper() ;
380    /*
381      Scan continuous and integer variables to see if continuous
382      are cover or network with integral rhs.
383    */
384    double continuousMultiplier = 1.0;
385    double * coeffMultiplier = NULL;
386    double largestObj = 0.0;
387    double smallestObj = COIN_DBL_MAX;
388    {
389        const double *rowLower = getRowLower() ;
390        const double *rowUpper = getRowUpper() ;
391        int numberRows = solver_->getNumRows() ;
392        double * rhs = new double [numberRows];
393        memset(rhs, 0, numberRows*sizeof(double));
394        int iColumn;
395        int numberColumns = solver_->getNumCols() ;
396        // Column copy of matrix
397        bool allPlusOnes = true;
398        bool allOnes = true;
399        int problemType = -1;
400        const double * element = solver_->getMatrixByCol()->getElements();
401        const int * row = solver_->getMatrixByCol()->getIndices();
402        const CoinBigIndex * columnStart = solver_->getMatrixByCol()->getVectorStarts();
403        const int * columnLength = solver_->getMatrixByCol()->getVectorLengths();
404        int numberInteger = 0;
405        int numberIntegerObj = 0;
406        int numberGeneralIntegerObj = 0;
407        int numberIntegerWeight = 0;
408        int numberContinuousObj = 0;
409        double cost = COIN_DBL_MAX;
410        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
411            if (upper[iColumn] == lower[iColumn]) {
412                CoinBigIndex start = columnStart[iColumn];
413                CoinBigIndex end = start + columnLength[iColumn];
414                for (CoinBigIndex j = start; j < end; j++) {
415                    int iRow = row[j];
416                    rhs[iRow] += lower[iColumn] * element[j];
417                }
418            } else {
419                double objValue = objective[iColumn];
420                if (solver_->isInteger(iColumn))
421                    numberInteger++;
422                if (objValue) {
423                    if (!solver_->isInteger(iColumn)) {
424                        numberContinuousObj++;
425                    } else {
426                        largestObj = CoinMax(largestObj, fabs(objValue));
427                        smallestObj = CoinMin(smallestObj, fabs(objValue));
428                        numberIntegerObj++;
429                        if (cost == COIN_DBL_MAX)
430                            cost = objValue;
431                        else if (cost != objValue)
432                            cost = -COIN_DBL_MAX;
433                        int gap = static_cast<int> (upper[iColumn] - lower[iColumn]);
434                        if (gap > 1) {
435                            numberGeneralIntegerObj++;
436                            numberIntegerWeight += gap;
437                        }
438                    }
439                }
440            }
441        }
442        int iType = 0;
443        if (!numberContinuousObj && numberIntegerObj <= 5 && numberIntegerWeight <= 100 &&
444                numberIntegerObj*3 < numberObjects_ && !parentModel_ && solver_->getNumRows() > 100)
445          iType = 1 + 4 + ((moreSpecialOptions_&536870912)==0) ? 2 : 0;
446        else if (!numberContinuousObj && numberIntegerObj <= 100 &&
447                 numberIntegerObj*5 < numberObjects_ && numberIntegerWeight <= 100 &&
448                 !parentModel_ &&
449                 solver_->getNumRows() > 100 && cost != -COIN_DBL_MAX)
450          iType = 4 + ((moreSpecialOptions_&536870912)==0) ? 2 : 0;
451        else if (!numberContinuousObj && numberIntegerObj <= 100 &&
452                 numberIntegerObj*5 < numberObjects_ &&
453                 !parentModel_ &&
454                 solver_->getNumRows() > 100 && cost != -COIN_DBL_MAX)
455            iType = 8;
456        int iTest = getMaximumNodes();
457        if (iTest >= 987654320 && iTest < 987654330 && numberObjects_ && !parentModel_) {
458            iType = iTest - 987654320;
459            printf("Testing %d integer variables out of %d objects (%d integer) have cost of %g - %d continuous\n",
460                   numberIntegerObj, numberObjects_, numberInteger, cost, numberContinuousObj);
461            if (iType == 9)
462                exit(77);
463            if (numberContinuousObj)
464                iType = 0;
465        }
466
467        //if (!numberContinuousObj&&(numberIntegerObj<=5||cost!=-COIN_DBL_MAX)&&
468        //numberIntegerObj*3<numberObjects_&&!parentModel_&&solver_->getNumRows()>100) {
469        if (iType) {
470            /*
471            A) put high priority on (if none)
472            B) create artificial objective (if clp)
473            */
474            int iPriority = -1;
475            for (int i = 0; i < numberObjects_; i++) {
476                int k = object_[i]->priority();
477                if (iPriority == -1)
478                    iPriority = k;
479                else if (iPriority != k)
480                    iPriority = -2;
481            }
482            bool branchOnSatisfied = ((iType & 1) != 0);
483            bool createFake = ((iType & 2) != 0);
484            bool randomCost = ((iType & 4) != 0);
485            if (iPriority >= 0) {
486                char general[200];
487                if (cost == -COIN_DBL_MAX) {
488                    sprintf(general, "%d integer variables out of %d objects (%d integer) have costs - high priority",
489                            numberIntegerObj, numberObjects_, numberInteger);
490                } else if (cost == COIN_DBL_MAX) {
491                    sprintf(general, "No integer variables out of %d objects (%d integer) have costs",
492                            numberObjects_, numberInteger);
493                    branchOnSatisfied = false;
494                } else {
495                    sprintf(general, "%d integer variables out of %d objects (%d integer) have cost of %g - high priority",
496                            numberIntegerObj, numberObjects_, numberInteger, cost);
497                }
498                messageHandler()->message(CBC_GENERAL,
499                                          messages())
500                << general << CoinMessageEol ;
501                sprintf(general, "branch on satisfied %c create fake objective %c random cost %c",
502                        branchOnSatisfied ? 'Y' : 'N',
503                        createFake ? 'Y' : 'N',
504                        randomCost ? 'Y' : 'N');
505                messageHandler()->message(CBC_GENERAL,
506                                          messages())
507                << general << CoinMessageEol ;
508                // switch off clp type branching
509                // no ? fastNodeDepth_ = -1;
510                int highPriority = (branchOnSatisfied) ? -999 : 100;
511                for (int i = 0; i < numberObjects_; i++) {
512                    CbcSimpleInteger * thisOne = dynamic_cast <CbcSimpleInteger *> (object_[i]);
513                    object_[i]->setPriority(1000);
514                    if (thisOne) {
515                        int iColumn = thisOne->columnNumber();
516                        if (objective[iColumn])
517                            thisOne->setPriority(highPriority);
518                    }
519                }
520            }
521#ifdef COIN_HAS_CLP
522            OsiClpSolverInterface * clpSolver
523            = dynamic_cast<OsiClpSolverInterface *> (solver_);
524            if (clpSolver && createFake) {
525                // Create artificial objective to be used when all else fixed
526                int numberColumns = clpSolver->getNumCols();
527                double * fakeObj = new double [numberColumns];
528                // Column copy
529                const CoinPackedMatrix  * matrixByCol = clpSolver->getMatrixByCol();
530                //const double * element = matrixByCol.getElements();
531                //const int * row = matrixByCol.getIndices();
532                //const CoinBigIndex * columnStart = matrixByCol.getVectorStarts();
533                const int * columnLength = matrixByCol->getVectorLengths();
534                const double * solution = clpSolver->getColSolution();
535#ifdef JJF_ZERO
536                int nAtBound = 0;
537                for (int i = 0; i < numberColumns; i++) {
538                    double lowerValue = lower[i];
539                    double upperValue = upper[i];
540                    if (clpSolver->isInteger(i)) {
541                        double lowerValue = lower[i];
542                        double upperValue = upper[i];
543                        double value = solution[i];
544                        if (value < lowerValue + 1.0e-6 ||
545                                value > upperValue - 1.0e-6)
546                            nAtBound++;
547                    }
548                }
549#endif
550                /*
551                  Generate a random objective function for problems where the given objective
552                  function is not terribly useful. (Nearly feasible, single integer variable,
553                  that sort of thing.
554                */
555                CoinDrand48(true, 1234567);
556                for (int i = 0; i < numberColumns; i++) {
557                    double lowerValue = lower[i];
558                    double upperValue = upper[i];
559                    double value = (randomCost) ? ceil((CoinDrand48() + 0.5) * 1000)
560                                   : i + 1 + columnLength[i] * 1000;
561                    value *= 0.001;
562                    //value += columnLength[i];
563                    if (lowerValue > -1.0e5 || upperValue < 1.0e5) {
564                        if (fabs(lowerValue) > fabs(upperValue))
565                            value = - value;
566                        if (clpSolver->isInteger(i)) {
567                            double solValue = solution[i];
568                            // Better to add in 0.5 or 1.0??
569                            if (solValue < lowerValue + 1.0e-6)
570                                value = fabs(value) + 0.5; //fabs(value*1.5);
571                            else if (solValue > upperValue - 1.0e-6)
572                                value = -fabs(value) - 0.5; //-fabs(value*1.5);
573                        }
574                    } else {
575                        value = 0.0;
576                    }
577                    fakeObj[i] = value;
578                }
579                // pass to solver
580                clpSolver->setFakeObjective(fakeObj);
581                delete [] fakeObj;
582            }
583#endif
584        } else if (largestObj < smallestObj*5.0 && !parentModel_ &&
585                   !numberContinuousObj &&
586                   !numberGeneralIntegerObj &&
587                   numberIntegerObj*2 < numberColumns) {
588            // up priorities on costed
589            int iPriority = -1;
590            for (int i = 0; i < numberObjects_; i++) {
591                int k = object_[i]->priority();
592                if (iPriority == -1)
593                    iPriority = k;
594                else if (iPriority != k)
595                    iPriority = -2;
596            }
597            if (iPriority >= 100) {
598#ifdef CLP_INVESTIGATE
599                printf("Setting variables with obj to high priority\n");
600#endif
601                for (int i = 0; i < numberObjects_; i++) {
602                    CbcSimpleInteger * obj =
603                        dynamic_cast <CbcSimpleInteger *>(object_[i]) ;
604                    if (obj) {
605                        int iColumn = obj->columnNumber();
606                        if (objective[iColumn])
607                            object_[i]->setPriority(iPriority - 1);
608                    }
609                }
610            }
611        }
612        int iRow;
613        for (iRow = 0; iRow < numberRows; iRow++) {
614            if (rowLower[iRow] > -1.0e20 &&
615                    fabs(rowLower[iRow] - rhs[iRow] - floor(rowLower[iRow] - rhs[iRow] + 0.5)) > 1.0e-10) {
616                continuousMultiplier = 0.0;
617                break;
618            }
619            if (rowUpper[iRow] < 1.0e20 &&
620                    fabs(rowUpper[iRow] - rhs[iRow] - floor(rowUpper[iRow] - rhs[iRow] + 0.5)) > 1.0e-10) {
621                continuousMultiplier = 0.0;
622                break;
623            }
624            // set rhs to limiting value
625            if (rowLower[iRow] != rowUpper[iRow]) {
626                if (rowLower[iRow] > -1.0e20) {
627                    if (rowUpper[iRow] < 1.0e20) {
628                        // no good
629                        continuousMultiplier = 0.0;
630                        break;
631                    } else {
632                        rhs[iRow] = rowLower[iRow] - rhs[iRow];
633                        if (problemType < 0)
634                            problemType = 3; // set cover
635                        else if (problemType != 3)
636                            problemType = 4;
637                    }
638                } else {
639                    rhs[iRow] = rowUpper[iRow] - rhs[iRow];
640                    if (problemType < 0)
641                        problemType = 1; // set partitioning <=
642                    else if (problemType != 1)
643                        problemType = 4;
644                }
645            } else {
646                rhs[iRow] = rowUpper[iRow] - rhs[iRow];
647                if (problemType < 0)
648                    problemType = 3; // set partitioning ==
649                else if (problemType != 2)
650                    problemType = 2;
651            }
652            if (fabs(rhs[iRow] - 1.0) > 1.0e-12)
653                problemType = 4;
654        }
655        if (continuousMultiplier) {
656            // 1 network, 2 cover, 4 negative cover
657            int possible = 7;
658            bool unitRhs = true;
659            // See which rows could be set cover
660            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
661                if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
662                    CoinBigIndex start = columnStart[iColumn];
663                    CoinBigIndex end = start + columnLength[iColumn];
664                    for (CoinBigIndex j = start; j < end; j++) {
665                        double value = element[j];
666                        if (value == 1.0) {
667                        } else if (value == -1.0) {
668                            rhs[row[j]] = -0.5;
669                            allPlusOnes = false;
670                        } else {
671                            rhs[row[j]] = -COIN_DBL_MAX;
672                            allOnes = false;
673                        }
674                    }
675                }
676            }
677            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
678                if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
679                    if (!isInteger(iColumn)) {
680                        CoinBigIndex start = columnStart[iColumn];
681                        CoinBigIndex end = start + columnLength[iColumn];
682                        double rhsValue = 0.0;
683                        // 1 all ones, -1 all -1s, 2 all +- 1, 3 no good
684                        int type = 0;
685                        for (CoinBigIndex j = start; j < end; j++) {
686                            double value = element[j];
687                            if (fabs(value) != 1.0) {
688                                type = 3;
689                                break;
690                            } else if (value == 1.0) {
691                                if (!type)
692                                    type = 1;
693                                else if (type != 1)
694                                    type = 2;
695                            } else {
696                                if (!type)
697                                    type = -1;
698                                else if (type != -1)
699                                    type = 2;
700                            }
701                            int iRow = row[j];
702                            if (rhs[iRow] == -COIN_DBL_MAX) {
703                                type = 3;
704                                break;
705                            } else if (rhs[iRow] == -0.5) {
706                                // different values
707                                unitRhs = false;
708                            } else if (rhsValue) {
709                                if (rhsValue != rhs[iRow])
710                                    unitRhs = false;
711                            } else {
712                                rhsValue = rhs[iRow];
713                            }
714                        }
715                        // if no elements OK
716                        if (type == 3) {
717                            // no good
718                            possible = 0;
719                            break;
720                        } else if (type == 2) {
721                            if (end - start > 2) {
722                                // no good
723                                possible = 0;
724                                break;
725                            } else {
726                                // only network
727                                possible &= 1;
728                                if (!possible)
729                                    break;
730                            }
731                        } else if (type == 1) {
732                            // only cover
733                            possible &= 2;
734                            if (!possible)
735                                break;
736                        } else if (type == -1) {
737                            // only negative cover
738                            possible &= 4;
739                            if (!possible)
740                                break;
741                        }
742                    }
743                }
744            }
745            if ((possible == 2 || possible == 4) && !unitRhs) {
746#if COIN_DEVELOP>1
747                printf("XXXXXX Continuous all +1 but different rhs\n");
748#endif
749                possible = 0;
750            }
751            // may be all integer
752            if (possible != 7) {
753                if (!possible)
754                    continuousMultiplier = 0.0;
755                else if (possible == 1)
756                    continuousMultiplier = 1.0;
757                else
758                    continuousMultiplier = 0.0; // 0.5 was incorrect;
759#if COIN_DEVELOP>1
760                if (continuousMultiplier)
761                    printf("XXXXXX multiplier of %g\n", continuousMultiplier);
762#endif
763                if (continuousMultiplier == 0.5) {
764                    coeffMultiplier = new double [numberColumns];
765                    bool allOne = true;
766                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
767                        coeffMultiplier[iColumn] = 1.0;
768                        if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
769                            if (!isInteger(iColumn)) {
770                                CoinBigIndex start = columnStart[iColumn];
771                                int iRow = row[start];
772                                double value = rhs[iRow];
773                                assert (value >= 0.0);
774                                if (value != 0.0 && value != 1.0)
775                                    allOne = false;
776                                coeffMultiplier[iColumn] = 0.5 * value;
777                            }
778                        }
779                    }
780                    if (allOne) {
781                        // back to old way
782                        delete [] coeffMultiplier;
783                        coeffMultiplier = NULL;
784                    }
785                }
786            } else {
787                // all integer
788                problemType_ = problemType;
789#if COIN_DEVELOP>1
790                printf("Problem type is %d\n", problemType_);
791#endif
792            }
793        }
794
795        // But try again
796        if (continuousMultiplier < 1.0) {
797            memset(rhs, 0, numberRows*sizeof(double));
798            int * count = new int [numberRows];
799            memset(count, 0, numberRows*sizeof(int));
800            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
801                CoinBigIndex start = columnStart[iColumn];
802                CoinBigIndex end = start + columnLength[iColumn];
803                if (upper[iColumn] == lower[iColumn]) {
804                    for (CoinBigIndex j = start; j < end; j++) {
805                        int iRow = row[j];
806                        rhs[iRow] += lower[iColumn] * element[j];
807                    }
808                } else if (solver_->isInteger(iColumn)) {
809                    for (CoinBigIndex j = start; j < end; j++) {
810                        int iRow = row[j];
811                        if (fabs(element[j] - floor(element[j] + 0.5)) > 1.0e-10)
812                            rhs[iRow]  = COIN_DBL_MAX;
813                    }
814                } else {
815                    for (CoinBigIndex j = start; j < end; j++) {
816                        int iRow = row[j];
817                        count[iRow]++;
818                        if (fabs(element[j]) != 1.0)
819                            rhs[iRow]  = COIN_DBL_MAX;
820                    }
821                }
822            }
823            // now look at continuous
824            bool allGood = true;
825            double direction = solver_->getObjSense() ;
826            int numberObj = 0;
827            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
828                if (upper[iColumn] > lower[iColumn]) {
829                    double objValue = objective[iColumn] * direction;
830                    if (objValue && !solver_->isInteger(iColumn)) {
831                        numberObj++;
832                        CoinBigIndex start = columnStart[iColumn];
833                        CoinBigIndex end = start + columnLength[iColumn];
834                        if (objValue > 0.0) {
835                            // wants to be as low as possible
836                            if (lower[iColumn] < -1.0e10 || fabs(lower[iColumn] - floor(lower[iColumn] + 0.5)) > 1.0e-10) {
837                                allGood = false;
838                                break;
839                            } else if (upper[iColumn] < 1.0e10 && fabs(upper[iColumn] - floor(upper[iColumn] + 0.5)) > 1.0e-10) {
840                                allGood = false;
841                                break;
842                            }
843                            bool singletonRow = true;
844                            bool equality = false;
845                            for (CoinBigIndex j = start; j < end; j++) {
846                                int iRow = row[j];
847                                if (count[iRow] > 1)
848                                    singletonRow = false;
849                                else if (rowLower[iRow] == rowUpper[iRow])
850                                    equality = true;
851                                double rhsValue = rhs[iRow];
852                                double lowerValue = rowLower[iRow];
853                                double upperValue = rowUpper[iRow];
854                                if (rhsValue < 1.0e20) {
855                                    if (lowerValue > -1.0e20)
856                                        lowerValue -= rhsValue;
857                                    if (upperValue < 1.0e20)
858                                        upperValue -= rhsValue;
859                                }
860                                if (fabs(rhsValue) > 1.0e20 || fabs(rhsValue - floor(rhsValue + 0.5)) > 1.0e-10
861                                        || fabs(element[j]) != 1.0) {
862                                    // no good
863                                    allGood = false;
864                                    break;
865                                }
866                                if (element[j] > 0.0) {
867                                    if (lowerValue > -1.0e20 && fabs(lowerValue - floor(lowerValue + 0.5)) > 1.0e-10) {
868                                        // no good
869                                        allGood = false;
870                                        break;
871                                    }
872                                } else {
873                                    if (upperValue < 1.0e20 && fabs(upperValue - floor(upperValue + 0.5)) > 1.0e-10) {
874                                        // no good
875                                        allGood = false;
876                                        break;
877                                    }
878                                }
879                            }
880                            if (!singletonRow && end > start + 1 && !equality)
881                                allGood = false;
882                            if (!allGood)
883                                break;
884                        } else {
885                            // wants to be as high as possible
886                            if (upper[iColumn] > 1.0e10 || fabs(upper[iColumn] - floor(upper[iColumn] + 0.5)) > 1.0e-10) {
887                                allGood = false;
888                                break;
889                            } else if (lower[iColumn] > -1.0e10 && fabs(lower[iColumn] - floor(lower[iColumn] + 0.5)) > 1.0e-10) {
890                                allGood = false;
891                                break;
892                            }
893                            bool singletonRow = true;
894                            bool equality = false;
895                            for (CoinBigIndex j = start; j < end; j++) {
896                                int iRow = row[j];
897                                if (count[iRow] > 1)
898                                    singletonRow = false;
899                                else if (rowLower[iRow] == rowUpper[iRow])
900                                    equality = true;
901                                double rhsValue = rhs[iRow];
902                                double lowerValue = rowLower[iRow];
903                                double upperValue = rowUpper[iRow];
904                                if (rhsValue < 1.0e20) {
905                                    if (lowerValue > -1.0e20)
906                                        lowerValue -= rhsValue;
907                                    if (upperValue < 1.0e20)
908                                        upperValue -= rhsValue;
909                                }
910                                if (fabs(rhsValue) > 1.0e20 || fabs(rhsValue - floor(rhsValue + 0.5)) > 1.0e-10
911                                        || fabs(element[j]) != 1.0) {
912                                    // no good
913                                    allGood = false;
914                                    break;
915                                }
916                                if (element[j] < 0.0) {
917                                    if (lowerValue > -1.0e20 && fabs(lowerValue - floor(lowerValue + 0.5)) > 1.0e-10) {
918                                        // no good
919                                        allGood = false;
920                                        break;
921                                    }
922                                } else {
923                                    if (upperValue < 1.0e20 && fabs(upperValue - floor(upperValue + 0.5)) > 1.0e-10) {
924                                        // no good
925                                        allGood = false;
926                                        break;
927                                    }
928                                }
929                            }
930                            if (!singletonRow && end > start + 1 && !equality)
931                                allGood = false;
932                            if (!allGood)
933                                break;
934                        }
935                    }
936                }
937            }
938            delete [] count;
939            if (allGood) {
940#if COIN_DEVELOP>1
941                if (numberObj)
942                    printf("YYYY analysis says all continuous with costs will be integer\n");
943#endif
944                continuousMultiplier = 1.0;
945            }
946        }
947        delete [] rhs;
948    }
949    /*
950      Take a first scan to see if there are unfixed continuous variables in the
951      objective.  If so, the minimum objective change could be arbitrarily small.
952      Also pick off the maximum coefficient of an unfixed integer variable.
953
954      If the objective is found to contain only integer variables, set the
955      fathoming discipline to strict.
956    */
957    double maximumCost = 0.0 ;
958    //double trueIncrement=0.0;
959    int iColumn ;
960    int numberColumns = getNumCols() ;
961    double scaleFactor = 1.0; // due to rhs etc
962    /*
963      Original model did not have integer bounds.
964    */
965    if ((specialOptions_&65536) == 0) {
966        /* be on safe side (later look carefully as may be able to
967           to get 0.5 say if bounds are multiples of 0.5 */
968        for (iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
969            if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
970                double value;
971                value = fabs(lower[iColumn]);
972                if (floor(value + 0.5) != value) {
973                    scaleFactor = CoinMin(scaleFactor, 0.5);
974                    if (floor(2.0*value + 0.5) != 2.0*value) {
975                        scaleFactor = CoinMin(scaleFactor, 0.25);
976                        if (floor(4.0*value + 0.5) != 4.0*value) {
977                            scaleFactor = 0.0;
978                        }
979                    }
980                }
981                value = fabs(upper[iColumn]);
982                if (floor(value + 0.5) != value) {
983                    scaleFactor = CoinMin(scaleFactor, 0.5);
984                    if (floor(2.0*value + 0.5) != 2.0*value) {
985                        scaleFactor = CoinMin(scaleFactor, 0.25);
986                        if (floor(4.0*value + 0.5) != 4.0*value) {
987                            scaleFactor = 0.0;
988                        }
989                    }
990                }
991            }
992        }
993    }
994    bool possibleMultiple = continuousMultiplier != 0.0 && scaleFactor != 0.0 ;
995    if (possibleMultiple) {
996        for (iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
997            if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
998                maximumCost = CoinMax(maximumCost, fabs(objective[iColumn])) ;
999            }
1000        }
1001    }
1002    setIntParam(CbcModel::CbcFathomDiscipline, possibleMultiple) ;
1003    /*
1004      If a nontrivial increment is possible, try and figure it out. We're looking
1005      for gcd(c<j>) for all c<j> that are coefficients of unfixed integer
1006      variables. Since the c<j> might not be integers, try and inflate them
1007      sufficiently that they look like integers (and we'll deflate the gcd
1008      later).
1009
1010      2520.0 is used as it is a nice multiple of 2,3,5,7
1011    */
1012    if (possibleMultiple && maximumCost) {
1013        int increment = 0 ;
1014        double multiplier = 2520.0 ;
1015        while (10.0*multiplier*maximumCost < 1.0e8)
1016            multiplier *= 10.0 ;
1017        int bigIntegers = 0; // Count of large costs which are integer
1018        for (iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
1019            if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
1020                double objValue = fabs(objective[iColumn]);
1021                if (!isInteger(iColumn)) {
1022                    if (!coeffMultiplier)
1023                        objValue *= continuousMultiplier;
1024                    else
1025                        objValue *= coeffMultiplier[iColumn];
1026                }
1027                if (objValue) {
1028                    double value = objValue * multiplier ;
1029                    if (value < 2.1e9) {
1030                        int nearest = static_cast<int> (floor(value + 0.5)) ;
1031                        if (fabs(value - floor(value + 0.5)) > 1.0e-8) {
1032                            increment = 0 ;
1033                            break ;
1034                        } else if (!increment) {
1035                            increment = nearest ;
1036                        } else {
1037                            increment = gcd(increment, nearest) ;
1038                        }
1039                    } else {
1040                        // large value - may still be multiple of 1.0
1041                        if (fabs(objValue - floor(objValue + 0.5)) > 1.0e-8) {
1042                            increment = 0;
1043                            break;
1044                        } else {
1045                            bigIntegers++;
1046                        }
1047                    }
1048                }
1049            }
1050        }
1051        delete [] coeffMultiplier;
1052        /*
1053          If the increment beats the current value for objective change, install it.
1054        */
1055        if (increment) {
1056            double value = increment ;
1057            double cutoff = getDblParam(CbcModel::CbcCutoffIncrement) ;
1058            if (bigIntegers) {
1059                // allow for 1.0
1060                increment = gcd(increment, static_cast<int> (multiplier));
1061                value = increment;
1062            }
1063            value /= multiplier ;
1064            value *= scaleFactor;
1065            //trueIncrement=CoinMax(cutoff,value);;
1066            if (value*0.999 > cutoff) {
1067                messageHandler()->message(CBC_INTEGERINCREMENT,
1068                                          messages())
1069                << value << CoinMessageEol ;
1070                setDblParam(CbcModel::CbcCutoffIncrement, CoinMax(value*0.999,value-1.0e-4)) ;
1071            }
1072        }
1073    }
1074
1075    return ;
1076}
1077
1078/*
1079saveModel called (carved out of) BranchandBound
1080*/
1081void CbcModel::saveModel(OsiSolverInterface * saveSolver, double * checkCutoffForRestart, bool * feasible)
1082{
1083    if (saveSolver && (specialOptions_&32768) != 0) {
1084        // See if worth trying reduction
1085        *checkCutoffForRestart = getCutoff();
1086        bool tryNewSearch = solverCharacteristics_->reducedCostsAccurate() &&
1087                            (*checkCutoffForRestart < 1.0e20);
1088        int numberColumns = getNumCols();
1089        if (tryNewSearch) {
1090#ifdef CLP_INVESTIGATE
1091            printf("after %d nodes, cutoff %g - looking\n",
1092                   numberNodes_, getCutoff());
1093#endif
1094            saveSolver->resolve();
1095            double direction = saveSolver->getObjSense() ;
1096            double gap = *checkCutoffForRestart - saveSolver->getObjValue() * direction ;
1097            double tolerance;
1098            saveSolver->getDblParam(OsiDualTolerance, tolerance) ;
1099            if (gap <= 0.0)
1100                gap = tolerance;
1101            gap += 100.0 * tolerance;
1102            double integerTolerance = getDblParam(CbcIntegerTolerance) ;
1103
1104            const double *lower = saveSolver->getColLower() ;
1105            const double *upper = saveSolver->getColUpper() ;
1106            const double *solution = saveSolver->getColSolution() ;
1107            const double *reducedCost = saveSolver->getReducedCost() ;
1108
1109            int numberFixed = 0 ;
1110            int numberFixed2 = 0;
1111            for (int i = 0 ; i < numberIntegers_ ; i++) {
1112                int iColumn = integerVariable_[i] ;
1113                double djValue = direction * reducedCost[iColumn] ;
1114                if (upper[iColumn] - lower[iColumn] > integerTolerance) {
1115                    if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) {
1116                        saveSolver->setColUpper(iColumn, lower[iColumn]) ;
1117                        numberFixed++ ;
1118                    } else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) {
1119                        saveSolver->setColLower(iColumn, upper[iColumn]) ;
1120                        numberFixed++ ;
1121                    }
1122                } else {
1123                    numberFixed2++;
1124                }
1125            }
1126#ifdef COIN_DEVELOP
1127            /*
1128              We're debugging. (specialOptions 1)
1129            */
1130            if ((specialOptions_&1) != 0) {
1131                const OsiRowCutDebugger *debugger = saveSolver->getRowCutDebugger() ;
1132                if (debugger) {
1133                    printf("Contains optimal\n") ;
1134                    OsiSolverInterface * temp = saveSolver->clone();
1135                    const double * solution = debugger->optimalSolution();
1136                    const double *lower = temp->getColLower() ;
1137                    const double *upper = temp->getColUpper() ;
1138                    int n = temp->getNumCols();
1139                    for (int i = 0; i < n; i++) {
1140                        if (temp->isInteger(i)) {
1141                            double value = floor(solution[i] + 0.5);
1142                            assert (value >= lower[i] && value <= upper[i]);
1143                            temp->setColLower(i, value);
1144                            temp->setColUpper(i, value);
1145                        }
1146                    }
1147                    temp->writeMps("reduced_fix");
1148                    delete temp;
1149                    saveSolver->writeMps("reduced");
1150                } else {
1151                    abort();
1152                }
1153            }
1154            printf("Restart could fix %d integers (%d already fixed)\n",
1155                   numberFixed + numberFixed2, numberFixed2);
1156#endif
1157            numberFixed += numberFixed2;
1158            if (numberFixed*20 < numberColumns)
1159                tryNewSearch = false;
1160        }
1161        if (tryNewSearch) {
1162            // back to solver without cuts?
1163            OsiSolverInterface * solver2 = continuousSolver_->clone();
1164            const double *lower = saveSolver->getColLower() ;
1165            const double *upper = saveSolver->getColUpper() ;
1166            for (int i = 0 ; i < numberIntegers_ ; i++) {
1167                int iColumn = integerVariable_[i] ;
1168                solver2->setColLower(iColumn, lower[iColumn]);
1169                solver2->setColUpper(iColumn, upper[iColumn]);
1170            }
1171            // swap
1172            delete saveSolver;
1173            saveSolver = solver2;
1174            double * newSolution = new double[numberColumns];
1175            double objectiveValue = *checkCutoffForRestart;
1176            CbcSerendipity heuristic(*this);
1177            if (bestSolution_)
1178                heuristic.setInputSolution(bestSolution_, bestObjective_);
1179            heuristic.setFractionSmall(0.9);
1180            heuristic.setFeasibilityPumpOptions(1008013);
1181            // Use numberNodes to say how many are original rows
1182            heuristic.setNumberNodes(continuousSolver_->getNumRows());
1183#ifdef COIN_DEVELOP
1184            if (continuousSolver_->getNumRows() <
1185                    saveSolver->getNumRows())
1186                printf("%d rows added ZZZZZ\n",
1187                       solver_->getNumRows() - continuousSolver_->getNumRows());
1188#endif
1189            int returnCode = heuristic.smallBranchAndBound(saveSolver,
1190                             -1, newSolution,
1191                             objectiveValue,
1192                             *checkCutoffForRestart, "Reduce");
1193            if (returnCode < 0) {
1194#ifdef COIN_DEVELOP
1195                printf("Restart - not small enough to do search after fixing\n");
1196#endif
1197                delete [] newSolution;
1198            } else {
1199                if ((returnCode&1) != 0) {
1200                    // increment number of solutions so other heuristics can test
1201                    numberSolutions_++;
1202                    numberHeuristicSolutions_++;
1203                    lastHeuristic_ = NULL;
1204                    setBestSolution(CBC_ROUNDING, objectiveValue, newSolution) ;
1205                }
1206                delete [] newSolution;
1207                *feasible = false; // stop search
1208            }
1209#if 0 // probably not needed def CBC_THREAD
1210            if (master_) {
1211                lockThread();
1212                if (parallelMode() > 0) {
1213                    while (master_->waitForThreadsInTree(0)) {
1214                        lockThread();
1215                        double dummyBest;
1216                        tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
1217                        //unlockThread();
1218                    }
1219                }
1220                master_->waitForThreadsInTree(2);
1221                delete master_;
1222                master_ = NULL;
1223                masterThread_ = NULL;
1224            }
1225#endif
1226        }
1227    }
1228}
1229/*
1230Adds integers, called from BranchandBound()
1231*/
1232void CbcModel::AddIntegers()
1233{
1234    int numberColumns = continuousSolver_->getNumCols();
1235    int numberRows = continuousSolver_->getNumRows();
1236    int * del = new int [CoinMax(numberColumns, numberRows)];
1237    int * original = new int [numberColumns];
1238    char * possibleRow = new char [numberRows];
1239    {
1240        const CoinPackedMatrix * rowCopy = continuousSolver_->getMatrixByRow();
1241        const int * column = rowCopy->getIndices();
1242        const int * rowLength = rowCopy->getVectorLengths();
1243        const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
1244        const double * rowLower = continuousSolver_->getRowLower();
1245        const double * rowUpper = continuousSolver_->getRowUpper();
1246        const double * element = rowCopy->getElements();
1247        for (int i = 0; i < numberRows; i++) {
1248            int nLeft = 0;
1249            bool possible = false;
1250            if (rowLower[i] < -1.0e20) {
1251                double value = rowUpper[i];
1252                if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1253                    possible = true;
1254            } else if (rowUpper[i] > 1.0e20) {
1255                double value = rowLower[i];
1256                if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1257                    possible = true;
1258            } else {
1259                double value = rowUpper[i];
1260                if (rowLower[i] == rowUpper[i] &&
1261                        fabs(value - floor(value + 0.5)) < 1.0e-8)
1262                    possible = true;
1263            }
1264            double allSame = (possible) ? 0.0 : -1.0;
1265            for (CoinBigIndex j = rowStart[i];
1266                    j < rowStart[i] + rowLength[i]; j++) {
1267                int iColumn = column[j];
1268                if (continuousSolver_->isInteger(iColumn)) {
1269                    if (fabs(element[j]) != 1.0)
1270                        possible = false;
1271                } else {
1272                    nLeft++;
1273                    if (!allSame) {
1274                      allSame = fabs(element[j]);
1275                    } else if (allSame>0.0) {
1276                      if (allSame!=fabs(element[j]))
1277                        allSame = -1.0;
1278                    }
1279                }
1280            }
1281            if (nLeft == rowLength[i] && allSame > 0.0)
1282                possibleRow[i] = 2;
1283            else if (possible || !nLeft)
1284                possibleRow[i] = 1;
1285            else
1286                possibleRow[i] = 0;
1287        }
1288    }
1289    int nDel = 0;
1290    for (int i = 0; i < numberColumns; i++) {
1291        original[i] = i;
1292        if (continuousSolver_->isInteger(i))
1293            del[nDel++] = i;
1294    }
1295    int nExtra = 0;
1296    OsiSolverInterface * copy1 = continuousSolver_->clone();
1297    int nPass = 0;
1298    while (nDel && nPass < 10) {
1299        nPass++;
1300        OsiSolverInterface * copy2 = copy1->clone();
1301        int nLeft = 0;
1302        for (int i = 0; i < nDel; i++)
1303            original[del[i]] = -1;
1304        for (int i = 0; i < numberColumns; i++) {
1305            int kOrig = original[i];
1306            if (kOrig >= 0)
1307                original[nLeft++] = kOrig;
1308        }
1309        assert (nLeft == numberColumns - nDel);
1310        copy2->deleteCols(nDel, del);
1311        numberColumns = copy2->getNumCols();
1312        const CoinPackedMatrix * rowCopy = copy2->getMatrixByRow();
1313        numberRows = rowCopy->getNumRows();
1314        const int * column = rowCopy->getIndices();
1315        const int * rowLength = rowCopy->getVectorLengths();
1316        const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
1317        const double * rowLower = copy2->getRowLower();
1318        const double * rowUpper = copy2->getRowUpper();
1319        const double * element = rowCopy->getElements();
1320        const CoinPackedMatrix * columnCopy = copy2->getMatrixByCol();
1321        const int * columnLength = columnCopy->getVectorLengths();
1322        nDel = 0;
1323        // Could do gcd stuff on ones with costs
1324        for (int i = 0; i < numberRows; i++) {
1325            if (!rowLength[i]) {
1326                del[nDel++] = i;
1327            } else if (possibleRow[i]) {
1328                if (rowLength[i] == 1) {
1329                    int k = rowStart[i];
1330                    int iColumn = column[k];
1331                    if (!copy2->isInteger(iColumn)) {
1332                        double mult = 1.0 / fabs(element[k]);
1333                        if (rowLower[i] < -1.0e20) {
1334                            // treat rhs as multiple of 1 unless elements all same
1335                            double value = ((possibleRow[i]==2) ? rowUpper[i] : 1.0) * mult;
1336                            if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
1337                                del[nDel++] = i;
1338                                if (columnLength[iColumn] == 1) {
1339                                    copy2->setInteger(iColumn);
1340                                    int kOrig = original[iColumn];
1341                                    setOptionalInteger(kOrig);
1342                                }
1343                            }
1344                        } else if (rowUpper[i] > 1.0e20) {
1345                            // treat rhs as multiple of 1 unless elements all same
1346                            double value = ((possibleRow[i]==2) ? rowLower[i] : 1.0) * mult;
1347                            if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
1348                                del[nDel++] = i;
1349                                if (columnLength[iColumn] == 1) {
1350                                    copy2->setInteger(iColumn);
1351                                    int kOrig = original[iColumn];
1352                                    setOptionalInteger(kOrig);
1353                                }
1354                            }
1355                        } else {
1356                            // treat rhs as multiple of 1 unless elements all same
1357                            double value = ((possibleRow[i]==2) ? rowUpper[i] : 1.0) * mult;
1358                            if (rowLower[i] == rowUpper[i] &&
1359                                    fabs(value - floor(value + 0.5)) < 1.0e-8) {
1360                                del[nDel++] = i;
1361                                copy2->setInteger(iColumn);
1362                                int kOrig = original[iColumn];
1363                                setOptionalInteger(kOrig);
1364                            }
1365                        }
1366                    }
1367                } else {
1368                    // only if all singletons
1369                    bool possible = false;
1370                    if (rowLower[i] < -1.0e20) {
1371                        double value = rowUpper[i];
1372                        if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1373                            possible = true;
1374                    } else if (rowUpper[i] > 1.0e20) {
1375                        double value = rowLower[i];
1376                        if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1377                            possible = true;
1378                    } else {
1379                        double value = rowUpper[i];
1380                        if (rowLower[i] == rowUpper[i] &&
1381                                fabs(value - floor(value + 0.5)) < 1.0e-8)
1382                            possible = true;
1383                    }
1384                    if (possible) {
1385                        for (CoinBigIndex j = rowStart[i];
1386                                j < rowStart[i] + rowLength[i]; j++) {
1387                            int iColumn = column[j];
1388                            if (columnLength[iColumn] != 1 || fabs(element[j]) != 1.0) {
1389                                possible = false;
1390                                break;
1391                            }
1392                        }
1393                        if (possible) {
1394                            for (CoinBigIndex j = rowStart[i];
1395                                    j < rowStart[i] + rowLength[i]; j++) {
1396                                int iColumn = column[j];
1397                                if (!copy2->isInteger(iColumn)) {
1398                                    copy2->setInteger(iColumn);
1399                                    int kOrig = original[iColumn];
1400                                    setOptionalInteger(kOrig);
1401                                }
1402                            }
1403                            del[nDel++] = i;
1404                        }
1405                    }
1406                }
1407            }
1408        }
1409        if (nDel) {
1410            copy2->deleteRows(nDel, del);
1411            // pack down possible
1412            int n=0;
1413            for (int i=0;i<nDel;i++)
1414              possibleRow[del[i]]=-1;
1415            for (int i=0;i<numberRows;i++) {
1416              if (possibleRow[i]>=0) 
1417                possibleRow[n++]=possibleRow[i];
1418            }
1419        }
1420        if (nDel != numberRows) {
1421            nDel = 0;
1422            for (int i = 0; i < numberColumns; i++) {
1423                if (copy2->isInteger(i)) {
1424                    del[nDel++] = i;
1425                    nExtra++;
1426                }
1427            }
1428        } else {
1429            nDel = 0;
1430        }
1431        delete copy1;
1432        copy1 = copy2->clone();
1433        delete copy2;
1434    }
1435    // See if what's left is a network
1436    bool couldBeNetwork = false;
1437    if (copy1->getNumRows() && copy1->getNumCols()) {
1438#ifdef COIN_HAS_CLP
1439        OsiClpSolverInterface * clpSolver
1440        = dynamic_cast<OsiClpSolverInterface *> (copy1);
1441        if (false && clpSolver) {
1442            numberRows = clpSolver->getNumRows();
1443            char * rotate = new char[numberRows];
1444            int n = clpSolver->getModelPtr()->findNetwork(rotate, 1.0);
1445            delete [] rotate;
1446#ifdef CLP_INVESTIGATE
1447            printf("INTA network %d rows out of %d\n", n, numberRows);
1448#endif
1449            if (CoinAbs(n) == numberRows) {
1450                couldBeNetwork = true;
1451                for (int i = 0; i < numberRows; i++) {
1452                    if (!possibleRow[i]) {
1453                        couldBeNetwork = false;
1454#ifdef CLP_INVESTIGATE
1455                        printf("but row %d is bad\n", i);
1456#endif
1457                        break;
1458                    }
1459                }
1460            }
1461        } else
1462#endif
1463        {
1464            numberColumns = copy1->getNumCols();
1465            numberRows = copy1->getNumRows();
1466            const double * rowLower = copy1->getRowLower();
1467            const double * rowUpper = copy1->getRowUpper();
1468            couldBeNetwork = true;
1469            for (int i = 0; i < numberRows; i++) {
1470                if (rowLower[i] > -1.0e20 && fabs(rowLower[i] - floor(rowLower[i] + 0.5)) > 1.0e-12) {
1471                    couldBeNetwork = false;
1472                    break;
1473                }
1474                if (rowUpper[i] < 1.0e20 && fabs(rowUpper[i] - floor(rowUpper[i] + 0.5)) > 1.0e-12) {
1475                    couldBeNetwork = false;
1476                    break;
1477                }
1478                if (possibleRow[i]==0) {
1479                    couldBeNetwork = false;
1480                    break;
1481                }
1482            }
1483            if (couldBeNetwork) {
1484                const CoinPackedMatrix  * matrixByCol = copy1->getMatrixByCol();
1485                const double * element = matrixByCol->getElements();
1486                //const int * row = matrixByCol->getIndices();
1487                const CoinBigIndex * columnStart = matrixByCol->getVectorStarts();
1488                const int * columnLength = matrixByCol->getVectorLengths();
1489                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1490                    CoinBigIndex start = columnStart[iColumn];
1491                    CoinBigIndex end = start + columnLength[iColumn];
1492                    if (end > start + 2) {
1493                        couldBeNetwork = false;
1494                        break;
1495                    }
1496                    int type = 0;
1497                    for (CoinBigIndex j = start; j < end; j++) {
1498                        double value = element[j];
1499                        if (fabs(value) != 1.0) {
1500                            couldBeNetwork = false;
1501                            break;
1502                        } else if (value == 1.0) {
1503                            if ((type&1) == 0)
1504                                type |= 1;
1505                            else
1506                                type = 7;
1507                        } else if (value == -1.0) {
1508                            if ((type&2) == 0)
1509                                type |= 2;
1510                            else
1511                                type = 7;
1512                        }
1513                    }
1514                    if (type > 3) {
1515                        couldBeNetwork = false;
1516                        break;
1517                    }
1518                }
1519            }
1520        }
1521    }
1522    if (couldBeNetwork) {
1523        for (int i = 0; i < numberColumns; i++)
1524            setOptionalInteger(original[i]);
1525    }
1526    if (nExtra || couldBeNetwork) {
1527        numberColumns = copy1->getNumCols();
1528        numberRows = copy1->getNumRows();
1529        if (!numberColumns || !numberRows) {
1530            int numberColumns = solver_->getNumCols();
1531            for (int i = 0; i < numberColumns; i++)
1532                assert(solver_->isInteger(i));
1533        }
1534#ifdef CLP_INVESTIGATE
1535        if (couldBeNetwork || nExtra)
1536            printf("INTA %d extra integers, %d left%s\n", nExtra,
1537                   numberColumns,
1538                   couldBeNetwork ? ", all network" : "");
1539#endif
1540        findIntegers(true, 2);
1541        convertToDynamic();
1542    }
1543#ifdef CLP_INVESTIGATE
1544    if (!couldBeNetwork && copy1->getNumCols() &&
1545            copy1->getNumRows()) {
1546        printf("INTA %d rows and %d columns remain\n",
1547               copy1->getNumRows(), copy1->getNumCols());
1548        if (copy1->getNumCols() < 200) {
1549            copy1->writeMps("moreint");
1550            printf("INTA Written remainder to moreint.mps.gz %d rows %d cols\n",
1551                   copy1->getNumRows(), copy1->getNumCols());
1552        }
1553    }
1554#endif
1555    delete copy1;
1556    delete [] del;
1557    delete [] original;
1558    delete [] possibleRow;
1559    // double check increment
1560    analyzeObjective();
1561}
1562/**
1563  \todo
1564  Normally, it looks like we enter here from command dispatch in the main
1565  routine, after calling the solver for an initial solution
1566  (CbcModel::initialSolve, which simply calls the solver's initialSolve
1567  routine.) The first thing we do is call resolve. Presumably there are
1568  circumstances where this is nontrivial? There's also a call from
1569  CbcModel::originalModel (tied up with integer presolve), which should be
1570  checked.
1571
1572*/
1573
1574/*
1575  The overall flow can be divided into three stages:
1576    * Prep: Check that the lp relaxation remains feasible at the root. If so,
1577      do all the setup for B&C.
1578    * Process the root node: Generate cuts, apply heuristics, and in general do
1579      the best we can to resolve the problem without B&C.
1580    * Do B&C search until we hit a limit or exhaust the search tree.
1581
1582  Keep in mind that in general there is no node in the search tree that
1583  corresponds to the active subproblem. The active subproblem is represented
1584  by the current state of the model,  of the solver, and of the constraint
1585  system held by the solver.
1586*/
1587#ifdef COIN_HAS_CPX
1588#include "OsiCpxSolverInterface.hpp"
1589#include "cplex.h"
1590#endif
1591void CbcModel::branchAndBound(int doStatistics)
1592
1593{
1594  if (!parentModel_)
1595    /*
1596      Capture a time stamp before we start (unless set).
1597    */
1598    if (!dblParam_[CbcStartSeconds]) {
1599      if (!useElapsedTime())
1600        dblParam_[CbcStartSeconds] = CoinCpuTime();
1601      else
1602        dblParam_[CbcStartSeconds] = CoinGetTimeOfDay();
1603    }
1604    dblParam_[CbcSmallestChange] = COIN_DBL_MAX;
1605    dblParam_[CbcSumChange] = 0.0;
1606    dblParam_[CbcLargestChange] = 0.0;
1607    intParam_[CbcNumberBranches] = 0;
1608    // Force minimization !!!!
1609    bool flipObjective = (solver_->getObjSense()<0.0);
1610    if (flipObjective)
1611      flipModel();
1612    dblParam_[CbcOptimizationDirection] = 1.0; // was solver_->getObjSense();
1613    strongInfo_[0] = 0;
1614    strongInfo_[1] = 0;
1615    strongInfo_[2] = 0;
1616    strongInfo_[3] = 0;
1617    strongInfo_[4] = 0;
1618    strongInfo_[5] = 0;
1619    strongInfo_[6] = 0;
1620    numberStrongIterations_ = 0;
1621    currentNode_ = NULL;
1622    // See if should do cuts old way
1623    if (parallelMode() < 0) {
1624        specialOptions_ |= 4096 + 8192;
1625    } else if (parallelMode() > 0) {
1626        specialOptions_ |= 4096;
1627    }
1628    int saveMoreSpecialOptions = moreSpecialOptions_;
1629    if (dynamic_cast<CbcTreeLocal *> (tree_))
1630        specialOptions_ |= 4096 + 8192;
1631#ifdef COIN_HAS_CLP
1632    {
1633        OsiClpSolverInterface * clpSolver
1634        = dynamic_cast<OsiClpSolverInterface *> (solver_);
1635        if (clpSolver) {
1636            // pass in disaster handler
1637            CbcDisasterHandler handler(this);
1638            clpSolver->passInDisasterHandler(&handler);
1639            // Initialise solvers seed (unless users says not)
1640            if ((specialOptions_&4194304)==0)
1641              clpSolver->getModelPtr()->setRandomSeed(1234567);
1642#ifdef JJF_ZERO
1643            // reduce factorization frequency
1644            int frequency = clpSolver->getModelPtr()->factorizationFrequency();
1645            clpSolver->getModelPtr()->setFactorizationFrequency(CoinMin(frequency, 120));
1646#endif
1647        }
1648    }
1649#endif
1650    // original solver (only set if pre-processing)
1651    OsiSolverInterface * originalSolver = NULL;
1652    int numberOriginalObjects = numberObjects_;
1653    OsiObject ** originalObject = NULL;
1654    // Save whether there were any objects
1655    bool noObjects = (numberObjects_ == 0);
1656    // Set up strategies
1657    /*
1658      See if the user has supplied a strategy object and deal with it if present.
1659      The call to setupOther will set numberStrong_ and numberBeforeTrust_, and
1660      perform integer preprocessing, if requested.
1661
1662      We need to hang on to a pointer to solver_. setupOther will assign a
1663      preprocessed solver to model, but will instruct assignSolver not to trash the
1664      existing one.
1665    */
1666    if (strategy_) {
1667        // May do preprocessing
1668        originalSolver = solver_;
1669        strategy_->setupOther(*this);
1670        if (strategy_->preProcessState()) {
1671            // pre-processing done
1672            if (strategy_->preProcessState() < 0) {
1673                // infeasible
1674                handler_->message(CBC_INFEAS, messages_) << CoinMessageEol ;
1675                status_ = 0 ;
1676                secondaryStatus_ = 1;
1677                originalContinuousObjective_ = COIN_DBL_MAX;
1678                if (flipObjective)
1679                  flipModel();
1680                return ;
1681            } else if (numberObjects_ && object_) {
1682                numberOriginalObjects = numberObjects_;
1683                // redo sequence
1684                numberIntegers_ = 0;
1685                int numberColumns = getNumCols();
1686                int nOrig = originalSolver->getNumCols();
1687                CglPreProcess * process = strategy_->process();
1688                assert (process);
1689                const int * originalColumns = process->originalColumns();
1690                // allow for cliques etc
1691                nOrig = CoinMax(nOrig, originalColumns[numberColumns-1] + 1);
1692                // try and redo debugger
1693                OsiRowCutDebugger * debugger = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
1694                if (debugger) {
1695                  if (numberColumns<=debugger->numberColumns())
1696                    debugger->redoSolution(numberColumns, originalColumns);
1697                  else
1698                    debugger=NULL; // no idea how to handle (SOS?)
1699                }
1700                // User-provided solution might have been best. Synchronise.
1701                if (bestSolution_) {
1702                    // need to redo - in case no better found in BAB
1703                    // just get integer part right
1704                    for (int i = 0; i < numberColumns; i++) {
1705                        int jColumn = originalColumns[i];
1706                        bestSolution_[i] = bestSolution_[jColumn];
1707                    }
1708                }
1709                originalObject = object_;
1710                // object number or -1
1711                int * temp = new int[nOrig];
1712                int iColumn;
1713                for (iColumn = 0; iColumn < nOrig; iColumn++)
1714                    temp[iColumn] = -1;
1715                int iObject;
1716                int nNonInt = 0;
1717                for (iObject = 0; iObject < numberOriginalObjects; iObject++) {
1718                    iColumn = originalObject[iObject]->columnNumber();
1719                    if (iColumn < 0) {
1720                        nNonInt++;
1721                    } else {
1722                        temp[iColumn] = iObject;
1723                    }
1724                }
1725                int numberNewIntegers = 0;
1726                int numberOldIntegers = 0;
1727                int numberOldOther = 0;
1728                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1729                    int jColumn = originalColumns[iColumn];
1730                    if (temp[jColumn] >= 0) {
1731                        int iObject = temp[jColumn];
1732                        CbcSimpleInteger * obj =
1733                            dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ;
1734                        if (obj)
1735                            numberOldIntegers++;
1736                        else
1737                            numberOldOther++;
1738                    } else if (isInteger(iColumn)) {
1739                        numberNewIntegers++;
1740                    }
1741                }
1742                /*
1743                  Allocate an array to hold the indices of the integer variables.
1744                  Make a large enough array for all objects
1745                */
1746                numberObjects_ = numberNewIntegers + numberOldIntegers + numberOldOther + nNonInt;
1747                object_ = new OsiObject * [numberObjects_];
1748                delete [] integerVariable_;
1749                integerVariable_ = new int [numberNewIntegers+numberOldIntegers];
1750                /*
1751                  Walk the variables again, filling in the indices and creating objects for
1752                  the integer variables. Initially, the objects hold the index and upper &
1753                  lower bounds.
1754                */
1755                numberIntegers_ = 0;
1756                int n = originalColumns[numberColumns-1] + 1;
1757                int * backward = new int[n];
1758                int i;
1759                for ( i = 0; i < n; i++)
1760                    backward[i] = -1;
1761                for (i = 0; i < numberColumns; i++)
1762                    backward[originalColumns[i]] = i;
1763                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1764                    int jColumn = originalColumns[iColumn];
1765                    if (temp[jColumn] >= 0) {
1766                        int iObject = temp[jColumn];
1767                        CbcSimpleInteger * obj =
1768                            dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ;
1769                        if (obj) {
1770                            object_[numberIntegers_] = originalObject[iObject]->clone();
1771                            // redo ids etc
1772                            //object_[numberIntegers_]->resetSequenceEtc(numberColumns,originalColumns);
1773                            object_[numberIntegers_]->resetSequenceEtc(numberColumns, backward);
1774                            integerVariable_[numberIntegers_++] = iColumn;
1775                        }
1776                    } else if (isInteger(iColumn)) {
1777                        object_[numberIntegers_] =
1778                            new CbcSimpleInteger(this, iColumn);
1779                        integerVariable_[numberIntegers_++] = iColumn;
1780                    }
1781                }
1782                delete [] backward;
1783                numberObjects_ = numberIntegers_;
1784                // Now append other column stuff
1785                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1786                    int jColumn = originalColumns[iColumn];
1787                    if (temp[jColumn] >= 0) {
1788                        int iObject = temp[jColumn];
1789                        CbcSimpleInteger * obj =
1790                            dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ;
1791                        if (!obj) {
1792                            object_[numberObjects_] = originalObject[iObject]->clone();
1793                            // redo ids etc
1794                            CbcObject * obj =
1795                                dynamic_cast <CbcObject *>(object_[numberObjects_]) ;
1796                            assert (obj);
1797                            obj->redoSequenceEtc(this, numberColumns, originalColumns);
1798                            numberObjects_++;
1799                        }
1800                    }
1801                }
1802                // now append non column stuff
1803                for (iObject = 0; iObject < numberOriginalObjects; iObject++) {
1804                    iColumn = originalObject[iObject]->columnNumber();
1805                    if (iColumn < 0) {
1806                        // already has column numbers changed
1807                        object_[numberObjects_] = originalObject[iObject]->clone();
1808#ifdef JJF_ZERO
1809                        // redo ids etc
1810                        CbcObject * obj =
1811                            dynamic_cast <CbcObject *>(object_[numberObjects_]) ;
1812                        assert (obj);
1813                        obj->redoSequenceEtc(this, numberColumns, originalColumns);
1814#endif
1815                        numberObjects_++;
1816                    }
1817                }
1818                delete [] temp;
1819                if (!numberObjects_)
1820                    handler_->message(CBC_NOINT, messages_) << CoinMessageEol ;
1821            } else {
1822                int numberColumns = getNumCols();
1823                CglPreProcess * process = strategy_->process();
1824                assert (process);
1825                const int * originalColumns = process->originalColumns();
1826                // try and redo debugger
1827                OsiRowCutDebugger * debugger = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
1828                if (debugger)
1829                    debugger->redoSolution(numberColumns, originalColumns);
1830            }
1831        } else {
1832            //no preprocessing
1833            originalSolver = NULL;
1834        }
1835        strategy_->setupCutGenerators(*this);
1836        strategy_->setupHeuristics(*this);
1837        // Set strategy print level to models
1838        strategy_->setupPrinting(*this, handler_->logLevel());
1839    }
1840    eventHappened_ = false;
1841    CbcEventHandler *eventHandler = getEventHandler() ;
1842    if (eventHandler)
1843        eventHandler->setModel(this);
1844#define CLIQUE_ANALYSIS
1845#ifdef CLIQUE_ANALYSIS
1846    // set up for probing
1847    // If we're doing clever stuff with cliques, additional info here.
1848    if (!parentModel_)
1849        probingInfo_ = new CglTreeProbingInfo(solver_);
1850    else
1851        probingInfo_ = NULL;
1852#else
1853    probingInfo_ = NULL;
1854#endif
1855
1856    // Try for dominated columns
1857    if ((specialOptions_&64) != 0) {
1858        CglDuplicateRow dupcuts(solver_);
1859        dupcuts.setMode(2);
1860        CglStored * storedCuts = dupcuts.outDuplicates(solver_);
1861        if (storedCuts) {
1862          COIN_DETAIL_PRINT(printf("adding dup cuts\n"));
1863            addCutGenerator(storedCuts, 1, "StoredCuts from dominated",
1864                            true, false, false, -200);
1865        }
1866    }
1867    if (!nodeCompare_)
1868        nodeCompare_ = new CbcCompareDefault();;
1869    // See if hot start wanted
1870    CbcCompareBase * saveCompare = NULL;
1871    // User supplied hotstart. Adapt for preprocessing.
1872    if (hotstartSolution_) {
1873        if (strategy_ && strategy_->preProcessState() > 0) {
1874            CglPreProcess * process = strategy_->process();
1875            assert (process);
1876            int n = solver_->getNumCols();
1877            const int * originalColumns = process->originalColumns();
1878            // columns should be in order ... but
1879            double * tempS = new double[n];
1880            for (int i = 0; i < n; i++) {
1881                int iColumn = originalColumns[i];
1882                tempS[i] = hotstartSolution_[iColumn];
1883            }
1884            delete [] hotstartSolution_;
1885            hotstartSolution_ = tempS;
1886            if (hotstartPriorities_) {
1887                int * tempP = new int [n];
1888                for (int i = 0; i < n; i++) {
1889                    int iColumn = originalColumns[i];
1890                    tempP[i] = hotstartPriorities_[iColumn];
1891                }
1892                delete [] hotstartPriorities_;
1893                hotstartPriorities_ = tempP;
1894            }
1895        }
1896        saveCompare = nodeCompare_;
1897        // depth first
1898        nodeCompare_ = new CbcCompareDepth();
1899    }
1900    if (!problemFeasibility_)
1901        problemFeasibility_ = new CbcFeasibilityBase();
1902# ifdef CBC_DEBUG
1903    std::string problemName ;
1904    solver_->getStrParam(OsiProbName, problemName) ;
1905    printf("Problem name - %s\n", problemName.c_str()) ;
1906    solver_->setHintParam(OsiDoReducePrint, false, OsiHintDo, 0) ;
1907# endif
1908    /*
1909      Assume we're done, and see if we're proven wrong.
1910    */
1911    status_ = 0 ;
1912    secondaryStatus_ = 0;
1913    phase_ = 0;
1914    /*
1915      Scan the variables, noting the integer variables. Create an
1916      CbcSimpleInteger object for each integer variable.
1917    */
1918    findIntegers(false) ;
1919    // Say not dynamic pseudo costs
1920    ownership_ &= ~0x40000000;
1921    // If dynamic pseudo costs then do
1922    if (numberBeforeTrust_)
1923        convertToDynamic();
1924    // Set up char array to say if integer (speed)
1925    delete [] integerInfo_;
1926    {
1927        int n = solver_->getNumCols();
1928        integerInfo_ = new char [n];
1929        for (int i = 0; i < n; i++) {
1930            if (solver_->isInteger(i))
1931                integerInfo_[i] = 1;
1932            else
1933                integerInfo_[i] = 0;
1934        }
1935    }
1936    if (preferredWay_) {
1937        // set all unset ones
1938        for (int iObject = 0 ; iObject < numberObjects_ ; iObject++) {
1939            CbcObject * obj =
1940                dynamic_cast <CbcObject *>(object_[iObject]) ;
1941            if (obj && !obj->preferredWay())
1942                obj->setPreferredWay(preferredWay_);
1943        }
1944    }
1945    /*
1946      Ensure that objects on the lists of OsiObjects, heuristics, and cut
1947      generators attached to this model all refer to this model.
1948    */
1949    synchronizeModel() ;
1950    if (!solverCharacteristics_) {
1951        OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
1952        if (solverCharacteristics) {
1953            solverCharacteristics_ = solverCharacteristics;
1954        } else {
1955            // replace in solver
1956            OsiBabSolver defaultC;
1957            solver_->setAuxiliaryInfo(&defaultC);
1958            solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
1959        }
1960    }
1961
1962    solverCharacteristics_->setSolver(solver_);
1963    // Set so we can tell we are in initial phase in resolve
1964    continuousObjective_ = -COIN_DBL_MAX ;
1965    /*
1966      Solve the relaxation.
1967
1968      Apparently there are circumstances where this will be non-trivial --- i.e.,
1969      we've done something since initialSolve that's trashed the solution to the
1970      continuous relaxation.
1971    */
1972    /* Tell solver we are in Branch and Cut
1973       Could use last parameter for subtle differences */
1974    solver_->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
1975#ifdef COIN_HAS_CLP
1976    {
1977        OsiClpSolverInterface * clpSolver
1978        = dynamic_cast<OsiClpSolverInterface *> (solver_);
1979        if (clpSolver) {
1980            ClpSimplex * clpSimplex = clpSolver->getModelPtr();
1981            if ((specialOptions_&32) == 0) {
1982                // take off names
1983                clpSimplex->dropNames();
1984            }
1985            // no crunch if mostly continuous
1986            if ((clpSolver->specialOptions()&(1 + 8)) != (1 + 8)) {
1987                int numberColumns = solver_->getNumCols();
1988                if (numberColumns > 1000 && numberIntegers_*4 < numberColumns)
1989                    clpSolver->setSpecialOptions(clpSolver->specialOptions()&(~1));
1990            }
1991            //#define NO_CRUNCH
1992#ifdef NO_CRUNCH
1993            printf("TEMP switching off crunch\n");
1994            int iOpt = clpSolver->specialOptions();
1995            iOpt &= ~1;
1996            iOpt |= 65536;
1997            clpSolver->setSpecialOptions(iOpt);
1998#endif
1999        }
2000    }
2001#endif
2002    bool feasible;
2003    numberSolves_ = 0 ;
2004    // If NLP then we assume already solved outside branchAndbound
2005    if (!solverCharacteristics_->solverType() || solverCharacteristics_->solverType() == 4) {
2006        feasible = resolve(NULL, 0) != 0 ;
2007    } else {
2008        // pick up given status
2009        feasible = (solver_->isProvenOptimal() &&
2010                    !solver_->isDualObjectiveLimitReached()) ;
2011    }
2012    if (problemFeasibility_->feasible(this, 0) < 0) {
2013        feasible = false; // pretend infeasible
2014    }
2015    numberSavedSolutions_ = 0;
2016    int saveNumberStrong = numberStrong_;
2017    int saveNumberBeforeTrust = numberBeforeTrust_;
2018    /*
2019      If the linear relaxation of the root is infeasible, bail out now. Otherwise,
2020      continue with processing the root node.
2021    */
2022    if (!feasible) {
2023        status_ = 0 ;
2024        if (!solver_->isProvenDualInfeasible()) {
2025            handler_->message(CBC_INFEAS, messages_) << CoinMessageEol ;
2026            secondaryStatus_ = 1;
2027        } else {
2028            handler_->message(CBC_UNBOUNDED, messages_) << CoinMessageEol ;
2029            secondaryStatus_ = 7;
2030        }
2031        originalContinuousObjective_ = COIN_DBL_MAX;
2032        solverCharacteristics_ = NULL;
2033        if (flipObjective)
2034          flipModel();
2035        return ;
2036    } else if (!numberObjects_ && (!strategy_ || strategy_->preProcessState() <= 0)) {
2037        // nothing to do
2038        if (flipObjective)
2039          flipModel();
2040        solverCharacteristics_ = NULL;
2041        bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
2042        int numberColumns = solver_->getNumCols();
2043        delete [] bestSolution_;
2044        bestSolution_ = new double[numberColumns];
2045        CoinCopyN(solver_->getColSolution(), numberColumns, bestSolution_);
2046        return ;
2047    }
2048    /*
2049      See if we're using the Osi side of the branching hierarchy. If so, either
2050      convert existing CbcObjects to OsiObjects, or generate them fresh. In the
2051      first case, CbcModel owns the objects on the object_ list. In the second
2052      case, the solver holds the objects and object_ simply points to the
2053      solver's list.
2054
2055      080417 The conversion code here (the block protected by `if (obj)') cannot
2056      possibly be correct. On the Osi side, descent is OsiObject -> OsiObject2 ->
2057      all other Osi object classes. On the Cbc side, it's OsiObject -> CbcObject
2058      -> all other Cbc object classes. It's structurally impossible for any Osi
2059      object to descend from CbcObject. The only thing I can see is that this is
2060      really dead code, and object detection is now handled from the Osi side.
2061    */
2062    // Convert to Osi if wanted
2063    bool useOsiBranching = false;
2064    //OsiBranchingInformation * persistentInfo = NULL;
2065    if (branchingMethod_ && branchingMethod_->chooseMethod()) {
2066        useOsiBranching = true;
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    double bestValue = 0.0 ;
3875    if (newNode) {
3876        bestValue = newNode->objectiveValue();
3877        if (newNode->branchingObject()) {
3878            newNode->initializeInfo() ;
3879            tree_->push(newNode) ;
3880            // save pointer to root node - so can pick up bounds
3881            if (!topOfTree_)
3882              topOfTree_ = dynamic_cast<CbcFullNodeInfo *>(newNode->nodeInfo()) ;
3883            if (statistics_) {
3884                if (numberNodes2_ == maximumStatistics_) {
3885                    maximumStatistics_ = 2 * maximumStatistics_;
3886                    CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_];
3887                    memset(temp, 0, maximumStatistics_*sizeof(CbcStatistics *));
3888                    memcpy(temp, statistics_, numberNodes2_*sizeof(CbcStatistics *));
3889                    delete [] statistics_;
3890                    statistics_ = temp;
3891                }
3892                assert (!statistics_[numberNodes2_]);
3893                statistics_[numberNodes2_] = new CbcStatistics(newNode, this);
3894            }
3895            numberNodes2_++;
3896#     ifdef CHECK_NODE
3897            printf("Node %x on tree\n", newNode) ;
3898#     endif
3899        } else {
3900            // continuous is integer
3901            double objectiveValue = newNode->objectiveValue();
3902            setBestSolution(CBC_SOLUTION, objectiveValue,
3903                            solver_->getColSolution()) ;
3904            delete newNode ;
3905            newNode = NULL ;
3906        }
3907    }
3908
3909    if (printFrequency_ <= 0) {
3910        printFrequency_ = 1000 ;
3911        if (getNumCols() > 2000)
3912            printFrequency_ = 100 ;
3913    }
3914    /*
3915      It is possible that strong branching fixes one variable and then the code goes round
3916      again and again.  This can take too long.  So we need to warn user - just once.
3917    */
3918    numberLongStrong_ = 0;
3919    CbcNode * createdNode = NULL;
3920#ifdef CBC_THREAD
3921    if ((specialOptions_&2048) != 0)
3922        numberThreads_ = 0;
3923    if (numberThreads_ ) {
3924        nodeCompare_->sayThreaded(); // need to use addresses
3925        master_ = new CbcBaseModel(*this,
3926                                   (parallelMode() < -1) ? 1 : 0);
3927        masterThread_ = master_->masterThread();
3928    }
3929#endif
3930#ifdef COIN_HAS_CLP
3931    {
3932        OsiClpSolverInterface * clpSolver
3933        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3934        if (clpSolver && !parentModel_) {
3935            clpSolver->computeLargestAway();
3936        }
3937    }
3938#endif
3939    /*
3940      At last, the actual branch-and-cut search loop, which will iterate until
3941      the live set is empty or we hit some limit (integrality gap, time, node
3942      count, etc.). The overall flow is to rebuild a subproblem, reoptimise using
3943      solveWithCuts(), choose a branching pattern with chooseBranch(), and finally
3944      add the node to the live set.
3945
3946      The first action is to winnow the live set to remove nodes which are worse
3947      than the current objective cutoff.
3948    */
3949    if (solver_->getRowCutDebuggerAlways()) {
3950        OsiRowCutDebugger * debuggerX = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
3951        const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger() ;
3952        if (!debugger) {
3953            // infeasible!!
3954            printf("before search\n");
3955            const double * lower = solver_->getColLower();
3956            const double * upper = solver_->getColUpper();
3957            const double * solution = debuggerX->optimalSolution();
3958            int numberColumns = solver_->getNumCols();
3959            for (int i = 0; i < numberColumns; i++) {
3960                if (solver_->isInteger(i)) {
3961                    if (solution[i] < lower[i] - 1.0e-6 || solution[i] > upper[i] + 1.0e-6)
3962                        printf("**** ");
3963                    printf("%d %g <= %g <= %g\n",
3964                           i, lower[i], solution[i], upper[i]);
3965                }
3966            }
3967            //abort();
3968        }
3969    }
3970    {
3971        // may be able to change cutoff now
3972        double cutoff = getCutoff();
3973        double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
3974        if (cutoff > bestObjective_ - increment) {
3975            cutoff = bestObjective_ - increment ;
3976            setCutoff(cutoff) ;
3977        }
3978    }
3979#ifdef CBC_THREAD
3980    bool goneParallel = false;
3981#endif
3982#define MAX_DEL_NODE 1
3983    CbcNode * delNode[MAX_DEL_NODE+1];
3984    int nDeleteNode = 0;
3985    // For Printing etc when parallel
3986    int lastEvery1000 = 0;
3987    int lastPrintEvery = 0;
3988    int numberConsecutiveInfeasible = 0;
3989#define PERTURB_IN_FATHOM
3990#ifdef PERTURB_IN_FATHOM
3991    // allow in fathom
3992    if ((moreSpecialOptions_& 262144) != 0)
3993      specialOptions_ |= 131072;
3994#endif
3995    while (true) {
3996        lockThread();
3997#ifdef COIN_HAS_CLP
3998        // See if we want dantzig row choice
3999        goToDantzig(100, savePivotMethod);
4000#endif
4001        if (tree_->empty()) {
4002#ifdef CBC_THREAD
4003            if (parallelMode() > 0 && master_) {
4004                int anyLeft = master_->waitForThreadsInTree(0);
4005                if (!anyLeft) {
4006                    master_->stopThreads(-1);
4007                    break;
4008                }
4009            } else {
4010                break;
4011            }
4012#else
4013            break;
4014#endif
4015        } else {
4016            unlockThread();
4017        }
4018        // If done 100 nodes see if worth trying reduction
4019        if (numberNodes_ == 50 || numberNodes_ == 100) {
4020#ifdef COIN_HAS_CLP
4021            OsiClpSolverInterface * clpSolver
4022            = dynamic_cast<OsiClpSolverInterface *> (solver_);
4023            if (clpSolver && ((specialOptions_&131072) == 0) && true) {
4024                ClpSimplex * simplex = clpSolver->getModelPtr();
4025                int perturbation = simplex->perturbation();
4026#ifdef CLP_INVESTIGATE
4027                printf("Testing its n,s %d %d solves n,s %d %d - pert %d\n",
4028                       numberIterations_, saveNumberIterations,
4029                       numberSolves_, saveNumberSolves, perturbation);
4030#endif
4031                if (perturbation == 50 && (numberIterations_ - saveNumberIterations) <
4032                        8*(numberSolves_ - saveNumberSolves)) {
4033                    // switch off perturbation
4034                    simplex->setPerturbation(100);
4035#ifdef CLP_INVESTIGATE
4036                    printf("Perturbation switched off\n");
4037#endif
4038                }
4039            }
4040#endif
4041            /*
4042              Decide if we want to do a restart.
4043            */
4044            if (saveSolver && (specialOptions_&(512 + 32768)) != 0) {
4045                bool tryNewSearch = solverCharacteristics_->reducedCostsAccurate() &&
4046                                    (getCutoff() < 1.0e20 && getCutoff() < checkCutoffForRestart);
4047                int numberColumns = getNumCols();
4048                if (tryNewSearch) {
4049                    checkCutoffForRestart = getCutoff() ;
4050#ifdef CLP_INVESTIGATE
4051                    printf("after %d nodes, cutoff %g - looking\n",
4052                           numberNodes_, getCutoff());
4053#endif
4054                    saveSolver->resolve();
4055                    double direction = saveSolver->getObjSense() ;
4056                    double gap = checkCutoffForRestart - saveSolver->getObjValue() * direction ;
4057                    double tolerance;
4058                    saveSolver->getDblParam(OsiDualTolerance, tolerance) ;
4059                    if (gap <= 0.0)
4060                        gap = tolerance;
4061                    gap += 100.0 * tolerance;
4062                    double integerTolerance = getDblParam(CbcIntegerTolerance) ;
4063
4064                    const double *lower = saveSolver->getColLower() ;
4065                    const double *upper = saveSolver->getColUpper() ;
4066                    const double *solution = saveSolver->getColSolution() ;
4067                    const double *reducedCost = saveSolver->getReducedCost() ;
4068
4069                    int numberFixed = 0 ;
4070                    int numberFixed2 = 0;
4071#ifdef COIN_DEVELOP
4072                    printf("gap %g\n", gap);
4073#endif
4074                    for (int i = 0 ; i < numberIntegers_ ; i++) {
4075                        int iColumn = integerVariable_[i] ;
4076                        double djValue = direction * reducedCost[iColumn] ;
4077                        if (upper[iColumn] - lower[iColumn] > integerTolerance) {
4078                            if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) {
4079                                //printf("%d to lb on dj of %g - bounds %g %g\n",
4080                                //     iColumn,djValue,lower[iColumn],upper[iColumn]);
4081                                saveSolver->setColUpper(iColumn, lower[iColumn]) ;
4082                                numberFixed++ ;
4083                            } else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) {
4084                                //printf("%d to ub on dj of %g - bounds %g %g\n",
4085                                //     iColumn,djValue,lower[iColumn],upper[iColumn]);
4086                                saveSolver->setColLower(iColumn, upper[iColumn]) ;
4087                                numberFixed++ ;
4088                            }
4089                        } else {
4090                            //printf("%d has dj of %g - already fixed to %g\n",
4091                            //     iColumn,djValue,lower[iColumn]);
4092                            numberFixed2++;
4093                        }
4094                    }
4095#ifdef COIN_DEVELOP
4096                    if ((specialOptions_&1) != 0) {
4097                        const OsiRowCutDebugger *debugger = saveSolver->getRowCutDebugger() ;
4098                        if (debugger) {
4099                            printf("Contains optimal\n") ;
4100                            OsiSolverInterface * temp = saveSolver->clone();
4101                            const double * solution = debugger->optimalSolution();
4102                            const double *lower = temp->getColLower() ;
4103                            const double *upper = temp->getColUpper() ;
4104                            int n = temp->getNumCols();
4105                            for (int i = 0; i < n; i++) {
4106                                if (temp->isInteger(i)) {
4107                                    double value = floor(solution[i] + 0.5);
4108                                    assert (value >= lower[i] && value <= upper[i]);
4109                                    temp->setColLower(i, value);
4110                                    temp->setColUpper(i, value);
4111                                }
4112                            }
4113                            temp->writeMps("reduced_fix");
4114                            delete temp;
4115                            saveSolver->writeMps("reduced");
4116                        } else {
4117                            abort();
4118                        }
4119                    }
4120                    printf("Restart could fix %d integers (%d already fixed)\n",
4121                           numberFixed + numberFixed2, numberFixed2);
4122#endif
4123                    numberFixed += numberFixed2;
4124                    if (numberFixed*10 < numberColumns && numberFixed*4 < numberIntegers_)
4125                        tryNewSearch = false;
4126                }
4127                if (tryNewSearch) {
4128                    // back to solver without cuts?
4129                    OsiSolverInterface * solver2 = saveSolver->clone();
4130                    const double *lower = saveSolver->getColLower() ;
4131                    const double *upper = saveSolver->getColUpper() ;
4132                    for (int i = 0 ; i < numberIntegers_ ; i++) {
4133                        int iColumn = integerVariable_[i] ;
4134                        solver2->setColLower(iColumn, lower[iColumn]);
4135                        solver2->setColUpper(iColumn, upper[iColumn]);
4136                    }
4137                    // swap
4138                    delete saveSolver;
4139                    saveSolver = solver2;
4140                    double * newSolution = new double[numberColumns];
4141                    double objectiveValue = checkCutoffForRestart;
4142                    // Save the best solution so far.
4143                    CbcSerendipity heuristic(*this);
4144                    if (bestSolution_)
4145                        heuristic.setInputSolution(bestSolution_, bestObjective_);
4146                    // Magic number
4147                    heuristic.setFractionSmall(0.8);
4148                    // `pumpTune' to stand-alone solver for explanations.
4149                    heuristic.setFeasibilityPumpOptions(1008013);
4150                    // Use numberNodes to say how many are original rows
4151                    heuristic.setNumberNodes(continuousSolver_->getNumRows());
4152#ifdef COIN_DEVELOP
4153                    if (continuousSolver_->getNumRows() <
4154                            solver_->getNumRows())
4155                        printf("%d rows added ZZZZZ\n",
4156                               solver_->getNumRows() - continuousSolver_->getNumRows());
4157#endif
4158                    int returnCode = heuristic.smallBranchAndBound(saveSolver,
4159                                     -1, newSolution,
4160                                     objectiveValue,
4161                                     checkCutoffForRestart, "Reduce");
4162                    if (returnCode < 0) {
4163#ifdef COIN_DEVELOP
4164                        printf("Restart - not small enough to do search after fixing\n");
4165#endif
4166                        delete [] newSolution;
4167                    } else {
4168                        // 1 for sol'n, 2 for finished, 3 for both
4169                        if ((returnCode&1) != 0) {
4170                            // increment number of solutions so other heuristics can test
4171                            numberSolutions_++;
4172                            numberHeuristicSolutions_++;
4173                            lastHeuristic_ = NULL;
4174                            setBestSolution(CBC_ROUNDING, objectiveValue, newSolution) ;
4175                        }
4176                        delete [] newSolution;
4177#ifdef CBC_THREAD
4178                        if (master_) {
4179                            lockThread();
4180                            if (parallelMode() > 0) {
4181                                while (master_->waitForThreadsInTree(0)) {
4182                                    lockThread();
4183                                    double dummyBest;
4184                                    tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4185                                    //unlockThread();
4186                                }
4187                            } else {
4188                                double dummyBest;
4189                                tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4190                            }
4191                            master_->waitForThreadsInTree(2);
4192                            delete master_;
4193                            master_ = NULL;
4194                            masterThread_ = NULL;
4195                        }
4196#endif
4197                        if (tree_->size()) {
4198                            double dummyBest;
4199                            tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4200                        }
4201                        break;
4202                    }
4203                }
4204                delete saveSolver;
4205                saveSolver = NULL;
4206            }
4207        }
4208        /*
4209          Check for abort on limits: node count, solution count, time, integrality gap.
4210        */
4211        if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
4212                numberSolutions_ < intParam_[CbcMaxNumSol] &&
4213                !maximumSecondsReached() &&
4214                !stoppedOnGap_ && !eventHappened_ && (maximumNumberIterations_ < 0 ||
4215                                                      numberIterations_ < maximumNumberIterations_))) {
4216            // out of loop
4217            break;
4218        }
4219#ifdef BONMIN
4220        assert(!solverCharacteristics_->solutionAddsCuts() || solverCharacteristics_->mipFeasible());
4221#endif
4222// Sets percentage of time when we try diving. Diving requires a bit of heap reorganisation, because
4223// we need to replace the comparison function to dive, and that requires reordering to retain the
4224// heap property.
4225#define DIVE_WHEN 1000
4226#define DIVE_STOP 2000
4227        int kNode = numberNodes_ % 4000;
4228        if (numberNodes_<100000 && kNode>DIVE_WHEN && kNode <= DIVE_STOP) {
4229            if (!parallelMode()) {
4230                if (kNode == DIVE_WHEN + 1 || numberConsecutiveInfeasible > 1) {
4231                    CbcCompareDefault * compare = dynamic_cast<CbcCompareDefault *>
4232                                                  (nodeCompare_);
4233                    // Don't interfere if user has replaced the compare function.
4234                    if (compare) {
4235                        //printf("Redoing tree\n");
4236                        compare->startDive(this);
4237                        numberConsecutiveInfeasible = 0;
4238                    }
4239                }
4240            }
4241        }
4242        // replace current cutoff?
4243        if (cutoff > getCutoff()) {
4244            double newCutoff = getCutoff();
4245            if (analyzeResults_) {
4246                // see if we could fix any (more)
4247                int n = 0;
4248                double * newLower = analyzeResults_;
4249                double * objLower = newLower + numberIntegers_;
4250                double * newUpper = objLower + numberIntegers_;
4251                double * objUpper = newUpper + numberIntegers_;
4252                for (int i = 0; i < numberIntegers_; i++) {
4253                    if (objLower[i] > newCutoff) {
4254                        n++;
4255                        if (objUpper[i] > newCutoff) {
4256                            newCutoff = -COIN_DBL_MAX;
4257                            break;
4258                        }
4259                    } else if (objUpper[i] > newCutoff) {
4260                        n++;
4261                    }
4262                }
4263                if (newCutoff == -COIN_DBL_MAX) {
4264                  COIN_DETAIL_PRINT(printf("Root analysis says finished\n"));
4265                } else if (n > numberFixedNow_) {
4266                  COIN_DETAIL_PRINT(printf("%d more fixed by analysis - now %d\n", n - numberFixedNow_, n));
4267                    numberFixedNow_ = n;
4268                }
4269            }
4270            if (eventHandler) {
4271                if (!eventHandler->event(CbcEventHandler::solution)) {
4272                    eventHappened_ = true; // exit
4273                }
4274                newCutoff = getCutoff();
4275            }
4276            lockThread();
4277            /*
4278              Clean the tree to reflect the new solution, then see if the
4279              node comparison predicate wants to make any changes. If so,
4280              call setComparison for the side effect of rebuilding the heap.
4281            */
4282            tree_->cleanTree(this,newCutoff,bestPossibleObjective_) ;
4283            if (nodeCompare_->newSolution(this) ||
4284                nodeCompare_->newSolution(this,continuousObjective_,
4285                                          continuousInfeasibilities_)) {
4286              tree_->setComparison(*nodeCompare_) ;
4287            }
4288            if (tree_->empty()) {
4289                continue;
4290            }
4291            unlockThread();
4292        }
4293        cutoff = getCutoff() ;
4294        /*
4295            Periodic activities: Opportunities to
4296            + tweak the nodeCompare criteria,
4297            + check if we've closed the integrality gap enough to quit,
4298            + print a summary line to let the user know we're working
4299        */
4300        if (numberNodes_ >= lastEvery1000) {
4301            lockThread();
4302#ifdef COIN_HAS_CLP
4303            // See if we want dantzig row choice
4304            goToDantzig(1000, savePivotMethod);
4305#endif
4306            lastEvery1000 = numberNodes_ + 1000;
4307            bool redoTree = nodeCompare_->every1000Nodes(this, numberNodes_) ;
4308#ifdef CHECK_CUT_SIZE
4309            verifyCutSize (tree_, *this);
4310#endif
4311            // redo tree if requested
4312            if (redoTree)
4313                tree_->setComparison(*nodeCompare_) ;
4314            unlockThread();
4315        }
4316        // Had hotstart before, now switched off
4317        if (saveCompare && !hotstartSolution_) {
4318            // hotstart switched off
4319            delete nodeCompare_; // off depth first
4320            nodeCompare_ = saveCompare;
4321            saveCompare = NULL;
4322            // redo tree
4323            lockThread();
4324            tree_->setComparison(*nodeCompare_) ;
4325            unlockThread();
4326        }
4327        if (numberNodes_ >= lastPrintEvery) {
4328            lastPrintEvery = numberNodes_ + printFrequency_;
4329            lockThread();
4330            int nNodes = tree_->size() ;
4331
4332            //MODIF PIERRE
4333            bestPossibleObjective_ = tree_->getBestPossibleObjective();
4334            unlockThread();
4335#ifdef CLP_INVESTIGATE
4336            if (getCutoff() < 1.0e20) {
4337                if (fabs(getCutoff() - (bestObjective_ - getCutoffIncrement())) > 1.0e-6 &&
4338                        !parentModel_)
4339                    printf("model cutoff in status %g, best %g, increment %g\n",
4340                           getCutoff(), bestObjective_, getCutoffIncrement());
4341                assert (getCutoff() < bestObjective_ - getCutoffIncrement() +
4342                        1.0e-6 + 1.0e-10*fabs(bestObjective_));
4343            }
4344#endif
4345            if (!intParam_[CbcPrinting]) {
4346                messageHandler()->message(CBC_STATUS, messages())
4347                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4348                << getCurrentSeconds()
4349                << CoinMessageEol ;
4350            } else if (intParam_[CbcPrinting] == 1) {
4351                messageHandler()->message(CBC_STATUS2, messages())
4352                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4353                << lastDepth << lastUnsatisfied << numberIterations_
4354                << getCurrentSeconds()
4355                << CoinMessageEol ;
4356            } else if (!numberExtraIterations_) {
4357                messageHandler()->message(CBC_STATUS2, messages())
4358                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4359                << lastDepth << lastUnsatisfied << numberIterations_
4360                << getCurrentSeconds()
4361                << CoinMessageEol ;
4362            } else {
4363                messageHandler()->message(CBC_STATUS3, messages())
4364                << numberNodes_ << numberExtraNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
4365                << lastDepth << lastUnsatisfied << numberIterations_ << numberExtraIterations_
4366                << getCurrentSeconds()
4367                << CoinMessageEol ;
4368            }
4369            if (eventHandler && !eventHandler->event(CbcEventHandler::treeStatus)) {
4370                eventHappened_ = true; // exit
4371            }
4372        }
4373        // See if can stop on gap
4374        if(canStopOnGap()) {
4375            stoppedOnGap_ = true ;
4376        }
4377
4378#ifdef CHECK_NODE_FULL
4379        verifyTreeNodes(tree_, *this) ;
4380#   endif
4381#   ifdef CHECK_CUT_COUNTS
4382        verifyCutCounts(tree_, *this) ;
4383#   endif
4384        /*
4385          Now we come to the meat of the loop. To create the active subproblem, we'll
4386          pop the most promising node in the live set, rebuild the subproblem it
4387          represents, and then execute the current arm of the branch to create the
4388          active subproblem.
4389        */
4390        CbcNode * node = NULL;
4391#ifdef CBC_THREAD
4392        if (!parallelMode() || parallelMode() == -1) {
4393#endif
4394            node = tree_->bestNode(cutoff) ;
4395            // Possible one on tree worse than cutoff
4396            // Weird comparison function can leave ineligible nodes on tree
4397            if (!node || node->objectiveValue() > cutoff)
4398                continue;
4399            // Do main work of solving node here
4400            doOneNode(this, node, createdNode);
4401#ifdef JJF_ZERO
4402            if (node) {
4403              if (createdNode) {
4404                printf("Node %d depth %d, created %d depth %d\n",
4405                       node->nodeNumber(), node->depth(),
4406                       createdNode->nodeNumber(), createdNode->depth());
4407              } else {
4408                printf("Node %d depth %d,  no created node\n",
4409                       node->nodeNumber(), node->depth());
4410              }
4411            } else if (createdNode) {
4412              printf("Node exhausted, created %d depth %d\n",
4413                     createdNode->nodeNumber(), createdNode->depth());
4414            } else {
4415              printf("Node exhausted,  no created node\n");
4416              numberConsecutiveInfeasible = 2;
4417            }
4418#endif
4419            //if (createdNode)
4420            //numberConsecutiveInfeasible=0;
4421            //else
4422            //numberConsecutiveInfeasible++;
4423#ifdef CBC_THREAD
4424        } else if (parallelMode() > 0) {
4425            //lockThread();
4426            //node = tree_->bestNode(cutoff) ;
4427            // Possible one on tree worse than cutoff
4428            if (true || !node || node->objectiveValue() > cutoff) {
4429                assert (master_);
4430                if (master_) {
4431                    int anyLeft = master_->waitForThreadsInTree(1);
4432                    // may need to go round again
4433                    if (anyLeft) {
4434                        continue;
4435                    } else {
4436                        master_->stopThreads(-1);
4437                    }
4438                }
4439            }
4440            //unlockThread();
4441        } else {
4442            // Deterministic parallel
4443            if (tree_->size() < CoinMax(numberThreads_, 8) && !goneParallel) {
4444                node = tree_->bestNode(cutoff) ;
4445                // Possible one on tree worse than cutoff
4446                if (!node || node->objectiveValue() > cutoff)
4447                    continue;
4448                // Do main work of solving node here
4449                doOneNode(this, node, createdNode);
4450                assert (createdNode);
4451                if (!createdNode->active()) {
4452                    delete createdNode;
4453                    createdNode = NULL;
4454                } else {
4455                    // Say one more pointing to this
4456                    node->nodeInfo()->increment() ;
4457                    tree_->push(createdNode) ;
4458                }
4459                if (node->active()) {
4460                    assert (node->nodeInfo());
4461                    if (node->nodeInfo()->numberBranchesLeft()) {
4462                        tree_->push(node) ;
4463                    } else {
4464                        node->setActive(false);
4465                    }
4466                } else {
4467                    if (node->nodeInfo()) {
4468                        if (!node->nodeInfo()->numberBranchesLeft())
4469                            node->nodeInfo()->allBranchesGone(); // can clean up
4470                        // So will delete underlying stuff
4471                        node->setActive(true);
4472                    }
4473                    delNode[nDeleteNode++] = node;
4474                    node = NULL;
4475                }
4476                if (nDeleteNode >= MAX_DEL_NODE) {
4477                    for (int i = 0; i < nDeleteNode; i++) {
4478                        //printf("trying to del %d %x\n",i,delNode[i]);
4479                        delete delNode[i];
4480                        //printf("done to del %d %x\n",i,delNode[i]);
4481                    }
4482                    nDeleteNode = 0;
4483                }
4484            } else {
4485                // Split and solve
4486                master_->deterministicParallel();
4487                goneParallel = true;
4488            }
4489        }
4490#endif
4491    }
4492    if (nDeleteNode) {
4493        for (int i = 0; i < nDeleteNode; i++) {
4494            delete delNode[i];
4495        }
4496        nDeleteNode = 0;
4497    }
4498#ifdef CBC_THREAD
4499    if (master_) {
4500        master_->stopThreads(-1);
4501        master_->waitForThreadsInTree(2);
4502        delete master_;
4503        master_ = NULL;
4504        masterThread_ = NULL;
4505        // adjust time to allow for children on some systems
4506        //dblParam_[CbcStartSeconds] -= CoinCpuTimeJustChildren();
4507    }
4508#endif
4509    /*
4510      End of the non-abort actions. The next block of code is executed if we've
4511      aborted because we hit one of the limits. Clean up by deleting the live set
4512      and break out of the node processing loop. Note that on an abort, node may
4513      have been pushed back onto the tree for further processing, in which case
4514      it'll be deleted in cleanTree. We need to check.
4515    */
4516    if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
4517          numberSolutions_ < intParam_[CbcMaxNumSol] &&
4518          !maximumSecondsReached() &&
4519          !stoppedOnGap_ &&
4520          !eventHappened_ &&
4521          (maximumNumberIterations_ < 0 || numberIterations_ < maximumNumberIterations_))
4522         ) {
4523        if (tree_->size()) {
4524            double dummyBest;
4525            tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
4526        }
4527        delete nextRowCut_;
4528        /* order is important here:
4529         * maximumSecondsReached() should be checked before eventHappened_ and
4530         * isNodeLimitReached() should be checked after eventHappened_
4531         * reason is, that at timelimit, eventHappened_ is set to true to make Cbc stop fast
4532         *   and if Ctrl+C is hit, then the nodelimit is set to -1 to make Cbc stop
4533         */
4534        if (stoppedOnGap_) {
4535            messageHandler()->message(CBC_GAP, messages())
4536            << bestObjective_ - bestPossibleObjective_
4537            << dblParam_[CbcAllowableGap]
4538            << dblParam_[CbcAllowableFractionGap]*100.0
4539            << CoinMessageEol ;
4540            secondaryStatus_ = 2;
4541            status_ = 0 ;
4542        } else if (maximumSecondsReached()) {
4543            handler_->message(CBC_MAXTIME, messages_) << CoinMessageEol ;
4544            secondaryStatus_ = 4;
4545            status_ = 1 ;
4546        } else if (numberSolutions_ >= intParam_[CbcMaxNumSol]) {
4547            handler_->message(CBC_MAXSOLS, messages_) << CoinMessageEol ;
4548            secondaryStatus_ = 6;
4549            status_ = 1 ;
4550        } else if (isNodeLimitReached()) {
4551            handler_->message(CBC_MAXNODES, messages_) << CoinMessageEol ;
4552            secondaryStatus_ = 3;
4553            status_ = 1 ;
4554        } else if (maximumNumberIterations_ >= 0 && numberIterations_ >= maximumNumberIterations_) {
4555            handler_->message(CBC_MAXITERS, messages_) << CoinMessageEol ;
4556            secondaryStatus_ = 8;
4557            status_ = 1 ;
4558        } else {
4559            handler_->message(CBC_EVENT, messages_) << CoinMessageEol ;
4560            secondaryStatus_ = 5;
4561            status_ = 5 ;
4562        }
4563    }
4564    /*
4565      That's it, we've exhausted the search tree, or broken out of the loop because
4566      we hit some limit on evaluation.
4567
4568      We may have got an intelligent tree so give it one more chance
4569    */
4570    // Tell solver we are not in Branch and Cut
4571    solver_->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo, NULL) ;
4572    tree_->endSearch();
4573    //  If we did any sub trees - did we give up on any?
4574    if ( numberStoppedSubTrees_)
4575        status_ = 1;
4576    numberNodes_ += numberExtraNodes_;
4577    numberIterations_ += numberExtraIterations_;
4578    if (eventHandler) {
4579        eventHandler->event(CbcEventHandler::endSearch);
4580    }
4581    if (!status_) {
4582        // Set best possible unless stopped on gap
4583        if (secondaryStatus_ != 2)
4584            bestPossibleObjective_ = bestObjective_;
4585        handler_->message(CBC_END_GOOD, messages_)
4586        << bestObjective_ << numberIterations_ << numberNodes_ << getCurrentSeconds()
4587        << CoinMessageEol ;
4588    } else {
4589        handler_->message(CBC_END, messages_)
4590        << bestObjective_ << bestPossibleObjective_
4591        << numberIterations_ << numberNodes_ << getCurrentSeconds()
4592        << CoinMessageEol ;
4593    }
4594    if (numberStrongIterations_)
4595        handler_->message(CBC_STRONG_STATS, messages_)
4596        << strongInfo_[0] << numberStrongIterations_ << strongInfo_[2]
4597        << strongInfo_[1] << CoinMessageEol ;
4598    if (!numberExtraNodes_)
4599        handler_->message(CBC_OTHER_STATS, messages_)
4600        << maximumDepthActual_
4601        << numberDJFixed_ << CoinMessageEol ;
4602    else
4603        handler_->message(CBC_OTHER_STATS2, messages_)
4604        << maximumDepthActual_
4605        << numberDJFixed_ << numberExtraNodes_ << numberExtraIterations_
4606        << CoinMessageEol ;
4607    if (doStatistics == 100) {
4608        for (int i = 0; i < numberObjects_; i++) {
4609            CbcSimpleIntegerDynamicPseudoCost * obj =
4610                dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object_[i]) ;
4611            if (obj)
4612                obj->print();
4613        }
4614    }
4615    if (statistics_) {
4616        // report in some way
4617        int * lookup = new int[numberObjects_];
4618        int i;
4619        for (i = 0; i < numberObjects_; i++)
4620            lookup[i] = -1;
4621        bool goodIds = false; //true;
4622        for (i = 0; i < numberObjects_; i++) {
4623            int iColumn = object_[i]->columnNumber();
4624            if (iColumn >= 0 && iColumn < numberColumns) {
4625                if (lookup[i] == -1) {
4626                    lookup[i] = iColumn;
4627                } else {
4628                    goodIds = false;
4629                    break;
4630                }
4631            } else {
4632                goodIds = false;
4633                break;
4634            }
4635        }
4636        if (!goodIds) {
4637            delete [] lookup;
4638            lookup = NULL;
4639        }
4640        if (doStatistics >= 3) {
4641            printf("  node parent depth column   value                    obj      inf\n");
4642            for ( i = 0; i < numberNodes2_; i++) {
4643                statistics_[i]->print(lookup);
4644            }
4645        }
4646        if (doStatistics > 1) {
4647            // Find last solution
4648            int k;
4649            for (k = numberNodes2_ - 1; k >= 0; k--) {
4650                if (statistics_[k]->endingObjective() != COIN_DBL_MAX &&
4651                        !statistics_[k]->endingInfeasibility())
4652                    break;
4653            }
4654            if (k >= 0) {
4655                int depth = statistics_[k]->depth();
4656                int * which = new int[depth+1];
4657                for (i = depth; i >= 0; i--) {
4658                    which[i] = k;
4659                    k = statistics_[k]->parentNode();
4660                }
4661                printf("  node parent depth column   value                    obj      inf\n");
4662                for (i = 0; i <= depth; i++) {
4663                    statistics_[which[i]]->print(lookup);
4664                }
4665                delete [] which;
4666            }
4667        }
4668        // now summary
4669        int maxDepth = 0;
4670        double averageSolutionDepth = 0.0;
4671        int numberSolutions = 0;
4672        double averageCutoffDepth = 0.0;
4673        double averageSolvedDepth = 0.0;
4674        int numberCutoff = 0;
4675        int numberDown = 0;
4676        int numberFirstDown = 0;
4677        double averageInfDown = 0.0;
4678        double averageObjDown = 0.0;
4679        int numberCutoffDown = 0;
4680        int numberUp = 0;
4681        int numberFirstUp = 0;
4682        double averageInfUp = 0.0;
4683        double averageObjUp = 0.0;
4684        int numberCutoffUp = 0;
4685        double averageNumberIterations1 = 0.0;
4686        double averageValue = 0.0;
4687        for ( i = 0; i < numberNodes2_; i++) {
4688            int depth =  statistics_[i]->depth();
4689            int way =  statistics_[i]->way();
4690            double value = statistics_[i]->value();
4691            double startingObjective =  statistics_[i]->startingObjective();
4692            int startingInfeasibility = statistics_[i]->startingInfeasibility();
4693            double endingObjective = statistics_[i]->endingObjective();
4694            int endingInfeasibility = statistics_[i]->endingInfeasibility();
4695            maxDepth = CoinMax(depth, maxDepth);
4696            // Only for completed
4697            averageNumberIterations1 += statistics_[i]->numberIterations();
4698            averageValue += value;
4699            if (endingObjective != COIN_DBL_MAX && !endingInfeasibility) {
4700                numberSolutions++;
4701                averageSolutionDepth += depth;
4702            }
4703            if (endingObjective == COIN_DBL_MAX) {
4704                numberCutoff++;
4705                averageCutoffDepth += depth;
4706                if (way < 0) {
4707                    numberDown++;
4708                    numberCutoffDown++;
4709                    if (way == -1)
4710                        numberFirstDown++;
4711                } else {
4712                    numberUp++;
4713                    numberCutoffUp++;
4714                    if (way == 1)
4715                        numberFirstUp++;
4716                }
4717            } else {
4718                averageSolvedDepth += depth;
4719                if (way < 0) {
4720                    numberDown++;
4721                    averageInfDown += startingInfeasibility - endingInfeasibility;
4722                    averageObjDown += endingObjective - startingObjective;
4723                    if (way == -1)
4724                        numberFirstDown++;
4725                } else {
4726                    numberUp++;
4727                    averageInfUp += startingInfeasibility - endingInfeasibility;
4728                    averageObjUp += endingObjective - startingObjective;
4729                    if (way == 1)
4730                        numberFirstUp++;
4731                }
4732            }
4733        }
4734        // Now print
4735        if (numberSolutions)
4736            averageSolutionDepth /= static_cast<double> (numberSolutions);
4737        int numberSolved = numberNodes2_ - numberCutoff;
4738        double averageNumberIterations2 = numberIterations_ - averageNumberIterations1
4739                                          - numberIterationsAtContinuous;
4740        if (numberCutoff) {
4741            averageCutoffDepth /= static_cast<double> (numberCutoff);
4742            averageNumberIterations2 /= static_cast<double> (numberCutoff);
4743        }
4744        if (numberNodes2_)
4745            averageValue /= static_cast<double> (numberNodes2_);
4746        if (numberSolved) {
4747            averageNumberIterations1 /= static_cast<double> (numberSolved);
4748            averageSolvedDepth /= static_cast<double> (numberSolved);
4749        }
4750        printf("%d solution(s) were found (by branching) at an average depth of %g\n",
4751               numberSolutions, averageSolutionDepth);
4752        printf("average value of variable being branched on was %g\n",
4753               averageValue);
4754        printf("%d nodes were cutoff at an average depth of %g with iteration count of %g\n",
4755               numberCutoff, averageCutoffDepth, averageNumberIterations2);
4756        printf("%d nodes were solved at an average depth of %g with iteration count of %g\n",
4757               numberSolved, averageSolvedDepth, averageNumberIterations1);
4758        if (numberDown) {
4759            averageInfDown /= static_cast<double> (numberDown);
4760            averageObjDown /= static_cast<double> (numberDown);
4761        }
4762        printf("Down %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
4763               numberDown, numberFirstDown, numberDown - numberFirstDown, numberCutoffDown,
4764               averageInfDown, averageObjDown);
4765        if (numberUp) {
4766            averageInfUp /= static_cast<double> (numberUp);
4767            averageObjUp /= static_cast<double> (numberUp);
4768        }
4769        printf("Up %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
4770               numberUp, numberFirstUp, numberUp - numberFirstUp, numberCutoffUp,
4771               averageInfUp, averageObjUp);
4772        for ( i = 0; i < numberNodes2_; i++)
4773            delete statistics_[i];
4774        delete [] statistics_;
4775        statistics_ = NULL;
4776        maximumStatistics_ = 0;
4777        delete [] lookup;
4778    }
4779    /*
4780      If we think we have a solution, restore and confirm it with a call to
4781      setBestSolution().  We need to reset the cutoff value so as not to fathom
4782      the solution on bounds.  Note that calling setBestSolution( ..., true)
4783      leaves the continuousSolver_ bounds vectors fixed at the solution value.
4784
4785      Running resolve() here is a failsafe --- setBestSolution has already
4786      reoptimised using the continuousSolver_. If for some reason we fail to
4787      prove optimality, run the problem again after instructing the solver to
4788      tell us more.
4789
4790      If all looks good, replace solver_ with continuousSolver_, so that the
4791      outside world will be able to obtain information about the solution using
4792      public methods.
4793
4794      Don't replace if we are trying to save cuts
4795    */
4796    if (bestSolution_ && (solverCharacteristics_->solverType() < 2 || solverCharacteristics_->solverType() == 4) && 
4797        ((specialOptions_&8388608)==0||(specialOptions_&2048)!=0)) {
4798        setCutoff(1.0e50) ; // As best solution should be worse than cutoff
4799        // change cutoff as constraint if wanted
4800        if (cutoffRowNumber_>=0) {
4801          if (solver_->getNumRows()>cutoffRowNumber_)
4802            solver_->setRowUpper(cutoffRowNumber_,1.0e50);
4803        }
4804        // also in continuousSolver_
4805        if (continuousSolver_) {
4806          // Solvers know about direction
4807          double direction = solver_->getObjSense();
4808          continuousSolver_->setDblParam(OsiDualObjectiveLimit, 1.0e50*direction);
4809        }
4810        phase_ = 5;
4811        double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
4812        if ((specialOptions_&4) == 0)
4813            bestObjective_ += 100.0 * increment + 1.0e-3; // only set if we are going to solve
4814        setBestSolution(CBC_END_SOLUTION, bestObjective_, bestSolution_, 1) ;
4815        continuousSolver_->resolve() ;
4816        if (!continuousSolver_->isProvenOptimal()) {
4817            continuousSolver_->messageHandler()->setLogLevel(2) ;
4818            continuousSolver_->initialSolve() ;
4819        }
4820        delete solver_ ;
4821        // above deletes solverCharacteristics_
4822        solverCharacteristics_ = NULL;
4823        solver_ = continuousSolver_ ;
4824        setPointers(solver_);
4825        continuousSolver_ = NULL ;
4826    }
4827    /*
4828      Clean up dangling objects. continuousSolver_ may already be toast.
4829    */
4830    delete lastws ;
4831    if (saveObjects) {
4832        for (int i = 0; i < numberObjects_; i++)
4833            delete saveObjects[i];
4834        delete [] saveObjects;
4835    }
4836    numberStrong_ = saveNumberStrong;
4837    numberBeforeTrust_ = saveNumberBeforeTrust;
4838    delete [] whichGenerator_ ;
4839    whichGenerator_ = NULL;
4840    delete [] lowerBefore ;
4841    delete [] upperBefore ;
4842    delete [] walkback_ ;
4843    walkback_ = NULL ;
4844    delete [] lastNodeInfo_ ;
4845    lastNodeInfo_ = NULL;
4846    delete [] lastNumberCuts_ ;
4847    lastNumberCuts_ = NULL;
4848    delete [] lastCut_;
4849    lastCut_ = NULL;
4850    delete [] addedCuts_ ;
4851    addedCuts_ = NULL ;
4852    //delete persistentInfo;
4853    // Get rid of characteristics
4854    solverCharacteristics_ = NULL;
4855    if (continuousSolver_) {
4856        delete continuousSolver_ ;
4857        continuousSolver_ = NULL ;
4858    }
4859    /*
4860      Destroy global cuts by replacing with an empty OsiCuts object.
4861    */
4862    globalCuts_ = CbcRowCuts() ;
4863    delete globalConflictCuts_;
4864    globalConflictCuts_=NULL;
4865    if (!bestSolution_ && (specialOptions_&8388608)==0) {
4866        // make sure lp solver is infeasible
4867        int numberColumns = solver_->getNumCols();
4868        const double * columnLower = solver_->getColLower();
4869        int iColumn;
4870        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4871            if (solver_->isInteger(iColumn))
4872                solver_->setColUpper(iColumn, columnLower[iColumn]);
4873        }
4874        solver_->initialSolve();
4875    }
4876#ifdef COIN_HAS_CLP
4877    {
4878        OsiClpSolverInterface * clpSolver
4879        = dynamic_cast<OsiClpSolverInterface *> (solver_);
4880        if (clpSolver) {
4881            // Possible restore of pivot method
4882            if (savePivotMethod) {
4883                // model may have changed
4884                savePivotMethod->setModel(NULL);
4885                clpSolver->getModelPtr()->setDualRowPivotAlgorithm(*savePivotMethod);
4886                delete savePivotMethod;
4887            }
4888            clpSolver->setLargestAway(-1.0);
4889        }
4890    }
4891#endif
4892    if ((fastNodeDepth_ >= 1000000 || (moreSpecialOptions_&33554432)!=0)
4893         && !parentModel_) {
4894      // delete object off end
4895      delete object_[numberObjects_];
4896      if ((moreSpecialOptions_&33554432)==0)
4897        fastNodeDepth_ -= 1000000;
4898    }
4899    delete saveSolver;
4900    // Undo preprocessing performed during BaB.
4901    if (strategy_ && strategy_->preProcessState() > 0) {
4902        // undo preprocessing
4903        CglPreProcess * process = strategy_->process();
4904        assert (process);
4905        int n = originalSolver->getNumCols();
4906        if (bestSolution_) {
4907            delete [] bestSolution_;
4908            bestSolution_ = new double [n];
4909            process->postProcess(*solver_);
4910        }
4911        strategy_->deletePreProcess();
4912        // Solution now back in originalSolver
4913        delete solver_;
4914        solver_ = originalSolver;
4915        if (bestSolution_) {
4916            bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
4917            memcpy(bestSolution_, solver_->getColSolution(), n*sizeof(double));
4918        }
4919        // put back original objects if there were any
4920        if (originalObject) {
4921            int iColumn;
4922            assert (ownObjects_);
4923            for (iColumn = 0; iColumn < numberObjects_; iColumn++)
4924                delete object_[iColumn];
4925            delete [] object_;
4926            numberObjects_ = numberOriginalObjects;
4927            object_ = originalObject;
4928            delete [] integerVariable_;
4929            numberIntegers_ = 0;
4930            for (iColumn = 0; iColumn < n; iColumn++) {
4931                if (solver_->isInteger(iColumn))
4932                    numberIntegers_++;
4933            }
4934            integerVariable_ = new int[numberIntegers_];
4935            numberIntegers_ = 0;
4936            for (iColumn = 0; iColumn < n; iColumn++) {
4937                if (solver_->isInteger(iColumn))
4938                    integerVariable_[numberIntegers_++] = iColumn;
4939            }
4940        }
4941    }
4942    if (flipObjective)
4943      flipModel();
4944#ifdef COIN_HAS_CLP
4945    {
4946        OsiClpSolverInterface * clpSolver
4947        = dynamic_cast<OsiClpSolverInterface *> (solver_);
4948        if (clpSolver)
4949            clpSolver->setFakeObjective(reinterpret_cast<double *> (NULL));
4950    }
4951#endif
4952    moreSpecialOptions_ = saveMoreSpecialOptions;
4953    return ;
4954}
4955
4956
4957// Solve the initial LP relaxation
4958void
4959CbcModel::initialSolve()
4960{
4961    assert (solver_);
4962    // Double check optimization directions line up
4963    dblParam_[CbcOptimizationDirection] = solver_->getObjSense();
4964    // Check if bounds are all integral (as may get messed up later)
4965    checkModel();
4966    if (!solverCharacteristics_) {
4967        OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
4968        if (solverCharacteristics) {
4969            solverCharacteristics_ = solverCharacteristics;
4970        } else {
4971            // replace in solver
4972            OsiBabSolver defaultC;
4973            solver_->setAuxiliaryInfo(&defaultC);
4974            solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
4975        }
4976    }
4977    solverCharacteristics_->setSolver(solver_);
4978    solver_->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
4979    solver_->initialSolve();
4980    solver_->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo, NULL) ;
4981    if (!solver_->isProvenOptimal())
4982        solver_->resolve();
4983    // But set up so Jon Lee will be happy
4984    status_ = -1;
4985    secondaryStatus_ = -1;
4986    originalContinuousObjective_ = solver_->getObjValue() * solver_->getObjSense();
4987    bestPossibleObjective_ = originalContinuousObjective_;
4988    delete [] continuousSolution_;
4989    continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
4990                                          solver_->getNumCols());
4991    setPointers(solver_);
4992    solverCharacteristics_ = NULL;
4993}
4994
4995/*! \brief Get an empty basis object
4996
4997  Return an empty CoinWarmStartBasis object with the requested capacity,
4998  appropriate for the current solver. The object is cloned from the object
4999  cached as emptyWarmStart_. If there is no cached object, the routine
5000  queries the solver for a warm start object, empties it, and caches the
5001  result.
5002*/
5003
5004CoinWarmStartBasis *CbcModel::getEmptyBasis (int ns, int na) const
5005
5006{
5007    CoinWarmStartBasis *emptyBasis ;
5008    /*
5009      Acquire an empty basis object, if we don't yet have one.
5010    */
5011    if (emptyWarmStart_ == 0) {
5012        if (solver_ == 0) {
5013            throw CoinError("Cannot construct basis without solver!",
5014                            "getEmptyBasis", "CbcModel") ;
5015        }
5016        emptyBasis =
5017            dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart()) ;
5018        if (emptyBasis == 0) {
5019            throw CoinError(
5020                "Solver does not appear to use a basis-oriented warm start.",
5021                "getEmptyBasis", "CbcModel") ;
5022        }
5023        emptyBasis->setSize(0, 0) ;
5024        emptyWarmStart_ = dynamic_cast<CoinWarmStart *>(emptyBasis) ;
5025    }
5026    /*
5027      Clone the empty basis object, resize it as requested, and return.
5028    */
5029    emptyBasis = dynamic_cast<CoinWarmStartBasis *>(emptyWarmStart_->clone()) ;
5030    assert(emptyBasis) ;
5031    if (ns != 0 || na != 0) emptyBasis->setSize(ns, na) ;
5032
5033    return (emptyBasis) ;
5034}
5035
5036
5037/** Default Constructor
5038
5039  Creates an empty model without an associated solver.
5040*/
5041CbcModel::CbcModel()
5042
5043        :
5044        solver_(NULL),
5045        ownership_(0x80000000),
5046        continuousSolver_(NULL),
5047        referenceSolver_(NULL),
5048        defaultHandler_(true),
5049        emptyWarmStart_(NULL),
5050        bestObjective_(COIN_DBL_MAX),
5051        bestPossibleObjective_(COIN_DBL_MAX),
5052        sumChangeObjective1_(0.0),
5053        sumChangeObjective2_(0.0),
5054        bestSolution_(NULL),
5055        savedSolutions_(NULL),
5056        currentSolution_(NULL),
5057        testSolution_(NULL),
5058        globalConflictCuts_(NULL),
5059        minimumDrop_(1.0e-4),
5060        numberSolutions_(0),
5061        numberSavedSolutions_(0),
5062        maximumSavedSolutions_(0),
5063        stateOfSearch_(0),
5064        whenCuts_(-1),
5065        hotstartSolution_(NULL),
5066        hotstartPriorities_(NULL),
5067        numberHeuristicSolutions_(0),
5068        numberNodes_(0),
5069        numberNodes2_(0),
5070        numberIterations_(0),
5071        numberSolves_(0),
5072        status_(-1),
5073        secondaryStatus_(-1),
5074        numberIntegers_(0),
5075        numberRowsAtContinuous_(0),
5076        cutoffRowNumber_(-1),
5077        maximumNumberCuts_(0),
5078        phase_(0),
5079        currentNumberCuts_(0),
5080        maximumDepth_(0),
5081        walkback_(NULL),
5082        lastNodeInfo_(NULL),
5083        lastCut_(NULL),
5084        lastDepth_(0),
5085        lastNumberCuts2_(0),
5086        maximumCuts_(0),
5087        lastNumberCuts_(NULL),
5088        addedCuts_(NULL),
5089        nextRowCut_(NULL),
5090        currentNode_(NULL),
5091        integerVariable_(NULL),
5092        integerInfo_(NULL),
5093        continuousSolution_(NULL),
5094        usedInSolution_(NULL),
5095        specialOptions_(0),
5096        moreSpecialOptions_(0),
5097        topOfTree_(NULL),
5098        subTreeModel_(NULL),
5099        heuristicModel_(NULL),
5100        numberStoppedSubTrees_(0),
5101        presolve_(0),
5102        numberStrong_(5),
5103        numberBeforeTrust_(10),
5104        numberPenalties_(20),
5105        stopNumberIterations_(-1),
5106        penaltyScaleFactor_(3.0),
5107        numberAnalyzeIterations_(0),
5108        analyzeResults_(NULL),
5109        numberInfeasibleNodes_(0),
5110        problemType_(0),
5111        printFrequency_(0),
5112        numberCutGenerators_(0),
5113        generator_(NULL),
5114        virginGenerator_(NULL),
5115        numberHeuristics_(0),
5116        heuristic_(NULL),
5117        lastHeuristic_(NULL),
5118        fastNodeDepth_(-1),
5119        eventHandler_(NULL),
5120        numberObjects_(0),
5121        object_(NULL),
5122        ownObjects_(true),
5123        originalColumns_(NULL),
5124        howOftenGlobalScan_(3),
5125        numberGlobalViolations_(0),
5126        numberExtraIterations_(0),
5127        numberExtraNodes_(0),
5128        continuousObjective_(COIN_DBL_MAX),
5129        originalContinuousObjective_(COIN_DBL_MAX),
5130        continuousInfeasibilities_(COIN_INT_MAX),
5131        maximumCutPassesAtRoot_(20),
5132        maximumCutPasses_(10),
5133        preferredWay_(0),
5134        currentPassNumber_(0),
5135        maximumWhich_(1000),
5136        maximumRows_(0),
5137        randomSeed_(-1),
5138        multipleRootTries_(0),
5139        currentDepth_(0),
5140        whichGenerator_(NULL),
5141        maximumStatistics_(0),
5142        statistics_(NULL),
5143        maximumDepthActual_(0),
5144        numberDJFixed_(0.0),
5145        probingInfo_(NULL),
5146        numberFixedAtRoot_(0),
5147        numberFixedNow_(0),
5148        stoppedOnGap_(false),
5149        eventHappened_(false),
5150        numberLongStrong_(0),
5151        numberOldActiveCuts_(0),
5152        numberNewCuts_(0),
5153        searchStrategy_(-1),
5154        strongStrategy_(0),
5155        numberStrongIterations_(0),
5156        resolveAfterTakeOffCuts_(true),
5157        maximumNumberIterations_(-1),
5158        continuousPriority_(COIN_INT_MAX),
5159        numberUpdateItems_(0),
5160        maximumNumberUpdateItems_(0),
5161        updateItems_(NULL),
5162        storedRowCuts_(NULL),
5163        numberThreads_(0),
5164        threadMode_(0),
5165        master_(NULL),
5166        masterThread_(NULL)
5167{
5168    memset(intParam_, 0, sizeof(intParam_));
5169    intParam_[CbcMaxNumNode] = 2147483647;
5170    intParam_[CbcMaxNumSol] = 9999999;
5171
5172    memset(dblParam_, 0, sizeof(dblParam_));
5173    dblParam_[CbcIntegerTolerance] = 1e-6;
5174    dblParam_[CbcCutoffIncrement] = 1e-5;
5175    dblParam_[CbcAllowableGap] = 1.0e-10;
5176    dblParam_[CbcMaximumSeconds] = 1.0e100;
5177    dblParam_[CbcCurrentCutoff] = 1.0e100;
5178    dblParam_[CbcOptimizationDirection] = 1.0;
5179    dblParam_[CbcCurrentObjectiveValue] = 1.0e100;
5180    dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100;
5181    strongInfo_[0] = 0;
5182    strongInfo_[1] = 0;
5183    strongInfo_[2] = 0;
5184    strongInfo_[3] = 0;
5185    strongInfo_[4] = 0;
5186    strongInfo_[5] = 0;
5187    strongInfo_[6] = 0;
5188    solverCharacteristics_ = NULL;
5189    nodeCompare_ = new CbcCompareDefault();;
5190    problemFeasibility_ = new CbcFeasibilityBase();
5191    tree_ = new CbcTree();
5192    branchingMethod_ = NULL;
5193    cutModifier_ = NULL;
5194    strategy_ = NULL;
5195    parentModel_ = NULL;
5196    cbcColLower_ = NULL;
5197    cbcColUpper_ = NULL;
5198    cbcRowLower_ = NULL;
5199    cbcRowUpper_ = NULL;
5200    cbcColSolution_ = NULL;
5201    cbcRowPrice_ = NULL;
5202    cbcReducedCost_ = NULL;
5203    cbcRowActivity_ = NULL;
5204    appData_ = NULL;
5205    handler_ = new CoinMessageHandler();
5206    handler_->setLogLevel(2);
5207    messages_ = CbcMessage();
5208    //eventHandler_ = new CbcEventHandler() ;
5209}
5210
5211/** Constructor from solver.
5212
5213  Creates a model complete with a clone of the solver passed as a parameter.
5214*/
5215
5216CbcModel::CbcModel(const OsiSolverInterface &rhs)
5217        :
5218        continuousSolver_(NULL),
5219        referenceSolver_(NULL),
5220        defaultHandler_(true),
5221        emptyWarmStart_(NULL),
5222        bestObjective_(COIN_DBL_MAX),
5223        bestPossibleObjective_(COIN_DBL_MAX),
5224        sumChangeObjective1_(0.0),
5225        sumChangeObjective2_(0.0),
5226        globalConflictCuts_(NULL),
5227        minimumDrop_(1.0e-4),
5228        numberSolutions_(0),
5229        numberSavedSolutions_(0),
5230        maximumSavedSolutions_(0),
5231        stateOfSearch_(0),
5232        whenCuts_(-1),
5233        hotstartSolution_(NULL),
5234        hotstartPriorities_(NULL),
5235        numberHeuristicSolutions_(0),
5236        numberNodes_(0),
5237        numberNodes2_(0),
5238        numberIterations_(0),
5239        numberSolves_(0),
5240        status_(-1),
5241        secondaryStatus_(-1),
5242        numberRowsAtContinuous_(0),
5243        cutoffRowNumber_(-1),
5244        maximumNumberCuts_(0),
5245        phase_(0),
5246        currentNumberCuts_(0),
5247        maximumDepth_(0),
5248        walkback_(NULL),
5249        lastNodeInfo_(NULL),
5250        lastCut_(NULL),
5251        lastDepth_(0),
5252        lastNumberCuts2_(0),
5253        maximumCuts_(0),
5254        lastNumberCuts_(NULL),
5255        addedCuts_(NULL),
5256        nextRowCut_(NULL),
5257        currentNode_(NULL),
5258        integerInfo_(NULL