source: stable/2.6/Cbc/src/CbcModel.cpp @ 1653

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

allow cleaner exit when integers all preprocessed away (and compiler error in gcc 4.6)

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