source: stable/2.8/Cbc/src/CbcModel.cpp @ 1910

Last change on this file since 1910 was 1910, checked in by tkr, 6 years ago

Merging r1909 from trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 715.0 KB
RevLine 
[1271]1/* $Id: CbcModel.cpp 1910 2013-04-25 19:46:48Z tkr $ */
[2]2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
[1573]4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
[2]6#if defined(_MSC_VER)
7// Turn off compiler warning about long names
8#  pragma warning(disable:4786)
9#endif
[311]10
[325]11#include "CbcConfig.h"
[311]12
[2]13#include <string>
14//#define CBC_DEBUG 1
15//#define CHECK_CUT_COUNTS
[1412]16//#define CHECK_NODE
[2]17//#define CHECK_NODE_FULL
[264]18//#define NODE_LOG
[640]19//#define GLOBAL_CUTS_JUST_POINTERS
[1271]20#ifdef CGL_DEBUG_GOMORY
21extern int gomory_try;
[687]22#endif
[2]23#include <cassert>
24#include <cmath>
25#include <cfloat>
[311]26#ifdef COIN_HAS_CLP
[208]27// include Presolve from Clp
28#include "ClpPresolve.hpp"
29#include "OsiClpSolverInterface.hpp"
[931]30#include "ClpNode.hpp"
31#include "ClpDualRowDantzig.hpp"
32#include "ClpSimplexPrimal.hpp"
[277]33#endif
34
35#include "CbcEventHandler.hpp"
[2]36
37#include "OsiSolverInterface.hpp"
[264]38#include "OsiAuxInfo.hpp"
[222]39#include "OsiSolverBranch.hpp"
[640]40#include "OsiChooseVariable.hpp"
[2]41#include "CoinWarmStartBasis.hpp"
42#include "CoinPackedMatrix.hpp"
[140]43#include "CoinHelperFunctions.hpp"
[2]44#include "CbcBranchActual.hpp"
[135]45#include "CbcBranchDynamic.hpp"
[2]46#include "CbcHeuristic.hpp"
[687]47#include "CbcHeuristicFPump.hpp"
[1315]48#include "CbcHeuristicRINS.hpp"
[1271]49#include "CbcHeuristicDive.hpp"
[2]50#include "CbcModel.hpp"
[186]51#include "CbcTreeLocal.hpp"
[150]52#include "CbcStatistics.hpp"
[97]53#include "CbcStrategy.hpp"
[2]54#include "CbcMessage.hpp"
55#include "OsiRowCut.hpp"
56#include "OsiColCut.hpp"
57#include "OsiRowCutDebugger.hpp"
58#include "OsiCuts.hpp"
59#include "CbcCountRowCut.hpp"
60#include "CbcCutGenerator.hpp"
[165]61#include "CbcFeasibilityBase.hpp"
[687]62#include "CbcFathom.hpp"
[1839]63#include "CbcFullNodeInfo.hpp"
[2]64// include Probing
65#include "CglProbing.hpp"
[1271]66#include "CglGomory.hpp"
67#include "CglTwomir.hpp"
[222]68// include preprocessing
69#include "CglPreProcess.hpp"
[640]70#include "CglDuplicateRow.hpp"
71#include "CglStored.hpp"
[833]72#include "CglClique.hpp"
[2]73
74#include "CoinTime.hpp"
[640]75#include "CoinMpsIO.hpp"
[2]76
77#include "CbcCompareActual.hpp"
78#include "CbcTree.hpp"
[1409]79// This may be dummy
80#include "CbcThread.hpp"
[2]81/* Various functions local to CbcModel.cpp */
82
[1844]83static void * doRootCbcThread(void * voidInfo);
84
[2]85namespace {
86
87//-------------------------------------------------------------------
[1286]88// Returns the greatest common denominator of two
89// positive integers, a and b, found using Euclid's algorithm
[2]90//-------------------------------------------------------------------
[1286]91static int gcd(int a, int b)
[2]92{
[1286]93    int remainder = -1;
94    // make sure a<=b (will always remain so)
95    if (a > b) {
96        // Swap a and b
97        int temp = a;
98        a = b;
99        b = temp;
[2]100    }
[1286]101    // if zero then gcd is nonzero (zero may occur in rhs of packed)
102    if (!a) {
103        if (b) {
104            return b;
105        } else {
106            printf("**** gcd given two zeros!!\n");
107            abort();
108        }
109    }
110    while (remainder) {
111        remainder = b % a;
112        b = a;
113        a = remainder;
114    }
115    return b;
[2]116}
117
[931]118
[2]119
120#ifdef CHECK_NODE_FULL
121
122/*
123  Routine to verify that tree linkage is correct. The invariant that is tested
124  is
125
126  reference count = (number of actual references) + (number of branches left)
127
128  The routine builds a set of paired arrays, info and count, by traversing the
129  tree. Each CbcNodeInfo is recorded in info, and the number of times it is
130  referenced (via the parent field) is recorded in count. Then a final check is
131  made to see if the numberPointingToThis_ field agrees.
132*/
133
134void verifyTreeNodes (const CbcTree * branchingTree, const CbcModel &model)
135
[1286]136{
137    if (model.getNodeCount() == 661) return;
138    printf("*** CHECKING tree after %d nodes\n", model.getNodeCount()) ;
139
140    int j ;
141    int nNodes = branchingTree->size() ;
[2]142# define MAXINFO 1000
[1286]143    int *count = new int [MAXINFO] ;
144    CbcNodeInfo **info = new CbcNodeInfo*[MAXINFO] ;
145    int nInfo = 0 ;
146    /*
147      Collect all CbcNodeInfo objects in info, by starting from each live node and
148      traversing back to the root. Nodes in the live set should have unexplored
149      branches remaining.
[2]150
[1286]151      TODO: The `while (nodeInfo)' loop could be made to break on reaching a
152        common ancester (nodeInfo is found in info[k]). Alternatively, the
153        check could change to signal an error if nodeInfo is not found above a
154        common ancestor.
155    */
156    for (j = 0 ; j < nNodes ; j++) {
157        CbcNode *node = branchingTree->nodePointer(j) ;
158        if (!node)
159            continue;
160        CbcNodeInfo *nodeInfo = node->nodeInfo() ;
161        int change = node->nodeInfo()->numberBranchesLeft() ;
162        assert(change) ;
163        while (nodeInfo) {
164            int k ;
165            for (k = 0 ; k < nInfo ; k++) {
166                if (nodeInfo == info[k]) break ;
167            }
168            if (k == nInfo) {
169                assert(nInfo < MAXINFO) ;
170                nInfo++ ;
171                info[k] = nodeInfo ;
172                count[k] = 0 ;
173            }
174            nodeInfo = nodeInfo->parent() ;
175        }
176    }
177    /*
178      Walk the info array. For each nodeInfo, look up its parent in info and
179      increment the corresponding count.
180    */
181    for (j = 0 ; j < nInfo ; j++) {
182        CbcNodeInfo *nodeInfo = info[j] ;
183        nodeInfo = nodeInfo->parent() ;
184        if (nodeInfo) {
185            int k ;
186            for (k = 0 ; k < nInfo ; k++) {
187                if (nodeInfo == info[k]) break ;
188            }
189            assert (k < nInfo) ;
190            count[k]++ ;
191        }
192    }
193    /*
194      Walk the info array one more time and check that the invariant holds. The
195      number of references (numberPointingToThis()) should equal the sum of the
196      number of actual references (held in count[]) plus the number of potential
197      references (unexplored branches, numberBranchesLeft()).
198    */
199    for (j = 0; j < nInfo; j++) {
200        CbcNodeInfo * nodeInfo = info[j] ;
201        if (nodeInfo) {
202            int k ;
203            for (k = 0; k < nInfo; k++)
204                if (nodeInfo == info[k])
205                    break ;
206            printf("Nodeinfo %x - %d left, %d count\n",
207                   nodeInfo,
208                   nodeInfo->numberBranchesLeft(),
209                   nodeInfo->numberPointingToThis()) ;
210            assert(nodeInfo->numberPointingToThis() ==
211                   count[k] + nodeInfo->numberBranchesLeft()) ;
212        }
213    }
[2]214
[1286]215    delete [] count ;
216    delete [] info ;
[2]217
[1286]218    return ;
219}
220
[36]221#endif  /* CHECK_NODE_FULL */
[2]222
[931]223
[2]224
225#ifdef CHECK_CUT_COUNTS
226
227/*
228  Routine to verify that cut reference counts are correct.
229*/
230void verifyCutCounts (const CbcTree * branchingTree, CbcModel &model)
231
[1286]232{
233    printf("*** CHECKING cuts after %d nodes\n", model.getNodeCount()) ;
[2]234
[1286]235    int j ;
236    int nNodes = branchingTree->size() ;
[2]237
[1286]238    /*
239      cut.tempNumber_ exists for the purpose of doing this verification. Clear it
240      in all cuts. We traverse the tree by starting from each live node and working
241      back to the root. At each CbcNodeInfo, check for cuts.
242    */
243    for (j = 0 ; j < nNodes ; j++) {
244        CbcNode *node = branchingTree->nodePointer(j) ;
245        CbcNodeInfo * nodeInfo = node->nodeInfo() ;
246        assert (node->nodeInfo()->numberBranchesLeft()) ;
247        while (nodeInfo) {
248            int k ;
249            for (k = 0 ; k < nodeInfo->numberCuts() ; k++) {
250                CbcCountRowCut *cut = nodeInfo->cuts()[k] ;
251                if (cut) cut->tempNumber_ = 0;
252            }
253            nodeInfo = nodeInfo->parent() ;
254        }
255    }
256    /*
257      Walk the live set again, this time collecting the list of cuts in use at each
258      node. addCuts1 will collect the cuts in model.addedCuts_. Take into account
259      that when we recreate the basis for a node, we compress out the slack cuts.
260    */
261    for (j = 0 ; j < nNodes ; j++) {
262        CoinWarmStartBasis *debugws = model.getEmptyBasis() ;
263        CbcNode *node = branchingTree->nodePointer(j) ;
264        CbcNodeInfo *nodeInfo = node->nodeInfo();
265        int change = node->nodeInfo()->numberBranchesLeft() ;
266        printf("Node %d %x (info %x) var %d way %d obj %g", j, node,
267               node->nodeInfo(), node->columnNumber(), node->way(),
268               node->objectiveValue()) ;
[2]269
[1286]270        model.addCuts1(node, debugws) ;
[2]271
[1286]272        int i ;
273        int numberRowsAtContinuous = model.numberRowsAtContinuous() ;
274        CbcCountRowCut **addedCuts = model.addedCuts() ;
275        for (i = 0 ; i < model.currentNumberCuts() ; i++) {
276            CoinWarmStartBasis::Status status =
277                debugws->getArtifStatus(i + numberRowsAtContinuous) ;
278            if (status != CoinWarmStartBasis::basic && addedCuts[i]) {
279                addedCuts[i]->tempNumber_ += change ;
280            }
281        }
[2]282
[1286]283        while (nodeInfo) {
284            nodeInfo = nodeInfo->parent() ;
285            if (nodeInfo) printf(" -> %x", nodeInfo);
286        }
287        printf("\n") ;
288        delete debugws ;
289    }
290    /*
291      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.
[2]292
[1286]293      TODO: Rewrite to check and print mismatch only when tempNumber_ == 0?
294    */
295    for (j = 0 ; j < nNodes ; j++) {
296        CbcNode *node = branchingTree->nodePointer(j) ;
297        CbcNodeInfo *nodeInfo = node->nodeInfo();
298        while (nodeInfo) {
299            int k ;
300            for (k = 0 ; k < nodeInfo->numberCuts() ; k++) {
301                CbcCountRowCut *cut = nodeInfo->cuts()[k] ;
302                if (cut && cut->tempNumber_ >= 0) {
303                    if (cut->tempNumber_ != cut->numberPointingToThis())
304                        printf("mismatch %x %d %x %d %d\n", nodeInfo, k,
305                               cut, cut->tempNumber_, cut->numberPointingToThis()) ;
306                    else
307                        printf("   match %x %d %x %d %d\n", nodeInfo, k,
308                               cut, cut->tempNumber_, cut->numberPointingToThis()) ;
309                    cut->tempNumber_ = -1 ;
310                }
311            }
312            nodeInfo = nodeInfo->parent() ;
313        }
314    }
[2]315
[1286]316    return ;
317}
[2]318
319#endif /* CHECK_CUT_COUNTS */
320
[931]321
[640]322#ifdef CHECK_CUT_SIZE
323
324/*
325  Routine to verify that cut reference counts are correct.
326*/
327void verifyCutSize (const CbcTree * branchingTree, CbcModel &model)
[1286]328{
[640]329
[1286]330    int j ;
331    int nNodes = branchingTree->size() ;
332    int totalCuts = 0;
[640]333
[1286]334    /*
335      cut.tempNumber_ exists for the purpose of doing this verification. Clear it
336      in all cuts. We traverse the tree by starting from each live node and working
337      back to the root. At each CbcNodeInfo, check for cuts.
338    */
339    for (j = 0 ; j < nNodes ; j++) {
340        CbcNode *node = branchingTree->nodePointer(j) ;
341        CbcNodeInfo * nodeInfo = node->nodeInfo() ;
342        assert (node->nodeInfo()->numberBranchesLeft()) ;
343        while (nodeInfo) {
344            totalCuts += nodeInfo->numberCuts();
345            nodeInfo = nodeInfo->parent() ;
346        }
[640]347    }
[1286]348    printf("*** CHECKING cuts (size) after %d nodes - %d cuts\n", model.getNodeCount(), totalCuts) ;
349    return ;
[2]350}
351
[640]352#endif /* CHECK_CUT_SIZE */
353
354}
355
[1286]356/* End unnamed namespace for CbcModel.cpp */
[2]357
[1845]358
[1286]359void
[2]360CbcModel::analyzeObjective ()
361/*
362  Try to find a minimum change in the objective function. The first scan
363  checks that there are no continuous variables with non-zero coefficients,
364  and grabs the largest objective coefficient associated with an unfixed
365  integer variable. The second scan attempts to scale up the objective
366  coefficients to a point where they are sufficiently close to integer that
367  we can pretend they are integer, and calculate a gcd over the coefficients
368  of interest. This will be the minimum increment for the scaled coefficients.
369  The final action is to scale the increment back for the original coefficients
370  and install it, if it's better than the existing value.
371
372  John's note: We could do better than this.
373
374  John's second note - apologies for changing s to z
375*/
[1286]376{
377    const double *objective = getObjCoefficients() ;
378    const double *lower = getColLower() ;
379    const double *upper = getColUpper() ;
380    /*
381      Scan continuous and integer variables to see if continuous
382      are cover or network with integral rhs.
383    */
384    double continuousMultiplier = 1.0;
385    double * coeffMultiplier = NULL;
386    double largestObj = 0.0;
387    double smallestObj = COIN_DBL_MAX;
388    {
389        const double *rowLower = getRowLower() ;
390        const double *rowUpper = getRowUpper() ;
391        int numberRows = solver_->getNumRows() ;
392        double * rhs = new double [numberRows];
393        memset(rhs, 0, numberRows*sizeof(double));
394        int iColumn;
395        int numberColumns = solver_->getNumCols() ;
396        // Column copy of matrix
397        int problemType = -1;
398        const double * element = solver_->getMatrixByCol()->getElements();
399        const int * row = solver_->getMatrixByCol()->getIndices();
400        const CoinBigIndex * columnStart = solver_->getMatrixByCol()->getVectorStarts();
401        const int * columnLength = solver_->getMatrixByCol()->getVectorLengths();
402        int numberInteger = 0;
403        int numberIntegerObj = 0;
404        int numberGeneralIntegerObj = 0;
405        int numberIntegerWeight = 0;
406        int numberContinuousObj = 0;
407        double cost = COIN_DBL_MAX;
408        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
409            if (upper[iColumn] == lower[iColumn]) {
410                CoinBigIndex start = columnStart[iColumn];
411                CoinBigIndex end = start + columnLength[iColumn];
412                for (CoinBigIndex j = start; j < end; j++) {
413                    int iRow = row[j];
414                    rhs[iRow] += lower[iColumn] * element[j];
415                }
416            } else {
417                double objValue = objective[iColumn];
418                if (solver_->isInteger(iColumn))
419                    numberInteger++;
420                if (objValue) {
421                    if (!solver_->isInteger(iColumn)) {
422                        numberContinuousObj++;
423                    } else {
424                        largestObj = CoinMax(largestObj, fabs(objValue));
425                        smallestObj = CoinMin(smallestObj, fabs(objValue));
426                        numberIntegerObj++;
427                        if (cost == COIN_DBL_MAX)
428                            cost = objValue;
429                        else if (cost != objValue)
430                            cost = -COIN_DBL_MAX;
431                        int gap = static_cast<int> (upper[iColumn] - lower[iColumn]);
432                        if (gap > 1) {
433                            numberGeneralIntegerObj++;
434                            numberIntegerWeight += gap;
435                        }
436                    }
437                }
438            }
439        }
440        int iType = 0;
441        if (!numberContinuousObj && numberIntegerObj <= 5 && numberIntegerWeight <= 100 &&
442                numberIntegerObj*3 < numberObjects_ && !parentModel_ && solver_->getNumRows() > 100)
[1883]443          iType = 1 + 4 + ((moreSpecialOptions_&536870912)==0) ? 2 : 0;
[1286]444        else if (!numberContinuousObj && numberIntegerObj <= 100 &&
445                 numberIntegerObj*5 < numberObjects_ && numberIntegerWeight <= 100 &&
446                 !parentModel_ &&
447                 solver_->getNumRows() > 100 && cost != -COIN_DBL_MAX)
[1883]448          iType = 4 + ((moreSpecialOptions_&536870912)==0) ? 2 : 0;
[1286]449        else if (!numberContinuousObj && numberIntegerObj <= 100 &&
450                 numberIntegerObj*5 < numberObjects_ &&
451                 !parentModel_ &&
452                 solver_->getNumRows() > 100 && cost != -COIN_DBL_MAX)
453            iType = 8;
454        int iTest = getMaximumNodes();
455        if (iTest >= 987654320 && iTest < 987654330 && numberObjects_ && !parentModel_) {
456            iType = iTest - 987654320;
457            printf("Testing %d integer variables out of %d objects (%d integer) have cost of %g - %d continuous\n",
458                   numberIntegerObj, numberObjects_, numberInteger, cost, numberContinuousObj);
459            if (iType == 9)
460                exit(77);
461            if (numberContinuousObj)
462                iType = 0;
463        }
464
465        //if (!numberContinuousObj&&(numberIntegerObj<=5||cost!=-COIN_DBL_MAX)&&
466        //numberIntegerObj*3<numberObjects_&&!parentModel_&&solver_->getNumRows()>100) {
467        if (iType) {
468            /*
469            A) put high priority on (if none)
470            B) create artificial objective (if clp)
471            */
472            int iPriority = -1;
473            for (int i = 0; i < numberObjects_; i++) {
474                int k = object_[i]->priority();
475                if (iPriority == -1)
476                    iPriority = k;
477                else if (iPriority != k)
478                    iPriority = -2;
479            }
480            bool branchOnSatisfied = ((iType & 1) != 0);
481            bool createFake = ((iType & 2) != 0);
482            bool randomCost = ((iType & 4) != 0);
483            if (iPriority >= 0) {
484                char general[200];
485                if (cost == -COIN_DBL_MAX) {
486                    sprintf(general, "%d integer variables out of %d objects (%d integer) have costs - high priority",
487                            numberIntegerObj, numberObjects_, numberInteger);
488                } else if (cost == COIN_DBL_MAX) {
489                    sprintf(general, "No integer variables out of %d objects (%d integer) have costs",
490                            numberObjects_, numberInteger);
491                    branchOnSatisfied = false;
492                } else {
493                    sprintf(general, "%d integer variables out of %d objects (%d integer) have cost of %g - high priority",
494                            numberIntegerObj, numberObjects_, numberInteger, cost);
495                }
496                messageHandler()->message(CBC_GENERAL,
497                                          messages())
498                << general << CoinMessageEol ;
499                sprintf(general, "branch on satisfied %c create fake objective %c random cost %c",
500                        branchOnSatisfied ? 'Y' : 'N',
501                        createFake ? 'Y' : 'N',
502                        randomCost ? 'Y' : 'N');
503                messageHandler()->message(CBC_GENERAL,
504                                          messages())
505                << general << CoinMessageEol ;
506                // switch off clp type branching
[1883]507                // no ? fastNodeDepth_ = -1;
[1286]508                int highPriority = (branchOnSatisfied) ? -999 : 100;
509                for (int i = 0; i < numberObjects_; i++) {
510                    CbcSimpleInteger * thisOne = dynamic_cast <CbcSimpleInteger *> (object_[i]);
511                    object_[i]->setPriority(1000);
512                    if (thisOne) {
513                        int iColumn = thisOne->columnNumber();
514                        if (objective[iColumn])
515                            thisOne->setPriority(highPriority);
516                    }
517                }
518            }
[1088]519#ifdef COIN_HAS_CLP
[1286]520            OsiClpSolverInterface * clpSolver
521            = dynamic_cast<OsiClpSolverInterface *> (solver_);
522            if (clpSolver && createFake) {
523                // Create artificial objective to be used when all else fixed
524                int numberColumns = clpSolver->getNumCols();
525                double * fakeObj = new double [numberColumns];
526                // Column copy
527                const CoinPackedMatrix  * matrixByCol = clpSolver->getMatrixByCol();
528                //const double * element = matrixByCol.getElements();
529                //const int * row = matrixByCol.getIndices();
530                //const CoinBigIndex * columnStart = matrixByCol.getVectorStarts();
531                const int * columnLength = matrixByCol->getVectorLengths();
532                const double * solution = clpSolver->getColSolution();
[1393]533#ifdef JJF_ZERO
[1286]534                int nAtBound = 0;
535                for (int i = 0; i < numberColumns; i++) {
536                    double lowerValue = lower[i];
537                    double upperValue = upper[i];
538                    if (clpSolver->isInteger(i)) {
539                        double lowerValue = lower[i];
540                        double upperValue = upper[i];
541                        double value = solution[i];
542                        if (value < lowerValue + 1.0e-6 ||
543                                value > upperValue - 1.0e-6)
544                            nAtBound++;
545                    }
546                }
[1271]547#endif
[1409]548                /*
549                  Generate a random objective function for problems where the given objective
550                  function is not terribly useful. (Nearly feasible, single integer variable,
551                  that sort of thing.
552                */
[1286]553                CoinDrand48(true, 1234567);
554                for (int i = 0; i < numberColumns; i++) {
555                    double lowerValue = lower[i];
556                    double upperValue = upper[i];
557                    double value = (randomCost) ? ceil((CoinDrand48() + 0.5) * 1000)
558                                   : i + 1 + columnLength[i] * 1000;
559                    value *= 0.001;
560                    //value += columnLength[i];
561                    if (lowerValue > -1.0e5 || upperValue < 1.0e5) {
562                        if (fabs(lowerValue) > fabs(upperValue))
563                            value = - value;
564                        if (clpSolver->isInteger(i)) {
565                            double solValue = solution[i];
566                            // Better to add in 0.5 or 1.0??
567                            if (solValue < lowerValue + 1.0e-6)
568                                value = fabs(value) + 0.5; //fabs(value*1.5);
569                            else if (solValue > upperValue - 1.0e-6)
570                                value = -fabs(value) - 0.5; //-fabs(value*1.5);
571                        }
572                    } else {
573                        value = 0.0;
574                    }
575                    fakeObj[i] = value;
576                }
577                // pass to solver
578                clpSolver->setFakeObjective(fakeObj);
579                delete [] fakeObj;
580            }
[1088]581#endif
[1286]582        } else if (largestObj < smallestObj*5.0 && !parentModel_ &&
583                   !numberContinuousObj &&
584                   !numberGeneralIntegerObj &&
585                   numberIntegerObj*2 < numberColumns) {
586            // up priorities on costed
587            int iPriority = -1;
588            for (int i = 0; i < numberObjects_; i++) {
589                int k = object_[i]->priority();
590                if (iPriority == -1)
591                    iPriority = k;
592                else if (iPriority != k)
593                    iPriority = -2;
594            }
595            if (iPriority >= 100) {
[1271]596#ifdef CLP_INVESTIGATE
[1286]597                printf("Setting variables with obj to high priority\n");
[1271]598#endif
[1286]599                for (int i = 0; i < numberObjects_; i++) {
600                    CbcSimpleInteger * obj =
601                        dynamic_cast <CbcSimpleInteger *>(object_[i]) ;
602                    if (obj) {
603                        int iColumn = obj->columnNumber();
604                        if (objective[iColumn])
605                            object_[i]->setPriority(iPriority - 1);
606                    }
607                }
608            }
609        }
610        int iRow;
611        for (iRow = 0; iRow < numberRows; iRow++) {
612            if (rowLower[iRow] > -1.0e20 &&
613                    fabs(rowLower[iRow] - rhs[iRow] - floor(rowLower[iRow] - rhs[iRow] + 0.5)) > 1.0e-10) {
614                continuousMultiplier = 0.0;
615                break;
616            }
617            if (rowUpper[iRow] < 1.0e20 &&
618                    fabs(rowUpper[iRow] - rhs[iRow] - floor(rowUpper[iRow] - rhs[iRow] + 0.5)) > 1.0e-10) {
619                continuousMultiplier = 0.0;
620                break;
621            }
622            // set rhs to limiting value
623            if (rowLower[iRow] != rowUpper[iRow]) {
624                if (rowLower[iRow] > -1.0e20) {
625                    if (rowUpper[iRow] < 1.0e20) {
626                        // no good
627                        continuousMultiplier = 0.0;
628                        break;
629                    } else {
630                        rhs[iRow] = rowLower[iRow] - rhs[iRow];
631                        if (problemType < 0)
632                            problemType = 3; // set cover
633                        else if (problemType != 3)
634                            problemType = 4;
635                    }
636                } else {
637                    rhs[iRow] = rowUpper[iRow] - rhs[iRow];
638                    if (problemType < 0)
639                        problemType = 1; // set partitioning <=
640                    else if (problemType != 1)
641                        problemType = 4;
642                }
643            } else {
644                rhs[iRow] = rowUpper[iRow] - rhs[iRow];
645                if (problemType < 0)
646                    problemType = 3; // set partitioning ==
647                else if (problemType != 2)
648                    problemType = 2;
649            }
650            if (fabs(rhs[iRow] - 1.0) > 1.0e-12)
651                problemType = 4;
652        }
653        if (continuousMultiplier) {
654            // 1 network, 2 cover, 4 negative cover
655            int possible = 7;
656            bool unitRhs = true;
657            // See which rows could be set cover
658            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
659                if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
660                    CoinBigIndex start = columnStart[iColumn];
661                    CoinBigIndex end = start + columnLength[iColumn];
662                    for (CoinBigIndex j = start; j < end; j++) {
663                        double value = element[j];
664                        if (value == 1.0) {
665                        } else if (value == -1.0) {
666                            rhs[row[j]] = -0.5;
667                        } else {
668                            rhs[row[j]] = -COIN_DBL_MAX;
669                        }
670                    }
671                }
672            }
673            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
674                if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
675                    if (!isInteger(iColumn)) {
676                        CoinBigIndex start = columnStart[iColumn];
677                        CoinBigIndex end = start + columnLength[iColumn];
678                        double rhsValue = 0.0;
679                        // 1 all ones, -1 all -1s, 2 all +- 1, 3 no good
680                        int type = 0;
681                        for (CoinBigIndex j = start; j < end; j++) {
682                            double value = element[j];
683                            if (fabs(value) != 1.0) {
684                                type = 3;
685                                break;
686                            } else if (value == 1.0) {
687                                if (!type)
688                                    type = 1;
689                                else if (type != 1)
690                                    type = 2;
691                            } else {
692                                if (!type)
693                                    type = -1;
694                                else if (type != -1)
695                                    type = 2;
696                            }
697                            int iRow = row[j];
698                            if (rhs[iRow] == -COIN_DBL_MAX) {
699                                type = 3;
700                                break;
701                            } else if (rhs[iRow] == -0.5) {
702                                // different values
703                                unitRhs = false;
704                            } else if (rhsValue) {
705                                if (rhsValue != rhs[iRow])
706                                    unitRhs = false;
707                            } else {
708                                rhsValue = rhs[iRow];
709                            }
710                        }
711                        // if no elements OK
712                        if (type == 3) {
713                            // no good
714                            possible = 0;
715                            break;
716                        } else if (type == 2) {
717                            if (end - start > 2) {
718                                // no good
719                                possible = 0;
720                                break;
721                            } else {
722                                // only network
723                                possible &= 1;
724                                if (!possible)
725                                    break;
726                            }
727                        } else if (type == 1) {
728                            // only cover
729                            possible &= 2;
730                            if (!possible)
731                                break;
732                        } else if (type == -1) {
733                            // only negative cover
734                            possible &= 4;
735                            if (!possible)
736                                break;
737                        }
738                    }
739                }
740            }
741            if ((possible == 2 || possible == 4) && !unitRhs) {
[931]742#if COIN_DEVELOP>1
[1286]743                printf("XXXXXX Continuous all +1 but different rhs\n");
[810]744#endif
[1286]745                possible = 0;
746            }
747            // may be all integer
748            if (possible != 7) {
749                if (!possible)
750                    continuousMultiplier = 0.0;
751                else if (possible == 1)
752                    continuousMultiplier = 1.0;
753                else
754                    continuousMultiplier = 0.0; // 0.5 was incorrect;
[931]755#if COIN_DEVELOP>1
[1286]756                if (continuousMultiplier)
757                    printf("XXXXXX multiplier of %g\n", continuousMultiplier);
[810]758#endif
[1286]759                if (continuousMultiplier == 0.5) {
760                    coeffMultiplier = new double [numberColumns];
761                    bool allOne = true;
762                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
763                        coeffMultiplier[iColumn] = 1.0;
764                        if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
765                            if (!isInteger(iColumn)) {
766                                CoinBigIndex start = columnStart[iColumn];
767                                int iRow = row[start];
768                                double value = rhs[iRow];
769                                assert (value >= 0.0);
770                                if (value != 0.0 && value != 1.0)
771                                    allOne = false;
772                                coeffMultiplier[iColumn] = 0.5 * value;
773                            }
774                        }
775                    }
776                    if (allOne) {
777                        // back to old way
778                        delete [] coeffMultiplier;
779                        coeffMultiplier = NULL;
780                    }
781                }
782            } else {
783                // all integer
784                problemType_ = problemType;
[931]785#if COIN_DEVELOP>1
[1286]786                printf("Problem type is %d\n", problemType_);
[814]787#endif
[1286]788            }
789        }
790
791        // But try again
792        if (continuousMultiplier < 1.0) {
793            memset(rhs, 0, numberRows*sizeof(double));
794            int * count = new int [numberRows];
795            memset(count, 0, numberRows*sizeof(int));
796            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
797                CoinBigIndex start = columnStart[iColumn];
798                CoinBigIndex end = start + columnLength[iColumn];
799                if (upper[iColumn] == lower[iColumn]) {
800                    for (CoinBigIndex j = start; j < end; j++) {
801                        int iRow = row[j];
802                        rhs[iRow] += lower[iColumn] * element[j];
803                    }
804                } else if (solver_->isInteger(iColumn)) {
805                    for (CoinBigIndex j = start; j < end; j++) {
806                        int iRow = row[j];
807                        if (fabs(element[j] - floor(element[j] + 0.5)) > 1.0e-10)
808                            rhs[iRow]  = COIN_DBL_MAX;
809                    }
810                } else {
811                    for (CoinBigIndex j = start; j < end; j++) {
812                        int iRow = row[j];
813                        count[iRow]++;
814                        if (fabs(element[j]) != 1.0)
815                            rhs[iRow]  = COIN_DBL_MAX;
816                    }
817                }
818            }
819            // now look at continuous
820            bool allGood = true;
821            double direction = solver_->getObjSense() ;
822            int numberObj = 0;
823            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
824                if (upper[iColumn] > lower[iColumn]) {
825                    double objValue = objective[iColumn] * direction;
826                    if (objValue && !solver_->isInteger(iColumn)) {
827                        numberObj++;
828                        CoinBigIndex start = columnStart[iColumn];
829                        CoinBigIndex end = start + columnLength[iColumn];
830                        if (objValue > 0.0) {
831                            // wants to be as low as possible
832                            if (lower[iColumn] < -1.0e10 || fabs(lower[iColumn] - floor(lower[iColumn] + 0.5)) > 1.0e-10) {
833                                allGood = false;
834                                break;
835                            } else if (upper[iColumn] < 1.0e10 && fabs(upper[iColumn] - floor(upper[iColumn] + 0.5)) > 1.0e-10) {
836                                allGood = false;
837                                break;
838                            }
839                            bool singletonRow = true;
840                            bool equality = false;
841                            for (CoinBigIndex j = start; j < end; j++) {
842                                int iRow = row[j];
843                                if (count[iRow] > 1)
844                                    singletonRow = false;
845                                else if (rowLower[iRow] == rowUpper[iRow])
846                                    equality = true;
847                                double rhsValue = rhs[iRow];
848                                double lowerValue = rowLower[iRow];
849                                double upperValue = rowUpper[iRow];
850                                if (rhsValue < 1.0e20) {
851                                    if (lowerValue > -1.0e20)
852                                        lowerValue -= rhsValue;
853                                    if (upperValue < 1.0e20)
854                                        upperValue -= rhsValue;
855                                }
856                                if (fabs(rhsValue) > 1.0e20 || fabs(rhsValue - floor(rhsValue + 0.5)) > 1.0e-10
857                                        || fabs(element[j]) != 1.0) {
858                                    // no good
859                                    allGood = false;
860                                    break;
861                                }
862                                if (element[j] > 0.0) {
863                                    if (lowerValue > -1.0e20 && fabs(lowerValue - floor(lowerValue + 0.5)) > 1.0e-10) {
864                                        // no good
865                                        allGood = false;
866                                        break;
867                                    }
868                                } else {
869                                    if (upperValue < 1.0e20 && fabs(upperValue - floor(upperValue + 0.5)) > 1.0e-10) {
870                                        // no good
871                                        allGood = false;
872                                        break;
873                                    }
874                                }
875                            }
876                            if (!singletonRow && end > start + 1 && !equality)
877                                allGood = false;
878                            if (!allGood)
879                                break;
880                        } else {
881                            // wants to be as high as possible
882                            if (upper[iColumn] > 1.0e10 || fabs(upper[iColumn] - floor(upper[iColumn] + 0.5)) > 1.0e-10) {
883                                allGood = false;
884                                break;
885                            } else if (lower[iColumn] > -1.0e10 && fabs(lower[iColumn] - floor(lower[iColumn] + 0.5)) > 1.0e-10) {
886                                allGood = false;
887                                break;
888                            }
889                            bool singletonRow = true;
890                            bool equality = false;
891                            for (CoinBigIndex j = start; j < end; j++) {
892                                int iRow = row[j];
893                                if (count[iRow] > 1)
894                                    singletonRow = false;
895                                else if (rowLower[iRow] == rowUpper[iRow])
896                                    equality = true;
897                                double rhsValue = rhs[iRow];
898                                double lowerValue = rowLower[iRow];
899                                double upperValue = rowUpper[iRow];
900                                if (rhsValue < 1.0e20) {
901                                    if (lowerValue > -1.0e20)
902                                        lowerValue -= rhsValue;
903                                    if (upperValue < 1.0e20)
904                                        upperValue -= rhsValue;
905                                }
906                                if (fabs(rhsValue) > 1.0e20 || fabs(rhsValue - floor(rhsValue + 0.5)) > 1.0e-10
907                                        || fabs(element[j]) != 1.0) {
908                                    // no good
909                                    allGood = false;
910                                    break;
911                                }
912                                if (element[j] < 0.0) {
913                                    if (lowerValue > -1.0e20 && fabs(lowerValue - floor(lowerValue + 0.5)) > 1.0e-10) {
914                                        // no good
915                                        allGood = false;
916                                        break;
917                                    }
918                                } else {
919                                    if (upperValue < 1.0e20 && fabs(upperValue - floor(upperValue + 0.5)) > 1.0e-10) {
920                                        // no good
921                                        allGood = false;
922                                        break;
923                                    }
924                                }
925                            }
926                            if (!singletonRow && end > start + 1 && !equality)
927                                allGood = false;
928                            if (!allGood)
929                                break;
930                        }
931                    }
932                }
933            }
934            delete [] count;
935            if (allGood) {
[931]936#if COIN_DEVELOP>1
[1286]937                if (numberObj)
938                    printf("YYYY analysis says all continuous with costs will be integer\n");
[837]939#endif
[1286]940                continuousMultiplier = 1.0;
941            }
942        }
943        delete [] rhs;
[833]944    }
[1286]945    /*
946      Take a first scan to see if there are unfixed continuous variables in the
947      objective.  If so, the minimum objective change could be arbitrarily small.
948      Also pick off the maximum coefficient of an unfixed integer variable.
[2]949
[1286]950      If the objective is found to contain only integer variables, set the
951      fathoming discipline to strict.
952    */
953    double maximumCost = 0.0 ;
954    //double trueIncrement=0.0;
955    int iColumn ;
956    int numberColumns = getNumCols() ;
957    double scaleFactor = 1.0; // due to rhs etc
[1409]958    /*
959      Original model did not have integer bounds.
960    */
[1286]961    if ((specialOptions_&65536) == 0) {
962        /* be on safe side (later look carefully as may be able to
963           to get 0.5 say if bounds are multiples of 0.5 */
964        for (iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
965            if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
966                double value;
967                value = fabs(lower[iColumn]);
968                if (floor(value + 0.5) != value) {
969                    scaleFactor = CoinMin(scaleFactor, 0.5);
970                    if (floor(2.0*value + 0.5) != 2.0*value) {
971                        scaleFactor = CoinMin(scaleFactor, 0.25);
972                        if (floor(4.0*value + 0.5) != 4.0*value) {
973                            scaleFactor = 0.0;
974                        }
975                    }
976                }
977                value = fabs(upper[iColumn]);
978                if (floor(value + 0.5) != value) {
979                    scaleFactor = CoinMin(scaleFactor, 0.5);
980                    if (floor(2.0*value + 0.5) != 2.0*value) {
981                        scaleFactor = CoinMin(scaleFactor, 0.25);
982                        if (floor(4.0*value + 0.5) != 4.0*value) {
983                            scaleFactor = 0.0;
984                        }
985                    }
986                }
987            }
988        }
989    }
990    bool possibleMultiple = continuousMultiplier != 0.0 && scaleFactor != 0.0 ;
991    if (possibleMultiple) {
992        for (iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
993            if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
994                maximumCost = CoinMax(maximumCost, fabs(objective[iColumn])) ;
995            }
996        }
997    }
998    setIntParam(CbcModel::CbcFathomDiscipline, possibleMultiple) ;
999    /*
1000      If a nontrivial increment is possible, try and figure it out. We're looking
1001      for gcd(c<j>) for all c<j> that are coefficients of unfixed integer
1002      variables. Since the c<j> might not be integers, try and inflate them
1003      sufficiently that they look like integers (and we'll deflate the gcd
1004      later).
[2]1005
[1286]1006      2520.0 is used as it is a nice multiple of 2,3,5,7
1007    */
1008    if (possibleMultiple && maximumCost) {
1009        int increment = 0 ;
1010        double multiplier = 2520.0 ;
1011        while (10.0*multiplier*maximumCost < 1.0e8)
1012            multiplier *= 10.0 ;
1013        int bigIntegers = 0; // Count of large costs which are integer
1014        for (iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
1015            if (upper[iColumn] > lower[iColumn] + 1.0e-8) {
1016                double objValue = fabs(objective[iColumn]);
1017                if (!isInteger(iColumn)) {
1018                    if (!coeffMultiplier)
1019                        objValue *= continuousMultiplier;
1020                    else
1021                        objValue *= coeffMultiplier[iColumn];
1022                }
1023                if (objValue) {
1024                    double value = objValue * multiplier ;
1025                    if (value < 2.1e9) {
1026                        int nearest = static_cast<int> (floor(value + 0.5)) ;
1027                        if (fabs(value - floor(value + 0.5)) > 1.0e-8) {
1028                            increment = 0 ;
1029                            break ;
1030                        } else if (!increment) {
1031                            increment = nearest ;
1032                        } else {
1033                            increment = gcd(increment, nearest) ;
1034                        }
1035                    } else {
1036                        // large value - may still be multiple of 1.0
1037                        if (fabs(objValue - floor(objValue + 0.5)) > 1.0e-8) {
1038                            increment = 0;
1039                            break;
1040                        } else {
1041                            bigIntegers++;
1042                        }
1043                    }
1044                }
1045            }
1046        }
1047        delete [] coeffMultiplier;
1048        /*
1049          If the increment beats the current value for objective change, install it.
1050        */
1051        if (increment) {
1052            double value = increment ;
1053            double cutoff = getDblParam(CbcModel::CbcCutoffIncrement) ;
1054            if (bigIntegers) {
1055                // allow for 1.0
1056                increment = gcd(increment, static_cast<int> (multiplier));
1057                value = increment;
1058            }
1059            value /= multiplier ;
1060            value *= scaleFactor;
1061            //trueIncrement=CoinMax(cutoff,value);;
1062            if (value*0.999 > cutoff) {
[1910]1063                if (messageHandler()->logLevel() > 2){
1064                    messageHandler()->message(CBC_INTEGERINCREMENT,
1065                                              messages())
1066                       << value << CoinMessageEol ;
1067                }
1068                setDblParam(CbcModel::CbcCutoffIncrement, CoinMax(value*0.999,value-1.0e-4)) ;
[1286]1069            }
1070        }
[474]1071    }
[2]1072
[1286]1073    return ;
[2]1074}
1075
[1332]1076/*
1077saveModel called (carved out of) BranchandBound
1078*/
[1330]1079void CbcModel::saveModel(OsiSolverInterface * saveSolver, double * checkCutoffForRestart, bool * feasible)
1080{
1081    if (saveSolver && (specialOptions_&32768) != 0) {
1082        // See if worth trying reduction
1083        *checkCutoffForRestart = getCutoff();
1084        bool tryNewSearch = solverCharacteristics_->reducedCostsAccurate() &&
1085                            (*checkCutoffForRestart < 1.0e20);
1086        int numberColumns = getNumCols();
1087        if (tryNewSearch) {
1088#ifdef CLP_INVESTIGATE
1089            printf("after %d nodes, cutoff %g - looking\n",
1090                   numberNodes_, getCutoff());
1091#endif
1092            saveSolver->resolve();
1093            double direction = saveSolver->getObjSense() ;
1094            double gap = *checkCutoffForRestart - saveSolver->getObjValue() * direction ;
1095            double tolerance;
1096            saveSolver->getDblParam(OsiDualTolerance, tolerance) ;
1097            if (gap <= 0.0)
1098                gap = tolerance;
1099            gap += 100.0 * tolerance;
1100            double integerTolerance = getDblParam(CbcIntegerTolerance) ;
1101
1102            const double *lower = saveSolver->getColLower() ;
1103            const double *upper = saveSolver->getColUpper() ;
1104            const double *solution = saveSolver->getColSolution() ;
1105            const double *reducedCost = saveSolver->getReducedCost() ;
1106
1107            int numberFixed = 0 ;
1108            int numberFixed2 = 0;
1109            for (int i = 0 ; i < numberIntegers_ ; i++) {
1110                int iColumn = integerVariable_[i] ;
1111                double djValue = direction * reducedCost[iColumn] ;
1112                if (upper[iColumn] - lower[iColumn] > integerTolerance) {
1113                    if (solution[iColumn] < lower[iColumn] + integerTolerance && djValue > gap) {
1114                        saveSolver->setColUpper(iColumn, lower[iColumn]) ;
1115                        numberFixed++ ;
1116                    } else if (solution[iColumn] > upper[iColumn] - integerTolerance && -djValue > gap) {
1117                        saveSolver->setColLower(iColumn, upper[iColumn]) ;
1118                        numberFixed++ ;
1119                    }
1120                } else {
1121                    numberFixed2++;
1122                }
1123            }
1124#ifdef COIN_DEVELOP
[1409]1125            /*
1126              We're debugging. (specialOptions 1)
1127            */
[1330]1128            if ((specialOptions_&1) != 0) {
1129                const OsiRowCutDebugger *debugger = saveSolver->getRowCutDebugger() ;
1130                if (debugger) {
1131                    printf("Contains optimal\n") ;
1132                    OsiSolverInterface * temp = saveSolver->clone();
1133                    const double * solution = debugger->optimalSolution();
1134                    const double *lower = temp->getColLower() ;
1135                    const double *upper = temp->getColUpper() ;
1136                    int n = temp->getNumCols();
1137                    for (int i = 0; i < n; i++) {
1138                        if (temp->isInteger(i)) {
1139                            double value = floor(solution[i] + 0.5);
1140                            assert (value >= lower[i] && value <= upper[i]);
1141                            temp->setColLower(i, value);
1142                            temp->setColUpper(i, value);
1143                        }
1144                    }
1145                    temp->writeMps("reduced_fix");
1146                    delete temp;
1147                    saveSolver->writeMps("reduced");
1148                } else {
1149                    abort();
1150                }
1151            }
1152            printf("Restart could fix %d integers (%d already fixed)\n",
1153                   numberFixed + numberFixed2, numberFixed2);
1154#endif
1155            numberFixed += numberFixed2;
1156            if (numberFixed*20 < numberColumns)
1157                tryNewSearch = false;
1158        }
1159        if (tryNewSearch) {
1160            // back to solver without cuts?
1161            OsiSolverInterface * solver2 = continuousSolver_->clone();
1162            const double *lower = saveSolver->getColLower() ;
1163            const double *upper = saveSolver->getColUpper() ;
1164            for (int i = 0 ; i < numberIntegers_ ; i++) {
1165                int iColumn = integerVariable_[i] ;
1166                solver2->setColLower(iColumn, lower[iColumn]);
1167                solver2->setColUpper(iColumn, upper[iColumn]);
1168            }
1169            // swap
1170            delete saveSolver;
1171            saveSolver = solver2;
1172            double * newSolution = new double[numberColumns];
1173            double objectiveValue = *checkCutoffForRestart;
1174            CbcSerendipity heuristic(*this);
1175            if (bestSolution_)
1176                heuristic.setInputSolution(bestSolution_, bestObjective_);
1177            heuristic.setFractionSmall(0.9);
1178            heuristic.setFeasibilityPumpOptions(1008013);
1179            // Use numberNodes to say how many are original rows
1180            heuristic.setNumberNodes(continuousSolver_->getNumRows());
1181#ifdef COIN_DEVELOP
1182            if (continuousSolver_->getNumRows() <
1183                    saveSolver->getNumRows())
1184                printf("%d rows added ZZZZZ\n",
1185                       solver_->getNumRows() - continuousSolver_->getNumRows());
1186#endif
1187            int returnCode = heuristic.smallBranchAndBound(saveSolver,
1188                             -1, newSolution,
1189                             objectiveValue,
1190                             *checkCutoffForRestart, "Reduce");
1191            if (returnCode < 0) {
1192#ifdef COIN_DEVELOP
1193                printf("Restart - not small enough to do search after fixing\n");
1194#endif
1195                delete [] newSolution;
1196            } else {
1197                if ((returnCode&1) != 0) {
1198                    // increment number of solutions so other heuristics can test
1199                    numberSolutions_++;
1200                    numberHeuristicSolutions_++;
1201                    lastHeuristic_ = NULL;
1202                    setBestSolution(CBC_ROUNDING, objectiveValue, newSolution) ;
1203                }
1204                delete [] newSolution;
1205                *feasible = false; // stop search
1206            }
[1412]1207#if 0 // probably not needed def CBC_THREAD
1208            if (master_) {
1209                lockThread();
1210                if (parallelMode() > 0) {
1211                    while (master_->waitForThreadsInTree(0)) {
1212                        lockThread();
1213                        double dummyBest;
1214                        tree_->cleanTree(this, -COIN_DBL_MAX, dummyBest) ;
1215                        //unlockThread();
1216                    }
1217                }
1218                master_->waitForThreadsInTree(2);
1219                delete master_;
1220                master_ = NULL;
1221                masterThread_ = NULL;
1222            }
1223#endif
[1330]1224        }
1225    }
1226}
1227/*
1228Adds integers, called from BranchandBound()
1229*/
1230void CbcModel::AddIntegers()
1231{
1232    int numberColumns = continuousSolver_->getNumCols();
1233    int numberRows = continuousSolver_->getNumRows();
1234    int * del = new int [CoinMax(numberColumns, numberRows)];
1235    int * original = new int [numberColumns];
1236    char * possibleRow = new char [numberRows];
1237    {
1238        const CoinPackedMatrix * rowCopy = continuousSolver_->getMatrixByRow();
1239        const int * column = rowCopy->getIndices();
1240        const int * rowLength = rowCopy->getVectorLengths();
1241        const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
1242        const double * rowLower = continuousSolver_->getRowLower();
1243        const double * rowUpper = continuousSolver_->getRowUpper();
1244        const double * element = rowCopy->getElements();
1245        for (int i = 0; i < numberRows; i++) {
1246            int nLeft = 0;
1247            bool possible = false;
1248            if (rowLower[i] < -1.0e20) {
1249                double value = rowUpper[i];
1250                if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1251                    possible = true;
1252            } else if (rowUpper[i] > 1.0e20) {
1253                double value = rowLower[i];
1254                if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1255                    possible = true;
1256            } else {
1257                double value = rowUpper[i];
1258                if (rowLower[i] == rowUpper[i] &&
1259                        fabs(value - floor(value + 0.5)) < 1.0e-8)
1260                    possible = true;
1261            }
[1639]1262            double allSame = (possible) ? 0.0 : -1.0;
[1330]1263            for (CoinBigIndex j = rowStart[i];
1264                    j < rowStart[i] + rowLength[i]; j++) {
1265                int iColumn = column[j];
1266                if (continuousSolver_->isInteger(iColumn)) {
1267                    if (fabs(element[j]) != 1.0)
1268                        possible = false;
1269                } else {
1270                    nLeft++;
[1639]1271                    if (!allSame) {
1272                      allSame = fabs(element[j]);
1273                    } else if (allSame>0.0) {
1274                      if (allSame!=fabs(element[j]))
1275                        allSame = -1.0;
1276                    }
[1330]1277                }
1278            }
[1639]1279            if (nLeft == rowLength[i] && allSame > 0.0)
1280                possibleRow[i] = 2;
1281            else if (possible || !nLeft)
[1330]1282                possibleRow[i] = 1;
1283            else
1284                possibleRow[i] = 0;
1285        }
1286    }
1287    int nDel = 0;
1288    for (int i = 0; i < numberColumns; i++) {
1289        original[i] = i;
1290        if (continuousSolver_->isInteger(i))
1291            del[nDel++] = i;
1292    }
1293    int nExtra = 0;
1294    OsiSolverInterface * copy1 = continuousSolver_->clone();
1295    int nPass = 0;
1296    while (nDel && nPass < 10) {
1297        nPass++;
1298        OsiSolverInterface * copy2 = copy1->clone();
1299        int nLeft = 0;
1300        for (int i = 0; i < nDel; i++)
1301            original[del[i]] = -1;
1302        for (int i = 0; i < numberColumns; i++) {
1303            int kOrig = original[i];
1304            if (kOrig >= 0)
1305                original[nLeft++] = kOrig;
1306        }
1307        assert (nLeft == numberColumns - nDel);
1308        copy2->deleteCols(nDel, del);
1309        numberColumns = copy2->getNumCols();
1310        const CoinPackedMatrix * rowCopy = copy2->getMatrixByRow();
1311        numberRows = rowCopy->getNumRows();
1312        const int * column = rowCopy->getIndices();
1313        const int * rowLength = rowCopy->getVectorLengths();
1314        const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
1315        const double * rowLower = copy2->getRowLower();
1316        const double * rowUpper = copy2->getRowUpper();
1317        const double * element = rowCopy->getElements();
1318        const CoinPackedMatrix * columnCopy = copy2->getMatrixByCol();
1319        const int * columnLength = columnCopy->getVectorLengths();
1320        nDel = 0;
1321        // Could do gcd stuff on ones with costs
1322        for (int i = 0; i < numberRows; i++) {
1323            if (!rowLength[i]) {
1324                del[nDel++] = i;
1325            } else if (possibleRow[i]) {
1326                if (rowLength[i] == 1) {
1327                    int k = rowStart[i];
1328                    int iColumn = column[k];
1329                    if (!copy2->isInteger(iColumn)) {
1330                        double mult = 1.0 / fabs(element[k]);
1331                        if (rowLower[i] < -1.0e20) {
[1639]1332                            // treat rhs as multiple of 1 unless elements all same
1333                            double value = ((possibleRow[i]==2) ? rowUpper[i] : 1.0) * mult;
[1330]1334                            if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
1335                                del[nDel++] = i;
1336                                if (columnLength[iColumn] == 1) {
1337                                    copy2->setInteger(iColumn);
1338                                    int kOrig = original[iColumn];
1339                                    setOptionalInteger(kOrig);
1340                                }
1341                            }
1342                        } else if (rowUpper[i] > 1.0e20) {
[1639]1343                            // treat rhs as multiple of 1 unless elements all same
1344                            double value = ((possibleRow[i]==2) ? rowLower[i] : 1.0) * mult;
[1330]1345                            if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
1346                                del[nDel++] = i;
1347                                if (columnLength[iColumn] == 1) {
1348                                    copy2->setInteger(iColumn);
1349                                    int kOrig = original[iColumn];
1350                                    setOptionalInteger(kOrig);
1351                                }
1352                            }
1353                        } else {
[1639]1354                            // treat rhs as multiple of 1 unless elements all same
1355                            double value = ((possibleRow[i]==2) ? rowUpper[i] : 1.0) * mult;
[1330]1356                            if (rowLower[i] == rowUpper[i] &&
1357                                    fabs(value - floor(value + 0.5)) < 1.0e-8) {
1358                                del[nDel++] = i;
1359                                copy2->setInteger(iColumn);
1360                                int kOrig = original[iColumn];
1361                                setOptionalInteger(kOrig);
1362                            }
1363                        }
1364                    }
1365                } else {
1366                    // only if all singletons
1367                    bool possible = false;
1368                    if (rowLower[i] < -1.0e20) {
1369                        double value = rowUpper[i];
1370                        if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1371                            possible = true;
1372                    } else if (rowUpper[i] > 1.0e20) {
1373                        double value = rowLower[i];
1374                        if (fabs(value - floor(value + 0.5)) < 1.0e-8)
1375                            possible = true;
1376                    } else {
1377                        double value = rowUpper[i];
1378                        if (rowLower[i] == rowUpper[i] &&
1379                                fabs(value - floor(value + 0.5)) < 1.0e-8)
1380                            possible = true;
1381                    }
1382                    if (possible) {
1383                        for (CoinBigIndex j = rowStart[i];
1384                                j < rowStart[i] + rowLength[i]; j++) {
1385                            int iColumn = column[j];
1386                            if (columnLength[iColumn] != 1 || fabs(element[j]) != 1.0) {
1387                                possible = false;
1388                                break;
1389                            }
1390                        }
1391                        if (possible) {
1392                            for (CoinBigIndex j = rowStart[i];
1393                                    j < rowStart[i] + rowLength[i]; j++) {
1394                                int iColumn = column[j];
1395                                if (!copy2->isInteger(iColumn)) {
1396                                    copy2->setInteger(iColumn);
1397                                    int kOrig = original[iColumn];
1398                                    setOptionalInteger(kOrig);
1399                                }
1400                            }
1401                            del[nDel++] = i;
1402                        }
1403                    }
1404                }
1405            }
1406        }
1407        if (nDel) {
1408            copy2->deleteRows(nDel, del);
[1773]1409            // pack down possible
1410            int n=0;
1411            for (int i=0;i<nDel;i++)
1412              possibleRow[del[i]]=-1;
1413            for (int i=0;i<numberRows;i++) {
1414              if (possibleRow[i]>=0) 
1415                possibleRow[n++]=possibleRow[i];
1416            }
[1330]1417        }
1418        if (nDel != numberRows) {
1419            nDel = 0;
1420            for (int i = 0; i < numberColumns; i++) {
1421                if (copy2->isInteger(i)) {
1422                    del[nDel++] = i;
1423                    nExtra++;
1424                }
1425            }
1426        } else {
1427            nDel = 0;
1428        }
1429        delete copy1;
1430        copy1 = copy2->clone();
1431        delete copy2;
1432    }
1433    // See if what's left is a network
1434    bool couldBeNetwork = false;
1435    if (copy1->getNumRows() && copy1->getNumCols()) {
1436#ifdef COIN_HAS_CLP
1437        OsiClpSolverInterface * clpSolver
1438        = dynamic_cast<OsiClpSolverInterface *> (copy1);
1439        if (false && clpSolver) {
1440            numberRows = clpSolver->getNumRows();
1441            char * rotate = new char[numberRows];
1442            int n = clpSolver->getModelPtr()->findNetwork(rotate, 1.0);
1443            delete [] rotate;
1444#ifdef CLP_INVESTIGATE
1445            printf("INTA network %d rows out of %d\n", n, numberRows);
1446#endif
1447            if (CoinAbs(n) == numberRows) {
1448                couldBeNetwork = true;
1449                for (int i = 0; i < numberRows; i++) {
1450                    if (!possibleRow[i]) {
1451                        couldBeNetwork = false;
1452#ifdef CLP_INVESTIGATE
1453                        printf("but row %d is bad\n", i);
1454#endif
1455                        break;
1456                    }
1457                }
1458            }
1459        } else
1460#endif
1461        {
1462            numberColumns = copy1->getNumCols();
1463            numberRows = copy1->getNumRows();
1464            const double * rowLower = copy1->getRowLower();
1465            const double * rowUpper = copy1->getRowUpper();
1466            couldBeNetwork = true;
1467            for (int i = 0; i < numberRows; i++) {
1468                if (rowLower[i] > -1.0e20 && fabs(rowLower[i] - floor(rowLower[i] + 0.5)) > 1.0e-12) {
1469                    couldBeNetwork = false;
1470                    break;
1471                }
1472                if (rowUpper[i] < 1.0e20 && fabs(rowUpper[i] - floor(rowUpper[i] + 0.5)) > 1.0e-12) {
1473                    couldBeNetwork = false;
1474                    break;
1475                }
[1773]1476                if (possibleRow[i]==0) {
1477                    couldBeNetwork = false;
1478                    break;
1479                }
[1330]1480            }
1481            if (couldBeNetwork) {
1482                const CoinPackedMatrix  * matrixByCol = copy1->getMatrixByCol();
1483                const double * element = matrixByCol->getElements();
1484                //const int * row = matrixByCol->getIndices();
1485                const CoinBigIndex * columnStart = matrixByCol->getVectorStarts();
1486                const int * columnLength = matrixByCol->getVectorLengths();
1487                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1488                    CoinBigIndex start = columnStart[iColumn];
1489                    CoinBigIndex end = start + columnLength[iColumn];
1490                    if (end > start + 2) {
1491                        couldBeNetwork = false;
1492                        break;
1493                    }
1494                    int type = 0;
1495                    for (CoinBigIndex j = start; j < end; j++) {
1496                        double value = element[j];
1497                        if (fabs(value) != 1.0) {
1498                            couldBeNetwork = false;
1499                            break;
1500                        } else if (value == 1.0) {
1501                            if ((type&1) == 0)
1502                                type |= 1;
1503                            else
1504                                type = 7;
1505                        } else if (value == -1.0) {
1506                            if ((type&2) == 0)
1507                                type |= 2;
1508                            else
1509                                type = 7;
1510                        }
1511                    }
1512                    if (type > 3) {
1513                        couldBeNetwork = false;
1514                        break;
1515                    }
1516                }
1517            }
1518        }
1519    }
1520    if (couldBeNetwork) {
1521        for (int i = 0; i < numberColumns; i++)
1522            setOptionalInteger(original[i]);
1523    }
1524    if (nExtra || couldBeNetwork) {
1525        numberColumns = copy1->getNumCols();
1526        numberRows = copy1->getNumRows();
1527        if (!numberColumns || !numberRows) {
1528            int numberColumns = solver_->getNumCols();
1529            for (int i = 0; i < numberColumns; i++)
1530                assert(solver_->isInteger(i));
1531        }
1532#ifdef CLP_INVESTIGATE
1533        if (couldBeNetwork || nExtra)
1534            printf("INTA %d extra integers, %d left%s\n", nExtra,
1535                   numberColumns,
1536                   couldBeNetwork ? ", all network" : "");
1537#endif
1538        findIntegers(true, 2);
1539        convertToDynamic();
1540    }
1541#ifdef CLP_INVESTIGATE
1542    if (!couldBeNetwork && copy1->getNumCols() &&
1543            copy1->getNumRows()) {
1544        printf("INTA %d rows and %d columns remain\n",
1545               copy1->getNumRows(), copy1->getNumCols());
1546        if (copy1->getNumCols() < 200) {
1547            copy1->writeMps("moreint");
1548            printf("INTA Written remainder to moreint.mps.gz %d rows %d cols\n",
1549                   copy1->getNumRows(), copy1->getNumCols());
1550        }
1551    }
1552#endif
1553    delete copy1;
1554    delete [] del;
1555    delete [] original;
1556    delete [] possibleRow;
1557    // double check increment
1558    analyzeObjective();
1559}
[2]1560/**
1561  \todo
1562  Normally, it looks like we enter here from command dispatch in the main
1563  routine, after calling the solver for an initial solution
1564  (CbcModel::initialSolve, which simply calls the solver's initialSolve
1565  routine.) The first thing we do is call resolve. Presumably there are
1566  circumstances where this is nontrivial? There's also a call from
1567  CbcModel::originalModel (tied up with integer presolve), which should be
1568  checked.
1569
1570*/
1571
1572/*
1573  The overall flow can be divided into three stages:
1574    * Prep: Check that the lp relaxation remains feasible at the root. If so,
1575      do all the setup for B&C.
1576    * Process the root node: Generate cuts, apply heuristics, and in general do
1577      the best we can to resolve the problem without B&C.
1578    * Do B&C search until we hit a limit or exhaust the search tree.
[1286]1579
[2]1580  Keep in mind that in general there is no node in the search tree that
1581  corresponds to the active subproblem. The active subproblem is represented
1582  by the current state of the model,  of the solver, and of the constraint
1583  system held by the solver.
1584*/
[1132]1585#ifdef COIN_HAS_CPX
1586#include "OsiCpxSolverInterface.hpp"
[1271]1587#include "cplex.h"
[1132]1588#endif
[1286]1589void CbcModel::branchAndBound(int doStatistics)
[2]1590
1591{
[1883]1592  if (!parentModel_)
[1286]1593    /*
[1623]1594      Capture a time stamp before we start (unless set).
[1286]1595    */
[1623]1596    if (!dblParam_[CbcStartSeconds]) {
1597      if (!useElapsedTime())
1598        dblParam_[CbcStartSeconds] = CoinCpuTime();
1599      else
1600        dblParam_[CbcStartSeconds] = CoinGetTimeOfDay();
1601    }
[1286]1602    dblParam_[CbcSmallestChange] = COIN_DBL_MAX;
1603    dblParam_[CbcSumChange] = 0.0;
1604    dblParam_[CbcLargestChange] = 0.0;
1605    intParam_[CbcNumberBranches] = 0;
[1798]1606    // Force minimization !!!!
1607    bool flipObjective = (solver_->getObjSense()<0.0);
1608    if (flipObjective)
1609      flipModel();
1610    dblParam_[CbcOptimizationDirection] = 1.0; // was solver_->getObjSense();
[1286]1611    strongInfo_[0] = 0;
1612    strongInfo_[1] = 0;
1613    strongInfo_[2] = 0;
1614    strongInfo_[3] = 0;
1615    strongInfo_[4] = 0;
1616    strongInfo_[5] = 0;
1617    strongInfo_[6] = 0;
1618    numberStrongIterations_ = 0;
1619    currentNode_ = NULL;
1620    // See if should do cuts old way
1621    if (parallelMode() < 0) {
1622        specialOptions_ |= 4096 + 8192;
1623    } else if (parallelMode() > 0) {
1624        specialOptions_ |= 4096;
1625    }
1626    int saveMoreSpecialOptions = moreSpecialOptions_;
1627    if (dynamic_cast<CbcTreeLocal *> (tree_))
1628        specialOptions_ |= 4096 + 8192;
[838]1629#ifdef COIN_HAS_CLP
[1286]1630    {
1631        OsiClpSolverInterface * clpSolver
1632        = dynamic_cast<OsiClpSolverInterface *> (solver_);
1633        if (clpSolver) {
1634            // pass in disaster handler
1635            CbcDisasterHandler handler(this);
1636            clpSolver->passInDisasterHandler(&handler);
[1798]1637            // Initialise solvers seed (unless users says not)
1638            if ((specialOptions_&4194304)==0)
1639              clpSolver->getModelPtr()->setRandomSeed(1234567);
[1393]1640#ifdef JJF_ZERO
[1286]1641            // reduce factorization frequency
1642            int frequency = clpSolver->getModelPtr()->factorizationFrequency();
1643            clpSolver->getModelPtr()->setFactorizationFrequency(CoinMin(frequency, 120));
[1016]1644#endif
[1286]1645        }
1646    }
[838]1647#endif
[1286]1648    // original solver (only set if pre-processing)
1649    OsiSolverInterface * originalSolver = NULL;
1650    int numberOriginalObjects = numberObjects_;
1651    OsiObject ** originalObject = NULL;
1652    // Save whether there were any objects
1653    bool noObjects = (numberObjects_ == 0);
1654    // Set up strategies
[1409]1655    /*
1656      See if the user has supplied a strategy object and deal with it if present.
1657      The call to setupOther will set numberStrong_ and numberBeforeTrust_, and
1658      perform integer preprocessing, if requested.
[1364]1659
[1409]1660      We need to hang on to a pointer to solver_. setupOther will assign a
1661      preprocessed solver to model, but will instruct assignSolver not to trash the
1662      existing one.
1663    */
[1357]1664    if (strategy_) {
1665        // May do preprocessing
1666        originalSolver = solver_;
1667        strategy_->setupOther(*this);
1668        if (strategy_->preProcessState()) {
1669            // pre-processing done
1670            if (strategy_->preProcessState() < 0) {
[1883]1671                // infeasible (or unbounded)
[1357]1672                status_ = 0 ;
[1883]1673                if (!solver_->isProvenDualInfeasible()) {
1674                  handler_->message(CBC_INFEAS, messages_) << CoinMessageEol ;
1675                  secondaryStatus_ = 1;
1676                } else {
1677                  handler_->message(CBC_UNBOUNDED, 
1678                                    messages_) << CoinMessageEol ;
1679                  secondaryStatus_ = 7;
1680                }
[1357]1681                originalContinuousObjective_ = COIN_DBL_MAX;
[1798]1682                if (flipObjective)
1683                  flipModel();
[1357]1684                return ;
1685            } else if (numberObjects_ && object_) {
1686                numberOriginalObjects = numberObjects_;
1687                // redo sequence
1688                numberIntegers_ = 0;
1689                int numberColumns = getNumCols();
1690                int nOrig = originalSolver->getNumCols();
1691                CglPreProcess * process = strategy_->process();
1692                assert (process);
1693                const int * originalColumns = process->originalColumns();
1694                // allow for cliques etc
1695                nOrig = CoinMax(nOrig, originalColumns[numberColumns-1] + 1);
1696                // try and redo debugger
1697                OsiRowCutDebugger * debugger = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
[1839]1698                if (debugger) {
1699                  if (numberColumns<=debugger->numberColumns())
[1357]1700                    debugger->redoSolution(numberColumns, originalColumns);
[1839]1701                  else
1702                    debugger=NULL; // no idea how to handle (SOS?)
1703                }
[1409]1704                // User-provided solution might have been best. Synchronise.
[1357]1705                if (bestSolution_) {
1706                    // need to redo - in case no better found in BAB
1707                    // just get integer part right
1708                    for (int i = 0; i < numberColumns; i++) {
1709                        int jColumn = originalColumns[i];
1710                        bestSolution_[i] = bestSolution_[jColumn];
1711                    }
1712                }
1713                originalObject = object_;
1714                // object number or -1
1715                int * temp = new int[nOrig];
1716                int iColumn;
1717                for (iColumn = 0; iColumn < nOrig; iColumn++)
1718                    temp[iColumn] = -1;
1719                int iObject;
1720                int nNonInt = 0;
1721                for (iObject = 0; iObject < numberOriginalObjects; iObject++) {
1722                    iColumn = originalObject[iObject]->columnNumber();
1723                    if (iColumn < 0) {
1724                        nNonInt++;
1725                    } else {
1726                        temp[iColumn] = iObject;
1727                    }
1728                }
1729                int numberNewIntegers = 0;
1730                int numberOldIntegers = 0;
1731                int numberOldOther = 0;
1732                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1733                    int jColumn = originalColumns[iColumn];
1734                    if (temp[jColumn] >= 0) {
1735                        int iObject = temp[jColumn];
1736                        CbcSimpleInteger * obj =
1737                            dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ;
1738                        if (obj)
1739                            numberOldIntegers++;
1740                        else
1741                            numberOldOther++;
1742                    } else if (isInteger(iColumn)) {
1743                        numberNewIntegers++;
1744                    }
1745                }
1746                /*
1747                  Allocate an array to hold the indices of the integer variables.
1748                  Make a large enough array for all objects
1749                */
1750                numberObjects_ = numberNewIntegers + numberOldIntegers + numberOldOther + nNonInt;
1751                object_ = new OsiObject * [numberObjects_];
1752                delete [] integerVariable_;
1753                integerVariable_ = new int [numberNewIntegers+numberOldIntegers];
1754                /*
1755                  Walk the variables again, filling in the indices and creating objects for
1756                  the integer variables. Initially, the objects hold the index and upper &
1757                  lower bounds.
1758                */
1759                numberIntegers_ = 0;
1760                int n = originalColumns[numberColumns-1] + 1;
1761                int * backward = new int[n];
1762                int i;
1763                for ( i = 0; i < n; i++)
1764                    backward[i] = -1;
1765                for (i = 0; i < numberColumns; i++)
1766                    backward[originalColumns[i]] = i;
1767                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1768                    int jColumn = originalColumns[iColumn];
1769                    if (temp[jColumn] >= 0) {
1770                        int iObject = temp[jColumn];
1771                        CbcSimpleInteger * obj =
1772                            dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ;
1773                        if (obj) {
1774                            object_[numberIntegers_] = originalObject[iObject]->clone();
1775                            // redo ids etc
1776                            //object_[numberIntegers_]->resetSequenceEtc(numberColumns,originalColumns);
1777                            object_[numberIntegers_]->resetSequenceEtc(numberColumns, backward);
1778                            integerVariable_[numberIntegers_++] = iColumn;
1779                        }
1780                    } else if (isInteger(iColumn)) {
1781                        object_[numberIntegers_] =
1782                            new CbcSimpleInteger(this, iColumn);
1783                        integerVariable_[numberIntegers_++] = iColumn;
1784                    }
1785                }
1786                delete [] backward;
1787                numberObjects_ = numberIntegers_;
1788                // Now append other column stuff
1789                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1790                    int jColumn = originalColumns[iColumn];
1791                    if (temp[jColumn] >= 0) {
1792                        int iObject = temp[jColumn];
1793                        CbcSimpleInteger * obj =
1794                            dynamic_cast <CbcSimpleInteger *>(originalObject[iObject]) ;
1795                        if (!obj) {
1796                            object_[numberObjects_] = originalObject[iObject]->clone();
1797                            // redo ids etc
1798                            CbcObject * obj =
1799                                dynamic_cast <CbcObject *>(object_[numberObjects_]) ;
1800                            assert (obj);
1801                            obj->redoSequenceEtc(this, numberColumns, originalColumns);
1802                            numberObjects_++;
1803                        }
1804                    }
1805                }
1806                // now append non column stuff
1807                for (iObject = 0; iObject < numberOriginalObjects; iObject++) {
1808                    iColumn = originalObject[iObject]->columnNumber();
1809                    if (iColumn < 0) {
1810                        // already has column numbers changed
1811                        object_[numberObjects_] = originalObject[iObject]->clone();
[1393]1812#ifdef JJF_ZERO
[1357]1813                        // redo ids etc
1814                        CbcObject * obj =
1815                            dynamic_cast <CbcObject *>(object_[numberObjects_]) ;
1816                        assert (obj);
1817                        obj->redoSequenceEtc(this, numberColumns, originalColumns);
1818#endif
1819                        numberObjects_++;
1820                    }
1821                }
1822                delete [] temp;
1823                if (!numberObjects_)
1824                    handler_->message(CBC_NOINT, messages_) << CoinMessageEol ;
1825            } else {
1826                int numberColumns = getNumCols();
1827                CglPreProcess * process = strategy_->process();
1828                assert (process);
1829                const int * originalColumns = process->originalColumns();
1830                // try and redo debugger
1831                OsiRowCutDebugger * debugger = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
1832                if (debugger)
1833                    debugger->redoSolution(numberColumns, originalColumns);
1834            }
1835        } else {
1836            //no preprocessing
1837            originalSolver = NULL;
1838        }
1839        strategy_->setupCutGenerators(*this);
1840        strategy_->setupHeuristics(*this);
1841        // Set strategy print level to models
1842        strategy_->setupPrinting(*this, handler_->logLevel());
1843    }
[1286]1844    eventHappened_ = false;
1845    CbcEventHandler *eventHandler = getEventHandler() ;
1846    if (eventHandler)
1847        eventHandler->setModel(this);
[1016]1848#define CLIQUE_ANALYSIS
[833]1849#ifdef CLIQUE_ANALYSIS
[1286]1850    // set up for probing
[1387]1851    // If we're doing clever stuff with cliques, additional info here.
[1286]1852    if (!parentModel_)
1853        probingInfo_ = new CglTreeProbingInfo(solver_);
1854    else
1855        probingInfo_ = NULL;
[833]1856#else
[1286]1857    probingInfo_ = NULL;
[833]1858#endif
[640]1859
[1286]1860    // Try for dominated columns
1861    if ((specialOptions_&64) != 0) {
1862        CglDuplicateRow dupcuts(solver_);
1863        dupcuts.setMode(2);
1864        CglStored * storedCuts = dupcuts.outDuplicates(solver_);
1865        if (storedCuts) {
[1641]1866          COIN_DETAIL_PRINT(printf("adding dup cuts\n"));
[1286]1867            addCutGenerator(storedCuts, 1, "StoredCuts from dominated",
1868                            true, false, false, -200);
1869        }
[1271]1870    }
[1286]1871    if (!nodeCompare_)
1872        nodeCompare_ = new CbcCompareDefault();;
1873    // See if hot start wanted
1874    CbcCompareBase * saveCompare = NULL;
[1387]1875    // User supplied hotstart. Adapt for preprocessing.
[1286]1876    if (hotstartSolution_) {
1877        if (strategy_ && strategy_->preProcessState() > 0) {
1878            CglPreProcess * process = strategy_->process();
1879            assert (process);
1880            int n = solver_->getNumCols();
1881            const int * originalColumns = process->originalColumns();
1882            // columns should be in order ... but
1883            double * tempS = new double[n];
1884            for (int i = 0; i < n; i++) {
1885                int iColumn = originalColumns[i];
1886                tempS[i] = hotstartSolution_[iColumn];
1887            }
1888            delete [] hotstartSolution_;
1889            hotstartSolution_ = tempS;
1890            if (hotstartPriorities_) {
1891                int * tempP = new int [n];
1892                for (int i = 0; i < n; i++) {
1893                    int iColumn = originalColumns[i];
1894                    tempP[i] = hotstartPriorities_[iColumn];
1895                }
1896                delete [] hotstartPriorities_;
1897                hotstartPriorities_ = tempP;
1898            }
[231]1899        }
[1286]1900        saveCompare = nodeCompare_;
1901        // depth first
1902        nodeCompare_ = new CbcCompareDepth();
[231]1903    }
[1286]1904    if (!problemFeasibility_)
1905        problemFeasibility_ = new CbcFeasibilityBase();
[2]1906# ifdef CBC_DEBUG
[1286]1907    std::string problemName ;
1908    solver_->getStrParam(OsiProbName, problemName) ;
1909    printf("Problem name - %s\n", problemName.c_str()) ;
1910    solver_->setHintParam(OsiDoReducePrint, false, OsiHintDo, 0) ;
[2]1911# endif
[1286]1912    /*
1913      Assume we're done, and see if we're proven wrong.
1914    */
1915    status_ = 0 ;
1916    secondaryStatus_ = 0;
1917    phase_ = 0;
1918    /*
1919      Scan the variables, noting the integer variables. Create an
1920      CbcSimpleInteger object for each integer variable.
1921    */
1922    findIntegers(false) ;
1923    // Say not dynamic pseudo costs
1924    ownership_ &= ~0x40000000;
1925    // If dynamic pseudo costs then do
1926    if (numberBeforeTrust_)
1927        convertToDynamic();
[1387]1928    // Set up char array to say if integer (speed)
[1286]1929    delete [] integerInfo_;
1930    {
1931        int n = solver_->getNumCols();
1932        integerInfo_ = new char [n];
1933        for (int i = 0; i < n; i++) {
1934            if (solver_->isInteger(i))
1935                integerInfo_[i] = 1;
1936            else
1937                integerInfo_[i] = 0;
1938        }
[197]1939    }
[1286]1940    if (preferredWay_) {
1941        // set all unset ones
1942        for (int iObject = 0 ; iObject < numberObjects_ ; iObject++) {
1943            CbcObject * obj =
1944                dynamic_cast <CbcObject *>(object_[iObject]) ;
1945            if (obj && !obj->preferredWay())
1946                obj->setPreferredWay(preferredWay_);
1947        }
[640]1948    }
[1286]1949    /*
1950      Ensure that objects on the lists of OsiObjects, heuristics, and cut
1951      generators attached to this model all refer to this model.
1952    */
1953    synchronizeModel() ;
1954    if (!solverCharacteristics_) {
1955        OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
1956        if (solverCharacteristics) {
1957            solverCharacteristics_ = solverCharacteristics;
1958        } else {
1959            // replace in solver
1960            OsiBabSolver defaultC;
1961            solver_->setAuxiliaryInfo(&defaultC);
1962            solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
1963        }
[395]1964    }
[862]1965
[1286]1966    solverCharacteristics_->setSolver(solver_);
1967    // Set so we can tell we are in initial phase in resolve
1968    continuousObjective_ = -COIN_DBL_MAX ;
1969    /*
1970      Solve the relaxation.
[2]1971
[1286]1972      Apparently there are circumstances where this will be non-trivial --- i.e.,
1973      we've done something since initialSolve that's trashed the solution to the
1974      continuous relaxation.
1975    */
1976    /* Tell solver we are in Branch and Cut
1977       Could use last parameter for subtle differences */
1978    solver_->setHintParam(OsiDoInBranchAndCut, true, OsiHintDo, NULL) ;
[1053]1979#ifdef COIN_HAS_CLP
[1286]1980    {
1981        OsiClpSolverInterface * clpSolver
1982        = dynamic_cast<OsiClpSolverInterface *> (solver_);
1983        if (clpSolver) {
1984            ClpSimplex * clpSimplex = clpSolver->getModelPtr();
1985            if ((specialOptions_&32) == 0) {
1986                // take off names
1987                clpSimplex->dropNames();
1988            }
1989            // no crunch if mostly continuous
1990            if ((clpSolver->specialOptions()&(1 + 8)) != (1 + 8)) {
1991                int numberColumns = solver_->getNumCols();
1992                if (numberColumns > 1000 && numberIntegers_*4 < numberColumns)
1993                    clpSolver->setSpecialOptions(clpSolver->specialOptions()&(~1));
1994            }
1995            //#define NO_CRUNCH
[1053]1996#ifdef NO_CRUNCH
[1286]1997            printf("TEMP switching off crunch\n");
1998            int iOpt = clpSolver->specialOptions();
1999            iOpt &= ~1;
2000            iOpt |= 65536;
2001            clpSolver->setSpecialOptions(iOpt);
[1053]2002#endif
[1286]2003        }
[1053]2004    }
2005#endif
[1286]2006    bool feasible;
[1656]2007    numberSolves_ = 0 ;
[1286]2008    // If NLP then we assume already solved outside branchAndbound
2009    if (!solverCharacteristics_->solverType() || solverCharacteristics_->solverType() == 4) {
2010        feasible = resolve(NULL, 0) != 0 ;
[480]2011    } else {
[1286]2012        // pick up given status
2013        feasible = (solver_->isProvenOptimal() &&
2014                    !solver_->isDualObjectiveLimitReached()) ;
[480]2015    }
[1286]2016    if (problemFeasibility_->feasible(this, 0) < 0) {
2017        feasible = false; // pretend infeasible
[640]2018    }
[1286]2019    numberSavedSolutions_ = 0;
2020    int saveNumberStrong = numberStrong_;
2021    int saveNumberBeforeTrust = numberBeforeTrust_;
2022    /*
2023      If the linear relaxation of the root is infeasible, bail out now. Otherwise,
2024      continue with processing the root node.
2025    */
2026    if (!feasible) {
2027        status_ = 0 ;
2028        if (!solver_->isProvenDualInfeasible()) {
2029            handler_->message(CBC_INFEAS, messages_) << CoinMessageEol ;
2030            secondaryStatus_ = 1;
2031        } else {
2032            handler_->message(CBC_UNBOUNDED, messages_) << CoinMessageEol ;
2033            secondaryStatus_ = 7;
2034        }
2035        originalContinuousObjective_ = COIN_DBL_MAX;
2036        solverCharacteristics_ = NULL;
[1798]2037        if (flipObjective)
2038          flipModel();
[1286]2039        return ;
[1654]2040    } else if (!numberObjects_ && (!strategy_ || strategy_->preProcessState() <= 0)) {
[1286]2041        // nothing to do
[1798]2042        if (flipObjective)
2043          flipModel();
[1286]2044        solverCharacteristics_ = NULL;
2045        bestObjective_ = solver_->getObjValue() * solver_->getObjSense();
2046        int numberColumns = solver_->getNumCols();
2047        delete [] bestSolution_;
2048        bestSolution_ = new double[numberColumns];
2049        CoinCopyN(solver_->getColSolution(), numberColumns, bestSolution_);
2050        return ;
[640]2051    }
[1409]2052    /*
2053      See if we're using the Osi side of the branching hierarchy. If so, either
2054      convert existing CbcObjects to OsiObjects, or generate them fresh. In the
2055      first case, CbcModel owns the objects on the object_ list. In the second
2056      case, the solver holds the objects and object_ simply points to the
2057      solver's list.
[1364]2058
[1409]2059      080417 The conversion code here (the block protected by `if (obj)') cannot
2060      possibly be correct. On the Osi side, descent is OsiObject -> OsiObject2 ->
2061      all other Osi object classes. On the Cbc side, it's OsiObject -> CbcObject
2062      -> all other Cbc object classes. It's structurally impossible for any Osi
2063      object to descend from CbcObject. The only thing I can see is that this is
2064      really dead code, and object detection is now handled from the Osi side.
2065    */
[1286]2066    // Convert to Osi if wanted
2067    //OsiBranchingInformation * persistentInfo = NULL;
2068    if (branchingMethod_ && branchingMethod_->chooseMethod()) {
2069        //persistentInfo = new OsiBranchingInformation(solver_);
2070        if (numberOriginalObjects) {
2071            for (int iObject = 0 ; iObject < numberObjects_ ; iObject++) {
2072                CbcObject * obj =
2073                    dynamic_cast <CbcObject *>(object_[iObject]) ;
2074                if (obj) {
2075                    CbcSimpleInteger * obj2 =
2076                        dynamic_cast <CbcSimpleInteger *>(obj) ;
2077                    if (obj2) {
2078                        // back to Osi land
2079                        object_[iObject] = obj2->osiObject();
2080                        delete obj;
2081                    } else {
2082                        OsiSimpleInteger * obj3 =
2083                            dynamic_cast <OsiSimpleInteger *>(obj) ;
2084                        if (!obj3) {
2085                            OsiSOS * obj4 =
2086                                dynamic_cast <OsiSOS *>(obj) ;
2087                            if (!obj4) {
2088                                CbcSOS * obj5 =
2089                                    dynamic_cast <CbcSOS *>(obj) ;
2090                                if (obj5) {
2091                                    // back to Osi land
2092                                    object_[iObject] = obj5->osiObject(solver_);
2093                                } else {
2094                                    printf("Code up CbcObject type in Osi land\n");
2095                                    abort();
2096                                }
2097                            }
2098                        }
2099                    }
2100                }
2101            }
2102            // and add to solver
2103            //if (!solver_->numberObjects()) {
2104            solver_->addObjects(numberObjects_, object_);
2105            //} else {
2106            //if (solver_->numberObjects()!=numberOriginalObjects) {
2107            //printf("should have trapped that solver has objects before\n");
2108            //abort();
2109            //}
2110            //}
2111        } else {
[1409]2112            /*
2113              As of 080104, findIntegersAndSOS is misleading --- the default OSI
2114              implementation finds only integers.
2115            */
[1286]2116            // do from solver
2117            deleteObjects(false);
2118            solver_->findIntegersAndSOS(false);
2119            numberObjects_ = solver_->numberObjects();
2120            object_ = solver_->objects();
2121            ownObjects_ = false;
2122        }
2123        branchingMethod_->chooseMethod()->setSolver(solver_);
[640]2124    }
[1387]2125    // take off heuristics if have to (some do not work with SOS, for example)
2126    // object should know what's safe.
[1286]2127    {
2128        int numberOdd = 0;
2129        int numberSOS = 0;
2130        for (int i = 0; i < numberObjects_; i++) {
2131            if (!object_[i]->canDoHeuristics())
2132                numberOdd++;
2133            CbcSOS * obj =
2134                dynamic_cast <CbcSOS *>(object_[i]) ;
2135            if (obj)
2136                numberSOS++;
2137        }
2138        if (numberOdd) {
2139            if (numberHeuristics_) {
2140                int k = 0;
2141                for (int i = 0; i < numberHeuristics_; i++) {
2142                    if (!heuristic_[i]->canDealWithOdd())
2143                        delete heuristic_[i];
2144                    else
2145                        heuristic_[k++] = heuristic_[i];
2146                }
2147                if (!k) {
2148                    delete [] heuristic_;
2149                    heuristic_ = NULL;
2150                }
2151                numberHeuristics_ = k;
2152                handler_->message(CBC_HEURISTICS_OFF, messages_) << numberOdd << CoinMessageEol ;
2153            }
[1883]2154            // If odd switch off AddIntegers
2155            specialOptions_ &= ~65536;
[1286]2156        } else if (numberSOS) {
2157            specialOptions_ |= 128; // say can do SOS in dynamic mode
2158            // switch off fast nodes for now
2159            fastNodeDepth_ = -1;
[1883]2160            moreSpecialOptions_ &= ~33554432; // no diving
[1286]2161        }
2162        if (numberThreads_ > 0) {
2163            // switch off fast nodes for now
2164            fastNodeDepth_ = -1;
2165        }
[931]2166    }
[1286]2167    // Save objective (just so user can access it)
[1569]2168    originalContinuousObjective_ = solver_->getObjValue()* solver_->getObjSense();
[1286]2169    bestPossibleObjective_ = originalContinuousObjective_;
2170    sumChangeObjective1_ = 0.0;
2171    sumChangeObjective2_ = 0.0;
2172    /*
2173      OsiRowCutDebugger knows an optimal answer for a subset of MIP problems.
2174      Assuming it recognises the problem, when called upon it will check a cut to
2175      see if it cuts off the optimal answer.
2176    */
2177    // If debugger exists set specialOptions_ bit
2178    if (solver_->getRowCutDebuggerAlways()) {
2179        specialOptions_ |= 1;
2180    }
[56]2181
2182# ifdef CBC_DEBUG
[1286]2183    if ((specialOptions_&1) == 0)
2184        solver_->activateRowCutDebugger(problemName.c_str()) ;
2185    if (solver_->getRowCutDebuggerAlways())
2186        specialOptions_ |= 1;
[2]2187# endif
2188
[1286]2189    /*
2190      Begin setup to process a feasible root node.
2191    */
2192    bestObjective_ = CoinMin(bestObjective_, 1.0e50) ;
2193    if (!bestSolution_) {
2194        numberSolutions_ = 0 ;
2195        numberHeuristicSolutions_ = 0 ;
2196    }
2197    stateOfSearch_ = 0;
2198    // Everything is minimization
2199    {
2200        // needed to sync cutoffs
2201        double value ;
2202        solver_->getDblParam(OsiDualObjectiveLimit, value) ;
2203        dblParam_[CbcCurrentCutoff] = value * solver_->getObjSense();
2204    }
2205    double cutoff = getCutoff() ;
2206    double direction = solver_->getObjSense() ;
2207    dblParam_[CbcOptimizationDirection] = direction;
2208    if (cutoff < 1.0e20 && direction < 0.0)
2209        messageHandler()->message(CBC_CUTOFF_WARNING1,
2210                                  messages())
2211        << cutoff << -cutoff << CoinMessageEol ;
2212    if (cutoff > bestObjective_)
2213        cutoff = bestObjective_ ;
2214    setCutoff(cutoff) ;
2215    /*
2216      We probably already have a current solution, but just in case ...
2217    */
2218    int numberColumns = getNumCols() ;
2219    if (!currentSolution_)
2220        currentSolution_ = new double[numberColumns] ;
2221    testSolution_ = currentSolution_;
2222    /*
2223      Create a copy of the solver, thus capturing the original (root node)
2224      constraint system (aka the continuous system).
2225    */
2226    continuousSolver_ = solver_->clone() ;
[152]2227
[1839]2228    // add cutoff as constraint if wanted
[1883]2229    if (cutoffRowNumber_==-2) {
[1839]2230      if (!parentModel_) {
2231        int numberColumns=solver_->getNumCols();
2232        double * obj = CoinCopyOfArray(solver_->getObjCoefficients(),numberColumns);
2233        int * indices = new int [numberColumns];
2234        int n=0;
2235        for (int i=0;i<numberColumns;i++) {
2236          if (obj[i]) {
2237            indices[n]=i;
2238            obj[n++]=obj[i];
2239          }
2240        }
2241        if (n) {
2242          double cutoff=getCutoff();
[1883]2243          // relax a little bit
2244          cutoff += 1.0e-4;
[1839]2245          double offset;
2246          solver_->getDblParam(OsiObjOffset, offset);
[1883]2247          cutoffRowNumber_ = solver_->getNumRows();
[1839]2248          solver_->addRow(n,indices,obj,-COIN_DBL_MAX,CoinMin(cutoff,1.0e25)+offset);
2249        } else {
2250          // no objective!
[1883]2251          cutoffRowNumber_ = -1;
[1839]2252        }
2253        delete [] indices;
2254        delete [] obj;
2255      } else {
2256        // switch off
[1883]2257        cutoffRowNumber_ = -1;
[1839]2258      }
2259    }
[1286]2260    numberRowsAtContinuous_ = getNumRows() ;
2261    solver_->saveBaseModel();
2262    /*
2263      Check the objective to see if we can deduce a nontrivial increment. If
2264      it's better than the current value for CbcCutoffIncrement, it'll be
2265      installed.
2266    */
2267    if (solverCharacteristics_->reducedCostsAccurate())
2268        analyzeObjective() ;
2269    {
2270        // may be able to change cutoff now
2271        double cutoff = getCutoff();
2272        double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
2273        if (cutoff > bestObjective_ - increment) {
2274            cutoff = bestObjective_ - increment ;
2275            setCutoff(cutoff) ;
2276        }
[1271]2277    }
[931]2278#ifdef COIN_HAS_CLP
[1286]2279    // Possible save of pivot method
2280    ClpDualRowPivot * savePivotMethod = NULL;
2281    {
2282        // pass tolerance and increment to solver
2283        OsiClpSolverInterface * clpSolver
2284        = dynamic_cast<OsiClpSolverInterface *> (solver_);
2285        if (clpSolver)
2286            clpSolver->setStuff(getIntegerTolerance(), getCutoffIncrement());
[1656]2287#ifdef CLP_RESOLVE
2288        if ((moreSpecialOptions_&1048576)!=0&&!parentModel_&&clpSolver) {
2289          resolveClp(clpSolver,0);
2290        }
2291#endif
[1286]2292    }
[931]2293#endif
[1286]2294    /*
2295      Set up for cut generation. addedCuts_ holds the cuts which are relevant for
2296      the active subproblem. whichGenerator will be used to record the generator
2297      that produced a given cut.
2298    */
2299    maximumWhich_ = 1000 ;
2300    delete [] whichGenerator_;
2301    whichGenerator_ = new int[maximumWhich_] ;
2302    memset(whichGenerator_, 0, maximumWhich_*sizeof(int));
2303    maximumNumberCuts_ = 0 ;
2304    currentNumberCuts_ = 0 ;
2305    delete [] addedCuts_ ;
2306    addedCuts_ = NULL ;
2307    OsiObject ** saveObjects = NULL;
2308    maximumRows_ = numberRowsAtContinuous_;
2309    currentDepth_ = 0;
2310    workingBasis_.resize(maximumRows_, numberColumns);
2311    /*
2312      Set up an empty heap and associated data structures to hold the live set
2313      (problems which require further exploration).
2314    */
2315    CbcCompareDefault * compareActual
2316    = dynamic_cast<CbcCompareDefault *> (nodeCompare_);
2317    if (compareActual) {
2318        compareActual->setBestPossible(direction*solver_->getObjValue());
2319        compareActual->setCutoff(getCutoff());
[1393]2320#ifdef JJF_ZERO
[1357]2321        if (false && !numberThreads_ && !parentModel_) {
[1286]2322            printf("CbcTreeArray ? threads ? parentArray\n");
2323            // Setup new style tree
2324            delete tree_;
2325            tree_ = new CbcTreeArray();
2326        }
[1341]2327#endif
[1132]2328    }
[1286]2329    tree_->setComparison(*nodeCompare_) ;
2330    /*
2331      Used to record the path from a node to the root of the search tree, so that
2332      we can then traverse from the root to the node when restoring a subproblem.
2333    */
2334    maximumDepth_ = 10 ;
2335    delete [] walkback_ ;
2336    walkback_ = new CbcNodeInfo * [maximumDepth_] ;
2337    lastDepth_ = 0;
2338    delete [] lastNodeInfo_ ;
2339    lastNodeInfo_ = new CbcNodeInfo * [maximumDepth_] ;
2340    delete [] lastNumberCuts_ ;
2341    lastNumberCuts_ = new int [maximumDepth_] ;
2342    maximumCuts_ = 100;
2343    lastNumberCuts2_ = 0;
2344    delete [] lastCut_;
2345    lastCut_ = new const OsiRowCut * [maximumCuts_];
2346    /*
2347      Used to generate bound edits for CbcPartialNodeInfo.
2348    */
2349    double * lowerBefore = new double [numberColumns] ;
2350    double * upperBefore = new double [numberColumns] ;
2351    /*
[1409]2352    Set up to run heuristics and generate cuts at the root node. The heavy
2353    lifting is hidden inside the calls to doHeuristicsAtRoot and solveWithCuts.
[2]2354
[1409]2355    To start, tell cut generators they can be a bit more aggressive at the
2356    root node.
[1364]2357
[1409]2358    QUESTION: phase_ = 0 is documented as `initial solve', phase = 1 as `solve
2359        with cuts at root'. Is phase_ = 1 the correct indication when
2360        doHeurisiticsAtRoot is called to run heuristics outside of the main
2361        cut / heurisitc / reoptimise loop in solveWithCuts?
[1364]2362
[1286]2363      Generate cuts at the root node and reoptimise. solveWithCuts does the heavy
2364      lifting. It will iterate a generate/reoptimise loop (including reduced cost
2365      fixing) until no cuts are generated, the change in objective falls off,  or
2366      the limit on the number of rounds of cut generation is exceeded.
[2]2367
[1286]2368      At the end of all this, any cuts will be recorded in cuts and also
2369      installed in the solver's constraint system. We'll have reoptimised, and
2370      removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
2371      adjusted accordingly).
[2]2372
[1286]2373      Tell cut generators they can be a bit more aggressive at root node
2374
2375      TODO: Why don't we make a copy of the solution after solveWithCuts?
2376      TODO: If numberUnsatisfied == 0, don't we have a solution?
2377    */
2378    phase_ = 1;
2379    int iCutGenerator;
2380    for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
[1409]2381        // If parallel switch off global cuts
2382        if (numberThreads_) {
2383            generator_[iCutGenerator]->setGlobalCuts(false);
2384            generator_[iCutGenerator]->setGlobalCutsAtRoot(false);
2385        }
[1286]2386        CglCutGenerator * generator = generator_[iCutGenerator]->generator();
2387        generator->setAggressiveness(generator->getAggressiveness() + 100);
[1839]2388        if (!generator->canDoGlobalCuts())
2389          generator->setGlobalCuts(false);
[1271]2390    }
[1286]2391    OsiCuts cuts ;
2392    int anyAction = -1 ;
2393    numberOldActiveCuts_ = 0 ;
2394    numberNewCuts_ = 0 ;
2395    // Array to mark solution
2396    delete [] usedInSolution_;
2397    usedInSolution_ = new int[numberColumns];
2398    CoinZeroN(usedInSolution_, numberColumns);
2399    /*
2400      For printing totals and for CbcNode (numberNodes_)
2401    */
2402    numberIterations_ = 0 ;
2403    numberNodes_ = 0 ;
2404    numberNodes2_ = 0 ;
2405    maximumStatistics_ = 0;
2406    maximumDepthActual_ = 0;
2407    numberDJFixed_ = 0.0;
[1315]2408    if (!parentModel_) {
2409        if ((specialOptions_&262144) != 0) {
2410            // create empty stored cuts
2411            //storedRowCuts_ = new CglStored(solver_->getNumCols());
2412        } else if ((specialOptions_&524288) != 0 && storedRowCuts_) {
2413            // tighten and set best solution
2414            // A) tight bounds on integer variables
[1409]2415            /*
2416                storedRowCuts_ are coming in from outside, probably for nonlinear.
2417              John was unsure about origin.
2418            */
[1315]2419            const double * lower = solver_->getColLower();
2420            const double * upper = solver_->getColUpper();
2421            const double * tightLower = storedRowCuts_->tightLower();
2422            const double * tightUpper = storedRowCuts_->tightUpper();
2423            int nTightened = 0;
2424            for (int i = 0; i < numberIntegers_; i++) {
2425                int iColumn = integerVariable_[i];
2426                if (tightLower[iColumn] > lower[iColumn]) {
2427                    nTightened++;
2428                    solver_->setColLower(iColumn, tightLower[iColumn]);
2429                }
2430                if (tightUpper[iColumn] < upper[iColumn]) {
2431                    nTightened++;
2432                    solver_->setColUpper(iColumn, tightUpper[iColumn]);
2433                }
2434            }
2435            if (nTightened)
[1641]2436              COIN_DETAIL_PRINT(printf("%d tightened by alternate cuts\n", nTightened));
[1315]2437            if (storedRowCuts_->bestObjective() < bestObjective_) {
2438                // B) best solution
2439                double objValue = storedRowCuts_->bestObjective();
2440                setBestSolution(CBC_SOLUTION, objValue,
2441                                storedRowCuts_->bestSolution()) ;
2442                // Do heuristics
2443                // Allow RINS
2444                for (int i = 0; i < numberHeuristics_; i++) {
2445                    CbcHeuristicRINS * rins
2446                    = dynamic_cast<CbcHeuristicRINS *> (heuristic_[i]);
2447                    if (rins) {
2448                        rins->setLastNode(-100);
2449                    }
2450                }
2451            }
2452        }
2453    }
[1409]2454    /*
2455      Run heuristics at the root. This is the only opportunity to run FPump; it
2456      will be removed from the heuristics list by doHeuristicsAtRoot.
2457    */
[1839]2458    // See if multiple runs wanted
2459    CbcModel ** rootModels=NULL;
2460    if (!parentModel_&&multipleRootTries_%100) {
2461      double rootTimeCpu=CoinCpuTime();
2462      double startTimeRoot=CoinGetTimeOfDay();
2463      int numberRootThreads=1;
2464      /* undocumented fine tuning
2465         aabbcc where cc is number of tries
2466         bb if nonzero is number of threads
2467         aa if nonzero just do heuristics
2468      */
2469      int numberModels = multipleRootTries_%100;
2470#ifdef CBC_THREAD
2471      numberRootThreads = (multipleRootTries_/100)%100;
2472      if (!numberRootThreads)
2473        numberRootThreads=numberModels;
2474#endif
2475      int otherOptions = (multipleRootTries_/10000)%100;
2476      rootModels = new CbcModel * [numberModels];
2477      unsigned int newSeed = randomSeed_;
2478      if (newSeed==0) {
2479        double time = fabs(CoinGetTimeOfDay());
2480        while (time>=COIN_INT_MAX)
2481          time *= 0.5;
2482        newSeed = static_cast<unsigned int>(time);
2483      } else if (newSeed<0) {
2484        newSeed = 123456789;
2485#ifdef COIN_HAS_CLP
2486        OsiClpSolverInterface * clpSolver
2487          = dynamic_cast<OsiClpSolverInterface *> (solver_);
2488        if (clpSolver) {
2489          newSeed += clpSolver->getModelPtr()->randomNumberGenerator()->getSeed();
2490        }
2491#endif
2492      }
2493      CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis *> (solver_->getEmptyWarmStart());
[1883]2494      int numberIterations=0;
[1839]2495      for (int i=0;i<numberModels;i++) { 
2496        rootModels[i]=new CbcModel(*this);
2497        rootModels[i]->setNumberThreads(0);
2498        rootModels[i]->setMaximumNodes(otherOptions ? -1 : 0);
2499        rootModels[i]->setRandomSeed(newSeed+10000000*i);
2500        rootModels[i]->randomNumberGenerator()->setSeed(newSeed+50000000*i);
2501        rootModels[i]->setMultipleRootTries(0);
2502        // use seed
2503        rootModels[i]->setSpecialOptions(specialOptions_ |(4194304|8388608));
[1883]2504        rootModels[i]->setMoreSpecialOptions(moreSpecialOptions_ & 
2505                                             (~134217728));
[1839]2506        rootModels[i]->solver_->setWarmStart(basis);
2507#ifdef COIN_HAS_CLP
2508        OsiClpSolverInterface * clpSolver
2509          = dynamic_cast<OsiClpSolverInterface *> (rootModels[i]->solver_);
2510        if (clpSolver) {
[1883]2511          ClpSimplex * simplex = clpSolver->getModelPtr();
2512          if (defaultHandler_)
2513            simplex->setDefaultMessageHandler();
2514          simplex->setRandomSeed(newSeed+20000000*i);
2515          simplex->allSlackBasis();
2516          int logLevel=simplex->logLevel();
2517          if (logLevel==1)
2518            simplex->setLogLevel(0);
2519          if (i!=0) {
2520            double random = simplex->randomNumberGenerator()->randomDouble();
2521            int bias = static_cast<int>(random*(numberIterations/4));
2522            simplex->setMaximumIterations(numberIterations/2+bias);
2523            simplex->primal();
2524            simplex->setMaximumIterations(COIN_INT_MAX);
2525            simplex->dual();
2526          } else {
2527            simplex->primal();
2528            numberIterations=simplex->numberIterations();
2529          }
2530          simplex->setLogLevel(logLevel);
2531          clpSolver->setWarmStart(NULL);
[1839]2532        }
2533#endif
[1883]2534        for (int j=0;j<numberHeuristics_;j++)
2535          rootModels[i]->heuristic_[j]->setSeed(rootModels[i]->heuristic_[j]->getSeed()+100000000*i);
2536        for (int j=0;j<numberCutGenerators_;j++)
2537          rootModels[i]->generator_[j]->generator()->refreshSolver(rootModels[i]->solver_);
[1839]2538      }
2539      delete basis;
2540#ifdef CBC_THREAD
2541      if (numberRootThreads==1) {
2542#endif
2543        for (int iModel=0;iModel<numberModels;iModel++) {
2544          doRootCbcThread(rootModels[iModel]);
2545          // see if solved at root node
2546          if (rootModels[iModel]->getMaximumNodes()) {
2547            feasible=false;
2548            break;
2549          }
2550        }
2551#ifdef CBC_THREAD
2552      } else {
2553        Coin_pthread_t * threadId = new Coin_pthread_t [numberRootThreads];
2554        for (int kModel=0;kModel<numberModels;kModel+=numberRootThreads) {
2555          bool finished=false;
2556          for (int iModel=kModel;iModel<CoinMin(numberModels,kModel+numberRootThreads);iModel++) {
2557            pthread_create(&(threadId[iModel-kModel].thr), NULL, 
2558                           doRootCbcThread,
2559                           rootModels[iModel]);
2560          }
2561          // wait
2562          for (int iModel=kModel;iModel<CoinMin(numberModels,kModel+numberRootThreads);iModel++) {
2563            pthread_join(threadId[iModel-kModel].thr, NULL);
2564          }
2565          // see if solved at root node
2566          for (int iModel=kModel;iModel<CoinMin(numberModels,kModel+numberRootThreads);iModel++) {
2567            if (rootModels[iModel]->getMaximumNodes())
2568              finished=true;
2569          }
2570          if (finished) {
2571            feasible=false;
2572            break;
2573          }
2574        }
2575        delete [] threadId;
2576      }
2577#endif
2578      // sort solutions
2579      int * which = new int [numberModels];
2580      double * value = new double [numberModels];
2581      int numberSolutions=0;
2582      for (int iModel=0;iModel<numberModels;iModel++) {
2583        if (rootModels[iModel]->bestSolution()) {
2584          which[numberSolutions]=iModel;
2585          value[numberSolutions++]=
2586            -rootModels[iModel]->getMinimizationObjValue();
2587        }
2588      }
2589      char general[100];
2590      rootTimeCpu=CoinCpuTime()-rootTimeCpu;
2591      if (numberRootThreads==1)
2592        sprintf(general,"Multiple root solvers took a total of %.2f seconds\n",
2593                rootTimeCpu);
2594      else
2595        sprintf(general,"Multiple root solvers took a total of %.2f seconds (%.2f elapsed)\n",
2596                rootTimeCpu,CoinGetTimeOfDay()-startTimeRoot);
2597      messageHandler()->message(CBC_GENERAL,
2598                                messages())
2599        << general << CoinMessageEol ;
2600      CoinSort_2(value,value+numberSolutions,which);
2601      // to get name
2602      CbcHeuristicRINS dummyHeuristic;
2603      dummyHeuristic.setHeuristicName("Multiple root solvers");
2604      lastHeuristic_=&dummyHeuristic;
2605      for (int i=0;i<numberSolutions;i++) {
2606        double objValue = - value[i];
2607        if (objValue<getCutoff()) {
2608          int iModel=which[i];
2609          setBestSolution(CBC_ROUNDING,objValue,
2610                          rootModels[iModel]->bestSolution());
2611        }
2612      }
2613      lastHeuristic_=NULL;
2614      delete [] which;
2615      delete [] value;
2616    }
[1286]2617    // Do heuristics
[1839]2618    if (numberObjects_&&!rootModels)
[1654]2619        doHeuristicsAtRoot();
[1682]2620    if (solverCharacteristics_->solutionAddsCuts()) {
2621      // With some heuristics solver needs a resolve here
2622      solver_->resolve();
2623      if(!isProvenOptimal()){
2624        solver_->initialSolve();
2625      }
[1663]2626    }
[1409]2627    /*
2628      Grepping through the code, it would appear that this is a command line
2629      debugging hook.  There's no obvious place in the code where this is set to
2630      a negative value.
[1387]2631
[1409]2632      User hook, says John.
2633    */
[1802]2634    if ( intParam_[CbcMaxNumNode] < 0
2635      ||numberSolutions_>=getMaximumSolutions())
2636      eventHappened_ = true; // stop as fast as possible
[1286]2637    stoppedOnGap_ = false ;
2638    // See if can stop on gap
2639    bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
[1717]2640    if(canStopOnGap()) {
[1286]2641        if (bestPossibleObjective_ < getCutoff())
2642            stoppedOnGap_ = true ;
2643        feasible = false;
2644        //eventHappened_=true; // stop as fast as possible
[1271]2645    }
[1409]2646    /*
2647      Set up for statistics collection, if requested. Standard values are
2648      documented in CbcModel.hpp. The magic number 100 will trigger a dump of
2649      CbcSimpleIntegerDynamicPseudoCost objects (no others). Looks like another
2650      command line debugging hook.
2651    */
[1286]2652    statistics_ = NULL;
2653    // Do on switch
2654    if (doStatistics > 0 && doStatistics <= 100) {
2655        maximumStatistics_ = 10000;
2656        statistics_ = new CbcStatistics * [maximumStatistics_];
2657        memset(statistics_, 0, maximumStatistics_*sizeof(CbcStatistics *));
[1271]2658    }
[1286]2659    // See if we can add integers
[1330]2660    if (noObjects && numberIntegers_ < solver_->getNumCols() && (specialOptions_&65536) != 0 && !parentModel_)
[1357]2661        AddIntegers();
[1409]2662    /*
2663      Do an initial round of cut generation for the root node. Depending on the
2664      type of underlying solver, we may want to do this even if the initial query
2665      to the objects indicates they're satisfied.
[1330]2666
[1409]2667      solveWithCuts does the heavy lifting. It will iterate a generate/reoptimise
2668      loop (including reduced cost fixing) until no cuts are generated, the
2669      change in objective falls off,  or the limit on the number of rounds of cut
2670      generation is exceeded.
[1364]2671
[1409]2672      At the end of all this, any cuts will be recorded in cuts and also
2673      installed in the solver's constraint system. We'll have reoptimised, and
2674      removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
2675      adjusted accordingly).
2676    */
[1330]2677    int iObject ;
2678    int preferredWay ;
2679    int numberUnsatisfied = 0 ;
2680    delete [] currentSolution_;
2681    currentSolution_ = new double [numberColumns];
2682    testSolution_ = currentSolution_;
2683    memcpy(currentSolution_, solver_->getColSolution(),
2684           numberColumns*sizeof(double)) ;
2685    // point to useful information
2686    OsiBranchingInformation usefulInfo = usefulInformation();
2687
2688    for (iObject = 0 ; iObject < numberObjects_ ; iObject++) {
2689        double infeasibility =
2690            object_[iObject]->infeasibility(&usefulInfo, preferredWay) ;
2691        if (infeasibility ) numberUnsatisfied++ ;
2692    }
2693    // replace solverType
[1883]2694    double * tightBounds = NULL;
[1330]2695    if (solverCharacteristics_->tryCuts())  {
2696
2697        if (numberUnsatisfied)   {
2698            // User event
2699            if (!eventHappened_ && feasible) {
[1839]2700                if (rootModels) {
[1883]2701                  // for fixings
2702                  int numberColumns=solver_->getNumCols();
2703                  tightBounds = new double [2*numberColumns];
2704                  {
2705                    const double * lower = solver_->getColLower();
2706                    const double * upper = solver_->getColUpper();
2707                    for (int i=0;i<numberColumns;i++) {
2708                      tightBounds[2*i+0]=lower[i];
2709                      tightBounds[2*i+1]=upper[i];
2710                    }
2711                  }
[1839]2712                  int numberModels = multipleRootTries_%100;
2713                  const OsiSolverInterface ** solvers = new 
2714                    const OsiSolverInterface * [numberModels];
2715                  int numberRows=continuousSolver_->getNumRows();
2716                  int maxCuts=0;
2717                  for (int i=0;i<numberModels;i++) {
2718                    solvers[i]=rootModels[i]->solver();
[1883]2719                    const double * lower = solvers[i]->getColLower();
2720                    const double * upper = solvers[i]->getColUpper();
2721                    for (int j=0;j<numberColumns;j++) {
2722                      tightBounds[2*j+0]=CoinMax(lower[j],tightBounds[2*j+0]);
2723                      tightBounds[2*j+1]=CoinMin(upper[j],tightBounds[2*j+1]);
2724                    }
[1839]2725                    int numberRows2=solvers[i]->getNumRows();
2726                    assert (numberRows2>=numberRows);
2727                    maxCuts += numberRows2-numberRows;
2728                    // accumulate statistics
2729                    for (int j=0;j<numberCutGenerators_;j++) {
2730                      generator_[j]->addStatistics(rootModels[i]->cutGenerator(j));
2731                    }
2732                  }
2733                  for (int j=0;j<numberCutGenerators_;j++) {
2734                    generator_[j]->scaleBackStatistics(numberModels);
2735                  }
[1883]2736                  //CbcRowCuts rowCut(maxCuts);
2737                  const OsiRowCutDebugger *debugger = NULL;
2738                  if ((specialOptions_&1) != 0) 
2739                    debugger = solver_->getRowCutDebugger() ;
[1839]2740                  for (int iModel=0;iModel<numberModels;iModel++) {
2741                    int numberRows2=solvers[iModel]->getNumRows();
2742                    const CoinPackedMatrix * rowCopy = solvers[iModel]->getMatrixByRow();
2743                    const int * rowLength = rowCopy->getVectorLengths(); 
2744                    const double * elements = rowCopy->getElements();
2745                    const int * column = rowCopy->getIndices();
2746                    const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
2747                    const double * rowLower = solvers[iModel]->getRowLower();
2748                    const double * rowUpper = solvers[iModel]->getRowUpper();
2749                    for (int iRow=numberRows;iRow<numberRows2;iRow++) {
2750                      OsiRowCut rc;
2751                      rc.setLb(rowLower[iRow]);
2752                      rc.setUb(rowUpper[iRow]);
2753                      CoinBigIndex start = rowStart[iRow];
2754                      rc.setRow(rowLength[iRow],column+start,elements+start,false);
[1883]2755                      if (debugger)
2756                        CoinAssert (!debugger->invalidCut(rc));
2757                      globalCuts_.addCutIfNotDuplicate(rc);
[1839]2758                    }
[1883]2759                    //int cutsAdded=globalCuts_.numberCuts()-numberCuts;
[1839]2760                    //numberCuts += cutsAdded;
2761                    //printf("Model %d gave %d cuts (out of %d possible)\n",
2762                    //     iModel,cutsAdded,numberRows2-numberRows);
2763                  }
[1883]2764                  // normally replace global cuts
2765                  //if (!globalCuts_.())
2766                  //globalCuts_=rowCutrowCut.addCuts(globalCuts_);
2767                  //rowCut.addCuts(globalCuts_);
2768                  int nTightened=0;
2769                  bool feasible=true;
2770                  {
2771                    double tolerance=1.0e-5;
2772                    const double * lower = solver_->getColLower();
2773                    const double * upper = solver_->getColUpper();
2774                    for (int i=0;i<numberColumns;i++) {
2775                      if (tightBounds[2*i+0]>tightBounds[2*i+1]) {
2776                        feasible=false;
2777                        printf("Bad bounds on %d %g,%g was %g,%g\n",
2778                               i,tightBounds[2*i+0],tightBounds[2*i+1],
2779                               lower[i],upper[i]);
2780                      }
2781                      //int k=0;
2782                      double oldLower=lower[i];
2783                      double oldUpper=upper[i];
2784                      if (tightBounds[2*i+0]>oldLower+tolerance) {
2785                        nTightened++;
2786                        //k++;
2787                        solver_->setColLower(i,tightBounds[2*i+0]);
2788                      }
2789                      if (tightBounds[2*i+1]<oldUpper-tolerance) {
2790                        nTightened++;
2791                        //k++;
2792                        solver_->setColUpper(i,tightBounds[2*i+1]);
2793                      }
2794                      //if (k)
2795                      //printf("new bounds on %d %g,%g was %g,%g\n",
2796                      //       i,tightBounds[2*i+0],tightBounds[2*i+1],
2797                      //       oldLower,oldUpper);
2798                    }
2799                    if (!feasible)
2800                      abort(); // deal with later
2801                  }
2802                  delete [] tightBounds;
2803                  tightBounds=NULL;
[1839]2804                  char printBuffer[200];
2805                  sprintf(printBuffer,"%d solvers added %d different cuts out of pool of %d",
[1883]2806                          numberModels,globalCuts_.sizeRowCuts(),maxCuts);
[1839]2807                  messageHandler()->message(CBC_GENERAL, messages())
2808                    << printBuffer << CoinMessageEol ;
[1883]2809                  if (nTightened) {
2810                    sprintf(printBuffer,"%d bounds were tightened",
2811                          nTightened);
2812                    messageHandler()->message(CBC_GENERAL, messages())
2813                      << printBuffer << CoinMessageEol ;
2814                  }
[1839]2815                  delete [] solvers;
2816                }
[1883]2817                if (!parentModel_&&(moreSpecialOptions_&67108864) != 0) {
2818                  // load cuts from file
2819                  FILE * fp = fopen("global.cuts","rb");
2820                  if (fp) {
2821                    size_t nRead;
2822                    int numberColumns=solver_->getNumCols();
2823                    int numCols;
2824                    nRead = fread(&numCols, sizeof(int), 1, fp);
2825                    if (nRead != 1)
2826                      throw("Error in fread");
2827                    if (numberColumns!=numCols) {
2828                      printf("Mismatch on columns %d %d\n",numberColumns,numCols);
2829                      fclose(fp);
2830                    } else {
2831                      // If rootModel just do some
2832                      double threshold=-1.0;
2833                      if (!multipleRootTries_)
2834                        threshold=0.5;
2835                      int initialCuts=0;
2836                      int initialGlobal = globalCuts_.sizeRowCuts();
2837                      double * elements = new double [numberColumns+2];
2838                      int * indices = new int [numberColumns];
2839                      int numberEntries=1;
2840                      while (numberEntries>0) {
2841                        nRead = fread(&numberEntries, sizeof(int), 1, fp);
2842                        if (nRead != 1)
2843                          throw("Error in fread");
2844                        double randomNumber=randomNumberGenerator_.randomDouble();
2845                        if (numberEntries>0) {
2846                          initialCuts++;
2847                          nRead = fread(elements, sizeof(double), numberEntries+2, fp);
2848                          if (nRead != static_cast<size_t>(numberEntries+2))
2849                            throw("Error in fread");
2850                          nRead = fread(indices, sizeof(int), numberEntries, fp);
2851                          if (nRead != static_cast<size_t>(numberEntries))
2852                            throw("Error in fread");
2853                          if (randomNumber>threshold) {
2854                            OsiRowCut rc;
2855                            rc.setLb(elements[numberEntries]);
2856                            rc.setUb(elements[numberEntries+1]);
2857                            rc.setRow(numberEntries,indices,elements,
2858                                      false);
2859                            rc.setGloballyValidAsInteger(2);
2860                            globalCuts_.addCutIfNotDuplicate(rc) ;
2861                          } 
2862                        }
2863                      }
2864                      fclose(fp);
2865                      // fixes
2866                      int nTightened=0;
2867                      fp = fopen("global.fix","rb");
2868                      if (fp) {
2869                        nRead = fread(indices, sizeof(int), 2, fp);
2870                        if (nRead != 2)
2871                          throw("Error in fread");
2872                        if (numberColumns!=indices[0]) {
2873                          printf("Mismatch on columns %d %d\n",numberColumns,
2874                                 indices[0]);
2875                        } else {
2876                          indices[0]=1;
2877                          while (indices[0]>=0) {
2878                            nRead = fread(indices, sizeof(int), 2, fp);
2879                            if (nRead != 2)
2880                              throw("Error in fread");
2881                            int iColumn=indices[0];
2882                            if (iColumn>=0) {
2883                              nTightened++;
2884                              nRead = fread(elements, sizeof(double), 4, fp);
2885                              if (nRead != 4)
2886                                throw("Error in fread");
2887                              solver_->setColLower(iColumn,elements[0]);
2888                              solver_->setColUpper(iColumn,elements[1]);
2889                            } 
2890                          }
2891                        }
2892                      }
2893                      if (fp)
2894                        fclose(fp);
2895                      char printBuffer[200];
2896                      sprintf(printBuffer,"%d cuts read in of which %d were unique, %d bounds tightened",
2897                             initialCuts,
2898                             globalCuts_.sizeRowCuts()-initialGlobal,nTightened); 
2899                      messageHandler()->message(CBC_GENERAL, messages())
2900                        << printBuffer << CoinMessageEol ;
2901                      delete [] elements;
2902                      delete [] indices;
2903                    }
2904                  }
2905                }
[1330]2906                feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
2907                                         NULL);
[1883]2908                if (multipleRootTries_&&
2909                    (moreSpecialOptions_&134217728)!=0) {
2910                  FILE * fp=NULL;
2911                  size_t nRead;
2912                  int numberColumns=solver_->getNumCols();
2913                  int initialCuts=0;
2914                  if ((moreSpecialOptions_&134217728)!=0) {
2915                    // append so go down to end
2916                    fp = fopen("global.cuts","r+b");
2917                    if (fp) {
2918                      int numCols;
2919                      nRead = fread(&numCols, sizeof(int), 1, fp);
2920                      if (nRead != 1)
2921                        throw("Error in fread");
2922                      if (numberColumns!=numCols) {
2923                        printf("Mismatch on columns %d %d\n",numberColumns,numCols);
2924                        fclose(fp);
2925                        fp=NULL;
2926                      }
2927                    }
2928                  }
2929                  double * elements = new double [numberColumns+2];
2930                  int * indices = new int [numberColumns];
2931                  if (fp) {
2932                    int numberEntries=1;
2933                    while (numberEntries>0) {
2934                      fpos_t position;
2935                      fgetpos(fp, &position);
2936                      nRead = fread(&numberEntries, sizeof(int), 1, fp);
2937                      if (nRead != 1)
2938                        throw("Error in fread");
2939                      if (numberEntries>0) {
2940                        initialCuts++;
2941                        nRead = fread(elements, sizeof(double), numberEntries+2, fp);
2942                        if (nRead != static_cast<size_t>(numberEntries+2))
2943                          throw("Error in fread");
2944                        nRead = fread(indices, sizeof(int), numberEntries, fp);
2945                        if (nRead != static_cast<size_t>(numberEntries))
2946                          throw("Error in fread");
2947                      } else {
2948                        // end
2949                        fsetpos(fp, &position);
2950                      }
2951                    }
2952                  } else {
2953                    fp = fopen("global.cuts","wb");
2954                    size_t nWrite;
2955                    nWrite=fwrite(&numberColumns,sizeof(int),1,fp);
2956                    if (nWrite != 1)
2957                      throw("Error in fwrite");
2958                  }
2959                  size_t nWrite;
2960                  // now append binding cuts
2961                  int numberC=continuousSolver_->getNumRows();
2962                  int numberRows=solver_->getNumRows();
2963                  printf("Saving %d cuts (up from %d)\n",
2964                         initialCuts+numberRows-numberC,initialCuts);
2965                  const double * rowLower = solver_->getRowLower();
2966                  const double * rowUpper = solver_->getRowUpper();
2967                  // Row copy
2968                  CoinPackedMatrix matrixByRow(*solver_->getMatrixByRow());
2969                  const double * elementByRow = matrixByRow.getElements();
2970                  const int * column = matrixByRow.getIndices();
2971                  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
2972                  const int * rowLength = matrixByRow.getVectorLengths();
2973                  for (int iRow=numberC;iRow<numberRows;iRow++) {
2974                    int n=rowLength[iRow];
2975                    assert (n);
2976                    CoinBigIndex start=rowStart[iRow];
2977                    memcpy(elements,elementByRow+start,n*sizeof(double));
2978                    memcpy(indices,column+start,n*sizeof(int));
2979                    elements[n]=rowLower[iRow];
2980                    elements[n+1]=rowUpper[iRow];
2981                    nWrite=fwrite(&n,sizeof(int),1,fp);
2982                    if (nWrite != 1)
2983                      throw("Error in fwrite");
2984                    nWrite=fwrite(elements,sizeof(double),n+2,fp);
2985                    if (nWrite != static_cast<size_t>(n+2))
2986                      throw("Error in fwrite");
2987                    nWrite=fwrite(indices,sizeof(int),n,fp);
2988                    if (nWrite != static_cast<size_t>(n))
2989                      throw("Error in fwrite");
2990                  }
2991                  // eof marker
2992                  int eofMarker=-1;
2993                  nWrite=fwrite(&eofMarker,sizeof(int),1,fp);
2994                  if (nWrite != 1)
2995                    throw("Error in fwrite");
2996                  fclose(fp);
2997                  // do tighter bounds (? later extra to original columns)
2998                  int nTightened=0;
2999                  const double * lower = solver_->getColLower();
3000                  const double * upper = solver_->getColUpper();
3001                  const double * originalLower = continuousSolver_->getColLower();
3002                  const double * originalUpper = continuousSolver_->getColUpper();
3003                  double tolerance=1.0e-5;
3004                  for (int i=0;i<numberColumns;i++) {
3005                    if (lower[i]>originalLower[i]+tolerance) {
3006                      nTightened++;
3007                    }
3008                    if (upper[i]<originalUpper[i]-tolerance) {
3009                      nTightened++;
3010                    }
3011                  }
3012                  if (nTightened) {
3013                    fp = fopen("global.fix","wb");
3014                    size_t nWrite;
3015                    indices[0]=numberColumns;
3016                    if (originalColumns_)
3017                      indices[1]=COIN_INT_MAX;
3018                    else
3019                      indices[1]=-1;
3020                    nWrite=fwrite(indices,sizeof(int),2,fp);
3021                    if (nWrite != 2)
3022                      throw("Error in fwrite");
3023                    for (int i=0;i<numberColumns;i++) {
3024                      int nTightened=0;
3025                      if (lower[i]>originalLower[i]+tolerance) {
3026                        nTightened++;
3027                      }
3028                      if (upper[i]<originalUpper[i]-tolerance) {
3029                        nTightened++;
3030                      }
3031                      if (nTightened) {
3032                        indices[0]=i;
3033                        if (originalColumns_)
3034                          indices[1]=originalColumns_[i];
3035                        elements[0]=lower[i];
3036                        elements[1]=upper[i];
3037                        elements[2]=originalLower[i];
3038                        elements[3]=originalUpper[i];
3039                        nWrite=fwrite(indices,sizeof(int),2,fp);
3040                        if (nWrite != 2)
3041                          throw("Error in fwrite");
3042                        nWrite=fwrite(elements,sizeof(double),4,fp);
3043                        if (nWrite != 4)
3044                          throw("Error in fwrite");
3045                      }
3046                    }
3047                    // eof marker
3048                    indices[0]=-1;
3049                    nWrite=fwrite(indices,sizeof(int),2,fp);
3050                    if (nWrite != 2)
3051                      throw("Error in fwrite");
3052                    fclose(fp);
3053                  }
3054                  delete [] elements;
3055                  delete [] indices;
3056                }
[1330]3057                if ((specialOptions_&524288) != 0 && !parentModel_
3058                        && storedRowCuts_) {
3059                    if (feasible) {
3060                        /* pick up stuff and try again
3061                        add cuts, maybe keep around
3062                        do best solution and if so new heuristics
3063                        obviously tighten bounds
3064                        */
3065                        // A and B probably done on entry
3066                        // A) tight bounds on integer variables
3067                        const double * lower = solver_->getColLower();
3068                        const double * upper = solver_->getColUpper();
3069                        const double * tightLower = storedRowCuts_->tightLower();
3070                        const double * tightUpper = storedRowCuts_->tightUpper();
3071                        int nTightened = 0;
3072                        for (int i = 0; i < numberIntegers_; i++) {
3073                            int iColumn = integerVariable_[i];
3074                            if (tightLower[iColumn] > lower[iColumn]) {
3075                                nTightened++;
3076                                solver_->setColLower(iColumn, tightLower[iColumn]);
[1286]3077                            }
[1330]3078                            if (tightUpper[iColumn] < upper[iColumn]) {
3079                                nTightened++;
3080                                solver_->setColUpper(iColumn, tightUpper[iColumn]);
3081                            }
[1286]3082                        }
[1330]3083                        if (nTightened)
[1641]3084                          COIN_DETAIL_PRINT(printf("%d tightened by alternate cuts\n", nTightened));
[1330]3085                        if (storedRowCuts_->bestObjective() < bestObjective_) {
3086                            // B) best solution
3087                            double objValue = storedRowCuts_->bestObjective();
3088                            setBestSolution(CBC_SOLUTION, objValue,
3089                                            storedRowCuts_->bestSolution()) ;
3090                            // Do heuristics
3091                            // Allow RINS
3092                            for (int i = 0; i < numberHeuristics_; i++) {
3093                                CbcHeuristicRINS * rins
3094                                = dynamic_cast<CbcHeuristicRINS *> (heuristic_[i]);
3095                                if (rins) {
3096                                    rins->setLastNode(-100);
[1286]3097                                }
3098                            }
[1330]3099                            doHeuristicsAtRoot();
[1286]3100                        }
[1393]3101#ifdef JJF_ZERO
[1330]3102                        int nCuts = storedRowCuts_->sizeRowCuts();
3103                        // add to global list
3104                        for (int i = 0; i < nCuts; i++) {
3105                            OsiRowCut newCut(*storedRowCuts_->rowCutPointer(i));
3106                            newCut.setGloballyValidAsInteger(2);
3107                            newCut.mutableRow().setTestForDuplicateIndex(false);
3108                            globalCuts_.insert(newCut) ;
[1286]3109                        }
[1330]3110#else
3111                        addCutGenerator(storedRowCuts_, -99, "Stored from previous run",
3112                                        true, false, false, -200);
[1271]3113#endif
[1330]3114                        // Set cuts as active
3115                        delete [] addedCuts_ ;
3116                        maximumNumberCuts_ = cuts.sizeRowCuts();
3117                        if (maximumNumberCuts_) {
3118                            addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
3119                        } else {
3120                            addedCuts_ = NULL;
[1286]3121                        }
[1330]3122                        for (int i = 0; i < maximumNumberCuts_; i++)
3123                            addedCuts_[i] = new CbcCountRowCut(*cuts.rowCutPtr(i),
3124                                                               NULL, -1, -1, 2);
[1641]3125                        COIN_DETAIL_PRINT(printf("size %d\n", cuts.sizeRowCuts()));
[1330]3126                        cuts = OsiCuts();
3127                        currentNumberCuts_ = maximumNumberCuts_;
3128                        feasible = solveWithCuts(cuts, maximumCutPassesAtRoot_,
3129                                                 NULL);
3130                        for (int i = 0; i < maximumNumberCuts_; i++)
3131                            delete addedCuts_[i];
[1286]3132                    }
[1330]3133                    delete storedRowCuts_;
3134                    storedRowCuts_ = NULL;
[1286]3135                }
[1330]3136            } else {
3137                feasible = false;
[1286]3138            }
[1330]3139        }       else if (solverCharacteristics_->solutionAddsCuts() ||
3140                   solverCharacteristics_->alwaysTryCutsAtRootNode()) {
3141            // may generate cuts and turn the solution
3142            //to an infeasible one
[1764]3143            feasible = solveWithCuts(cuts, 2,
[1330]3144                                     NULL);
[1286]3145        }
[1271]3146    }
[1839]3147    if (rootModels) {
3148      int numberModels = multipleRootTries_%100;
3149      for (int i=0;i<numberModels;i++)
3150        delete rootModels[i];
3151      delete [] rootModels;
3152    }
[1330]3153    // check extra info on feasibility
3154    if (!solverCharacteristics_->mipFeasible())
3155        feasible = false;
[1839]3156    // If max nodes==0 - don't do strong branching
3157    if (!getMaximumNodes()) {
3158      if (feasible)
3159        feasible=false;
3160      else
3161        setMaximumNodes(1); //allow to stop on success
3162    }
3163    topOfTree_=NULL;
[1656]3164#ifdef CLP_RESOLVE
3165    if ((moreSpecialOptions_&2097152)!=0&&!parentModel_&&feasible) {
3166      OsiClpSolverInterface * clpSolver
3167        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3168      if (clpSolver)
3169        resolveClp(clpSolver,0);
3170    }
3171#endif
[1286]3172    // make cut generators less aggressive
3173    for (iCutGenerator = 0; iCutGenerator < numberCutGenerators_; iCutGenerator++) {
3174        CglCutGenerator * generator = generator_[iCutGenerator]->generator();
3175        generator->setAggressiveness(generator->getAggressiveness() - 100);
[1132]3176    }
[1286]3177    currentNumberCuts_ = numberNewCuts_ ;
[1682]3178    if (solverCharacteristics_->solutionAddsCuts()) {
3179      // With some heuristics solver needs a resolve here (don't know if this is bug in heuristics)
3180      solver_->resolve();
3181      if(!isProvenOptimal()){
3182        solver_->initialSolve();
3183      }
3184    }
[1286]3185    // See if can stop on gap
3186    bestPossibleObjective_ = solver_->getObjValue() * solver_->getObjSense();
[1717]3187    if(canStopOnGap()) {
[1286]3188        if (bestPossibleObjective_ < getCutoff())
3189            stoppedOnGap_ = true ;
3190        feasible = false;
[1132]3191    }
[1286]3192    // User event
3193    if (eventHappened_)
3194        feasible = false;
3195#if defined(COIN_HAS_CLP)&&defined(COIN_HAS_CPX)
[1409]3196    /*
3197      This is the notion of using Cbc stuff to get going, then calling cplex to
3198      finish off.
3199    */
[1286]3200    if (feasible && (specialOptions_&16384) != 0 && fastNodeDepth_ == -2 && !parentModel_) {
3201        // Use Cplex to do search!
3202        double time1 = CoinCpuTime();
3203        OsiClpSolverInterface * clpSolver
3204        = dynamic_cast<OsiClpSolverInterface *> (solver_);
3205        OsiCpxSolverInterface cpxSolver;
3206        double direction = clpSolver->getObjSense();
3207        cpxSolver.setObjSense(direction);
3208        // load up cplex
3209        const CoinPackedMatrix * matrix = continuousSolver_->getMatrixByCol();
3210        const double * rowLower = continuousSolver_->getRowLower();
3211        const double * rowUpper = continuousSolver_->getRowUpper();
3212        const double * columnLower = continuousSolver_->getColLower();
3213        const double * columnUpper = continuousSolver_->getColUpper();
3214        const double * objective = continuousSolver_->getObjCoefficients();
3215        cpxSolver.loadProblem(*matrix, columnLower, columnUpper,
3216                              objective, rowLower, rowUpper);
3217        double * setSol = new double [numberIntegers_];
3218        int * setVar = new int [numberIntegers_];
3219        // cplex doesn't know about objective offset
3220        double offset = clpSolver->getModelPtr()->objectiveOffset();
3221        for (int i = 0; i < numberIntegers_; i++) {
3222            int iColumn = integerVariable_[i];
3223            cpxSolver.setInteger(iColumn);
3224            if (bestSolution_) {
3225                setSol[i] = bestSolution_[iColumn];
3226                setVar[i] = iColumn;
3227            }
3228        }
3229        CPXENVptr env = cpxSolver.getEnvironmentPtr();
3230        CPXLPptr lpPtr = cpxSolver.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL);
3231        cpxSolver.switchToMIP();
3232        if (bestSolution_) {
3233            CPXcopymipstart(env, lpPtr, numberIntegers_, setVar, setSol);
3234        }
3235        if (clpSolver->getNumRows() > continuousSolver_->getNumRows() && false) {
3236            // add cuts
3237            const CoinPackedMatrix * matrix = clpSolver->getMatrixByRow();
3238            const double * rhs = clpSolver->getRightHandSide();
3239            const char * rowSense = clpSolver->getRowSense();
3240            const double * elementByRow = matrix->getElements();
3241            const int * column = matrix->getIndices();
3242            const CoinBigIndex * rowStart = matrix->getVectorStarts();
3243            const int * rowLength = matrix->getVectorLengths();
3244            int nStart = continuousSolver_->getNumRows();
3245            int nRows = clpSolver->getNumRows();
3246            int size = rowStart[nRows-1] + rowLength[nRows-1] -
3247                       rowStart[nStart];
3248            int nAdd = 0;
3249            double * rmatval = new double [size];
3250            int * rmatind = new int [size];
3251            int * rmatbeg = new int [nRows-nStart+1];
3252            size = 0;
3253            rmatbeg[0] = 0;
3254            for (int i = nStart; i < nRows; i++) {
3255                for (int k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
3256                    rmatind[size] = column[k];
3257                    rmatval[size++] = elementByRow[k];
3258                }
3259                nAdd++;
3260                rmatbeg[nAdd] = size;
3261            }
3262            CPXaddlazyconstraints(env, lpPtr, nAdd, size,
3263                                  rhs, rowSense, rmatbeg,
3264                                  rmatind, rmatval, NULL);
3265            CPXsetintparam( env, CPX_PARAM_REDUCE,
3266                            // CPX_PREREDUCE_NOPRIMALORDUAL (0)
3267                            CPX_PREREDUCE_PRIMALONLY);
3268        }
3269        if (getCutoff() < 1.0e50) {
3270            double useCutoff = getCutoff() + offset;
3271            if (bestObjective_ < 1.0e50)
3272                useCutoff = bestObjective_ + offset + 1.0e-7;
3273            cpxSolver.setDblParam(OsiDualObjectiveLimit, useCutoff*
3274                                  direction);
3275            if ( direction > 0.0 )
3276                CPXsetdblparam( env, CPX_PARAM_CUTUP, useCutoff ) ; // min
3277            else
3278                CPXsetdblparam( env, CPX_PARAM_CUTLO, useCutoff ) ; // max
3279        }
3280        CPXsetdblparam(env, CPX_PARAM_EPGAP, dblParam_[CbcAllowableFractionGap]);
3281        delete [] setSol;
3282        delete [] setVar;
3283        char printBuffer[200];
3284        if (offset) {
3285            sprintf(printBuffer, "Add %g to all Cplex messages for true objective",
3286                    -offset);
3287            messageHandler()->message(CBC_GENERAL, messages())
3288            << printBuffer << CoinMessageEol ;
3289            cpxSolver.setDblParam(OsiObjOffset, offset);
3290        }
3291        cpxSolver.branchAndBound();
3292        double timeTaken = CoinCpuTime() - time1;
3293        sprintf(printBuffer, "Cplex took %g seconds",
3294                timeTaken);
3295        messageHandler()->message(CBC_GENERAL, messages())
3296        << printBuffer << CoinMessageEol ;
3297        numberExtraNodes_ = CPXgetnodecnt(env, lpPtr);
3298        numberExtraIterations_ = CPXgetmipitcnt(env, lpPtr);
3299        double value = cpxSolver.getObjValue() * direction;
3300        if (cpxSolver.isProvenOptimal() && value <= getCutoff()) {
3301            feasible = true;
3302            clpSolver->setWarmStart(NULL);
3303            // try and do solution
3304            double * newSolution =
3305                CoinCopyOfArray(cpxSolver.getColSolution(),
3306                                getNumCols());
3307            setBestSolution(CBC_STRONGSOL, value, newSolution) ;
3308            delete [] newSolution;
3309        }
3310        feasible = false;
[1132]3311    }
[1286]3312#endif
[1883]3313    if (!parentModel_&&(moreSpecialOptions_&268435456) != 0) {
3314      // try idiotic idea
3315      CbcObject * obj = new CbcIdiotBranch(this);
3316      obj->setPriority(1); // temp
3317      addObjects(1, &obj);
3318      delete obj;
3319    }
3320   
[1409]3321    /*
3322      A hook to use clp to quickly explore some part of the tree.
3323    */
[1286]3324    if (fastNodeDepth_ == 1000 &&/*!parentModel_*/(specialOptions_&2048) == 0) {
3325        fastNodeDepth_ = -1;
3326        CbcObject * obj =
3327            new CbcFollowOn(this);
3328        obj->setPriority(1);
3329        addObjects(1, &obj);
[1883]3330        delete obj;
[1132]3331    }
[1286]3332    int saveNumberSolves = numberSolves_;