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

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

allow sos in lp files and fix odd SOS branching

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