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

Last change on this file since 1332 was 1332, checked in by EdwinStraver, 11 years ago

Further breakouts from branchAndBound()

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