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

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

allow end cuts and lagomory

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 661.4 KB
Line 
1/* $Id: CbcModel.cpp 1656 2011-05-31 08:10:30Z 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    /*
2395      Grepping through the code, it would appear that this is a command line
2396      debugging hook.  There's no obvious place in the code where this is set to
2397      a negative value.
2398
2399      User hook, says John.
2400    */
2401    if ( intParam_[CbcMaxNumNode] < 0)
2402        eventHappened_ = true; // stop as fast as possible
2403    stoppedOnGap_ = false ;
2404    // See if can stop on gap
2405    bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
2406    double testGap = CoinMax(dblParam_[CbcAllowableGap],
2407                             CoinMax(fabs(bestObjective_), fabs(bestPossibleObjective_))
2408                             * dblParam_[CbcAllowableFractionGap]);
2409    if (bestObjective_ - bestPossibleObjective_ < testGap && getCutoffIncrement() >= 0.0) {
2410        if (bestPossibleObjective_ < getCutoff())
2411            stoppedOnGap_ = true ;
2412        feasible = false;
2413        //eventHappened_=true; // stop as fast as possible
2414    }
2415    /*
2416      Set up for statistics collection, if requested. Standard values are
2417      documented in CbcModel.hpp. The magic number 100 will trigger a dump of
2418      CbcSimpleIntegerDynamicPseudoCost objects (no others). Looks like another
2419      command line debugging hook.
2420    */
2421    statistics_ = NULL;
2422    // Do on switch
2423    if (doStatistics > 0 && doStatistics <= 100) {
2424        maximumStatistics_ = 10000;
2425        statistics_ = new CbcStatistics * [maximumStatistics_];
2426        memset(statistics_, 0, maximumStatistics_*sizeof(CbcStatistics *));
2427    }
2428    // See if we can add integers
2429    if (noObjects && numberIntegers_ < solver_->getNumCols() && (specialOptions_&65536) != 0 && !parentModel_)
2430        AddIntegers();
2431    /*
2432      Do an initial round of cut generation for the root node. Depending on the
2433      type of underlying solver, we may want to do this even if the initial query
2434      to the objects indicates they're satisfied.
2435
2436      solveWithCuts does the heavy lifting. It will iterate a generate/reoptimise
2437      loop (including reduced cost fixing) until no cuts are generated, the
2438      change in objective falls off,  or the limit on the number of rounds of cut
2439      generation is exceeded.
2440
2441      At the end of all this, any cuts will be recorded in cuts and also
2442      installed in the solver's constraint system. We'll have reoptimised, and
2443      removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
2444      adjusted accordingly).
2445    */
2446    int iObject ;
2447    int preferredWay ;
2448    int numberUnsatisfied = 0 ;
2449    delete [] currentSolution_;
2450    currentSolution_ = new double [numberColumns];
2451    testSolution_ = currentSolution_;
2452    memcpy(currentSolution_, solver_->getColSolution(),
2453           numberColumns*sizeof(double)) ;
2454    // point to useful information
2455    OsiBranchingInformation usefulInfo = usefulInformation();
2456
2457    for (iObject = 0 ; iObject < numberObjects_ ; iObject++) {
2458        double infeasibility =
2459            object_[iObject]->infeasibility(&usefulInfo, preferredWay) ;
2460        if (infeasibility ) numberUnsatisfied++ ;
2461    }
2462    // replace solverType
2463    if (solverCharacteristics_->tryCuts())  {
2464
2465        if (numberUnsatisfied)   {
2466            // User event
2467            if (!eventHappened_ && feasible) {
2468                feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
2469                                         NULL);
2470                if ((specialOptions_&524288) != 0 && !parentModel_
2471                        && storedRowCuts_) {
2472                    if (feasible) {
2473                        /* pick up stuff and try again
2474                        add cuts, maybe keep around
2475                        do best solution and if so new heuristics
2476                        obviously tighten bounds
2477                        */
2478                        // A and B probably done on entry
2479                        // A) tight bounds on integer variables
2480                        const double * lower = solver_->getColLower();
2481                        const double * upper = solver_->getColUpper();
2482                        const double * tightLower = storedRowCuts_->tightLower();
2483                        const double * tightUpper = storedRowCuts_->tightUpper();
2484                        int nTightened = 0;
2485                        for (int i = 0; i < numberIntegers_; i++) {
2486                            int iColumn = integerVariable_[i];
2487                            if (tightLower[iColumn] > lower[iColumn]) {
2488                                nTightened++;
2489                                solver_->setColLower(iColumn, tightLower[iColumn]);
2490                            }
2491                            if (tightUpper[iColumn] < upper[iColumn]) {
2492                                nTightened++;
2493                                solver_->setColUpper(iColumn, tightUpper[iColumn]);
2494                            }
2495                        }
2496                        if (nTightened)
2497                          COIN_DETAIL_PRINT(printf("%d tightened by alternate cuts\n", nTightened));
2498                        if (storedRowCuts_->bestObjective() < bestObjective_) {
2499                            // B) best solution
2500                            double objValue = storedRowCuts_->bestObjective();
2501                            setBestSolution(CBC_SOLUTION, objValue,
2502                                            storedRowCuts_->bestSolution()) ;
2503                            // Do heuristics
2504                            // Allow RINS
2505                            for (int i = 0; i < numberHeuristics_; i++) {
2506                                CbcHeuristicRINS * rins
2507                                = dynamic_cast<CbcHeuristicRINS *> (heuristic_[i]);
2508                                if (rins) {
2509                                    rins->setLastNode(-100);
2510                                }
2511                            }
2512                            doHeuristicsAtRoot();
2513                        }
2514#ifdef JJF_ZERO
2515                        int nCuts = storedRowCuts_->sizeRowCuts();
2516                        // add to global list
2517                        for (int i = 0; i < nCuts; i++) {
2518                            OsiRowCut newCut(*storedRowCuts_->rowCutPointer(i));
2519                            newCut.setGloballyValidAsInteger(2);
2520                            newCut.mutableRow().setTestForDuplicateIndex(false);
2521                            globalCuts_.insert(newCut) ;
2522                        }
2523#else
2524                        addCutGenerator(storedRowCuts_, -99, "Stored from previous run",
2525                                        true, false, false, -200);
2526#endif
2527                        // Set cuts as active
2528                        delete [] addedCuts_ ;
2529                        maximumNumberCuts_ = cuts.sizeRowCuts();
2530                        if (maximumNumberCuts_) {
2531                            addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
2532                        } else {
2533                            addedCuts_ = NULL;
2534                        }
2535                        for (int i = 0; i < maximumNumberCuts_; i++)
2536                            addedCuts_[i] = new CbcCountRowCut(*cuts.rowCutPtr(i),
2537                                                               NULL, -1, -1, 2);
2538                        COIN_DETAIL_PRINT(printf("size %d\n", cuts.sizeRowCuts()));
2539                        cuts = OsiCuts();
2540                        currentNumberCuts_ = maximumNumberCuts_;
2541                        feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
2542                                                 NULL);
2543                        for (int i = 0; i < maximumNumberCuts_; i++)
2544                            delete addedCuts_[i];
2545                    }
2546                    delete storedRowCuts_;
2547                    storedRowCuts_ = NULL;
2548                }
2549            } else {
2550                feasible = false;
2551            }
2552        }       else if (solverCharacteristics_->solutionAddsCuts() ||
2553                   solverCharacteristics_->alwaysTryCutsAtRootNode()) {
2554            // may generate cuts and turn the solution
2555            //to an infeasible one
2556            feasible = solveWithCuts(cuts, 1,
2557                                     NULL);
2558        }
2559    }
2560    // check extra info on feasibility
2561    if (!solverCharacteristics_->mipFeasible())
2562        feasible = false;
2563
2564#ifdef CLP_RESOLVE
2565    if ((moreSpecialOptions_&2097152)!=0&&!parentModel_&&feasible) {
2566      OsiClpSolverInterface * clpSolver
2567        = dynamic_cast<OsiClpSolverInterface *> (solver_);
2568      if (clpSolver)
2569        resolveClp(clpSolver,0);
2570    }
2571#endif
2572    // make cut generators less aggressive
2573    for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
2574        CglCutGenerator * generator = generator_[iCutGenerator]->generator();
2575        generator->setAggressiveness(generator->getAggressiveness() - 100);
2576    }
2577    currentNumberCuts_ = numberNewCuts_ ;
2578    // See if can stop on gap
2579    bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
2580    testGap = CoinMax(dblParam_[CbcAllowableGap],
2581                      CoinMax(fabs(bestObjective_), fabs(bestPossibleObjective_))
2582                      * dblParam_[CbcAllowableFractionGap]);
2583    if (bestObjective_ - bestPossibleObjective_ < testGap && getCutoffIncrement() >= 0.0) {
2584        if (bestPossibleObjective_ < getCutoff())
2585            stoppedOnGap_ = true ;
2586        feasible = false;
2587    }
2588    // User event
2589    if (eventHappened_)
2590        feasible = false;
2591#if defined(COIN_HAS_CLP)&&defined(COIN_HAS_CPX)
2592    /*
2593      This is the notion of using Cbc stuff to get going, then calling cplex to
2594      finish off.
2595    */
2596    if (feasible && (specialOptions_&16384) != 0 && fastNodeDepth_ == -2 && !parentModel_) {
2597        // Use Cplex to do search!
2598        double time1 = CoinCpuTime();
2599        OsiClpSolverInterface * clpSolver
2600        = dynamic_cast<OsiClpSolverInterface *> (solver_);
2601        OsiCpxSolverInterface cpxSolver;
2602        double direction = clpSolver->getObjSense();
2603        cpxSolver.setObjSense(direction);
2604        // load up cplex
2605        const CoinPackedMatrix * matrix = continuousSolver_->getMatrixByCol();
2606        const double * rowLower = continuousSolver_->getRowLower();
2607        const double * rowUpper = continuousSolver_->getRowUpper();
2608        const double * columnLower = continuousSolver_->getColLower();
2609        const double * columnUpper = continuousSolver_->getColUpper();
2610        const double * objective = continuousSolver_->getObjCoefficients();
2611        cpxSolver.loadProblem(*matrix, columnLower, columnUpper,
2612                              objective, rowLower, rowUpper);
2613        double * setSol = new double [numberIntegers_];
2614        int * setVar = new int [numberIntegers_];
2615        // cplex doesn't know about objective offset
2616        double offset = clpSolver->getModelPtr()->objectiveOffset();
2617        for (int i = 0; i < numberIntegers_; i++) {
2618            int iColumn = integerVariable_[i];
2619            cpxSolver.setInteger(iColumn);
2620            if (bestSolution_) {
2621                setSol[i] = bestSolution_[iColumn];
2622                setVar[i] = iColumn;
2623            }
2624        }
2625        CPXENVptr env = cpxSolver.getEnvironmentPtr();
2626        CPXLPptr lpPtr = cpxSolver.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL);
2627        cpxSolver.switchToMIP();
2628        if (bestSolution_) {
2629            CPXcopymipstart(env, lpPtr, numberIntegers_, setVar, setSol);
2630        }
2631        if (clpSolver->getNumRows() > continuousSolver_->getNumRows() && false) {
2632            // add cuts
2633            const CoinPackedMatrix * matrix = clpSolver->getMatrixByRow();
2634            const double * rhs = clpSolver->getRightHandSide();
2635            const char * rowSense = clpSolver->getRowSense();
2636            const double * elementByRow = matrix->getElements();
2637            const int * column = matrix->getIndices();
2638            const CoinBigIndex * rowStart = matrix->getVectorStarts();
2639            const int * rowLength = matrix->getVectorLengths();
2640            int nStart = continuousSolver_->getNumRows();
2641            int nRows = clpSolver->getNumRows();
2642            int size = rowStart[nRows-1] + rowLength[nRows-1] -
2643                       rowStart[nStart];
2644            int nAdd = 0;
2645            double * rmatval = new double [size];
2646            int * rmatind = new int [size];
2647            int * rmatbeg = new int [nRows-nStart+1];
2648            size = 0;
2649            rmatbeg[0] = 0;
2650            for (int i = nStart; i < nRows; i++) {
2651                for (int k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2652                    rmatind[size] = column[k];
2653                    rmatval[size++] = elementByRow[k];
2654                }
2655                nAdd++;
2656                rmatbeg[nAdd] = size;
2657            }
2658            CPXaddlazyconstraints(env, lpPtr, nAdd, size,
2659                                  rhs, rowSense, rmatbeg,
2660                                  rmatind, rmatval, NULL);
2661            CPXsetintparam( env, CPX_PARAM_REDUCE,
2662                            // CPX_PREREDUCE_NOPRIMALORDUAL (0)
2663                            CPX_PREREDUCE_PRIMALONLY);
2664        }
2665        if (getCutoff() < 1.0e50) {
2666            double useCutoff = getCutoff() + offset;
2667            if (bestObjective_ < 1.0e50)
2668                useCutoff = bestObjective_ + offset + 1.0e-7;
2669            cpxSolver.setDblParam(OsiDualObjectiveLimit, useCutoff*
2670                                  direction);
2671            if ( direction > 0.0 )
2672                CPXsetdblparam( env, CPX_PARAM_CUTUP, useCutoff ) ; // min
2673            else
2674                CPXsetdblparam( env, CPX_PARAM_CUTLO, useCutoff ) ; // max
2675        }
2676        CPXsetdblparam(env, CPX_PARAM_EPGAP, dblParam_[CbcAllowableFractionGap]);
2677        delete [] setSol;
2678        delete [] setVar;
2679        char printBuffer[200];
2680        if (offset) {
2681            sprintf(printBuffer, "Add %g to all Cplex messages for true objective",
2682                    -offset);
2683            messageHandler()->message(CBC_GENERAL, messages())
2684            << printBuffer << CoinMessageEol ;
2685            cpxSolver.setDblParam(OsiObjOffset, offset);
2686        }
2687        cpxSolver.branchAndBound();
2688        double timeTaken = CoinCpuTime() - time1;
2689        sprintf(printBuffer, "Cplex took %g seconds",
2690                timeTaken);
2691        messageHandler()->message(CBC_GENERAL, messages())
2692        << printBuffer << CoinMessageEol ;
2693        numberExtraNodes_ = CPXgetnodecnt(env, lpPtr);
2694        numberExtraIterations_ = CPXgetmipitcnt(env, lpPtr);
2695        double value = cpxSolver.getObjValue() * direction;
2696        if (cpxSolver.isProvenOptimal() && value <= getCutoff()) {
2697            feasible = true;
2698            clpSolver->setWarmStart(NULL);
2699            // try and do solution
2700            double * newSolution =
2701                CoinCopyOfArray(cpxSolver.getColSolution(),
2702                                getNumCols());
2703            setBestSolution(CBC_STRONGSOL, value, newSolution) ;
2704            delete [] newSolution;
2705        }
2706        feasible = false;
2707    }
2708#endif
2709    /*
2710      A hook to use clp to quickly explore some part of the tree.
2711    */
2712    if (fastNodeDepth_ == 1000 &&/*!parentModel_*/(specialOptions_&2048) == 0) {
2713        fastNodeDepth_ = -1;
2714        CbcObject * obj =
2715            new CbcFollowOn(this);
2716        obj->setPriority(1);
2717        addObjects(1, &obj);
2718    }
2719    int saveNumberSolves = numberSolves_;
2720    int saveNumberIterations = numberIterations_;
2721    if (fastNodeDepth_ >= 0 &&/*!parentModel_*/(specialOptions_&2048) == 0) {
2722        // add in a general depth object doClp
2723        int type = (fastNodeDepth_ <= 100) ? fastNodeDepth_ : -(fastNodeDepth_ - 100);
2724        CbcObject * obj =
2725            new CbcGeneralDepth(this, type);
2726        addObjects(1, &obj);
2727        // mark as done
2728        fastNodeDepth_ += 1000000;
2729        delete obj;
2730        // fake number of objects
2731        numberObjects_--;
2732        if (parallelMode() < -1) {
2733            // But make sure position is correct
2734            OsiObject * obj2 = object_[numberObjects_];
2735            obj = dynamic_cast<CbcObject *> (obj2);
2736            assert (obj);
2737            obj->setPosition(numberObjects_);
2738        }
2739    }
2740#ifdef COIN_HAS_CLP
2741#ifdef NO_CRUNCH
2742    if (true) {
2743        OsiClpSolverInterface * clpSolver
2744        = dynamic_cast<OsiClpSolverInterface *> (solver_);
2745        if (clpSolver && !parentModel_) {
2746            ClpSimplex * clpSimplex = clpSolver->getModelPtr();
2747            clpSimplex->setSpecialOptions(clpSimplex->specialOptions() | 131072);
2748            //clpSimplex->startPermanentArrays();
2749            clpSimplex->setPersistenceFlag(2);
2750        }
2751    }
2752#endif
2753#endif
2754// Save copy of solver
2755    OsiSolverInterface * saveSolver = NULL;
2756    if (!parentModel_ && (specialOptions_&(512 + 32768)) != 0)
2757        saveSolver = solver_->clone();
2758    double checkCutoffForRestart = 1.0e100;
2759    saveModel(saveSolver, &checkCutoffForRestart, &feasible);
2760    if ((specialOptions_&262144) != 0 && !parentModel_) {
2761        // Save stuff and return!
2762        storedRowCuts_->saveStuff(bestObjective_, bestSolution_,
2763                                  solver_->getColLower(),
2764                                  solver_->getColUpper());
2765        delete [] lowerBefore;
2766        delete [] upperBefore;
2767        delete saveSolver;
2768        return;
2769    }
2770    /*
2771      We've taken the continuous relaxation as far as we can. Time to branch.
2772      The first order of business is to actually create a node. chooseBranch
2773      currently uses strong branching to evaluate branch object candidates,
2774      unless forced back to simple branching. If chooseBranch concludes that a
2775      branching candidate is monotone (anyAction == -1) or infeasible (anyAction
2776      == -2) when forced to integer values, it returns here immediately.
2777
2778      Monotone variables trigger a call to resolve(). If the problem remains
2779      feasible, try again to choose a branching variable. At the end of the loop,
2780      resolved == true indicates that some variables were fixed.
2781
2782      Loss of feasibility will result in the deletion of newNode.
2783    */
2784
2785    bool resolved = false ;
2786    CbcNode *newNode = NULL ;
2787    numberFixedAtRoot_ = 0;
2788    numberFixedNow_ = 0;
2789    int numberIterationsAtContinuous = numberIterations_;
2790    //solverCharacteristics_->setSolver(solver_);
2791    if (feasible) {
2792        //#define HOTSTART -1
2793#if HOTSTART<0
2794        if (bestSolution_ && !parentModel_ && !hotstartSolution_ &&
2795                (moreSpecialOptions_&1024) != 0) {
2796            // Set priorities so only branch on ones we need to
2797            // use djs and see if only few branches needed
2798#ifndef NDEBUG
2799            double integerTolerance = getIntegerTolerance() ;
2800#endif
2801            bool possible = true;
2802            const double * saveLower = continuousSolver_->getColLower();
2803            const double * saveUpper = continuousSolver_->getColUpper();
2804            for (int i = 0; i < numberObjects_; i++) {
2805                const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object_[i]);
2806                if (thisOne) {
2807                    int iColumn = thisOne->columnNumber();
2808                    if (saveUpper[iColumn] > saveLower[iColumn] + 1.5) {
2809                        possible = false;
2810                        break;
2811                    }
2812                } else {
2813                    possible = false;
2814                    break;
2815                }
2816            }
2817            if (possible) {
2818                OsiSolverInterface * solver = continuousSolver_->clone();
2819                int numberColumns = solver->getNumCols();
2820                for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
2821                    double value = bestSolution_[iColumn] ;
2822                    value = CoinMax(value, saveLower[iColumn]) ;
2823                    value = CoinMin(value, saveUpper[iColumn]) ;
2824                    value = floor(value + 0.5);
2825                    if (solver->isInteger(iColumn)) {
2826                        solver->setColLower(iColumn, value);
2827                        solver->setColUpper(iColumn, value);
2828                    }
2829                }
2830                solver->setHintParam(OsiDoDualInResolve, false, OsiHintTry);
2831                // objlim and all slack
2832                double direction = solver->getObjSense();
2833                solver->setDblParam(OsiDualObjectiveLimit, 1.0e50*direction);
2834                CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis *> (solver->getEmptyWarmStart());
2835                solver->setWarmStart(basis);
2836                delete basis;
2837                bool changed = true;
2838                hotstartPriorities_ = new int [numberColumns];
2839                for (int iColumn = 0; iColumn < numberColumns; iColumn++)
2840                    hotstartPriorities_[iColumn] = 1;
2841                while (changed) {
2842                    changed = false;
2843                    solver->resolve();
2844                    if (!solver->isProvenOptimal()) {
2845                        possible = false;
2846                        break;
2847                    }
2848                    const double * dj = solver->getReducedCost();
2849                    const double * colLower = solver->getColLower();
2850                    const double * colUpper = solver->getColUpper();
2851                    const double * solution = solver->getColSolution();
2852                    int nAtLbNatural = 0;
2853                    int nAtUbNatural = 0;
2854                    int nZeroDj = 0;
2855                    int nForced = 0;
2856                    for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
2857                        double value = solution[iColumn] ;
2858                        value = CoinMax(value, saveLower[iColumn]) ;
2859                        value = CoinMin(value, saveUpper[iColumn]) ;
2860                        if (solver->isInteger(iColumn)) {
2861                            assert(fabs(value - solution[iColumn]) <= integerTolerance) ;
2862                            if (hotstartPriorities_[iColumn] == 1) {
2863                                if (dj[iColumn] < -1.0e-6) {
2864                                    // negative dj
2865                                    if (saveUpper[iColumn] == colUpper[iColumn]) {
2866                                        nAtUbNatural++;
2867                                        hotstartPriorities_[iColumn] = 2;
2868                                        solver->setColLower(iColumn, saveLower[iColumn]);
2869                                        solver->setColUpper(iColumn, saveUpper[iColumn]);
2870                                    } else {
2871                                        nForced++;
2872                                    }
2873                                } else if (dj[iColumn] > 1.0e-6) {
2874                                    // positive dj
2875                                    if (saveLower[iColumn] == colLower[iColumn]) {
2876                                        nAtLbNatural++;
2877                                        hotstartPriorities_[iColumn] = 2;
2878                                        solver->setColLower(iColumn, saveLower[iColumn]);
2879                                        solver->setColUpper(iColumn, saveUpper[iColumn]);
2880                                    } else {
2881                                        nForced++;
2882                                    }
2883                                } else {
2884                                    // zero dj
2885                                    nZeroDj++;
2886                                }
2887                            }
2888                        }
2889                    }
2890#ifdef CLP_INVESTIGATE
2891                    printf("%d forced, %d naturally at lower, %d at upper - %d zero dj\n",
2892                           nForced, nAtLbNatural, nAtUbNatural, nZeroDj);
2893#endif
2894                    if (nAtLbNatural || nAtUbNatural) {
2895                        changed = true;
2896                    } else {
2897                        if (nForced + nZeroDj > 50 ||
2898                                (nForced + nZeroDj)*4 > numberIntegers_)
2899                            possible = false;
2900                    }
2901                }
2902                delete solver;
2903            }
2904            if (possible) {
2905                setHotstartSolution(bestSolution_);
2906                if (!saveCompare) {
2907                    // create depth first comparison
2908                    saveCompare = nodeCompare_;
2909                    // depth first
2910                    nodeCompare_ = new CbcCompareDepth();
2911                    tree_->setComparison(*nodeCompare_) ;
2912                }
2913            } else {
2914                delete [] hotstartPriorities_;
2915                hotstartPriorities_ = NULL;
2916            }
2917        }
2918#endif
2919#if HOTSTART>0
2920        if (hotstartSolution_ && !hotstartPriorities_) {
2921            // Set up hot start
2922            OsiSolverInterface * solver = solver_->clone();
2923            double direction = solver_->getObjSense() ;
2924            int numberColumns = solver->getNumCols();
2925            double * saveLower = CoinCopyOfArray(solver->getColLower(), numberColumns);
2926            double * saveUpper = CoinCopyOfArray(solver->getColUpper(), numberColumns);
2927            // move solution
2928            solver->setColSolution(hotstartSolution_);
2929            // point to useful information
2930            const double * saveSolution = testSolution_;
2931            testSolution_ = solver->getColSolution();
2932            OsiBranchingInformation usefulInfo = usefulInformation();
2933            testSolution_ = saveSolution;
2934            /*
2935            Run through the objects and use feasibleRegion() to set variable bounds
2936            so as to fix the variables specified in the objects at their value in this
2937            solution. Since the object list contains (at least) one object for every
2938            integer variable, this has the effect of fixing all integer variables.
2939            */
2940            for (int i = 0; i < numberObjects_; i++)
2941                object_[i]->feasibleRegion(solver, &usefulInfo);
2942            solver->resolve();
2943            assert (solver->isProvenOptimal());
2944            double gap = CoinMax((solver->getObjValue() - solver_->getObjValue()) * direction, 0.0) ;
2945            const double * dj = solver->getReducedCost();
2946            const double * colLower = solver->getColLower();
2947            const double * colUpper = solver->getColUpper();
2948            const double * solution = solver->getColSolution();
2949            int nAtLbNatural = 0;
2950            int nAtUbNatural = 0;
2951            int nAtLbNaturalZero = 0;
2952            int nAtUbNaturalZero = 0;
2953            int nAtLbFixed = 0;
2954            int nAtUbFixed = 0;
2955            int nAtOther = 0;
2956            int nAtOtherNatural = 0;
2957            int nNotNeeded = 0;
2958            delete [] hotstartSolution_;
2959            hotstartSolution_ = new double [numberColumns];
2960            delete [] hotstartPriorities_;
2961            hotstartPriorities_ = new int [numberColumns];
2962            int * order = (int *) saveUpper;
2963            int nFix = 0;
2964            double bestRatio = COIN_DBL_MAX;
2965            for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
2966                double value = solution[iColumn] ;
2967                value = CoinMax(value, saveLower[iColumn]) ;
2968                value = CoinMin(value, saveUpper[iColumn]) ;
2969                double sortValue = COIN_DBL_MAX;
2970                if (solver->isInteger(iColumn)) {
2971                    assert(fabs(value - solution[iColumn]) <= 1.0e-5) ;
2972                    double value2 = floor(value + 0.5);
2973                    if (dj[iColumn] < -1.0e-6) {
2974                        // negative dj
2975                        //assert (value2==colUpper[iColumn]);
2976                        if (saveUpper[iColumn] == colUpper[iColumn]) {
2977                            nAtUbNatural++;
2978                            sortValue = 0.0;
2979                            double value = -dj[iColumn];
2980                            if (value > gap)
2981                                nFix++;
2982                            else if (gap < value*bestRatio)
2983                                bestRatio = gap / value;
2984                            if (saveLower[iColumn] != colLower[iColumn]) {
2985                                nNotNeeded++;
2986                                sortValue = 1.0e20;
2987                            }
2988                        } else if (saveLower[iColumn] == colUpper[iColumn]) {
2989                            nAtLbFixed++;
2990                            sortValue = dj[iColumn];
2991                        } else {
2992                            nAtOther++;
2993                            sortValue = 0.0;
2994                            if (saveLower[iColumn] != colLower[iColumn] &&
2995                                    saveUpper[iColumn] != colUpper[iColumn]) {
2996                                nNotNeeded++;
2997                                sortValue = 1.0e20;
2998                            }
2999                        }
3000                    } else if (dj[iColumn] > 1.0e-6) {
3001                        // positive dj
3002                        //assert (value2==colLower[iColumn]);
3003                        if (saveLower[iColumn] == colLower[iColumn]) {
3004                            nAtLbNatural++;
3005                            sortValue = 0.0;
3006                            double value = dj[iColumn];
3007                            if (value > gap)
3008                                nFix++;
3009                            else if (gap < value*bestRatio)
3010                                bestRatio = gap / value;
3011                            if (saveUpper[iColumn] != colUpper[iColumn]) {
3012                                nNotNeeded++;
3013                                sortValue = 1.0e20;
3014                            }
3015                        } else if (saveUpper[iColumn] == colLower[iColumn]) {
3016                            nAtUbFixed++;
3017                            sortValue = -dj[iColumn];
3018                        } else {
3019                            nAtOther++;
3020                            sortValue = 0.0;
3021                            if (saveLower[iColumn] != colLower[iColumn] &&
3022                                    saveUpper[iColumn] != colUpper[iColumn]) {
3023                                nNotNeeded++;
3024                                sortValue = 1.0e20;
3025                            }
3026                        }
3027                    } else {
3028                        // zero dj
3029                        if (value2 == saveUpper[iColumn]) {
3030                            nAtUbNaturalZero++;
3031                            sortValue = 0.0;
3032                            if (saveLower[iColumn] != colLower[iColumn]) {
3033                                nNotNeeded++;
3034                                sortValue = 1.0e20;
3035                            }
3036                        } else if (value2 == saveLower[iColumn]) {
3037                            nAtLbNaturalZero++;
3038                            sortValue = 0.0;
3039                        } else {
3040                            nAtOtherNatural++;
3041                            sortValue = 0.0;
3042                            if (saveLower[iColumn] != colLower[iColumn] &&
3043                                    saveUpper[iColumn] != colUpper[iColumn]) {
3044                                nNotNeeded++;
3045                                sortValue = 1.0e20;
3046                            }
3047                        }
3048                    }
3049#if HOTSTART==3
3050                    sortValue = -fabs(dj[iColumn]);
3051#endif
3052                }
3053                hotstartSolution_[iColumn] = value ;
3054                saveLower[iColumn] = sortValue;
3055                order[iColumn] = iColumn;
3056            }
3057            COIN_DETAIL_PRINT(printf("** can fix %d columns - best ratio for others is %g on gap of %g\n",
3058                                     nFix, bestRatio, gap));
3059            int nNeg = 0;
3060            CoinSort_2(saveLower, saveLower + numberColumns, order);
3061            for (int i = 0; i < numberColumns; i++) {
3062                if (saveLower[i] < 0.0) {
3063                    nNeg++;
3064#if HOTSTART==2||HOTSTART==3
3065                    // swap sign ?
3066                    saveLower[i] = -saveLower[i];
3067#endif
3068                }
3069            }
3070            CoinSort_2(saveLower, saveLower + nNeg, order);
3071            for (int i = 0; i < numberColumns; i++) {
3072#if HOTSTART==1
3073                hotstartPriorities_[order[i]] = 100;
3074#else
3075                hotstartPriorities_[order[i]] = -(i + 1);
3076#endif
3077            }
3078            COIN_DETAIL_PRINT(printf("nAtLbNat %d,nAtUbNat %d,nAtLbNatZero %d,nAtUbNatZero %d,nAtLbFixed %d,nAtUbFixed %d,nAtOther %d,nAtOtherNat %d, useless %d %d\n",
3079                   nAtLbNatural,
3080                   nAtUbNatural,
3081                   nAtLbNaturalZero,
3082                   nAtUbNaturalZero,
3083                   nAtLbFixed,
3084                   nAtUbFixed,
3085                   nAtOther,
3086                                     nAtOtherNatural, nNotNeeded, nNeg));
3087            delete [] saveLower;
3088            delete [] saveUpper;
3089            if (!saveCompare) {
3090                // create depth first comparison
3091                saveCompare = nodeCompare_;
3092                // depth first
3093                nodeCompare_ = new CbcCompareDepth();
3094                tree_->setComparison(*nodeCompare_) ;
3095            }
3096        }
3097#endif
3098        newNode = new CbcNode ;
3099        // Set objective value (not so obvious if NLP etc)
3100        setObjectiveValue(newNode, NULL);
3101        anyAction = -1 ;
3102        // To make depth available we may need a fake node
3103        CbcNode fakeNode;
3104        if (!currentNode_) {
3105            // Not true if sub trees assert (!numberNodes_);
3106            currentNode_ = &fakeNode;
3107        }
3108        phase_ = 3;
3109        // only allow 1000 passes
3110        int numberPassesLeft = 1000;
3111        // This is first crude step
3112        if (numberAnalyzeIterations_) {
3113            delete [] analyzeResults_;
3114            analyzeResults_ = new double [4*numberIntegers_];
3115            numberFixedAtRoot_ = newNode->analyze(this, analyzeResults_);
3116            if (numberFixedAtRoot_ > 0) {
3117              COIN_DETAIL_PRINT(printf("%d fixed by analysis\n", numberFixedAtRoot_));
3118                setPointers(solver_);
3119                numberFixedNow_ = numberFixedAtRoot_;
3120            } else if (numberFixedAtRoot_ < 0) {
3121              COIN_DETAIL_PRINT(printf("analysis found to be infeasible\n"));
3122                anyAction = -2;
3123                delete newNode ;
3124                newNode = NULL ;
3125                feasible = false ;
3126            }
3127        }
3128        OsiSolverBranch * branches = NULL;
3129        anyAction = chooseBranch(newNode, numberPassesLeft, NULL, cuts, resolved,
3130                                 NULL, NULL, NULL, branches);
3131        if (anyAction == -2 || newNode->objectiveValue() >= cutoff) {
3132            if (anyAction != -2) {
3133                // zap parent nodeInfo
3134#ifdef COIN_DEVELOP
3135                printf("zapping CbcNodeInfo %x\n", reinterpret_cast<int>(newNode->nodeInfo()->parent()));
3136#endif
3137                if (newNode->nodeInfo())
3138                    newNode->nodeInfo()->nullParent();
3139            }
3140            delete newNode ;
3141            newNode = NULL ;
3142            feasible = false ;
3143        }
3144    }
3145    if (newNode && probingInfo_) {
3146        int number01 = probingInfo_->numberIntegers();
3147        //const fixEntry * entry = probingInfo_->fixEntries();
3148        const int * toZero = probingInfo_->toZero();
3149        //const int * toOne = probingInfo_->toOne();
3150        //const int * integerVariable = probingInfo_->integerVariable();
3151        if (toZero[number01]) {
3152            CglTreeProbingInfo info(*probingInfo_);
3153#ifdef JJF_ZERO
3154            /*
3155              Marginal idea. Further exploration probably good. Build some extra
3156              cliques from probing info. Not quite worth the effort?
3157            */
3158            OsiSolverInterface * fake = info.analyze(*solver_, 1);
3159            if (fake) {
3160                fake->writeMps("fake");
3161                CglFakeClique cliqueGen(fake);
3162                cliqueGen.setStarCliqueReport(false);
3163                cliqueGen.setRowCliqueReport(false);
3164                cliqueGen.setMinViolation(0.1);
3165                addCutGenerator(&cliqueGen, 1, "Fake cliques", true, false, false, -200);
3166                generator_[numberCutGenerators_-1]->setTiming(true);
3167            }
3168#endif
3169            if (probingInfo_->packDown()) {
3170#ifdef CLP_INVESTIGATE
3171                printf("%d implications on %d 0-1\n", toZero[number01], number01);
3172#endif
3173                // Create a cut generator that remembers implications discovered at root.
3174                CglImplication implication(probingInfo_);
3175                addCutGenerator(&implication, 1, "ImplicationCuts", true, false, false, -200);
3176            } else {
3177                delete probingInfo_;
3178                probingInfo_ = NULL;
3179            }
3180        } else {
3181            delete probingInfo_;
3182
3183            probingInfo_ = NULL;
3184        }
3185    }
3186    /*
3187      At this point, the root subproblem is infeasible or fathomed by bound
3188      (newNode == NULL), or we're live with an objective value that satisfies the
3189      current objective cutoff.
3190    */
3191    assert (!newNode || newNode->objectiveValue() <= cutoff) ;
3192    // Save address of root node as we don't want to delete it
3193    // initialize for print out
3194    int lastDepth = 0;
3195    int lastUnsatisfied = 0;
3196    if (newNode)
3197        lastUnsatisfied = newNode->numberUnsatisfied();
3198    /*
3199      The common case is that the lp relaxation is feasible but doesn't satisfy
3200      integrality (i.e., newNode->branchingObject(), indicating we've been able to
3201      select a branching variable). Remove any cuts that have gone slack due to
3202      forcing monotone variables. Then tack on an CbcFullNodeInfo object and full
3203      basis (via createInfo()) and stash the new cuts in the nodeInfo (via
3204      addCuts()). If, by some miracle, we have an integral solution at the root
3205      (newNode->branchingObject() is NULL), takeOffCuts() will ensure that the solver holds
3206      a valid solution for use by setBestSolution().
3207    */
3208    CoinWarmStartBasis *lastws = NULL ;
3209    if (feasible && newNode->branchingObject()) {
3210        if (resolved) {
3211            takeOffCuts(cuts, false, NULL) ;
3212#     ifdef CHECK_CUT_COUNTS
3213            { printf("Number of rows after chooseBranch fix (root)"
3214                         "(active only) %d\n",
3215                         numberRowsAtContinuous_ + numberNewCuts_ + numberOldActiveCuts_) ;
3216                const CoinWarmStartBasis* debugws =
3217                    dynamic_cast <const CoinWarmStartBasis*>(solver_->getWarmStart()) ;
3218                debugws->print() ;
3219                delete debugws ;
3220            }
3221#     endif
3222        }
3223        //newNode->createInfo(this,NULL,NULL,NULL,NULL,0,0) ;
3224        //newNode->nodeInfo()->addCuts(cuts,
3225        //                       newNode->numberBranches(),whichGenerator_) ;
3226        if (lastws) delete lastws ;
3227        lastws = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
3228    }
3229    /*
3230      Continuous data to be used later
3231    */
3232    continuousObjective_ = solver_->getObjValue() * solver_->getObjSense();
3233    continuousInfeasibilities_ = 0 ;
3234    if (newNode) {
3235        continuousObjective_ = newNode->objectiveValue() ;
3236        delete [] continuousSolution_;
3237        continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
3238                                              numberColumns);
3239        continuousInfeasibilities_ = newNode->numberUnsatisfied() ;
3240    }
3241    /*
3242      Bound may have changed so reset in objects
3243    */
3244    {
3245        int i ;
3246        for (i = 0; i < numberObjects_; i++)
3247            object_[i]->resetBounds(solver_) ;
3248    }
3249    /*
3250      Feasible? Then we should have either a live node prepped for future
3251      expansion (indicated by variable() >= 0), or (miracle of miracles) an
3252      integral solution at the root node.
3253
3254      initializeInfo sets the reference counts in the nodeInfo object.  Since
3255      this node is still live, push it onto the heap that holds the live set.
3256    */
3257    double bestValue = 0.0 ;
3258    if (newNode) {
3259        bestValue = newNode->objectiveValue();
3260        if (newNode->branchingObject()) {
3261            newNode->initializeInfo() ;
3262            tree_->push(newNode) ;
3263            if (statistics_) {
3264                if (numberNodes2_ == maximumStatistics_) {
3265                    maximumStatistics_ = 2 * maximumStatistics_;
3266                    CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_];
3267                    memset(temp, 0, maximumStatistics_*sizeof(CbcStatistics *));
3268                    memcpy(temp, statistics_, numberNodes2_*sizeof(CbcStatistics *));
3269                    delete [] statistics_;
3270                    statistics_ = temp;
3271                }
3272                assert (!statistics_[numberNodes2_]);
3273                statistics_[numberNodes2_] = new CbcStatistics(newNode, this);
3274            }
3275            numberNodes2_++;
3276#     ifdef CHECK_NODE
3277            printf("Node %x on tree\n", newNode) ;
3278#     endif
3279        } else {
3280            // continuous is integer
3281            double objectiveValue = newNode->objectiveValue();
3282            setBestSolution(CBC_SOLUTION, objectiveValue,
3283                            solver_->getColSolution()) ;
3284            delete newNode ;
3285            newNode = NULL ;
3286        }
3287    }
3288
3289    if (printFrequency_ <= 0) {
3290        printFrequency_ = 1000 ;
3291        if (getNumCols() > 2000)
3292            printFrequency_ = 100 ;
3293    }
3294    /*
3295      It is possible that strong branching fixes one variable and then the code goes round
3296      again and again.  This can take too long.  So we need to warn user - just once.
3297    */
3298    numberLongStrong_ = 0;
3299    CbcNode * createdNode = NULL;
3300#ifdef CBC_THREAD
3301    if ((specialOptions_&2048) != 0)
3302        numberThreads_ = 0;
3303    if (numberThreads_ ) {
3304        nodeCompare_->sayThreaded(); // need to use addresses
3305        master_ = new CbcBaseModel(*this,
3306                                   (parallelMode() < -1) ? 1 : 0);
3307        masterThread_ = master_->masterThread();
3308    }
3309#endif
3310#ifdef COIN_HAS_CLP
3311    {
3312        OsiClpSolverInterface * clpSolver
3313        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3314        if (clpSolver && !parentModel_) {
3315            clpSolver->computeLargestAway();
3316        }
3317    }
3318#endif
3319    /*
3320      At last, the actual branch-and-cut search loop, which will iterate until
3321      the live set is empty or we hit some limit (integrality gap, time, node
3322      count, etc.). The overall flow is to rebuild a subproblem, reoptimise using
3323      solveWithCuts(), choose a branching pattern with chooseBranch(), and finally
3324      add the node to the live set.
3325
3326      The first action is to winnow the live set to remove nodes which are worse
3327      than the current objective cutoff.
3328    */
3329    if (solver_->getRowCutDebuggerAlways()) {
3330        OsiRowCutDebugger * debuggerX = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
3331        const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger() ;
3332        if (!debugger) {
3333            // infeasible!!
3334            printf("before search\n");
3335            const double * lower = solver_->getColLower();
3336            const double * upper = solver_->getColUpper();
3337            const double * solution = debuggerX->optimalSolution();
3338            int numberColumns = solver_->getNumCols();
3339            for (int i = 0; i < numberColumns; i++) {
3340                if (solver_->isInteger(i)) {
3341                    if (solution[i] < lower[i] - 1.0e-6 || solution[i] > upper[i] + 1.0e-6)
3342                        printf("**** ");
3343                    printf("%d %g <= %g <= %g\n",
3344                           i, lower[i], solution[i], upper[i]);
3345                }
3346            }
3347            //abort();
3348        }
3349    }
3350    {
3351        // may be able to change cutoff now
3352        double cutoff = getCutoff();
3353        double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
3354        if (cutoff > bestObjective_ - increment) {
3355            cutoff = bestObjective_ - increment ;
3356            setCutoff(cutoff) ;
3357        }
3358    }
3359#ifdef CBC_THREAD
3360    bool goneParallel = false;
3361#endif
3362#define MAX_DEL_NODE 1
3363    CbcNode * delNode[MAX_DEL_NODE+1];
3364    int nDeleteNode = 0;
3365    // For Printing etc when parallel
3366    int lastEvery1000 = 0;
3367    int lastPrintEvery = 0;
3368    int numberConsecutiveInfeasible = 0;
3369#define PERTURB_IN_FATHOM
3370#ifdef PERTURB_IN_FATHOM
3371    // allow in fathom
3372    if ((moreSpecialOptions_& 262144) != 0)
3373      specialOptions_ |= 131072;
3374#endif
3375    while (true) {
3376        lockThread();
3377#ifdef COIN_HAS_CLP
3378        // See if we want dantzig row choice
3379        goToDantzig(100, savePivotMethod);
3380#endif
3381        if (tree_->empty()) {
3382#ifdef CBC_THREAD
3383            if (parallelMode() > 0 && master_) {
3384                int anyLeft = master_->waitForThreadsInTree(0);
3385                if (!anyLeft) {
3386                    master_->stopThreads(-1);
3387                    break;
3388                }
3389            } else {
3390                break;
3391            }
3392#else
3393            break;
3394#endif
3395        } else {
3396            unlockThread();
3397        }
3398        // If done 100 nodes see if worth trying reduction
3399        if (numberNodes_ == 50 || numberNodes_ == 100) {
3400#ifdef COIN_HAS_CLP
3401            OsiClpSolverInterface * clpSolver
3402            = dynamic_cast<OsiClpSolverInterface *> (solver_);
3403            if (clpSolver && ((specialOptions_&131072) == 0) && true) {
3404                ClpSimplex * simplex = clpSolver->getModelPtr();
3405                int perturbation = simplex->perturbation();
3406#ifdef CLP_INVESTIGATE
3407                printf("Testing its n,s %d %d solves n,s %d %d - pert %d\n",
3408                       numberIterations_, saveNumberIterations,
3409                       numberSolves_, saveNumberSolves, perturbation);
3410#endif
3411                if (perturbation == 50 && (numberIterations_ - saveNumberIterations) <
3412                        8*(numberSolves_ - saveNumberSolves)) {
3413                    // switch off perturbation
3414                    simplex->setPerturbation(100);
3415#ifdef CLP_INVESTIGATE
3416                    printf("Perturbation switched off\n");
3417#endif
3418                }
3419            }
3420#endif
3421            /*
3422              Decide if we want to do a restart.
3423            */
3424            if (saveSolver) {
3425                bool tryNewSearch = solverCharacteristics_->reducedCostsAccurate() &&
3426                                    (getCutoff() < 1.0e20 && getCutoff() < checkCutoffForRestart);
3427                int numberColumns = getNumCols();
3428                if (tryNewSearch) {
3429                    checkCutoffForRestart = getCutoff() ;
3430#ifdef CLP_INVESTIGATE
3431                    printf("after %d nodes, cutoff %g - looking\n",
3432                           numberNodes_, getCutoff());
3433#endif
3434                    saveSolver->resolve();
3435                    double direction = saveSolver->getObjSense() ;
3436                    double gap = checkCutoffForRestart - saveSolver->getObjValue() * direction ;
3437                    double tolerance;
3438                    saveSolver->getDblParam(OsiDualTolerance, tolerance) ;
3439                    if (gap <= 0.0)
3440                        gap = tolerance;
3441                    gap += 100.0 * tolerance;
3442                    double integerTolerance = getDblParam(CbcIntegerTolerance) ;
3443
3444                    const double *lower = saveSolver->getColLower() ;
3445                    const double *upper = saveSolver->getColUpper() ;
3446                    const double *solution = saveSolver->getColSolution() ;
3447                    const double *reducedCost = saveSolver->getReducedCost() ;
3448
3449                    int numberFixed = 0 ;
3450                    int numberFixed2 = 0;
3451#ifdef COIN_DEVELOP
3452                    printf("gap %g\n", gap);
3453#endif
3454                    for (int i = 0 ; i < numberIntegers_ ; i++) {
3455                        int iColumn = integerVariable_[i] ;
3456                        double djValue = direction * reducedCost[iColumn] ;
3457                        if (upper[iColumn] - lower[iColumn] > integerTolerance) {
3458                            if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) {
3459                                //printf("%d to lb on dj of %g - bounds %g %g\n",
3460                                //     iColumn,djValue,lower[iColumn],upper[iColumn]);
3461                                saveSolver->setColUpper(iColumn, lower[iColumn]) ;
3462                                numberFixed++ ;
3463                            } else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) {
3464                                //printf("%d to ub on dj of %g - bounds %g %g\n",
3465                                //     iColumn,djValue,lower[iColumn],upper[iColumn]);
3466                                saveSolver->setColLower(iColumn, upper[iColumn]) ;
3467                                numberFixed++ ;
3468                            }
3469                        } else {
3470                            //printf("%d has dj of %g - already fixed to %g\n",
3471                            //     iColumn,djValue,lower[iColumn]);
3472                            numberFixed2++;
3473                        }
3474                    }
3475#ifdef COIN_DEVELOP
3476                    if ((specialOptions_&1) != 0) {
3477                        const OsiRowCutDebugger *debugger = saveSolver->getRowCutDebugger() ;
3478                        if (debugger) {
3479                            printf("Contains optimal\n") ;
3480                            OsiSolverInterface * temp = saveSolver->clone();
3481                            const double * solution = debugger->optimalSolution();
3482                            const double *lower = temp->getColLower() ;
3483                            const double *upper = temp->getColUpper() ;
3484                            int n = temp->getNumCols();
3485                            for (int i = 0; i < n; i++) {
3486                                if (temp->isInteger(i)) {
3487                                    double value = floor(solution[i] + 0.5);
3488                                    assert (value >= lower[i] && value <= upper[i]);
3489                                    temp->setColLower(i, value);
3490                                    temp->setColUpper(i, value);
3491                                }
3492                            }
3493                            temp->writeMps("reduced_fix");
3494                            delete temp;
3495                            saveSolver->writeMps("reduced");
3496                        } else {
3497                            abort();
3498                        }
3499                    }
3500                    printf("Restart could fix %d integers (%d already fixed)\n",
3501                           numberFixed + numberFixed2, numberFixed2);
3502#endif
3503                    numberFixed += numberFixed2;
3504                    if (numberFixed*10 < numberColumns && numberFixed*4 < numberIntegers_)
3505                        tryNewSearch = false;
3506                }
3507                if (tryNewSearch) {
3508                    // back to solver without cuts?
3509                    OsiSolverInterface * solver2 = saveSolver->clone();
3510                    const double *lower = saveSolver->getColLower() ;
3511                    const double *upper = saveSolver->getColUpper() ;
3512                    for (int i = 0 ; i < numberIntegers_ ; i++) {
3513                        int iColumn = integerVariable_[i] ;
3514                        solver2->setColLower(iColumn, lower[iColumn]);
3515                        solver2->setColUpper(iColumn, upper[iColumn]);
3516                    }
3517                    // swap
3518                    delete saveSolver;
3519                    saveSolver = solver2;
3520                    double * newSolution = new double[numberColumns];
3521                    double objectiveValue = checkCutoffForRestart;
3522                    // Save the best solution so far.
3523                    CbcSerendipity heuristic(*this);
3524                    if (bestSolution_)
3525                        heuristic.setInputSolution(bestSolution_, bestObjective_);
3526                    // Magic number
3527                    heuristic.setFractionSmall(0.8);
3528                    // `pumpTune' to stand-alone solver for explanations.
3529                    heuristic.setFeasibilityPumpOptions(1008013);
3530                    // Use numberNodes to say how many are original rows
3531                    heuristic.setNumberNodes(continuousSolver_->getNumRows());
3532#ifdef COIN_DEVELOP
3533                    if (continuousSolver_->getNumRows() <
3534                            solver_->getNumRows())
3535                        printf("%d rows added ZZZZZ\n",
3536                               solver_->getNumRows() - continuousSolver_->getNumRows());
3537#endif
3538                    int returnCode = heuristic.smallBranchAndBound(saveSolver,
3539                                     -1, newSolution,
3540                                     objectiveValue,
3541                                     checkCutoffForRestart, "Reduce");
3542                    if (returnCode < 0) {
3543#ifdef COIN_DEVELOP
3544                        printf("Restart - not small enough to do search after fixing\n");
3545#endif
3546                        delete [] newSolution;
3547                    } else {
3548                        // 1 for sol'n, 2 for finished, 3 for both
3549                        if ((returnCode&1) != 0) {
3550                            // increment number of solutions so other heuristics can test
3551                            numberSolutions_++;
3552                            numberHeuristicSolutions_++;
3553                            lastHeuristic_ = NULL;
3554                            setBestSolution(CBC_ROUNDING, objectiveValue, newSolution) ;
3555                        }
3556                        delete [] newSolution;
3557#ifdef CBC_THREAD
3558                        if (master_) {
3559                            lockThread();
3560                            if (parallelMode() > 0) {
3561                                while (master_->waitForThreadsInTree(0)) {
3562                                    lockThread();
3563                                    double dummyBest;
3564                                    tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
3565                                    //unlockThread();
3566                                }
3567                            } else {
3568                                double dummyBest;
3569                                tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
3570                            }
3571                            master_->waitForThreadsInTree(2);
3572                            delete master_;
3573                            master_ = NULL;
3574                            masterThread_ = NULL;
3575                        }
3576#endif
3577                        if (tree_->size()) {
3578                            double dummyBest;
3579                            tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
3580                        }
3581                        break;
3582                    }
3583                }
3584                delete saveSolver;
3585                saveSolver = NULL;
3586            }
3587        }
3588        /*
3589          Check for abort on limits: node count, solution count, time, integrality gap.
3590        */
3591        if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
3592                numberSolutions_ < intParam_[CbcMaxNumSol] &&
3593                !maximumSecondsReached() &&
3594                !stoppedOnGap_ && !eventHappened_ && (maximumNumberIterations_ < 0 ||
3595                                                      numberIterations_ < maximumNumberIterations_))) {
3596            // out of loop
3597            break;
3598        }
3599#ifdef BONMIN
3600        assert(!solverCharacteristics_->solutionAddsCuts() || solverCharacteristics_->mipFeasible());
3601#endif
3602// Sets percentage of time when we try diving. Diving requires a bit of heap reorganisation, because
3603// we need to replace the comparison function to dive, and that requires reordering to retain the
3604// heap property.
3605#define DIVE_WHEN 1000
3606#define DIVE_STOP 2000
3607        int kNode = numberNodes_ % 4000;
3608        if (numberNodes_<100000 && kNode>DIVE_WHEN && kNode <= DIVE_STOP) {
3609            if (!parallelMode()) {
3610                if (kNode == DIVE_WHEN + 1 || numberConsecutiveInfeasible > 1) {
3611                    CbcCompareDefault * compare = dynamic_cast<CbcCompareDefault *>
3612                                                  (nodeCompare_);
3613                    // Don't interfere if user has replaced the compare function.
3614                    if (compare) {
3615                        //printf("Redoing tree\n");
3616                        compare->startDive(this);
3617                        numberConsecutiveInfeasible = 0;
3618                    }
3619                }
3620            }
3621        }
3622        // replace current cutoff?
3623        if (cutoff > getCutoff()) {
3624            double newCutoff = getCutoff();
3625            if (analyzeResults_) {
3626                // see if we could fix any (more)
3627                int n = 0;
3628                double * newLower = analyzeResults_;
3629                double * objLower = newLower + numberIntegers_;
3630                double * newUpper = objLower + numberIntegers_;
3631                double * objUpper = newUpper + numberIntegers_;
3632                for (int i = 0; i < numberIntegers_; i++) {
3633                    if (objLower[i] > newCutoff) {
3634                        n++;
3635                        if (objUpper[i] > newCutoff) {
3636                            newCutoff = -COIN_DBL_MAX;
3637                            break;
3638                        }
3639                    } else if (objUpper[i] > newCutoff) {
3640                        n++;
3641                    }
3642                }
3643                if (newCutoff == -COIN_DBL_MAX) {
3644                  COIN_DETAIL_PRINT(printf("Root analysis says finished\n"));
3645                } else if (n > numberFixedNow_) {
3646                  COIN_DETAIL_PRINT(printf("%d more fixed by analysis - now %d\n", n - numberFixedNow_, n));
3647                    numberFixedNow_ = n;
3648                }
3649            }
3650            if (eventHandler) {
3651                if (!eventHandler->event(CbcEventHandler::solution)) {
3652                    eventHappened_ = true; // exit
3653                }
3654                newCutoff = getCutoff();
3655            }
3656            lockThread();
3657            /*
3658              Clean the tree to reflect the new solution, then see if the
3659              node comparison predicate wants to make any changes. If so,
3660              call setComparison for the side effect of rebuilding the heap.
3661            */
3662            tree_->cleanTree(this,newCutoff,bestPossibleObjective_) ;
3663            if (nodeCompare_->newSolution(this) ||
3664                nodeCompare_->newSolution(this,continuousObjective_,
3665                                          continuousInfeasibilities_)) {
3666              tree_->setComparison(*nodeCompare_) ;
3667            }
3668            if (tree_->empty()) {
3669                continue;
3670            }
3671            unlockThread();
3672        }
3673        cutoff = getCutoff() ;
3674        /*
3675            Periodic activities: Opportunities to
3676            + tweak the nodeCompare criteria,
3677            + check if we've closed the integrality gap enough to quit,
3678            + print a summary line to let the user know we're working
3679        */
3680        if (numberNodes_ >= lastEvery1000) {
3681            lockThread();
3682#ifdef COIN_HAS_CLP
3683            // See if we want dantzig row choice
3684            goToDantzig(1000, savePivotMethod);
3685#endif
3686            lastEvery1000 = numberNodes_ + 1000;
3687            bool redoTree = nodeCompare_->every1000Nodes(this, numberNodes_) ;
3688#ifdef CHECK_CUT_SIZE
3689            verifyCutSize (tree_, *this);
3690#endif
3691            // redo tree if requested
3692            if (redoTree)
3693                tree_->setComparison(*nodeCompare_) ;
3694            unlockThread();
3695        }
3696        // Had hotstart before, now switched off
3697        if (saveCompare && !hotstartSolution_) {
3698            // hotstart switched off
3699            delete nodeCompare_; // off depth first
3700            nodeCompare_ = saveCompare;
3701            saveCompare = NULL;
3702            // redo tree
3703            lockThread();
3704            tree_->setComparison(*nodeCompare_) ;
3705            unlockThread();
3706        }
3707        if (numberNodes_ >= lastPrintEvery) {
3708            lastPrintEvery = numberNodes_ + printFrequency_;
3709            lockThread();
3710            int nNodes = tree_->size() ;
3711
3712            //MODIF PIERRE
3713            bestPossibleObjective_ = tree_->getBestPossibleObjective();
3714            unlockThread();
3715#ifdef CLP_INVESTIGATE
3716            if (getCutoff() < 1.0e20) {
3717                if (fabs(getCutoff() - (bestObjective_ - getCutoffIncrement())) > 1.0e-6 &&
3718                        !parentModel_)
3719                    printf("model cutoff in status %g, best %g, increment %g\n",
3720                           getCutoff(), bestObjective_, getCutoffIncrement());
3721                assert (getCutoff() < bestObjective_ - getCutoffIncrement() +
3722                        1.0e-6 + 1.0e-10*fabs(bestObjective_));
3723            }
3724#endif
3725            if (!intParam_[CbcPrinting]) {
3726                messageHandler()->message(CBC_STATUS, messages())
3727                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
3728                << getCurrentSeconds()
3729                << CoinMessageEol ;
3730            } else if (intParam_[CbcPrinting] == 1) {
3731                messageHandler()->message(CBC_STATUS2, messages())
3732                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
3733                << lastDepth << lastUnsatisfied << numberIterations_
3734                << getCurrentSeconds()
3735                << CoinMessageEol ;
3736            } else if (!numberExtraIterations_) {
3737                messageHandler()->message(CBC_STATUS2, messages())
3738                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
3739                << lastDepth << lastUnsatisfied << numberIterations_
3740                << getCurrentSeconds()
3741                << CoinMessageEol ;
3742            } else {
3743                messageHandler()->message(CBC_STATUS3, messages())
3744                << numberNodes_ << numberExtraNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
3745                << lastDepth << lastUnsatisfied << numberIterations_ << numberExtraIterations_
3746                << getCurrentSeconds()
3747                << CoinMessageEol ;
3748            }
3749            if (eventHandler && !eventHandler->event(CbcEventHandler::treeStatus)) {
3750                eventHappened_ = true; // exit
3751            }
3752        }
3753        // See if can stop on gap
3754        double testGap = CoinMax(dblParam_[CbcAllowableGap],
3755                                 CoinMax(fabs(bestObjective_), fabs(bestPossibleObjective_))
3756                                 * dblParam_[CbcAllowableFractionGap]);
3757        if (bestObjective_ - bestPossibleObjective_ < testGap && getCutoffIncrement() >= 0.0) {
3758            stoppedOnGap_ = true ;
3759        }
3760
3761#ifdef CHECK_NODE_FULL
3762        verifyTreeNodes(tree_, *this) ;
3763#   endif
3764#   ifdef CHECK_CUT_COUNTS
3765        verifyCutCounts(tree_, *this) ;
3766#   endif
3767        /*
3768          Now we come to the meat of the loop. To create the active subproblem, we'll
3769          pop the most promising node in the live set, rebuild the subproblem it
3770          represents, and then execute the current arm of the branch to create the
3771          active subproblem.
3772        */
3773        CbcNode * node = NULL;
3774#ifdef CBC_THREAD
3775        if (!parallelMode() || parallelMode() == -1) {
3776#endif
3777            node = tree_->bestNode(cutoff) ;
3778            // Possible one on tree worse than cutoff
3779            // Weird comparison function can leave ineligible nodes on tree
3780            if (!node || node->objectiveValue() > cutoff)
3781                continue;
3782            // Do main work of solving node here
3783            doOneNode(this, node, createdNode);
3784#ifdef JJF_ZERO
3785            if (node) {
3786                if (createdNode) {
3787                    printf("Node %d depth %d, created %d depth %d\n",
3788                           node->nodeNumber(), node->depth(),
3789                           createdNode->nodeNumber(), createdNode->depth());
3790                } else {
3791                    printf("Node %d depth %d,  no created node\n",
3792                           node->nodeNumber(), node->depth());
3793                }
3794            } else if (createdNode) {
3795                printf("Node exhausted, created %d depth %d\n",
3796                       createdNode->nodeNumber(), createdNode->depth());
3797            } else {
3798                printf("Node exhausted,  no created node\n");
3799                numberConsecutiveInfeasible = 2;
3800            }
3801#endif
3802            //if (createdNode)
3803            //numberConsecutiveInfeasible=0;
3804            //else
3805            //numberConsecutiveInfeasible++;
3806#ifdef CBC_THREAD
3807        } else if (parallelMode() > 0) {
3808            //lockThread();
3809            //node = tree_->bestNode(cutoff) ;
3810            // Possible one on tree worse than cutoff
3811            if (true || !node || node->objectiveValue() > cutoff) {
3812                assert (master_);
3813                if (master_) {
3814                    int anyLeft = master_->waitForThreadsInTree(1);
3815                    // may need to go round again
3816                    if (anyLeft) {
3817                        continue;
3818                    } else {
3819                        master_->stopThreads(-1);
3820                    }
3821                }
3822            }
3823            //unlockThread();
3824        } else {
3825            // Deterministic parallel
3826            if (tree_->size() < CoinMax(numberThreads_, 8) && !goneParallel) {
3827                node = tree_->bestNode(cutoff) ;
3828                // Possible one on tree worse than cutoff
3829                if (!node || node->objectiveValue() > cutoff)
3830                    continue;
3831                // Do main work of solving node here
3832                doOneNode(this, node, createdNode);
3833                assert (createdNode);
3834                if (!createdNode->active()) {
3835                    delete createdNode;
3836                    createdNode = NULL;
3837                } else {
3838                    // Say one more pointing to this
3839                    node->nodeInfo()->increment() ;
3840                    tree_->push(createdNode) ;
3841                }
3842                if (node->active()) {
3843                    assert (node->nodeInfo());
3844                    if (node->nodeInfo()->numberBranchesLeft()) {
3845                        tree_->push(node) ;
3846                    } else {
3847                        node->setActive(false);
3848                    }
3849                } else {
3850                    if (node->nodeInfo()) {
3851                        if (!node->nodeInfo()->numberBranchesLeft())
3852                            node->nodeInfo()->allBranchesGone(); // can clean up
3853                        // So will delete underlying stuff
3854                        node->setActive(true);
3855                    }
3856                    delNode[nDeleteNode++] = node;
3857                    node = NULL;
3858                }
3859                if (nDeleteNode >= MAX_DEL_NODE) {
3860                    for (int i = 0; i < nDeleteNode; i++) {
3861                        //printf("trying to del %d %x\n",i,delNode[i]);
3862                        delete delNode[i];
3863                        //printf("done to del %d %x\n",i,delNode[i]);
3864                    }
3865                    nDeleteNode = 0;
3866                }
3867            } else {
3868                // Split and solve
3869                master_->deterministicParallel();
3870                goneParallel = true;
3871            }
3872        }
3873#endif
3874    }
3875    if (nDeleteNode) {
3876        for (int i = 0; i < nDeleteNode; i++) {
3877            delete delNode[i];
3878        }
3879        nDeleteNode = 0;
3880    }
3881#ifdef CBC_THREAD
3882    if (master_) {
3883        master_->stopThreads(-1);
3884        master_->waitForThreadsInTree(2);
3885        delete master_;
3886        master_ = NULL;
3887        masterThread_ = NULL;
3888        // adjust time to allow for children on some systems
3889        //dblParam_[CbcStartSeconds] -= CoinCpuTimeJustChildren();
3890    }
3891#endif
3892    /*
3893      End of the non-abort actions. The next block of code is executed if we've
3894      aborted because we hit one of the limits. Clean up by deleting the live set
3895      and break out of the node processing loop. Note that on an abort, node may
3896      have been pushed back onto the tree for further processing, in which case
3897      it'll be deleted in cleanTree. We need to check.
3898    */
3899    if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
3900            numberSolutions_ < intParam_[CbcMaxNumSol] &&
3901            !maximumSecondsReached() &&
3902            !stoppedOnGap_ && !eventHappened_ && (maximumNumberIterations_ < 0 ||
3903                                                  numberIterations_ < maximumNumberIterations_))) {
3904        if (tree_->size()) {
3905            double dummyBest;
3906            tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
3907        }
3908        delete nextRowCut_;
3909        if (stoppedOnGap_) {
3910            messageHandler()->message(CBC_GAP, messages())
3911            << bestObjective_ - bestPossibleObjective_
3912            << dblParam_[CbcAllowableGap]
3913            << dblParam_[CbcAllowableFractionGap]*100.0
3914            << CoinMessageEol ;
3915            secondaryStatus_ = 2;
3916            status_ = 0 ;
3917        } else if (isNodeLimitReached()) {
3918            handler_->message(CBC_MAXNODES, messages_) << CoinMessageEol ;
3919            secondaryStatus_ = 3;
3920            status_ = 1 ;
3921        } else if (maximumSecondsReached()) {
3922            handler_->message(CBC_MAXTIME, messages_) << CoinMessageEol ;
3923            secondaryStatus_ = 4;
3924            status_ = 1 ;
3925        } else if (eventHappened_) {
3926            handler_->message(CBC_EVENT, messages_) << CoinMessageEol ;
3927            secondaryStatus_ = 5;
3928            status_ = 5 ;
3929        } else {
3930            handler_->message(CBC_MAXSOLS, messages_) << CoinMessageEol ;
3931            secondaryStatus_ = 6;
3932            status_ = 1 ;
3933        }
3934    }
3935    /*
3936      That's it, we've exhausted the search tree, or broken out of the loop because
3937      we hit some limit on evaluation.
3938
3939      We may have got an intelligent tree so give it one more chance
3940    */
3941    // Tell solver we are not in Branch and Cut
3942    solver_->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo, NULL) ;
3943    tree_->endSearch();
3944    //  If we did any sub trees - did we give up on any?
3945    if ( numberStoppedSubTrees_)
3946        status_ = 1;
3947    numberNodes_ += numberExtraNodes_;
3948    numberIterations_ += numberExtraIterations_;
3949    if (eventHandler) {
3950        eventHandler->event(CbcEventHandler::endSearch);
3951    }
3952    if (!status_) {
3953        // Set best possible unless stopped on gap
3954        if (secondaryStatus_ != 2)
3955            bestPossibleObjective_ = bestObjective_;
3956        handler_->message(CBC_END_GOOD, messages_)
3957        << bestObjective_ << numberIterations_ << numberNodes_ << getCurrentSeconds()
3958        << CoinMessageEol ;
3959    } else {
3960        handler_->message(CBC_END, messages_)
3961        << bestObjective_ << bestPossibleObjective_
3962        << numberIterations_ << numberNodes_ << getCurrentSeconds()
3963        << CoinMessageEol ;
3964    }
3965    if (numberStrongIterations_)
3966        handler_->message(CBC_STRONG_STATS, messages_)
3967        << strongInfo_[0] << numberStrongIterations_ << strongInfo_[2]
3968        << strongInfo_[1] << CoinMessageEol ;
3969    if (!numberExtraNodes_)
3970        handler_->message(CBC_OTHER_STATS, messages_)
3971        << maximumDepthActual_
3972        << numberDJFixed_ << CoinMessageEol ;
3973    else
3974        handler_->message(CBC_OTHER_STATS2, messages_)
3975        << maximumDepthActual_
3976        << numberDJFixed_ << numberExtraNodes_ << numberExtraIterations_
3977        << CoinMessageEol ;
3978    if (doStatistics == 100) {
3979        for (int i = 0; i < numberObjects_; i++) {
3980            CbcSimpleIntegerDynamicPseudoCost * obj =
3981                dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object_[i]) ;
3982            if (obj)
3983                obj->print();
3984        }
3985    }
3986    if (statistics_) {
3987        // report in some way
3988        int * lookup = new int[numberObjects_];
3989        int i;
3990        for (i = 0; i < numberObjects_; i++)
3991            lookup[i] = -1;
3992        bool goodIds = false; //true;
3993        for (i = 0; i < numberObjects_; i++) {
3994            int iColumn = object_[i]->columnNumber();
3995            if (iColumn >= 0 && iColumn < numberColumns) {
3996                if (lookup[i] == -1) {
3997                    lookup[i] = iColumn;
3998                } else {
3999                    goodIds = false;
4000                    break;
4001                }
4002            } else {
4003                goodIds = false;
4004                break;
4005            }
4006        }
4007        if (!goodIds) {
4008            delete [] lookup;
4009            lookup = NULL;
4010        }
4011        if (doStatistics >= 3) {
4012            printf("  node parent depth column   value                    obj      inf\n");
4013            for ( i = 0; i < numberNodes2_; i++) {
4014                statistics_[i]->print(lookup);
4015            }
4016        }
4017        if (doStatistics > 1) {
4018            // Find last solution
4019            int k;
4020            for (k = numberNodes2_ - 1; k >= 0; k--) {
4021                if (statistics_[k]->endingObjective() != COIN_DBL_MAX &&
4022                        !statistics_[k]->endingInfeasibility())
4023                    break;
4024            }
4025            if (k >= 0) {
4026                int depth = statistics_[k]->depth();
4027                int * which = new int[depth+1];
4028                for (i = depth; i >= 0; i--) {
4029                    which[i] = k;
4030                    k = statistics_[k]->parentNode();
4031                }
4032                printf("  node parent depth column   value                    obj      inf\n");
4033                for (i = 0; i <= depth; i++) {
4034                    statistics_[which[i]]->print(lookup);
4035                }
4036                delete [] which;
4037            }
4038        }
4039        // now summary
4040        int maxDepth = 0;
4041        double averageSolutionDepth = 0.0;
4042        int numberSolutions = 0;
4043        double averageCutoffDepth = 0.0;
4044        double averageSolvedDepth = 0.0;
4045        int numberCutoff = 0;
4046        int numberDown = 0;
4047        int numberFirstDown = 0;
4048        double averageInfDown = 0.0;
4049        double averageObjDown = 0.0;
4050        int numberCutoffDown = 0;
4051        int numberUp = 0;
4052        int numberFirstUp = 0;
4053        double averageInfUp = 0.0;
4054        double averageObjUp = 0.0;
4055        int numberCutoffUp = 0;
4056        double averageNumberIterations1 = 0.0;
4057        double averageValue = 0.0;
4058        for ( i = 0; i < numberNodes2_; i++) {
4059            int depth =  statistics_[i]->depth();
4060            int way =  statistics_[i]->way();
4061            double value = statistics_[i]->value();
4062            double startingObjective =  statistics_[i]->startingObjective();
4063            int startingInfeasibility = statistics_[i]->startingInfeasibility();
4064            double endingObjective = statistics_[i]->endingObjective();
4065            int endingInfeasibility = statistics_[i]->endingInfeasibility();
4066            maxDepth = CoinMax(depth, maxDepth);
4067            // Only for completed
4068            averageNumberIterations1 += statistics_[i]->numberIterations();
4069            averageValue += value;
4070            if (endingObjective != COIN_DBL_MAX && !endingInfeasibility) {
4071                numberSolutions++;
4072                averageSolutionDepth += depth;
4073            }
4074            if (endingObjective == COIN_DBL_MAX) {
4075                numberCutoff++;
4076                averageCutoffDepth += depth;
4077                if (way < 0) {
4078                    numberDown++;
4079                    numberCutoffDown++;
4080                    if (way == -1)
4081                        numberFirstDown++;
4082                } else {
4083                    numberUp++;
4084                    numberCutoffUp++;
4085                    if (way == 1)
4086                        numberFirstUp++;
4087                }
4088            } else {
4089                averageSolvedDepth += depth;
4090                if (way < 0) {
4091                    numberDown++;
4092                    averageInfDown += startingInfeasibility - endingInfeasibility;
4093                    averageObjDown += endingObjective - startingObjective;
4094                    if (way == -1)
4095                        numberFirstDown++;
4096                } else {
4097                    numberUp++;
4098                    averageInfUp += startingInfeasibility - endingInfeasibility;
4099                    averageObjUp += endingObjective - startingObjective;
4100                    if (way == 1)
4101                        numberFirstUp++;
4102                }
4103            }
4104        }
4105        // Now print
4106        if (numberSolutions)
4107            averageSolutionDepth /= static_cast<double> (numberSolutions);
4108        int numberSolved = numberNodes2_ - numberCutoff;
4109        double averageNumberIterations2 = numberIterations_ - averageNumberIterations1
4110                                          - numberIterationsAtContinuous;
4111        if (numberCutoff) {
4112            averageCutoffDepth /= static_cast<double> (numberCutoff);
4113            averageNumberIterations2 /= static_cast<double> (numberCutoff);
4114        }
4115        if (numberNodes2_)
4116            averageValue /= static_cast<double> (numberNodes2_);
4117        if (numberSolved) {
4118            averageNumberIterations1 /= static_cast<double> (numberSolved);
4119            averageSolvedDepth /= static_cast<double> (numberSolved);
4120        }
4121        printf("%d solution(s) were found (by branching) at an average depth of %g\n",
4122               numberSolutions, averageSolutionDepth);
4123        printf("average value of variable being branched on was %g\n",
4124               averageValue);
4125        printf("%d nodes were cutoff at an average depth of %g with iteration count of %g\n",
4126               numberCutoff, averageCutoffDepth, averageNumberIterations2);
4127        printf("%d nodes were solved at an average depth of %g with iteration count of %g\n",
4128               numberSolved, averageSolvedDepth, averageNumberIterations1);
4129        if (numberDown) {
4130            averageInfDown /= static_cast<double> (numberDown);
4131            averageObjDown /= static_cast<double> (numberDown);
4132        }
4133        printf("Down %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
4134               numberDown, numberFirstDown, numberDown - numberFirstDown, numberCutoffDown,
4135               averageInfDown, averageObjDown);
4136        if (numberUp) {
4137            averageInfUp /= static_cast<double> (numberUp);
4138            averageObjUp /= static_cast<double> (numberUp);
4139        }
4140        printf("Up %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
4141               numberUp, numberFirstUp, numberUp - numberFirstUp, numberCutoffUp,
4142               averageInfUp, averageObjUp);
4143        for ( i = 0; i < numberNodes2_; i++)
4144            delete statistics_[i];
4145        delete [] statistics_;
4146        statistics_ = NULL;
4147        maximumStatistics_ = 0;
4148        delete [] lookup;
4149    }
4150    /*
4151      If we think we have a solution, restore and confirm it with a call to
4152      setBestSolution().  We need to reset the cutoff value so as not to fathom
4153      the solution on bounds.  Note that calling setBestSolution( ..., true)
4154      leaves the continuousSolver_ bounds vectors fixed at the solution value.
4155
4156      Running resolve() here is a failsafe --- setBestSolution has already
4157      reoptimised using the continuousSolver_. If for some reason we fail to
4158      prove optimality, run the problem again after instructing the solver to
4159      tell us more.
4160
4161      If all looks good, replace solver_ with continuousSolver_, so that the
4162      outside world will be able to obtain information about the solution using
4163      public methods.
4164    */
4165    if (bestSolution_ && (solverCharacteristics_->solverType() < 2 || solverCharacteristics_->solverType() == 4)) {
4166        setCutoff(1.0e50) ; // As best solution should be worse than cutoff
4167        phase_ = 5;
4168        double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
4169        if ((specialOptions_&4) == 0)
4170            bestObjective_ += 100.0 * increment + 1.0e-3; // only set if we are going to solve
4171        setBestSolution(CBC_END_SOLUTION, bestObjective_, bestSolution_, 1) ;
4172        continuousSolver_->resolve() ;
4173        if (!continuousSolver_->isProvenOptimal()) {
4174            continuousSolver_->messageHandler()->setLogLevel(2) ;
4175            continuousSolver_->initialSolve() ;
4176        }
4177        delete solver_ ;
4178        // above deletes solverCharacteristics_
4179        solverCharacteristics_ = NULL;
4180        solver_ = continuousSolver_ ;
4181        setPointers(solver_);
4182        continuousSolver_ = NULL ;
4183    }
4184    /*
4185      Clean up dangling objects. continuousSolver_ may already be toast.
4186    */
4187    delete lastws ;
4188    if (saveObjects) {
4189        for (int i = 0; i < numberObjects_; i++)
4190            delete saveObjects[i];
4191        delete [] saveObjects;
4192    }
4193    numberStrong_ = saveNumberStrong;
4194    numberBeforeTrust_ = saveNumberBeforeTrust;
4195    delete [] whichGenerator_ ;
4196    whichGenerator_ = NULL;
4197    delete [] lowerBefore ;
4198    delete [] upperBefore ;
4199    delete [] walkback_ ;
4200    walkback_ = NULL ;
4201    delete [] lastNodeInfo_ ;
4202    lastNodeInfo_ = NULL;
4203    delete [] lastNumberCuts_ ;
4204    lastNumberCuts_ = NULL;
4205    delete [] lastCut_;
4206    lastCut_ = NULL;
4207    delete [] addedCuts_ ;
4208    addedCuts_ = NULL ;
4209    //delete persistentInfo;
4210    // Get rid of characteristics
4211    solverCharacteristics_ = NULL;
4212    if (continuousSolver_) {
4213        delete continuousSolver_ ;
4214        continuousSolver_ = NULL ;
4215    }
4216    /*
4217      Destroy global cuts by replacing with an empty OsiCuts object.
4218    */
4219    globalCuts_ = OsiCuts() ;
4220    if (!bestSolution_) {
4221        // make sure lp solver is infeasible
4222        int numberColumns = solver_->getNumCols();
4223        const double * columnLower = solver_->getColLower();
4224        int iColumn;
4225        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4226            if (solver_->isInteger(iColumn))
4227                solver_->setColUpper(iColumn, columnLower[iColumn]);
4228        }
4229        solver_->initialSolve();
4230    }
4231#ifdef COIN_HAS_CLP
4232    {
4233        OsiClpSolverInterface * clpSolver
4234        = dynamic_cast<OsiClpSolverInterface *> (solver_);
4235        if (clpSolver) {
4236            // Possible restore of pivot method
4237            if (savePivotMethod) {
4238                // model may have changed
4239                savePivotMethod->setModel(NULL);
4240                clpSolver->getModelPtr()->setDualRowPivotAlgorithm(*savePivotMethod);
4241                delete savePivotMethod;
4242            }
4243            clpSolver->setLargestAway(-1.0);
4244        }
4245    }
4246#endif
4247    if (fastNodeDepth_ >= 1000000 && !parentModel_) {
4248        // delete object off end
4249        delete object_[numberObjects_];
4250        fastNodeDepth_ -= 1000000;
4251    }
4252    delete saveSolver;
4253    // Undo preprocessing performed during BaB.
4254    if (strategy_ && strategy_->preProcessState() > 0) {
4255        // undo preprocessing
4256        CglPreProcess * process = strategy_->process();
4257        assert (process);
4258        int n = originalSolver->getNumCols();
4259        if (bestSolution_) {
4260            delete [] bestSolution_;
4261            bestSolution_ = new double [n];
4262            process->postProcess(*solver_);
4263        }
4264        strategy_->deletePreProcess();
4265        // Solution now back in originalSolver
4266        delete solver_;
4267        solver_ = originalSolver;
4268        if (bestSolution_) {
4269            bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
4270            memcpy(bestSolution_, solver_->getColSolution(), n*sizeof(double));
4271        }
4272        // put back original objects if there were any
4273        if (originalObject) {
4274            int iColumn;
4275            assert (ownObjects_);
4276            for (iColumn = 0; iColumn < numberObjects_; iColumn++)
4277                delete object_[iColumn];
4278            delete [] object_;
4279            numberObjects_ = numberOriginalObjects;
4280            object_ = originalObject;
4281            delete [] integerVariable_;
4282            numberIntegers_ = 0;
4283            for (iColumn = 0; iColumn < n; iColumn++) {
4284                if (solver_->isInteger(iColumn))
4285                    numberIntegers_++;
4286            }
4287            integerVariable_ = new int[numberIntegers_];
4288            numberIntegers_ = 0;
4289            for (iColumn = 0; iColumn < n; iColumn++) {
4290                if (solver_->isInteger(iColumn))
4291                    integerVariable_[numberIntegers_++] = iColumn;
4292            }
4293        }
4294    }
4295#ifdef COIN_HAS_CLP
4296    {
4297        OsiClpSolverInterface * clpSolver
4298        = dynamic_cast<OsiClpSolverInterface *> (solver_);
4299        if (clpSolver)
4300            clpSolver->setFakeObjective(reinterpret_cast<double *> (NULL));
4301    }
4302#endif
4303    moreSpecialOptions_ = saveMoreSpecialOptions;
4304    return ;
4305}
4306
4307
4308// Solve the initial LP relaxation
4309void
4310CbcModel::initialSolve()
4311{
4312    assert (solver_);
4313    // Double check optimization directions line up
4314    dblParam_[CbcOptimizationDirection] = solver_->getObjSense();
4315    // Check if bounds are all integral (as may get messed up later)
4316    checkModel();
4317    if (!solverCharacteristics_) {
4318        OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
4319        if (solverCharacteristics) {
4320            solverCharacteristics_ = solverCharacteristics;
4321        } else {
4322            // replace in solver
4323            OsiBabSolver defaultC;
4324            solver_->setAuxiliaryInfo(&defaultC);
4325            solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
4326        }
4327    }
4328    solverCharacteristics_->setSolver(solver_);
4329    solver_->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
4330    solver_->initialSolve();
4331    solver_->setHintParam(OsiDoInBranchAndCut, false, OsiHintDo, NULL) ;
4332    if (!solver_->isProvenOptimal())
4333        solver_->resolve();
4334    // But set up so Jon Lee will be happy
4335    status_ = -1;
4336    secondaryStatus_ = -1;
4337    originalContinuousObjective_ = solver_->getObjValue() * solver_->getObjSense();
4338    bestPossibleObjective_ = originalContinuousObjective_;
4339    delete [] continuousSolution_;
4340    continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
4341                                          solver_->getNumCols());
4342    setPointers(solver_);
4343    solverCharacteristics_ = NULL;
4344}
4345
4346/*! \brief Get an empty basis object
4347
4348  Return an empty CoinWarmStartBasis object with the requested capacity,
4349  appropriate for the current solver. The object is cloned from the object
4350  cached as emptyWarmStart_. If there is no cached object, the routine
4351  queries the solver for a warm start object, empties it, and caches the
4352  result.
4353*/
4354
4355CoinWarmStartBasis *CbcModel::getEmptyBasis (int ns, int na) const
4356
4357{
4358    CoinWarmStartBasis *emptyBasis ;
4359    /*
4360      Acquire an empty basis object, if we don't yet have one.
4361    */
4362    if (emptyWarmStart_ == 0) {
4363        if (solver_ == 0) {
4364            throw CoinError("Cannot construct basis without solver!",
4365                            "getEmptyBasis", "CbcModel") ;
4366        }
4367        emptyBasis =
4368            dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart()) ;
4369        if (emptyBasis == 0) {
4370            throw CoinError(
4371                "Solver does not appear to use a basis-oriented warm start.",
4372                "getEmptyBasis", "CbcModel") ;
4373        }
4374        emptyBasis->setSize(0, 0) ;
4375        emptyWarmStart_ = dynamic_cast<CoinWarmStart *>(emptyBasis) ;
4376    }
4377    /*
4378      Clone the empty basis object, resize it as requested, and return.
4379    */
4380    emptyBasis = dynamic_cast<CoinWarmStartBasis *>(emptyWarmStart_->clone()) ;
4381    assert(emptyBasis) ;
4382    if (ns != 0 || na != 0) emptyBasis->setSize(ns, na) ;
4383
4384    return (emptyBasis) ;
4385}
4386
4387
4388/** Default Constructor
4389
4390  Creates an empty model without an associated solver.
4391*/
4392CbcModel::CbcModel()
4393
4394        :
4395        solver_(NULL),
4396        ownership_(0x80000000),
4397        continuousSolver_(NULL),
4398        referenceSolver_(NULL),
4399        defaultHandler_(true),
4400        emptyWarmStart_(NULL),
4401        bestObjective_(COIN_DBL_MAX),
4402        bestPossibleObjective_(COIN_DBL_MAX),
4403        sumChangeObjective1_(0.0),
4404        sumChangeObjective2_(0.0),
4405        bestSolution_(NULL),
4406        savedSolutions_(NULL),
4407        currentSolution_(NULL),
4408        testSolution_(NULL),
4409        minimumDrop_(1.0e-4),
4410        numberSolutions_(0),
4411        numberSavedSolutions_(0),
4412        maximumSavedSolutions_(0),
4413        stateOfSearch_(0),
4414        whenCuts_(-1),
4415        hotstartSolution_(NULL),
4416        hotstartPriorities_(NULL),
4417        numberHeuristicSolutions_(0),
4418        numberNodes_(0),
4419        numberNodes2_(0),
4420        numberIterations_(0),
4421        numberSolves_(0),
4422        status_(-1),
4423        secondaryStatus_(-1),
4424        numberIntegers_(0),
4425        numberRowsAtContinuous_(0),
4426        maximumNumberCuts_(0),
4427        phase_(0),
4428        currentNumberCuts_(0),
4429        maximumDepth_(0),
4430        walkback_(NULL),
4431        lastNodeInfo_(NULL),
4432        lastCut_(NULL),
4433        lastDepth_(0),
4434        lastNumberCuts2_(0),
4435        maximumCuts_(0),
4436        lastNumberCuts_(NULL),
4437        addedCuts_(NULL),
4438        nextRowCut_(NULL),
4439        currentNode_(NULL),
4440        integerVariable_(NULL),
4441        integerInfo_(NULL),
4442        continuousSolution_(NULL),
4443        usedInSolution_(NULL),
4444        specialOptions_(0),
4445        moreSpecialOptions_(0),
4446        subTreeModel_(NULL),
4447        numberStoppedSubTrees_(0),
4448        presolve_(0),
4449        numberStrong_(5),
4450        numberBeforeTrust_(10),
4451        numberPenalties_(20),
4452        stopNumberIterations_(-1),
4453        penaltyScaleFactor_(3.0),
4454        numberAnalyzeIterations_(0),
4455        analyzeResults_(NULL),
4456        numberInfeasibleNodes_(0),
4457        problemType_(0),
4458        printFrequency_(0),
4459        numberCutGenerators_(0),
4460        generator_(NULL),
4461        virginGenerator_(NULL),
4462        numberHeuristics_(0),
4463        heuristic_(NULL),
4464        lastHeuristic_(NULL),
4465# ifdef COIN_HAS_CLP
4466        fastNodeDepth_(-1),
4467#endif
4468        eventHandler_(NULL),
4469        numberObjects_(0),
4470        object_(NULL),
4471        ownObjects_(true),
4472        originalColumns_(NULL),
4473        howOftenGlobalScan_(1),
4474        numberGlobalViolations_(0),
4475        numberExtraIterations_(0),
4476        numberExtraNodes_(0),
4477        continuousObjective_(COIN_DBL_MAX),
4478        originalContinuousObjective_(COIN_DBL_MAX),
4479        continuousInfeasibilities_(COIN_INT_MAX),
4480        maximumCutPassesAtRoot_(20),
4481        maximumCutPasses_(10),
4482        preferredWay_(0),
4483        currentPassNumber_(0),
4484        maximumWhich_(1000),
4485        maximumRows_(0),
4486        currentDepth_(0),
4487        whichGenerator_(NULL),
4488        maximumStatistics_(0),
4489        statistics_(NULL),
4490        maximumDepthActual_(0),
4491        numberDJFixed_(0.0),
4492        probingInfo_(NULL),
4493        numberFixedAtRoot_(0),
4494        numberFixedNow_(0),
4495        stoppedOnGap_(false),
4496        eventHappened_(false),
4497        numberLongStrong_(0),
4498        numberOldActiveCuts_(0),
4499        numberNewCuts_(0),
4500        searchStrategy_(-1),
4501        numberStrongIterations_(0),
4502        resolveAfterTakeOffCuts_(true),
4503        maximumNumberIterations_(-1),
4504        continuousPriority_(COIN_INT_MAX),
4505        numberUpdateItems_(0),
4506        maximumNumberUpdateItems_(0),
4507        updateItems_(NULL),
4508        storedRowCuts_(NULL),
4509        numberThreads_(0),
4510        threadMode_(0),
4511        master_(NULL),
4512        masterThread_(NULL)
4513{
4514    memset(intParam_, 0, sizeof(intParam_));
4515    intParam_[CbcMaxNumNode] = 2147483647;
4516    intParam_[CbcMaxNumSol] = 9999999;
4517
4518    memset(dblParam_, 0, sizeof(dblParam_));
4519    dblParam_[CbcIntegerTolerance] = 1e-6;
4520    dblParam_[CbcCutoffIncrement] = 1e-5;
4521    dblParam_[CbcAllowableGap] = 1.0e-10;
4522    dblParam_[CbcMaximumSeconds] = 1.0e100;
4523    dblParam_[CbcCurrentCutoff] = 1.0e100;
4524    dblParam_[CbcOptimizationDirection] = 1.0;
4525    dblParam_[CbcCurrentObjectiveValue] = 1.0e100;
4526    dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100;
4527    strongInfo_[0] = 0;
4528    strongInfo_[1] = 0;
4529    strongInfo_[2] = 0;
4530    strongInfo_[3] = 0;
4531    strongInfo_[4] = 0;
4532    strongInfo_[5] = 0;
4533    strongInfo_[6] = 0;
4534    solverCharacteristics_ = NULL;
4535    nodeCompare_ = new CbcCompareDefault();;
4536    problemFeasibility_ = new CbcFeasibilityBase();
4537    tree_ = new CbcTree();
4538    branchingMethod_ = NULL;
4539    cutModifier_ = NULL;
4540    strategy_ = NULL;
4541    parentModel_ = NULL;
4542    cbcColLower_ = NULL;
4543    cbcColUpper_ = NULL;
4544    cbcRowLower_ = NULL;
4545    cbcRowUpper_ = NULL;
4546    cbcColSolution_ = NULL;
4547    cbcRowPrice_ = NULL;
4548    cbcReducedCost_ = NULL;
4549    cbcRowActivity_ = NULL;
4550    appData_ = NULL;
4551    handler_ = new CoinMessageHandler();
4552    handler_->setLogLevel(2);
4553    messages_ = CbcMessage();
4554    //eventHandler_ = new CbcEventHandler() ;
4555}
4556
4557/** Constructor from solver.
4558
4559  Creates a model complete with a clone of the solver passed as a parameter.
4560*/
4561
4562CbcModel::CbcModel(const OsiSolverInterface &rhs)
4563        :
4564        continuousSolver_(NULL),
4565        referenceSolver_(NULL),
4566        defaultHandler_(true),
4567        emptyWarmStart_(NULL),
4568        bestObjective_(COIN_DBL_MAX),
4569        bestPossibleObjective_(COIN_DBL_MAX),
4570        sumChangeObjective1_(0.0),
4571        sumChangeObjective2_(0.0),
4572        minimumDrop_(1.0e-4),
4573        numberSolutions_(0),
4574        numberSavedSolutions_(0),
4575        maximumSavedSolutions_(0),
4576        stateOfSearch_(0),
4577        whenCuts_(-1),
4578        hotstartSolution_(NULL),
4579        hotstartPriorities_(NULL),
4580        numberHeuristicSolutions_(0),
4581        numberNodes_(0),
4582        numberNodes2_(0),
4583        numberIterations_(0),
4584        numberSolves_(0),
4585        status_(-1),
4586        secondaryStatus_(-1),
4587        numberRowsAtContinuous_(0),
4588        maximumNumberCuts_(0),
4589        phase_(0),
4590        currentNumberCuts_(0),
4591        maximumDepth_(0),
4592        walkback_(NULL),
4593        lastNodeInfo_(NULL),
4594        lastCut_(NULL),
4595        lastDepth_(0),
4596        lastNumberCuts2_(0),
4597        maximumCuts_(0),
4598        lastNumberCuts_(NULL),
4599        addedCuts_(NULL),
4600        nextRowCut_(NULL),
4601        currentNode_(NULL),
4602        integerInfo_(NULL),
4603        specialOptions_(0),
4604        moreSpecialOptions_(0),
4605        subTreeModel_(NULL),
4606        numberStoppedSubTrees_(0),
4607        presolve_(0),
4608        numberStrong_(5),
4609        numberBeforeTrust_(10),
4610        numberPenalties_(20),
4611        stopNumberIterations_(-1),
4612        penaltyScaleFactor_(3.0),
4613        numberAnalyzeIterations_(0),
4614        analyzeResults_(NULL),
4615        numberInfeasibleNodes_(0),
4616        problemType_(0),
4617        printFrequency_(0),
4618        numberCutGenerators_(0),
4619        generator_(NULL),
4620        virginGenerator_(NULL),
4621        numberHeuristics_(0),
4622        heuristic_(NULL),
4623        lastHeuristic_(NULL),
4624# ifdef COIN_HAS_CLP
4625        fastNodeDepth_(-1),
4626#endif
4627        eventHandler_(NULL),
4628        numberObjects_(0),
4629        object_(NULL),
4630        ownObjects_(true),
4631        originalColumns_(NULL),
4632        howOftenGlobalScan_(1),
4633        numberGlobalViolations_(0),
4634        numberExtraIterations_(0),
4635        numberExtraNodes_(0),
4636        continuousObjective_(COIN_DBL_MAX),
4637        originalContinuousObjective_(COIN_DBL_MAX),
4638        continuousInfeasibilities_(COIN_INT_MAX),
4639        maximumCutPassesAtRoot_(20),
4640        maximumCutPasses_(10),
4641        preferredWay_(0),
4642        currentPassNumber_(0),
4643        maximumWhich_(1000),
4644        maximumRows_(0),
4645        currentDepth_(0),
4646        whichGenerator_(NULL),
4647        maximumStatistics_(0),
4648        statistics_(NULL),
4649        maximumDepthActual_(0),
4650        numberDJFixed_(0.0),
4651        probingInfo_(NULL),
4652        numberFixedAtRoot_(0),
4653        numberFixedNow_(0),
4654        stoppedOnGap_(false),
4655        eventHappened_(false),
4656        numberLongStrong_(0),
4657        numberOldActiveCuts_(0),
4658        numberNewCuts_(0),
4659        searchStrategy_(-1),
4660        numberStrongIterations_(0),
4661        resolveAfterTakeOffCuts_(true),
4662        maximumNumberIterations_(-1),
4663        continuousPriority_(COIN_INT_MAX),
4664        numberUpdateItems_(0),
4665        maximumNumberUpdateItems_(0),
4666        updateItems_(NULL),
4667        storedRowCuts_(NULL),
4668        numberThreads_(0),
4669        threadMode_(0),
4670        master_(NULL),
4671        masterThread_(NULL)
4672{
4673    memset(intParam_, 0, sizeof(intParam_));
4674    intParam_[CbcMaxNumNode] = 2147483647;
4675    intParam_[CbcMaxNumSol] = 9999999;
4676
4677    memset(dblParam_, 0, sizeof(dblParam_));
4678    dblParam_[CbcIntegerTolerance] = 1e-6;
4679    dblParam_[CbcCutoffIncrement] = 1e-5;
4680    dblParam_[CbcAllowableGap] = 1.0e-10;
4681    dblParam_[CbcMaximumSeconds] = 1.0e100;
4682    dblParam_[CbcCurrentCutoff] = 1.0e100;
4683    dblParam_[CbcOptimizationDirection] = 1.0;
4684    dblParam_[CbcCurrentObjectiveValue] = 1.0e100;
4685    dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100;
4686    strongInfo_[0] = 0;
4687    strongInfo_[1] = 0;
4688    strongInfo_[2] = 0;
4689    strongInfo_[3] = 0;
4690    strongInfo_[4] = 0;
4691    strongInfo_[5] = 0;
4692    strongInfo_[6] = 0;
4693    solverCharacteristics_ = NULL;
4694    nodeCompare_ = new CbcCompareDefault();;
4695    problemFeasibility_ = new CbcFeasibilityBase();
4696    tree_ = new CbcTree();
4697    branchingMethod_ = NULL;
4698    cutModifier_ = NULL;
4699    strategy_ = NULL;
4700    parentModel_ = NULL;
4701    appData_ = NULL;
4702    solver_ = rhs.clone();
4703    handler_ = new CoinMessageHandler();
4704    if (!solver_->defaultHandler()&&
4705        solver_->messageHandler()->logLevel(0)!=-1000)
4706      passInMessageHandler(solver_->messageHandler());
4707    handler_->setLogLevel(2);
4708    messages_ = CbcMessage();
4709    //eventHandler_ = new CbcEventHandler() ;
4710    referenceSolver_ = solver_->clone();
4711    ownership_ = 0x80000000;
4712    cbcColLower_ = NULL;
4713    cbcColUpper_ = NULL;
4714    cbcRowLower_ = NULL;
4715    cbcRowUpper_ = NULL;
4716    cbcColSolution_ = NULL;
4717    cbcRowPrice_ = NULL;
4718    cbcReducedCost_ = NULL;
4719    cbcRowActivity_ = NULL;
4720
4721    // Initialize solution and integer variable vectors
4722    bestSolution_ = NULL; // to say no solution found
4723    savedSolutions_ = NULL;
4724    numberIntegers_ = 0;
4725    int numberColumns = solver_->getNumCols();
4726    int iColumn;
4727    if (numberColumns) {
4728        // Space for current solution
4729        currentSolution_ = new double[numberColumns];
4730        continuousSolution_ = new double[numberColumns];
4731        usedInSolution_ = new int[numberColumns];
4732        CoinZeroN(usedInSolution_, numberColumns);
4733        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4734            if ( solver_->isInteger(iColumn))
4735                numberIntegers_++;
4736        }
4737    } else {
4738        // empty model
4739        currentSolution_ = NULL;
4740        continuousSolution_ = NULL;
4741        usedInSolution_ = NULL;
4742    }
4743    testSolution_ = currentSolution_;
4744    if (numberIntegers_) {
4745        integerVariable_ = new int [numberIntegers_];
4746        numberIntegers_ = 0;
4747        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4748            if ( solver_->isInteger(iColumn))
4749                integerVariable_[numberIntegers_++] = iColumn;
4750        }
4751    } else {
4752        integerVariable_ = NULL;
4753    }
4754}
4755
4756/*
4757  Assign a solver to the model (model assumes ownership)
4758
4759  The integer variable vector is initialized if it's not already present.
4760  If deleteSolver then current solver deleted (if model owned)
4761
4762  Assuming ownership matches usage in OsiSolverInterface
4763  (cf. assignProblem, loadProblem).
4764
4765  TODO: What to do about solver parameters? A simple copy likely won't do it,
4766        because the SI must push the settings into the underlying solver. In
4767        the context of switching solvers in cbc, this means that command line
4768        settings will get lost. Stash the command line somewhere and reread it
4769        here, maybe?
4770
4771  TODO: More generally, how much state should be transferred from the old
4772        solver to the new solver? Best perhaps to see how usage develops.
4773        What's done here mimics the CbcModel(OsiSolverInterface) constructor.
4774*/
4775void
4776CbcModel::assignSolver(OsiSolverInterface *&solver, bool deleteSolver)
4777
4778{
4779    // resize best solution if exists
4780    if (bestSolution_ && solver && solver_) {
4781        int nOld = solver_->getNumCols();
4782        int nNew = solver->getNumCols();
4783        if (nNew > nOld) {
4784            double * temp = new double[nNew];
4785            memcpy(temp, bestSolution_, nOld*sizeof(double));
4786            memset(temp + nOld, 0, (nNew - nOld)*sizeof(double));
4787            delete [] bestSolution_;
4788            bestSolution_ = temp;
4789        }
4790    }
4791    // Keep the current message level for solver (if solver exists)
4792    if (solver_)
4793        solver->messageHandler()->setLogLevel(solver_->messageHandler()->logLevel()) ;
4794
4795    if (modelOwnsSolver() && deleteSolver) {
4796        solverCharacteristics_ = NULL;
4797        delete solver_ ;
4798    }
4799    solver_ = solver;
4800    solver = NULL ;
4801    setModelOwnsSolver(true) ;
4802    /*
4803      Basis information is solver-specific.
4804    */
4805    if (emptyWarmStart_) {
4806        delete emptyWarmStart_  ;
4807        emptyWarmStart_ = 0 ;
4808    }
4809    bestSolutionBasis_ = CoinWarmStartBasis();
4810    /*
4811      Initialize integer variable vector.
4812    */
4813    numberIntegers_ = 0;
4814    int numberColumns = solver_->getNumCols();
4815    int iColumn;
4816    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4817        if ( solver_->isInteger(iColumn))
4818            numberIntegers_++;
4819    }
4820    delete [] integerVariable_;
4821    if (numberIntegers_) {
4822        integerVariable_ = new int [numberIntegers_];
4823        numberIntegers_ = 0;
4824        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4825            if ( solver_->isInteger(iColumn))
4826                integerVariable_[numberIntegers_++] = iColumn;
4827        }
4828    } else {
4829        integerVariable_ = NULL;
4830    }
4831
4832    return ;
4833}
4834
4835// Copy constructor.
4836
4837CbcModel::CbcModel(const CbcModel & rhs, bool cloneHandler)
4838        :
4839        continuousSolver_(NULL),
4840        referenceSolver_(NULL),
4841        defaultHandler_(rhs.defaultHandler_),
4842        emptyWarmStart_(NULL),
4843        bestObjective_(rhs.bestObjective_),
4844        bestPossibleObjective_(rhs.bestPossibleObjective_),
4845        sumChangeObjective1_(rhs.sumChangeObjective1_),
4846        sumChangeObjective2_(rhs.sumChangeObjective2_),
4847        minimumDrop_(rhs.minimumDrop_),
4848        numberSolutions_(rhs.numberSolutions_),
4849        numberSavedSolutions_(rhs.numberSavedSolutions_),
4850        maximumSavedSolutions_(rhs.maximumSavedSolutions_),
4851        stateOfSearch_(rhs.stateOfSearch_),
4852        whenCuts_(rhs.whenCuts_),
4853        numberHeuristicSolutions_(rhs.numberHeuristicSolutions_),
4854        numberNodes_(rhs.numberNodes_),
4855        numberNodes2_(rhs.numberNodes2_),
4856        numberIterations_(rhs.numberIterations_),
4857        numberSolves_(rhs.numberSolves_),
4858        status_(rhs.status_),
4859        secondaryStatus_(rhs.secondaryStatus_),
4860        specialOptions_(rhs.specialOptions_),
4861        moreSpecialOptions_(rhs.moreSpecialOptions_),
4862        subTreeModel_(rhs.subTreeModel_),
4863        numberStoppedSubTrees_(rhs.numberStoppedSubTrees_),
4864        presolve_(rhs.presolve_),
4865        numberStrong_(rhs.numberStrong_),
4866        numberBeforeTrust_(rhs.numberBeforeTrust_),
4867        numberPenalties_(rhs.numberPenalties_),
4868        stopNumberIterations_(rhs.stopNumberIterations_),
4869        penaltyScaleFactor_(rhs.penaltyScaleFactor_),
4870        numberAnalyzeIterations_(rhs.numberAnalyzeIterations_),
4871        analyzeResults_(NULL),
4872        numberInfeasibleNodes_(rhs.numberInfeasibleNodes_),
4873        problemType_(rhs.problemType_),
4874        printFrequency_(rhs.printFrequency_),
4875# ifdef COIN_HAS_CLP
4876        fastNodeDepth_(rhs.fastNodeDepth_),
4877#endif
4878        howOftenGlobalScan_(rhs.howOftenGlobalScan_),
4879        numberGlobalViolations_(rhs.numberGlobalViolations_),
4880        numberExtraIterations_(rhs.numberExtraIterations_),
4881        numberExtraNodes_(rhs.numberExtraNodes_),
4882        continuousObjective_(rhs.continuousObjective_),
4883        originalContinuousObjective_(rhs.originalContinuousObjective_),
4884        continuousInfeasibilities_(rhs.continuousInfeasibilities_),
4885        maximumCutPassesAtRoot_(rhs.maximumCutPassesAtRoot_),
4886        maximumCutPasses_( rhs.maximumCutPasses_),
4887        preferredWay_(rhs.preferredWay_),
4888        currentPassNumber_(rhs.currentPassNumber_),
4889        maximumWhich_(rhs.maximumWhich_),
4890        maximumRows_(0),
4891        currentDepth_(0),
4892        whichGenerator_(NULL),
4893        maximumStatistics_(0),
4894        statistics_(NULL),
4895        maximumDepthActual_(0),
4896        numberDJFixed_(0.0),
4897        probingInfo_(NULL),
4898        numberFixedAtRoot_(rhs.numberFixedAtRoot_),
4899        numberFixedNow_(rhs.numberFixedNow_),
4900        stoppedOnGap_(rhs.stoppedOnGap_),
4901        eventHappened_(rhs.eventHappened_),
4902        numberLongStrong_(rhs.numberLongStrong_),
4903        numberOldActiveCuts_(rhs.numberOldActiveCuts_),
4904        numberNewCuts_(rhs.numberNewCuts_),
4905        searchStrategy_(rhs.searchStrategy_),
4906        numberStrongIterations_(rhs.numberStrongIterations_),
4907        resolveAfterTakeOffCuts_(rhs.resolveAfterTakeOffCuts_),
4908        maximumNumberIterations_(rhs.maximumNumberIterations_),
4909        continuousPriority_(rhs.continuousPriority_),
4910        numberUpdateItems_(rhs.numberUpdateItems_),
4911        maximumNumberUpdateItems_(rhs.maximumNumberUpdateItems_),
4912        updateItems_(NULL),
4913        storedRowCuts_(NULL),
4914        numberThreads_(rhs.numberThreads_),
4915        threadMode_(rhs.threadMode_),
4916        master_(NULL),
4917        masterThread_(NULL)
4918{
4919    memcpy(intParam_, rhs.intParam_, sizeof(intParam_));
4920    memcpy(dblParam_, rhs.dblParam_, sizeof(dblParam_));
4921    strongInfo_[0] = rhs.strongInfo_[0];
4922    strongInfo_[1] = rhs.strongInfo_[1];
4923    strongInfo_[2] = rhs.strongInfo_[2];
4924    strongInfo_[3] = rhs.strongInfo_[3];
4925    strongInfo_[4] = rhs.strongInfo_[4];
4926    strongInfo_[5] = rhs.strongInfo_[5];
4927    strongInfo_[6] = rhs.strongInfo_[6];
4928    solverCharacteristics_ = NULL;
4929    if (rhs.emptyWarmStart_) emptyWarmStart_ = rhs.emptyWarmStart_->clone() ;
4930    if (defaultHandler_ || cloneHandler) {
4931        handler_ = new CoinMessageHandler();
4932        handler_->setLogLevel(2);
4933    } else {
4934        handler_ = rhs.handler_;
4935    }
4936    messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
4937    numberCutGenerators_ = rhs.numberCutGenerators_;
4938    if (numberCutGenerators_) {
4939        generator_ = new CbcCutGenerator * [numberCutGenerators_];
4940        virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_];
4941        int i;
4942        for (i = 0; i < numberCutGenerators_; i++) {
4943            generator_[i] = new CbcCutGenerator(*rhs.generator_[i]);
4944            virginGenerator_[i] = new CbcCutGenerator(*rhs.virginGenerator_[i]);
4945        }
4946    } else {
4947        generator_ = NULL;
4948        virginGenerator_ = NULL;
4949    }
4950    globalCuts_ = rhs.globalCuts_;
4951    numberHeuristics_ = rhs.numberHeuristics_;
4952    if (numberHeuristics_) {
4953        heuristic_ = new CbcHeuristic * [numberHeuristics_];
4954        int i;
4955        for (i = 0; i < numberHeuristics_; i++) {
4956            heuristic_[i] = rhs.heuristic_[i]->clone();
4957        }
4958    } else {
4959        heuristic_ = NULL;
4960    }
4961    lastHeuristic_ = NULL;
4962    if (rhs.eventHandler_) {
4963        eventHandler_ = rhs.eventHandler_->clone() ;
4964    } else {
4965        eventHandler_ = NULL ;
4966    }
4967    ownObjects_ = rhs.ownObjects_;
4968    if (ownObjects_) {
4969        numberObjects_ = rhs.numberObjects_;
4970        if (numberObjects_) {
4971            object_ = new OsiObject * [numberObjects_];
4972            int i;
4973            for (i = 0; i < numberObjects_; i++) {
4974                object_[i] = (rhs.object_[i])->clone();
4975                CbcObject * obj = dynamic_cast <CbcObject *>(object_[i]) ;
4976                // Could be OsiObjects
4977                if (obj)
4978                    obj->setModel(this);
4979            }
4980        } else {
4981            object_ = NULL;
4982        }
4983    } else {
4984        // assume will be redone
4985        numberObjects_ = 0;
4986        object_ = NULL;
4987    }
4988    if (rhs.referenceSolver_)
4989        referenceSolver_ = rhs.referenceSolver_->clone();
4990    else
4991        referenceSolver_ = NULL;
4992    solver_ = rhs.solver_->clone();
4993    if (rhs.originalColumns_) {
4994        int numberColumns = solver_->getNumCols();
4995        originalColumns_ = new int [numberColumns];
4996        memcpy(originalColumns_, rhs.originalColumns_, numberColumns*sizeof(int));
4997    } else {
4998        originalColumns_ = NULL;
4999    }
5000    if (maximumNumberUpdateItems_) {
5001        updateItems_ = new CbcObjectUpdateData [maximumNumberUpdateItems_];
5002        for (int i = 0; i < maximumNumberUpdateItems_; i++)
5003            updateItems_[i] = rhs.updateItems_[i];
5004    }
5005    if (maximumWhich_ && rhs.whichGenerator_)
5006        whichGenerator_ = CoinCopyOfArray(rhs.whichGenerator_, maximumWhich_);
5007    nodeCompare_ = rhs.nodeCompare_->clone();
5008    problemFeasibility_ = rhs.problemFeasibility_->clone();
5009    tree_ = rhs.tree_->clone();
5010    if (rhs.branchingMethod_)
5011        branchingMethod_ = rhs.branchingMethod_->clone();
5012    else
5013        branchingMethod_ = NULL;
5014    if (rhs.cutModifier_)
5015        cutModifier_ = rhs.cutModifier_->clone();
5016    else
5017        cutModifier_ = NULL;
5018    cbcColLower_ = NULL;
5019    cbcColUpper_ = NULL;
5020    cbcRowLower_ = NULL;
5021    cbcRowUpper_ = NULL;
5022    cbcColSolution_ = NULL;
5023    cbcRowPrice_ = NULL;
5024    cbcReducedCost_ = NULL;
5025    cbcRowActivity_ = NULL;
5026    if (rhs.strategy_)
5027        strategy_ = rhs.strategy_->clone();
5028    else
5029        strategy_ = NULL;
5030    parentModel_ = rhs.parentModel_;
5031    appData_ = rhs.appData_;
5032    messages_ = rhs.messages_;
5033    ownership_ = rhs.ownership_ | 0x80000000;
5034    messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
5035    numberIntegers_ = rhs.numberIntegers_;
5036    randomNumberGenerator_ = rhs.randomNumberGenerator_;
5037    if (numberIntegers_) {
5038        integerVariable_ = new int [numberIntegers_];
5039        memcpy(integerVariable_, rhs.integerVariable_, numberIntegers_*sizeof(int));
5040        integerInfo_ = CoinCopyOfArray(rhs.integerInfo_, solver_->getNumCols());
5041    } else {
5042        integerVariable_ = NULL;
5043        integerInfo_ = NULL;
5044    }
5045    if (rhs.hotstartSolution_) {
5046        int numberColumns = solver_->getNumCols();
5047        hotstartSolution_ = CoinCopyOfArray(rhs.hotstartSolution_, numberColumns);
5048        hotstartPriorities_ = CoinCopyOfArray(rhs.hotstartPriorities_, numberColumns);
5049    } else {
5050        hotstartSolution_ = NULL;
5051        hotstartPriorities_ = NULL;
5052    }
5053    if (rhs.bestSolution_) {
5054        int numberColumns = solver_->getNumCols();
5055        bestSolution_ = new double[numberColumns];
5056        memcpy(bestSolution_, rhs.bestSolution_, numberColumns*sizeof(double));
5057    } else {
5058        bestSolution_ = NULL;
5059    }
5060    int numberColumns = solver_->getNumCols();
5061    if (maximumSavedSolutions_ && rhs.savedSolutions_) {
5062        savedSolutions_ = new double * [maximumSavedSolutions_];
5063        for (int i = 0; i < maximumSavedSolutions_; i++)
5064            savedSolutions_[i] = CoinCopyOfArray(rhs.savedSolutions_[i], numberColumns + 2);
5065    } else {
5066        savedSolutions_ = NULL;
5067    }
5068    // Space for current solution
5069    currentSolution_ = new double[numberColumns];
5070    continuousSolution_ = new double[numberColumns];
5071    usedInSolution_ = new int[numberColumns];
5072    CoinZeroN(usedInSolution_, numberColumns);
5073    testSolution_ = currentSolution_;
5074    numberRowsAtContinuous_ = rhs.numberRowsAtContinuous_;
5075    maximumNumberCuts_ = rhs.maximumNumberCuts_;
5076    phase_ = rhs.phase_;
5077    currentNumberCuts_ = rhs.currentNumberCuts_;
5078    maximumDepth_ = rhs.maximumDepth_;
5079    // These are only used as temporary arrays so need not be filled
5080    if (maximumNumberCuts_) {
5081        addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
5082    } else {
5083        addedCuts_ = NULL;
5084    }
5085    bestSolutionBasis_ = rhs.bestSolutionBasis_;
5086    nextRowCut_ = NULL;
5087    currentNode_ = NULL;
5088    if (maximumDepth_) {
5089        walkback_ = new CbcNodeInfo * [maximumDepth_];
5090        lastNodeInfo_ = new CbcNodeInfo * [maximumDepth_] ;
5091        lastNumberCuts_ = new int [maximumDepth_] ;
5092    } else {
5093        walkback_ = NULL;
5094        lastNodeInfo_ = NULL;
5095        lastNumberCuts_ = NULL;
5096    }
5097    maximumCuts_ = rhs.maximumCuts_;
5098    if (maximumCuts_) {
5099        lastCut_ = new const OsiRowCut * [maximumCuts_] ;
5100    } else {
5101        lastCut_ = NULL;
5102    }
5103    synchronizeModel();
5104    if (cloneHandler && !defaultHandler_) {
5105        delete handler_;
5106        CoinMessageHandler * handler = rhs.handler_->clone();
5107        passInMessageHandler(handler);
5108    }
5109}
5110
5111// Assignment operator
5112CbcModel &
5113CbcModel::operator=(const CbcModel & rhs)
5114{
5115    if (this != &rhs) {
5116        if (modelOwnsSolver()) {
5117            solverCharacteristics_ = NULL;
5118            delete solver_;
5119            solver_ = NULL;
5120        }
5121        gutsOfDestructor();
5122        if (defaultHandler_) {
5123            delete handler_;
5124            handler_ = NULL;
5125        }
5126        defaultHandler_ = rhs.defaultHandler_;
5127        if (defaultHandler_) {
5128            handler_ = new CoinMessageHandler();
5129            handler_->setLogLevel(2);
5130        } else {
5131            handler_ = rhs.handler_;
5132        }
5133        messages_ = rhs.messages_;
5134        messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
5135        if (rhs.solver_) {
5136            solver_ = rhs.solver_->clone() ;
5137        } else {
5138            solver_ = 0 ;
5139        }
5140        ownership_ = 0x80000000;
5141        delete continuousSolver_ ;
5142        if (rhs.continuousSolver_) {
5143            continuousSolver_ = rhs.continuousSolver_->clone() ;
5144        } else {
5145            continuousSolver_ = 0 ;
5146        }
5147        delete referenceSolver_;
5148        if (rhs.referenceSolver_) {
5149            referenceSolver_ = rhs.referenceSolver_->clone() ;
5150        } else {
5151            referenceSolver_ = NULL ;
5152        }
5153
5154        delete emptyWarmStart_ ;
5155        if (rhs.emptyWarmStart_) {
5156            emptyWarmStart_ = rhs.emptyWarmStart_->clone() ;
5157        } else {
5158            emptyWarmStart_ = 0 ;
5159        }
5160
5161        bestObjective_ = rhs.bestObjective_;
5162        bestPossibleObjective_ = rhs.bestPossibleObjective_;
5163        sumChangeObjective1_ = rhs.sumChangeObjective1_;
5164        sumChangeObjective2_ = rhs.sumChangeObjective2_;
5165        delete [] bestSolution_;
5166        if (rhs.bestSolution_) {
5167            int numberColumns = rhs.getNumCols();
5168            bestSolution_ = new double[numberColumns];
5169            memcpy(bestSolution_, rhs.bestSolution_, numberColumns*sizeof(double));
5170        } else {
5171            bestSolution_ = NULL;
5172        }
5173        for (int i = 0; i < maximumSavedSolutions_; i++)
5174            delete [] savedSolutions_[i];
5175        delete [] savedSolutions_;
5176        savedSolutions_ = NULL;
5177        int numberColumns = rhs.getNumCols();
5178        if (numberColumns) {
5179            // Space for current solution
5180            currentSolution_ = new double[numberColumns];
5181            continuousSolution_ = new double[numberColumns];
5182            usedInSolution_ = new int[numberColumns];
5183            CoinZeroN(usedInSolution_, numberColumns);
5184        } else {
5185            currentSolution_ = NULL;
5186            continuousSolution_ = NULL;
5187            usedInSolution_ = NULL;
5188        }
5189        if (maximumSavedSolutions_) {
5190            savedSolutions_ = new double * [maximumSavedSolutions_];
5191            for (int i = 0; i < maximumSavedSolutions_; i++)
5192                savedSolutions_[i] = CoinCopyOfArray(rhs.savedSolutions_[i], numberColumns + 2);
5193        } else {
5194            savedSolutions_ = NULL;
5195        }
5196        testSolution_ = currentSolution_;
5197        minimumDrop_ = rhs.minimumDrop_;
5198        numberSolutions_ = rhs.numberSolutions_;
5199        numberSavedSolutions_ = rhs.numberSavedSolutions_;
5200        maximumSavedSolutions_ = rhs.maximumSavedSolutions_;
5201        stateOfSearch_ = rhs.stateOfSearch_;
5202        whenCuts_ = rhs.whenCuts_;
5203        numberHeuristicSolutions_ = rhs.numberHeuristicSolutions_;
5204        numberNodes_ = rhs.numberNodes_;
5205        numberNodes2_ = rhs.numberNodes2_;
5206        numberIterations_ = rhs.numberIterations_;
5207        numberSolves_ = rhs.numberSolves_;
5208        status_ = rhs.status_;
5209        secondaryStatus_ = rhs.secondaryStatus_;
5210        specialOptions_ = rhs.specialOptions_;
5211        moreSpecialOptions_ = rhs.moreSpecialOptions_;
5212        subTreeModel_ = rhs.subTreeModel_;
5213        numberStoppedSubTrees_ = rhs.numberStoppedSubTrees_;
5214        presolve_ = rhs.presolve_;
5215        numberStrong_ = rhs.numberStrong_;
5216        numberBeforeTrust_ = rhs.numberBeforeTrust_;
5217        numberPenalties_ = rhs.numberPenalties_;
5218        stopNumberIterations_ = rhs.stopNumberIterations_;
5219        penaltyScaleFactor_ = rhs.penaltyScaleFactor_;
5220        numberAnalyzeIterations_ = rhs.numberAnalyzeIterations_;
5221        delete [] analyzeResults_;
5222        analyzeResults_ = NULL;
5223        numberInfeasibleNodes_ = rhs.numberInfeasibleNodes_;
5224        problemType_ = rhs.problemType_;
5225        printFrequency_ = rhs.printFrequency_;
5226        howOftenGlobalScan_ = rhs.howOftenGlobalScan_;
5227        numberGlobalViolations_ = rhs.numberGlobalViolations_;
5228        numberExtraIterations_ = rhs.numberExtraIterations_;
5229        numberExtraNodes_ = rhs.numberExtraNodes_;
5230        continuousObjective_ = rhs.continuousObjective_;
5231        originalContinuousObjective_ = rhs.originalContinuousObjective_;
5232        continuousInfeasibilities_ = rhs.continuousInfeasibilities_;
5233        maximumCutPassesAtRoot_ = rhs.maximumCutPassesAtRoot_;
5234        maximumCutPasses_ = rhs.maximumCutPasses_;
5235        preferredWay_ = rhs.preferredWay_;
5236        currentPassNumber_ = rhs.currentPassNumber_;
5237        memcpy(intParam_, rhs.intParam_, sizeof(intParam_));
5238        memcpy(dblParam_, rhs.dblParam_, sizeof(dblParam_));
5239        globalCuts_ = rhs.globalCuts_;
5240        int i;
5241        for (i = 0; i < numberCutGenerators_; i++) {
5242            delete generator_[i];
5243            delete virginGenerator_[i];
5244        }
5245        delete [] generator_;
5246        delete [] virginGenerator_;
5247        delete [] heuristic_;
5248        maximumWhich_ = rhs.maximumWhich_;
5249        delete [] whichGenerator_;
5250        whichGenerator_ = NULL;
5251        if (maximumWhich_ && rhs.whichGenerator_)
5252            whichGenerator_ = CoinCopyOfArray(rhs.whichGenerator_, maximumWhich_);
5253        maximumRows_ = 0;
5254        currentDepth_ = 0;
5255        randomNumberGenerator_ = rhs.randomNumberGenerator_;
5256        workingBasis_ = CoinWarmStartBasis();
5257        for (i = 0; i < maximumStatistics_; i++)
5258            delete statistics_[i];
5259        delete [] statistics_;
5260        maximumStatistics_ = 0;
5261        statistics_ = NULL;
5262        delete probingInfo_;
5263        probingInfo_ = NULL;
5264        numberFixedAtRoot_ = rhs.numberFixedAtRoot_;
5265        numberFixedNow_ = rhs.numberFixedNow_;
5266        stoppedOnGap_ = rhs.stoppedOnGap_;
5267        eventHappened_ = rhs.eventHappened_;
5268        numberLongStrong_ = rhs.numberLongStrong_;
5269        numberOldActiveCuts_ = rhs.numberOldActiveCuts_;
5270        numberNewCuts_ = rhs.numberNewCuts_;
5271        resolveAfterTakeOffCuts_ = rhs.resolveAfterTakeOffCuts_;
5272        maximumNumberIterations_ = rhs.maximumNumberIterations_;
5273        continuousPriority_ = rhs.continuousPriority_;
5274        numberUpdateItems_ = rhs.numberUpdateItems_;
5275        maximumNumberUpdateItems_ = rhs.maximumNumberUpdateItems_;
5276        delete [] updateItems_;
5277        if (maximumNumberUpdateItems_) {
5278            updateItems_ = new CbcObjectUpdateData [maximumNumberUpdateItems_];
5279            for (i = 0; i < maximumNumberUpdateItems_; i++)
5280                updateItems_[i] = rhs.updateItems_[i];
5281        } else {
5282            updateItems_ = NULL;
5283        }
5284        numberThreads_ = rhs.numberThreads_;
5285        threadMode_ = rhs.threadMode_;
5286        delete master_;
5287        master_ = NULL;
5288        masterThread_ = NULL;
5289        searchStrategy_ = rhs.searchStrategy_;
5290        numberStrongIterations_ = rhs.numberStrongIterations_;
5291        strongInfo_[0] = rhs.strongInfo_[0];
5292        strongInfo_[1] = rhs.strongInfo_[1];
5293        strongInfo_[2] = rhs.strongInfo_[2];
5294        strongInfo_[3] = rhs.strongInfo_[3];
5295        strongInfo_[4] = rhs.strongInfo_[4];
5296        strongInfo_[5] = rhs.strongInfo_[5];
5297        strongInfo_[6] = rhs.strongInfo_[6];
5298        solverCharacteristics_ = NULL;
5299        lastHeuristic_ = NULL;
5300        numberCutGenerators_ = rhs.numberCutGenerators_;
5301        if (numberCutGenerators_) {
5302            generator_ = new CbcCutGenerator * [numberCutGenerators_];
5303            virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_];
5304            int i;
5305            for (i = 0; i < numberCutGenerators_; i++) {
5306                generator_[i] = new CbcCutGenerator(*rhs.generator_[i]);
5307                virginGenerator_[i] = new CbcCutGenerator(*rhs.virginGenerator_[i]);
5308            }
5309        } else {
5310            generator_ = NULL;
5311            virginGenerator_ = NULL;
5312        }
5313        numberHeuristics_ = rhs.numberHeuristics_;
5314        if (numberHeuristics_) {
5315            heuristic_ = new CbcHeuristic * [numberHeuristics_];
5316            memcpy(heuristic_, rhs.heuristic_,
5317                   numberHeuristics_*sizeof(CbcHeuristic *));
5318        } else {
5319            heuristic_ = NULL;
5320        }
5321        lastHeuristic_ = NULL;
5322        if (eventHandler_)
5323            delete eventHandler_ ;
5324        if (rhs.eventHandler_) {
5325            eventHandler_ = rhs.eventHandler_->clone() ;
5326        } else {
5327            eventHandler_ = NULL ;
5328        }
5329# ifdef COIN_HAS_CLP
5330        fastNodeDepth_ = rhs.fastNodeDepth_;
5331#endif
5332        if (ownObjects_) {
5333            for (i = 0; i < numberObjects_; i++)
5334                delete object_[i];
5335            delete [] object_;
5336            numberObjects_ = rhs.numberObjects_;
5337            if (numberObjects_) {
5338                object_ = new OsiObject * [numberObjects_];
5339                int i;
5340                for (i = 0; i < numberObjects_; i++)
5341                    object_[i] = (rhs.object_[i])->clone();
5342            } else {
5343                object_ = NULL;
5344            }
5345        } else {
5346            // assume will be redone
5347            numberObjects_ = 0;
5348            object_ = NULL;
5349        }
5350        delete [] originalColumns_;
5351        if (rhs.originalColumns_) {
5352            int numberColumns = rhs.getNumCols();
5353            originalColumns_ = new int [numberColumns];
5354            memcpy(originalColumns_, rhs.originalColumns_, numberColumns*sizeof(int));
5355        } else {
5356            originalColumns_ = NULL;
5357        }
5358        nodeCompare_ = rhs.nodeCompare_->clone();
5359        problemFeasibility_ = rhs.problemFeasibility_->clone();
5360        delete tree_;
5361        tree_ = rhs.tree_->clone();
5362        if (rhs.branchingMethod_)
5363            branchingMethod_ = rhs.branchingMethod_->clone();
5364        else
5365            branchingMethod_ = NULL;
5366        if (rhs.cutModifier_)
5367            cutModifier_ = rhs.cutModifier_->clone();
5368        else
5369            cutModifier_ = NULL;
5370        delete strategy_;
5371        if (rhs.strategy_)
5372            strategy_ = rhs.strategy_->clone();
5373        else
5374            strategy_ = NULL;
5375        parentModel_ = rhs.parentModel_;
5376        appData_ = rhs.appData_;
5377
5378        delete [] integerVariable_;
5379        numberIntegers_ = rhs.numberIntegers_;
5380        if (numberIntegers_) {
5381            integerVariable_ = new int [numberIntegers_];
5382            memcpy(integerVariable_, rhs.integerVariable_,
5383                   numbe