source: stable/2.7/Cbc/src/CbcModel.cpp @ 1681

Last change on this file since 1681 was 1681, checked in by forrest, 8 years ago

change from ifdef COIN_HAS_BONMIB to runtime test

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