Changeset 1943 for trunk


Ignore:
Timestamp:
Jul 21, 2013 5:05:45 AM (6 years ago)
Author:
forrest
Message:

more options, copy statistics structure analysis
start coding of "switch" variables i.e. badly scaled ints or hi/lo
changes to allow more influence on small branch and bound
changes to get correct printout with threads

Location:
trunk/Cbc/src
Files:
24 edited

Legend:

Unmodified
Added
Removed
  • trunk/Cbc/src/CbcEventHandler.cpp

    r1573 r1943  
    102102}
    103103
     104//-------------------------------------------------------------------
     105// event() -- return the action for an event.
     106//-------------------------------------------------------------------
     107
     108CbcEventHandler::CbcAction CbcEventHandler::event(CbcEvent event,
     109                                                  void * /*data*/)
     110/*
     111  If an event/action map exists and contains an entry for the event, return it.
     112  Otherwise return the default action.
     113*/
     114{
     115    if (eaMap_ != 0) {
     116        eaMapPair::iterator entry = eaMap_->find(event) ;
     117        if (entry != eaMap_->end()) {
     118            return (entry->second) ;
     119        } else {
     120            return (dfltAction_) ;
     121        }
     122    } else {
     123        return (dfltAction_) ;
     124    }
     125}
     126
  • trunk/Cbc/src/CbcEventHandler.hpp

    r1632 r1943  
    9999        /*! After failed heuristic. */
    100100        afterHeuristic,
     101        /*! On entry to small branch and bound. */
     102        smallBranchAndBound,
    101103        /*! End of search. */
    102104        endSearch
     
    121123        addCuts,
    122124        /*! Pretend solution never happened. */
    123         killSolution
     125        killSolution,
     126        /*! Take action on modified data. */
     127        takeAction
    124128
    125129    } ;
     
    140144    */
    141145    virtual CbcAction event(CbcEvent whichEvent) ;
     146
     147    /*! \brief Return the action to be taken for an event - and modify data.
     148
     149      Return the action that should be taken in response to the event passed as
     150      the parameter. The default implementation simply reads a return code
     151      from a map.
     152    */
     153    virtual CbcAction event(CbcEvent whichEvent, void * data) ;
    142154
    143155    //@}
  • trunk/Cbc/src/CbcFullNodeInfo.cpp

    r1899 r1943  
    130130{
    131131    OsiSolverInterface *solver = model->solver() ;
    132 
     132    // may be end game
     133    if (!active_)
     134      return;
    133135    // branch - do bounds
    134136    assert (active_ == 7 || active_ == 15);
  • trunk/Cbc/src/CbcHeuristic.cpp

    r1886 r1943  
    2424#include "CbcHeuristic.hpp"
    2525#include "CbcHeuristicFPump.hpp"
     26#include "CbcEventHandler.hpp"
    2627#include "CbcStrategy.hpp"
    2728#include "CglPreProcess.hpp"
     
    148149        numberNodesDone_(0),
    149150        inputSolution_(NULL)
    150 {}
     151{
     152}
    151153
    152154void
     
    407409                    if ((numCouldRun_ % howOften_) == 0 &&
    408410                            numberSolutionsFound_*howOften_ < numCouldRun_) {
     411                      //#define COIN_DEVELOP
    409412#ifdef COIN_DEVELOP
    410413                        int old = howOften_;
     
    664667}
    665668
    666 
     669//static int saveModel=0;
    667670// Do mini branch and bound (return 1 if solution)
    668671int
     
    671674                                  double cutoff, std::string name) const
    672675{
     676  CbcEventHandler *eventHandler = model_->getEventHandler() ;
     677  // Use this fraction
     678  double fractionSmall = fractionSmall_;
     679  int maximumSolutions =  model_->getMaximumSolutions();
     680  int iterationMultiplier = 100;
     681  if (eventHandler) {
     682    typedef struct {
     683      double fractionSmall;
     684      double spareDouble[3];
     685      OsiSolverInterface * solver;
     686      void * sparePointer[2];
     687      int numberNodes;
     688      int maximumSolutions;
     689      int iterationMultiplier;
     690      int howOften;
     691      int spareInt[3];
     692    } SmallMod;
     693    SmallMod temp;
     694    temp.solver=solver;
     695    temp.fractionSmall=fractionSmall;
     696    temp.numberNodes=numberNodes;
     697    temp.iterationMultiplier=iterationMultiplier;
     698    temp.howOften=howOften_;
     699    temp.maximumSolutions=maximumSolutions;
     700    CbcEventHandler::CbcAction status =
     701      eventHandler->event(CbcEventHandler::smallBranchAndBound,
     702                          &temp);
     703    if (status==CbcEventHandler::killSolution)
     704      return -1;
     705    if (status==CbcEventHandler::takeAction) {
     706      fractionSmall=temp.fractionSmall;
     707      numberNodes=temp.numberNodes;
     708      iterationMultiplier=temp.iterationMultiplier;
     709      howOften_=temp.howOften;
     710      maximumSolutions=temp.maximumSolutions;
     711    }
     712  }
     713#if 0
     714  if (saveModel || model_->getMaximumSolutions()==100) {
     715    printf("writing model\n");
     716    solver->writeMpsNative("before.mps", NULL, NULL, 2, 1);
     717  }
     718#endif
    673719    // size before
    674720    int shiftRows = 0;
     
    681727           name.c_str(), solver->getNumRows(), solver->getNumCols());
    682728#endif
    683     // Use this fraction
    684     double fractionSmall = fractionSmall_;
    685729    double before = 2 * numberRowsStart + numberColumnsStart;
    686730    if (before > 40000.0) {
     
    729773    //assert ((saveModelOptions&2048) == 0);
    730774    model_->setSpecialOptions(saveModelOptions | 2048);
    731     {
     775    if (fractionSmall<1.0) {
    732776        int saveLogLevel = solver->messageHandler()->logLevel();
    733777        if (saveLogLevel == 1)
     
    864908        OsiSolverInterface * solver2 = NULL;
    865909        if ((model_->moreSpecialOptions()&65536)!=0)
    866           process.setOptions(2+4+8); // no cuts
     910          process.setOptions(2+4+8+16); // no cuts
     911        else
     912          process.setOptions(16); // no complicated dupcol stuff
    867913        /* Do not try and produce equality cliques and
    868914           do up to 2 passes (normally) 5 if restart */
    869915        int numberPasses = 2;
     916        if ((model_->moreSpecialOptions2()&16)!=0) {
     917          // quick
     918          process.setOptions(2+4+8+16); // no cuts
     919          numberPasses = 1;
     920        }
    870921        if (numberNodes < 0) {
    871922          numberPasses = 5;
     
    9561007                    if ((saveModelOptions&2048) == 0)
    9571008                      model.setMoreSpecialOptions(model_->moreSpecialOptions());
     1009                      model.setMoreSpecialOptions2(model_->moreSpecialOptions2());
    9581010                    // off conflict analysis
    9591011                    model.setMoreSpecialOptions(model.moreSpecialOptions()&~4194304);
     
    10361088#endif
    10371089                model.setParentModel(*model_);
    1038                 model.setMaximumSolutions(model_->getMaximumSolutions());
     1090                model.setMaximumSolutions(maximumSolutions);
    10391091                model.setOriginalColumns(process.originalColumns());
    10401092                model.setSearchStrategy(-1);
     
    12081260                        setCutAndHeuristicOptions(model);
    12091261                        // not too many iterations
    1210                         model.setMaximumNumberIterations(100*(numberNodes + 10));
     1262                        model.setMaximumNumberIterations(iterationMultiplier*(numberNodes + 10));
    12111263                        // Not fast stuff
    12121264                        model.setFastNodeDepth(-1);
     
    12281280                    int saveOptions = model_->specialOptions();
    12291281                    model_->setSpecialOptions(saveOptions|1048576);
     1282                    // and switch off debugger
     1283                    model.setSpecialOptions(model.specialOptions()&(~1));
     1284#if 0 //def COIN_HAS_CLP
     1285                    OsiClpSolverInterface * clpSolver
     1286                      = dynamic_cast<OsiClpSolverInterface *> (model.solver());
     1287                    if (clpSolver)
     1288                      clpSolver->zapDebugger();
     1289#endif
    12301290                    model.branchAndBound();
    12311291                    model_->setHeuristicModel(NULL);
     
    13331393    } else {
    13341394        returnCode = 2; // infeasible finished
     1395        printf("Infeasible on initial solve\n");
    13351396    }
    13361397    model_->setSpecialOptions(saveModelOptions);
     
    18081869    double primalTolerance;
    18091870    solver->getDblParam(OsiPrimalTolerance, primalTolerance);
     1871    double useTolerance = primalTolerance;
    18101872
    18111873    int numberRows = matrix_.getNumRows();
     
    18661928        int iColumn = integerVariable[i];
    18671929        double value = newSolution[iColumn];
     1930        double thisTolerance = integerTolerance;
    18681931        if (fabs(floor(value + 0.5) - value) > integerTolerance) {
    18691932            double below = floor(value);
     
    19351998                            double thisCost = -direction * objective[iColumn] * distance;
    19361999                            if (integerType[iColumn]) {
    1937                                 distance = ceil(distance - primalTolerance);
    1938                                 if (currentValue - distance >= lowerValue - primalTolerance) {
    1939                                     if (absInfeasibility - distance*absElement < -gap - primalTolerance)
     2000                                distance = ceil(distance - useTolerance);
     2001                                if (currentValue - distance >= lowerValue - useTolerance) {
     2002                                    if (absInfeasibility - distance*absElement < -gap - useTolerance)
    19402003                                        thisCost = 1.0e100; // no good
    19412004                                    else
     
    19612024                            if (integerType[iColumn]) {
    19622025                                distance = ceil(distance - 1.0e-7);
    1963                                 assert (currentValue - distance <= upperValue + primalTolerance);
    1964                                 if (absInfeasibility - distance*absElement < -gap - primalTolerance)
     2026                                assert (currentValue - distance <= upperValue + useTolerance);
     2027                                if (absInfeasibility - distance*absElement < -gap - useTolerance)
    19652028                                    thisCost = 1.0e100; // no good
    19662029                                else
     
    21392202        // and now all if improving
    21402203        double lastChange = penaltyChange ? 1.0 : 0.0;
    2141         while (lastChange > 1.0e-2) {
     2204        int numberPasses=0;
     2205        while (lastChange > 1.0e-2 && numberPasses < 1000) {
    21422206            lastChange = 0;
     2207            numberPasses++;
    21432208            for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    21442209                bool isInteger = (integerType[iColumn] != 0);
     
    23552420                    bool good = true;
    23562421                    double newValue = newSolution[iColumn] + move;
    2357                     if (newValue < lower[iColumn] - primalTolerance ||
    2358                             newValue > upper[iColumn] + primalTolerance) {
     2422                    if (newValue < lower[iColumn] - useTolerance ||
     2423                            newValue > upper[iColumn] + useTolerance) {
    23592424                        move = 0.0;
    23602425                    } else {
  • trunk/Cbc/src/CbcHeuristic.hpp

    r1880 r1943  
    332332
    333333    /// How often to do (code can change)
    334     int howOften_;
     334    mutable int howOften_;
    335335    /// How much to increase how often
    336336    double decayFactor_;
  • trunk/Cbc/src/CbcHeuristicDINS.cpp

    r1899 r1943  
    181181                values_[i] = NULL;
    182182        } else {
    183             assert (numberIntegers == numberIntegers_);
     183            assert (numberIntegers >= numberIntegers_);
    184184        }
    185185        // move solutions (0 will be most recent)
     
    193193        }
    194194        int i;
    195         for (i = 0; i < numberIntegers; i++) {
     195        for (i = 0; i < numberIntegers_; i++) {
    196196            int iColumn = integerVariable[i];
    197197            double value = bestSolution[iColumn];
     
    205205        OsiSolverInterface * solver = model_->solver();
    206206
    207         int numberIntegers = model_->numberIntegers();
     207        //int numberIntegers = model_->numberIntegers();
    208208        const int * integerVariable = model_->integerVariable();
    209209
     
    223223            const double * continuousSolution = newSolver->getColSolution();
    224224            // Space for added constraint
    225             double * element = new double [numberIntegers];
    226             int * column = new int [numberIntegers];
     225            double * element = new double [numberIntegers_];
     226            int * column = new int [numberIntegers_];
    227227            int i;
    228228            int nFix = 0;
     
    233233            double bias = localSpace;
    234234            int okSame = numberKeptSolutions_ - 1;
    235             for (i = 0; i < numberIntegers; i++) {
     235            for (i = 0; i < numberIntegers_; i++) {
    236236                int iColumn = integerVariable[i];
    237237                const OsiObject * object = model_->object(i);
     
    333333            << generalPrint
    334334            << CoinMessageEol;
    335             if (nFix > numberIntegers / 10) {
     335            if (nFix > numberIntegers_ / 10) {
    336336#ifdef JJF_ZERO
    337337                newSolver->initialSolve();
    338338                printf("obj %g\n", newSolver->getObjValue());
    339                 for (i = 0; i < numberIntegers; i++) {
     339                for (i = 0; i < numberIntegers_; i++) {
    340340                    int iColumn = integerVariable[i];
    341341                    printf("%d new bounds %g %g - solutions %g %g\n",
  • trunk/Cbc/src/CbcHeuristicFPump.cpp

    r1880 r1943  
    2626#include "CoinTime.hpp"
    2727#include "CbcEventHandler.hpp"
    28 
     28#ifdef SWITCH_VARIABLES
     29#include "CbcSimpleIntegerDynamicPseudoCost.hpp"
     30#endif
    2931
    3032// Default Constructor
     
    829831                            if (returnCode && newSolutionValue < saveValue)
    830832                                numberBandBsolutions++;
     833                        } else if (numberColumns>numberIntegersOrig) {
     834                          // relax continuous
     835                          solver->resolve();
     836                          if (solver->isProvenOptimal()) {
     837                            memcpy(newSolution,solver->getColSolution(),
     838                                   numberColumns);
     839                            newSolutionValue = -saveOffset;
     840                            for (  i = 0 ; i < numberColumns ; i++ )
     841                              newSolutionValue += saveObjective[i] * newSolution[i];
     842                            newSolutionValue *= direction;
     843                            sprintf(pumpPrint, "Relaxing continuous gives %g", newSolutionValue);
     844                          } else {
     845                            sprintf(pumpPrint,"Infeasible when relaxing continuous!\n");
     846                          }
     847                          model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
     848                            << pumpPrint
     849                            << CoinMessageEol;
    831850                        }
    832851                    }
     
    10091028                    double newValue = 0.0;
    10101029                    if (newSolution[iColumn] < lower[iColumn] + primalTolerance) {
    1011                         newValue = costValue + scaleFactor * saveObjective[iColumn];
     1030                      newValue = costValue + scaleFactor * saveObjective[iColumn];
    10121031                    } else {
    1013                         if (newSolution[iColumn] > upper[iColumn] - primalTolerance) {
    1014                             newValue = -costValue + scaleFactor * saveObjective[iColumn];
    1015                         }
     1032                      if (newSolution[iColumn] > upper[iColumn] - primalTolerance) {
     1033                        newValue = -costValue + scaleFactor * saveObjective[iColumn];
     1034                      }
    10161035                    }
    10171036#ifdef RAND_RAND
     
    12361255                    {
    12371256                        const double * newSolution = solver->getColSolution();
    1238                         for (  i = 0 ; i < numberColumns ; i++ ) {
    1239                             if (solver->isInteger(i)) {
    1240                                 double value = newSolution[i];
     1257                        for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++ ) {
     1258                            if (solver->isInteger(iColumn)) {
     1259                                double value = newSolution[iColumn];
    12411260                                double nearest = floor(value + 0.5);
    12421261                                newSumInfeas += fabs(value - nearest);
    1243                                 if (fabs(value - nearest) > 1.0e-6)
     1262                                if (fabs(value - nearest) > 1.0e-6) {
    12441263                                    newNumberInfeas++;
    1245                             }
    1246                             newTrueSolutionValue += saveObjective[i] * newSolution[i];
     1264                                }
     1265                            }
     1266                            newTrueSolutionValue += saveObjective[iColumn] * newSolution[iColumn];
    12471267                        }
    12481268                        newTrueSolutionValue *= direction;
     
    23642384    const double * columnUpper = solver->getColUpper();
    23652385    // Check if valid with current solution (allow for 0.99999999s)
     2386    double newSumInfeas = 0.0;
     2387    int newNumberInfeas = 0;
    23662388    for (i = 0; i < numberIntegers; i++) {
    23672389        int iColumn = integerVariable[i];
    23682390        double value = solution[iColumn];
    23692391        double round = floor(value + 0.5);
    2370         if (fabs(value - round) > primalTolerance)
    2371             break;
     2392        if (fabs(value - round) > primalTolerance) {
     2393          newSumInfeas += fabs(value-round);
     2394          newNumberInfeas++;
     2395        }
    23722396    }
    2373     if (i == numberIntegers) {
     2397    if (!newNumberInfeas) {
    23742398        // may be able to use solution even if 0.99999's
    23752399        double * saveLower = CoinCopyOfArray(columnLower, numberColumns);
  • trunk/Cbc/src/CbcHeuristicLocal.cpp

    r1839 r1943  
    936936        : CbcHeuristic()
    937937{
     938    increment_ = 0.01;
    938939    feasibilityPump_ = NULL;
    939940    numberSolutions_ = 0;
     
    948949        : CbcHeuristic(model)
    949950{
     951    increment_ = 0.01;
    950952    feasibilityPump_ = NULL;
    951953    numberSolutions_ = 0;
     
    987989  numberSolutions_(rhs.numberSolutions_)
    988990{
     991    increment_ = rhs.increment_;
    989992    feasibilityPump_ = NULL;
    990993    if (model_ && rhs.used_) {
     
    10041007    if (this != &rhs) {
    10051008        CbcHeuristic::operator=(rhs);
     1009        increment_ = rhs.increment_;
    10061010        numberSolutions_ = rhs.numberSolutions_;
    10071011        delete [] used_;
     
    10951099  double cutoff=model_->getCutoff();
    10961100  assert (cutoff<1.0e20);
    1097   if (model_->getCutoffIncrement()<1.0e-4)
    1098     cutoff -= 0.01;
     1101  if (model_->getCutoffIncrement()<1.0e-4) {
     1102    cutoff -= increment_;
     1103  }
    10991104  double offset;
    11001105  newSolver->getDblParam(OsiObjOffset, offset);
     
    11851190            numberNodesDone_,numberNodes_,
    11861191            numberIncrease,sumIncrease,numberDecrease,sumDecrease);
     1192    if (!numberIncrease&&!numberDecrease) {
     1193      // somehow tolerances are such that we can slip through
     1194      // change for next time
     1195      increment_ += CoinMax(increment_,fabs(solutionValue+offset)*1.0e-10);
     1196    }
    11871197  } else {
    11881198    sprintf(proxPrint,"Proximity search ran %d nodes - no new solution",
  • trunk/Cbc/src/CbcHeuristicLocal.hpp

    r1802 r1943  
    125125    virtual int solution(double & objectiveValue,
    126126                         double * newSolution);
    127 
     127    /// Set extra increment
     128    inline void setIncrement(double value)
     129    { increment_ = value;}
    128130    /// Used array so we can set
    129131    inline int * used() const {
     
    133135protected:
    134136    // Data
    135     // Copy of Feasibility pump
     137    /// Increment to use if no change
     138    double increment_;
     139    /// Copy of Feasibility pump
    136140    CbcHeuristicFPump * feasibilityPump_;
    137     // Number of solutions so we only do after new solution
     141    /// Number of solutions so we only do after new solution
    138142    int numberSolutions_;
    139143    /// Whether a variable has been in a solution (also when)
  • trunk/Cbc/src/CbcHeuristicRINS.cpp

    r1899 r1943  
    325325#endif
    326326            }
    327             //printf("%d integers have same value\n",nFix);
     327            if (solutionValue==-COIN_DBL_MAX) {
     328              // return fixings in betterSolution
     329              const double * colLower = newSolver->getColLower();
     330              const double * colUpper = newSolver->getColUpper();
     331              for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     332                if (colLower[iColumn]==colUpper[iColumn])
     333                  betterSolution[iColumn]=colLower[iColumn];
     334                else
     335                  betterSolution[iColumn]=COIN_DBL_MAX;
     336              }
     337              delete newSolver;
     338              return 0;
     339            }
     340            //printf("RINS %d integers have same value\n",nFix);
    328341            returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
    329342                                             model_->getCutoff(), "CbcHeuristicRINS");
  • trunk/Cbc/src/CbcMessage.cpp

    r1915 r1943  
    5454    {CBC_OTHER_STATS, 35, 1, "Maximum depth %d, %g variables fixed on reduced cost"},
    5555    {CBC_HEURISTICS_OFF, 36, 1, "Heuristics switched off as %d branching objects are of wrong type"},
    56     {CBC_STATUS2, 37, 1, "%d nodes, %d on tree, best %g - possible %g depth %d unsat %d its %d (%.2f seconds)"},
     56    {CBC_STATUS2, 37, 1, "%d nodes, %d on tree, best %g - possible %g depth %d unsat %d value %g its %d (%.2f seconds)"},
    5757    {CBC_FPUMP1, 38, 1, "%s"},
    5858    {CBC_FPUMP2, 39, 2, "%s"},
  • trunk/Cbc/src/CbcMipStartIO.cpp

    r1899 r1943  
    1111#include <OsiSolverInterface.hpp>
    1212#include "CbcMessage.hpp"
     13#include "CbcHeuristic.hpp"
    1314#include <CbcModel.hpp>
    1415#include "CbcMipStartIO.hpp"
     
    2324
    2425   for ( size_t i=0 ; i<l ; ++i )
    25       if (!(isdigit(str[i])||(str[i]=='.')))
     26     if (!(isdigit(str[i])||(str[i]=='.')||(str[i]=='-')||(str[i]=='e')))
    2627         return false;
    2728
     
    113114   int notFound = 0;
    114115   char colNotFound[256] = "";
     116   int nContinuousFixed = 0;
     117#ifndef JUST_FIX_INTEGER
     118#define JUST_FIX_INTEGER 0
     119#endif
     120#if JUST_FIX_INTEGER > 1
     121   // all not mentioned are at zero
     122   for ( int i=0 ; (i<lp->getNumCols()) ; ++i )
     123     {
     124       if (lp->isInteger(i))
     125         lp->setColBounds( i, 0.0, 0.0 );
     126     }
     127#endif
    115128   for ( int i=0 ; (i<(int)colValues.size()) ; ++i )
    116129   {
     
    126139         const int idx = mIt->second;
    127140         double v = colValues[i].second;
     141#if JUST_FIX_INTEGER
     142         if (!lp->isInteger(idx))
     143           continue;
     144#endif
    128145         if (v<1e-8)
    129146            v = 0.0;
    130147         if (lp->isInteger(idx))  // just to avoid small
    131148            v = floor( v+0.5 );   // fractional garbage
     149         else
     150           nContinuousFixed++;
    132151         lp->setColBounds( idx, v, v );
    133152         ++fixed;
     
    148167        << printLine << CoinMessageEol;
    149168   }
    150 
     169#if JUST_FIX_INTEGER
     170   lp->setHintParam(OsiDoPresolveInInitial, true, OsiHintDo) ;
     171#endif
    151172   lp->initialSolve();
    152173   if (!lp->isProvenOptimal())
     
    154175      model->messageHandler()->message(CBC_GENERAL, model->messages())
    155176        << "Warning: mipstart values could not be used to build a solution." << CoinMessageEol;
    156       status = 1;
    157       goto TERMINATE;
     177      if (nContinuousFixed) {
     178        model->messageHandler()->message(CBC_GENERAL, model->messages())
     179          << "Trying just fixing integer variables." << CoinMessageEol;
     180        int numberColumns = lp->getNumCols();
     181        const double * oldLower = model->solver()->getColLower();
     182        const double * oldUpper = model->solver()->getColUpper();
     183        for ( int i=0 ; i<numberColumns ; ++i ) {
     184          if (!lp->isInteger(i)) {
     185            lp->setColLower(i,oldLower[i]);
     186            lp->setColUpper(i,oldUpper[i]);
     187          }
     188        }
     189        lp->initialSolve();
     190        if (!lp->isProvenOptimal())
     191          model->messageHandler()->message(CBC_GENERAL, model->messages())
     192            << "Still no good." << CoinMessageEol;
     193      }
     194      if (!lp->isProvenOptimal()) {
     195        status = 1;
     196        goto TERMINATE;
     197      }
    158198   }
    159199
     
    165205        << printLine << CoinMessageEol;
    166206      double start = CoinCpuTime();
     207#if 1
     208      CbcSerendipity heuristic(*model);
     209      heuristic.setFractionSmall(2.0);
     210      heuristic.setFeasibilityPumpOptions(1008013);
     211      int returnCode = heuristic.smallBranchAndBound(lp,
     212                                                     1000, sol,
     213                                                     compObj,
     214                                                     model->getCutoff(),
     215                                                     "ReduceInMIPStart");
     216      if ((returnCode&1) != 0) {
     217         sprintf( printLine,"Mini branch and bound defined values for remaining variables in %.2f seconds.",
     218                  CoinCpuTime()-start);
     219         model->messageHandler()->message(CBC_GENERAL, model->messages())
     220           << printLine << CoinMessageEol;
     221         foundIntegerSol = true;
     222         obj = compObj;
     223      }
     224#else
    167225      CbcModel babModel( *lp );
    168226      babModel.setLogLevel( 0 );
     
    180238         obj = compObj = babModel.getObjValue();
    181239      }
     240#endif
    182241      else
    183242      {
     
    200259      model->messageHandler()->message(CBC_GENERAL, model->messages())
    201260        << printLine << CoinMessageEol;
     261      {
     262        int numberColumns=lp->getNumCols();
     263        double largestInfeasibility = 0.0;
     264        double primalTolerance ;
     265        double offset;
     266        lp->getDblParam(OsiObjOffset, offset);
     267        lp->getDblParam(OsiPrimalTolerance, primalTolerance) ;
     268        const double *objective = lp->getObjCoefficients() ;
     269        const double * rowLower = lp->getRowLower() ;
     270        const double * rowUpper = lp->getRowUpper() ;
     271        const double * columnLower = lp->getColLower() ;
     272        const double * columnUpper = lp->getColUpper() ;
     273        int numberRows = lp->getNumRows() ;
     274        double *rowActivity = new double[numberRows] ;
     275        memset(rowActivity, 0, numberRows*sizeof(double)) ;
     276        double *rowSum = new double[numberRows] ;
     277        memset(rowSum, 0, numberRows*sizeof(double)) ;
     278        const double * element = lp->getMatrixByCol()->getElements();
     279        const int * row = lp->getMatrixByCol()->getIndices();
     280        const CoinBigIndex * columnStart = lp->getMatrixByCol()->getVectorStarts();
     281        const int * columnLength = lp->getMatrixByCol()->getVectorLengths();
     282        const CoinPackedMatrix * rowCopy = lp->getMatrixByRow();
     283        const int * column = rowCopy->getIndices();
     284        const int * rowLength = rowCopy->getVectorLengths();
     285        const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     286        const double * elementByRow = rowCopy->getElements();
     287        double objValue=-offset;
     288        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     289          double value = sol[iColumn];
     290          if (lp->isInteger(iColumn))
     291            assert (fabs(value-floor(value+0.5))<1.0e-6);
     292          objValue += value*objective[iColumn];
     293          if (value>columnUpper[iColumn]) {
     294            if (value-columnUpper[iColumn]>1.0e-8)
     295              printf("column %d has value %.12g above %.12g\n",iColumn,value,columnUpper[iColumn]);
     296            value=columnUpper[iColumn];
     297          } else if (value<columnLower[iColumn]) {
     298            if (value-columnLower[iColumn]<-1.0e-8)
     299              printf("column %d has value %.12g below %.12g\n",iColumn,value,columnLower[iColumn]);
     300            value=columnLower[iColumn];
     301          }
     302          if (value) {
     303            CoinBigIndex start = columnStart[iColumn];
     304            CoinBigIndex end = start + columnLength[iColumn];
     305            for (CoinBigIndex j = start; j < end; j++) {
     306              int iRow = row[j];
     307              if (fabs(value)<1.0e-6&&fabs(value*element[j])>1.0e-5)
     308                printf("Column %d row %d value %.8g element %g %s\n",
     309                       iColumn,iRow,value,element[j],lp->isInteger(iColumn) ? "integer" : "");
     310              rowActivity[iRow] += value * element[j];
     311              rowSum[iRow] += fabs(value * element[j]);
     312            }
     313          }
     314        }
     315        for (int i = 0 ; i < numberRows ; i++) {
     316#if 0 //def CLP_INVESTIGATE
     317          double inf;
     318          inf = rowLower[i] - rowActivity[i];
     319          if (inf > primalTolerance)
     320            printf("Row %d inf %g sum %g %g <= %g <= %g\n",
     321                   i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
     322          inf = rowActivity[i] - rowUpper[i];
     323          if (inf > primalTolerance)
     324            printf("Row %d inf %g sum %g %g <= %g <= %g\n",
     325                   i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
     326#endif
     327          double infeasibility = CoinMax(rowActivity[i]-rowUpper[i],
     328                                         rowLower[i]-rowActivity[i]);
     329          // but allow for errors
     330          double factor = CoinMax(1.0,rowSum[i]*1.0e-3);
     331          if (infeasibility>largestInfeasibility*factor) {
     332            largestInfeasibility = infeasibility/factor;
     333            printf("Cinf of %g on row %d sum %g scaled %g\n",
     334                   infeasibility,i,rowSum[i],largestInfeasibility);
     335            if (infeasibility>1.0e10) {
     336              for (CoinBigIndex j=rowStart[i];
     337                   j<rowStart[i]+rowLength[i];j++) {
     338                printf("col %d element %g\n",
     339                       column[j],elementByRow[j]);
     340              }
     341            }
     342          }
     343        }
     344        delete [] rowActivity ;
     345        delete [] rowSum;
     346        if (largestInfeasibility > 10.0*primalTolerance)
     347          printf("Clargest infeasibility is %g - obj %g\n", largestInfeasibility,objValue);
     348        else
     349          printf("Cfeasible (%g) - obj %g\n", largestInfeasibility,objValue);
     350      }
    202351      for ( int i=0 ; (i<lp->getNumCols()) ; ++i )
    203352      {
     353#if 0
    204354         if (sol[i]<1e-8)
    205355            sol[i] = 0.0;
     
    207357            if (lp->isInteger(i))
    208358               sol[i] = floor( sol[i]+0.5 );
    209       }
     359#else
     360         if (lp->isInteger(i)) {
     361           if (fabs(sol[i] - floor( sol[i]+0.5 ))>1.0e-8)
     362             printf("bad sol for %d - %.12g\n",i,sol[i]);
     363           sol[i] = floor( sol[i]+0.5 );
     364         }
     365#endif
     366      }
     367      {
     368        int numberColumns=lp->getNumCols();
     369        double largestInfeasibility = 0.0;
     370        double primalTolerance ;
     371        double offset;
     372        lp->getDblParam(OsiObjOffset, offset);
     373        lp->getDblParam(OsiPrimalTolerance, primalTolerance) ;
     374        const double *objective = lp->getObjCoefficients() ;
     375        const double * rowLower = lp->getRowLower() ;
     376        const double * rowUpper = lp->getRowUpper() ;
     377        const double * columnLower = lp->getColLower() ;
     378        const double * columnUpper = lp->getColUpper() ;
     379        int numberRows = lp->getNumRows() ;
     380        double *rowActivity = new double[numberRows] ;
     381        memset(rowActivity, 0, numberRows*sizeof(double)) ;
     382        double *rowSum = new double[numberRows] ;
     383        memset(rowSum, 0, numberRows*sizeof(double)) ;
     384        const double * element = lp->getMatrixByCol()->getElements();
     385        const int * row = lp->getMatrixByCol()->getIndices();
     386        const CoinBigIndex * columnStart = lp->getMatrixByCol()->getVectorStarts();
     387        const int * columnLength = lp->getMatrixByCol()->getVectorLengths();
     388        const CoinPackedMatrix * rowCopy = lp->getMatrixByRow();
     389        const int * column = rowCopy->getIndices();
     390        const int * rowLength = rowCopy->getVectorLengths();
     391        const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     392        const double * elementByRow = rowCopy->getElements();
     393        double objValue=-offset;
     394        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     395          double value = sol[iColumn];
     396          if (lp->isInteger(iColumn))
     397            assert (fabs(value-floor(value+0.5))<1.0e-6);
     398          objValue += value*objective[iColumn];
     399          if (value>columnUpper[iColumn]) {
     400            if (value-columnUpper[iColumn]>1.0e-8)
     401              printf("column %d has value %.12g above %.12g\n",iColumn,value,columnUpper[iColumn]);
     402            value=columnUpper[iColumn];
     403          } else if (value<columnLower[iColumn]) {
     404            if (value-columnLower[iColumn]<-1.0e-8)
     405              printf("column %d has value %.12g below %.12g\n",iColumn,value,columnLower[iColumn]);
     406            value=columnLower[iColumn];
     407          }
     408          if (value) {
     409            CoinBigIndex start = columnStart[iColumn];
     410            CoinBigIndex end = start + columnLength[iColumn];
     411            for (CoinBigIndex j = start; j < end; j++) {
     412              int iRow = row[j];
     413              rowActivity[iRow] += value * element[j];
     414              rowSum[iRow] += fabs(value * element[j]);
     415            }
     416          }
     417        }
     418        for (int i = 0 ; i < numberRows ; i++) {
     419#if 0 //def CLP_INVESTIGATE
     420          double inf;
     421          inf = rowLower[i] - rowActivity[i];
     422          if (inf > primalTolerance)
     423            printf("Row %d inf %g sum %g %g <= %g <= %g\n",
     424                   i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
     425          inf = rowActivity[i] - rowUpper[i];
     426          if (inf > primalTolerance)
     427            printf("Row %d inf %g sum %g %g <= %g <= %g\n",
     428                   i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
     429#endif
     430          double infeasibility = CoinMax(rowActivity[i]-rowUpper[i],
     431                                         rowLower[i]-rowActivity[i]);
     432          // but allow for errors
     433          double factor = CoinMax(1.0,rowSum[i]*1.0e-3);
     434          if (infeasibility>largestInfeasibility*factor) {
     435            largestInfeasibility = infeasibility/factor;
     436            printf("Dinf of %g on row %d sum %g scaled %g\n",
     437                   infeasibility,i,rowSum[i],largestInfeasibility);
     438            if (infeasibility>1.0e10) {
     439              for (CoinBigIndex j=rowStart[i];
     440                   j<rowStart[i]+rowLength[i];j++) {
     441                printf("col %d element %g\n",
     442                       column[j],elementByRow[j]);
     443              }
     444            }
     445          }
     446        }
     447        delete [] rowActivity ;
     448        delete [] rowSum;
     449        if (largestInfeasibility > 10.0*primalTolerance)
     450          printf("Dlargest infeasibility is %g - obj %g\n", largestInfeasibility,objValue);
     451        else
     452          printf("Dfeasible (%g) - obj %g\n", largestInfeasibility,objValue);
     453      }
     454#if JUST_FIX_INTEGER
     455      const double * oldLower = model->solver()->getColLower();
     456      const double * oldUpper = model->solver()->getColUpper();
     457      const double * dj = lp->getReducedCost();
     458      int nNaturalLB=0;
     459      int nMaybeLB=0;
     460      int nForcedLB=0;
     461      int nNaturalUB=0;
     462      int nMaybeUB=0;
     463      int nForcedUB=0;
     464      int nOther=0;
     465      for ( int i=0 ; i<lp->getNumCols() ; ++i ) {
     466        if (lp->isInteger(i)) {
     467          if (sol[i]==oldLower[i]) {
     468            if (dj[i]>1.0e-5)
     469              nNaturalLB++;
     470            else if (dj[i]<-1.0e-5)
     471              nForcedLB++;
     472            else
     473              nMaybeLB++;
     474          } else if (sol[i]==oldUpper[i]) {
     475            if (dj[i]<-1.0e-5)
     476              nNaturalUB++;
     477            else if (dj[i]>1.0e-5)
     478              nForcedUB++;
     479            else
     480              nMaybeUB++;
     481          } else {
     482            nOther++;
     483          }
     484        }
     485      }
     486      printf("%d other, LB %d natural, %d neutral, %d forced, UB %d natural, %d neutral, %d forced\n",
     487             nOther,nNaturalLB,nMaybeLB,nForcedLB,
     488             nNaturalUB,nMaybeUB,nForcedUB=0);
     489#endif
    210490   }
    211491
  • trunk/Cbc/src/CbcModel.cpp

    r1924 r1943  
    24532453        }
    24542454    }
     2455#ifdef SWITCH_VARIABLES
     2456    // see if any switching variables
     2457    if (numberIntegers_<solver_->getNumCols())
     2458      findSwitching();
     2459#endif
    24552460    /*
    24562461      Run heuristics at the root. This is the only opportunity to run FPump; it
     
    34063411    numberFixedAtRoot_ = 0;
    34073412    numberFixedNow_ = 0;
     3413    if (!parentModel_&&(moreSpecialOptions2_&2)!=0) {
     3414#ifdef COIN_HAS_CLP
     3415      OsiClpSolverInterface * clpSolver
     3416        = dynamic_cast<OsiClpSolverInterface *> (solver_);
     3417      if (clpSolver) {
     3418        if (getCutoff()>1.0e20) {
     3419          printf("Zapping costs\n");
     3420          int numberColumns=solver_->getNumCols();
     3421          double * zeroCost = new double [numberColumns];
     3422          // could make small random
     3423          memset(zeroCost,0,numberColumns*sizeof(double));
     3424          solver_->setObjective(zeroCost);
     3425          double objValue = solver_->getObjValue();
     3426          solver_->setDblParam(OsiObjOffset,-objValue);
     3427          clpSolver->getModelPtr()->setObjectiveValue(objValue);
     3428          delete [] zeroCost;
     3429        } else {
     3430          moreSpecialOptions2_ &= ~2;
     3431        }
     3432      } else {
     3433#endif
     3434          moreSpecialOptions2_ &= ~2;
     3435#ifdef COIN_HAS_CLP
     3436      }
     3437#endif
     3438    }
    34083439    int numberIterationsAtContinuous = numberIterations_;
    34093440    //solverCharacteristics_->setSolver(solver_);
     
    38113842    assert (!newNode || newNode->objectiveValue() <= cutoff) ;
    38123843    // Save address of root node as we don't want to delete it
    3813     // initialize for print out
    3814     int lastDepth = 0;
    3815     int lastUnsatisfied = 0;
    3816     if (newNode)
    3817         lastUnsatisfied = newNode->numberUnsatisfied();
    38183844    /*
    38193845      The common case is that the lp relaxation is feasible but doesn't satisfy
     
    43564382                messageHandler()->message(CBC_STATUS2, messages())
    43574383                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
    4358                 << lastDepth << lastUnsatisfied << numberIterations_
     4384                << tree_->lastDepth() << tree_->lastUnsatisfied()
     4385                << tree_->lastObjective() << numberIterations_
    43594386                << getCurrentSeconds()
    43604387                << CoinMessageEol ;
     
    43624389                messageHandler()->message(CBC_STATUS2, messages())
    43634390                << numberNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
    4364                 << lastDepth << lastUnsatisfied << numberIterations_
     4391                << tree_->lastDepth() << tree_->lastUnsatisfied() << numberIterations_
    43654392                << getCurrentSeconds()
    43664393                << CoinMessageEol ;
     
    43684395                messageHandler()->message(CBC_STATUS3, messages())
    43694396                << numberNodes_ << numberExtraNodes_ << nNodes << bestObjective_ << bestPossibleObjective_
    4370                 << lastDepth << lastUnsatisfied << numberIterations_ << numberExtraIterations_
     4397                << tree_->lastDepth() << tree_->lastUnsatisfied() << numberIterations_ << numberExtraIterations_
    43714398                << getCurrentSeconds()
    43724399                << CoinMessageEol ;
     
    44464473        } else {
    44474474            // Deterministic parallel
    4448             if (tree_->size() < CoinMax(numberThreads_, 8) && !goneParallel) {
     4475          if ((tree_->size() < CoinMax(numberThreads_, 8)||
     4476               hotstartSolution_) && !goneParallel) {
    44494477                node = tree_->bestNode(cutoff) ;
    44504478                // Possible one on tree worse than cutoff
     
    51015129        specialOptions_(0),
    51025130        moreSpecialOptions_(0),
     5131        moreSpecialOptions2_(0),
    51035132        topOfTree_(NULL),
    51045133        subTreeModel_(NULL),
     
    52655294        specialOptions_(0),
    52665295        moreSpecialOptions_(0),
     5296        moreSpecialOptions2_(0),
    52675297        topOfTree_(NULL),
    52685298        subTreeModel_(NULL),
     
    55535583        specialOptions_(rhs.specialOptions_),
    55545584        moreSpecialOptions_(rhs.moreSpecialOptions_),
     5585        moreSpecialOptions2_(rhs.moreSpecialOptions2_),
    55555586        topOfTree_(NULL),
    55565587        subTreeModel_(rhs.subTreeModel_),
     
    59075938        specialOptions_ = rhs.specialOptions_;
    59085939        moreSpecialOptions_ = rhs.moreSpecialOptions_;
     5940        moreSpecialOptions2_ = rhs.moreSpecialOptions2_;
    59095941        subTreeModel_ = rhs.subTreeModel_;
    59105942        heuristicModel_ = NULL;
     
    63636395    specialOptions_ = rhs.specialOptions_;
    63646396    moreSpecialOptions_ = rhs.moreSpecialOptions_;
     6397    moreSpecialOptions2_ = rhs.moreSpecialOptions2_;
    63656398    numberStrong_ = rhs.numberStrong_;
    63666399    numberBeforeTrust_ = rhs.numberBeforeTrust_;
     
    72607293    }
    72617294    numberDJFixed_ += numberFixed - numberTightened;
     7295#ifdef SWITCH_VARIABLES
     7296    if (numberFixed)
     7297      fixAssociated(NULL,0);
     7298#endif
    72627299    return numberFixed;
    72637300}
     
    1099411031            // treat as if will cost what it says up
    1099511032            double upCost = costValue;
     11033#ifndef BRANCH_BREAKEVEN
     11034#define BRANCH_BREAKEVEN 0.3
     11035#else
     11036            preferredWay=1;
     11037#endif
    1099611038            // and balance at breakeven of 0.3
    10997             double downCost = (0.7 * upCost) / 0.3;
     11039            double downCost = ((1.0-BRANCH_BREAKEVEN) * upCost) / BRANCH_BREAKEVEN;
    1099811040            if (obj1a) {
    1099911041                upCost = obj1a->upPseudoCost();
     
    1103111073        branchingMethod_ = new CbcBranchDynamicDecision();
    1103211074    }
     11075#ifdef SWITCH_VARIABLES
     11076    // see if any switching variables
     11077    if (numberIntegers_<solver_->getNumCols())
     11078      findSwitching();
     11079#endif
    1103311080    synchronizeNumberBeforeTrust();
    1103411081}
     11082#ifdef SWITCH_VARIABLES
     11083// Convert Dynamic to Switching
     11084int
     11085CbcModel::findSwitching()
     11086{
     11087  if ((moreSpecialOptions2_&1)==0)
     11088    return 0;
     11089  const CoinPackedMatrix * rowCopy = solver_->getMatrixByRow();
     11090  const int * column = rowCopy->getIndices();
     11091  const int * rowLength = rowCopy->getVectorLengths();
     11092  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     11093  const double * rowLower = solver_->getRowLower();
     11094  const double * rowUpper = solver_->getRowUpper();
     11095  const double * columnLower = solver_->getColLower();
     11096  const double * columnUpper = solver_->getColUpper();
     11097  const double * element = rowCopy->getElements();
     11098  //const double * element = solver_->getMatrixByCol()->getElements();
     11099  const int * row = solver_->getMatrixByCol()->getIndices();
     11100  const CoinBigIndex * columnStart = solver_->getMatrixByCol()->getVectorStarts();
     11101  const int * columnLength = solver_->getMatrixByCol()->getVectorLengths();
     11102  int numberRows = solver_->getNumRows();
     11103  int numberColumns = solver_->getNumCols();
     11104  int * sort = new int[2*numberRows+2+numberColumns];
     11105  int * whichRow = sort+numberRows+1;
     11106  int * marked = whichRow+numberRows+1;
     11107  memset(marked,0,numberColumns*sizeof(int));
     11108  int nnSwitch=0;
     11109  int nnSwitchTotal=0;
     11110  int n2Switch=0;
     11111  double largeRatio1=1000.0;
     11112  double largeRatio2=100.0;
     11113  for (int i=0;i<numberIntegers_;i++) {
     11114    int iColumn = integerVariable_[i];
     11115    if (columnLower[iColumn]||columnUpper[iColumn]!=1.0)
     11116      continue;
     11117    if (!dynamic_cast <CbcSimpleInteger *> (object_[i]))
     11118      continue;
     11119    int nAdd=0;
     11120    bool takeThis=false;
     11121    CoinBigIndex start = columnStart[iColumn];
     11122    CoinBigIndex end = start + columnLength[iColumn];
     11123    for (CoinBigIndex j = start; j < end; j++) {
     11124      int iRow = row[j];
     11125      if (rowLength[iRow]!=2) {
     11126        continue;
     11127      }
     11128      // for now just 0.0 in rhs
     11129      if (!rowLower[iRow]) {
     11130        if (rowUpper[iRow]!=COIN_DBL_MAX)
     11131          continue;
     11132      } else if (rowLower[iRow]!=-COIN_DBL_MAX) {
     11133        continue;
     11134      } else if (rowUpper[iRow]) {
     11135        continue;
     11136      }
     11137      CoinBigIndex k = rowStart[iRow];
     11138      double bValue, cValue;
     11139      int cColumn;
     11140      if (column[k]==iColumn) {
     11141        bValue=element[k];
     11142        cValue=element[k+1];
     11143        cColumn=column[k+1];
     11144      } else {
     11145        bValue=element[k+1];
     11146        cValue=element[k];
     11147        cColumn=column[k];
     11148      }
     11149      if (solver_->isInteger(cColumn))
     11150        continue;
     11151      if (columnLower[cColumn]<0.0)
     11152        continue;
     11153      if (bValue*cValue>0.0)
     11154        continue;
     11155      if (fabs(bValue)>largeRatio1*fabs(cValue))
     11156        takeThis=true;
     11157      // add to list
     11158      whichRow[nAdd]=iRow;
     11159      sort[nAdd++]=cColumn;
     11160    }
     11161    if (nAdd) {
     11162      n2Switch++;
     11163      CoinSort_2(sort,sort+nAdd,whichRow);
     11164      int last=sort[0];
     11165      for (int k=1;k<nAdd;k++) {
     11166        if (sort[k]==last)
     11167          takeThis=true;
     11168        else
     11169          last=sort[k];
     11170      }
     11171      if (takeThis) {
     11172        int last=sort[0];
     11173        marked[last]++;
     11174        for (int k=1;k<nAdd;k++) {
     11175          if (sort[k]!=last) {
     11176            last=sort[k];
     11177            marked[last]++;
     11178          }
     11179        }
     11180        //printf("Column %d has %d other columns\n",iColumn,nAdd);
     11181        sort[nAdd]=COIN_INT_MAX;
     11182        whichRow[nAdd]=COIN_INT_MAX;
     11183        CbcSimpleIntegerDynamicPseudoCost * thisOne =
     11184          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *> (object_[i]);
     11185        if (thisOne) {
     11186          assert(iColumn == thisOne->columnNumber());
     11187          object_[i]=new CbcSwitchingBinary(thisOne,nAdd,sort,whichRow);
     11188          delete thisOne;
     11189        } else {
     11190          CbcSimpleInteger * thisOne =
     11191            dynamic_cast <CbcSimpleInteger *> (object_[i]);
     11192          assert (thisOne);
     11193          assert(iColumn == thisOne->columnNumber());
     11194          CbcSimpleIntegerDynamicPseudoCost tempObj(this,iColumn,0.1);
     11195          object_[i]=new CbcSwitchingBinary(&tempObj,nAdd,sort,whichRow);
     11196          delete thisOne;
     11197        }
     11198      }
     11199    }
     11200    // see if there is an interesting row
     11201    for (CoinBigIndex j = start; j < end; j++) {
     11202      int iRow = row[j];
     11203      // for now just 0.0 in rhs
     11204      if (!rowLower[iRow]) {
     11205        if (rowUpper[iRow]!=COIN_DBL_MAX)
     11206          continue;
     11207      } else if (rowLower[iRow]!=-COIN_DBL_MAX) {
     11208        continue;
     11209      } else if (rowUpper[iRow]) {
     11210        continue;
     11211      }
     11212      int nOther=0;
     11213      double bEl=0.0;
     11214      double cMax=-COIN_DBL_MAX;
     11215      double cMin=COIN_DBL_MAX;
     11216      for (CoinBigIndex k = rowStart[iRow];
     11217           k<rowStart[iRow]+rowLength[iRow];k++) {
     11218        int jColumn = column[k];
     11219        if (jColumn==iColumn) {
     11220          bEl=element[k];
     11221        } else {
     11222          sort[nOther++]=jColumn;
     11223          if (solver_->isInteger(jColumn)) {
     11224            cMin=-1.0;
     11225            cMax=1.0;
     11226            break;
     11227          } else {
     11228            cMax=CoinMax(cMax,element[k]);
     11229            cMin=CoinMin(cMin,element[k]);
     11230            if (columnLower[jColumn]<0.0) {
     11231              cMin=-1.0;
     11232              cMax=1.0;
     11233              break;
     11234            }
     11235          }
     11236        }
     11237      }
     11238      double largestC = CoinMax(fabs(cMin),fabs(cMax));
     11239      if (((cMin>0.0&&bEl<0.0&&!rowUpper[iRow])||
     11240           (cMin<0.0&&bEl>0.0&&!rowLower[iRow]))&&cMin*cMax>0.0&&
     11241          fabs(bEl)>largeRatio2*largestC) {
     11242        // forces to zero
     11243        CbcSwitchingBinary * object =
     11244          dynamic_cast <CbcSwitchingBinary *> (object_[i]);
     11245        if (!object) {
     11246          // create empty one
     11247          CbcSimpleIntegerDynamicPseudoCost * thisOne =
     11248            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *> (object_[i]);
     11249          if (thisOne) {
     11250            assert(iColumn == thisOne->columnNumber());
     11251            object=new CbcSwitchingBinary(thisOne,0,sort,whichRow);
     11252            delete thisOne;
     11253          } else {
     11254            CbcSimpleInteger * thisOne =
     11255              dynamic_cast <CbcSimpleInteger *> (object_[i]);
     11256            assert (thisOne);
     11257            assert(iColumn == thisOne->columnNumber());
     11258            CbcSimpleIntegerDynamicPseudoCost tempObj(this,iColumn,0.1);
     11259            object=new CbcSwitchingBinary(&tempObj,0,sort,whichRow);
     11260            delete thisOne;
     11261          }
     11262          object_[i]=object;
     11263        }
     11264        object->addZeroSwitches(nOther,sort);
     11265        nnSwitch++;
     11266        nnSwitchTotal+=nOther;
     11267      }
     11268    }
     11269  }
     11270  if (n2Switch+nnSwitch) {
     11271    if (handler_->logLevel()>2)
     11272      printf("%d two switch variables - %d multi (total multi %d)\n",
     11273             n2Switch,nnSwitch,nnSwitchTotal);
     11274    memset(whichRow,0,(numberRows+1)*sizeof(int));
     11275    for (int i=0;i<numberColumns;i++) {
     11276      whichRow[marked[i]]++;
     11277    }
     11278    if (handler_->logLevel()>2) {
     11279      for (int i=0;i<numberRows+1;i++) {
     11280        if (whichRow[i])
     11281          printf("%d variables have %d switches\n",whichRow[i],i);
     11282      }
     11283    }
     11284  }
     11285  delete [] sort;
     11286  // say switches exist
     11287  if (n2Switch+nnSwitch)
     11288    moreSpecialOptions2_|=4;
     11289  return n2Switch+nnSwitch;
     11290}
     11291// Fix associated variables
     11292int
     11293CbcModel::fixAssociated(OsiSolverInterface * solver,int cleanBasis)
     11294{
     11295  int nChanged=0;
     11296  if ((moreSpecialOptions2_&4)!=0) {
     11297    int n=-1;
     11298    while (n) {
     11299      n=0;
     11300      for (int i=0;i<numberObjects_;i++) {
     11301        CbcSwitchingBinary * object = dynamic_cast<CbcSwitchingBinary *> (object_[i]);
     11302        if (object) {
     11303          n += object->setAssociatedBounds(solver,cleanBasis);
     11304        }
     11305      }
     11306      nChanged+=n;
     11307    }
     11308  }
     11309  return nChanged;
     11310}
     11311/* Debug associated variables
     11312   printLevel - 1 summary if bad on fixed
     11313                2 summary if bad on satisfied
     11314                3 for individuals
     11315 */
     11316int
     11317CbcModel::checkAssociated(const OsiSolverInterface * solver,
     11318                          const double * solution,int printLevel)
     11319{
     11320  int nBad=0;
     11321  int nBadFixed=0;
     11322  if ((moreSpecialOptions2_&4)!=0) {
     11323    int nAt0=0;
     11324    int nAt1=0;
     11325    int nBetween=0;
     11326    for (int i=0;i<numberObjects_;i++) {
     11327      CbcSwitchingBinary * object = dynamic_cast<CbcSwitchingBinary *> (object_[i]);
     11328      if (object) {
     11329        int state[3];
     11330        nBad += object->checkAssociatedBounds(solver,solution,printLevel,state,
     11331                                              nBadFixed);
     11332        if (state[0]==0)
     11333          nBetween++;
     11334        else if (state[0]==-1)
     11335          nAt0++;
     11336        else
     11337          nAt1++;
     11338      }
     11339    }
     11340    if (handler_->logLevel()>2) {
     11341      if (printLevel>1||(printLevel==1&&nBadFixed)) {
     11342        printf("%d switches, %d at 0, %d at 1, %d between - %d bad values (%d when fixed)\n",
     11343               nBetween+nAt0+nAt1,nAt0,nAt1,nBetween,nBad,nBadFixed);
     11344        if (nBadFixed && printLevel!=3)
     11345          checkAssociated(solver,solution,3);
     11346      }
     11347    }
     11348  }
     11349  return nBad;
     11350}
     11351#endif
    1103511352// Set numberBeforeTrust in all objects
    1103611353void
     
    1139011707        memcpy(saveUpper, getColUpper(), numberColumns*sizeof(double));
    1139111708        memcpy(saveLower, getColLower(), numberColumns*sizeof(double));
     11709        //#define CLP_INVESTIGATE4
     11710#ifdef CLP_INVESTIGATE4
     11711        {
     11712          int nBad=checkAssociated(solver_,solver_->getColSolution(),1);
     11713          if (nBad)
     11714            checkAssociated(solver_,solver_->getColSolution(),3);
     11715          double largestInfeasibility = 0.0;
     11716          double primalTolerance ;
     11717          double offset;
     11718          solver_->getDblParam(OsiObjOffset, offset);
     11719          solver_->getDblParam(OsiPrimalTolerance, primalTolerance) ;
     11720          const double *objective = getObjCoefficients() ;
     11721          const double * rowLower = solver_->getRowLower() ;
     11722          const double * rowUpper = solver_->getRowUpper() ;
     11723          const double * columnLower = solver_->getColLower() ;
     11724          const double * columnUpper = solver_->getColUpper() ;
     11725          int numberRows = solver_->getNumRows() ;
     11726          double *rowActivity = new double[numberRows] ;
     11727          memset(rowActivity, 0, numberRows*sizeof(double)) ;
     11728          double *rowSum = new double[numberRows] ;
     11729          memset(rowSum, 0, numberRows*sizeof(double)) ;
     11730          int * marked = new int [numberColumns];
     11731          for (int i=0;i<numberColumns;i++)
     11732            marked[i]=-1;
     11733          for (int i=0;i<numberIntegers_;i++)
     11734            marked[integerVariable_[i]]=-2;
     11735          if ((moreSpecialOptions2_&4)!=0) {
     11736            for (int i=0;i<numberObjects_;i++) {
     11737              CbcSwitchingBinary * object = dynamic_cast<CbcSwitchingBinary *> (object_[i]);
     11738              if (object) {
     11739                int iColumn = object->columnNumber();
     11740                const int * other = object->otherVariable();
     11741                marked[iColumn]=-3-other[0];
     11742                int n=object->numberOther();
     11743                for (int k=0;k<n;k++)
     11744                  marked[other[k]]=iColumn;
     11745              }
     11746            }
     11747          }
     11748          const double * element = solver_->getMatrixByCol()->getElements();
     11749          const int * row = solver_->getMatrixByCol()->getIndices();
     11750          const CoinBigIndex * columnStart = solver_->getMatrixByCol()->getVectorStarts();
     11751          const int * columnLength = solver_->getMatrixByCol()->getVectorLengths();
     11752          const CoinPackedMatrix * rowCopy = solver_->getMatrixByRow();
     11753          const int * column = rowCopy->getIndices();
     11754          const int * rowLength = rowCopy->getVectorLengths();
     11755          const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     11756          const double * elementByRow = rowCopy->getElements();
     11757          double objValue=-offset;
     11758          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     11759            double value = solution[iColumn];
     11760            objValue += value*objective[iColumn];
     11761            if (value>columnUpper[iColumn]) {
     11762              if (value-columnUpper[iColumn]>1.0e-8)
     11763                printf("column %d has value %.12g above %.12g\n",iColumn,value,columnUpper[iColumn]);
     11764              value=columnUpper[iColumn];
     11765            } else if (value<columnLower[iColumn]) {
     11766              if (value-columnLower[iColumn]<-1.0e-8)
     11767                printf("column %d has value %.12g below %.12g\n",iColumn,value,columnLower[iColumn]);
     11768              value=columnLower[iColumn];
     11769            }
     11770            if (value) {
     11771              CoinBigIndex start = columnStart[iColumn];
     11772              CoinBigIndex end = start + columnLength[iColumn];
     11773              for (CoinBigIndex j = start; j < end; j++) {
     11774                int iRow = row[j];
     11775                rowActivity[iRow] += value * element[j];
     11776                rowSum[iRow] += fabs(value * element[j]);
     11777              }
     11778            }
     11779          }
     11780          for (int i = 0 ; i < numberRows ; i++) {
     11781#if 0 //def CLP_INVESTIGATE
     11782            double inf;
     11783            inf = rowLower[i] - rowActivity[i];
     11784            if (inf > primalTolerance)
     11785              printf("Row %d inf %g sum %g %g <= %g <= %g\n",
     11786                     i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
     11787            inf = rowActivity[i] - rowUpper[i];
     11788            if (inf > primalTolerance)
     11789              printf("Row %d inf %g sum %g %g <= %g <= %g\n",
     11790                     i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
     11791#endif
     11792            double infeasibility = CoinMax(rowActivity[i]-rowUpper[i],
     11793                                           rowLower[i]-rowActivity[i]);
     11794            // but allow for errors
     11795            double factor = CoinMax(1.0,rowSum[i]*1.0e-3);
     11796            if (infeasibility>largestInfeasibility*factor) {
     11797              largestInfeasibility = infeasibility/factor;
     11798              printf("Ainf of %g on row %d sum %g scaled %g\n",
     11799                     infeasibility,i,rowSum[i],largestInfeasibility);
     11800              if (infeasibility>1.0e10) {
     11801                for (CoinBigIndex j=rowStart[i];
     11802                     j<rowStart[i]+rowLength[i];j++) {
     11803                  printf("col %d element %g marked %d\n",
     11804                         column[j],elementByRow[j],marked[column[j]]);
     11805                }
     11806              }
     11807            }
     11808          }
     11809          delete [] rowActivity ;
     11810          delete [] rowSum;
     11811          delete [] marked;
     11812          if (largestInfeasibility > 10.0*primalTolerance)
     11813            printf("Alargest infeasibility is %g - obj %g\n", largestInfeasibility,objValue);
     11814          else
     11815            printf("Afeasible (%g) - obj %g\n", largestInfeasibility,objValue);
     11816        }
     11817#endif
    1139211818        // point to useful information
    1139311819        OsiBranchingInformation usefulInfo = usefulInformation();
     
    1140211828        for (i = 0; i < numberObjects_; i++)
    1140311829            object_[i]->feasibleRegion(solver_, &usefulInfo);
     11830#ifdef CLP_INVESTIGATE4
     11831        {
     11832          int nBad=checkAssociated(solver_,solver_->getColSolution(),1);
     11833          if (nBad)
     11834            checkAssociated(solver_,solver_->getColSolution(),3);
     11835          double largestInfeasibility = 0.0;
     11836          double primalTolerance ;
     11837          double offset;
     11838          solver_->getDblParam(OsiObjOffset, offset);
     11839          solver_->getDblParam(OsiPrimalTolerance, primalTolerance) ;
     11840          const double *objective = getObjCoefficients() ;
     11841          const double * rowLower = solver_->getRowLower() ;
     11842          const double * rowUpper = solver_->getRowUpper() ;
     11843          const double * columnLower = solver_->getColLower() ;
     11844          const double * columnUpper = solver_->getColUpper() ;
     11845          int numberRows = solver_->getNumRows() ;
     11846          double *rowActivity = new double[numberRows] ;
     11847          memset(rowActivity, 0, numberRows*sizeof(double)) ;
     11848          double *rowSum = new double[numberRows] ;
     11849          memset(rowSum, 0, numberRows*sizeof(double)) ;
     11850          int * marked = new int [numberColumns];
     11851          for (int i=0;i<numberColumns;i++)
     11852            marked[i]=-1;
     11853          for (int i=0;i<numberIntegers_;i++)
     11854            marked[integerVariable_[i]]=-2;
     11855          if ((moreSpecialOptions2_&4)!=0) {
     11856            for (int i=0;i<numberObjects_;i++) {
     11857              CbcSwitchingBinary * object = dynamic_cast<CbcSwitchingBinary *> (object_[i]);
     11858              if (object) {
     11859                int iColumn = object->columnNumber();
     11860                const int * other = object->otherVariable();
     11861                marked[iColumn]=-3-other[0];
     11862                int n=object->numberOther();
     11863                for (int k=0;k<n;k++)
     11864                  marked[other[k]]=iColumn;
     11865              }
     11866            }
     11867          }
     11868          const double * element = solver_->getMatrixByCol()->getElements();
     11869          const int * row = solver_->getMatrixByCol()->getIndices();
     11870          const CoinBigIndex * columnStart = solver_->getMatrixByCol()->getVectorStarts();
     11871          const int * columnLength = solver_->getMatrixByCol()->getVectorLengths();
     11872          const CoinPackedMatrix * rowCopy = solver_->getMatrixByRow();
     11873          const int * column = rowCopy->getIndices();
     11874          const int * rowLength = rowCopy->getVectorLengths();
     11875          const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     11876          const double * elementByRow = rowCopy->getElements();
     11877          double objValue=-offset;
     11878          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     11879            double value = solution[iColumn];
     11880            objValue += value*objective[iColumn];
     11881            if (value>columnUpper[iColumn]) {
     11882              if (value-columnUpper[iColumn]>1.0e-8)
     11883                printf("column %d has value %.12g above %.12g\n",iColumn,value,columnUpper[iColumn]);
     11884              value=columnUpper[iColumn];
     11885            } else if (value<columnLower[iColumn]) {
     11886              if (value-columnLower[iColumn]<-1.0e-8)
     11887                printf("column %d has value %.12g below %.12g\n",iColumn,value,columnLower[iColumn]);
     11888              value=columnLower[iColumn];
     11889            }
     11890            if (value) {
     11891              CoinBigIndex start = columnStart[iColumn];
     11892              CoinBigIndex end = start + columnLength[iColumn];
     11893              for (CoinBigIndex j = start; j < end; j++) {
     11894                int iRow = row[j];
     11895                rowActivity[iRow] += value * element[j];
     11896                rowSum[iRow] += fabs(value * element[j]);
     11897              }
     11898            }
     11899          }
     11900          for (int i = 0 ; i < numberRows ; i++) {
     11901#if 0 //def CLP_INVESTIGATE
     11902            double inf;
     11903            inf = rowLower[i] - rowActivity[i];
     11904            if (inf > primalTolerance)
     11905              printf("Row %d inf %g sum %g %g <= %g <= %g\n",
     11906                     i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
     11907            inf = rowActivity[i] - rowUpper[i];
     11908            if (inf > primalTolerance)
     11909              printf("Row %d inf %g sum %g %g <= %g <= %g\n",
     11910                     i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
     11911#endif
     11912            double infeasibility = CoinMax(rowActivity[i]-rowUpper[i],
     11913                                           rowLower[i]-rowActivity[i]);
     11914            // but allow for errors
     11915            double factor = CoinMax(1.0,rowSum[i]*1.0e-3);
     11916            if (infeasibility>largestInfeasibility*factor) {
     11917              largestInfeasibility = infeasibility/factor;
     11918              printf("inf of %g on row %d sum %g scaled %g\n",
     11919                     infeasibility,i,rowSum[i],largestInfeasibility);
     11920              if (infeasibility>1.0e10) {
     11921                for (CoinBigIndex j=rowStart[i];
     11922                     j<rowStart[i]+rowLength[i];j++) {
     11923                  printf("col %d element %g marked %d\n",
     11924                         column[j],elementByRow[j],marked[column[j]]);
     11925                }
     11926              }
     11927            }
     11928          }
     11929          delete [] rowActivity ;
     11930          delete [] rowSum;
     11931          delete [] marked;
     11932          if (largestInfeasibility > 10.0*primalTolerance)
     11933            printf("Largest infeasibility is %g - obj %g\n", largestInfeasibility,objValue);
     11934          else
     11935            printf("Feasible (%g) - obj %g\n", largestInfeasibility,objValue);
     11936        }
     11937#endif
    1140411938        // If relaxed then leave bounds on basic variables
    1140511939        if (fixVariables == -1 && (specialOptions_&16) == 0) {
     
    1142211956        }
    1142311957        // We can switch off check
    11424         if ((specialOptions_&4) == 0) {
     11958        if ((specialOptions_&4) == 0 && (moreSpecialOptions2_&10) != 8) {
    1142511959            if ((specialOptions_&2) == 0 && solverCharacteristics_->warmStart()) {
    1142611960                /*
     
    1144811982            solver_->setHintParam(OsiDoDualInInitial, true, OsiHintTry);
    1144911983            solver_->initialSolve();
     11984#ifdef SWITCH_VARIABLES
     11985            if (solver_->isProvenOptimal()) {
     11986              int nBad=checkAssociated(solver_,solver_->getColSolution(),1);
     11987              if (nBad)
     11988                checkAssociated(solver_,solver_->getColSolution(),3);
     11989            }
     11990#endif
    1145011991#ifdef JJF_ZERO
    1145111992            if (solver_->isProvenOptimal()) {
    11452                 solver_->writeMps("feasible");
     11993              solver_->writeMpsNative("feasible.mps",NULL,NULL,2);
     11994#ifdef COIN_HAS_CLP
     11995              OsiClpSolverInterface * clpSolver
     11996                = dynamic_cast<OsiClpSolverInterface *> (solver_);
     11997              if (clpSolver ) {
     11998                clpSolver->getModelPtr()->writeBasis("feasible.bas",true);
     11999              }
     12000#endif
    1145312001                printf("XXXXXXXXXXXX - saving feasible\n");
    1145412002            }
     
    1162312171                const CoinBigIndex * columnStart = solver_->getMatrixByCol()->getVectorStarts();
    1162412172                const int * columnLength = solver_->getMatrixByCol()->getVectorLengths();
    11625                 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    11626                     double value = solution[iColumn];
     12173                double offset;
     12174                solver_->getDblParam(OsiObjOffset, offset);
     12175                double objValue=-offset;
     12176                const double *objective = getObjCoefficients() ;
     12177                for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     12178                  double value = solution[iColumn];
     12179                  objValue += value*objective[iColumn];
    1162712180                    if (value) {
    1162812181                        CoinBigIndex start = columnStart[iColumn];
     
    1164712200                               i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
    1164812201#endif
    11649                     largestInfeasibility = CoinMax(largestInfeasibility,
    11650                                                    rowLower[i] - rowActivity[i]);
    11651                     largestInfeasibility = CoinMax(largestInfeasibility,
    11652                                                    rowActivity[i] - rowUpper[i]);
     12202                    double infeasibility = CoinMax(rowActivity[i]-rowUpper[i],
     12203                                                   rowLower[i]-rowActivity[i]);
     12204                    // but allow for errors
     12205                    double factor = CoinMax(1.0,rowSum[i]*1.0e-3);
     12206                    if (infeasibility>largestInfeasibility*factor) {
     12207                      largestInfeasibility = infeasibility/factor;
     12208                      //printf("inf of %g on row %d sum %g scaled %g\n",
     12209                      //     infeasibility,i,rowSum[i],largestInfeasibility);
     12210                    }
    1165312211                }
    1165412212                delete [] rowActivity ;
    1165512213                delete [] rowSum;
     12214                if (handler_->logLevel()>2) {
     12215                  if (largestInfeasibility > 10.0*primalTolerance)
     12216                    printf("BLargest infeasibility is %g - obj %g (%g)\n", largestInfeasibility,objValue,objectiveValue);
     12217                  else
     12218                    printf("BFeasible (%g) - obj %g %g\n", largestInfeasibility,objValue,objectiveValue);
     12219                }
     12220                //if (fabs(objValue-objectiveValue)>1.0e-7*fabs(objectiveValue)) {
     12221                //printf("Bad obj values\n");
     12222                objectiveValue = objValue;
     12223                //}
    1165612224#ifdef CLP_INVESTIGATE
    1165712225                if (largestInfeasibility > 10.0*primalTolerance)
     
    1186712435        assert(basis != NULL);
    1186812436        objectiveValue = checkSolution(cutoff, solution, fixVariables, objectiveValue);
    11869         if (saveObjectiveValue + 1.0e-3 < objectiveValue) {
     12437        if (cutoff>1.0e40&&objectiveValue<1.0e10)
     12438          saveObjectiveValue = objectiveValue; // take anyway
     12439        if (saveObjectiveValue + 1.0e-3 +1.0e-7*fabs(saveObjectiveValue)
     12440            < objectiveValue) {
    1187012441#if COIN_DEVELOP>1
    1187112442            printf("First try at solution had objective %.16g, rechecked as %.16g\n",
     
    1199612567                    cutoff -= 2.0e-5;
    1199712568            }
     12569            if (!parentModel_&&(moreSpecialOptions2_&2)!=0) {
     12570              // put back objective
     12571              solver_->setObjective(continuousSolver_->getObjCoefficients());
     12572              double offset;
     12573              continuousSolver_->getDblParam(OsiObjOffset,offset);
     12574              solver_->setDblParam(OsiObjOffset,offset);
     12575              moreSpecialOptions2_ &= ~2;
     12576            }
    1199812577            // This is not correct - that way cutoff can go up if maximization
    1199912578            //double direction = solver_->getObjSense();
     
    1254513124                    } else {
    1254613125                        // relax a bit
    12547                         value = CoinMin(saveUpper, value + 1.0e-5 * (fabs(saveUpper) + 1));
     13126                        value = CoinMin(saveUpper, value + 1.0e-8 * (fabs(saveUpper) + 1));
    1254813127                    }
    1254913128                    if (value - saveLower < 1.0e-7)
     
    1255913138                    } else {
    1256013139                        // relax a bit
    12561                         value = CoinMax(saveLower, value - 1.0e-5 * (fabs(saveLower) + 1));
     13140                        value = CoinMax(saveLower, value - 1.0e-8 * (fabs(saveLower) + 1));
    1256213141                    }
    1256313142                    if (saveUpper - value < 1.0e-7)
     
    1260713186                                    newLower = CoinMax(lower[jColumn],
    1260813187                                                       newLower
    12609                                                        - 1.0e-5 * (fabs(lower[jColumn]) + 1));
     13188                                                       - 1.0e-8 * (fabs(lower[jColumn]) + 1));
    1261013189                                    newUpper = CoinMin(upper[jColumn],
    1261113190                                                       newUpper
    12612                                                        + 1.0e-5 * (fabs(upper[jColumn]) + 1));
     13191                                                       + 1.0e-8 * (fabs(upper[jColumn]) + 1));
    1261313192                                }
    1261413193                                solver->setColLower(jColumn, newLower);
     
    1314513724    if (probingInfo_ && currentDepth_ > 0) {
    1314613725        int nFix = probingInfo_->fixColumns(*solver);
     13726#ifdef SWITCH_VARIABLES
     13727        if (nFix>0)
     13728          fixAssociated(solver_,0);
     13729#endif
    1314713730        if (nFix < 0) {
    1314813731#ifdef COIN_HAS_CLP
     
    1324213825#else
    1324313826    solver->resolve();
     13827#endif
     13828#ifdef SWITCH_VARIABLES
     13829    if (solver_->isProvenOptimal()) {
     13830      int nBad=checkAssociated(solver_,solver_->getColSolution(),0);
     13831      if (nBad)
     13832        checkAssociated(solver_,solver_->getColSolution(),1);
     13833    }
    1324413834#endif
    1324513835    return solver->isProvenOptimal() ? 1 : 0;
     
    1473215322    int save2 = maximumDepth_;
    1473315323    int retCode = addCuts(node, lastws, numberFixedNow_ > numberFixedAtRoot_);
     15324#ifdef SWITCH_VARIABLES
     15325    fixAssociated(solver_,0);
     15326#endif
    1473415327    //if (save1<maximumNumberCuts_) {
    1473515328    // increased
     
    1571416307                            }
    1571516308                        } else if (ifSol < 0)   { // just returning an estimate
    15716                             estValue = CoinMin(heurValue, estValue) ;
     16309                            estValue = heurValue; //CoinMin(heurValue, estValue) ;
    1571716310                            heurValue = saveValue ;
    1571816311                        }
     
    1604716640        }
    1604816641    }
    16049 #ifdef CLP_INVESTIGATE4
     16642#ifdef CLP_INVESTIGATE5
    1605016643    if (priority)
    1605116644        printf("Before fathom %d not trusted out of %d\n",
     
    1631416907        OsiClpSolverInterface * clpSolver
    1631516908        = dynamic_cast<OsiClpSolverInterface *> (solver_);
    16316         if (clpSolver && numberNodes_ >= numberNodes && numberNodes_ < 2*numberNodes) {
     16909        if (clpSolver && numberNodes_ >= numberNodes && numberNodes_ < 2*numberNodes  && clpSolver->getNumRows() < 10000) {
    1631716910            if (numberIterations_ < (numberSolves_ + numberNodes_)*10) {
    1631816911                //if (numberIterations_<numberNodes_*20) {
  • trunk/Cbc/src/CbcModel.hpp

    r1927 r1943  
    513513    void findIntegers(bool startAgain, int type = 0);
    514514
     515#ifdef SWITCH_VARIABLES
     516    /// Convert Dynamic to Switching
     517    int findSwitching();
     518    /// Fix associated variables
     519    int fixAssociated(OsiSolverInterface * solver,int cleanBasis);
     520    /// Debug associated variables
     521    int checkAssociated(const OsiSolverInterface * solver,
     522                        const double * solution, int printLevel);
     523#endif
    515524    //@}
    516525
     
    757766    inline int getCurrentPassNumber() const {
    758767        return currentPassNumber_;
     768    }
     769    /** Set current cut pass number in this round of cuts.
     770        (1 is first pass) */
     771    inline void setCurrentPassNumber(int value) {
     772        currentPassNumber_ = value;
    759773    }
    760774
     
    17701784        22 bit (4194304) - do not initialize random seed in solver (user has)
    17711785        23 bit (8388608) - leave solver_ with cuts
     1786        24 bit (16777216) - just get feasible if no cutoff
    17721787    */
    17731788    inline void setSpecialOptions(int value) {
     
    18361851    inline int moreSpecialOptions() const {
    18371852        return moreSpecialOptions_;
     1853    }
     1854    /** Set more more special options
     1855        0 bit (1) - find switching variables
     1856        1 bit (2) - using fake objective until solution
     1857        2 bit (4) - switching variables exist
     1858        3 bit (8) - skip most of setBestSolution checks
     1859        4 bit (16) - very lightweight preprocessing in smallB&B
     1860        5 bit (32) - event handler needs to be cloned when parallel
     1861    */
     1862    inline void setMoreSpecialOptions2(int value) {
     1863        moreSpecialOptions2_ = value;
     1864    }
     1865    /// Get more special options2
     1866    inline int moreSpecialOptions2() const {
     1867        return moreSpecialOptions2_;
    18381868    }
    18391869    /// Set cutoff as constraint
     
    25692599        18 bit (262144) - donor CbcModel
    25702600        19 bit (524288) - recipient CbcModel
     2601        20 bit (1048576) - waiting for sub model to return
     2602        22 bit (4194304) - do not initialize random seed in solver (user has)
     2603        23 bit (8388608) - leave solver_ with cuts
     2604        24 bit (16777216) - just get feasible if no cutoff
    25712605    */
    25722606    int specialOptions_;
     
    25862620    */
    25872621    int moreSpecialOptions_;
     2622    /** More more special options
     2623        0 bit (1) - find switching variables
     2624        1 bit (2) - using fake objective until solution
     2625        2 bit (4) - switching variables exist
     2626        3 bit (8) - skip most of setBestSolution checks
     2627        4 bit (16) - very lightweight preprocessing in smallB&B
     2628        5 bit (32) - event handler needs to be cloned when parallel
     2629    */
     2630    int moreSpecialOptions2_;
    25882631    /// User node comparison function
    25892632    CbcCompareBase * nodeCompare_;
  • trunk/Cbc/src/CbcNode.cpp

    r1886 r1943  
    15381538    // Set guessed solution value
    15391539    guessedObjectiveValue_ = objectiveValue_ + estimatedDegradation;
     1540    //printf ("Node %d depth %d unsatisfied %d sum %g obj %g guess %g\n",
     1541    //      model->getNodeCount(),depth_,numberUnsatisfied_,
     1542    //      sumInfeasibilities_,objectiveValue_,guessedObjectiveValue_);
    15401543    /*
    15411544      Cleanup, then we're outta here.
     
    36693672    // Set guessed solution value
    36703673    guessedObjectiveValue_ = objectiveValue_ + estimatedDegradation;
     3674    int kColumn=-1;
     3675    if (branch_) {
     3676      CbcObject * obj = (dynamic_cast<CbcBranchingObject *>(branch_))->object();
     3677      CbcSimpleInteger * branchObj =
     3678        dynamic_cast <CbcSimpleInteger *>(obj) ;
     3679      if (branchObj) {
     3680        kColumn=branchObj->columnNumber();
     3681      }
     3682    }
     3683    if (model->logLevel()>1)
     3684    printf ("Node %d depth %d unsatisfied %d sum %g obj %g guess %g branching on %d\n",
     3685            model->getNodeCount(),depth_,numberUnsatisfied_,
     3686            sumInfeasibilities_,objectiveValue_,guessedObjectiveValue_,
     3687            kColumn);
    36713688#ifdef DO_ALL_AT_ROOT
    36723689    if (strongType) {
  • trunk/Cbc/src/CbcNodeInfo.cpp

    r1899 r1943  
    507507{
    508508    active_ &= (~mode);
    509 }
    510 
     509    if (mode==7) {
     510      for (int i = 0; i < numberCuts_; i++) {
     511        delete cuts_[i];
     512        cuts_[i] = NULL;
     513      }
     514      delete [] cuts_;
     515      cuts_=NULL;
     516      numberCuts_=0;
     517    }
     518}
     519
  • trunk/Cbc/src/CbcSimpleInteger.cpp

    r1899 r1943  
    2222#include "CbcMessage.hpp"
    2323#include "CbcSimpleInteger.hpp"
     24#include "CbcSimpleIntegerDynamicPseudoCost.hpp"
    2425#include "CbcBranchActual.hpp"
    2526#include "CoinSort.hpp"
     
    145146    solver->setColLower(columnNumber_, newValue);
    146147    solver->setColUpper(columnNumber_, newValue);
     148#ifdef SWITCH_VARIABLES
     149    const CbcSwitchingBinary * sObject = dynamic_cast<const CbcSwitchingBinary *> (this);
     150    if (sObject)
     151      sObject->setAssociatedBounds(solver,1);
     152#endif
    147153    return fabs(value - newValue);
    148154}
     
    309315    up_[0] = ceil(value_);
    310316    up_[1] = model_->getColUpper()[iColumn];
     317    // fix extreme cases
     318    if (up_[0]==1.0)
     319      down_[1]=0.0;
     320    if (down_[1]==0.0)
     321      up_[0]=1.0;
    311322}
    312323// Useful constructor for fixing
     
    555566    if (nlb < olb + 1.0e-8 && nub > oub - 1.0e-8 && false)
    556567        printf("bad null change for column %d - bounds %g,%g\n", iColumn, olb, oub);
     568#endif
     569#ifdef SWITCH_VARIABLES
     570    if (model_->logLevel()>2)
     571      printf("for column %d - old bounds %g,%g - new %g,%g\n", iColumn, olb, oub,
     572             nlb,nub);
     573    CbcSwitchingBinary * sObject = dynamic_cast<CbcSwitchingBinary *> (originalCbcObject_);
     574    if (sObject)
     575      sObject->setAssociatedBounds();
     576    //(dynamic_cast<CbcSimpleInteger *>(originalCbcObject_))->setAssociatedBounds();
    557577#endif
    558578    return 0.0;
  • trunk/Cbc/src/CbcSimpleInteger.hpp

    r1899 r1943  
    210210    virtual CbcBranchingObject * createCbcBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) ;
    211211    /// Fills in a created branching object
    212     void fillCreateBranch(CbcIntegerBranchingObject * branching, const OsiBranchingInformation * info, int way) ;
     212  /*virtual*/ void fillCreateBranch(CbcIntegerBranchingObject * branching, const OsiBranchingInformation * info, int way) ;
    213213
    214214    using CbcObject::solverBranch ;
  • trunk/Cbc/src/CbcSimpleIntegerDynamicPseudoCost.cpp

    r1899 r1943  
    2424#include "CoinError.hpp"
    2525#include "CbcSimpleIntegerDynamicPseudoCost.hpp"
     26#ifdef COIN_HAS_CLP
     27#include "OsiClpSolverInterface.hpp"
     28#endif
    2629#ifdef COIN_DEVELOP
    2730typedef struct {
     
    840843    if (!info->hotstartSolution_ && priority_ != -999) {
    841844#ifndef NDEBUG
     845#ifndef SWITCH_VARIABLES
    842846        double nearest = floor(value + 0.5);
    843847        assert (fabs(value - nearest) > info->integerTolerance_);
     848#endif
    844849#endif
    845850    } else if (info->hotstartSolution_) {
     
    10161021            if (fabs(value - nearest) > integerTolerance)
    10171022                unsatisfied++;
     1023#ifdef SWITCH_VARIABLES
     1024            const CbcSwitchingBinary * sObject = dynamic_cast<const CbcSwitchingBinary *> (this);
     1025            if (sObject) {
     1026              int state[3],nBadFixed;
     1027              unsatisfied +=
     1028                sObject->checkAssociatedBounds(solver,solution,0,
     1029                                               state,nBadFixed);
     1030            }
     1031#endif
    10181032        }
    10191033    }
     
    13981412    return CbcCompareRanges(thisBd, otherBd, replaceIfOverlap);
    13991413}
    1400 
     1414#ifdef SWITCH_VARIABLES
     1415/** Default Constructor
     1416
     1417  Equivalent to an unspecified binary variable.
     1418*/
     1419CbcSwitchingBinary::CbcSwitchingBinary ()
     1420        : CbcSimpleIntegerDynamicPseudoCost(),
     1421          zeroLowerBound_(NULL),
     1422          oneLowerBound_(NULL),
     1423          zeroUpperBound_(NULL),
     1424          oneUpperBound_(NULL),
     1425          otherVariable_(NULL),
     1426          numberOther_(0),
     1427          type_(0)
     1428{
     1429}
     1430
     1431/** Useful constructor
     1432*/
     1433CbcSwitchingBinary::CbcSwitchingBinary (CbcSimpleIntegerDynamicPseudoCost * oldObject,
     1434                                        int nOdd,const int * other, const int * otherRow)
     1435  : CbcSimpleIntegerDynamicPseudoCost(*oldObject),
     1436    zeroLowerBound_(NULL),
     1437    oneLowerBound_(NULL),
     1438    zeroUpperBound_(NULL),
     1439    oneUpperBound_(NULL),
     1440    otherVariable_(NULL),
     1441    numberOther_(0),
     1442    type_(0)
     1443{
     1444  if (nOdd)
     1445    type_=2;
     1446  const CoinPackedMatrix * rowCopy = model_->solver()->getMatrixByRow();
     1447  const int * column = rowCopy->getIndices();
     1448  //const int * rowLength = rowCopy->getVectorLengths();
     1449  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     1450  //const double * rowLower = model_->solver()->getRowLower();
     1451  const double * rowUpper = model_->solver()->getRowUpper();
     1452  const double * columnLower = model_->solver()->getColLower();
     1453  const double * columnUpper = model_->solver()->getColUpper();
     1454  const double * element = rowCopy->getElements();
     1455  int last = other[0];
     1456  int nPair=0;
     1457  int nInGroup=1;
     1458  for (int i=1;i<=nOdd;i++) {
     1459    if (other[i]==last) {
     1460      nInGroup++;
     1461    } else {
     1462      if (nInGroup>2 && model_->logLevel()>2)
     1463        printf("%d in group for column %d - some redundancy\n",
     1464               nInGroup,columnNumber_);
     1465      nPair++;
     1466      last=other[i];
     1467      nInGroup=1;
     1468    }
     1469  }
     1470  zeroLowerBound_ = new double [4*nPair];
     1471  oneLowerBound_ = zeroLowerBound_+nPair;
     1472  zeroUpperBound_ = oneLowerBound_+nPair;
     1473  oneUpperBound_ = zeroUpperBound_+nPair;
     1474  otherVariable_ = new int [nPair];
     1475  numberOther_ = nPair;
     1476  if (nPair>1&&model_->logLevel()>2)
     1477    printf("%d pairs for column %d\n",
     1478               nPair,columnNumber_);
     1479  // Now fill
     1480  last = other[0];
     1481  nPair=0;
     1482  int rows[20];
     1483  rows[0]= otherRow[0];
     1484  nInGroup=1;
     1485  for (int i=1;i<=nOdd;i++) {
     1486    if (other[i]==last) {
     1487      rows[nInGroup++]=otherRow[i];
     1488    } else {
     1489      double newLowerZero=0.0;
     1490      double newUpperZero=COIN_DBL_MAX;
     1491      double newLowerOne=0.0;
     1492      double newUpperOne=COIN_DBL_MAX;
     1493      int cColumn=-1;
     1494      for (int j=0;j<nInGroup;j++) {
     1495        int iRow = rows[j];
     1496        CoinBigIndex k = rowStart[iRow];
     1497        double bValue, cValue;
     1498        if (column[k]==columnNumber_) {
     1499          bValue=element[k];
     1500          cValue=element[k+1];
     1501          cColumn=column[k+1];
     1502        } else {
     1503          bValue=element[k+1];
     1504          cValue=element[k];
     1505          cColumn=column[k];
     1506        }
     1507        if (rowUpper[iRow]) {
     1508          // G row - convert to L
     1509          bValue=-bValue;
     1510          cValue=-cValue;
     1511        }
     1512        if (bValue>0.0) {
     1513          // binary*abs(bValue) <= continuous*abs(cValue);
     1514          newLowerOne = -bValue/cValue;
     1515        } else {
     1516          // binary*abs(bValue) >= continuous*abs(cValue);
     1517          newUpperOne = -bValue/cValue;
     1518          newUpperZero = 0.0;
     1519        }
     1520      }
     1521      zeroLowerBound_[nPair]=newLowerZero;
     1522      oneLowerBound_[nPair]=newLowerOne;
     1523      zeroUpperBound_[nPair]=newUpperZero;
     1524      oneUpperBound_[nPair]=newUpperOne;
     1525      // make current bounds tight
     1526      double newLower = CoinMin(newLowerZero,newLowerOne);
     1527      if (newLower>columnLower[cColumn])
     1528        model_->solver()->setColLower(cColumn,newLower);
     1529      double newUpper = CoinMax(newUpperZero,newUpperOne);
     1530      if (newUpper<columnUpper[cColumn])
     1531        model_->solver()->setColUpper(cColumn,newUpper);
     1532      otherVariable_[nPair++]=cColumn;
     1533      last=other[i];
     1534      rows[0] = otherRow[i];
     1535      nInGroup=1;
     1536    }
     1537  }
     1538}
     1539// Copy constructor
     1540CbcSwitchingBinary::CbcSwitchingBinary ( const CbcSwitchingBinary & rhs)
     1541        : CbcSimpleIntegerDynamicPseudoCost(rhs),
     1542          numberOther_(rhs.numberOther_),
     1543          type_(rhs.type_)
     1544{
     1545  zeroLowerBound_ = CoinCopyOfArray(rhs.zeroLowerBound_,4*numberOther_);
     1546  oneLowerBound_ = zeroLowerBound_+numberOther_;
     1547  zeroUpperBound_ = oneLowerBound_+numberOther_;
     1548  oneUpperBound_ = zeroUpperBound_+numberOther_;
     1549  otherVariable_ = CoinCopyOfArray(rhs.otherVariable_,numberOther_);
     1550}
     1551
     1552// Clone
     1553CbcObject *
     1554CbcSwitchingBinary::clone() const
     1555{
     1556    return new CbcSwitchingBinary(*this);
     1557}
     1558
     1559// Assignment operator
     1560CbcSwitchingBinary &
     1561CbcSwitchingBinary::operator=( const CbcSwitchingBinary & rhs)
     1562{
     1563    if (this != &rhs) {
     1564        CbcSimpleIntegerDynamicPseudoCost::operator=(rhs);
     1565        numberOther_=rhs.numberOther_;
     1566        type_ = rhs.type_;
     1567        delete [] zeroLowerBound_;
     1568        delete [] otherVariable_;
     1569        zeroLowerBound_ = CoinCopyOfArray(rhs.zeroLowerBound_,4*numberOther_);
     1570        oneLowerBound_ = zeroLowerBound_+numberOther_;
     1571        zeroUpperBound_ = oneLowerBound_+numberOther_;
     1572        oneUpperBound_ = zeroUpperBound_+numberOther_;
     1573        otherVariable_ = CoinCopyOfArray(rhs.otherVariable_,numberOther_);
     1574    }
     1575    return *this;
     1576}
     1577
     1578// Destructor
     1579CbcSwitchingBinary::~CbcSwitchingBinary ()
     1580{
     1581  delete [] zeroLowerBound_;
     1582  delete [] otherVariable_;
     1583}
     1584// Add in zero switches
     1585void
     1586CbcSwitchingBinary::addZeroSwitches(int nAdd,const int * columns)
     1587{
     1588  type_ |= 1;
     1589  int nNew = numberOther_+nAdd;
     1590  double * bounds = new double[4*nNew];
     1591  int * other = new int [nNew];
     1592  memcpy(other,otherVariable_,numberOther_*sizeof(int));
     1593  delete [] otherVariable_;
     1594  otherVariable_=other;
     1595  memcpy(bounds,zeroLowerBound_,numberOther_*sizeof(double));
     1596  memcpy(bounds+nNew,oneLowerBound_,numberOther_*sizeof(double));
     1597  memcpy(bounds+2*nNew,zeroUpperBound_,numberOther_*sizeof(double));
     1598  memcpy(bounds+3*nNew,oneUpperBound_,numberOther_*sizeof(double));
     1599  delete [] zeroLowerBound_;
     1600  zeroLowerBound_ = bounds;
     1601  oneLowerBound_ = zeroLowerBound_+nNew;
     1602  zeroUpperBound_ = oneLowerBound_+nNew;
     1603  oneUpperBound_ = zeroUpperBound_+nNew;
     1604  for (int i=0;i<nAdd;i++) {
     1605    zeroLowerBound_[numberOther_]=0.0;
     1606    oneLowerBound_[numberOther_]=0.0;
     1607    zeroUpperBound_[numberOther_]=0.0;
     1608    oneUpperBound_[numberOther_]=COIN_DBL_MAX;
     1609    otherVariable_[numberOther_++]=columns[i];
     1610  }
     1611}
     1612// Same - returns true if contents match(ish)
     1613bool
     1614CbcSwitchingBinary::same(const CbcSwitchingBinary * otherObject) const
     1615{
     1616  bool okay = CbcSimpleIntegerDynamicPseudoCost::same(otherObject);
     1617    return okay;
     1618}
     1619double
     1620CbcSwitchingBinary::infeasibility(const OsiBranchingInformation * info,
     1621        int &preferredWay) const
     1622{
     1623    assert (downDynamicPseudoCost_ > 1.0e-40 && upDynamicPseudoCost_ > 1.0e-40);
     1624    double * solution = const_cast<double *>(model_->testSolution());
     1625    const double * lower = model_->getCbcColLower();
     1626    const double * upper = model_->getCbcColUpper();
     1627    double saveValue = solution[columnNumber_];
     1628    if (!lower[columnNumber_]&&upper[columnNumber_]==1.0) {
     1629      double integerTolerance =
     1630        model_->getDblParam(CbcModel::CbcIntegerTolerance);
     1631      if (saveValue<integerTolerance) {
     1632        // check others OK
     1633        bool allGood=true;
     1634        double tolerance;
     1635        model_->solver()->getDblParam(OsiPrimalTolerance, tolerance) ;
     1636        for (int i=0;i<numberOther_;i++) {
     1637          int otherColumn = otherVariable_[i];
     1638          double value = solution[otherColumn];
     1639          if (value<zeroLowerBound_[i]-tolerance||
     1640              value>zeroUpperBound_[i]+tolerance)
     1641            allGood=false;
     1642        }
     1643        if (!allGood)
     1644          solution[columnNumber_]=2.0*integerTolerance;
     1645      } else if (saveValue>1.0-integerTolerance) {
     1646        // check others OK
     1647        bool allGood=true;
     1648        double tolerance;
     1649        model_->solver()->getDblParam(OsiPrimalTolerance, tolerance) ;
     1650        for (int i=0;i<numberOther_;i++) {
     1651          int otherColumn = otherVariable_[i];
     1652          double value = solution[otherColumn];
     1653          if (value<oneLowerBound_[i]-tolerance||
     1654              value>oneUpperBound_[i]+tolerance)
     1655            allGood=false;
     1656        }
     1657        if (!allGood)
     1658          solution[columnNumber_]=1.0-2.0*integerTolerance;
     1659      }
     1660    }
     1661    double inf = CbcSimpleIntegerDynamicPseudoCost::infeasibility(info,preferredWay);
     1662    solution[columnNumber_]=saveValue;
     1663    return inf;
     1664}
     1665// Set associated bounds
     1666int
     1667CbcSwitchingBinary::setAssociatedBounds(OsiSolverInterface * solver,
     1668                                        int cleanBasis) const
     1669{
     1670  if (!solver)
     1671    solver = model_->solver();
     1672#ifdef COIN_HAS_CLP
     1673  OsiClpSolverInterface * clpSolver
     1674    = dynamic_cast<OsiClpSolverInterface *> (solver);
     1675  if (cleanBasis!=1)
     1676    clpSolver=NULL;
     1677#endif
     1678  const double * columnLower = solver->getColLower();
     1679  const double * columnUpper = solver->getColUpper();
     1680  int nChanged=0;
     1681  if (!columnUpper[columnNumber_]) {
     1682#ifdef COIN_HAS_CLP
     1683    if (clpSolver)
     1684      clpSolver->setColumnStatus(columnNumber_,ClpSimplex::isFixed);
     1685#endif
     1686    for (int i=0;i<numberOther_;i++) {
     1687      int otherColumn = otherVariable_[i];
     1688      if (zeroLowerBound_[i]>columnLower[otherColumn]) {
     1689        solver->setColLower(otherColumn,zeroLowerBound_[i]);
     1690        nChanged++;
     1691      }
     1692      if (zeroUpperBound_[i]<columnUpper[otherColumn]) {
     1693        solver->setColUpper(otherColumn,zeroUpperBound_[i]);
     1694#ifdef COIN_DEVELOP
     1695        const double * solution = solver->getColSolution();
     1696        double value = solution[otherColumn];
     1697        if (value - zeroUpperBound_[i] > 1.0e-5 && model_->logLevel()>1)
     1698          printf("value for continuous %d %g - above %g - switch %d is %.12g (ub 0)\n",
     1699                 otherColumn, value, zeroUpperBound_[i],columnNumber_,solution[columnNumber_]);
     1700#endif
     1701        nChanged++;
     1702      }
     1703    }
     1704  } else if (columnLower[columnNumber_]==1.0) {
     1705#ifdef COIN_HAS_CLP
     1706    if (clpSolver)
     1707      clpSolver->setColumnStatus(columnNumber_,ClpSimplex::isFixed);
     1708#endif
     1709    for (int i=0;i<numberOther_;i++) {
     1710      int otherColumn = otherVariable_[i];
     1711      if (oneLowerBound_[i]>columnLower[otherColumn]) {
     1712        solver->setColLower(otherColumn,oneLowerBound_[i]);
     1713        nChanged++;
     1714      }
     1715      if (oneUpperBound_[i]<columnUpper[otherColumn]) {
     1716        solver->setColUpper(otherColumn,oneUpperBound_[i]);
     1717        nChanged++;
     1718      }
     1719    }
     1720  } else if (cleanBasis>=2) {
     1721    // if all OK then can fix
     1722    int state[3];
     1723    int nBadFixed;
     1724    const double * solution = solver->getColSolution();
     1725    if (!checkAssociatedBounds(solver,solution,
     1726                               0,state,nBadFixed)) {
     1727      const double *reducedCost = solver->getReducedCost() ;
     1728      double good=true;
     1729      double integerTolerance =
     1730        model_->getDblParam(CbcModel::CbcIntegerTolerance);
     1731      if (solution[columnNumber_]<integerTolerance) {
     1732        if (cleanBasis==2||reducedCost[columnNumber_]>1.0e-6)
     1733          solver->setColUpper(columnNumber_,0.0);
     1734        else
     1735          good=false;
     1736      } else if (solution[columnNumber_]>1.0-integerTolerance) {
     1737        if (cleanBasis==2||reducedCost[columnNumber_]<-1.0e-6)
     1738          solver->setColLower(columnNumber_,1.0);
     1739        else
     1740          good=false;
     1741      }
     1742      if (good)
     1743        nChanged=setAssociatedBounds(solver,0);
     1744    }
     1745  } else {
     1746    // see if any continuous bounds force binary
     1747    for (int i=0;i<numberOther_;i++) {
     1748      int otherColumn = otherVariable_[i];
     1749      if (columnLower[otherColumn]>zeroUpperBound_[i]) {
     1750        // can't be zero
     1751        solver->setColLower(columnNumber_,1.0);
     1752        nChanged++;
     1753      } else if (columnLower[otherColumn]>oneUpperBound_[i]) {
     1754        // can't be one
     1755        solver->setColUpper(columnNumber_,0.0);
     1756        nChanged++;
     1757      }
     1758      if (columnUpper[otherColumn]<zeroLowerBound_[i]) {
     1759        // can't be zero
     1760        solver->setColLower(columnNumber_,1.0);
     1761        nChanged++;
     1762      } else if (columnUpper[otherColumn]<oneLowerBound_[i]) {
     1763        // can't be one
     1764        solver->setColUpper(columnNumber_,0.0);
     1765        nChanged++;
     1766      }
     1767    }
     1768  }
     1769  return nChanged;
     1770}
     1771// Check associated bounds
     1772int
     1773CbcSwitchingBinary::checkAssociatedBounds(const OsiSolverInterface * solver,
     1774                                          const double * solution, int printLevel,
     1775                                          int state[3], int & nBadFixed) const
     1776{
     1777  state[0] = 0;
     1778  int nBad=0;
     1779  if (!solver)
     1780    solver = model_->solver();
     1781  double tolerance;
     1782  solver->getDblParam(OsiPrimalTolerance, tolerance) ;
     1783  const double * columnLower = solver->getColLower();
     1784  const double * columnUpper = solver->getColUpper();
     1785  double integerTolerance =
     1786    model_->getDblParam(CbcModel::CbcIntegerTolerance);
     1787  bool printIt = printLevel>2 && model_->logLevel()>1;
     1788  if (solution[columnNumber_]<integerTolerance) {
     1789    state[0] = -1;
     1790    for (int i=0;i<numberOther_;i++) {
     1791      int otherColumn = otherVariable_[i];
     1792      if (zeroLowerBound_[i]>solution[otherColumn]+tolerance*5.0) {
     1793        nBad++;
     1794        if (columnUpper[columnNumber_]==0.0) {
     1795          nBadFixed++;
     1796          //printIt=true;
     1797        }
     1798        if (printIt)
     1799          printf("switch %d at zero, other %d at %.12g below bound of %.12g\n",
     1800                 columnNumber_,otherColumn,solution[otherColumn],zeroLowerBound_[i]);
     1801      }
     1802      if (zeroUpperBound_[i]<solution[otherColumn]-tolerance*5.0) {
     1803        nBad++;
     1804        if (columnUpper[columnNumber_]==0.0) {
     1805          nBadFixed++;
     1806          //printIt=true;
     1807        }
     1808        if (printIt)
     1809          printf("switch %d at zero, other %d at %.12g above bound of %.12g\n",
     1810                 columnNumber_,otherColumn,solution[otherColumn],zeroUpperBound_[i]);
     1811      }
     1812    }
     1813  } else if (solution[columnNumber_]>1.0-integerTolerance) {
     1814    state[0] = 1;
     1815    for (int i=0;i<numberOther_;i++) {
     1816      int otherColumn = otherVariable_[i];
     1817      if (oneLowerBound_[i]>solution[otherColumn]+tolerance*5.0) {
     1818        nBad++;
     1819        if (columnLower[columnNumber_]==1.0) {
     1820          nBadFixed++;
     1821          //printIt=true;
     1822        }
     1823        if (printIt)
     1824          printf("switch %d at one, other %d at %.12g below bound of %.12g\n",
     1825                 columnNumber_,otherColumn,solution[otherColumn],oneLowerBound_[i]);
     1826      }
     1827      if (oneUpperBound_[i]<solution[otherColumn]-tolerance*5.0) {
     1828        nBad++;
     1829        if (columnLower[columnNumber_]==1.0) {
     1830          nBadFixed++;
     1831          //printIt=true;
     1832        }
     1833        if (printIt)
     1834          printf("switch %d at one, other %d at %.12g above bound of %.12g\n",
     1835                 columnNumber_,otherColumn,solution[otherColumn],oneUpperBound_[i]);
     1836      }
     1837    }
     1838  } else {
     1839    // in between - compute tight variables
     1840    state[1]=0;
     1841    state[2]=0;
     1842    // for now just compute ones away from bounds
     1843    for (int i=0;i<numberOther_;i++) {
     1844      int otherColumn = otherVariable_[i];
     1845      double otherValue = solution[otherColumn];
     1846      if (otherValue>columnLower[otherColumn]+tolerance&&
     1847          otherValue<columnUpper[otherColumn]-tolerance)
     1848        state[1]++;
     1849    }
     1850  }
     1851  return nBad;
     1852}
     1853#endif
  • trunk/Cbc/src/CbcSimpleIntegerDynamicPseudoCost.hpp

    r1899 r1943  
    7272
    7373    /// Fills in a created branching object
    74     void fillCreateBranch(CbcIntegerBranchingObject * branching, const OsiBranchingInformation * info, int way) ;
     74    // void fillCreateBranch(CbcIntegerBranchingObject * branching, const OsiBranchingInformation * info, int way) ;
    7575
    7676
     
    462462    double changeInGuessed_;
    463463};
     464#ifdef SWITCH_VARIABLES
     465/** Define a single integer class but with associated switched variable
     466    So Binary variable switches on/off a continuous variable
     467    designed for badly scaled problems
     468 */
     469
     470
     471class CbcSwitchingBinary : public CbcSimpleIntegerDynamicPseudoCost {
     472
     473public:
     474
     475    // Default Constructor
     476    CbcSwitchingBinary ();
     477
     478    // Useful constructor
     479    CbcSwitchingBinary (CbcSimpleIntegerDynamicPseudoCost * oldObject,
     480                        int nOdd,const int * other, const int * otherRow);
     481
     482
     483    // Copy constructor
     484    CbcSwitchingBinary ( const CbcSwitchingBinary &);
     485
     486    /// Clone
     487    virtual CbcObject * clone() const;
     488
     489    // Assignment operator
     490    CbcSwitchingBinary & operator=( const CbcSwitchingBinary& rhs);
     491
     492    // Destructor
     493    virtual ~CbcSwitchingBinary ();
     494
     495    /// Add in zero switches
     496    void addZeroSwitches(int nAdd,const int * columns);
     497    /// Infeasibility - large is 0.5
     498    virtual double infeasibility(const OsiBranchingInformation * info,
     499                                 int &preferredWay) const;
     500
     501    /// Same - returns true if contents match(ish)
     502    bool same(const CbcSwitchingBinary * obj) const;
     503    /// Set associated bounds
     504    virtual int setAssociatedBounds(OsiSolverInterface * solver=NULL,
     505                          int cleanBasis=0) const;
     506    /// Check associated bounds
     507    int checkAssociatedBounds(const OsiSolverInterface * solver,const double * solution,
     508                              int printLevel, int state[3], int & nBadFixed) const;
     509    /// Lower bound when binary zero
     510    inline const double * zeroLowerBound() const
     511    { return zeroLowerBound_; }
     512    /// Lower bound when binary one
     513    inline const double * oneLowerBound() const
     514    { return oneLowerBound_; }
     515    /// Upper bound when binary zero
     516    inline const double * zeroUpperBound() const
     517    { return zeroUpperBound_; }
     518    /// Upper bound when binary one
     519    inline const double * oneUpperBound() const
     520    { return oneUpperBound_; }
     521    /** Continuous variable -
     522    */
     523    inline const int * otherVariable() const
     524    { return otherVariable_;}
     525    /// Number of other variables
     526    inline int numberOther() const
     527    { return numberOther_;}
     528    /** Type
     529        1 - single switch
     530        2 - double switch
     531        3 - both
     532    */
     533    inline int type() const
     534    { return type_;}
     535protected:
     536    /// data
     537
     538    /// Lower bound when binary zero
     539    double * zeroLowerBound_;
     540    /// Lower bound when binary one
     541    double * oneLowerBound_;
     542    /// Upper bound when binary zero
     543    double * zeroUpperBound_;
     544    /// Upper bound when binary one
     545    double * oneUpperBound_;
     546    /** Continuous variable -
     547    */
     548    int * otherVariable_;
     549    /// Number of other variables
     550    int numberOther_;
     551    /** Type
     552        1 - single switch
     553        2 - double switch
     554        3 - both
     555    */
     556    int type_;
     557};
    464558#endif
    465 
     559#endif
     560
  • trunk/Cbc/src/CbcSolver.cpp

    r1937 r1943  
    38763876                                          process.setOptions(2+4+8); // no cuts
    38773877                                        cbcPreProcessPointer = & process;
     3878                                        int saveOptions = osiclp->getModelPtr()->moreSpecialOptions();
     3879                                        if ((model_.specialOptions()&16777216)!=0&&
     3880                                            model_.getCutoff()>1.0e30) {
     3881                                          osiclp->getModelPtr()->setMoreSpecialOptions(saveOptions|262144);
     3882                                        }
    38783883                                        solver2 = process.preProcessNonDefault(*saveSolver, translate[preProcess], numberPasses,
    38793884                                                                               tunePreProcess);
     
    38813886                                          saveSolver->writeMps("before");*/
    38823887                                        osiclp->getModelPtr()->setPerturbation(savePerturbation);
     3888                                        osiclp->getModelPtr()->setMoreSpecialOptions(saveOptions);
    38833889                                    }
    38843890#elif CBC_OTHER_SOLVER==1
     
    49534959                                   double obj;
    49544960                                   int status = computeCompleteSolution( babModel_, colNames, mipStart, &x[0], obj );
    4955                                    if (!status)
     4961                                   if (!status) {
    49564962                                     babModel_->setBestSolution( &x[0], static_cast<int>(x.size()), obj, false );
     4963                                     babModel_->setSolutionCount(1);
     4964                                   }
    49574965                                }
    49584966
     
    61716179                                if (biLinearProblem)
    61726180                                  babModel_->setSpecialOptions(babModel_->specialOptions() &(~(512|32768)));
     6181                                babModel_->setMoreSpecialOptions2(parameters_[whichParam(CBC_PARAM_INT_MOREMOREMIPOPTIONS, numberParameters_, parameters_)].intValue());
    61736182                                babModel_->branchAndBound(statistics);
    61746183                                if (truncateColumns<babModel_->solver()->getNumCols()) {
     
    94329441    delete [] numberExact;
    94339442}
     9443static void sortOnOther(int * column,
     9444                 const CoinBigIndex * rowStart,
     9445                 int * order,
     9446                 int * other,
     9447                 int nRow,
     9448                 int nInRow,
     9449                 int where)
     9450{
     9451     if (nRow < 2 || where >= nInRow)
     9452          return;
     9453     // do initial sort
     9454     int kRow;
     9455     int iRow;
     9456     for ( kRow = 0; kRow < nRow; kRow++) {
     9457          iRow = order[kRow];
     9458          other[kRow] = column[rowStart[iRow] + where];
     9459     }
     9460     CoinSort_2(other, other + nRow, order);
     9461     int first = 0;
     9462     iRow = order[0];
     9463     int firstC = column[rowStart[iRow] + where];
     9464     kRow = 1;
     9465     while (kRow < nRow) {
     9466          int lastC = 9999999;;
     9467          for (; kRow < nRow + 1; kRow++) {
     9468               if (kRow < nRow) {
     9469                    iRow = order[kRow];
     9470                    lastC = column[rowStart[iRow] + where];
     9471               } else {
     9472                    lastC = 9999999;
     9473               }
     9474               if (lastC > firstC)
     9475                    break;
     9476          }
     9477          // sort
     9478          sortOnOther(column, rowStart, order + first, other, kRow - first,
     9479                      nInRow, where + 1);
     9480          firstC = lastC;
     9481          first = kRow;
     9482     }
     9483}
    94349484static void statistics(ClpSimplex * originalModel, ClpSimplex * model)
    94359485{
    9436     int numberColumns = originalModel->numberColumns();
    9437     const char * integerInformation  = originalModel->integerInformation();
    9438     const double * columnLower = originalModel->columnLower();
    9439     const double * columnUpper = originalModel->columnUpper();
    9440     int numberIntegers = 0;
    9441     int numberBinary = 0;
    9442     int iRow, iColumn;
    9443     if (integerInformation) {
    9444         for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    9445             if (integerInformation[iColumn]) {
    9446                 if (columnUpper[iColumn] > columnLower[iColumn]) {
    9447                     numberIntegers++;
    9448                     if (columnUpper[iColumn] == 0.0 && columnLower[iColumn] == 1)
    9449                         numberBinary++;
    9450                 }
    9451             }
    9452         }
    9453     }
    9454     numberColumns = model->numberColumns();
    9455     int numberRows = model->numberRows();
    9456     columnLower = model->columnLower();
    9457     columnUpper = model->columnUpper();
    9458     const double * rowLower = model->rowLower();
    9459     const double * rowUpper = model->rowUpper();
    9460     const double * objective = model->objective();
    9461     CoinPackedMatrix * matrix = model->matrix();
    9462     CoinBigIndex numberElements = matrix->getNumElements();
    9463     const int * columnLength = matrix->getVectorLengths();
    9464     //const CoinBigIndex * columnStart = matrix->getVectorStarts();
    9465     const double * elementByColumn = matrix->getElements();
    9466     int * number = new int[numberRows+1];
    9467     memset(number, 0, (numberRows + 1)*sizeof(int));
    9468     int numberObjSingletons = 0;
    9469     /* cType
    9470        0 0/inf, 1 0/up, 2 lo/inf, 3 lo/up, 4 free, 5 fix, 6 -inf/0, 7 -inf/up,
    9471        8 0/1
    9472     */
    9473     int cType[9];
    9474     std::string cName[] = {"0.0->inf,", "0.0->up,", "lo->inf,", "lo->up,", "free,", "fixed,", "-inf->0.0,",
    9475                            "-inf->up,", "0.0->1.0"
    9476                           };
    9477     int nObjective = 0;
    9478     memset(cType, 0, sizeof(cType));
    9479     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    9480         int length = columnLength[iColumn];
    9481         if (length == 1 && objective[iColumn])
    9482             numberObjSingletons++;
    9483         number[length]++;
    9484         if (objective[iColumn])
    9485             nObjective++;
    9486         if (columnLower[iColumn] > -1.0e20) {
    9487             if (columnLower[iColumn] == 0.0) {
    9488                 if (columnUpper[iColumn] > 1.0e20)
    9489                     cType[0]++;
    9490                 else if (columnUpper[iColumn] == 1.0)
    9491                     cType[8]++;
    9492                 else if (columnUpper[iColumn] == 0.0)
    9493                     cType[5]++;
    9494                 else
    9495                     cType[1]++;
    9496             } else {
    9497                 if (columnUpper[iColumn] > 1.0e20)
    9498                     cType[2]++;
    9499                 else if (columnUpper[iColumn] == columnLower[iColumn])
    9500                     cType[5]++;
    9501                 else
    9502                     cType[3]++;
    9503             }
    9504         } else {
    9505             if (columnUpper[iColumn] > 1.0e20)
    9506                 cType[4]++;
    9507             else if (columnUpper[iColumn] == 0.0)
    9508                 cType[6]++;
    9509             else
    9510                 cType[7]++;
    9511         }
    9512     }
    9513     /* rType
    9514        0 E 0, 1 E 1, 2 E -1, 3 E other, 4 G 0, 5 G 1, 6 G other,
    9515        7 L 0,  8 L 1, 9 L other, 10 Range 0/1, 11 Range other, 12 free
    9516     */
    9517     int rType[13];
    9518     std::string rName[] = {"E 0.0,", "E 1.0,", "E -1.0,", "E other,", "G 0.0,", "G 1.0,", "G other,",
    9519                            "L 0.0,", "L 1.0,", "L other,", "Range 0.0->1.0,", "Range other,", "Free"
    9520                           };
    9521     memset(rType, 0, sizeof(rType));
    9522     for (iRow = 0; iRow < numberRows; iRow++) {
    9523         if (rowLower[iRow] > -1.0e20) {
    9524             if (rowLower[iRow] == 0.0) {
    9525                 if (rowUpper[iRow] > 1.0e20)
    9526                     rType[4]++;
    9527                 else if (rowUpper[iRow] == 1.0)
    9528                     rType[10]++;
    9529                 else if (rowUpper[iRow] == 0.0)
    9530                     rType[0]++;
    9531                 else
    9532                     rType[11]++;
    9533             } else if (rowLower[iRow] == 1.0) {
    9534                 if (rowUpper[iRow] > 1.0e20)
    9535                     rType[5]++;
    9536                 else if (rowUpper[iRow] == rowLower[iRow])
    9537                     rType[1]++;
    9538                 else
    9539                     rType[11]++;
    9540             } else if (rowLower[iRow] == -1.0) {
    9541                 if (rowUpper[iRow] > 1.0e20)
    9542                     rType[6]++;
    9543                 else if (rowUpper[iRow] == rowLower[iRow])
    9544                     rType[2]++;
    9545                 else
    9546                     rType[11]++;
    9547             } else {
    9548                 if (rowUpper[iRow] > 1.0e20)
    9549                     rType[6]++;
    9550                 else if (rowUpper[iRow] == rowLower[iRow])
    9551                     rType[3]++;
    9552                 else
    9553                     rType[11]++;
    9554             }
    9555         } else {
    9556             if (rowUpper[iRow] > 1.0e20)
    9557                 rType[12]++;
    9558             else if (rowUpper[iRow] == 0.0)
    9559                 rType[7]++;
    9560             else if (rowUpper[iRow] == 1.0)
    9561                 rType[8]++;
    9562             else
    9563                 rType[9]++;
    9564         }
    9565     }
    9566     // Basic statistics
    9567     printf("\n\nProblem has %d rows, %d columns (%d with objective) and %d elements\n",
    9568            numberRows, numberColumns, nObjective, numberElements);
    9569     if (number[0] + number[1]) {
    9570         printf("There are ");
    9571         if (numberObjSingletons)
    9572             printf("%d singletons with objective ", numberObjSingletons);
    9573         int numberNoObj = number[1] - numberObjSingletons;
    9574         if (numberNoObj)
    9575             printf("%d singletons with no objective ", numberNoObj);
    9576         if (number[0])
    9577             printf("** %d columns have no entries", number[0]);
    9578         printf("\n");
    9579     }
    9580     printf("Column breakdown:\n");
    9581     int k;
    9582     for (k = 0; k < static_cast<int> (sizeof(cType) / sizeof(int)); k++) {
    9583         printf("%d of type %s ", cType[k], cName[k].c_str());
    9584         if (((k + 1) % 3) == 0)
    9585             printf("\n");
    9586     }
    9587     if ((k % 3) != 0)
    9588         printf("\n");
    9589     printf("Row breakdown:\n");
    9590     for (k = 0; k < static_cast<int> (sizeof(rType) / sizeof(int)); k++) {
    9591         printf("%d of type %s ", rType[k], rName[k].c_str());
    9592         if (((k + 1) % 3) == 0)
    9593             printf("\n");
    9594     }
    9595     if ((k % 3) != 0)
    9596         printf("\n");
    9597     if (model->logLevel() < 2)
    9598         return ;
    9599     int kMax = model->logLevel() > 3 ? 1000000 : 10;
    9600     k = 0;
    9601     for (iRow = 1; iRow <= numberRows; iRow++) {
    9602         if (number[iRow]) {
    9603             k++;
    9604             printf("%d columns have %d entries\n", number[iRow], iRow);
    9605             if (k == kMax)
    9606                 break;
    9607         }
    9608     }
    9609     if (k < numberRows) {
    9610         int kk = k;
    9611         k = 0;
    9612         for (iRow = numberRows; iRow >= 1; iRow--) {
    9613             if (number[iRow]) {
    9614                 k++;
    9615                 if (k == kMax)
     9486     int numberColumns = originalModel->numberColumns();
     9487     const char * integerInformation  = originalModel->integerInformation();
     9488     const double * columnLower = originalModel->columnLower();
     9489     const double * columnUpper = originalModel->columnUpper();
     9490     int numberIntegers = 0;
     9491     int numberBinary = 0;
     9492     int iRow, iColumn;
     9493     if (integerInformation) {
     9494          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     9495               if (integerInformation[iColumn]) {
     9496                    if (columnUpper[iColumn] > columnLower[iColumn]) {
     9497                         numberIntegers++;
     9498                         if (columnLower[iColumn] == 0.0 && columnUpper[iColumn] == 1)
     9499                              numberBinary++;
     9500                    }
     9501               }
     9502          }
     9503          printf("Original problem has %d integers (%d of which binary)\n",
     9504                 numberIntegers,numberBinary);
     9505     }
     9506     numberColumns = model->numberColumns();
     9507     int numberRows = model->numberRows();
     9508     columnLower = model->columnLower();
     9509     columnUpper = model->columnUpper();
     9510     const double * rowLower = model->rowLower();
     9511     const double * rowUpper = model->rowUpper();
     9512     const double * objective = model->objective();
     9513     if (model->integerInformation()) {
     9514       const char * integerInformation  = model->integerInformation();
     9515       int numberIntegers = 0;
     9516       int numberBinary = 0;
     9517       double * obj = new double [numberColumns];
     9518       int * which = new int [numberColumns];
     9519       for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     9520         if (columnUpper[iColumn] > columnLower[iColumn]) {
     9521           if (integerInformation[iColumn]) {
     9522             numberIntegers++;
     9523             if (columnLower[iColumn] == 0.0 && columnUpper[iColumn] == 1)
     9524               numberBinary++;
     9525           }
     9526         }
     9527       }
     9528       if(numberColumns != originalModel->numberColumns())
     9529         printf("Presolved problem has %d integers (%d of which binary)\n",
     9530                numberIntegers,numberBinary);
     9531       for (int ifInt=0;ifInt<2;ifInt++) {
     9532         for (int ifAbs=0;ifAbs<2;ifAbs++) {
     9533           int numberSort=0;
     9534           int numberZero=0;
     9535           int numberDifferentObj=0;
     9536           for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     9537             if (columnUpper[iColumn] > columnLower[iColumn]) {
     9538               if (!ifInt||integerInformation[iColumn]) {
     9539                 obj[numberSort]=(ifAbs) ? fabs(objective[iColumn]) :
     9540                   objective[iColumn];
     9541                 which[numberSort++]=iColumn;
     9542                 if (!objective[iColumn])
     9543                   numberZero++;
     9544               }
     9545             }
     9546           }
     9547           CoinSort_2(obj,obj+numberSort,which);
     9548           double last=obj[0];
     9549           for (int jColumn = 1; jColumn < numberSort; jColumn++) {
     9550             if (fabs(obj[jColumn]-last)>1.0e-12) {
     9551               numberDifferentObj++;
     9552               last=obj[jColumn];
     9553             }
     9554           }
     9555           numberDifferentObj++;
     9556           printf("==== ");
     9557           if (ifInt)
     9558             printf("for integers ");
     9559           if (!ifAbs)
     9560             printf("%d zero objective ",numberZero);
     9561           else
     9562             printf("absolute objective values ");
     9563           printf("%d different\n",numberDifferentObj);
     9564           bool saveModel=false;
     9565           int target=model->logLevel();
     9566           if (target>10000) {
     9567             if (ifInt&&!ifAbs)
     9568               saveModel=true;
     9569             target-=10000;
     9570           }
     9571
     9572           if (target<=100)
     9573             target=12;
     9574           else
     9575             target-=100;
     9576           if (numberDifferentObj<target) {
     9577             int iLast=0;
     9578             double last=obj[0];
     9579             for (int jColumn = 1; jColumn < numberSort; jColumn++) {
     9580               if (fabs(obj[jColumn]-last)>1.0e-12) {
     9581                 printf("%d variables have objective of %g\n",
     9582                        jColumn-iLast,last);
     9583                 iLast=jColumn;
     9584                 last=obj[jColumn];
     9585               }
     9586             }
     9587             printf("%d variables have objective of %g\n",
     9588                    numberSort-iLast,last);
     9589             if (saveModel) {
     9590               int spaceNeeded=numberSort+numberDifferentObj;
     9591               int * columnAdd = new int[spaceNeeded+numberDifferentObj+1];
     9592               double * elementAdd = new double[spaceNeeded];
     9593               int * rowAdd = new int[2*numberDifferentObj+1];
     9594               int * newIsInteger = rowAdd+numberDifferentObj+1;
     9595               double * objectiveNew = new double[3*numberDifferentObj];
     9596               double * lowerNew = objectiveNew+numberDifferentObj;
     9597               double * upperNew = lowerNew+numberDifferentObj;
     9598               memset(columnAdd+spaceNeeded,0,
     9599                      (numberDifferentObj+1)*sizeof(int));
     9600               ClpSimplex tempModel=*model;
     9601               int iLast=0;
     9602               double last=obj[0];
     9603               numberDifferentObj=0;
     9604               int numberElements=0;
     9605               rowAdd[0]=0;
     9606               double * objective = tempModel.objective();
     9607               for (int jColumn = 1; jColumn < numberSort+1; jColumn++) {
     9608                 if (jColumn==numberSort||fabs(obj[jColumn]-last)>1.0e-12) {
     9609                   // not if just one
     9610                   if (jColumn-iLast>1) {
     9611                     bool allInteger=integerInformation!=NULL;
     9612                     int iColumn=which[iLast];
     9613                     objectiveNew[numberDifferentObj]=objective[iColumn];
     9614                     double lower=0.0;
     9615                     double upper=0.0;
     9616                     for (int kColumn=iLast;kColumn<jColumn;kColumn++) {
     9617                       iColumn=which[kColumn];
     9618                       objective[iColumn]=0.0;
     9619                       double lowerValue=columnLower[iColumn];
     9620                       double upperValue=columnUpper[iColumn];
     9621                       double elementValue=-1.0;
     9622                       if (objectiveNew[numberDifferentObj]*objective[iColumn]<0.0) {
     9623                         lowerValue=-columnUpper[iColumn];
     9624                         upperValue=-columnLower[iColumn];
     9625                         elementValue=1.0;
     9626                       }
     9627                       columnAdd[numberElements]=iColumn;
     9628                       elementAdd[numberElements++]=elementValue;
     9629                       if (integerInformation&&!integerInformation[iColumn])
     9630                         allInteger=false;
     9631                       if (lower!=-COIN_DBL_MAX) {
     9632                         if (lowerValue!=-COIN_DBL_MAX)
     9633                           lower += lowerValue;
     9634                         else
     9635                           lower=-COIN_DBL_MAX;
     9636                       }
     9637                       if (upper!=COIN_DBL_MAX) {
     9638                         if (upperValue!=COIN_DBL_MAX)
     9639                           upper += upperValue;
     9640                         else
     9641                           upper=COIN_DBL_MAX;
     9642                       }
     9643                     }
     9644                     columnAdd[numberElements]=numberColumns+numberDifferentObj;
     9645                     elementAdd[numberElements++]=1.0;
     9646                     newIsInteger[numberDifferentObj]= (allInteger) ? 1 : 0;
     9647                     lowerNew[numberDifferentObj]=lower;
     9648                     upperNew[numberDifferentObj]=upper;
     9649                     numberDifferentObj++;
     9650                     rowAdd[numberDifferentObj]=numberElements;
     9651                   }
     9652                   iLast=jColumn;
     9653                   last=obj[jColumn];
     9654                 }
     9655               }
     9656               // add columns
     9657               tempModel.addColumns(numberDifferentObj, lowerNew, upperNew,
     9658                                    objectiveNew,
     9659                                    columnAdd+spaceNeeded, NULL, NULL);
     9660               // add constraints and make integer if all integer in group
     9661               for (int iObj=0; iObj < numberDifferentObj; iObj++) {
     9662                 lowerNew[iObj]=0.0;
     9663                 upperNew[iObj]=0.0;
     9664                 if (newIsInteger[iObj])
     9665                   tempModel.setInteger(numberColumns+iObj);
     9666               }
     9667               tempModel.addRows(numberDifferentObj, lowerNew, upperNew,
     9668                                 rowAdd,columnAdd,elementAdd);
     9669               delete [] columnAdd;
     9670               delete [] elementAdd;
     9671               delete [] rowAdd;
     9672               delete [] objectiveNew;
     9673               // save
     9674               std::string tempName = model->problemName();
     9675               if (ifInt)
     9676                 tempName += "_int";
     9677               if (ifAbs)
     9678                 tempName += "_abs";
     9679               tempName += ".mps";
     9680               tempModel.writeMps(tempName.c_str());
     9681             }
     9682           }
     9683         }
     9684       }
     9685       delete [] which;
     9686       delete [] obj;
     9687       printf("===== end objective counts\n");
     9688     }
     9689     CoinPackedMatrix * matrix = model->matrix();
     9690     CoinBigIndex numberElements = matrix->getNumElements();
     9691     const int * columnLength = matrix->getVectorLengths();
     9692     //const CoinBigIndex * columnStart = matrix->getVectorStarts();
     9693     const double * elementByColumn = matrix->getElements();
     9694     int * number = new int[numberRows+1];
     9695     memset(number, 0, (numberRows + 1)*sizeof(int));
     9696     int numberObjSingletons = 0;
     9697     /* cType
     9698        0 0/inf, 1 0/up, 2 lo/inf, 3 lo/up, 4 free, 5 fix, 6 -inf/0, 7 -inf/up,
     9699        8 0/1
     9700     */
     9701     int cType[9];
     9702     std::string cName[] = {"0.0->inf,", "0.0->up,", "lo->inf,", "lo->up,", "free,", "fixed,", "-inf->0.0,",
     9703                            "-inf->up,", "0.0->1.0"
     9704                           };
     9705     int nObjective = 0;
     9706     memset(cType, 0, sizeof(cType));
     9707     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     9708          int length = columnLength[iColumn];
     9709          if (length == 1 && objective[iColumn])
     9710               numberObjSingletons++;
     9711          number[length]++;
     9712          if (objective[iColumn])
     9713               nObjective++;
     9714          if (columnLower[iColumn] > -1.0e20) {
     9715               if (columnLower[iColumn] == 0.0) {
     9716                    if (columnUpper[iColumn] > 1.0e20)
     9717                         cType[0]++;
     9718                    else if (columnUpper[iColumn] == 1.0)
     9719                         cType[8]++;
     9720                    else if (columnUpper[iColumn] == 0.0)
     9721                         cType[5]++;
     9722                    else
     9723                         cType[1]++;
     9724               } else {
     9725                    if (columnUpper[iColumn] > 1.0e20)
     9726                         cType[2]++;
     9727                    else if (columnUpper[iColumn] == columnLower[iColumn])
     9728                         cType[5]++;
     9729                    else
     9730                         cType[3]++;
     9731               }
     9732          } else {
     9733               if (columnUpper[iColumn] > 1.0e20)
     9734                    cType[4]++;
     9735               else if (columnUpper[iColumn] == 0.0)
     9736                    cType[6]++;
     9737               else
     9738                    cType[7]++;
     9739          }
     9740     }
     9741     /* rType
     9742        0 E 0, 1 E 1, 2 E -1, 3 E other, 4 G 0, 5 G 1, 6 G other,
     9743        7 L 0,  8 L 1, 9 L other, 10 Range 0/1, 11 Range other, 12 free
     9744     */
     9745     int rType[13];
     9746     std::string rName[] = {"E 0.0,", "E 1.0,", "E -1.0,", "E other,", "G 0.0,", "G 1.0,", "G other,",
     9747                            "L 0.0,", "L 1.0,", "L other,", "Range 0.0->1.0,", "Range other,", "Free"
     9748                           };
     9749     memset(rType, 0, sizeof(rType));
     9750     for (iRow = 0; iRow < numberRows; iRow++) {
     9751          if (rowLower[iRow] > -1.0e20) {
     9752               if (rowLower[iRow] == 0.0) {
     9753                    if (rowUpper[iRow] > 1.0e20)
     9754                         rType[4]++;
     9755                    else if (rowUpper[iRow] == 1.0)
     9756                         rType[10]++;
     9757                    else if (rowUpper[iRow] == 0.0)
     9758                         rType[0]++;
     9759                    else
     9760                         rType[11]++;
     9761               } else if (rowLower[iRow] == 1.0) {
     9762                    if (rowUpper[iRow] > 1.0e20)
     9763                         rType[5]++;
     9764                    else if (rowUpper[iRow] == rowLower[iRow])
     9765                         rType[1]++;
     9766                    else
     9767                         rType[11]++;
     9768               } else if (rowLower[iRow] == -1.0) {
     9769                    if (rowUpper[iRow] > 1.0e20)
     9770                         rType[6]++;
     9771                    else if (rowUpper[iRow] == rowLower[iRow])
     9772                         rType[2]++;
     9773                    else
     9774                         rType[11]++;
     9775               } else {
     9776                    if (rowUpper[iRow] > 1.0e20)
     9777                         rType[6]++;
     9778                    else if (rowUpper[iRow] == rowLower[iRow])
     9779                         rType[3]++;
     9780                    else
     9781                         rType[11]++;
     9782               }
     9783          } else {
     9784               if (rowUpper[iRow] > 1.0e20)
     9785                    rType[12]++;
     9786               else if (rowUpper[iRow] == 0.0)
     9787                    rType[7]++;
     9788               else if (rowUpper[iRow] == 1.0)
     9789                    rType[8]++;
     9790               else
     9791                    rType[9]++;
     9792          }
     9793     }
     9794     // Basic statistics
     9795     printf("\n\nProblem has %d rows, %d columns (%d with objective) and %d elements\n",
     9796            numberRows, numberColumns, nObjective, numberElements);
     9797     if (number[0] + number[1]) {
     9798          printf("There are ");
     9799          if (numberObjSingletons)
     9800               printf("%d singletons with objective ", numberObjSingletons);
     9801          int numberNoObj = number[1] - numberObjSingletons;
     9802          if (numberNoObj)
     9803               printf("%d singletons with no objective ", numberNoObj);
     9804          if (number[0])
     9805               printf("** %d columns have no entries", number[0]);
     9806          printf("\n");
     9807     }
     9808     printf("Column breakdown:\n");
     9809     int k;
     9810     for (k = 0; k < static_cast<int> (sizeof(cType) / sizeof(int)); k++) {
     9811          printf("%d of type %s ", cType[k], cName[k].c_str());
     9812          if (((k + 1) % 3) == 0)
     9813               printf("\n");
     9814     }
     9815     if ((k % 3) != 0)
     9816          printf("\n");
     9817     printf("Row breakdown:\n");
     9818     for (k = 0; k < static_cast<int> (sizeof(rType) / sizeof(int)); k++) {
     9819          printf("%d of type %s ", rType[k], rName[k].c_str());
     9820          if (((k + 1) % 3) == 0)
     9821               printf("\n");
     9822     }
     9823     if ((k % 3) != 0)
     9824          printf("\n");
     9825     //#define SYM
     9826#ifndef SYM
     9827     if (model->logLevel() < 2)
     9828          return ;
     9829#endif
     9830     int kMax = model->logLevel() > 3 ? 1000000 : 10;
     9831     k = 0;
     9832     for (iRow = 1; iRow <= numberRows; iRow++) {
     9833          if (number[iRow]) {
     9834               k++;
     9835               printf("%d columns have %d entries\n", number[iRow], iRow);
     9836               if (k == kMax)
    96169837                    break;
    9617             }
    9618         }
    9619         if (k > kk) {
    9620             printf("\n    .........\n\n");
    9621             iRow = k;
    9622             k = 0;
    9623             for (; iRow < numberRows; iRow++) {
    9624                 if (number[iRow]) {
     9838          }
     9839     }
     9840     if (k < numberRows) {
     9841          int kk = k;
     9842          k = 0;
     9843          for (iRow = numberRows; iRow >= 1; iRow--) {
     9844               if (number[iRow]) {
    96259845                    k++;
    9626                     printf("%d columns have %d entries\n", number[iRow], iRow);
    96279846                    if (k == kMax)
    9628                         break;
    9629                 }
    9630             }
    9631         }
    9632     }
    9633     delete [] number;
    9634     printf("\n\n");
    9635     // get row copy
    9636     CoinPackedMatrix rowCopy = *matrix;
    9637     rowCopy.reverseOrdering();
    9638     //const int * column = rowCopy.getIndices();
    9639     const int * rowLength = rowCopy.getVectorLengths();
    9640     //const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
    9641     //const double * element = rowCopy.getElements();
    9642     number = new int[numberColumns+1];
    9643     memset(number, 0, (numberColumns + 1)*sizeof(int));
    9644     for (iRow = 0; iRow < numberRows; iRow++) {
    9645         int length = rowLength[iRow];
    9646         number[length]++;
    9647     }
    9648     if (number[0])
    9649         printf("** %d rows have no entries\n", number[0]);
    9650     k = 0;
    9651     for (iColumn = 1; iColumn <= numberColumns; iColumn++) {
    9652         if (number[iColumn]) {
    9653             k++;
    9654             printf("%d rows have %d entries\n", number[iColumn], iColumn);
    9655             if (k == kMax)
    9656                 break;
    9657         }
    9658     }
    9659     if (k < numberColumns) {
    9660         int kk = k;
    9661         k = 0;
    9662         for (iColumn = numberColumns; iColumn >= 1; iColumn--) {
    9663             if (number[iColumn]) {
    9664                 k++;
    9665                 if (k == kMax)
     9847                         break;
     9848               }
     9849          }
     9850          if (k > kk) {
     9851               printf("\n    .........\n\n");
     9852               iRow = k;
     9853               k = 0;
     9854               for (; iRow < numberRows; iRow++) {
     9855                    if (number[iRow]) {
     9856                         k++;
     9857                         printf("%d columns have %d entries\n", number[iRow], iRow);
     9858                         if (k == kMax)
     9859                              break;
     9860                    }
     9861               }
     9862          }
     9863     }
     9864     delete [] number;
     9865     printf("\n\n");
     9866     if (model->logLevel() == 63
     9867#ifdef SYM
     9868               || true
     9869#endif
     9870        ) {
     9871          // get column copy
     9872          CoinPackedMatrix columnCopy = *matrix;
     9873          const int * columnLength = columnCopy.getVectorLengths();
     9874          number = new int[numberRows+1];
     9875          memset(number, 0, (numberRows + 1)*sizeof(int));
     9876          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     9877               int length = columnLength[iColumn];
     9878               number[length]++;
     9879          }
     9880          k = 0;
     9881          for (iRow = 1; iRow <= numberRows; iRow++) {
     9882               if (number[iRow]) {
     9883                    k++;
     9884               }
     9885          }
     9886          int * row = columnCopy.getMutableIndices();
     9887          const CoinBigIndex * columnStart = columnCopy.getVectorStarts();
     9888          double * element = columnCopy.getMutableElements();
     9889          int * order = new int[numberColumns];
     9890          int * other = new int[numberColumns];
     9891          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     9892               int length = columnLength[iColumn];
     9893               order[iColumn] = iColumn;
     9894               other[iColumn] = length;
     9895               CoinBigIndex start = columnStart[iColumn];
     9896               CoinSort_2(row + start, row + start + length, element + start);
     9897          }
     9898          CoinSort_2(other, other + numberColumns, order);
     9899          int jColumn = number[0] + number[1];
     9900          for (iRow = 2; iRow <= numberRows; iRow++) {
     9901               if (number[iRow]) {
     9902                    printf("XX %d columns have %d entries\n", number[iRow], iRow);
     9903                    int kColumn = jColumn + number[iRow];
     9904                    sortOnOther(row, columnStart,
     9905                                order + jColumn, other, number[iRow], iRow, 0);
     9906                    // Now print etc
     9907                    if (iRow < 500000) {
     9908                         for (int lColumn = jColumn; lColumn < kColumn; lColumn++) {
     9909                              iColumn = order[lColumn];
     9910                              CoinBigIndex start = columnStart[iColumn];
     9911                              if (model->logLevel() == 63) {
     9912                                   printf("column %d %g <= ", iColumn, columnLower[iColumn]);
     9913                                   for (CoinBigIndex i = start; i < start + iRow; i++)
     9914                                        printf("( %d, %g) ", row[i], element[i]);
     9915                                   printf("<= %g\n", columnUpper[iColumn]);
     9916                              }
     9917                         }
     9918                    }
     9919                    jColumn = kColumn;
     9920               }
     9921          }
     9922          delete [] order;
     9923          delete [] other;
     9924          delete [] number;
     9925     }
     9926     // get row copy
     9927     CoinPackedMatrix rowCopy = *matrix;
     9928     rowCopy.reverseOrdering();
     9929     const int * rowLength = rowCopy.getVectorLengths();
     9930     number = new int[numberColumns+1];
     9931     memset(number, 0, (numberColumns + 1)*sizeof(int));
     9932     if (model->logLevel() > 3) {
     9933       // get column copy
     9934       CoinPackedMatrix columnCopy = *matrix;
     9935       const int * columnLength = columnCopy.getVectorLengths();
     9936       const int * row = columnCopy.getIndices();
     9937       const CoinBigIndex * columnStart = columnCopy.getVectorStarts();
     9938       const double * element = columnCopy.getElements();
     9939       const double * elementByRow = rowCopy.getElements();
     9940       const int * rowStart = rowCopy.getVectorStarts();
     9941       const int * column = rowCopy.getIndices();
     9942       int nPossibleZeroCost=0;
     9943       int nPossibleNonzeroCost=0;
     9944       for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     9945         int length = columnLength[iColumn];
     9946         if (columnLower[iColumn]<-1.0e30&&columnUpper[iColumn]>1.0e30) {
     9947           if (length==1) {
     9948             printf("Singleton free %d - cost %g\n",iColumn,objective[iColumn]);
     9949           } else if (length==2) {
     9950             int iRow0=row[columnStart[iColumn]];
     9951             int iRow1=row[columnStart[iColumn]+1];
     9952             double element0=element[columnStart[iColumn]];
     9953             double element1=element[columnStart[iColumn]+1];
     9954             int n0=rowLength[iRow0];
     9955             int n1=rowLength[iRow1];
     9956             printf("Doubleton free %d - cost %g - %g in %srow with %d entries and %g in %srow with %d entries\n",
     9957                    iColumn,objective[iColumn],element0,(rowLower[iRow0]==rowUpper[iRow0]) ? "==" : "",n0,
     9958                    element1,(rowLower[iRow1]==rowUpper[iRow1]) ? "==" : "",n1);
     9959
     9960           }
     9961         }
     9962         if (length==1) {
     9963           int iRow=row[columnStart[iColumn]];
     9964           double value=COIN_DBL_MAX;
     9965           for (int i=rowStart[iRow];i<rowStart[iRow]+rowLength[iRow];i++) {
     9966             int jColumn=column[i];
     9967             if (jColumn!=iColumn) {
     9968               if (value!=elementByRow[i]) {
     9969                 if (value==COIN_DBL_MAX) {
     9970                   value=elementByRow[i];
     9971                 } else {
     9972                   value = -COIN_DBL_MAX;
     9973                   break;
     9974                 }
     9975               }
     9976             }
     9977           }
     9978           if (!objective[iColumn]) {
     9979             if (model->logLevel() > 4)
     9980             printf("Singleton %d with no objective in row with %d elements - rhs %g,%g\n",iColumn,rowLength[iRow],rowLower[iRow],rowUpper[iRow]);
     9981             nPossibleZeroCost++;
     9982           } else if (value!=-COIN_DBL_MAX) {
     9983             if (model->logLevel() > 4)
     9984               printf("Singleton %d (%s) with objective in row %d (%s) with %d equal elements - rhs %g,%g\n",iColumn,model->getColumnName(iColumn).c_str(),
     9985                      iRow,model->getRowName(iRow).c_str(),
     9986                      rowLength[iRow],rowLower[iRow],rowUpper[iRow]);
     9987             nPossibleNonzeroCost++;
     9988           }
     9989         }
     9990       }
     9991       if (nPossibleZeroCost||nPossibleNonzeroCost)
     9992         printf("%d singletons with zero cost, %d with valid cost\n",
     9993                nPossibleZeroCost,nPossibleNonzeroCost);
     9994       // look for DW
     9995       int * blockStart = new int [2*(numberRows+numberColumns)+1+numberRows];
     9996       int * columnBlock = blockStart+numberRows;
     9997       int * nextColumn = columnBlock+numberColumns;
     9998       int * blockCount = nextColumn+numberColumns;
     9999       int * blockEls = blockCount+numberRows+1;
     10000       int direction[2]={-1,1};
     10001       int bestBreak=-1;
     10002       double bestValue=0.0;
     10003       int iPass=0;
     10004       int halfway=(numberRows+1)/2;
     10005       int firstMaster=-1;
     10006       int lastMaster=-2;
     10007       while (iPass<2) {
     10008         int increment=direction[iPass];
     10009         int start= increment>0 ? 0 : numberRows-1;
     10010         int stop=increment>0 ? numberRows : -1;
     10011         int numberBlocks=0;
     10012         int thisBestBreak=-1;
     10013         double thisBestValue=COIN_DBL_MAX;
     10014         int numberRowsDone=0;
     10015         int numberMarkedColumns=0;
     10016         int maximumBlockSize=0;
     10017         for (int i=0;i<numberRows+2*numberColumns;i++)
     10018           blockStart[i]=-1;
     10019         for (int i=0;i<numberRows+1;i++)
     10020           blockCount[i]=0;
     10021         for (int iRow=start;iRow!=stop;iRow+=increment) {
     10022           int iBlock = -1;
     10023           for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
     10024             int iColumn=column[j];
     10025             int whichColumnBlock=columnBlock[iColumn];
     10026             if (whichColumnBlock>=0) {
     10027               // column marked
     10028               if (iBlock<0) {
     10029                 // put row in that block
     10030                 iBlock=whichColumnBlock;
     10031               } else if (iBlock!=whichColumnBlock) {
     10032                 // merge
     10033                 blockCount[iBlock]+=blockCount[whichColumnBlock];
     10034                 blockCount[whichColumnBlock]=0;
     10035                 int jColumn=blockStart[whichColumnBlock];
     10036                 while (jColumn>=0) {
     10037                   columnBlock[jColumn]=iBlock;
     10038                   iColumn=jColumn;
     10039                   jColumn=nextColumn[jColumn];
     10040                 }
     10041                 nextColumn[iColumn]=blockStart[iBlock];
     10042                 blockStart[iBlock]=blockStart[whichColumnBlock];
     10043                 blockStart[whichColumnBlock]=-1;
     10044               }
     10045             }
     10046           }
     10047           int n=numberMarkedColumns;
     10048           if (iBlock<0) {
     10049             //new block
     10050             if (rowLength[iRow]) {
     10051               numberBlocks++;
     10052               iBlock=numberBlocks;
     10053               int jColumn=column[rowStart[iRow]];
     10054               columnBlock[jColumn]=iBlock;
     10055               blockStart[iBlock]=jColumn;
     10056               numberMarkedColumns++;
     10057               for (CoinBigIndex j=rowStart[iRow]+1;j<rowStart[iRow]+rowLength[iRow];j++) {
     10058                 int iColumn=column[j];
     10059                 columnBlock[iColumn]=iBlock;
     10060                 numberMarkedColumns++;
     10061                 nextColumn[jColumn]=iColumn;
     10062                 jColumn=iColumn;
     10063               }
     10064               blockCount[iBlock]=numberMarkedColumns-n;
     10065             } else {
     10066               // empty
     10067               iBlock=numberRows;
     10068             }
     10069           } else {
     10070             // put in existing block
     10071             int jColumn=blockStart[iBlock];
     10072             for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
     10073               int iColumn=column[j];
     10074               assert (columnBlock[iColumn]<0||columnBlock[iColumn]==iBlock);
     10075               if (columnBlock[iColumn]<0) {
     10076                 columnBlock[iColumn]=iBlock;
     10077                 numberMarkedColumns++;
     10078                 nextColumn[iColumn]=jColumn;
     10079                 jColumn=iColumn;
     10080               }
     10081             }
     10082             blockStart[iBlock]=jColumn;
     10083             blockCount[iBlock]+=numberMarkedColumns-n;
     10084           }
     10085           maximumBlockSize=CoinMax(maximumBlockSize,blockCount[iBlock]);
     10086           numberRowsDone++;
     10087           if (thisBestValue*numberRowsDone > maximumBlockSize&&numberRowsDone>halfway) {
     10088             thisBestBreak=iRow;
     10089             thisBestValue=static_cast<double>(maximumBlockSize)/
     10090               static_cast<double>(numberRowsDone);
     10091           }
     10092         }
     10093         if (thisBestBreak==stop)
     10094           thisBestValue=COIN_DBL_MAX;
     10095         iPass++;
     10096         if (iPass==1) {
     10097           bestBreak=thisBestBreak;
     10098           bestValue=thisBestValue;
     10099         } else {
     10100           if (bestValue<thisBestValue) {
     10101             firstMaster=0;
     10102             lastMaster=bestBreak;
     10103           } else {
     10104             firstMaster=thisBestBreak; // ? +1
     10105             lastMaster=numberRows;
     10106           }
     10107         }
     10108       }
     10109       if (firstMaster<lastMaster) {
     10110         printf("%d master rows %d <= < %d\n",lastMaster-firstMaster,
     10111                firstMaster,lastMaster);
     10112         for (int i=0;i<numberRows+2*numberColumns;i++)
     10113           blockStart[i]=-1;
     10114         for (int i=firstMaster;i<lastMaster;i++)
     10115           blockStart[i]=-2;
     10116         int firstRow=0;
     10117         int numberBlocks=-1;
     10118         while (true) {
     10119           for (;firstRow<numberRows;firstRow++) {
     10120             if (blockStart[firstRow]==-1)
     10121               break;
     10122           }
     10123           if (firstRow==numberRows)
     10124             break;
     10125           int nRows=0;
     10126           numberBlocks++;
     10127           int numberStack=1;
     10128           blockCount[0] = firstRow;
     10129           while (numberStack) {
     10130             int iRow=blockCount[--numberStack];
     10131             for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
     10132               int iColumn=column[j];
     10133               int iBlock=columnBlock[iColumn];
     10134               if (iBlock<0) {
     10135                 columnBlock[iColumn]=numberBlocks;
     10136                 for (CoinBigIndex k=columnStart[iColumn];
     10137                      k<columnStart[iColumn]+columnLength[iColumn];k++) {
     10138                   int jRow=row[k];
     10139                   int rowBlock=blockStart[jRow];
     10140                   if (rowBlock==-1) {
     10141                     nRows++;
     10142                     blockStart[jRow]=numberBlocks;
     10143                     blockCount[numberStack++]=jRow;
     10144                   }
     10145                 }
     10146               }
     10147             }
     10148           }
     10149           if (!nRows) {
     10150             // empty!!
     10151             numberBlocks--;
     10152           }
     10153           firstRow++;
     10154         }
     10155         // adjust
     10156         numberBlocks++;
     10157         for (int i=0;i<numberBlocks;i++) {
     10158           blockCount[i]=0;
     10159           nextColumn[i]=0;
     10160         }
     10161         int numberEmpty=0;
     10162         int numberMaster=0;
     10163         memset(blockEls,0,numberBlocks*sizeof(int));
     10164         for (int iRow = 0; iRow < numberRows; iRow++) {
     10165           int iBlock=blockStart[iRow];
     10166           if (iBlock>=0) {
     10167             blockCount[iBlock]++;
     10168             blockEls[iBlock]+=rowLength[iRow];
     10169           } else {
     10170             if (iBlock==-2)
     10171               numberMaster++;
     10172             else
     10173               numberEmpty++;
     10174           }
     10175         }
     10176         int numberEmptyColumns=0;
     10177         int numberMasterColumns=0;
     10178         for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     10179           int iBlock=columnBlock[iColumn];
     10180           if (iBlock>=0) {
     10181             nextColumn[iBlock]++;
     10182           } else {
     10183             if (columnLength[iColumn])
     10184               numberMasterColumns++;
     10185             else
     10186               numberEmptyColumns++;
     10187           }
     10188         }
     10189         int largestRows=0;
     10190         int largestColumns=0;
     10191         for (int i=0;i<numberBlocks;i++) {
     10192           if (blockCount[i]+nextColumn[i]>largestRows+largestColumns) {
     10193             largestRows=blockCount[i];
     10194             largestColumns=nextColumn[i];
     10195           }
     10196         }
     10197         bool useful=true;
     10198         if (numberMaster>halfway||largestRows*3>numberRows)
     10199           useful=false;
     10200         printf("%s %d blocks (largest %d,%d), %d master rows (%d empty) out of %d, %d master columns (%d empty) out of %d\n",
     10201                useful ? "**Useful" : "NoGood",
     10202                numberBlocks,largestRows,largestColumns,numberMaster,numberEmpty,numberRows,
     10203                numberMasterColumns,numberEmptyColumns,numberColumns);
     10204         FILE * fp=NULL;
     10205         bool justIntegers=true;
     10206         bool oneFile=true;
     10207         int logLevel = model->logLevel();
     10208         if (logLevel>19) {
     10209           logLevel-=2;
     10210           oneFile=true;
     10211           fp = fopen("fake.bnd","w");
     10212         }
     10213         if (logLevel==19)
     10214           justIntegers=false;
     10215         for (int i=0;i<numberBlocks;i++)
     10216           printf("Block %d has %d rows and %d columns (%d elements)\n",
     10217                  i,blockCount[i],nextColumn[i],blockEls[i]);
     10218         if (logLevel >= 17 && logLevel <= 21) {
     10219           int * whichRows=new int[numberRows+numberColumns];
     10220           int * whichColumns=whichRows+numberRows;
     10221           char name[20];
     10222           for (int iBlock=0;iBlock<numberBlocks;iBlock++) {
     10223             sprintf(name,"block%d.mps",iBlock);
     10224             int nRows=0;
     10225             for (int iRow=0;iRow<numberRows;iRow++) {
     10226               if (blockStart[iRow]==iBlock)
     10227                 whichRows[nRows++]=iRow;
     10228             }
     10229             int nColumns=0;
     10230             for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     10231               if (columnBlock[iColumn]==iBlock)
     10232                 whichColumns[nColumns++]=iColumn;
     10233             }
     10234             ClpSimplex subset(model,nRows,whichRows,nColumns,whichColumns);
     10235             for (int jRow=0;jRow<nRows;jRow++) {
     10236               int iRow = whichRows[jRow];
     10237               std::string name = model->getRowName(iRow);
     10238               subset.setRowName(jRow,name);
     10239             }
     10240             int nInteger=0;
     10241             for (int jColumn=0;jColumn<nColumns;jColumn++) {
     10242               int iColumn = whichColumns[jColumn];
     10243               if (model->isInteger(iColumn)) {
     10244                 subset.setInteger(jColumn);
     10245                 nInteger++;
     10246               }
     10247               std::string name = model->getColumnName(iColumn);
     10248               subset.setColumnName(jColumn,name);
     10249             }
     10250             if (logLevel == 17) {
     10251               subset.writeMps(name,0,1);
     10252             } else if (nInteger) {
     10253               OsiClpSolverInterface subset2(&subset);
     10254               CbcModel smallModel(subset2);
     10255               smallModel.branchAndBound();
     10256               const double * solution = smallModel.bestSolution();
     10257               if (solution) {
     10258                 if (!oneFile) {
     10259                   sprintf(name,"block%d.bnd",iBlock);
     10260                   fp = fopen(name,"w");
     10261                   assert (fp);
     10262                 }
     10263                 fprintf(fp,"BBB objective %g for block %d\n",
     10264                         smallModel.getObjValue(),iBlock);
     10265                 for (int jColumn=0;jColumn<nColumns;jColumn++) {
     10266                   if (subset.isInteger(jColumn)||
     10267                       !justIntegers)
     10268                     fprintf(fp," FX BOUND1    %.8s  %g\n",
     10269                             subset.getColumnName(jColumn).c_str(),
     10270                             solution[jColumn]);
     10271                 }
     10272                 if (!oneFile)
     10273                   fclose(fp);
     10274               } else {
     10275                 printf("***** Problem is infeasible\n");
     10276                 abort();
     10277               }
     10278             }
     10279           }
     10280           if (oneFile)
     10281             fclose(fp);
     10282           delete [] whichRows;
     10283         }
     10284       }
     10285       delete [] blockStart;
     10286     }
     10287     for (iRow = 0; iRow < numberRows; iRow++) {
     10288          int length = rowLength[iRow];
     10289          number[length]++;
     10290     }
     10291     if (number[0])
     10292          printf("** %d rows have no entries\n", number[0]);
     10293     k = 0;
     10294     for (iColumn = 1; iColumn <= numberColumns; iColumn++) {
     10295          if (number[iColumn]) {
     10296               k++;
     10297               printf("%d rows have %d entries\n", number[iColumn], iColumn);
     10298               if (k == kMax)
    966610299                    break;
    9667             }
    9668         }
    9669         if (k > kk) {
    9670             printf("\n    .........\n\n");
    9671             iColumn = k;
    9672             k = 0;
    9673             for (; iColumn < numberColumns; iColumn++) {
    9674                 if (number[iColumn]) {
     10300          }
     10301     }
     10302     if (k < numberColumns) {
     10303          int kk = k;
     10304          k = 0;
     10305          for (iColumn = numberColumns; iColumn >= 1; iColumn--) {
     10306               if (number[iColumn]) {
    967510307                    k++;
    9676                     printf("%d rows have %d entries\n", number[iColumn], iColumn);
    967710308                    if (k == kMax)
    9678                         break;
    9679                 }
    9680             }
    9681         }
    9682     }
    9683     delete [] number;
    9684     // Now do breakdown of ranges
    9685     breakdown("Elements", numberElements, elementByColumn);
    9686     breakdown("RowLower", numberRows, rowLower);
    9687     breakdown("RowUpper", numberRows, rowUpper);
    9688     breakdown("ColumnLower", numberColumns, columnLower);
    9689     breakdown("ColumnUpper", numberColumns, columnUpper);
    9690     breakdown("Objective", numberColumns, objective);
     10309                         break;
     10310               }
     10311          }
     10312          if (k > kk) {
     10313               printf("\n    .........\n\n");
     10314               iColumn = k;
     10315               k = 0;
     10316               for (; iColumn < numberColumns; iColumn++) {
     10317                    if (number[iColumn]) {
     10318                         k++;
     10319                         printf("%d rows have %d entries\n", number[iColumn], iColumn);
     10320                         if (k == kMax)
     10321                              break;
     10322                    }
     10323               }
     10324          }
     10325     }
     10326     if (model->logLevel() == 63
     10327#ifdef SYM
     10328               || true
     10329#endif
     10330        ) {
     10331          int * column = rowCopy.getMutableIndices();
     10332          const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
     10333          double * element = rowCopy.getMutableElements();
     10334          int * order = new int[numberRows];
     10335          int * other = new int[numberRows];
     10336          for (iRow = 0; iRow < numberRows; iRow++) {
     10337               int length = rowLength[iRow];
     10338               order[iRow] = iRow;
     10339               other[iRow] = length;
     10340               CoinBigIndex start = rowStart[iRow];
     10341               CoinSort_2(column + start, column + start + length, element + start);
     10342          }
     10343          CoinSort_2(other, other + numberRows, order);
     10344          int jRow = number[0] + number[1];
     10345          double * weight = new double[numberRows];
     10346          double * randomColumn = new double[numberColumns+1];
     10347          double * randomRow = new double [numberRows+1];
     10348          int * sortRow = new int [numberRows];
     10349          int * possibleRow = new int [numberRows];
     10350          int * backRow = new int [numberRows];
     10351          int * stackRow = new int [numberRows];
     10352          int * sortColumn = new int [numberColumns];
     10353          int * possibleColumn = new int [numberColumns];
     10354          int * backColumn = new int [numberColumns];
     10355          int * backColumn2 = new int [numberColumns];
     10356          int * mapRow = new int [numberRows];
     10357          int * mapColumn = new int [numberColumns];
     10358          int * stackColumn = new int [numberColumns];
     10359          double randomLower = CoinDrand48();
     10360          double randomUpper = CoinDrand48();
     10361          double randomInteger = CoinDrand48();
     10362          int * startAdd = new int[numberRows+1];
     10363          int * columnAdd = new int [2*numberElements];
     10364          double * elementAdd = new double[2*numberElements];
     10365          int nAddRows = 0;
     10366          startAdd[0] = 0;
     10367          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     10368               randomColumn[iColumn] = CoinDrand48();
     10369               backColumn2[iColumn] = -1;
     10370          }
     10371          for (iColumn = 2; iColumn <= numberColumns; iColumn++) {
     10372               if (number[iColumn]) {
     10373                    printf("XX %d rows have %d entries\n", number[iColumn], iColumn);
     10374                    int kRow = jRow + number[iColumn];
     10375                    sortOnOther(column, rowStart,
     10376                                order + jRow, other, number[iColumn], iColumn, 0);
     10377                    // Now print etc
     10378                    if (iColumn < 500000) {
     10379                         int nLook = 0;
     10380                         for (int lRow = jRow; lRow < kRow; lRow++) {
     10381                              iRow = order[lRow];
     10382                              CoinBigIndex start = rowStart[iRow];
     10383                              if (model->logLevel() == 63) {
     10384                                   printf("row %d %g <= ", iRow, rowLower[iRow]);
     10385                                   for (CoinBigIndex i = start; i < start + iColumn; i++)
     10386                                        printf("( %d, %g) ", column[i], element[i]);
     10387                                   printf("<= %g\n", rowUpper[iRow]);
     10388                              }
     10389                              int first = column[start];
     10390                              double sum = 0.0;
     10391                              for (CoinBigIndex i = start; i < start + iColumn; i++) {
     10392                                   int jColumn = column[i];
     10393                                   double value = element[i];
     10394                                   jColumn -= first;
     10395                                   assert (jColumn >= 0);
     10396                                   sum += value * randomColumn[jColumn];
     10397                              }
     10398                              if (rowLower[iRow] > -1.0e30 && rowLower[iRow])
     10399                                   sum += rowLower[iRow] * randomLower;
     10400                              else if (!rowLower[iRow])
     10401                                   sum += 1.234567e-7 * randomLower;
     10402                              if (rowUpper[iRow] < 1.0e30 && rowUpper[iRow])
     10403                                   sum += rowUpper[iRow] * randomUpper;
     10404                              else if (!rowUpper[iRow])
     10405                                   sum += 1.234567e-7 * randomUpper;
     10406                              sortRow[nLook] = iRow;
     10407                              randomRow[nLook++] = sum;
     10408                              // best way is to number unique elements and bounds and use
     10409                              if (fabs(sum) > 1.0e4)
     10410                                   sum *= 1.0e-6;
     10411                              weight[iRow] = sum;
     10412                         }
     10413                         assert (nLook <= numberRows);
     10414                         CoinSort_2(randomRow, randomRow + nLook, sortRow);
     10415                         randomRow[nLook] = COIN_DBL_MAX;
     10416                         double last = -COIN_DBL_MAX;
     10417                         int iLast = -1;
     10418                         for (int iLook = 0; iLook < nLook + 1; iLook++) {
     10419                              if (randomRow[iLook] > last) {
     10420                                   if (iLast >= 0) {
     10421                                        int n = iLook - iLast;
     10422                                        if (n > 1) {
     10423                                             //printf("%d rows possible?\n",n);
     10424                                        }
     10425                                   }
     10426                                   iLast = iLook;
     10427                                   last = randomRow[iLook];
     10428                              }
     10429                         }
     10430                    }
     10431                    jRow = kRow;
     10432               }
     10433          }
     10434          CoinPackedMatrix columnCopy = *matrix;
     10435          const int * columnLength = columnCopy.getVectorLengths();
     10436          const int * row = columnCopy.getIndices();
     10437          const CoinBigIndex * columnStart = columnCopy.getVectorStarts();
     10438          const double * elementByColumn = columnCopy.getElements();
     10439          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     10440               int length = columnLength[iColumn];
     10441               CoinBigIndex start = columnStart[iColumn];
     10442               double sum = objective[iColumn];
     10443               if (columnLower[iColumn] > -1.0e30 && columnLower[iColumn])
     10444                    sum += columnLower[iColumn] * randomLower;
     10445               else if (!columnLower[iColumn])
     10446                    sum += 1.234567e-7 * randomLower;
     10447               if (columnUpper[iColumn] < 1.0e30 && columnUpper[iColumn])
     10448                    sum += columnUpper[iColumn] * randomUpper;
     10449               else if (!columnUpper[iColumn])
     10450                    sum += 1.234567e-7 * randomUpper;
     10451               if (model->isInteger(iColumn))
     10452                    sum += 9.87654321e-6 * randomInteger;
     10453               for (CoinBigIndex i = start; i < start + length; i++) {
     10454                    int iRow = row[i];
     10455                    sum += elementByColumn[i] * weight[iRow];
     10456               }
     10457               sortColumn[iColumn] = iColumn;
     10458               randomColumn[iColumn] = sum;
     10459          }
     10460          {
     10461               CoinSort_2(randomColumn, randomColumn + numberColumns, sortColumn);
     10462               for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     10463                    int i = sortColumn[iColumn];
     10464                    backColumn[i] = iColumn;
     10465               }
     10466               randomColumn[numberColumns] = COIN_DBL_MAX;
     10467               double last = -COIN_DBL_MAX;
     10468               int iLast = -1;
     10469               for (int iLook = 0; iLook < numberColumns + 1; iLook++) {
     10470                    if (randomColumn[iLook] > last) {
     10471                         if (iLast >= 0) {
     10472                              int n = iLook - iLast;
     10473                              if (n > 1) {
     10474                                   //printf("%d columns possible?\n",n);
     10475                              }
     10476                              for (int i = iLast; i < iLook; i++) {
     10477                                   possibleColumn[sortColumn[i]] = n;
     10478                              }
     10479                         }
     10480                         iLast = iLook;
     10481                         last = randomColumn[iLook];
     10482                    }
     10483               }
     10484               for (iRow = 0; iRow < numberRows; iRow++) {
     10485                    CoinBigIndex start = rowStart[iRow];
     10486                    double sum = 0.0;
     10487                    int length = rowLength[iRow];
     10488                    for (CoinBigIndex i = start; i < start + length; i++) {
     10489                         int jColumn = column[i];
     10490                         double value = element[i];
     10491                         jColumn = backColumn[jColumn];
     10492                         sum += value * randomColumn[jColumn];
     10493                         //if (iColumn==23089||iRow==23729)
     10494                         //printf("row %d cola %d colb %d value %g rand %g sum %g\n",
     10495                         //   iRow,jColumn,column[i],value,randomColumn[jColumn],sum);
     10496                    }
     10497                    sortRow[iRow] = iRow;
     10498                    randomRow[iRow] = weight[iRow];
     10499                    randomRow[iRow] = sum;
     10500               }
     10501               CoinSort_2(randomRow, randomRow + numberRows, sortRow);
     10502               for (iRow = 0; iRow < numberRows; iRow++) {
     10503                    int i = sortRow[iRow];
     10504                    backRow[i] = iRow;
     10505               }
     10506               randomRow[numberRows] = COIN_DBL_MAX;
     10507               last = -COIN_DBL_MAX;
     10508               iLast = -1;
     10509               // Do backward indices from order
     10510               for (iRow = 0; iRow < numberRows; iRow++) {
     10511                    other[order[iRow]] = iRow;
     10512               }
     10513               for (int iLook = 0; iLook < numberRows + 1; iLook++) {
     10514                    if (randomRow[iLook] > last) {
     10515                         if (iLast >= 0) {
     10516                              int n = iLook - iLast;
     10517                              if (n > 1) {
     10518                                   //printf("%d rows possible?\n",n);
     10519                                   // Within group sort as for original "order"
     10520                                   for (int i = iLast; i < iLook; i++) {
     10521                                        int jRow = sortRow[i];
     10522                                        order[i] = other[jRow];
     10523                                   }
     10524                                   CoinSort_2(order + iLast, order + iLook, sortRow + iLast);
     10525                              }
     10526                              for (int i = iLast; i < iLook; i++) {
     10527                                   possibleRow[sortRow[i]] = n;
     10528                              }
     10529                         }
     10530                         iLast = iLook;
     10531                         last = randomRow[iLook];
     10532                    }
     10533               }
     10534               // Temp out
     10535               for (int iLook = 0; iLook < numberRows - 1000000; iLook++) {
     10536                    iRow = sortRow[iLook];
     10537                    CoinBigIndex start = rowStart[iRow];
     10538                    int length = rowLength[iRow];
     10539                    int numberPossible = possibleRow[iRow];
     10540                    for (CoinBigIndex i = start; i < start + length; i++) {
     10541                         int jColumn = column[i];
     10542                         if (possibleColumn[jColumn] != numberPossible)
     10543                              numberPossible = -1;
     10544                    }
     10545                    int n = numberPossible;
     10546                    if (numberPossible > 1) {
     10547                         //printf("pppppossible %d\n",numberPossible);
     10548                         for (int jLook = iLook + 1; jLook < iLook + numberPossible; jLook++) {
     10549                              int jRow = sortRow[jLook];
     10550                              CoinBigIndex start2 = rowStart[jRow];
     10551                              assert (numberPossible == possibleRow[jRow]);
     10552                              assert(length == rowLength[jRow]);
     10553                              for (CoinBigIndex i = start2; i < start2 + length; i++) {
     10554                                   int jColumn = column[i];
     10555                                   if (possibleColumn[jColumn] != numberPossible)
     10556                                        numberPossible = -1;
     10557                              }
     10558                         }
     10559                         if (numberPossible < 2) {
     10560                              // switch off
     10561                              for (int jLook = iLook; jLook < iLook + n; jLook++)
     10562                                   possibleRow[sortRow[jLook]] = -1;
     10563                         }
     10564                         // skip rest
     10565                         iLook += n - 1;
     10566                    } else {
     10567                         possibleRow[iRow] = -1;
     10568                    }
     10569               }
     10570               for (int iLook = 0; iLook < numberRows; iLook++) {
     10571                    iRow = sortRow[iLook];
     10572                    int numberPossible = possibleRow[iRow];
     10573                    // Only if any integers
     10574                    int numberIntegers = 0;
     10575                    CoinBigIndex start = rowStart[iRow];
     10576                    int length = rowLength[iRow];
     10577                    for (CoinBigIndex i = start; i < start + length; i++) {
     10578                         int jColumn = column[i];
     10579                         if (model->isInteger(jColumn))
     10580                              numberIntegers++;
     10581                    }
     10582                    if (numberPossible > 1 && !numberIntegers) {
     10583                         //printf("possible %d - but no integers\n",numberPossible);
     10584                    }
     10585                    if (numberPossible > 1 && (numberIntegers || false)) {
     10586                         //
     10587                         printf("possible %d - %d integers\n", numberPossible, numberIntegers);
     10588                         int lastLook = iLook;
     10589                         int nMapRow = -1;
     10590                         for (int jLook = iLook + 1; jLook < iLook + numberPossible; jLook++) {
     10591                              // stop if too many failures
     10592                              if (jLook > iLook + 10 && nMapRow < 0)
     10593                                   break;
     10594                              // Create identity mapping
     10595                              int i;
     10596                              for (i = 0; i < numberRows; i++)
     10597                                   mapRow[i] = i;
     10598                              for (i = 0; i < numberColumns; i++)
     10599                                   mapColumn[i] = i;
     10600                              int offset = jLook - iLook;
     10601                              int nStackC = 0;
     10602                              // build up row and column mapping
     10603                              int nStackR = 1;
     10604                              stackRow[0] = iLook;
     10605                              bool good = true;
     10606                              while (nStackR) {
     10607                                   nStackR--;
     10608                                   int look1 = stackRow[nStackR];
     10609                                   int look2 = look1 + offset;
     10610                                   assert (randomRow[look1] == randomRow[look2]);
     10611                                   int row1 = sortRow[look1];
     10612                                   int row2 = sortRow[look2];
     10613                                   assert (mapRow[row1] == row1);
     10614                                   assert (mapRow[row2] == row2);
     10615                                   mapRow[row1] = row2;
     10616                                   mapRow[row2] = row1;
     10617                                   CoinBigIndex start1 = rowStart[row1];
     10618                                   CoinBigIndex offset2 = rowStart[row2] - start1;
     10619                                   int length = rowLength[row1];
     10620                                   assert( length == rowLength[row2]);
     10621                                   for (CoinBigIndex i = start1; i < start1 + length; i++) {
     10622                                        int jColumn1 = column[i];
     10623                                        int jColumn2 = column[i+offset2];
     10624                                        if (randomColumn[backColumn[jColumn1]] !=
     10625                                                  randomColumn[backColumn[jColumn2]]) {
     10626                                             good = false;
     10627                                             break;
     10628                                        }
     10629                                        if (mapColumn[jColumn1] == jColumn1) {
     10630                                             // not touched
     10631                                             assert (mapColumn[jColumn2] == jColumn2);
     10632                                             if (jColumn1 != jColumn2) {
     10633                                                  // Put on stack
     10634                                                  mapColumn[jColumn1] = jColumn2;
     10635                                                  mapColumn[jColumn2] = jColumn1;
     10636                                                  stackColumn[nStackC++] = jColumn1;
     10637                                             }
     10638                                        } else {
     10639                                             if (mapColumn[jColumn1] != jColumn2 ||
     10640                                                       mapColumn[jColumn2] != jColumn1) {
     10641                                                  // bad
     10642                                                  good = false;
     10643                                                  printf("bad col\n");
     10644                                                  break;
     10645                                             }
     10646                                        }
     10647                                   }
     10648                                   if (!good)
     10649                                        break;
     10650                                   while (nStackC) {
     10651                                        nStackC--;
     10652                                        int iColumn = stackColumn[nStackC];
     10653                                        int iColumn2 = mapColumn[iColumn];
     10654                                        assert (iColumn != iColumn2);
     10655                                        int length = columnLength[iColumn];
     10656                                        assert (length == columnLength[iColumn2]);
     10657                                        CoinBigIndex start = columnStart[iColumn];
     10658                                        CoinBigIndex offset2 = columnStart[iColumn2] - start;
     10659                                        for (CoinBigIndex i = start; i < start + length; i++) {
     10660                                             int iRow = row[i];
     10661                                             int iRow2 = row[i+offset2];
     10662                                             if (mapRow[iRow] == iRow) {
     10663                                                  // First (but be careful)
     10664                                                  if (iRow != iRow2) {
     10665                                                       //mapRow[iRow]=iRow2;
     10666                                                       //mapRow[iRow2]=iRow;
     10667                                                       int iBack = backRow[iRow];
     10668                                                       int iBack2 = backRow[iRow2];
     10669                                                       if (randomRow[iBack] == randomRow[iBack2] &&
     10670                                                                 iBack2 - iBack == offset) {
     10671                                                            stackRow[nStackR++] = iBack;
     10672                                                       } else {
     10673                                                            //printf("randomRow diff - weights %g %g\n",
     10674                                                            //     weight[iRow],weight[iRow2]);
     10675                                                            // bad
     10676                                                            good = false;
     10677                                                            break;
     10678                                                       }
     10679                                                  }
     10680                                             } else {
     10681                                                  if (mapRow[iRow] != iRow2 ||
     10682                                                            mapRow[iRow2] != iRow) {
     10683                                                       // bad
     10684                                                       good = false;
     10685                                                       printf("bad row\n");
     10686                                                       break;
     10687                                                  }
     10688                                             }
     10689                                        }
     10690                                        if (!good)
     10691                                             break;
     10692                                   }
     10693                              }
     10694                              // then check OK
     10695                              if (good) {
     10696                                   for (iRow = 0; iRow < numberRows; iRow++) {
     10697                                        CoinBigIndex start = rowStart[iRow];
     10698                                        int length = rowLength[iRow];
     10699                                        if (mapRow[iRow] == iRow) {
     10700                                             for (CoinBigIndex i = start; i < start + length; i++) {
     10701                                                  int jColumn = column[i];
     10702                                                  backColumn2[jColumn] = i - start;
     10703                                             }
     10704                                             for (CoinBigIndex i = start; i < start + length; i++) {
     10705                                                  int jColumn = column[i];
     10706                                                  if (mapColumn[jColumn] != jColumn) {
     10707                                                       int jColumn2 = mapColumn[jColumn];
     10708                                                       CoinBigIndex i2 = backColumn2[jColumn2];
     10709                                                       if (i2 < 0) {
     10710                                                            good = false;
     10711                                                       } else if (element[i] != element[i2+start]) {
     10712                                                            good = false;
     10713                                                       }
     10714                                                  }
     10715                                             }
     10716                                             for (CoinBigIndex i = start; i < start + length; i++) {
     10717                                                  int jColumn = column[i];
     10718                                                  backColumn2[jColumn] = -1;
     10719                                             }
     10720                                        } else {
     10721                                             int row2 = mapRow[iRow];
     10722                                             assert (iRow = mapRow[row2]);
     10723                                             if (rowLower[iRow] != rowLower[row2] ||
     10724                                                       rowLower[row2] != rowLower[iRow])
     10725                                                  good = false;
     10726                                             CoinBigIndex offset2 = rowStart[row2] - start;
     10727                                             for (CoinBigIndex i = start; i < start + length; i++) {
     10728                                                  int jColumn = column[i];
     10729                                                  double value = element[i];
     10730                                                  int jColumn2 = column[i+offset2];
     10731                                                  double value2 = element[i+offset2];
     10732                                                  if (value != value2 || mapColumn[jColumn] != jColumn2 ||
     10733                                                            mapColumn[jColumn2] != jColumn)
     10734                                                       good = false;
     10735                                             }
     10736                                        }
     10737                                   }
     10738                                   if (good) {
     10739                                        // check rim
     10740                                        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     10741                                             if (mapColumn[iColumn] != iColumn) {
     10742                                                  int iColumn2 = mapColumn[iColumn];
     10743                                                  if (objective[iColumn] != objective[iColumn2])
     10744                                                       good = false;
     10745                                                  if (columnLower[iColumn] != columnLower[iColumn2])
     10746                                                       good = false;
     10747                                                  if (columnUpper[iColumn] != columnUpper[iColumn2])
     10748                                                       good = false;
     10749                                                  if (model->isInteger(iColumn) != model->isInteger(iColumn2))
     10750                                                       good = false;
     10751                                             }
     10752                                        }
     10753                                   }
     10754                                   if (good) {
     10755                                        // temp
     10756                                        if (nMapRow < 0) {
     10757                                             //const double * solution = model->primalColumnSolution();
     10758                                             // find mapped
     10759                                             int nMapColumn = 0;
     10760                                             for (int i = 0; i < numberColumns; i++) {
     10761                                                  if (mapColumn[i] > i)
     10762                                                       nMapColumn++;
     10763                                             }
     10764                                             nMapRow = 0;
     10765                                             int kRow = -1;
     10766                                             for (int i = 0; i < numberRows; i++) {
     10767                                                  if (mapRow[i] > i) {
     10768                                                       nMapRow++;
     10769                                                       kRow = i;
     10770                                                  }
     10771                                             }
     10772                                             printf("%d columns, %d rows\n", nMapColumn, nMapRow);
     10773                                             if (nMapRow == 1) {
     10774                                                  CoinBigIndex start = rowStart[kRow];
     10775                                                  int length = rowLength[kRow];
     10776                                                  printf("%g <= ", rowLower[kRow]);
     10777                                                  for (CoinBigIndex i = start; i < start + length; i++) {
     10778                                                       int jColumn = column[i];
     10779                                                       if (mapColumn[jColumn] != jColumn)
     10780                                                            printf("* ");
     10781                                                       printf("%d,%g ", jColumn, element[i]);
     10782                                                  }
     10783                                                  printf("<= %g\n", rowUpper[kRow]);
     10784                                             }
     10785                                        }
     10786                                        // temp
     10787                                        int row1 = sortRow[lastLook];
     10788                                        int row2 = sortRow[jLook];
     10789                                        lastLook = jLook;
     10790                                        CoinBigIndex start1 = rowStart[row1];
     10791                                        CoinBigIndex offset2 = rowStart[row2] - start1;
     10792                                        int length = rowLength[row1];
     10793                                        assert( length == rowLength[row2]);
     10794                                        CoinBigIndex put = startAdd[nAddRows];
     10795                                        double multiplier = length < 11 ? 2.0 : 1.125;
     10796                                        double value = 1.0;
     10797                                        for (CoinBigIndex i = start1; i < start1 + length; i++) {
     10798                                             int jColumn1 = column[i];
     10799                                             int jColumn2 = column[i+offset2];
     10800                                             columnAdd[put] = jColumn1;
     10801                                             elementAdd[put++] = value;
     10802                                             columnAdd[put] = jColumn2;
     10803                                             elementAdd[put++] = -value;
     10804                                             value *= multiplier;
     10805                                        }
     10806                                        nAddRows++;
     10807                                        startAdd[nAddRows] = put;
     10808                                   } else {
     10809                                        printf("ouch - did not check out as good\n");
     10810                                   }
     10811                              }
     10812                         }
     10813                         // skip rest
     10814                         iLook += numberPossible - 1;
     10815                    }
     10816               }
     10817          }
     10818          if (nAddRows) {
     10819               double * lower = new double [nAddRows];
     10820               double * upper = new double[nAddRows];
     10821               int i;
     10822               //const double * solution = model->primalColumnSolution();
     10823               for (i = 0; i < nAddRows; i++) {
     10824                    lower[i] = 0.0;
     10825                    upper[i] = COIN_DBL_MAX;
     10826               }
     10827               printf("Adding %d rows with %d elements\n", nAddRows,
     10828                      startAdd[nAddRows]);
     10829               //ClpSimplex newModel(*model);
     10830               //newModel.addRows(nAddRows,lower,upper,startAdd,columnAdd,elementAdd);
     10831               //newModel.writeMps("modified.mps");
     10832               delete [] lower;
     10833               delete [] upper;
     10834          }
     10835          delete [] startAdd;
     10836          delete [] columnAdd;
     10837          delete [] elementAdd;
     10838          delete [] order;
     10839          delete [] other;
     10840          delete [] randomColumn;
     10841          delete [] weight;
     10842          delete [] randomRow;
     10843          delete [] sortRow;
     10844          delete [] backRow;
     10845          delete [] possibleRow;
     10846          delete [] sortColumn;
     10847          delete [] backColumn;
     10848          delete [] backColumn2;
     10849          delete [] possibleColumn;
     10850          delete [] mapRow;
     10851          delete [] mapColumn;
     10852          delete [] stackRow;
     10853          delete [] stackColumn;
     10854     }
     10855     delete [] number;
     10856     // Now do breakdown of ranges
     10857     breakdown("Elements", numberElements, elementByColumn);
     10858     breakdown("RowLower", numberRows, rowLower);
     10859     breakdown("RowUpper", numberRows, rowUpper);
     10860     breakdown("ColumnLower", numberColumns, columnLower);
     10861     breakdown("ColumnUpper", numberColumns, columnUpper);
     10862     breakdown("Objective", numberColumns, objective);
    969110863}
    969210864
  • trunk/Cbc/src/CbcThread.cpp

    r1899 r1943  
    13171317        otherModel->numberSolutions_ = numberSolutions_;
    13181318        otherModel->numberHeuristicSolutions_ = numberHeuristicSolutions_;
    1319         otherModel->numberNodes_ = 1; //numberNodes_;
     1319        otherModel->numberNodes_ = numberNodes_;
    13201320        otherModel->numberIterations_ = numberIterations_;
    13211321#ifdef JJF_ZERO
     
    15301530            tree_ = NULL;
    15311531        }
     1532        if ((moreSpecialOptions2_&32)!=0)
     1533          delete eventHandler_;
    15321534        eventHandler_ = NULL;
    15331535        delete solverCharacteristics_;
     
    15411543    } else if (mode == -1) {
    15421544        delete eventHandler_;
    1543         eventHandler_ = baseModel->eventHandler_;
     1545        if ((moreSpecialOptions2_&32)==0||!baseModel->eventHandler_) {
     1546          eventHandler_ = baseModel->eventHandler_;
     1547        } else {
     1548          eventHandler_ = baseModel->eventHandler_->clone();
     1549          eventHandler_->setModel(this);
     1550        }
    15441551        assert (!statistics_);
    15451552        assert(baseModel->solverCharacteristics_);
  • trunk/Cbc/src/CbcTree.cpp

    r1813 r1943  
    363363{
    364364    x->setNodeNumber(maximumNodeNumber_);
     365    lastObjective_ = x->objectiveValue();
     366    lastDepth_ = x->depth();
     367    lastUnsatisfied_ = x->numberUnsatisfied();
    365368    maximumNodeNumber_++;
    366369#   if CBC_DEBUG_HEAP > 2
     
    567570        if (value >= cutoff || !node->active()) {
    568571            if (node) {
     572                if (cutoff<-1.0e30)
     573                  node->nodeInfo()->deactivate(7);
    569574                nodeArray[--kDelete] = node;
    570575                depth[kDelete] = node->depth();
  • trunk/Cbc/src/CbcTree.hpp

    r1813 r1943  
    163163    inline int * newBounds() const { return newBound_; }
    164164
     165    /// Last objective in branch-and-cut search tree
     166    inline double lastObjective() const {
     167        return lastObjective_;
     168    }
     169    /// Last depth in branch-and-cut search tree
     170    inline int lastDepth() const {
     171        return lastDepth_;
     172    }
     173    /// Last number of objects unsatisfied
     174    inline int lastUnsatisfied() const {
     175        return lastUnsatisfied_;
     176    }
    165177    /// Adds branching information to complete state
    166178    void addBranchingInformation(const CbcModel * model, const CbcNodeInfo * nodeInfo,
     
    182194    /// Storage vector for the heap
    183195    std::vector <CbcNode *> nodes_;
    184           /// Sort predicate for heap ordering.
     196    /// Sort predicate for heap ordering.
    185197    CbcCompare comparison_;
    186198    /// Maximum "node" number so far to split ties
     
    190202    /// Maximum size of variable list
    191203    int maximumBranching_;
     204    /// Objective of last node pushed on tree
     205    double lastObjective_;
     206    /// Depth of last node pushed on tree
     207    int lastDepth_;
     208    /// Number unsatisfied of last node pushed on tree
     209    int lastUnsatisfied_;
    192210    /** Integer variables branched or bounded
    193211        top bit set if new upper bound
Note: See TracChangeset for help on using the changeset viewer.