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

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

fix minor stuff like provenoptimal

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