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

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

fix multiple solvers in stable

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