source: branches/sandbox/Cbc/src/CbcModel.cpp @ 1357

Last change on this file since 1357 was 1357, checked in by coin, 12 years ago

run 'astyle -A4 -p' and dos2unix

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