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

Last change on this file since 1791 was 1791, checked in by stefan, 7 years ago

merge r1790 from stable/2.7 (correct message if iterlim reached) and introduce secondaryStatus 8 for stop at iteration limit

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 663.4 KB
Line 
1/* $Id: CbcModel.cpp 1791 2012-06-08 15:15:10Z stefan $ */
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_ &&
3919          !eventHappened_ &&
3920          (maximumNumberIterations_ < 0 || numberIterations_ < maximumNumberIterations_))
3921         ) {
3922        if (tree_->size()) {
3923            double dummyBest;
3924            tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
3925        }
3926        delete nextRowCut_;
3927        /* order is important here:
3928         * maximumSecondsReached() should be checked before eventHappened_ and
3929         * isNodeLimitReached() should be checked after eventHappened_
3930         * reason is, that at timelimit, eventHappened_ is set to true to make Cbc stop fast
3931         *   and if Ctrl+C is hit, then the nodelimit is set to -1 to make Cbc stop
3932         */
3933        if (stoppedOnGap_) {
3934            messageHandler()->message(CBC_GAP, messages())
3935            << bestObjective_ - bestPossibleObjective_
3936            << dblParam_[CbcAllowableGap]
3937            << dblParam_[CbcAllowableFractionGap]*100.0
3938            << CoinMessageEol ;
3939            secondaryStatus_ = 2;
3940            status_ = 0 ;
3941        } else if (maximumSecondsReached()) {
3942            handler_->message(CBC_MAXTIME, messages_) << CoinMessageEol ;
3943            secondaryStatus_ = 4;
3944            status_ = 1 ;
3945        } else if (eventHappened_) {
3946            handler_->message(CBC_EVENT, messages_) << CoinMessageEol ;
3947            secondaryStatus_ = 5;
3948            status_ = 5 ;
3949        } else if (isNodeLimitReached()) {
3950            handler_->message(CBC_MAXNODES, messages_) << CoinMessageEol ;
3951            secondaryStatus_ = 3;
3952            status_ = 1 ;
3953        } else if (maximumNumberIterations_ >= 0 && numberIterations_ >= maximumNumberIterations_) {
3954            handler_->message(CBC_MAXITERS, messages_) << CoinMessageEol ;
3955            secondaryStatus_ = 8;
3956            status_ = 1 ;
3957        } else {
3958            assert(numberSolutions_ >= intParam_[CbcMaxNumSol]);
3959            handler_->message(CBC_MAXSOLS, messages_) << CoinMessageEol ;
3960            secondaryStatus_ = 6;
3961            status_ = 1 ;
3962        }
3963    }
3964    /*
3965      That's it, we've exhausted the search tree, or broken out of the loop because
3966      we hit some limit on evaluation.
3967
3968      We may have got an intelligent tree so give it one more chance
3969    */
3970    // Tell solver we are not in Branch and Cut
3971    solver_->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo, NULL) ;
3972    tree_->endSearch();
3973    //  If we did any sub trees - did we give up on any?
3974    if ( numberStoppedSubTrees_)
3975        status_ = 1;
3976    numberNodes_ += numberExtraNodes_;
3977    numberIterations_ += numberExtraIterations_;
3978    if (eventHandler) {
3979        eventHandler->event(CbcEventHandler::endSearch);
3980    }
3981    if (!status_) {
3982        // Set best possible unless stopped on gap
3983        if (secondaryStatus_ != 2)
3984            bestPossibleObjective_ = bestObjective_;
3985        handler_->message(CBC_END_GOOD, messages_)
3986        << bestObjective_ << numberIterations_ << numberNodes_ << getCurrentSeconds()
3987        << CoinMessageEol ;
3988    } else {
3989        handler_->message(CBC_END, messages_)
3990        << bestObjective_ << bestPossibleObjective_
3991        << numberIterations_ << numberNodes_ << getCurrentSeconds()
3992        << CoinMessageEol ;
3993    }
3994    if (numberStrongIterations_)
3995        handler_->message(CBC_STRONG_STATS, messages_)
3996        << strongInfo_[0] << numberStrongIterations_ << strongInfo_[2]
3997        << strongInfo_[1] << CoinMessageEol ;
3998    if (!numberExtraNodes_)
3999        handler_->message(CBC_OTHER_STATS, messages_)
4000        << maximumDepthActual_
4001        << numberDJFixed_ << CoinMessageEol ;
4002    else
4003        handler_->message(CBC_OTHER_STATS2, messages_)
4004        << maximumDepthActual_
4005        << numberDJFixed_ << numberExtraNodes_ << numberExtraIterations_
4006        << CoinMessageEol ;
4007    if (doStatistics == 100) {
4008        for (int i = 0; i < numberObjects_; i++) {
4009            CbcSimpleIntegerDynamicPseudoCost * obj =
4010                dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object_[i]) ;
4011            if (obj)
4012                obj->print();
4013        }
4014    }
4015    if (statistics_) {
4016        // report in some way
4017        int * lookup = new int[numberObjects_];
4018        int i;
4019        for (i = 0; i < numberObjects_; i++)
4020            lookup[i] = -1;
4021        bool goodIds = false; //true;
4022        for (i = 0; i < numberObjects_; i++) {
4023            int iColumn = object_[i]->columnNumber();
4024            if (iColumn >= 0 && iColumn < numberColumns) {
4025                if (lookup[i] == -1) {
4026                    lookup[i] = iColumn;
4027                } else {
4028                    goodIds = false;
4029                    break;
4030                }
4031            } else {
4032                goodIds = false;
4033                break;
4034            }
4035        }
4036        if (!goodIds) {
4037            delete [] lookup;
4038            lookup = NULL;
4039        }
4040        if (doStatistics >= 3) {
4041            printf("  node parent depth column   value                    obj      inf\n");
4042            for ( i = 0; i < numberNodes2_; i++) {
4043                statistics_[i]->print(lookup);
4044            }
4045        }
4046        if (doStatistics > 1) {
4047            // Find last solution
4048            int k;
4049            for (k = numberNodes2_ - 1; k >= 0; k--) {
4050                if (statistics_[k]->endingObjective() != COIN_DBL_MAX &&
4051                        !statistics_[k]->endingInfeasibility())
4052                    break;
4053            }
4054            if (k >= 0) {
4055                int depth = statistics_[k]->depth();
4056                int * which = new int[depth+1];
4057                for (i = depth; i >= 0; i--) {
4058                    which[i] = k;
4059                    k = statistics_[k]->parentNode();
4060                }
4061                printf("  node parent depth column   value                    obj      inf\n");
4062                for (i = 0; i <= depth; i++) {
4063                    statistics_[which[i]]->print(lookup);
4064                }
4065                delete [] which;
4066            }
4067        }
4068        // now summary
4069        int maxDepth = 0;
4070        double averageSolutionDepth = 0.0;
4071        int numberSolutions = 0;
4072        double averageCutoffDepth = 0.0;
4073        double averageSolvedDepth = 0.0;
4074        int numberCutoff = 0;
4075        int numberDown = 0;
4076        int numberFirstDown = 0;
4077        double averageInfDown = 0.0;
4078        double averageObjDown = 0.0;
4079        int numberCutoffDown = 0;
4080        int numberUp = 0;
4081        int numberFirstUp = 0;
4082        double averageInfUp = 0.0;
4083        double averageObjUp = 0.0;
4084        int numberCutoffUp = 0;
4085        double averageNumberIterations1 = 0.0;
4086        double averageValue = 0.0;
4087        for ( i = 0; i < numberNodes2_; i++) {
4088            int depth =  statistics_[i]->depth();
4089            int way =  statistics_[i]->way();
4090            double value = statistics_[i]->value();
4091            double startingObjective =  statistics_[i]->startingObjective();
4092            int startingInfeasibility = statistics_[i]->startingInfeasibility();
4093            double endingObjective = statistics_[i]->endingObjective();
4094            int endingInfeasibility = statistics_[i]->endingInfeasibility();
4095            maxDepth = CoinMax(depth, maxDepth);
4096            // Only for completed
4097            averageNumberIterations1 += statistics_[i]->numberIterations();
4098            averageValue += value;
4099            if (endingObjective != COIN_DBL_MAX && !endingInfeasibility) {
4100                numberSolutions++;
4101                averageSolutionDepth += depth;
4102            }
4103            if (endingObjective == COIN_DBL_MAX) {
4104                numberCutoff++;
4105                averageCutoffDepth += depth;
4106                if (way < 0) {
4107                    numberDown++;
4108                    numberCutoffDown++;
4109                    if (way == -1)
4110                        numberFirstDown++;
4111                } else {
4112                    numberUp++;
4113                    numberCutoffUp++;
4114                    if (way == 1)
4115                        numberFirstUp++;
4116                }
4117            } else {
4118                averageSolvedDepth += depth;
4119                if (way < 0) {
4120                    numberDown++;
4121                    averageInfDown += startingInfeasibility - endingInfeasibility;
4122                    averageObjDown += endingObjective - startingObjective;
4123                    if (way == -1)
4124                        numberFirstDown++;
4125                } else {
4126                    numberUp++;
4127                    averageInfUp += startingInfeasibility - endingInfeasibility;
4128                    averageObjUp += endingObjective - startingObjective;
4129                    if (way == 1)
4130                        numberFirstUp++;
4131                }
4132            }
4133        }
4134        // Now print
4135        if (numberSolutions)
4136            averageSolutionDepth /= static_cast<double> (numberSolutions);
4137        int numberSolved = numberNodes2_ - numberCutoff;
4138        double averageNumberIterations2 = numberIterations_ - averageNumberIterations1
4139                                          - numberIterationsAtContinuous;
4140        if (numberCutoff) {
4141            averageCutoffDepth /= static_cast<double> (numberCutoff);
4142            averageNumberIterations2 /= static_cast<double> (numberCutoff);
4143        }
4144        if (numberNodes2_)
4145            averageValue /= static_cast<double> (numberNodes2_);
4146        if (numberSolved) {
4147            averageNumberIterations1 /= static_cast<double> (numberSolved);
4148            averageSolvedDepth /= static_cast<double> (numberSolved);
4149        }
4150        printf("%d solution(s) were found (by branching) at an average depth of %g\n",
4151               numberSolutions, averageSolutionDepth);
4152        printf("average value of variable being branched on was %g\n",
4153               averageValue);
4154        printf("%d nodes were cutoff at an average depth of %g with iteration count of %g\n",
4155               numberCutoff, averageCutoffDepth, averageNumberIterations2);
4156        printf("%d nodes were solved at an average depth of %g with iteration count of %g\n",
4157               numberSolved, averageSolvedDepth, averageNumberIterations1);
4158        if (numberDown) {
4159            averageInfDown /= static_cast<double> (numberDown);
4160            averageObjDown /= static_cast<double> (numberDown);
4161        }
4162        printf("Down %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
4163               numberDown, numberFirstDown, numberDown - numberFirstDown, numberCutoffDown,
4164               averageInfDown, averageObjDown);
4165        if (numberUp) {
4166            averageInfUp /= static_cast<double> (numberUp);
4167            averageObjUp /= static_cast<double> (numberUp);
4168        }
4169        printf("Up %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
4170               numberUp, numberFirstUp, numberUp - numberFirstUp, numberCutoffUp,
4171               averageInfUp, averageObjUp);
4172        for ( i = 0; i < numberNodes2_; i++)
4173            delete statistics_[i];
4174        delete [] statistics_;
4175        statistics_ = NULL;
4176        maximumStatistics_ = 0;
4177        delete [] lookup;
4178    }
4179    /*
4180      If we think we have a solution, restore and confirm it with a call to
4181      setBestSolution().  We need to reset the cutoff value so as not to fathom
4182      the solution on bounds.  Note that calling setBestSolution( ..., true)
4183      leaves the continuousSolver_ bounds vectors fixed at the solution value.
4184
4185      Running resolve() here is a failsafe --- setBestSolution has already
4186      reoptimised using the continuousSolver_. If for some reason we fail to
4187      prove optimality, run the problem again after instructing the solver to
4188      tell us more.
4189
4190      If all looks good, replace solver_ with continuousSolver_, so that the
4191      outside world will be able to obtain information about the solution using
4192      public methods.
4193    */
4194    if (bestSolution_ && (solverCharacteristics_->solverType() < 2 || solverCharacteristics_->solverType() == 4)) {
4195        setCutoff(1.0e50) ; // As best solution should be worse than cutoff
4196        // also in continuousSolver_
4197        if (continuousSolver_) {
4198          // Solvers know about direction
4199          double direction = solver_->getObjSense();
4200          continuousSolver_->setDblParam(OsiDualObjectiveLimit, 1.0e50*direction);
4201        }
4202        phase_ = 5;
4203        double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
4204        if ((specialOptions_&4) == 0)
4205            bestObjective_ += 100.0 * increment + 1.0e-3; // only set if we are going to solve
4206        setBestSolution(CBC_END_SOLUTION, bestObjective_, bestSolution_, 1) ;
4207        continuousSolver_->resolve() ;
4208        if (!continuousSolver_->isProvenOptimal()) {
4209            continuousSolver_->messageHandler()->setLogLevel(2) ;
4210            continuousSolver_->initialSolve() ;
4211        }
4212        delete solver_ ;
4213        // above deletes solverCharacteristics_
4214        solverCharacteristics_ = NULL;
4215        solver_ = continuousSolver_ ;
4216        setPointers(solver_);
4217        continuousSolver_ = NULL ;
4218    }
4219    /*
4220      Clean up dangling objects. continuousSolver_ may already be toast.
4221    */
4222    delete lastws ;
4223    if (saveObjects) {
4224        for (int i = 0; i < numberObjects_; i++)
4225            delete saveObjects[i];
4226        delete [] saveObjects;
4227    }
4228    numberStrong_ = saveNumberStrong;
4229    numberBeforeTrust_ = saveNumberBeforeTrust;
4230    delete [] whichGenerator_ ;
4231    whichGenerator_ = NULL;
4232    delete [] lowerBefore ;
4233    delete [] upperBefore ;
4234    delete [] walkback_ ;
4235    walkback_ = NULL ;
4236    delete [] lastNodeInfo_ ;
4237    lastNodeInfo_ = NULL;
4238    delete [] lastNumberCuts_ ;
4239    lastNumberCuts_ = NULL;
4240    delete [] lastCut_;
4241    lastCut_ = NULL;
4242    delete [] addedCuts_ ;
4243    addedCuts_ = NULL ;
4244    //delete persistentInfo;
4245    // Get rid of characteristics
4246    solverCharacteristics_ = NULL;
4247    if (continuousSolver_) {
4248        delete continuousSolver_ ;
4249        continuousSolver_ = NULL ;
4250    }
4251    /*
4252      Destroy global cuts by replacing with an empty OsiCuts object.
4253    */
4254    globalCuts_ = OsiCuts() ;
4255    if (!bestSolution_) {
4256        // make sure lp solver is infeasible
4257        int numberColumns = solver_->getNumCols();
4258        const double * columnLower = solver_->getColLower();
4259        int iColumn;
4260        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4261            if (solver_->isInteger(iColumn))
4262                solver_->setColUpper(iColumn, columnLower[iColumn]);
4263        }
4264        solver_->initialSolve();
4265    }
4266#ifdef COIN_HAS_CLP
4267    {
4268        OsiClpSolverInterface * clpSolver
4269        = dynamic_cast<OsiClpSolverInterface *> (solver_);
4270        if (clpSolver) {
4271            // Possible restore of pivot method
4272            if (savePivotMethod) {
4273                // model may have changed
4274                savePivotMethod->setModel(NULL);
4275                clpSolver->getModelPtr()->setDualRowPivotAlgorithm(*savePivotMethod);
4276                delete savePivotMethod;
4277            }
4278            clpSolver->setLargestAway(-1.0);
4279        }
4280    }
4281#endif
4282    if (fastNodeDepth_ >= 1000000 && !parentModel_) {
4283        // delete object off end
4284        delete object_[numberObjects_];
4285        fastNodeDepth_ -= 1000000;
4286    }
4287    delete saveSolver;
4288    // Undo preprocessing performed during BaB.
4289    if (strategy_ && strategy_->preProcessState() > 0) {
4290        // undo preprocessing
4291        CglPreProcess * process = strategy_->process();
4292        assert (process);
4293        int n = originalSolver->getNumCols();
4294        if (bestSolution_) {
4295            delete [] bestSolution_;
4296            bestSolution_ = new double [n];
4297            process->postProcess(*solver_);
4298        }
4299        strategy_->deletePreProcess();
4300        // Solution now back in originalSolver
4301        delete solver_;
4302        solver_ = originalSolver;
4303        if (bestSolution_) {
4304            bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
4305            memcpy(bestSolution_, solver_->getColSolution(), n*sizeof(double));
4306        }
4307        // put back original objects if there were any
4308        if (originalObject) {
4309            int iColumn;
4310            assert (ownObjects_);
4311            for (iColumn = 0; iColumn < numberObjects_; iColumn++)
4312                delete object_[iColumn];
4313            delete [] object_;
4314            numberObjects_ = numberOriginalObjects;
4315            object_ = originalObject;
4316            delete [] integerVariable_;
4317            numberIntegers_ = 0;
4318            for (iColumn = 0; iColumn < n; iColumn++) {
4319                if (solver_->isInteger(iColumn))
4320                    numberIntegers_++;
4321            }
4322            integerVariable_ = new int[numberIntegers_];
4323            numberIntegers_ = 0;
4324            for (iColumn = 0; iColumn < n; iColumn++) {
4325                if (solver_->isInteger(iColumn))
4326                    integerVariable_[numberIntegers_++] = iColumn;
4327            }
4328        }
4329    }
4330#ifdef COIN_HAS_CLP
4331    {
4332        OsiClpSolverInterface * clpSolver
4333        = dynamic_cast<OsiClpSolverInterface *> (solver_);
4334        if (clpSolver)
4335            clpSolver->setFakeObjective(reinterpret_cast<double *> (NULL));
4336    }
4337#endif
4338    moreSpecialOptions_ = saveMoreSpecialOptions;
4339    return ;
4340}
4341
4342
4343// Solve the initial LP relaxation
4344void
4345CbcModel::initialSolve()
4346{
4347    assert (solver_);
4348    // Double check optimization directions line up
4349    dblParam_[CbcOptimizationDirection] = solver_->getObjSense();
4350    // Check if bounds are all integral (as may get messed up later)
4351    checkModel();
4352    if (!solverCharacteristics_) {
4353        OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
4354        if (solverCharacteristics) {
4355            solverCharacteristics_ = solverCharacteristics;
4356        } else {
4357            // replace in solver
4358            OsiBabSolver defaultC;
4359            solver_->setAuxiliaryInfo(&defaultC);
4360            solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
4361        }
4362    }
4363    solverCharacteristics_->setSolver(solver_);
4364    solver_->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
4365    solver_->initialSolve();
4366    solver_->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo, NULL) ;
4367    if (!solver_->isProvenOptimal())
4368        solver_->resolve();
4369    // But set up so Jon Lee will be happy
4370    status_ = -1;
4371    secondaryStatus_ = -1;
4372    originalContinuousObjective_ = solver_->getObjValue() * solver_->getObjSense();
4373    bestPossibleObjective_ = originalContinuousObjective_;
4374    delete [] continuousSolution_;
4375    continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
4376                                          solver_->getNumCols());
4377    setPointers(solver_);
4378    solverCharacteristics_ = NULL;
4379}
4380
4381/*! \brief Get an empty basis object
4382
4383  Return an empty CoinWarmStartBasis object with the requested capacity,
4384  appropriate for the current solver. The object is cloned from the object
4385  cached as emptyWarmStart_. If there is no cached object, the routine
4386  queries the solver for a warm start object, empties it, and caches the
4387  result.
4388*/
4389
4390CoinWarmStartBasis *CbcModel::getEmptyBasis (int ns, int na) const
4391
4392{
4393    CoinWarmStartBasis *emptyBasis ;
4394    /*
4395      Acquire an empty basis object, if we don't yet have one.
4396    */
4397    if (emptyWarmStart_ == 0) {
4398        if (solver_ == 0) {
4399            throw CoinError("Cannot construct basis without solver!",
4400                            "getEmptyBasis", "CbcModel") ;
4401        }
4402        emptyBasis =
4403            dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart()) ;
4404        if (emptyBasis == 0) {
4405            throw CoinError(
4406                "Solver does not appear to use a basis-oriented warm start.",
4407                "getEmptyBasis", "CbcModel") ;
4408        }
4409        emptyBasis->setSize(0, 0) ;
4410        emptyWarmStart_ = dynamic_cast<CoinWarmStart *>(emptyBasis) ;
4411    }
4412    /*
4413      Clone the empty basis object, resize it as requested, and return.
4414    */
4415    emptyBasis = dynamic_cast<CoinWarmStartBasis *>(emptyWarmStart_->clone()) ;
4416    assert(emptyBasis) ;
4417    if (ns != 0 || na != 0) emptyBasis->setSize(ns, na) ;
4418
4419    return (emptyBasis) ;
4420}
4421
4422
4423/** Default Constructor
4424
4425  Creates an empty model without an associated solver.
4426*/
4427CbcModel::CbcModel()
4428
4429        :
4430        solver_(NULL),
4431        ownership_(0x80000000),
4432        continuousSolver_(NULL),
4433        referenceSolver_(NULL),
4434        defaultHandler_(true),
4435        emptyWarmStart_(NULL),
4436        bestObjective_(COIN_DBL_MAX),
4437        bestPossibleObjective_(COIN_DBL_MAX),
4438        sumChangeObjective1_(0.0),
4439        sumChangeObjective2_(0.0),
4440        bestSolution_(NULL),
4441        savedSolutions_(NULL),
4442        currentSolution_(NULL),
4443        testSolution_(NULL),
4444        minimumDrop_(1.0e-4),
4445        numberSolutions_(0),
4446        numberSavedSolutions_(0),
4447        maximumSavedSolutions_(0),
4448        stateOfSearch_(0),
4449        whenCuts_(-1),
4450        hotstartSolution_(NULL),
4451        hotstartPriorities_(NULL),
4452        numberHeuristicSolutions_(0),
4453        numberNodes_(0),
4454        numberNodes2_(0),
4455        numberIterations_(0),
4456        numberSolves_(0),
4457        status_(-1),
4458        secondaryStatus_(-1),
4459        numberIntegers_(0),
4460        numberRowsAtContinuous_(0),
4461        maximumNumberCuts_(0),
4462        phase_(0),
4463        currentNumberCuts_(0),
4464        maximumDepth_(0),
4465        walkback_(NULL),
4466        lastNodeInfo_(NULL),
4467        lastCut_(NULL),
4468        lastDepth_(0),
4469        lastNumberCuts2_(0),
4470        maximumCuts_(0),
4471        lastNumberCuts_(NULL),
4472        addedCuts_(NULL),
4473        nextRowCut_(NULL),
4474        currentNode_(NULL),
4475        integerVariable_(NULL),
4476        integerInfo_(NULL),
4477        continuousSolution_(NULL),
4478        usedInSolution_(NULL),
4479        specialOptions_(0),
4480        moreSpecialOptions_(0),
4481        subTreeModel_(NULL),
4482        numberStoppedSubTrees_(0),
4483        presolve_(0),
4484        numberStrong_(5),
4485        numberBeforeTrust_(10),
4486        numberPenalties_(20),
4487        stopNumberIterations_(-1),
4488        penaltyScaleFactor_(3.0),
4489        numberAnalyzeIterations_(0),
4490        analyzeResults_(NULL),
4491        numberInfeasibleNodes_(0),
4492        problemType_(0),
4493        printFrequency_(0),
4494        numberCutGenerators_(0),
4495        generator_(NULL),
4496        virginGenerator_(NULL),
4497        numberHeuristics_(0),
4498        heuristic_(NULL),
4499        lastHeuristic_(NULL),
4500        fastNodeDepth_(-1),
4501        eventHandler_(NULL),
4502        numberObjects_(0),
4503        object_(NULL),
4504        ownObjects_(true),
4505        originalColumns_(NULL),
4506        howOftenGlobalScan_(1),
4507        numberGlobalViolations_(0),
4508        numberExtraIterations_(0),
4509        numberExtraNodes_(0),
4510        continuousObjective_(COIN_DBL_MAX),
4511        originalContinuousObjective_(COIN_DBL_MAX),
4512        continuousInfeasibilities_(COIN_INT_MAX),
4513        maximumCutPassesAtRoot_(20),
4514        maximumCutPasses_(10),
4515        preferredWay_(0),
4516        currentPassNumber_(0),
4517        maximumWhich_(1000),
4518        maximumRows_(0),
4519        currentDepth_(0),
4520        whichGenerator_(NULL),
4521        maximumStatistics_(0),
4522        statistics_(NULL),
4523        maximumDepthActual_(0),
4524        numberDJFixed_(0.0),
4525        probingInfo_(NULL),
4526        numberFixedAtRoot_(0),
4527        numberFixedNow_(0),
4528        stoppedOnGap_(false),
4529        eventHappened_(false),
4530        numberLongStrong_(0),
4531        numberOldActiveCuts_(0),
4532        numberNewCuts_(0),
4533        searchStrategy_(-1),
4534        numberStrongIterations_(0),
4535        resolveAfterTakeOffCuts_(true),
4536        maximumNumberIterations_(-1),
4537        continuousPriority_(COIN_INT_MAX),
4538        numberUpdateItems_(0),
4539        maximumNumberUpdateItems_(0),
4540        updateItems_(NULL),
4541        storedRowCuts_(NULL),
4542        numberThreads_(0),
4543        threadMode_(0),
4544        master_(NULL),
4545        masterThread_(NULL)
4546{
4547    memset(intParam_, 0, sizeof(intParam_));
4548    intParam_[CbcMaxNumNode] = 2147483647;
4549    intParam_[CbcMaxNumSol] = 9999999;
4550
4551    memset(dblParam_, 0, sizeof(dblParam_));
4552    dblParam_[CbcIntegerTolerance] = 1e-6;
4553    dblParam_[CbcCutoffIncrement] = 1e-5;
4554    dblParam_[CbcAllowableGap] = 1.0e-10;
4555    dblParam_[CbcMaximumSeconds] = 1.0e100;
4556    dblParam_[CbcCurrentCutoff] = 1.0e100;
4557    dblParam_[CbcOptimizationDirection] = 1.0;
4558    dblParam_[CbcCurrentObjectiveValue] = 1.0e100;
4559    dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100;
4560    strongInfo_[0] = 0;
4561    strongInfo_[1] = 0;
4562    strongInfo_[2] = 0;
4563    strongInfo_[3] = 0;
4564    strongInfo_[4] = 0;
4565    strongInfo_[5] = 0;
4566    strongInfo_[6] = 0;
4567    solverCharacteristics_ = NULL;
4568    nodeCompare_ = new CbcCompareDefault();;
4569    problemFeasibility_ = new CbcFeasibilityBase();
4570    tree_ = new CbcTree();
4571    branchingMethod_ = NULL;
4572    cutModifier_ = NULL;
4573    strategy_ = NULL;
4574    parentModel_ = NULL;
4575    cbcColLower_ = NULL;
4576    cbcColUpper_ = NULL;
4577    cbcRowLower_ = NULL;
4578    cbcRowUpper_ = NULL;
4579    cbcColSolution_ = NULL;
4580    cbcRowPrice_ = NULL;
4581    cbcReducedCost_ = NULL;
4582    cbcRowActivity_ = NULL;
4583    appData_ = NULL;
4584    handler_ = new CoinMessageHandler();
4585    handler_->setLogLevel(2);
4586    messages_ = CbcMessage();
4587    //eventHandler_ = new CbcEventHandler() ;
4588}
4589
4590/** Constructor from solver.
4591
4592  Creates a model complete with a clone of the solver passed as a parameter.
4593*/
4594
4595CbcModel::CbcModel(const OsiSolverInterface &rhs)
4596        :
4597        continuousSolver_(NULL),
4598        referenceSolver_(NULL),
4599        defaultHandler_(true),
4600        emptyWarmStart_(NULL),
4601        bestObjective_(COIN_DBL_MAX),
4602        bestPossibleObjective_(COIN_DBL_MAX),
4603        sumChangeObjective1_(0.0),
4604        sumChangeObjective2_(0.0),
4605        minimumDrop_(1.0e-4),
4606        numberSolutions_(0),
4607        numberSavedSolutions_(0),
4608        maximumSavedSolutions_(0),
4609        stateOfSearch_(0),
4610        whenCuts_(-1),
4611        hotstartSolution_(NULL),
4612        hotstartPriorities_(NULL),
4613        numberHeuristicSolutions_(0),
4614        numberNodes_(0),
4615        numberNodes2_(0),
4616        numberIterations_(0),
4617        numberSolves_(0),
4618        status_(-1),
4619        secondaryStatus_(-1),
4620        numberRowsAtContinuous_(0),
4621        maximumNumberCuts_(0),
4622        phase_(0),
4623        currentNumberCuts_(0),
4624        maximumDepth_(0),
4625        walkback_(NULL),
4626        lastNodeInfo_(NULL),
4627        lastCut_(NULL),
4628        lastDepth_(0),
4629        lastNumberCuts2_(0),
4630        maximumCuts_(0),
4631        lastNumberCuts_(NULL),
4632        addedCuts_(NULL),
4633        nextRowCut_(NULL),
4634        currentNode_(NULL),
4635        integerInfo_(NULL),
4636        specialOptions_(0),
4637        moreSpecialOptions_(0),
4638        subTreeModel_(NULL),
4639        numberStoppedSubTrees_(0),
4640        presolve_(0),
4641        numberStrong_(5),
4642        numberBeforeTrust_(10),
4643        numberPenalties_(20),
4644        stopNumberIterations_(-1),
4645        penaltyScaleFactor_(3.0),
4646        numberAnalyzeIterations_(0),
4647        analyzeResults_(NULL),
4648        numberInfeasibleNodes_(0),
4649        problemType_(0),
4650        printFrequency_(0),
4651        numberCutGenerators_(0),
4652        generator_(NULL),
4653        virginGenerator_(NULL),
4654        numberHeuristics_(0),
4655        heuristic_(NULL),
4656        lastHeuristic_(NULL),
4657        fastNodeDepth_(-1),
4658        eventHandler_(NULL),
4659        numberObjects_(0),
4660        object_(NULL),
4661        ownObjects_(true),
4662        originalColumns_(NULL),
4663        howOftenGlobalScan_(1),
4664        numberGlobalViolations_(0),
4665        numberExtraIterations_(0),
4666        numberExtraNodes_(0),
4667        continuousObjective_(COIN_DBL_MAX),
4668        originalContinuousObjective_(COIN_DBL_MAX),
4669        continuousInfeasibilities_(COIN_INT_MAX),
4670        maximumCutPassesAtRoot_(20),
4671        maximumCutPasses_(10),
4672        preferredWay_(0),
4673        currentPassNumber_(0),
4674        maximumWhich_(1000),
4675        maximumRows_(0),
4676        currentDepth_(0),
4677        whichGenerator_(NULL),
4678        maximumStatistics_(0),
4679        statistics_(NULL),
4680        maximumDepthActual_(0),
4681        numberDJFixed_(0.0),
4682        probingInfo_(NULL),
4683        numberFixedAtRoot_(0),
4684        numberFixedNow_(0),
4685        stoppedOnGap_(false),
4686        eventHappened_(false),
4687        numberLongStrong_(0),
4688        numberOldActiveCuts_(0),
4689        numberNewCuts_(0),
4690        searchStrategy_(-1),
4691        numberStrongIterations_(0),
4692        resolveAfterTakeOffCuts_(true),
4693        maximumNumberIterations_(-1),
4694        continuousPriority_(COIN_INT_MAX),
4695        numberUpdateItems_(0),
4696        maximumNumberUpdateItems_(0),
4697        updateItems_(NULL),
4698        storedRowCuts_(NULL),
4699        numberThreads_(0),
4700        threadMode_(0),
4701        master_(NULL),
4702        masterThread_(NULL)
4703{
4704    memset(intParam_, 0, sizeof(intParam_));
4705    intParam_[CbcMaxNumNode] = 2147483647;
4706    intParam_[CbcMaxNumSol] = 9999999;
4707
4708    memset(dblParam_, 0, sizeof(dblParam_));
4709    dblParam_[CbcIntegerTolerance] = 1e-6;
4710    dblParam_[CbcCutoffIncrement] = 1e-5;
4711    dblParam_[CbcAllowableGap] = 1.0e-10;
4712    dblParam_[CbcMaximumSeconds] = 1.0e100;
4713    dblParam_[CbcCurrentCutoff] = 1.0e100;
4714    dblParam_[CbcOptimizationDirection] = 1.0;
4715    dblParam_[CbcCurrentObjectiveValue] = 1.0e100;
4716    dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100;
4717    strongInfo_[0] = 0;
4718    strongInfo_[1] = 0;
4719    strongInfo_[2] = 0;
4720    strongInfo_[3] = 0;
4721    strongInfo_[4] = 0;
4722    strongInfo_[5] = 0;
4723    strongInfo_[6] = 0;
4724    solverCharacteristics_ = NULL;
4725    nodeCompare_ = new CbcCompareDefault();;
4726    problemFeasibility_ = new CbcFeasibilityBase();
4727    tree_ = new CbcTree();
4728    branchingMethod_ = NULL;
4729    cutModifier_ = NULL;
4730    strategy_ = NULL;
4731    parentModel_ = NULL;
4732    appData_ = NULL;
4733    solver_ = rhs.clone();
4734    handler_ = new CoinMessageHandler();
4735    if (!solver_->defaultHandler()&&
4736        solver_->messageHandler()->logLevel(0)!=-1000)
4737      passInMessageHandler(solver_->messageHandler());
4738    handler_->setLogLevel(2);
4739    messages_ = CbcMessage();
4740    //eventHandler_ = new CbcEventHandler() ;
4741    referenceSolver_ = solver_->clone();
4742    ownership_ = 0x80000000;
4743    cbcColLower_ = NULL;
4744    cbcColUpper_ = NULL;
4745    cbcRowLower_ = NULL;
4746    cbcRowUpper_ = NULL;
4747    cbcColSolution_ = NULL;
4748    cbcRowPrice_ = NULL;
4749    cbcReducedCost_ = NULL;
4750    cbcRowActivity_ = NULL;
4751
4752    // Initialize solution and integer variable vectors
4753    bestSolution_ = NULL; // to say no solution found
4754    savedSolutions_ = NULL;
4755    numberIntegers_ = 0;
4756    int numberColumns = solver_->getNumCols();
4757    int iColumn;
4758    if (numberColumns) {
4759        // Space for current solution
4760        currentSolution_ = new double[numberColumns];
4761        continuousSolution_ = new double[numberColumns];
4762        usedInSolution_ = new int[numberColumns];
4763        CoinZeroN(usedInSolution_, numberColumns);
4764        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4765            if ( solver_->isInteger(iColumn))
4766                numberIntegers_++;
4767        }
4768    } else {
4769        // empty model
4770        currentSolution_ = NULL;
4771        continuousSolution_ = NULL;
4772        usedInSolution_ = NULL;
4773    }
4774    testSolution_ = currentSolution_;
4775    if (numberIntegers_) {
4776        integerVariable_ = new int [numberIntegers_];
4777        numberIntegers_ = 0;
4778        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4779            if ( solver_->isInteger(iColumn))
4780                integerVariable_[numberIntegers_++] = iColumn;
4781        }
4782    } else {
4783        integerVariable_ = NULL;
4784    }
4785}
4786
4787/*
4788  Assign a solver to the model (model assumes ownership)
4789
4790  The integer variable vector is initialized if it's not already present.
4791  If deleteSolver then current solver deleted (if model owned)
4792
4793  Assuming ownership matches usage in OsiSolverInterface
4794  (cf. assignProblem, loadProblem).
4795
4796  TODO: What to do about solver parameters? A simple copy likely won't do it,
4797        because the SI must push the settings into the underlying solver. In
4798        the context of switching solvers in cbc, this means that command line
4799        settings will get lost. Stash the command line somewhere and reread it
4800        here, maybe?
4801
4802  TODO: More generally, how much state should be transferred from the old
4803        solver to the new solver? Best perhaps to see how usage develops.
4804        What's done here mimics the CbcModel(OsiSolverInterface) constructor.
4805*/
4806void
4807CbcModel::assignSolver(OsiSolverInterface *&solver, bool deleteSolver)
4808
4809{
4810    // resize best solution if exists
4811    if (bestSolution_ && solver && solver_) {
4812        int nOld = solver_->getNumCols();
4813        int nNew = solver->getNumCols();
4814        if (nNew > nOld) {
4815            double * temp = new double[nNew];
4816            memcpy(temp, bestSolution_, nOld*sizeof(double));
4817            memset(temp + nOld, 0, (nNew - nOld)*sizeof(double));
4818            delete [] bestSolution_;
4819            bestSolution_ = temp;
4820        }
4821    }
4822    // Keep the current message level for solver (if solver exists)
4823    if (solver_)
4824        solver->messageHandler()->setLogLevel(solver_->messageHandler()->logLevel()) ;
4825
4826    if (modelOwnsSolver() && deleteSolver) {
4827        solverCharacteristics_ = NULL;
4828        delete solver_ ;
4829    }
4830    solver_ = solver;
4831    solver = NULL ;
4832    setModelOwnsSolver(true) ;
4833    /*
4834      Basis information is solver-specific.
4835    */
4836    if (emptyWarmStart_) {
4837        delete emptyWarmStart_  ;
4838        emptyWarmStart_ = 0 ;
4839    }
4840    bestSolutionBasis_ = CoinWarmStartBasis();
4841    /*
4842      Initialize integer variable vector.
4843    */
4844    numberIntegers_ = 0;
4845    int numberColumns = solver_->getNumCols();
4846    int iColumn;
4847    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4848        if ( solver_->isInteger(iColumn))
4849            numberIntegers_++;
4850    }
4851    delete [] integerVariable_;
4852    if (numberIntegers_) {
4853        integerVariable_ = new int [numberIntegers_];
4854        numberIntegers_ = 0;
4855        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4856            if ( solver_->isInteger(iColumn))
4857                integerVariable_[numberIntegers_++] = iColumn;
4858        }
4859    } else {
4860        integerVariable_ = NULL;
4861    }
4862
4863    return ;
4864}
4865
4866// Copy constructor.
4867
4868CbcModel::CbcModel(const CbcModel & rhs, bool cloneHandler)
4869        :
4870        continuousSolver_(NULL),
4871        referenceSolver_(NULL),
4872        defaultHandler_(rhs.defaultHandler_),
4873        emptyWarmStart_(NULL),
4874        bestObjective_(rhs.bestObjective_),
4875        bestPossibleObjective_(rhs.bestPossibleObjective_),
4876        sumChangeObjective1_(rhs.sumChangeObjective1_),
4877        sumChangeObjective2_(rhs.sumChangeObjective2_),
4878        minimumDrop_(rhs.minimumDrop_),
4879        numberSolutions_(rhs.numberSolutions_),
4880        numberSavedSolutions_(rhs.numberSavedSolutions_),
4881        maximumSavedSolutions_(rhs.maximumSavedSolutions_),
4882        stateOfSearch_(rhs.stateOfSearch_),
4883        whenCuts_(rhs.whenCuts_),
4884        numberHeuristicSolutions_(rhs.numberHeuristicSolutions_),
4885        numberNodes_(rhs.numberNodes_),
4886        numberNodes2_(rhs.numberNodes2_),
4887        numberIterations_(rhs.numberIterations_),
4888        numberSolves_(rhs.numberSolves_),
4889        status_(rhs.status_),
4890        secondaryStatus_(rhs.secondaryStatus_),
4891        specialOptions_(rhs.specialOptions_),
4892        moreSpecialOptions_(rhs.moreSpecialOptions_),
4893        subTreeModel_(rhs.subTreeModel_),
4894        numberStoppedSubTrees_(rhs.numberStoppedSubTrees_),
4895        presolve_(rhs.presolve_),
4896        numberStrong_(rhs.numberStrong_),
4897        numberBeforeTrust_(rhs.numberBeforeTrust_),
4898        numberPenalties_(rhs.numberPenalties_),
4899        stopNumberIterations_(rhs.stopNumberIterations_),
4900        penaltyScaleFactor_(rhs.penaltyScaleFactor_),
4901        numberAnalyzeIterations_(rhs.numberAnalyzeIterations_),
4902        analyzeResults_(NULL),
4903        numberInfeasibleNodes_(rhs.numberInfeasibleNodes_),
4904        problemType_(rhs.problemType_),
4905        printFrequency_(rhs.printFrequency_),
4906        fastNodeDepth_(rhs.fastNodeDepth_),
4907        howOftenGlobalScan_(rhs.howOftenGlobalScan_),
4908        numberGlobalViolations_(rhs.numberGlobalViolations_),
4909        numberExtraIterations_(rhs.numberExtraIterations_),
4910        numberExtraNodes_(rhs.numberExtraNodes_),
4911        continuousObjective_(rhs.continuousObjective_),
4912        originalContinuousObjective_(rhs.originalContinuousObjective_),
4913        continuousInfeasibilities_(rhs.continuousInfeasibilities_),
4914        maximumCutPassesAtRoot_(rhs.maximumCutPassesAtRoot_),
4915        maximumCutPasses_( rhs.maximumCutPasses_),
4916        preferredWay_(rhs.preferredWay_),
4917        currentPassNumber_(rhs.currentPassNumber_),
4918        maximumWhich_(rhs.maximumWhich_),
4919        maximumRows_(0),
4920        currentDepth_(0),
4921        whichGenerator_(NULL),
4922        maximumStatistics_(0),
4923        statistics_(NULL),
4924        maximumDepthActual_(0),
4925        numberDJFixed_(0.0),
4926        probingInfo_(NULL),
4927        numberFixedAtRoot_(rhs.numberFixedAtRoot_),
4928        numberFixedNow_(rhs.numberFixedNow_),
4929        stoppedOnGap_(rhs.stoppedOnGap_),
4930        eventHappened_(rhs.eventHappened_),
4931        numberLongStrong_(rhs.numberLongStrong_),
4932        numberOldActiveCuts_(rhs.numberOldActiveCuts_),
4933        numberNewCuts_(rhs.numberNewCuts_),
4934        searchStrategy_(rhs.searchStrategy_),
4935        numberStrongIterations_(rhs.numberStrongIterations_),
4936        resolveAfterTakeOffCuts_(rhs.resolveAfterTakeOffCuts_),
4937        maximumNumberIterations_(rhs.maximumNumberIterations_),
4938        continuousPriority_(rhs.continuousPriority_),
4939        numberUpdateItems_(rhs.numberUpdateItems_),
4940        maximumNumberUpdateItems_(rhs.maximumNumberUpdateItems_),
4941        updateItems_(NULL),
4942        storedRowCuts_(NULL),
4943        numberThreads_(rhs.numberThreads_),
4944        threadMode_(rhs.threadMode_),
4945        master_(NULL),
4946        masterThread_(NULL)
4947{
4948    memcpy(intParam_, rhs.intParam_, sizeof(intParam_));
4949    memcpy(dblParam_, rhs.dblParam_, sizeof(dblParam_));
4950    strongInfo_[0] = rhs.strongInfo_[0];
4951    strongInfo_[1] = rhs.strongInfo_[1];
4952    strongInfo_[2] = rhs.strongInfo_[2];
4953    strongInfo_[3] = rhs.strongInfo_[3];
4954    strongInfo_[4] = rhs.strongInfo_[4];
4955    strongInfo_[5] = rhs.strongInfo_[5];
4956    strongInfo_[6] = rhs.strongInfo_[6];
4957    solverCharacteristics_ = NULL;
4958    if (rhs.emptyWarmStart_) emptyWarmStart_ = rhs.emptyWarmStart_->clone() ;
4959    if (defaultHandler_ || cloneHandler) {
4960        handler_ = new CoinMessageHandler();
4961        handler_->setLogLevel(2);
4962    } else {
4963        handler_ = rhs.handler_;
4964    }
4965    messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
4966    numberCutGenerators_ = rhs.numberCutGenerators_;
4967    if (numberCutGenerators_) {
4968        generator_ = new CbcCutGenerator * [numberCutGenerators_];
4969        virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_];
4970        int i;
4971        for (i = 0; i < numberCutGenerators_; i++) {
4972            generator_[i] = new CbcCutGenerator(*rhs.generator_[i]);
4973            virginGenerator_[i] = new CbcCutGenerator(*rhs.virginGenerator_[i]);
4974        }
4975    } else {
4976        generator_ = NULL;
4977        virginGenerator_ = NULL;
4978    }
4979    globalCuts_ = rhs.globalCuts_;
4980    numberHeuristics_ = rhs.numberHeuristics_;
4981    if (numberHeuristics_) {
4982        heuristic_ = new CbcHeuristic * [numberHeuristics_];
4983        int i;
4984        for (i = 0; i < numberHeuristics_; i++) {
4985            heuristic_[i] = rhs.heuristic_[i]->clone();
4986        }
4987    } else {
4988        heuristic_ = NULL;
4989    }
4990    lastHeuristic_ = NULL;
4991    if (rhs.eventHandler_) {
4992        eventHandler_ = rhs.eventHandler_->clone() ;
4993    } else {
4994        eventHandler_ = NULL ;
4995    }
4996    ownObjects_ = rhs.ownObjects_;
4997    if (ownObjects_) {
4998        numberObjects_ = rhs.numberObjects_;
4999        if (numberObjects_) {
5000            object_ = new OsiObject * [numberObjects_];
5001            int i;
5002            for (i = 0; i < numberObjects_; i++) {
5003                object_[i] = (rhs.object_[i])->clone();
5004                CbcObject * obj = dynamic_cast <CbcObject *>(object_[i]) ;
5005                // Could be OsiObjects
5006                if (obj)
5007                    obj->setModel(this);
5008            }
5009        } else {
5010            object_ = NULL;
5011        }
5012    } else {
5013        // assume will be redone
5014        numberObjects_ = 0;
5015        object_ = NULL;
5016    }
5017    if (rhs.referenceSolver_)
5018        referenceSolver_ = rhs.referenceSolver_->clone();
5019    else
5020        referenceSolver_ = NULL;
5021    solver_ = rhs.solver_->clone();
5022    if (rhs.originalColumns_) {
5023        int numberColumns = solver_->getNumCols();
5024        originalColumns_ = new int [numberColumns];
5025        memcpy(originalColumns_, rhs.originalColumns_, numberColumns*sizeof(int));
5026    } else {
5027        originalColumns_ = NULL;
5028    }
5029    if (maximumNumberUpdateItems_) {
5030        updateItems_ = new CbcObjectUpdateData [maximumNumberUpdateItems_];
5031        for (int i = 0; i < maximumNumberUpdateItems_; i++)
5032            updateItems_[i] = rhs.updateItems_[i];
5033    }
5034    if (maximumWhich_ && rhs.whichGenerator_)
5035        whichGenerator_ = CoinCopyOfArray(rhs.whichGenerator_, maximumWhich_);
5036    nodeCompare_ = rhs.nodeCompare_->clone();
5037    problemFeasibility_ = rhs.problemFeasibility_->clone();
5038    tree_ = rhs.tree_->clone();
5039    if (rhs.branchingMethod_)
5040        branchingMethod_ = rhs.branchingMethod_->clone();
5041    else
5042        branchingMethod_ = NULL;
5043    if (rhs.cutModifier_)
5044        cutModifier_ = rhs.cutModifier_->clone();
5045    else
5046        cutModifier_ = NULL;
5047    cbcColLower_ = NULL;
5048    cbcColUpper_ = NULL;
5049    cbcRowLower_ = NULL;
5050    cbcRowUpper_ = NULL;
5051    cbcColSolution_ = NULL;
5052    cbcRowPrice_ = NULL;
5053    cbcReducedCost_ = NULL;
5054    cbcRowActivity_ = NULL;
5055    if (rhs.strategy_)
5056        strategy_ = rhs.strategy_->clone();
5057    else
5058        strategy_ = NULL;
5059    parentModel_ = rhs.parentModel_;
5060    appData_ = rhs.appData_;
5061    messages_ = rhs.messages_;
5062    ownership_ = rhs.ownership_ | 0x80000000;
5063    messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
5064    numberIntegers_ = rhs.numberIntegers_;
5065    randomNumberGenerator_ = rhs.randomNumberGenerator_;
5066    if (numberIntegers_) {
5067        integerVariable_ = new int [numberIntegers_];
5068        memcpy(integerVariable_, rhs.integerVariable_, numberIntegers_*sizeof(int));
5069        integerInfo_ = CoinCopyOfArray(rhs.integerInfo_, solver_->getNumCols());
5070    } else {
5071        integerVariable_ = NULL;
5072        integerInfo_ = NULL;
5073    }
5074    if (rhs.hotstartSolution_) {
5075        int numberColumns = solver_->getNumCols();
5076        hotstartSolution_ = CoinCopyOfArray(rhs.hotstartSolution_, numberColumns);
5077        hotstartPriorities_ = CoinCopyOfArray(rhs.hotstartPriorities_, numberColumns);
5078    } else {
5079        hotstartSolution_ = NULL;
5080        hotstartPriorities_ = NULL;
5081    }
5082    if (rhs.bestSolution_) {
5083        int numberColumns = solver_->getNumCols();
5084        bestSolution_ = new double[numberColumns];
5085        memcpy(bestSolution_, rhs.bestSolution_, numberColumns*sizeof(double));
5086    } else {
5087        bestSolution_ = NULL;
5088    }
5089    int numberColumns = solver_->getNumCols();
5090    if (maximumSavedSolutions_ && rhs.savedSolutions_) {
5091        savedSolutions_ = new double * [maximumSavedSolutions_];
5092        for (int i = 0; i < maximumSavedSolutions_; i++)
5093            savedSolutions_[i] = CoinCopyOfArray(rhs.savedSolutions_[i], numberColumns + 2);
5094    } else {
5095        savedSolutions_ = NULL;
5096    }
5097    // Space for current solution
5098    currentSolution_ = new double[numberColumns];
5099    continuousSolution_ = new double[numberColumns];
5100    usedInSolution_ = new int[numberColumns];
5101    CoinZeroN(usedInSolution_, numberColumns);
5102    testSolution_ = currentSolution_;
5103    numberRowsAtContinuous_ = rhs.numberRowsAtContinuous_;
5104    maximumNumberCuts_ = rhs.maximumNumberCuts_;
5105    phase_ = rhs.phase_;
5106    currentNumberCuts_ = rhs.currentNumberCuts_;
5107    maximumDepth_ = rhs.maximumDepth_;
5108    // These are only used as temporary arrays so need not be filled
5109    if (maximumNumberCuts_) {
5110        addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
5111    } else {
5112        addedCuts_ = NULL;
5113    }
5114    bestSolutionBasis_ = rhs.bestSolutionBasis_;
5115    nextRowCut_ = NULL;
5116    currentNode_ = NULL;
5117    if (maximumDepth_) {
5118        walkback_ = new CbcNodeInfo * [maximumDepth_];
5119        lastNodeInfo_ = new CbcNodeInfo * [maximumDepth_] ;
5120        lastNumberCuts_ = new int [maximumDepth_] ;
5121    } else {
5122        walkback_ = NULL;
5123        lastNodeInfo_ = NULL;
5124        lastNumberCuts_ = NULL;
5125    }
5126    maximumCuts_ = rhs.maximumCuts_;
5127    if (maximumCuts_) {
5128        lastCut_ = new const OsiRowCut * [maximumCuts_] ;
5129    } else {
5130        lastCut_ = NULL;
5131    }
5132    synchronizeModel();
5133    if (cloneHandler && !defaultHandler_) {
5134        delete handler_;
5135        CoinMessageHandler * handler = rhs.handler_->clone();
5136        passInMessageHandler(handler);
5137    }
5138}
5139
5140// Assignment operator
5141CbcModel &
5142CbcModel::operator=(const CbcModel & rhs)
5143{
5144    if (this != &rhs) {
5145        if (modelOwnsSolver()) {
5146            solverCharacteristics_ = NULL;
5147            delete solver_;
5148            solver_ = NULL;
5149        }
5150        gutsOfDestructor();
5151        if (defaultHandler_) {
5152            delete handler_;
5153            handler_ = NULL;
5154        }
5155        defaultHandler_ = rhs.defaultHandler_;
5156        if (defaultHandler_) {
5157            handler_ = new CoinMessageHandler();
5158            handler_->setLogLevel(2);
5159        } else {
5160            handler_ = rhs.handler_;
5161        }
5162        messages_ = rhs.messages_;
5163        messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
5164        if (rhs.solver_) {
5165            solver_ = rhs.solver_->clone() ;
5166        } else {
5167            solver_ = 0 ;
5168        }
5169        ownership_ = 0x80000000;
5170        delete continuousSolver_ ;
5171        if (rhs.continuousSolver_) {
5172            continuousSolver_ = rhs.continuousSolver_->clone() ;
5173        } else {
5174            continuousSolver_ = 0 ;
5175        }
5176        delete referenceSolver_;
5177        if (rhs.referenceSolver_) {
5178            referenceSolver_ = rhs.referenceSolver_->clone() ;
5179        } else {
5180            referenceSolver_ = NULL ;
5181        }
5182
5183        delete emptyWarmStart_ ;
5184        if (rhs.emptyWarmStart_) {
5185            emptyWarmStart_ = rhs.emptyWarmStart_->clone() ;
5186        } else {
5187            emptyWarmStart_ = 0 ;
5188        }
5189
5190        bestObjective_ = rhs.bestObjective_;
5191        bestPossibleObjective_ = rhs.bestPossibleObjective_;
5192        sumChangeObjective1_ = rhs.sumChangeObjective1_;
5193        sumChangeObjective2_ = rhs.sumChangeObjective2_;
5194        delete [] bestSolution_;
5195        if (rhs.bestSolution_) {
5196            int numberColumns = rhs.getNumCols();
5197            bestSolution_ = new double[numberColumns];
5198            memcpy(bestSolution_, rhs.bestSolution_, numberColumns*sizeof(double));
5199        } else {
5200            bestSolution_ = NULL;
5201        }
5202        for (int i = 0; i < maximumSavedSolutions_; i++)
5203            delete [] savedSolutions_[i];
5204        delete [] savedSolutions_;
5205        savedSolutions_ = NULL;
5206        int numberColumns = rhs.getNumCols();
5207        if (numberColumns) {
5208            // Space for current solution
5209            currentSolution_ = new double[numberColumns];
5210            continuousSolution_ = new double[numberColumns];
5211            usedInSolution_ = new int[numberColumns];
5212            CoinZeroN(usedInSolution_, numberColumns);
5213        } else {
5214            currentSolution_ = NULL;
5215            continuousSolution_ = NULL;
5216            usedInSolution_ = NULL;
5217        }
5218        if (maximumSavedSolutions_) {
5219            savedSolutions_ = new double * [maximumSavedSolutions_];
5220            for (int i = 0; i < maximumSavedSolutions_; i++)
5221                savedSolutions_[i] = CoinCopyOfArray(rhs.savedSolutions_[i], numberColumns + 2);
5222        } else {
5223            savedSolutions_ = NULL;
5224        }
5225        testSolution_ = currentSolution_;
5226        minimumDrop_ = rhs.minimumDrop_;
5227        numberSolutions_ = rhs.numberSolutions_;
5228        numberSavedSolutions_ = rhs.numberSavedSolutions_;
5229        maximumSavedSolutions_ = rhs.maximumSavedSolutions_;
5230        stateOfSearch_ = rhs.stateOfSearch_;
5231        whenCuts_ = rhs.whenCuts_;
5232        numberHeuristicSolutions_ = rhs.numberHeuristicSolutions_;
5233        numberNodes_ = rhs.numberNodes_;
5234        numberNodes2_ = rhs.numberNodes2_;
5235        numberIterations_ = rhs.numberIterations_;
5236        numberSolves_ = rhs.numberSolves_;
5237        status_ = rhs.status_;
5238        secondaryStatus_ = rhs.secondaryStatus_;
5239        specialOptions_ = rhs.specialOptions_;
5240        moreSpecialOptions_ = rhs.moreSpecialOptions_;
5241        subTreeModel_ = rhs.subTreeModel_;
5242        numberStoppedSubTrees_ = rhs.numberStoppedSubTrees_;
5243        presolve_ = rhs.presolve_;
5244        numberStrong_ = rhs.numberStrong_;
5245        numberBeforeTrust_ = rhs.numberBeforeTrust_;
5246        numberPenalties_ = rhs.numberPenalties_;
5247        stopNumberIterations_ = rhs.stopNumberIterations_;
5248        penaltyScaleFactor_ = rhs.penaltyScaleFactor_;
5249        numberAnalyzeIterations_ = rhs.numberAnalyzeIterations_;
5250        delete [] analyzeResults_;
5251        analyzeResults_ = NULL;
5252        numberInfeasibleNodes_ = rhs.numberInfeasibleNodes_;
5253        problemType_ = rhs.problemType_;
5254        printFrequency_ = rhs.printFrequency_;
5255        howOftenGlobalScan_ = rhs.howOftenGlobalScan_;
5256        numberGlobalViolations_ = rhs.numberGlobalViolations_;
5257        numberExtraIterations_ = rhs.numberExtraIterations_;
5258        numberExtraNodes_ = rhs.numberExtraNodes_;
5259        continuousObjective_ = rhs.continuousObjective_;
5260        originalContinuousObjective_ = rhs.originalContinuousObjective_;
5261        continuousInfeasibilities_ = rhs.continuousInfeasibilities_;
5262        maximumCutPassesAtRoot_ = rhs.maximumCutPassesAtRoot_;
5263        maximumCutPasses_ = rhs.maximumCutPasses_;
5264        preferredWay_ = rhs.preferredWay_;
5265        currentPassNumber_ = rhs.currentPassNumber_;
5266        memcpy(intParam_, rhs.intParam_, sizeof(intParam_));
5267        memcpy(dblParam_, rhs.dblParam_, sizeof(dblParam_));
5268        globalCuts_ = rhs.globalCuts_;
5269        int i;
5270        for (i = 0; i < numberCutGenerators_; i++) {
5271            delete generator_[i];
5272            delete virginGenerator_[i];
5273        }
5274        delete [] generator_;
5275        delete [] virginGenerator_;
5276        delete [] heuristic_;
5277        maximumWhich_ = rhs.maximumWhich_;
5278        delete [] whichGenerator_;
5279        whichGenerator_ = NULL;
5280        if (maximumWhich_ && rhs.whichGenerator_)
5281            whichGenerator_ = CoinCopyOfArray(rhs.whichGenerator_, maximumWhich_);
5282        maximumRows_ = 0;
5283        currentDepth_ = 0;
5284        randomNumberGenerator_ = rhs.randomNumberGenerator_;
5285        workingBasis_ = CoinWarmStartBasis();
5286        for (i = 0; i < maximumStatistics_; i++)
5287            delete statistics_[i];
5288        delete [] statistics_;
5289        maximumStatistics_ = 0;
5290        statistics_ = NULL;
5291        delete probingInfo_;
5292        probingInfo_ = NULL;
5293        numberFixedAtRoot_ = rhs.numberFixedAtRoot_;
5294        numberFixedNow_ = rhs.numberFixedNow_;
5295        stoppedOnGap_ = rhs.stoppedOnGap_;
5296        eventHappened_ = rhs.eventHappened_;
5297        numberLongStrong_ = rhs.numberLongStrong_;
5298        numberOldActiveCuts_ = rhs.numberOldActiveCuts_;
5299        numberNewCuts_ = rhs.numberNewCuts_;
5300        resolveAfterTakeOffCuts_ = rhs.resolveAfterTakeOffCuts_;
5301        maximumNumberIterations_ = rhs.maximumNumberIterations_;
5302        continuousPriority_ = rhs.continuousPriority_;
5303        numberUpdateItems_ = rhs.numberUpdateItems_;
5304        maximumNumberUpdateItems_ = rhs.maximumNumberUpdateItems_;
5305        delete [] updateItems_;
5306        if (maximumNumberUpdateItems_) {
5307            updateItems_ = new CbcObjectUpdateData [maximumNumberUpdateItems_];
5308            for (i = 0; i < maximumNumberUpdateItems_; i++)
5309                updateItems_[i] = rhs.updateItems_[i];
5310        } else {
5311            updateItems_ = NULL;
5312        }
5313        numberThreads_ = rhs.numberThreads_;
5314        threadMode_ = rhs.threadMode_;
5315        delete master_;
5316        master_ = NULL;
5317        masterThread_ = NULL;
5318        searchStrategy_ = rhs.searchStrategy_;
5319        numberStrongIterations_ = rhs.numberStrongIterations_;
5320        strongInfo_[0] = rhs.strongInfo_[0];
5321        strongInfo_[1] = rhs.strongInfo_[1];
5322        strongInfo_[2] = rhs.strongInfo_[2];
5323        strongInfo_[3] = rhs.strongInfo_[3];
5324        strongInfo_[4] = rhs.strongInfo_[4];
5325        strongInfo_[5] = rhs.strongInfo_[5];
5326        strongInfo_[6] = rhs.strongInfo_[6];
5327        solverCharacteristics_ = NULL;
5328        lastHeuristic_ = NULL;
5329        numberCutGenerators_ = rhs.numberCutGenerators_;
5330        if (numberCutGenerators_) {
5331            generator_ = new CbcCutGenerator * [numberCutGenerators_];
5332            virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_];
5333            int i;
5334            for (i = 0; i < numberCutGenerators_; i++) {
5335                generator_[i] = new CbcCutGenerator(*rhs.generator_[i]);
5336                virginGenerator_[i] = new CbcCutGenerator(*rhs.virginGenerator_[i]);
5337            }
5338        } else {
5339            generator_ = NULL;
5340            virginGenerator_ = NULL;
5341        }
5342        numberHeuristics_ = rhs.numberHeuristics_;
5343        if (numberHeuristics_) {
5344            heuristic_ = new CbcHeuristic * [numberHeuristics_];
5345            memcpy(heuristic_, rhs.heuristic_,
5346                   numberHeuristics_*sizeof(CbcHeuristic *));
5347        } else {
5348            heuristic_ = NULL;
5349        }
5350        lastHeuristic_ = NULL;
5351        if (eventHandler_)
5352            delete eventHandler_ ;
5353        if (rhs.eventHandler_) {
5354            eventHandler_ = rhs.eventHandler_->clone() ;
5355        } else {
5356            eventHandler_ = NULL ;
5357        }
5358        fastNodeDepth_ = rhs.fastNodeDepth_;
5359        if (ownObjects_) {
5360            for (i = 0; i < numberObjects_; i++)
5361                delete object_[i];
5362            delete [] object_;
5363            numberObjects_ = rhs.numberObjects_;
5364            if (numberObjects_) {
5365                object_ = new OsiObject * [numberObjects_];
5366                int i;
5367                for (i = 0; i < numberObjects_; i++)
5368                    object_[i] = (rhs.object_[i])->clone();
5369            } else {
5370                object_ = NULL;
5371            }
5372        } else {
5373            // assume will be redone
5374            numberObjects_ = 0;
5375            object_ = NULL;
5376        }
5377        delete [] originalColumns_;
5378        if (rhs.originalColumns_) {
5379            int numberColumns = rhs.getNumCols();
5380            originalColumns_ = new int [numberColumns];
5381            memcpy(originalColumns_, rhs.originalColumns_, numberColumns*sizeof(int));
5382        } else {
5383            originalColumns_ = NULL;
5384        }
5385        nodeCompare_