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

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

add Proximity heuristic (Fischetti and Monaci) - shouldn't break anything

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