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

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

to print out saying optimal in more cases

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