Changeset 687


Ignore:
Timestamp:
Jul 15, 2007 5:28:33 PM (12 years ago)
Author:
forrest
Message:

add threaded to trunk

Location:
trunk/Cbc/src
Files:
4 added
28 edited

Legend:

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

    r640 r687  
    179179{
    180180  return NULL;
     181}
     182/* Pass in information on branch just done and create CbcObjectUpdateData instance.
     183   If object does not need data then backward pointer will be NULL.
     184   Assumes can get information from solver */
     185CbcObjectUpdateData
     186CbcObject::createUpdateInformation(const OsiSolverInterface * solver,
     187                                                        const CbcNode * node,
     188                                                        const CbcBranchingObject * branchingObject)
     189{
     190  return CbcObjectUpdateData();
    181191}
    182192 
     
    318328  return *this;
    319329}
    320 
    321  
     330// Default constructor
     331CbcObjectUpdateData::CbcObjectUpdateData()
     332  : object_(NULL),
     333    way_(0),
     334    objectNumber_(-1),
     335    change_(0.0),
     336    status_(0),
     337    intDecrease_(0),
     338    branchingValue_(0.0)
     339{
     340}
     341
     342// Useful constructor
     343CbcObjectUpdateData::CbcObjectUpdateData (CbcObject * object,
     344                                          int way,
     345                                          double change,
     346                                          int status,
     347                                          int intDecrease,
     348                                          double branchingValue)
     349  : object_(object),
     350    way_(way),
     351    objectNumber_(-1),
     352    change_(change),
     353    status_(status),
     354    intDecrease_(intDecrease),
     355    branchingValue_(branchingValue)
     356{
     357}
     358
     359// Destructor
     360CbcObjectUpdateData::~CbcObjectUpdateData ()
     361{
     362}
     363
     364// Copy constructor
     365CbcObjectUpdateData::CbcObjectUpdateData ( const CbcObjectUpdateData & rhs)
     366  : object_(rhs.object_),
     367    way_(rhs.way_),
     368    objectNumber_(rhs.objectNumber_),
     369    change_(rhs.change_),
     370    status_(rhs.status_),
     371    intDecrease_(rhs.intDecrease_),
     372    branchingValue_(rhs.branchingValue_)
     373{
     374}
     375
     376// Assignment operator
     377CbcObjectUpdateData &
     378CbcObjectUpdateData::operator=( const CbcObjectUpdateData& rhs)
     379{
     380  if (this!=&rhs) {
     381    object_ = rhs.object_;
     382    way_ = rhs.way_;
     383    objectNumber_ = rhs.objectNumber_;
     384    change_ = rhs.change_;
     385    status_ = rhs.status_;
     386    intDecrease_ = rhs.intDecrease_;
     387    branchingValue_ = rhs.branchingValue_;
     388  }
     389  return *this;
     390}
     391
     392 
  • trunk/Cbc/src/CbcBranchBase.hpp

    r640 r687  
    1515class CbcBranchingObject;
    1616class OsiChooseVariable;
     17class CbcObjectUpdateData;
    1718
    1819//#############################################################################
     
    190191                            double tolerance) const;
    191192
     193  /** Pass in information on branch just done and create CbcObjectUpdateData instance.
     194      If object does not need data then backward pointer will be NULL.
     195      Assumes can get information from solver */
     196  virtual CbcObjectUpdateData createUpdateInformation(const OsiSolverInterface * solver,
     197                                                        const CbcNode * node,
     198                                                        const CbcBranchingObject * branchingObject);
     199
     200  /// Update object by CbcObjectUpdateData
     201  virtual void updateInformation(const CbcObjectUpdateData & data) {};
     202
    192203  /// Identifier (normally column number in matrix)
    193204  inline int id() const
     
    323334  {way_=way;};
    324335
     336   /// update model
     337  inline void setModel(CbcModel * model)
     338  { model_ = model;};
    325339  /// Return model
    326340  inline CbcModel * model() const
     
    489503protected:
    490504};
     505/*  This stores data so an object can be updated
     506 */
     507class CbcObjectUpdateData {
     508
     509public:
     510
     511  /// Default Constructor
     512  CbcObjectUpdateData ();
     513
     514  /// Useful constructor
     515  CbcObjectUpdateData (CbcObject * object,
     516                       int way,
     517                       double change,
     518                       int status,
     519                       int intDecrease_,
     520                       double branchingValue);
     521 
     522  /// Copy constructor
     523  CbcObjectUpdateData ( const CbcObjectUpdateData &);
     524   
     525  /// Assignment operator
     526  CbcObjectUpdateData & operator=( const CbcObjectUpdateData& rhs);
     527
     528  /// Destructor
     529  virtual ~CbcObjectUpdateData ();
     530
     531 
     532public:
     533  /// data
     534
     535  /// Object
     536  CbcObject * object_;
     537  /// Branch as defined by instance of CbcObject
     538  int way_;
     539  /// Object number
     540  int objectNumber_;
     541  /// Change in objective
     542  double change_;
     543  /// Status 0 Optimal, 1 infeasible, 2 unknown
     544  int status_;
     545  /// Decrease in number unsatisfied
     546  int intDecrease_;
     547  /// Branching value
     548  double branchingValue_;
     549
     550};
    491551
    492552#endif
  • trunk/Cbc/src/CbcBranchDynamic.cpp

    r640 r687  
    267267{
    268268}
     269// Copy some information i.e. just variable stuff
     270void
     271CbcSimpleIntegerDynamicPseudoCost::copySome(CbcSimpleIntegerDynamicPseudoCost * otherObject)
     272{
     273  downDynamicPseudoCost_=otherObject->downDynamicPseudoCost_;
     274  upDynamicPseudoCost_=otherObject->upDynamicPseudoCost_;
     275  sumDownCost_ = otherObject->sumDownCost_;
     276  sumUpCost_ = otherObject->sumUpCost_;
     277  sumDownChange_ = otherObject->sumDownChange_;
     278  sumUpChange_ = otherObject->sumUpChange_;
     279  sumDownCostSquared_ = otherObject->sumDownCostSquared_;
     280  sumUpCostSquared_ = otherObject->sumUpCostSquared_;
     281  sumDownDecrease_ = otherObject->sumDownDecrease_;
     282  sumUpDecrease_ = otherObject->sumUpDecrease_;
     283  lastDownCost_ = otherObject->lastDownCost_;
     284  lastUpCost_ = otherObject->lastUpCost_;
     285  lastDownDecrease_ = otherObject->lastDownDecrease_;
     286  lastUpDecrease_ = otherObject->lastUpDecrease_;
     287  numberTimesDown_ = otherObject->numberTimesDown_;
     288  numberTimesUp_ = otherObject->numberTimesUp_;
     289  numberTimesDownInfeasible_ = otherObject->numberTimesDownInfeasible_;
     290  numberTimesUpInfeasible_ = otherObject->numberTimesUpInfeasible_;
     291  numberTimesDownLocalFixed_ = otherObject->numberTimesDownLocalFixed_;
     292  numberTimesUpLocalFixed_ = otherObject->numberTimesUpLocalFixed_;
     293  numberTimesDownTotalFixed_ = otherObject->numberTimesDownTotalFixed_;
     294  numberTimesUpTotalFixed_ = otherObject->numberTimesUpTotalFixed_;
     295  numberTimesProbingTotal_ = otherObject->numberTimesProbingTotal_;
     296}
    269297// Creates a branching object
    270298CbcBranchingObject *
     
    646674  double downCost = CoinMax((value-below)*downDynamicPseudoCost_,0.0);
    647675  return downCost;
     676}
     677/* Pass in information on branch just done and create CbcObjectUpdateData instance.
     678   If object does not need data then backward pointer will be NULL.
     679   Assumes can get information from solver */
     680CbcObjectUpdateData
     681CbcSimpleIntegerDynamicPseudoCost::createUpdateInformation(const OsiSolverInterface * solver,
     682                                                           const CbcNode * node,
     683                                                           const CbcBranchingObject * branchingObject)
     684{
     685  double originalValue=node->objectiveValue();
     686  int originalUnsatisfied = node->numberUnsatisfied();
     687  double objectiveValue = solver->getObjValue()*solver->getObjSense();
     688  int unsatisfied=0;
     689  int i;
     690  //might be base model - doesn't matter
     691  int numberIntegers = model_->numberIntegers();;
     692  const double * solution = solver->getColSolution();
     693  //const double * lower = solver->getColLower();
     694  //const double * upper = solver->getColUpper();
     695  double change = CoinMax(0.0,objectiveValue-originalValue);
     696  int iStatus;
     697  if (solver->isProvenOptimal())
     698    iStatus=0; // optimal
     699  else if (solver->isIterationLimitReached()
     700           &&!solver->isDualObjectiveLimitReached())
     701    iStatus=2; // unknown
     702  else
     703    iStatus=1; // infeasible
     704
     705  bool feasible = iStatus!=1;
     706  if (feasible) {
     707    double integerTolerance =
     708      model_->getDblParam(CbcModel::CbcIntegerTolerance);
     709    const int * integerVariable = model_->integerVariable();
     710    for (i=0;i<numberIntegers;i++) {
     711      int j=integerVariable[i];
     712      double value = solution[j];
     713      double nearest = floor(value+0.5);
     714      if (fabs(value-nearest)>integerTolerance)
     715        unsatisfied++;
     716    }
     717  }
     718  int way = branchingObject->way();
     719  way = - way; // because after branch so moved on
     720  double value = branchingObject->value();
     721  return CbcObjectUpdateData (this, way,
     722                                  change, iStatus,
     723                                  originalUnsatisfied-unsatisfied,value);
     724}
     725// Update object by CbcObjectUpdateData
     726void
     727CbcSimpleIntegerDynamicPseudoCost::updateInformation(const CbcObjectUpdateData & data)
     728{
     729  bool feasible = data.status_!=1;
     730  int way = data.way_;
     731  double value = data.branchingValue_;
     732  double change = data.change_;
     733#define TYPE2 1
     734#define TYPERATIO 0.9
     735  if (way<0) {
     736    // down
     737    if (feasible) {
     738      //printf("(down change %g value down %g ",change,value-floor(value));
     739      incrementNumberTimesDown();
     740      addToSumDownChange(1.0e-30+value-floor(value));
     741      addToSumDownDecrease(data.intDecrease_);
     742#if TYPE2==0
     743      addToSumDownCost(change/(1.0e-30+(value-floor(value))));
     744      setDownDynamicPseudoCost(sumDownCost()/(double) numberTimesDown());
     745#elif TYPE2==1
     746      addToSumDownCost(change);
     747      setDownDynamicPseudoCost(sumDownCost()/sumDownChange());
     748#elif TYPE2==2
     749      addToSumDownCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(value-floor(value))));
     750      setDownDynamicPseudoCost(sumDownCost()*(TYPERATIO/sumDownChange()+(1.0-TYPERATIO)/(double) numberTimesDown()));
     751#endif
     752    } else {
     753      //printf("(down infeasible value down %g ",change,value-floor(value));
     754      incrementNumberTimesDownInfeasible();
     755#if INFEAS==2
     756      double distanceToCutoff=0.0;
     757      double objectiveValue = model->getCurrentMinimizationObjValue();
     758      distanceToCutoff =  model->getCutoff()  - originalValue;
     759      if (distanceToCutoff<1.0e20)
     760        change = distanceToCutoff*2.0;
     761      else
     762        change = downDynamicPseudoCost()*(value-floor(value))*10.0;
     763      change = CoinMax(1.0e-12*(1.0+fabs(originalValue)),change);
     764      incrementNumberTimesDown();
     765      addToSumDownChange(1.0e-30+value-floor(value));
     766      addToSumDownDecrease(data.intDecrease_);
     767#if TYPE2==0
     768      addToSumDownCost(change/(1.0e-30+(value-floor(value))));
     769      setDownDynamicPseudoCost(sumDownCost()/(double) numberTimesDown());
     770#elif TYPE2==1
     771      addToSumDownCost(change);
     772      setDownDynamicPseudoCost(sumDownCost()/sumDownChange());
     773#elif TYPE2==2
     774      addToSumDownCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(value-floor(value))));
     775      setDownDynamicPseudoCost(sumDownCost()*(TYPERATIO/sumDownChange()+(1.0-TYPERATIO)/(double) numberTimesDown()));
     776#endif
     777#endif
     778    }
     779  } else {
     780    // up
     781    if (feasible) {
     782      //printf("(up change %g value down %g ",change,ceil(value)-value);
     783      incrementNumberTimesUp();
     784      addToSumUpChange(1.0e-30+ceil(value)-value);
     785      addToSumUpDecrease(data.intDecrease_);
     786#if TYPE2==0
     787      addToSumUpCost(change/(1.0e-30+(ceil(value)-value)));
     788      setUpDynamicPseudoCost(sumUpCost()/(double) numberTimesUp());
     789#elif TYPE2==1
     790      addToSumUpCost(change);
     791      setUpDynamicPseudoCost(sumUpCost()/sumUpChange());
     792#elif TYPE2==2
     793      addToSumUpCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(ceil(value)-value)));
     794      setUpDynamicPseudoCost(sumUpCost()*(TYPERATIO/sumUpChange()+(1.0-TYPERATIO)/(double) numberTimesUp()));
     795#endif
     796    } else {
     797      //printf("(up infeasible value down %g ",change,ceil(value)-value);
     798      incrementNumberTimesUpInfeasible();
     799#if INFEAS==2
     800      double distanceToCutoff=0.0;
     801      double objectiveValue = model->getCurrentMinimizationObjValue();
     802      distanceToCutoff =  model->getCutoff()  - originalValue;
     803      if (distanceToCutoff<1.0e20)
     804        change = distanceToCutoff*2.0;
     805      else
     806        change = upDynamicPseudoCost()*(ceil(value)-value)*10.0;
     807      change = CoinMax(1.0e-12*(1.0+fabs(originalValue)),change);
     808      incrementNumberTimesUp();
     809      addToSumUpChange(1.0e-30+ceil(value)-value);
     810      addToSumUpDecrease(data.intDecrease_);
     811#if TYPE2==0
     812      addToSumUpCost(change/(1.0e-30+(ceil(value)-value)));
     813      setUpDynamicPseudoCost(sumUpCost()/(double) numberTimesUp());
     814#elif TYPE2==1
     815      addToSumUpCost(change);
     816      setUpDynamicPseudoCost(sumUpCost()/sumUpChange());
     817#elif TYPE2==2
     818      addToSumUpCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(ceil(value)-value)));
     819      setUpDynamicPseudoCost(sumUpCost()*(TYPERATIO/sumUpChange()+(1.0-TYPERATIO)/(double) numberTimesUp()));
     820#endif
     821#endif
     822    }
     823  }
     824  //print(1,0.5);
    648825}
    649826// Pass in probing information
     
    9521129  int way = object_->way();
    9531130  double value = object_->value();
    954 #define TYPE2 1
    955 #define TYPERATIO 0.9
     1131  //#define TYPE2 1
     1132  //#define TYPERATIO 0.9
    9561133  if (way<0) {
    9571134    // down
     
    10431220    }
    10441221  }
    1045   //object->print();
     1222  //object->print(1,0.5);
    10461223  delete object_;
    10471224  object_=NULL;
     
    10611238                            double changeDown, int numInfDown)
    10621239{
    1063   int stateOfSearch = thisOne->model()->stateOfSearch();
     1240  CbcModel * model = thisOne->model();
     1241  int stateOfSearch = model->stateOfSearch();
    10641242  int betterWay=0;
    10651243  double value=0.0;
     1244  if (!bestObject_) {
     1245    bestCriterion_=-1.0;
     1246    bestNumberUp_=INT_MAX;
     1247    bestNumberDown_=INT_MAX;
     1248  }
    10661249  if (stateOfSearch<=1&&thisOne->model()->currentNode()->depth()>10) {
    1067 #if 0
    1068     if (!bestObject_) {
    1069       bestNumberUp_=INT_MAX;
    1070       bestNumberDown_=INT_MAX;
    1071     }
     1250#define TRY_STUFF 1
     1251#ifdef TRY_STUFF
    10721252    // before solution - choose smallest number
    10731253    // could add in depth as well
     
    11081288    }
    11091289#else
    1110     if (!bestObject_) {
    1111       bestCriterion_=-1.0;
    1112     }
    11131290    // got a solution
    11141291    double minValue = CoinMin(changeDown,changeUp);
     
    11241301#endif
    11251302  } else {
    1126     if (!bestObject_) {
    1127       bestCriterion_=-1.0;
    1128     }
     1303#if TRY_STUFF > 1
     1304    // Get current number of infeasibilities, cutoff and current objective
     1305    CbcNode * node = model->currentNode();
     1306    int numberUnsatisfied = node->numberUnsatisfied();
     1307    double cutoff = model->getCutoff();
     1308    double objectiveValue = node->objectiveValue();
     1309#endif
    11291310    // got a solution
    11301311    double minValue = CoinMin(changeDown,changeUp);
    11311312    double maxValue = CoinMax(changeDown,changeUp);
    11321313    // Reduce
     1314#ifdef TRY_STUFF
     1315    //maxValue = CoinMin(maxValue,minValue*4.0);
     1316#else
    11331317    maxValue = CoinMin(maxValue,minValue*2.0);
     1318#endif
    11341319    value = WEIGHT_AFTER*minValue + (1.0-WEIGHT_AFTER)*maxValue;
    1135     if (value>bestCriterion_+1.0e-8) {
     1320    double useValue = value;
     1321    double useBest = bestCriterion_;
     1322#if TRY_STUFF>1
     1323    int thisNumber = CoinMin(numInfUp,numInfDown);
     1324    int bestNumber = CoinMin(bestNumberUp_,bestNumberDown_);
     1325    double distance = cutoff-objectiveValue;
     1326    assert (distance>=0.0);
     1327    if (useValue+0.1*distance>useBest&&useValue*1.1>useBest&&
     1328        useBest+0.1*distance>useValue&&useBest*1.1>useValue) {
     1329      // not much in it - look at unsatisfied
     1330      if (thisNumber<numberUnsatisfied||bestNumber<numberUnsatisfied) {
     1331        double perInteger = distance/ ((double) numberUnsatisfied);
     1332        useValue += thisNumber*perInteger;
     1333        useBest += bestNumber*perInteger;
     1334      }
     1335    }
     1336#endif
     1337    if (useValue>useBest+1.0e-8) {
    11361338      if (changeUp<=changeDown) {
    11371339        betterWay=1;
  • trunk/Cbc/src/CbcBranchDynamic.hpp

    r640 r687  
    6464  virtual CbcBranchingObject * createBranch(OsiSolverInterface * solver,
    6565                                            const OsiBranchingInformation * info, int way) ;
     66
     67  /** Pass in information on branch just done and create CbcObjectUpdateData instance.
     68      If object does not need data then backward pointer will be NULL.
     69      Assumes can get information from solver */
     70  virtual CbcObjectUpdateData createUpdateInformation(const OsiSolverInterface * solver,
     71                                                        const CbcNode * node,
     72                                                        const CbcBranchingObject * branchingObject);
     73  /// Update object by CbcObjectUpdateData
     74  virtual void updateInformation(const CbcObjectUpdateData & data) ;
     75  /// Copy some information i.e. just variable stuff
     76  void copySome(CbcSimpleIntegerDynamicPseudoCost * otherObject);
    6677
    6778  /** Create an OsiSolverBranch object
  • trunk/Cbc/src/CbcCompareBase.hpp

    r485 r687  
    1818*/
    1919#include "CbcNode.hpp"
     20#include "CbcConfig.h"
    2021
    2122class CbcModel;
     
    7677  {
    7778    assert (x);
     79    assert (y);
     80#ifndef CBC_THREAD
    7881    CbcNodeInfo * infoX = x->nodeInfo();
    7982    assert (infoX);
    8083    int nodeNumberX = infoX->nodeNumber();
    81     assert (y);
    8284    CbcNodeInfo * infoY = y->nodeInfo();
    8385    assert (infoY);
     
    8587    assert (nodeNumberX!=nodeNumberY);
    8688    return (nodeNumberX>nodeNumberY);
     89#else
     90    // doesn't work if threaded
     91    assert (x!=y);
     92    return (x>y);
     93#endif
    8794  };
    8895protected:
  • trunk/Cbc/src/CbcCountRowCut.cpp

    r310 r687  
    4141                                int whichGenerator)
    4242  : OsiRowCut(rhs),
    43   owner_(info),
    44   ownerCut_(whichOne),
     43    owner_(info),
     44    ownerCut_(whichOne),
    4545    numberPointingToThis_(0),
    46   whichCutGenerator_(whichGenerator)
     46    whichCutGenerator_(whichGenerator)
    4747{
    4848#ifdef CHECK_CUT_COUNTS
     
    5858  // Look at owner and delete
    5959  owner_->deleteCut(ownerCut_);
     60  ownerCut_=-1234567;
    6061}
    6162// Increment number of references
     
    6364CbcCountRowCut::increment(int change)
    6465{
     66  assert(ownerCut_!=-1234567);
    6567  numberPointingToThis_+=change;
    6668}
     
    7072CbcCountRowCut::decrement(int change)
    7173{
    72   assert(numberPointingToThis_>=change);
     74  assert(ownerCut_!=-1234567);
     75  //assert(numberPointingToThis_>=change);
     76  assert(numberPointingToThis_>=0);
     77  if(numberPointingToThis_<change) {
     78    assert(numberPointingToThis_>0);
     79    printf("negative cut count %d - %d\n",numberPointingToThis_, change);
     80    change = numberPointingToThis_;
     81  }
    7382  numberPointingToThis_-=change;
    7483  return numberPointingToThis_;
    7584}
     85
    7686// Set information
    7787void
  • trunk/Cbc/src/CbcCutGenerator.cpp

    r640 r687  
    8383  model_ = rhs.model_;
    8484  generator_=rhs.generator_->clone();
    85   generator_->refreshSolver(model_->solver());
     85  //generator_->refreshSolver(model_->solver());
    8686  whenCutGenerator_=rhs.whenCutGenerator_;
    8787  whenCutGeneratorInSub_ = rhs.whenCutGeneratorInSub_;
     
    160160*/
    161161bool
    162 CbcCutGenerator::generateCuts( OsiCuts & cs , bool fullScan, CbcNode * node)
     162CbcCutGenerator::generateCuts( OsiCuts & cs , bool fullScan, OsiSolverInterface * solver, CbcNode * node)
    163163{
    164164  int howOften = whenCutGenerator_;
     
    172172    howOften=1;
    173173  bool returnCode=false;
    174   OsiSolverInterface * solver = model_->solver();
     174  //OsiSolverInterface * solver = model_->solver();
    175175  int depth;
    176176  if (node)
     
    210210    if (!generator) {
    211211      // Pass across model information in case it could be useful
    212       //OsiSolverInterface * solver = model_->solver();
    213212      //void * saveData = solver->getApplicationData();
    214213      //solver->setApplicationData(model_);
     
    235234      int numberColumns = solver->getNumCols();
    236235      double primalTolerance = 1.0e-8;
    237 #if 0
    238       int numberChanged=0,ifCut=0;
    239       CoinPackedVector lbs;
    240       CoinPackedVector ubs;
    241       for (j=0;j<numberColumns;j++) {
    242         if (solver->isInteger(j)) {
    243           if (tightUpper[j]<upper[j]) {
    244             numberChanged++;
    245             assert (tightUpper[j]==floor(tightUpper[j]+0.5));
    246             ubs.insert(j,tightUpper[j]);
    247             if (tightUpper[j]<solution[j]-primalTolerance)
    248               ifCut=1;
    249           }
    250           if (tightLower[j]>lower[j]) {
    251             numberChanged++;
    252             assert (tightLower[j]==floor(tightLower[j]+0.5));
    253             lbs.insert(j,tightLower[j]);
    254             if (tightLower[j]>solution[j]+primalTolerance)
    255               ifCut=1;
    256           }
    257         } else {
    258           if (tightUpper[j]==tightLower[j]&&
    259               upper[j]>lower[j]) {
    260             // fix
    261             //solver->setColLower(j,tightLower[j]);
    262             //solver->setColUpper(j,tightUpper[j]);
    263             double value = tightUpper[j];
    264             numberChanged++;
    265             if (value<upper[j])
    266               ubs.insert(j,value);
    267             if (value>lower[j])
    268               lbs.insert(j,value);
    269           }
    270         }
    271       }
    272       if (numberChanged) {
    273         OsiColCut cc;
    274         cc.setUbs(ubs);
    275         cc.setLbs(lbs);
    276         if (ifCut) {
    277           cc.setEffectiveness(100.0);
    278         } else {
    279           cc.setEffectiveness(1.0e-5);
    280         }
    281         cs.insert(cc);
    282       }
    283       // need to resolve if some bounds changed
    284       returnCode = !solver->basisIsAvailable();
    285       assert (!returnCode);
    286 #else
    287236      const char * tightenBounds = generator->tightenBounds();
    288       for (j=0;j<numberColumns;j++) {
    289         if (solver->isInteger(j)) {
    290           if (tightUpper[j]<upper[j]) {
    291             double nearest = floor(tightUpper[j]+0.5);
    292             //assert (fabs(tightUpper[j]-nearest)<1.0e-5); may be infeasible
    293             solver->setColUpper(j,nearest);
    294             if (nearest<solution[j]-primalTolerance)
    295               returnCode=true;
    296           }
    297           if (tightLower[j]>lower[j]) {
    298             double nearest = floor(tightLower[j]+0.5);
    299             //assert (fabs(tightLower[j]-nearest)<1.0e-5); may be infeasible
    300             solver->setColLower(j,nearest);
    301             if (nearest>solution[j]+primalTolerance)
    302               returnCode=true;
    303           }
    304         } else {
    305           if (upper[j]>lower[j]) {
    306             if (tightUpper[j]==tightLower[j]) {
    307               // fix
    308               solver->setColLower(j,tightLower[j]);
    309               solver->setColUpper(j,tightUpper[j]);
    310               if (tightLower[j]>solution[j]+primalTolerance||
    311                   tightUpper[j]<solution[j]-primalTolerance)
    312                 returnCode=true;
    313             } else if (tightenBounds&&tightenBounds[j]) {
    314               solver->setColLower(j,CoinMax(tightLower[j],lower[j]));
    315               solver->setColUpper(j,CoinMin(tightUpper[j],upper[j]));
    316               if (tightLower[j]>solution[j]+primalTolerance||
    317                   tightUpper[j]<solution[j]-primalTolerance)
     237      if ((model_->getThreadMode()&2)==0) {
     238        for (j=0;j<numberColumns;j++) {
     239          if (solver->isInteger(j)) {
     240            if (tightUpper[j]<upper[j]) {
     241              double nearest = floor(tightUpper[j]+0.5);
     242              //assert (fabs(tightUpper[j]-nearest)<1.0e-5); may be infeasible
     243              solver->setColUpper(j,nearest);
     244              if (nearest<solution[j]-primalTolerance)
    318245                returnCode=true;
    319246            }
     247            if (tightLower[j]>lower[j]) {
     248              double nearest = floor(tightLower[j]+0.5);
     249              //assert (fabs(tightLower[j]-nearest)<1.0e-5); may be infeasible
     250              solver->setColLower(j,nearest);
     251              if (nearest>solution[j]+primalTolerance)
     252                returnCode=true;
     253            }
     254          } else {
     255            if (upper[j]>lower[j]) {
     256              if (tightUpper[j]==tightLower[j]) {
     257                // fix
     258                solver->setColLower(j,tightLower[j]);
     259                solver->setColUpper(j,tightUpper[j]);
     260                if (tightLower[j]>solution[j]+primalTolerance||
     261                    tightUpper[j]<solution[j]-primalTolerance)
     262                  returnCode=true;
     263              } else if (tightenBounds&&tightenBounds[j]) {
     264                solver->setColLower(j,CoinMax(tightLower[j],lower[j]));
     265                solver->setColUpper(j,CoinMin(tightUpper[j],upper[j]));
     266                if (tightLower[j]>solution[j]+primalTolerance||
     267                    tightUpper[j]<solution[j]-primalTolerance)
     268                  returnCode=true;
     269              }
     270            }
    320271          }
    321         }
     272        }
     273      } else {
     274        CoinPackedVector lbs;
     275        CoinPackedVector ubs;
     276        int numberChanged=0;
     277        bool ifCut=false;
     278        for (j=0;j<numberColumns;j++) {
     279          if (solver->isInteger(j)) {
     280            if (tightUpper[j]<upper[j]) {
     281              double nearest = floor(tightUpper[j]+0.5);
     282              //assert (fabs(tightUpper[j]-nearest)<1.0e-5); may be infeasible
     283              ubs.insert(j,nearest);
     284              numberChanged++;
     285              if (nearest<solution[j]-primalTolerance)
     286                ifCut=true;
     287            }
     288            if (tightLower[j]>lower[j]) {
     289              double nearest = floor(tightLower[j]+0.5);
     290              //assert (fabs(tightLower[j]-nearest)<1.0e-5); may be infeasible
     291              lbs.insert(j,nearest);
     292              numberChanged++;
     293              if (nearest>solution[j]+primalTolerance)
     294                ifCut=true;
     295            }
     296          } else {
     297            if (upper[j]>lower[j]) {
     298              if (tightUpper[j]==tightLower[j]) {
     299                // fix
     300                lbs.insert(j,tightLower[j]);
     301                ubs.insert(j,tightUpper[j]);
     302                if (tightLower[j]>solution[j]+primalTolerance||
     303                    tightUpper[j]<solution[j]-primalTolerance)
     304                  ifCut=true;
     305              } else if (tightenBounds&&tightenBounds[j]) {
     306                lbs.insert(j,CoinMax(tightLower[j],lower[j]));
     307                ubs.insert(j,CoinMin(tightUpper[j],upper[j]));
     308                if (tightLower[j]>solution[j]+primalTolerance||
     309                    tightUpper[j]<solution[j]-primalTolerance)
     310                  ifCut=true;
     311              }
     312            }
     313          }
     314        }
     315        if (numberChanged) {
     316          OsiColCut cc;
     317          cc.setUbs(ubs);
     318          cc.setLbs(lbs);
     319          if (ifCut) {
     320            cc.setEffectiveness(100.0);
     321          } else {
     322            cc.setEffectiveness(1.0e-5);
     323          }
     324          cs.insert(cc);
     325        }
    322326      }
    323327      //if (!solver->basisIsAvailable())
    324328      //returnCode=true;
    325 #endif
    326329#if 0
    327330      // Pass across info to pseudocosts
     
    352355      int k ;
    353356      int nOdd=0;
    354       const OsiSolverInterface * solver = model_->solver();
     357      //const OsiSolverInterface * solver = model_->solver();
    355358      for (k = numberRowCutsAfter-1;k>=numberRowCutsBefore;k--) {
    356359        OsiRowCut & thisCut = cs.rowCut(k) ;
     
    382385      for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
    383386        OsiRowCut thisCut = cs.rowCut(k) ;
     387        if (thisCut.lb()>thisCut.ub()||
     388            thisCut.lb()>1.0e8||
     389            thisCut.ub()<-1.0e8)
     390          printf("cut from %s has bounds %g and %g!\n",
     391                 generatorName_,thisCut.lb(),thisCut.ub());
    384392        if (thisCut.lb()<=thisCut.ub()) {
    385393          /* check size of elements.
  • trunk/Cbc/src/CbcCutGenerator.hpp

    r640 r687  
    6464    If node then can find out depth
    6565  */
    66   bool generateCuts( OsiCuts &cs, bool fullScan,CbcNode * node);
     66  bool generateCuts( OsiCuts &cs, bool fullScan,OsiSolverInterface * solver,CbcNode * node);
    6767  //@}
    6868
     
    177177  inline double timeInCutGenerator() const
    178178  { return timeInCutGenerator_;};
     179  inline void incrementTimeInCutGenerator(double value)
     180  { timeInCutGenerator_ += value;};
    179181  /// Get the \c CglCutGenerator corresponding to this \c CbcCutGenerator.
    180182  inline CglCutGenerator * generator() const
     
    237239  inline void setNumberActiveCutsAtRoot(int value)
    238240  { numberActiveCutsAtRoot_ = value;};
     241  /// Set model
     242  inline void setModel(CbcModel * model)
     243  { model_ = model;};
    239244  //@}
    240245 
  • trunk/Cbc/src/CbcFathom.cpp

    r310 r687  
    99#include <cfloat>
    1010
    11 #include "OsiSolverInterface.hpp"
     11#ifdef COIN_HAS_CLP
     12#include "OsiClpSolverInterface.hpp"
     13#endif
    1214#include "CbcModel.hpp"
    1315#include "CbcMessage.hpp"
     
    4547  model_ = model;
    4648}
     49#ifdef COIN_HAS_CLP
     50
     51//#############################################################################
     52// Constructors, destructors clone and assignment
     53//#############################################################################
     54
     55//-------------------------------------------------------------------
     56// Default Constructor
     57//-------------------------------------------------------------------
     58CbcOsiSolver::CbcOsiSolver ()
     59  : OsiClpSolverInterface()
     60{
     61  cbcModel_ = NULL;
     62}
     63//-------------------------------------------------------------------
     64// Clone
     65//-------------------------------------------------------------------
     66OsiSolverInterface *
     67CbcOsiSolver::clone(bool copyData) const
     68{
     69  assert (copyData);
     70  return new CbcOsiSolver(*this);
     71}
     72
     73
     74//-------------------------------------------------------------------
     75// Copy constructor
     76//-------------------------------------------------------------------
     77CbcOsiSolver::CbcOsiSolver (
     78                  const CbcOsiSolver & rhs)
     79  : OsiClpSolverInterface(rhs)
     80{
     81  cbcModel_ = rhs.cbcModel_;
     82}
     83
     84//-------------------------------------------------------------------
     85// Destructor
     86//-------------------------------------------------------------------
     87CbcOsiSolver::~CbcOsiSolver ()
     88{
     89}
     90
     91//-------------------------------------------------------------------
     92// Assignment operator
     93//-------------------------------------------------------------------
     94CbcOsiSolver &
     95CbcOsiSolver::operator=(const CbcOsiSolver& rhs)
     96{
     97  if (this != &rhs) {
     98    OsiClpSolverInterface::operator=(rhs);
     99    cbcModel_ = rhs.cbcModel_;
     100  }
     101  return *this;
     102}
     103#endif
  • trunk/Cbc/src/CbcFathom.hpp

    r36 r687  
    33#ifndef CbcFathom_H
    44#define CbcFathom_H
     5#include "CbcConfig.h"
    56
    67class CbcModel;
     
    6364 
    6465};
     66#ifdef COIN_HAS_CLP
     67#include "OsiClpSolverInterface.hpp"
    6568
     69//#############################################################################
     70
     71/**
     72   
     73This is for codes where solver needs to know about CbcModel
     74*/
     75
     76class CbcOsiSolver : public OsiClpSolverInterface {
     77 
     78public:
     79 
     80  /**@name Constructors and destructors */
     81  //@{
     82  /// Default Constructor
     83  CbcOsiSolver ();
     84 
     85  /// Clone
     86  virtual OsiSolverInterface * clone(bool copyData=true) const;
     87 
     88  /// Copy constructor
     89  CbcOsiSolver (const CbcOsiSolver &);
     90 
     91  /// Assignment operator
     92  CbcOsiSolver & operator=(const CbcOsiSolver& rhs);
     93 
     94  /// Destructor
     95  virtual ~CbcOsiSolver ();
     96 
     97  //@}
     98 
     99 
     100  /**@name Sets and Gets */
     101  //@{
     102  /// Set Cbc Model
     103  inline void setCbcModel(CbcModel * model)
     104  { cbcModel_=model;};
     105  /// Return Cbc Model
     106  inline CbcModel * cbcModel() const
     107  { return cbcModel_;};
     108  //@}
     109 
     110  //---------------------------------------------------------------------------
     111 
     112protected:
     113 
     114 
     115  /**@name Private member data */
     116  //@{
     117  /// Pointer back to CbcModel
     118  CbcModel * cbcModel_;
     119  //@}
     120};
    66121#endif
     122#endif
  • trunk/Cbc/src/CbcHeuristic.cpp

    r640 r687  
    178178      // probably faster to use a basis to get integer solutions
    179179      model.setSpecialOptions(2);
     180#ifdef CBC_THREAD
     181      if (model_->getNumberThreads()>0&&(model_->getThreadMode()&1)!=0) {
     182        // See if at root node
     183        bool atRoot = model_->getNodeCount()==0;
     184        int passNumber = model_->getCurrentPassNumber();
     185        if (atRoot&&passNumber==1)
     186          model.setNumberThreads(model_->getNumberThreads());
     187      }
     188#endif
    180189      model.branchAndBound();
    181190      if (logLevel>1)
     
    302311      (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3)))
    303312    return 0; // switched off
    304 
    305313  OsiSolverInterface * solver = model_->solver();
    306314  const double * lower = solver->getColLower();
     
    315323
    316324  int numberRows = matrix_.getNumRows();
    317 
     325  assert (numberRows<=solver->getNumRows());
    318326  int numberIntegers = model_->numberIntegers();
    319327  const int * integerVariable = model_->integerVariable();
     
    322330  double newSolutionValue = direction*solver->getObjValue();
    323331  int returnCode = 0;
    324 
    325332  // Column copy
    326333  const double * element = matrix_.getElements();
     
    893900      }
    894901    }
     902    // Just in case of some stupidity
     903    double objOffset=0.0;
     904    solver->getDblParam(OsiObjOffset,objOffset);
     905    newSolutionValue = -objOffset;
     906    for ( i=0 ; i<numberColumns ; i++ )
     907      newSolutionValue += objective[i]*newSolution[i];
     908    newSolutionValue *= direction;
     909    //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
    895910    if (newSolutionValue<solutionValue) {
    896911      // paranoid check
  • trunk/Cbc/src/CbcHeuristic.hpp

    r640 r687  
    7979  inline int numberNodes() const
    8080  { return numberNodes_;};
     81  /// Just set model - do not do anything else
     82  inline void setModelOnly(CbcModel * model)
     83  { model_ = model;};
     84 
    8185
    8286  /// Sets fraction of new(rows+columns)/old(rows+columns) before doing small branch and bound (default 1.0)
  • trunk/Cbc/src/CbcHeuristicFPump.cpp

    r640 r687  
    175175                         double * betterSolution)
    176176{
     177#define LEN_PRINT 132
     178  char pumpPrint[LEN_PRINT];
     179  pumpPrint[0]='\0';
    177180  if (!when()||(when()==1&&model_->phase()!=1))
    178181    return 0; // switched off
     
    364367      }
    365368      if (maximumTime_>0.0&&CoinCpuTime()>=startTime_+maximumTime_) break;
    366       numberPasses++;
    367369      memcpy(newSolution,solution,numberColumns*sizeof(double));
    368370      int flip;
    369       returnCode = rounds(solver,newSolution,saveObjective,numberIntegers,integerVariable,
    370                           roundExpensive_,defaultRounding_,&flip);
     371      returnCode = rounds(solver, newSolution,saveObjective,numberIntegers,integerVariable,
     372                          pumpPrint,numberPasses,roundExpensive_,defaultRounding_,&flip);
     373      numberPasses++;
    371374      if (returnCode) {
    372375        // SOLUTION IS INTEGER
     
    382385          newSolutionValue += saveObjective[i]*newSolution[i];
    383386        newSolutionValue *= direction;
    384         if (model_->logLevel())
    385           printf(" - solution found of %g",newSolutionValue);
     387        sprintf(pumpPrint+strlen(pumpPrint)," - solution found of %g",newSolutionValue);
    386388        newLineNeeded=false;
    387389        if (newSolutionValue<solutionValue) {
     
    416418            solutionFound=true;
    417419            if (general&&saveValue!=newSolutionValue)
    418               printf(" - cleaned solution of %g\n",solutionValue);
    419             else
    420               printf("\n");
     420              sprintf(pumpPrint+strlen(pumpPrint)," - cleaned solution of %g",solutionValue);
     421            if (pumpPrint[0]!='\0')
     422              model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
     423                << pumpPrint
     424                <<CoinMessageEol;
     425            pumpPrint[0]='\0';
    421426          } else {
    422           if (model_->logLevel())
    423             printf(" - not good enough after small branch and bound\n");
     427            sprintf(pumpPrint+strlen(pumpPrint)," - not good enough after mini branch and bound");
     428            model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
     429              << pumpPrint
     430              <<CoinMessageEol;
     431            pumpPrint[0]='\0';
    424432          }
    425433        } else {
    426           if (model_->logLevel())
    427             printf(" - not good enough\n");
     434          sprintf(pumpPrint+strlen(pumpPrint)," - not good enough");
     435          model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
     436            << pumpPrint
     437            <<CoinMessageEol;
     438          pumpPrint[0]='\0';
    428439          newLineNeeded=false;
    429440          returnCode=0;
     
    448459        if (matched || numberPasses%100 == 0) {
    449460          // perturbation
    450           if (model_->logLevel())
    451             printf("Perturbation applied");
     461          sprintf(pumpPrint+strlen(pumpPrint)," perturbation applied");
    452462          newLineNeeded=true;
    453463          for (i=0;i<numberIntegers;i++) {
     
    531541          memcpy(newSolution,solution,numberColumns*sizeof(double));
    532542          int flip;
    533           returnCode = rounds(solver,newSolution,saveObjective,numberIntegers,integerVariable,
    534                               roundExpensive_,defaultRounding_,&flip);
     543          returnCode = rounds(solver, newSolution,saveObjective,numberIntegers,integerVariable,
     544                              pumpPrint,numberPasses,roundExpensive_,defaultRounding_,&flip);
     545          numberPasses++;
    535546          if (returnCode) {
    536547            // solution - but may not be better
     
    540551              newSolutionValue += saveObjective[i]*newSolution[i];
    541552            newSolutionValue *= direction;
    542             if (model_->logLevel())
    543               printf(" - intermediate solution found of %g",newSolutionValue);
     553            sprintf(pumpPrint+strlen(pumpPrint)," - intermediate solution found of %g",newSolutionValue);
    544554            if (newSolutionValue<solutionValue) {
    545555              memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
     
    636646          }
    637647        }
    638         if (model_->logLevel())
    639           printf("\npass %3d: obj. %10.5f --> ", numberPasses+totalNumberPasses,solver->getObjValue());
     648        if (pumpPrint[0]!='\0')
     649          model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
     650            << pumpPrint
     651            <<CoinMessageEol;
     652        pumpPrint[0]='\0';
     653        sprintf(pumpPrint+strlen(pumpPrint),"Pass %3d: obj. %10.5f --> ", numberPasses+totalNumberPasses,solver->getObjValue());
    640654        if (closestSolution&&solver->getObjValue()<closestObjectiveValue) {
    641655          int i;
     
    656670      scaleFactor *= weightFactor_;
    657671    } // END WHILE
    658     if (newLineNeeded&&model_->logLevel())
    659       printf(" - no solution found\n");
     672    if (!solutionFound) {
     673      sprintf(pumpPrint+strlen(pumpPrint),"No solution found this major pass");
     674      model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
     675        << pumpPrint
     676        <<CoinMessageEol;
     677      pumpPrint[0]='\0';
     678    }
    660679    delete solver;
    661680    for ( j=0;j<NUMBER_OLD;j++)
     
    724743      newSolver->initialSolve();
    725744      if (!newSolver->isProvenOptimal()) {
    726         newSolver->writeMps("bad.mps");
     745        //newSolver->writeMps("bad.mps");
    727746        assert (newSolver->isProvenOptimal());
    728747        break;
    729748      }
    730       printf("%d integers at bound fixed and %d continuous",
     749      sprintf(pumpPrint+strlen(pumpPrint),"Before mini branch and bound, %d integers at bound fixed and %d continuous",
    731750             nFix,nFixC);
    732       if (nFixC2+nFixI==0)
    733         printf("\n");
    734       else
    735         printf(" of which %d were internal integer and %d internal continuous\n",
    736              nFixI,nFixC2);
     751      if (nFixC2+nFixI!=0)
     752        sprintf(pumpPrint+strlen(pumpPrint)," of which %d were internal integer and %d internal continuous",
     753                nFixI,nFixC2);
     754      model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
     755        << pumpPrint
     756        <<CoinMessageEol;
     757      pumpPrint[0]='\0';
    737758      double saveValue = newSolutionValue;
    738759      returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,newSolutionValue,
    739                                        newSolutionValue,"CbcHeuristicLocalAfterFPump");
     760                                       cutoff,"CbcHeuristicLocalAfterFPump");
    740761      if ((returnCode&2)!=0) {
    741762        // could add cut
     
    743764      }
    744765      if (returnCode) {
    745         printf("old sol of %g new of %g\n",saveValue,newSolutionValue);
     766        sprintf(pumpPrint+strlen(pumpPrint),"Mini branch and bound improved solution from %g to %g",saveValue,newSolutionValue);
     767        model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
     768          << pumpPrint
     769          <<CoinMessageEol;
     770        pumpPrint[0]='\0';
    746771        memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
    747772        if (fixContinuous) {
     
    763788            double value = newSolver->getObjValue()*newSolver->getObjSense();
    764789            if (value<newSolutionValue) {
    765               printf("freeing continuous gives a solution of %g\n", value);
     790              sprintf(pumpPrint+strlen(pumpPrint),"Freeing continuous variables gives a solution of %g", value);
     791              model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
     792                << pumpPrint
     793                <<CoinMessageEol;
     794              pumpPrint[0]='\0';
    766795              newSolutionValue=value;
    767796              memcpy(betterSolution,newSolver->getColSolution(),numberColumns*sizeof(double));
    768797            }
    769798          } else {
    770             newSolver->writeMps("bad3.mps");
     799            //newSolver->writeMps("bad3.mps");
    771800          }
    772801        }
     
    786815      double gap = relativeIncrement_*fabs(solutionValue);
    787816      cutoff -= CoinMax(CoinMax(gap,absoluteIncrement_),model_->getCutoffIncrement());
    788       printf("round again with cutoff of %g\n",cutoff);
     817      sprintf(pumpPrint+strlen(pumpPrint),"Round again with cutoff of %g",cutoff);
     818      model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
     819        << pumpPrint
     820        <<CoinMessageEol;
     821      pumpPrint[0]='\0';
    789822      if ((accumulate_&3)<2&&usedColumn)
    790823        memset(usedColumn,0,numberColumns);
    791       totalNumberPasses += numberPasses;
     824      totalNumberPasses += numberPasses-1;
    792825    } else {
    793826      break;
     
    817850    newSolver->addRow(numberIntegersOrig,closestSolution,
    818851                      lastSolution,-COIN_DBL_MAX,rhs+10.0);
    819     double saveValue = newSolutionValue;
     852    //double saveValue = newSolutionValue;
    820853    //newSolver->writeMps("sub");
    821854    int returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,newSolutionValue,
     
    857890                          const double * objective,
    858891                          int numberIntegers, const int * integerVariable,
     892                          char * pumpPrint, int & iter,
    859893                          bool roundExpensive, double downValue, int *flip)
    860894{
     
    865899  int i;
    866900
    867   static int iter = 0;
    868901  int numberColumns = model_->getNumCols();
    869902  // tmp contains the current obj coefficients
     
    905938
    906939  if (nnv > nn) nnv = nn;
    907   if (iter != 0&&model_->logLevel())
    908     printf("up = %5d , down = %5d", flip_up, flip_down); fflush(stdout);
     940  if (iter != 0)
     941    sprintf(pumpPrint+strlen(pumpPrint),"up = %5d , down = %5d", flip_up, flip_down);
    909942  *flip = flip_up + flip_down;
    910943  delete [] tmp;
     
    913946  const double * columnUpper = solver->getColUpper();
    914947  if (*flip == 0 && iter != 0) {
    915     if(model_->logLevel())
    916       printf(" -- rand = %4d (%4d) ", nnv, nn);
     948    sprintf(pumpPrint+strlen(pumpPrint)," -- rand = %4d (%4d) ", nnv, nn);
    917949     for (i = 0; i < nnv; i++) {
    918950       // was solution[list[i]] = 1. - solution[list[i]]; but does that work for 7>=x>=6
     
    930962     }
    931963     *flip = nnv;
    932   } else if (model_->logLevel()) {
    933     printf(" ");
     964  } else {
     965    //sprintf(pumpPrint+strlen(pumpPrint)," ");
    934966  }
    935967  delete [] list; delete [] val;
    936   iter++;
     968  //iter++;
    937969   
    938970  const double * rowLower = solver->getRowLower();
  • trunk/Cbc/src/CbcHeuristicFPump.hpp

    r640 r687  
    165165  int rounds(OsiSolverInterface * solver,double * solution, const double * objective,
    166166             int numberIntegers, const int * integerVariable,
     167             char * pumpPrint,int & passNumber,
    167168             bool roundExpensive=false,
    168169             double downValue=0.5, int *flip=0);
  • trunk/Cbc/src/CbcHeuristicGreedy.cpp

    r640 r687  
    764764      memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
    765765      solutionValue = newSolutionValue;
    766       model_->messageHandler()->message(CBC_HEURISTIC_SOLUTION,model_->messages())
    767         << newSolutionValue
    768         << "CbcHeuristicGreedy"<<model_->getCurrentSeconds()
    769         <<CoinMessageEol;
    770766      returnCode=1;
    771767    }
  • trunk/Cbc/src/CbcHeuristicLocal.cpp

    r640 r687  
    549549        returnCode=1;
    550550        solutionValue = newSolutionValue + bestChange;
    551         if (bestChange<-1.0e-1)
    552           model_->messageHandler()->message(CBC_HEURISTIC_SOLUTION,model_->messages())
    553             << solutionValue
    554             << "CbcHeuristicLocal"<<model_->getCurrentSeconds()
    555             <<CoinMessageEol;
    556551      } else {
    557552        // bad solution - should not happen so debug if see message
  • trunk/Cbc/src/CbcLinked.cpp

    r640 r687  
    33#include "CbcConfig.h"
    44
    5 #ifndef COIN_HAS_LINK
    6 #ifdef COIN_HAS_ASL
    7 #define COIN_HAS_LINK
    8 #endif
    9 #endif
    105#include "CoinTime.hpp"
    116
     
    138#include "CoinModel.hpp"
    149#include "ClpSimplex.hpp"
     10#ifdef COIN_HAS_ASL
     11#include "ClpAmplObjective.hpp"
     12#endif
    1513// returns jColumn (-2 if linear term, -1 if unknown) and coefficient
    1614static
     
    9290  return jColumn;
    9391}
    94 #ifdef COIN_HAS_LINK
     92#include "ClpQuadraticObjective.hpp"
    9593#include <cassert>
    9694#if defined(_MSC_VER)
     
    105103#include "ClpPackedMatrix.hpp"
    106104#include "CoinTime.hpp"
    107 #include "ClpQuadraticObjective.hpp"
    108105#include "CbcModel.hpp"
    109106#include "CbcCutGenerator.hpp"
     
    180177  //sprintf(temp,"cc%d",iPass);
    181178  //writeMps(temp);
     179  //writeMps("tight");
     180  //exit(33);
    182181  //printf("wrote cc%d\n",iPass);
    183182  OsiClpSolverInterface::initialSolve();
     
    185184  if (modelPtr_->status()==0&&(secondaryStatus==2||secondaryStatus==4))
    186185    modelPtr_->cleanup(1);
    187   if (!isProvenOptimal())
    188     writeMps("yy");
     186  //if (!isProvenOptimal())
     187  //writeMps("yy");
    189188  if (isProvenOptimal()&&quadraticModel_&&modelPtr_->numberColumns()==quadraticModel_->numberColumns()) {
    190189    // see if qp can get better solution
     
    233232          int numberGenerators = cbcModel_->numberCutGenerators();
    234233          int iGenerator;
     234          cbcModel_->lockThread();
    235235          for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
    236236            CbcCutGenerator * generator = cbcModel_->cutGenerator(iGenerator);
     
    263263            }
    264264          }
     265          cbcModel_->unlockThread();
    265266        }
    266267      }
     
    461462          bestSolution_ = CoinCopyOfArray(modelPtr_->getColSolution(),modelPtr_->getNumCols());
    462463          bestObjectiveValue_ = value;
     464          if (maxIts<=10000&&cbcModel_) {
     465            OsiSolverLink * solver2 = dynamic_cast<OsiSolverLink *> (cbcModel_->solver());
     466            assert (solver2);
     467            if (solver2!=this) {
     468              // in strong branching - need to store in original solver
     469              if (value<solver2->bestObjectiveValue_-1.0e-3) {
     470                delete [] solver2->bestSolution_;
     471                solver2->bestSolution_ = CoinCopyOfArray(bestSolution_,modelPtr_->getNumCols());
     472                solver2->bestObjectiveValue_ = value;
     473              }
     474            }
     475          }
    463476          // If model has stored then add cut (if convex)
    464477          if (cbcModel_&&(specialOptions2_&4)!=0&&quadraticModel_) {
     
    470483              CglStored * gen2 = dynamic_cast<CglStored *> (gen);
    471484              if (gen2) {
     485                cbcModel_->lockThread();
    472486                // add OA cut
    473                 double offset;
     487                double offset=0.0;
    474488                int numberColumns = quadraticModel_->numberColumns();
    475489                double * gradient = new double [numberColumns+1];
     
    492506                    int yColumn = obj->yColumn();
    493507                    if (xColumn!=yColumn) {
    494                       double coefficient = 2.0*obj->coefficient();
     508                      double coefficient = /* 2.0* */obj->coefficient();
    495509                      gradient[xColumn] += coefficient*solution[yColumn];
    496510                      gradient[yColumn] += coefficient*solution[xColumn];
     
    520534                delete [] gradient;
    521535                delete [] column;
     536                cbcModel_->unlockThread();
    522537                break;
    523538              }
     
    634649          delete [] lower2;
    635650          delete [] upper2;
    636           if (isProvenOptimal())
    637             writeMps("zz");
     651          //if (isProvenOptimal())
     652          //writeMps("zz");
    638653        }
    639654      }
     
    652667            CglTemporary * gen2 = dynamic_cast<CglTemporary *> (gen);
    653668            if (gen2) {
     669              cbcModel_->lockThread();
    654670              const double * solution = getColSolution();
    655671              // add OA cut
    656               double offset;
     672              double offset=0.0;
    657673              int numberColumns = quadraticModel_->numberColumns();
    658674              double * gradient = new double [numberColumns+1];
     
    677693                  int yColumn = obj->yColumn();
    678694                  if (xColumn!=yColumn) {
    679                     double coefficient = 2.0*obj->coefficient();
     695                    double coefficient = /* 2.0* */obj->coefficient();
    680696                    gradient[xColumn] += coefficient*solution[yColumn];
    681697                    gradient[yColumn] += coefficient*solution[xColumn];
     
    710726              delete [] gradient;
    711727              delete [] column;
     728              cbcModel_->unlockThread();
    712729              break;
    713730            }
     
    735752            //const double * columnLower = modelPtr_->columnLower();
    736753            //const double * columnUpper = modelPtr_->columnUpper();
     754            cbcModel_->lockThread();
    737755            for (int iNon=0;iNon<numberNonLinearRows_;iNon++) {
    738756              int iRow = rowNonLinear_[iNon];
     
    754772                int yColumn = obj->yColumn();
    755773                if (xColumn!=yColumn) {
    756                   double coefficient = 2.0*obj->coefficient();
     774                  double coefficient = /* 2.0* */obj->coefficient();
    757775                  gradient[xColumn] += coefficient*solution[yColumn];
    758776                  gradient[yColumn] += coefficient*solution[xColumn];
     
    791809                gen2->addCut(offset-1.0e-7,COIN_DBL_MAX,n,column,gradient);
    792810            }
     811            cbcModel_->unlockThread();
    793812            delete [] gradient;
    794813            delete [] column;
     
    812831//-------------------------------------------------------------------
    813832OsiSolverLink::OsiSolverLink ()
    814   : OsiClpSolverInterface()
     833  : CbcOsiSolver()
    815834{
    816835  gutsOfDestructor(true);
     
    893912*/
    894913OsiSolverLink::OsiSolverLink ( CoinModel & coinModel)
    895   : OsiClpSolverInterface()
     914  : CbcOsiSolver()
    896915{
    897916  gutsOfDestructor(true);
     
    899918}
    900919// need bounds
    901 static void fakeBounds(OsiSolverInterface * solver,int column,double maximumValue)
     920static void fakeBounds(OsiSolverInterface * solver,int column,double maximumValue,
     921                       CoinModel * model1, CoinModel * model2)
    902922{
    903923  double lo = solver->getColLower()[column];
    904   if (lo<-maximumValue)
     924  if (lo<-maximumValue) {
    905925    solver->setColLower(column,-maximumValue);
     926    if (model1)
     927      model1->setColLower(column,-maximumValue);
     928    if (model2)
     929      model2->setColLower(column,-maximumValue);
     930  }
    906931  double up = solver->getColUpper()[column];
    907   if (up>maximumValue)
     932  if (up>maximumValue) {
    908933    solver->setColUpper(column,maximumValue);
    909 }
    910 void OsiSolverLink::load ( CoinModel & coinModel, bool tightenBounds,int logLevel)
     934    if (model1)
     935      model1->setColUpper(column,maximumValue);
     936    if (model2)
     937      model2->setColUpper(column,maximumValue);
     938  }
     939}
     940void OsiSolverLink::load ( CoinModel & coinModelOriginal, bool tightenBounds,int logLevel)
    911941{
    912942  // first check and set up arrays
    913   int numberColumns = coinModel.numberColumns();
    914   int numberRows = coinModel.numberRows();
     943  int numberColumns = coinModelOriginal.numberColumns();
     944  int numberRows = coinModelOriginal.numberRows();
    915945  // List of nonlinear entries
    916946  int * which = new int[numberColumns];
     
    920950  int numberErrors=0;
    921951  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    922     CoinModelLink triple=coinModel.firstInColumn(iColumn);
     952    CoinModelLink triple=coinModelOriginal.firstInColumn(iColumn);
    923953    bool linear=true;
    924954    int n=0;
    925955    // See if quadratic objective
    926     const char * expr = coinModel.getColumnObjectiveAsString(iColumn);
     956    const char * expr = coinModelOriginal.getColumnObjectiveAsString(iColumn);
    927957    if (strcmp(expr,"Numeric")) {
    928958      linear=false;
     
    930960    while (triple.row()>=0) {
    931961      int iRow = triple.row();
    932       const char * expr = coinModel.getElementAsString(iRow,iColumn);
     962      const char * expr = coinModelOriginal.getElementAsString(iRow,iColumn);
    933963      if (strcmp(expr,"Numeric")) {
    934964        linear=false;
    935965      }
    936       triple=coinModel.next(triple);
     966      triple=coinModelOriginal.next(triple);
    937967      n++;
    938968    }
     
    944974  if (!numberVariables_) {
    945975    delete [] which;
    946     coinModel_ = coinModel;
     976    coinModel_ = coinModelOriginal;
    947977    int nInt=0;
    948978    for (iColumn=0;iColumn<numberColumns;iColumn++) {   
     
    951981    }
    952982    printf("There are %d integers\n",nInt);
    953     loadFromCoinModel(coinModel,true);
     983    loadFromCoinModel(coinModelOriginal,true);
    954984    OsiObject ** objects = new OsiObject * [nInt];
    955985    nInt=0;
     
    968998    return;
    969999  } else {
    970     coinModel_ = coinModel;
     1000    coinModel_ = coinModelOriginal;
    9711001    // arrays for tightening bounds
    9721002    int * freeRow = new int [numberRows];
     
    10041034          ifFirst=false;
    10051035        }
    1006         coinModel.setObjective(iColumn,linearTerm);
     1036        coinModelOriginal.setObjective(iColumn,linearTerm);
    10071037      }
    10081038    }
     
    10391069            ifFirst=false;
    10401070          }
    1041           coinModel.setElement(iRow,iColumn,linearTerm);
     1071          coinModelOriginal.setElement(iRow,iColumn,linearTerm);
    10421072        }
    10431073        triple=coinModel_.next(triple);
     
    10531083    }
    10541084    printf("There are %d bilinear and %d integers\n",nBi,nInt);
    1055     loadFromCoinModel(coinModel,true);
     1085    loadFromCoinModel(coinModelOriginal,true);
     1086    CoinModel coinModel = coinModelOriginal;
    10561087    if (tightenBounds&&numberColumns<100) {
    10571088      // first fake bounds
    10581089      for (iColumn=0;iColumn<numberColumns;iColumn++) {
    10591090        if (tryColumn[iColumn]) {
    1060           fakeBounds(this,iColumn,defaultBound_);
     1091          fakeBounds(this,iColumn,defaultBound_,&coinModel,&coinModel_);
    10611092        }
    10621093      }
     
    10911122            columnLower[iColumn]=value;
    10921123            coinModel_.setColumnLower(iColumn,value);
     1124            coinModel.setColumnLower(iColumn,value);
    10931125          }
    10941126          objective[iColumn]=-1.0;
     
    11021134            columnUpper[iColumn]=value;
    11031135            coinModel_.setColumnUpper(iColumn,value);
     1136            coinModel.setColumnUpper(iColumn,value);
    11041137          }
    11051138          objective[iColumn]=0.0;
     
    11271160      // add in objective as constraint
    11281161      objectiveVariable_= numberColumns;
    1129       objectiveRow_ = modelPtr_->numberRows();
    1130       addCol(0,NULL,NULL,-COIN_DBL_MAX,COIN_DBL_MAX,1.0);
     1162      objectiveRow_ = coinModel.numberRows();
     1163      coinModel.addColumn(0,NULL,NULL,-COIN_DBL_MAX,COIN_DBL_MAX,1.0);
    11311164      int * column = new int[numberColumns+1];
    11321165      double * element = new double[numberColumns+1];
    1133       double * objective = modelPtr_->objective();
     1166      double * objective = coinModel.objectiveArray();
    11341167      int n=0;
    1135       int starts[2]={0,0};
    11361168      for (int i=0;i<numberColumns;i++) {
    11371169        if (objective[i]) {
     
    11431175      column[n]=objectiveVariable_;
    11441176      element[n++]=-1.0;
    1145       starts[1]=n;
    1146       double offset = - modelPtr_->objectiveOffset();
     1177      double offset = - coinModel.objectiveOffset();
    11471178      assert (!offset); // get sign right if happens
    1148       modelPtr_->setObjectiveOffset(0.0);
     1179      coinModel.setObjectiveOffset(0.0);
    11491180      double lowerBound = -COIN_DBL_MAX;
    1150       addRows(1,starts,column,element,&lowerBound,&offset);
     1181      coinModel.addRow(n,column,element,lowerBound,offset);
    11511182      delete [] column;
    11521183      delete [] element;
     
    11591190    double * sort = new double [nBi];
    11601191    nBi=nInt;
     1192    const OsiObject ** justBi = const_cast<const OsiObject **> (objects+nInt);
    11611193    for (iColumn=0;iColumn<numberColumns;iColumn++) {
    11621194      if (quadraticObjective)
     
    11661198      if (strcmp(expr,"Numeric")) {
    11671199        // need bounds
    1168         fakeBounds(this,iColumn,defaultBound_);
     1200        fakeBounds(this,iColumn,defaultBound_,&coinModel,&coinModel_);
    11691201        // value*x*y
    11701202        char temp[20000];
     
    11791211            if (quadraticObjective) {
    11801212              columnQuadratic[numberQuadratic]=jColumn;
    1181               elementQuadratic[numberQuadratic++]=2.0*value; // convention
     1213              if (jColumn==iColumn)
     1214                elementQuadratic[numberQuadratic++]=2.0*value; // convention
     1215              else
     1216                elementQuadratic[numberQuadratic++]=1.0*value; // convention
    11821217            }
    11831218            // need bounds
    1184             fakeBounds(this,jColumn,defaultBound_);
     1219            fakeBounds(this,jColumn,defaultBound_,&coinModel,&coinModel_);
    11851220            double meshI = coinModel_.isInteger(iColumn) ? 1.0 : 0.0;
    11861221            if (meshI)
     
    12091244              meshJ=0.0;
    12101245            }
    1211             OsiBiLinear * newObj = new OsiBiLinear(this,iColumn,jColumn,objectiveRow_,value,meshI,meshJ,
    1212                                                    nBi-nInt,(const OsiObject **) (objects+nInt));
     1246            OsiBiLinear * newObj = new OsiBiLinear(&coinModel,iColumn,jColumn,objectiveRow_,value,meshI,meshJ,
     1247                                                   nBi-nInt,justBi);
    12131248            newObj->setPriority(biLinearPriority_);
    12141249            objects[nBi++] = newObj;
     
    12401275        if (strcmp("Numeric",el)) {
    12411276          // need bounds
    1242           fakeBounds(this,iColumn,defaultBound_);
     1277          fakeBounds(this,iColumn,defaultBound_,&coinModel,&coinModel_);
    12431278          // value*x*y
    12441279          char temp[20000];
     
    12521287            if (jColumn>=0) {
    12531288              // need bounds
    1254               fakeBounds(this,jColumn,defaultBound_);
     1289              fakeBounds(this,jColumn,defaultBound_,&coinModel,&coinModel_);
    12551290              double meshI = coinModel_.isInteger(iColumn) ? 1.0 : 0.0;
    12561291              if (meshI)
     
    12791314                meshJ=0.0;
    12801315              }
    1281               OsiBiLinear * newObj = new OsiBiLinear(this,iColumn,jColumn,iRow,value,meshI,meshJ,
    1282                                                      nBi-nInt,(const OsiObject **) (objects+nInt));
     1316              OsiBiLinear * newObj = new OsiBiLinear(&coinModel,iColumn,jColumn,iRow,value,meshI,meshJ,
     1317                                                     nBi-nInt,justBi);
    12831318              newObj->setPriority(biLinearPriority_);
    12841319              objects[nBi++] = newObj;
     
    13081343             stats[0],stats[1],stats[2],nBi-nInt-nDiff);
    13091344    }
     1345    // reload with all bilinear stuff
     1346    loadFromCoinModel(coinModel,true);
     1347    //exit(77);
    13101348    nInt=0;
    13111349    for (iColumn=0;iColumn<numberColumns;iColumn++) {   
     
    14581496    }
    14591497  }
     1498  delete [] which;
     1499  if ((specialOptions2_&16)!=0)
     1500    addTighterConstraints();
     1501}
     1502// Add reformulated bilinear constraints
     1503void
     1504OsiSolverLink::addTighterConstraints()
     1505{
     1506  // This is first attempt - for now get working on trimloss
     1507  int numberW=0;
     1508  int * xW = new int[numberObjects_];
     1509  int * yW = new int[numberObjects_];
     1510  // Points to firstlambda
     1511  int * wW = new int[numberObjects_];
     1512  // Coefficient
     1513  double * alphaW = new double[numberObjects_];
     1514  // Objects
     1515  OsiBiLinear ** objW = new OsiBiLinear * [numberObjects_];
     1516  int numberColumns = getNumCols();
     1517  int firstLambda=numberColumns;
     1518  // set up list (better to rethink and do properly as column ordered)
     1519  int * list = new int[numberColumns];
     1520  memset(list,0,numberColumns*sizeof(int));
     1521  int i;
     1522  for ( i =0;i<numberObjects_;i++) {
     1523    OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
     1524    if (obj) {
     1525      //obj->setBranchingStrategy(4); // ***** temp
     1526      objW[numberW]=obj;
     1527      xW[numberW]=obj->xColumn();
     1528      yW[numberW]=obj->yColumn();
     1529      list[xW[numberW]]=1;
     1530      list[yW[numberW]]=1;
     1531      wW[numberW]=obj->firstLambda();
     1532      firstLambda = CoinMin(firstLambda,obj->firstLambda());
     1533      alphaW[numberW]=obj->coefficient();
     1534      //assert (alphaW[numberW]==1.0); // fix when occurs
     1535      numberW++;
     1536    }
     1537  }
     1538  int nList = 0;
     1539  for (i=0;i<numberColumns;i++) {
     1540    if (list[i])
     1541      list[nList++]=i;
     1542  }
     1543  // set up mark array
     1544  char * mark = new char [firstLambda*firstLambda];
     1545  memset(mark,0,firstLambda*firstLambda);
     1546  for (i=0;i<numberW;i++) {
     1547    int x = xW[i];
     1548    int y = yW[i];
     1549    mark[x*firstLambda+y]=1;
     1550    mark[y*firstLambda+x]=1;
     1551  }
     1552  int numberRows2 = originalRowCopy_->getNumRows();
     1553  int * addColumn = new int [numberColumns];
     1554  double * addElement = new double [numberColumns];
     1555  int * addW = new int [numberColumns];
     1556  assert (objectiveRow_<0); // fix when occurs
     1557  for (int iRow=0;iRow<numberRows2;iRow++) {
     1558    for (int iList=0;iList<nList;iList++) {
     1559      int kColumn = list[iList];
     1560      const double * columnLower = getColLower();
     1561      //const double * columnUpper = getColUpper();
     1562      const double * rowLower = getRowLower();
     1563      const double * rowUpper = getRowUpper();
     1564      const CoinPackedMatrix * rowCopy = getMatrixByRow();
     1565      const double * element = rowCopy->getElements();
     1566      const int * column = rowCopy->getIndices();
     1567      const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     1568      const int * rowLength = rowCopy->getVectorLengths();
     1569      CoinBigIndex j;
     1570      int numberElements = rowLength[iRow];
     1571      int n=0;
     1572      for (j=rowStart[iRow];j<rowStart[iRow]+numberElements;j++) {
     1573        int iColumn = column[j];
     1574        if (iColumn>=firstLambda) {
     1575          // no good
     1576          n=-1;
     1577          break;
     1578        }
     1579        if (mark[iColumn*firstLambda+kColumn])
     1580          n++;
     1581      }
     1582      if (n==numberElements) {
     1583        printf("can add row %d\n",iRow);
     1584        assert (columnLower[kColumn]>=0); // might be able to fix
     1585        n=0;
     1586        for (j=rowStart[iRow];j<rowStart[iRow]+numberElements;j++) {
     1587          int xColumn=kColumn;
     1588          int yColumn = column[j];
     1589          int k;
     1590          for (k=0;k<numberW;k++) {
     1591            if ((xW[k]==yColumn&&yW[k]==xColumn)||
     1592                (yW[k]==yColumn&&xW[k]==xColumn))
     1593              break;
     1594          }
     1595          assert (k<numberW);
     1596          if (xW[k]!=xColumn) {
     1597            int temp=xColumn;
     1598            xColumn=yColumn;
     1599            yColumn=temp;
     1600          }
     1601          addW[n/4]=k;
     1602          int start = wW[k];
     1603          double value = element[j];
     1604          for (int kk=0;kk<4;kk++) {
     1605            // Dummy value
     1606            addElement[n]= value;
     1607            addColumn[n++]=start+kk;
     1608          }
     1609        }
     1610        addColumn[n++] = kColumn;
     1611        double lo = rowLower[iRow];
     1612        double up = rowUpper[iRow];
     1613        if (lo>-1.0e20) {
     1614          // and tell object
     1615          for (j=0;j<n-1;j+=4) {
     1616            int iObject = addW[j/4];
     1617            objW[iObject]->addExtraRow(matrix_->getNumRows(),addElement[j]);
     1618          }
     1619          addElement[n-1]=-lo;
     1620          if (lo==up)
     1621            addRow(n,addColumn,addElement,0.0,0.0);
     1622          else
     1623            addRow(n,addColumn,addElement,0.0,COIN_DBL_MAX);
     1624          matrix_->appendRow(n,addColumn,addElement);
     1625        }
     1626        if (up<1.0e20&&up>lo) {
     1627          // and tell object
     1628          for (j=0;j<n-1;j+=4) {
     1629            int iObject = addW[j/4];
     1630            objW[iObject]->addExtraRow(matrix_->getNumRows(),addElement[j]);
     1631          }
     1632          addElement[n-1]=-up;
     1633          addRow(n,addColumn,addElement,-COIN_DBL_MAX,0.0);
     1634          matrix_->appendRow(n,addColumn,addElement);
     1635        }
     1636      }
     1637    }
     1638  }
    14601639#if 0
    1461   //int fake[10]={3,4,2,5,5,3,1,11,2,0};
    1462   int fake[10]={11,5,5,4,3,3,2,2,1,0};
    1463   for (int kk=0;kk<10;kk++) {
    1464     setColUpper(kk,fake[kk]);
    1465     setColLower(kk,fake[kk]);
     1640  // possibly do bounds
     1641  for (int iColumn=0;iColumn<firstLambda;iColumn++) {
     1642    for (int iList=0;iList<nList;iList++) {
     1643      int kColumn = list[iList];
     1644      const double * columnLower = getColLower();
     1645      const double * columnUpper = getColUpper();
     1646      if (mark[iColumn*firstLambda+kColumn]) {
     1647        printf("can add column %d\n",iColumn);
     1648        assert (columnLower[kColumn]>=0); // might be able to fix
     1649        int xColumn=kColumn;
     1650        int yColumn = iColumn;
     1651        int k;
     1652        for (k=0;k<numberW;k++) {
     1653          if ((xW[k]==yColumn&&yW[k]==xColumn)||
     1654              (yW[k]==yColumn&&xW[k]==xColumn))
     1655            break;
     1656        }
     1657        assert (k<numberW);
     1658        if (xW[k]!=xColumn) {
     1659          int temp=xColumn;
     1660          xColumn=yColumn;
     1661          yColumn=temp;
     1662        }
     1663        int start = wW[k];
     1664        int n=0;
     1665        for (int kk=0;kk<4;kk++) {
     1666          // Dummy value
     1667          addElement[n]= 1.0e-19;
     1668          addColumn[n++]=start+kk;
     1669        }
     1670        // Tell object about this
     1671        objW[k]->addExtraRow(matrix_->getNumRows(),1.0);
     1672        addColumn[n++] = kColumn;
     1673        double lo = columnLower[iColumn];
     1674        double up = columnUpper[iColumn];
     1675        if (lo>-1.0e20) {
     1676          addElement[n-1]=-lo;
     1677          if (lo==up)
     1678            addRow(n,addColumn,addElement,0.0,0.0);
     1679          else
     1680            addRow(n,addColumn,addElement,0.0,COIN_DBL_MAX);
     1681          matrix_->appendRow(n,addColumn,addElement);
     1682        }
     1683        if (up<1.0e20&&up>lo) {
     1684          addElement[n-1]=-up;
     1685          addRow(n,addColumn,addElement,-COIN_DBL_MAX,0.0);
     1686          matrix_->appendRow(n,addColumn,addElement);
     1687        }
     1688      }
     1689    }
    14661690  }
    14671691#endif
    1468   delete [] which;
     1692  delete [] xW;
     1693  delete [] yW;
     1694  delete [] wW;
     1695  delete [] alphaW;
     1696  delete [] addColumn;
     1697  delete [] addElement;
     1698  delete [] addW;
     1699  delete [] mark;
     1700  delete [] list;
     1701  delete [] objW;
    14691702}
    14701703// Set all biLinear priorities on x-x variables
     
    14911724        objNew->setYOtherSatisfied(oldSatisfied);
    14921725        objNew->setYMeshSize(meshSize);
    1493         objNew->setXYSatisfied(0.5*meshSize);
     1726        objNew->setXYSatisfied(0.25*meshSize);
    14941727        objNew->setPriority(value);
    14951728        objNew->setBranchingStrategy(8);
     
    15021735  delete [] newObject;
    15031736}
     1737/* Set options and priority on all or some biLinear variables
     1738   1 - on I-I
     1739   2 - on I-x
     1740   4 - on x-x
     1741      or combinations.
     1742      -1 means leave (for priority value and strategy value)
     1743*/
     1744void
     1745OsiSolverLink::setBranchingStrategyOnVariables(int strategyValue, int priorityValue,
     1746                                               int mode)
     1747{
     1748  int i;
     1749  for ( i =0;i<numberObjects_;i++) {
     1750    OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
     1751    if (obj) {
     1752      bool change=false;
     1753      if (obj->xMeshSize()<1.0&&obj->yMeshSize()<1.0&&(mode&4)!=0)
     1754        change=true;
     1755      else if (((obj->xMeshSize()==1.0&&obj->yMeshSize()<1.0)||
     1756                (obj->xMeshSize()<1.0&&obj->yMeshSize()==1.0))&&(mode&2)!=0)
     1757        change=true;
     1758      else if (obj->xMeshSize()==1.0&&obj->yMeshSize()==1.0&&(mode&1)!=0)
     1759        change=true;
     1760      else if (obj->xMeshSize()>1.0||obj->yMeshSize()>1.0)
     1761        abort();
     1762      if (change) {
     1763        if (strategyValue>=0)
     1764          obj->setBranchingStrategy(strategyValue);
     1765        if (priorityValue>=0)
     1766          obj->setPriority(priorityValue);
     1767      }
     1768    }
     1769  }
     1770}
     1771
    15041772// Say convex (should work it out)
    15051773void
     
    17121980  }
    17131981  model.initialSolve();
    1714   model.writeMps("bad.mps");
     1982  //model.writeMps("bad.mps");
    17151983  // redo number of columns
    17161984  numberColumns = model.numberColumns();
     
    20792347      char temp[20];
    20802348      sprintf(temp,"pass%d.mps",iPass);
    2081       model.writeMps(temp);
     2349      //model.writeMps(temp);
    20822350#ifdef CLP_DEBUG
    20832351      if (model.status()) {
     
    21382406      }
    21392407      model.primal(1);
    2140       model.writeMps("xx.mps");
     2408      //model.writeMps("xx.mps");
    21412409      iPass--;
    21422410      goodMove=-1;
     
    22702538    cbcModel->cutGenerator(6)->setTiming(true);
    22712539    // For now - switch off most heuristics (because CglPreProcess is bad with QP)
    2272 #if 0   
     2540#if 1   
    22732541    CbcHeuristicFPump heuristicFPump(*cbcModel);
    22742542    heuristicFPump.setWhen(13);
     
    22912559   
    22922560    CbcRounding rounding(*cbcModel);
     2561    rounding.setHeuristicName("rounding");
    22932562    cbcModel->addHeuristic(&rounding);
    22942563   
     
    23292598    // if convex
    23302599    if ((specialOptions2()&4)!=0) {
     2600      if (cbcModel_)
     2601        cbcModel_->lockThread();
    23312602      // add OA cut
    23322603      double offset;
     
    23502621      delete [] gradient;
    23512622      delete [] column;
     2623      if (cbcModel_)
     2624        cbcModel_->unlockThread();
    23522625    }
    23532626    delete qp;
     
    23692642  double * solution = CoinCopyOfArray(temp->primalColumnSolution(),numberColumns);
    23702643  delete temp;
     2644  if (mode==0) {
     2645    return solution;
     2646  } else if (mode==2) {
     2647    const double * lower = getColLower();
     2648    const double * upper = getColUpper();
     2649    for (int iObject =0;iObject<numberObjects_;iObject++) {
     2650      OsiSimpleInteger * obj = dynamic_cast<OsiSimpleInteger *> (object_[iObject]);
     2651      if (obj&&(obj->priority()<biLinearPriority_||biLinearPriority_<=0)) {
     2652        int iColumn = obj->columnNumber();
     2653        double value = solution[iColumn];
     2654        value = floor(value+0.5);
     2655        if (fabs(value-solution[iColumn])>0.01) {
     2656          setColLower(iColumn,CoinMax(lower[iColumn],value-CoinMax(defaultBound_,0.0)));
     2657          setColUpper(iColumn,CoinMin(upper[iColumn],value+CoinMax(defaultBound_,1.0)));
     2658        } else {
     2659          // could fix to integer
     2660          setColLower(iColumn,CoinMax(lower[iColumn],value-CoinMax(defaultBound_,0.0)));
     2661          setColUpper(iColumn,CoinMin(upper[iColumn],value+CoinMax(defaultBound_,0.0)));
     2662        }
     2663      }
     2664    }
     2665    return solution;
     2666  }
    23712667  OsiClpSolverInterface newSolver;
    23722668  if (mode==1) {
     
    23762672    tempModel = coinModel_;
    23772673    // solve modified problem
     2674    char * mark = new char[numberColumns];
     2675    memset(mark,0,numberColumns);
    23782676    for (int iObject =0;iObject<numberObjects_;iObject++) {
    23792677      OsiSimpleInteger * obj = dynamic_cast<OsiSimpleInteger *> (object_[iObject]);
     
    23832681        value = ceil(value-1.0e-7);
    23842682        tempModel.associateElement(coinModel_.columnName(iColumn),value);
    2385       }
    2386     }
     2683        mark[iColumn]=1;
     2684      }
     2685      OsiBiLinear * objB = dynamic_cast<OsiBiLinear *> (object_[iObject]);
     2686      if (objB) {
     2687        // if one or both continuous then fix one
     2688        if (objB->xMeshSize()<1.0) {
     2689          int xColumn = objB->xColumn();
     2690          double value = solution[xColumn];
     2691          tempModel.associateElement(coinModel_.columnName(xColumn),value);
     2692          mark[xColumn]=1;
     2693        } else if (objB->yMeshSize()<1.0) {
     2694          int yColumn = objB->yColumn();
     2695          double value = solution[yColumn];
     2696          tempModel.associateElement(coinModel_.columnName(yColumn),value);
     2697          mark[yColumn]=1;
     2698        }
     2699      }
     2700    }
     2701    CoinModel * reOrdered = tempModel.reorder(mark);
     2702    assert (reOrdered);
     2703    tempModel=*reOrdered;
     2704    delete reOrdered;
     2705    delete [] mark;
    23872706    newSolver.loadFromCoinModel(tempModel,true);
    23882707    for (int iObject =0;iObject<numberObjects_;iObject++) {
     
    23942713        newSolver.setColLower(iColumn,value);
    23952714        newSolver.setColUpper(iColumn,value);
     2715      }
     2716      OsiBiLinear * objB = dynamic_cast<OsiBiLinear *> (object_[iObject]);
     2717      if (objB) {
     2718        // if one or both continuous then fix one
     2719        if (objB->xMeshSize()<1.0) {
     2720          int xColumn = objB->xColumn();
     2721          double value = solution[xColumn];
     2722          newSolver.setColLower(xColumn,value);
     2723          newSolver.setColUpper(xColumn,value);
     2724        } else if (objB->yMeshSize()<1.0) {
     2725          int yColumn = objB->yColumn();
     2726          double value = solution[yColumn];
     2727          newSolver.setColLower(yColumn,value);
     2728          newSolver.setColUpper(yColumn,value);
     2729        }
    23962730      }
    23972731    }
     
    24912825    clpModel->setLogLevel(saveLogLevel);
    24922826    returnCode=-1; // infeasible//std::cout<<"Problem is infeasible - tightenPrimalBounds!"<<std::endl;
    2493     clpModel->writeMps("infeas2.mps");
     2827    //clpModel->writeMps("infeas2.mps");
    24942828  } else {
    24952829    clpModel->setLogLevel(saveLogLevel);
     
    27103044{
    27113045  assert (copyData);
    2712   return new OsiSolverLink(*this);
     3046  OsiSolverLink * newModel = new OsiSolverLink(*this);
     3047  return newModel;
    27133048}
    27143049
     
    27193054OsiSolverLink::OsiSolverLink (
    27203055                  const OsiSolverLink & rhs)
    2721   : OsiClpSolverInterface(rhs)
     3056  : CbcOsiSolver(rhs)
    27223057{
    27233058  gutsOfDestructor(true);
     
    27433078  if (this != &rhs) {
    27443079    gutsOfDestructor();
    2745     OsiClpSolverInterface::operator=(rhs);
     3080    CbcOsiSolver::operator=(rhs);
    27463081    gutsOfCopy(rhs);
    27473082  }
     
    27713106  convex_ = NULL;
    27723107  whichNonLinear_ = NULL;
    2773   cbcModel_ = NULL;
    27743108  info_ = NULL;
    27753109  fixVariables_=NULL;
     
    28013135  biLinearPriority_ = rhs.biLinearPriority_;
    28023136  numberFix_ = rhs.numberFix_;
    2803   cbcModel_ = rhs.cbcModel_;
    28043137  if (numberVariables_) {
    28053138    if (rhs.matrix_)
     
    30943427      }
    30953428    }
    3096     newSolver.writeMps("xx");
     3429    //newSolver.writeMps("xx");
    30973430    CbcModel model(newSolver);
    30983431    // Now do requested saves and modifications
     
    33183651          setRowBounds(i,-COIN_DBL_MAX,COIN_DBL_MAX);
    33193652        initialSolve();
    3320         if (!isProvenOptimal())
    3321           getModelPtr()->writeMps("bad.mps");
     3653        //if (!isProvenOptimal())
     3654        //getModelPtr()->writeMps("bad.mps");
    33223655        if (isProvenOptimal()) {
    33233656          delete [] bestSolution_;
     
    43024635    xyRow_(-1),
    43034636    convexity_(-1),
     4637    numberExtraRows_(0),
     4638    multiplier_(NULL),
     4639    extraRow_(NULL),
    43044640    chosen_(-1)
    43054641{
     
    43304666    xyRow_(xyRow),
    43314667    convexity_(-1),
     4668    numberExtraRows_(0),
     4669    multiplier_(NULL),
     4670    extraRow_(NULL),
    43324671    chosen_(-1)
    43334672{
     
    45174856  }
    45184857}
     4858// Useful constructor
     4859OsiBiLinear::OsiBiLinear (CoinModel * coinModel, int xColumn,
     4860                          int yColumn, int xyRow, double coefficient,
     4861                          double xMesh, double yMesh,
     4862                          int numberExistingObjects,const OsiObject ** objects )
     4863  : OsiObject2(),
     4864    coefficient_(coefficient),
     4865    xMeshSize_(xMesh),
     4866    yMeshSize_(yMesh),
     4867    xSatisfied_(1.0e-6),
     4868    ySatisfied_(1.0e-6),
     4869    xOtherSatisfied_(0.0),
     4870    yOtherSatisfied_(0.0),
     4871    xySatisfied_(1.0e-6),
     4872    xyBranchValue_(0.0),
     4873    xColumn_(xColumn),
     4874    yColumn_(yColumn),
     4875    firstLambda_(-1),
     4876    branchingStrategy_(0),
     4877    boundType_(0),
     4878    xRow_(-1),
     4879    yRow_(-1),
     4880    xyRow_(xyRow),
     4881    convexity_(-1),
     4882    numberExtraRows_(0),
     4883    multiplier_(NULL),
     4884    extraRow_(NULL),
     4885    chosen_(-1)
     4886{
     4887  double columnLower[4];
     4888  double columnUpper[4];
     4889  double objective[4];
     4890  double rowLower[3];
     4891  double rowUpper[3];
     4892  CoinBigIndex starts[5];
     4893  int index[16];
     4894  double element[16];
     4895  int i;
     4896  starts[0]=0;
     4897  // rows
     4898  int numberRows = coinModel->numberRows();
     4899  // convexity
     4900  rowLower[0]=1.0;
     4901  rowUpper[0]=1.0;
     4902  convexity_ = numberRows;
     4903  starts[1]=0;
     4904  // x
     4905  rowLower[1]=0.0;
     4906  rowUpper[1]=0.0;
     4907  index[0]=xColumn_;
     4908  element[0]=-1.0;
     4909  xRow_ = numberRows+1;
     4910  starts[2]=1;
     4911  int nAdd=2;
     4912  if (xColumn_!=yColumn_) {
     4913    rowLower[2]=0.0;
     4914    rowUpper[2]=0.0;
     4915    index[1]=yColumn;
     4916    element[1]=-1.0;
     4917    nAdd=3;
     4918    yRow_ = numberRows+2;
     4919    starts[3]=2;
     4920  } else {
     4921    yRow_=-1;
     4922    branchingStrategy_=1;
     4923  }
     4924  // may be objective
     4925  assert (xyRow_>=-1);
     4926  for (i=0;i<nAdd;i++) {
     4927    CoinBigIndex iStart = starts[i];
     4928    coinModel->addRow(starts[i+1]-iStart,index+iStart,element+iStart,rowLower[i],rowUpper[i]);
     4929  }
     4930  int n=0;
     4931  // order is LxLy, LxUy, UxLy and UxUy
     4932  firstLambda_ = coinModel->numberColumns();
     4933  // bit sloppy as theoretically could be infeasible but otherwise need to do more work
     4934  double xB[2];
     4935  double yB[2];
     4936  const double * lower = coinModel->columnLowerArray();
     4937  const double * upper = coinModel->columnUpperArray();
     4938  xB[0]=lower[xColumn_];
     4939  xB[1]=upper[xColumn_];
     4940  yB[0]=lower[yColumn_];
     4941  yB[1]=upper[yColumn_];
     4942  if (xMeshSize_!=floor(xMeshSize_)) {
     4943    // not integral
     4944    xSatisfied_ = CoinMax(xSatisfied_,0.51*xMeshSize_);
     4945    if (!yMeshSize_) {
     4946      xySatisfied_ = CoinMax(xySatisfied_,xSatisfied_*CoinMax(fabs(yB[0]),fabs(yB[1])));
     4947    }
     4948  }
     4949  if (yMeshSize_!=floor(yMeshSize_)) {
     4950    // not integral
     4951    ySatisfied_ = CoinMax(ySatisfied_,0.51*yMeshSize_);
     4952    if (!xMeshSize_) {
     4953      xySatisfied_ = CoinMax(xySatisfied_,ySatisfied_*CoinMax(fabs(xB[0]),fabs(xB[1])));
     4954    }
     4955  }
     4956  // adjust
     4957  double distance;
     4958  double steps;
     4959  if (xMeshSize_) {
     4960    distance = xB[1]-xB[0];
     4961    steps = floor ((distance+0.5*xMeshSize_)/xMeshSize_);
     4962    distance = xB[0]+xMeshSize_*steps;
     4963    if (fabs(xB[1]-distance)>xSatisfied_) {
     4964      printf("bad x mesh %g %g %g -> %g\n",xB[0],xMeshSize_,xB[1],distance);
     4965      //double newValue = CoinMax(fabs(xB[1]-distance),xMeshSize_);
     4966      //printf("xSatisfied increased to %g\n",newValue);
     4967      //xSatisfied_ = newValue;
     4968      //xB[1]=distance;
     4969      //coinModel->setColUpper(xColumn_,distance);
     4970    }
     4971  }
     4972  if (yMeshSize_) {
     4973    distance = yB[1]-yB[0];
     4974    steps = floor ((distance+0.5*yMeshSize_)/yMeshSize_);
     4975    distance = yB[0]+yMeshSize_*steps;
     4976    if (fabs(yB[1]-distance)>ySatisfied_) {
     4977      printf("bad y mesh %g %g %g -> %g\n",yB[0],yMeshSize_,yB[1],distance);
     4978      //double newValue = CoinMax(fabs(yB[1]-distance),yMeshSize_);
     4979      //printf("ySatisfied increased to %g\n",newValue);
     4980      //ySatisfied_ = newValue;
     4981      //yB[1]=distance;
     4982      //coinModel->setColUpper(yColumn_,distance);
     4983    }
     4984  }
     4985  for (i=0;i<4;i++) {
     4986    double x = (i<2) ? xB[0] : xB[1];
     4987    double y = ((i&1)==0) ? yB[0] : yB[1];
     4988    columnLower[i]=0.0;
     4989    columnUpper[i]=2.0;
     4990    objective[i]=0.0;
     4991    double value;
     4992    // xy
     4993    value=coefficient_*x*y;
     4994    if (xyRow_>=0) {
     4995      if (fabs(value)<1.0e-19)
     4996        value = 1.0e-19;
     4997      element[n]=value;
     4998      index[n++]=xyRow_;
     4999    } else {
     5000      objective[i]=value;
     5001    }
     5002    // convexity
     5003    value=1.0;
     5004    element[n]=value;
     5005    index[n++]=0+numberRows;
     5006    // x
     5007    value=x;
     5008    if (fabs(value)<1.0e-19)
     5009      value = 1.0e-19;
     5010    element[n]=value;
     5011    index[n++]=1+numberRows;
     5012    if (xColumn_!=yColumn_) {
     5013      // y
     5014      value=y;
     5015      if (fabs(value)<1.0e-19)
     5016      value = 1.0e-19;
     5017      element[n]=value;
     5018      index[n++]=2+numberRows;
     5019    }
     5020    starts[i+1]=n;
     5021  }
     5022  for (i=0;i<4;i++) {
     5023    CoinBigIndex iStart = starts[i];
     5024    coinModel->addColumn(starts[i+1]-iStart,index+iStart,element+iStart,columnLower[i],
     5025                         columnUpper[i],objective[i]);
     5026  }
     5027  // At least one has to have a mesh
     5028  if (!xMeshSize_&&(!yMeshSize_||yRow_<0)) {
     5029    printf("one of x and y must have a mesh size\n");
     5030    abort();
     5031  } else if (yRow_>=0) {
     5032    if (!xMeshSize_)
     5033      branchingStrategy_ = 2;
     5034    else if (!yMeshSize_)
     5035      branchingStrategy_ = 1;
     5036  }
     5037  // Now add constraints to link in x and or y to existing ones.
     5038  bool xDone=false;
     5039  bool yDone=false;
     5040  // order is LxLy, LxUy, UxLy and UxUy
     5041  for (i=numberExistingObjects-1;i>=0;i--) {
     5042    const OsiObject * obj = objects[i];
     5043    const OsiBiLinear * obj2 =
     5044      dynamic_cast <const OsiBiLinear *>(obj) ;
     5045    if (obj2) {
     5046      if (xColumn_==obj2->xColumn_&&!xDone) {
     5047        // make sure y equal
     5048        double rhs=0.0;
     5049        CoinBigIndex starts[2];
     5050        int index[4];
     5051        double element[4]= {1.0,1.0,-1.0,-1.0};
     5052        starts[0]=0;
     5053        starts[1]=4;
     5054        index[0]=firstLambda_+0;
     5055        index[1]=firstLambda_+1;
     5056        index[2]=obj2->firstLambda_+0;
     5057        index[3]=obj2->firstLambda_+1;
     5058        coinModel->addRow(4,index,element,rhs,rhs);
     5059        xDone=true;
     5060      }
     5061      if (yColumn_==obj2->yColumn_&&yRow_>=0&&!yDone) {
     5062        // make sure x equal
     5063        double rhs=0.0;
     5064        CoinBigIndex starts[2];
     5065        int index[4];
     5066        double element[4]= {1.0,1.0,-1.0,-1.0};
     5067        starts[0]=0;
     5068        starts[1]=4;
     5069        index[0]=firstLambda_+0;
     5070        index[1]=firstLambda_+2;
     5071        index[2]=obj2->firstLambda_+0;
     5072        index[3]=obj2->firstLambda_+2;
     5073        coinModel->addRow(4,index,element,rhs,rhs);
     5074        yDone=true;
     5075      }
     5076    }
     5077  }
     5078}
    45195079
    45205080// Copy constructor
     
    45395099   xyRow_(rhs.xyRow_),
    45405100   convexity_(rhs.convexity_),
     5101   numberExtraRows_(rhs.numberExtraRows_),
     5102   multiplier_(NULL),
     5103   extraRow_(NULL),
    45415104   chosen_(rhs.chosen_)
    45425105{
     5106  if (numberExtraRows_) {
     5107    multiplier_ = CoinCopyOfArray(rhs.multiplier_,numberExtraRows_);
     5108    extraRow_ = CoinCopyOfArray(rhs.extraRow_,numberExtraRows_);
     5109  }
    45435110}
    45445111
     
    45745141    xyRow_ = rhs.xyRow_;
    45755142    convexity_ = rhs.convexity_;
     5143    numberExtraRows_ = rhs.numberExtraRows_;
     5144    delete [] multiplier_;
     5145    delete [] extraRow_;
     5146    if (numberExtraRows_) {
     5147      multiplier_ = CoinCopyOfArray(rhs.multiplier_,numberExtraRows_);
     5148      extraRow_ = CoinCopyOfArray(rhs.extraRow_,numberExtraRows_);
     5149    } else {
     5150      multiplier_ = NULL;
     5151      extraRow_ = NULL;
     5152    }
    45765153    chosen_ = rhs.chosen_;
    45775154  }
     
    45825159OsiBiLinear::~OsiBiLinear ()
    45835160{
    4584 }
    4585 
     5161  delete [] multiplier_;
     5162  delete [] extraRow_;
     5163}
     5164// Adds in data for extra row with variable coefficients
     5165void
     5166OsiBiLinear::addExtraRow(int row, double multiplier)
     5167{
     5168  int * tempI = new int [numberExtraRows_+1];
     5169  double * tempD = new double [numberExtraRows_+1];
     5170  memcpy(tempI,extraRow_,numberExtraRows_*sizeof(int));
     5171  memcpy(tempD,multiplier_,numberExtraRows_*sizeof(double));
     5172  tempI[numberExtraRows_]=row;
     5173  tempD[numberExtraRows_]=multiplier;
     5174  if (numberExtraRows_)
     5175    assert (row>tempI[numberExtraRows_-1]);
     5176  numberExtraRows_++;
     5177  delete [] extraRow_;
     5178  extraRow_ = tempI;
     5179  delete [] multiplier_;
     5180  multiplier_ = tempD;
     5181}
     5182static bool testCoarse=true;
    45865183// Infeasibility - large is 0.5
    45875184double
     
    46575254  double steps;
    46585255  bool xSatisfied;
    4659   double xNew;
     5256  double xNew=xB[0];
    46605257  if (xMeshSize_) {
    46615258    if (x<0.5*(xB[0]+xB[1])) {
     
    46735270    }
    46745271    // but if first coarse grid then only if gap small
    4675     if (false&&(branchingStrategy_&8)!=0&&xSatisfied&&
     5272    if (testCoarse&&(branchingStrategy_&8)!=0&&xSatisfied&&
    46765273        xB[1]-xB[0]>=xMeshSize_) {
    46775274      // but allow if fine grid would allow
     
    46875284  }
    46885285  bool ySatisfied;
    4689   double yNew;
     5286  double yNew=yB[0];
    46905287  if (yMeshSize_) {
    46915288    if (y<0.5*(yB[0]+yB[1])) {
     
    47035300    }
    47045301    // but if first coarse grid then only if gap small
    4705     if (false&&(branchingStrategy_&8)!=0&&ySatisfied&&
     5302    if (testCoarse&&(branchingStrategy_&8)!=0&&ySatisfied&&
    47065303        yB[1]-yB[0]>=yMeshSize_) {
    47075304      // but allow if fine grid would allow
     
    47655362    xyLambda /= coefficient_;
    47665363  }
     5364  if (0) {
     5365    // only true with positive values
     5366    // see if all convexification constraints OK with true
     5367    assert (xyTrue+1.0e-5>xB[0]*y+yB[0]*x - xB[0]*yB[0]);
     5368    assert (xyTrue+1.0e-5>xB[1]*y+yB[1]*x - xB[1]*yB[1]);
     5369    assert (xyTrue-1.0e-5<xB[1]*y+yB[0]*x - xB[1]*yB[0]);
     5370    assert (xyTrue-1.0e-5<xB[0]*y+yB[1]*x - xB[0]*yB[1]);
     5371    // see if all convexification constraints OK with lambda version
     5372#if 1
     5373    assert (xyLambda+1.0e-5>xB[0]*y+yB[0]*x - xB[0]*yB[0]);
     5374    assert (xyLambda+1.0e-5>xB[1]*y+yB[1]*x - xB[1]*yB[1]);
     5375    assert (xyLambda-1.0e-5<xB[1]*y+yB[0]*x - xB[1]*yB[0]);
     5376    assert (xyLambda-1.0e-5<xB[0]*y+yB[1]*x - xB[0]*yB[1]);
     5377#endif
     5378    // see if other bound stuff true
     5379    assert (xyLambda+1.0e-5>xB[0]*y);
     5380    assert (xyLambda+1.0e-5>yB[0]*x);
     5381    assert (xyLambda-1.0e-5<xB[1]*y);
     5382    assert (xyLambda-1.0e-5<yB[1]*x);
     5383#define SIZE 2
     5384    if (yColumn_==xColumn_+SIZE) {
     5385#if SIZE==6
     5386      double bMax = 2200.0;
     5387      double bMin = bMax - 100.0;
     5388      double b[] = {330.0,360.0,380.0,430.0,490.0,530.0};
     5389#elif SIZE==2
     5390      double bMax = 1900.0;
     5391      double bMin = bMax - 200.0;
     5392      double b[] = {460.0,570.0};
     5393#else
     5394      abort();
     5395#endif
     5396      double sum =0.0;
     5397      double sum2 =0.0;
     5398      int m=xColumn_;
     5399      double x = info->solution_[m];
     5400      double xB[2];
     5401      double yB[2];
     5402      xB[0]=info->lower_[m];
     5403      xB[1]=info->upper_[m];
     5404      for (int i=0;i<SIZE*SIZE;i+=SIZE) {
     5405        int n = i+SIZE+m;
     5406        double y = info->solution_[n];
     5407        yB[0]=info->lower_[n];
     5408        yB[1]=info->upper_[n];
     5409        int firstLambda=SIZE*SIZE+2*SIZE+4*i+4*m;
     5410        double xyLambda=0.0;
     5411        for (int j=0;j<4;j++) {
     5412          int iX = j>>1;
     5413          int iY = j&1;
     5414          xyLambda += xB[iX]*yB[iY]*info->solution_[firstLambda+j];
     5415        }
     5416        sum += xyLambda*b[i/SIZE];
     5417        double xyTrue = x*y;
     5418        sum2 += xyTrue*b[i/SIZE];
     5419      }
     5420      if (sum>bMax*x+1.0e-5||sum<bMin*x-1.0e-5) {
     5421        //if (sum<bMax*x+1.0e-5&&sum>bMin*x-1.0e-5) {
     5422        printf("bmin*x %g b*w %g bmax*x %g (true) %g\n",bMin*x,sum,bMax*x,sum2);
     5423        printf("m %d lb %g value %g up %g\n",
     5424               m,xB[0],x,xB[1]);
     5425        sum=0.0;
     5426        for (int i=0;i<SIZE*SIZE;i+=SIZE) {
     5427          int n = i+SIZE+m;
     5428          double y = info->solution_[n];
     5429          yB[0]=info->lower_[n];
     5430          yB[1]=info->upper_[n];
     5431          printf("n %d lb %g value %g up %g\n",
     5432                 n,yB[0],y,yB[1]);
     5433          int firstLambda=SIZE*SIZE+2*SIZE+4*i+m*4;
     5434          double xyLambda=0.0;
     5435          for (int j=0;j<4;j++) {
     5436            int iX = j>>1;
     5437            int iY = j&1;
     5438            xyLambda += xB[iX]*yB[iY]*info->solution_[firstLambda+j];
     5439            printf("j %d l %d new xylambda %g ",j,firstLambda+j,xyLambda);
     5440          }
     5441          sum += xyLambda*b[i/SIZE];
     5442          printf(" - sum now %g\n",sum);
     5443        }
     5444      }
     5445      if (sum2>bMax*x+1.0e-5||sum2<bMin*x-1.0e-5) {
     5446        printf("bmin*x %g b*x*y %g bmax*x %g (estimate) %g\n",bMin*x,sum2,bMax*x,sum);
     5447        printf("m %d lb %g value %g up %g\n",
     5448               m,xB[0],x,xB[1]);
     5449        sum2=0.0;
     5450        for (int i=0;i<SIZE*SIZE;i+=SIZE) {
     5451          int n = i+SIZE+m;
     5452          double y = info->solution_[n];
     5453          yB[0]=info->lower_[n];
     5454          yB[1]=info->upper_[n];
     5455          printf("n %d lb %g value %g up %g\n",
     5456                 n,yB[0],y,yB[1]);
     5457          double xyTrue = x*y;
     5458          sum2 += xyTrue*b[i/SIZE];
     5459          printf("xyTrue %g - sum now %g\n",xyTrue,sum2);
     5460        }
     5461      }
     5462    }
     5463  }
    47675464  // If pseudo shadow prices then see what would happen
    47685465  //double pseudoEstimate = 0.0;
     
    47815478      // == row so move x and y not xy
    47825479    }
     5480  }
     5481  if ((branchingStrategy_&16)!=0) {
     5482    // always treat as satisfied!!
     5483    xSatisfied=true;
     5484    ySatisfied=true;
     5485    xyTrue=xyLambda;
    47835486  }
    47845487  if ( !xSatisfied) {
     
    48525555          }
    48535556        }
     5557        if (testCoarse&&(branchingStrategy_&8)!=0&&xB[1]-xB[0]<1.0001*xSatisfied_&&
     5558            yB[1]-yB[0]<1.0001*ySatisfied_)
     5559          feasible=true;
    48545560        if (feasible) {
    48555561          if (xB[1]-xB[0]>=xSatisfied_&&xMeshSize_) {
     
    49385644  }
    49395645  whichWay=whichWay_;
     5646  //if (infeasibility_&&priority_==10)
     5647  //printf("x %d %g %g %g, y %d %g %g %g\n",xColumn_,xB[0],x,xB[1],yColumn_,yB[0],y,yB[1]);
    49405648  return infeasibility_;
    49415649}
    4942 
     5650// Sets infeasibility and other when pseudo shadow prices
     5651void
     5652OsiBiLinear::getPseudoShadow(const OsiBranchingInformation * info)
     5653{
     5654  // order is LxLy, LxUy, UxLy and UxUy
     5655  double xB[2];
     5656  double yB[2];
     5657  xB[0]=info->lower_[xColumn_];
     5658  xB[1]=info->upper_[xColumn_];
     5659  yB[0]=info->lower_[yColumn_];
     5660  yB[1]=info->upper_[yColumn_];
     5661  double x = info->solution_[xColumn_];
     5662  x = CoinMax(x,xB[0]);
     5663  x = CoinMin(x,xB[1]);
     5664  double y = info->solution_[yColumn_];
     5665  y = CoinMax(y,yB[0]);
     5666  y = CoinMin(y,yB[1]);
     5667  int j;
     5668  double xyTrue = x*y;
     5669  double xyLambda = 0.0;
     5670  if ((branchingStrategy_&4)==0) {
     5671    for (j=0;j<4;j++) {
     5672      int iX = j>>1;
     5673      int iY = j&1;
     5674      xyLambda += xB[iX]*yB[iY]*info->solution_[firstLambda_+j];
     5675    }
     5676  } else {
     5677    if (xyRow_>=0) {
     5678      const double * element = info->elementByColumn_;
     5679      const int * row = info->row_;
     5680      const CoinBigIndex * columnStart = info->columnStart_;
     5681      const int * columnLength = info->columnLength_;
     5682      for (j=0;j<4;j++) {
     5683        int iColumn = firstLambda_+j;
     5684        int iStart = columnStart[iColumn];
     5685        int iEnd = iStart + columnLength[iColumn];
     5686        int k=iStart;
     5687        double sol = info->solution_[iColumn];
     5688        for (;k<iEnd;k++) {
     5689          if (xyRow_==row[k])
     5690            xyLambda += element[k]*sol;
     5691        }
     5692      }
     5693    } else {
     5694      // objective
     5695      const double * objective = info->objective_;
     5696      for (j=0;j<4;j++) {
     5697        int iColumn = firstLambda_+j;
     5698        double sol = info->solution_[iColumn];
     5699        xyLambda += objective[iColumn]*sol;
     5700      }
     5701    }
     5702    xyLambda /= coefficient_;
     5703  }
     5704  assert (info->defaultDual_>=0.0);
     5705  // If we move to xy then we move by coefficient * (xyTrue-xyLambda) on row xyRow_
     5706  double movement = xyTrue-xyLambda;
     5707  infeasibility_=0.0;
     5708  const double * pi = info->pi_;
     5709  const double * activity = info->rowActivity_;
     5710  const double * lower = info->rowLower_;
     5711  const double * upper = info->rowUpper_;
     5712  double tolerance = info->primalTolerance_;
     5713  double direction = info->direction_;
     5714  if (xyRow_>=0) {
     5715    assert (!boundType_);
     5716    if (lower[xyRow_]<-1.0e20)
     5717      assert (pi[xyRow_]<=1.0e-3);
     5718    if (upper[xyRow_]>1.0e20)
     5719      assert (pi[xyRow_]>=-1.0e-3);
     5720    double valueP = pi[xyRow_]*direction;
     5721    // if move makes infeasible then make at least default
     5722    double newValue = activity[xyRow_] + movement*coefficient_;
     5723    if (newValue>upper[xyRow_]+tolerance||newValue<lower[xyRow_]-tolerance)
     5724      infeasibility_ += fabs(movement*coefficient_)*CoinMax(fabs(valueP),info->defaultDual_);
     5725  } else {
     5726    // objective
     5727    assert (movement>-1.0e-7);
     5728    infeasibility_ += movement;
     5729  }
     5730  for (int i=0;i<numberExtraRows_;i++) {
     5731    int iRow = extraRow_[i];
     5732    if (lower[iRow]<-1.0e20)
     5733      assert (pi[iRow]<=1.0e-3);
     5734    if (upper[iRow]>1.0e20)
     5735      assert (pi[iRow]>=-1.0e-3);
     5736    double valueP = pi[iRow]*direction;
     5737    // if move makes infeasible then make at least default
     5738    double newValue = activity[iRow] + movement*multiplier_[i];
     5739    if (newValue>upper[iRow]+tolerance||newValue<lower[iRow]-tolerance)
     5740      infeasibility_ += fabs(movement*multiplier_[i])*CoinMax(fabs(valueP),info->defaultDual_);
     5741  }
     5742  if (infeasibility_<1.0e-7)
     5743    infeasibility_=0.0;
     5744  otherInfeasibility_ = CoinMax(1.0e-12,infeasibility_*10.0);
     5745}
    49435746// This looks at solution and sets bounds to contain solution
    49445747double
     
    54766279  const int * row = matrix->getIndices();
    54776280  const CoinBigIndex * columnStart = matrix->getVectorStarts();
    5478   //const int * columnLength = matrix->getVectorLengths();
     6281  const int * columnLength = matrix->getVectorLengths();
    54796282  // order is LxLy, LxUy, UxLy and UxUy
    54806283  double xB[2];
     
    54976300    double y = yB[iY];
    54986301    CoinBigIndex k = columnStart[j+firstLambda_];
     6302    CoinBigIndex last = k+columnLength[j+firstLambda_];
    54996303    double value;
    55006304    // xy
     
    55316335      element[k++]=value;
    55326336      numberUpdated++;
     6337    }
     6338    // Do extra rows
     6339    for (int i=0;i<numberExtraRows_;i++) {
     6340      int iRow = extraRow_[i];
     6341      for (;k<last;k++) {
     6342        if (row[k]==iRow)
     6343          break;
     6344      }
     6345      assert (k<last);
     6346      element[k++] = x*y*multiplier_[i];
    55336347    }
    55346348  }
     
    60996913      if (newUp>upper[iRow]+tolerance||newUp<lower[iRow]-tolerance)
    61006914        u = CoinMax(u,info->defaultDual_);
    6101       upEstimate += u*upMovement;
     6915      upEstimate += u*upMovement*fabs(el2);
    61026916      // if down makes infeasible then make at least default
    61036917      double newDown = activity[iRow] - downMovement*el2;
    61046918      if (newDown>upper[iRow]+tolerance||newDown<lower[iRow]-tolerance)
    61056919        d = CoinMax(d,info->defaultDual_);
    6106       downEstimate += d*downMovement;
     6920      downEstimate += d*downMovement*fabs(el2);
    61076921    }
    61086922    if (downEstimate>=upEstimate) {
     
    66927506  return nelCreate;
    66937507}
    6694 #endif
    66957508#include "ClpConstraint.hpp"
    66967509#include "ClpConstraintLinear.hpp"
     
    67057518                    int mode)
    67067519{
     7520#ifdef COIN_HAS_ASL
     7521  // matrix etc will be changed
     7522  CoinModel coinModel2 = coinModel;
     7523  if (coinModel2.moreInfo()) {
     7524    // for now just ampl objective
     7525    ClpSimplex * model = new ClpSimplex();
     7526    model->loadProblem(coinModel2);
     7527    int numberConstraints;
     7528    ClpConstraint ** constraints=NULL;
     7529    int type = model->loadNonLinear(coinModel2.moreInfo(),
     7530                                    numberConstraints,constraints);
     7531    if (type==1||type==3) {
     7532      int returnCode = model->nonlinearSLP(numberPasses,deltaTolerance);
     7533      assert (!returnCode);
     7534    } else if (type==2||type==4) {
     7535      int returnCode = model->nonlinearSLP(numberConstraints,constraints,
     7536                                           numberPasses,deltaTolerance);
     7537      assert (!returnCode);
     7538    } else {
     7539      printf("error or linear - fix %d\n",type);
     7540    }
     7541    //exit(66);
     7542    return model;
     7543  }
    67077544  // first check and set up arrays
    67087545  int numberColumns = coinModel.numberColumns();
    67097546  int numberRows = coinModel.numberRows();
    67107547  // List of nonlinear rows
    6711   int * which = new int[numberRows+1];
     7548  int * which = new int[numberRows];
    67127549  bool testLinear=false;
    67137550  int numberConstraints=0;
    67147551  int iColumn;
    6715   // matrix etc will be changed
    6716   CoinModel coinModel2 = coinModel;
    67177552  bool linearObjective=true;
    67187553  int maximumQuadraticElements=0;
     
    67497584    for (iColumn=0;iColumn<numberColumns;iColumn++)
    67507585      coinModel2.setObjective(iColumn,0.0);
    6751     which[numberConstraints++]=-1;
    67527586  }
    67537587  int iRow;
     
    67987632  ClpSimplex * model = new ClpSimplex();
    67997633  // return if nothing
    6800   if (!numberConstraints) {
     7634  if (!numberConstraints&&linearObjective) {
    68017635    delete [] which;
    68027636    model->loadProblem(coinModel);
     
    68147648  int saveNumber=numberConstraints;
    68157649  numberConstraints=0;
     7650  ClpQuadraticObjective * quadObj = NULL;
    68167651  if (!linearObjective) {
    68177652    int numberQuadratic=0;
    6818     int largestColumn=-1;
    68197653    CoinZeroN(linearTerm,numberColumns);
    68207654    for (iColumn=0;iColumn<numberColumns;iColumn++) {
     
    68237657      const char * expr = coinModel.getColumnObjectiveAsString(iColumn);
    68247658      if (strcmp(expr,"Numeric")) {
    6825         largestColumn = CoinMax(largestColumn,iColumn);
    68267659        // value*x*y
    68277660        char temp[20000];
     
    68357668          if (jColumn>=0) {
    68367669            columnQuadratic[numberQuadratic]=jColumn;
    6837             elementQuadratic[numberQuadratic++]=2.0*value; // convention
    6838             largestColumn = CoinMax(largestColumn,jColumn);
     7670            if (jColumn!=iColumn)
     7671              elementQuadratic[numberQuadratic++]=1.0*value; // convention
     7672            else if (jColumn==iColumn)
     7673              elementQuadratic[numberQuadratic++]=2.0*value; // convention
    68397674          } else if (jColumn==-2) {
    68407675            linearTerm[iColumn] = value;
    6841             // and put in as row -1
    6842             columnQuadratic[numberQuadratic]=-1;
    6843             elementQuadratic[numberQuadratic++]=value;
    6844             largestColumn = CoinMax(largestColumn,iColumn);
    68457676          } else {
    68467677            printf("bad nonlinear term %s\n",temp);
     
    68527683        // linear part
    68537684        linearTerm[iColumn] = coinModel.getColumnObjective(iColumn);
    6854         // and put in as row -1
    6855         columnQuadratic[numberQuadratic]=-1;
    6856         elementQuadratic[numberQuadratic++]=linearTerm[iColumn];
    6857         if (linearTerm[iColumn])
    6858           largestColumn = CoinMax(largestColumn,iColumn);
    68597685      }
    68607686    }
    68617687    startQuadratic[numberColumns] = numberQuadratic;
    6862     // here we create ClpConstraint
    6863     constraints[numberConstraints++] = new ClpConstraintQuadratic(-1, largestColumn+1, numberColumns,
    6864                                                                   startQuadratic,columnQuadratic,elementQuadratic);
     7688    quadObj = new ClpQuadraticObjective(linearTerm, numberColumns,
     7689                                           startQuadratic,columnQuadratic,elementQuadratic);
    68657690  }
    68667691  int iConstraint;
     
    69497774  delete [] which;
    69507775  model->loadProblem(coinModel2);
    6951   int returnCode = model->nonlinearSLP(numberConstraints, constraints,
    6952                                        numberPasses,deltaTolerance);
    6953   // See if any integers
    6954   for (iConstraint=0;iConstraint<saveNumber;iConstraint++)
    6955     delete constraints[iConstraint];
    6956  
     7776  if (quadObj)
     7777    model->setObjective(quadObj);
     7778  delete quadObj;
     7779  int returnCode;
     7780  if (numberConstraints) {
     7781    returnCode = model->nonlinearSLP(numberConstraints, constraints,
     7782                                     numberPasses,deltaTolerance);
     7783    for (iConstraint=0;iConstraint<saveNumber;iConstraint++)
     7784      delete constraints[iConstraint];
     7785  } else {
     7786    returnCode = model->nonlinearSLP(numberPasses,deltaTolerance);
     7787  }
    69577788  delete [] constraints;
    69587789  assert (!returnCode);
    69597790  return model;
    6960 }
     7791#else
     7792  printf("loadNonLinear needs ampl\n");
     7793  abort();
     7794  return NULL;
     7795#endif
     7796}
     7797OsiChooseStrongSubset::OsiChooseStrongSubset() :
     7798  OsiChooseStrong(),
     7799  numberObjectsToUse_(0)
     7800{
     7801}
     7802
     7803OsiChooseStrongSubset::OsiChooseStrongSubset(const OsiSolverInterface * solver) :
     7804  OsiChooseStrong(solver),
     7805  numberObjectsToUse_(-1)
     7806{
     7807}
     7808
     7809OsiChooseStrongSubset::OsiChooseStrongSubset(const OsiChooseStrongSubset & rhs)
     7810  : OsiChooseStrong(rhs)
     7811
     7812  numberObjectsToUse_ = -1;
     7813}
     7814
     7815OsiChooseStrongSubset &
     7816OsiChooseStrongSubset::operator=(const OsiChooseStrongSubset & rhs)
     7817{
     7818  if (this != &rhs) {
     7819    OsiChooseStrong::operator=(rhs);
     7820    numberObjectsToUse_ = -1;
     7821  }
     7822  return *this;
     7823}
     7824
     7825
     7826OsiChooseStrongSubset::~OsiChooseStrongSubset ()
     7827{
     7828}
     7829
     7830// Clone
     7831OsiChooseVariable *
     7832OsiChooseStrongSubset::clone() const
     7833{
     7834  return new OsiChooseStrongSubset(*this);
     7835}
     7836// Initialize
     7837int
     7838OsiChooseStrongSubset::setupList ( OsiBranchingInformation *info, bool initialize)
     7839{
     7840  assert (solver_==info->solver_);
     7841  // Only has to work with Clp
     7842  OsiSolverInterface * solverA = const_cast<OsiSolverInterface *> (solver_);
     7843  OsiSolverLink * solver = dynamic_cast<OsiSolverLink *> (solverA);
     7844  assert (solver);
     7845  int numberObjects = solver->numberObjects();
     7846  if (numberObjects>numberObjects_) {
     7847    // redo useful arrays
     7848    delete [] upTotalChange_;
     7849    delete [] downTotalChange_;
     7850    delete [] upNumber_;
     7851    delete [] downNumber_;
     7852    numberObjects_ = solver->numberObjects();
     7853    upTotalChange_ = new double [numberObjects_];
     7854    downTotalChange_ = new double [numberObjects_];
     7855    upNumber_ = new int [numberObjects_];
     7856    downNumber_ = new int [numberObjects_];
     7857    CoinZeroN(upTotalChange_,numberObjects_);
     7858    CoinZeroN(downTotalChange_,numberObjects_);
     7859    CoinZeroN(upNumber_,numberObjects_);
     7860    CoinZeroN(downNumber_,numberObjects_);
     7861  }
     7862  if (numberObjectsToUse_<0) {
     7863    // Sort objects so bilinear at end
     7864    OsiObject ** sorted = new OsiObject * [numberObjects];
     7865    OsiObject ** objects = solver->objects();
     7866    numberObjects_=0;
     7867    int numberBiLinear=0;
     7868    int i;
     7869    for (i=0;i<numberObjects;i++) {
     7870      OsiObject * obj = objects[i];
     7871      OsiBiLinear * objB = dynamic_cast<OsiBiLinear *> (obj);
     7872      if (!objB)
     7873        objects[numberObjects_++]=obj;
     7874      else
     7875        sorted[numberBiLinear++]=obj;
     7876    }
     7877    numberObjectsToUse_ = numberObjects_;
     7878    for (i=0;i<numberBiLinear;i++)
     7879      objects[numberObjects_++]=sorted[i];
     7880    delete [] sorted;
     7881    // See if any master objects
     7882    for (i=0;i<numberObjectsToUse_;i++) {
     7883      OsiUsesBiLinear * obj = dynamic_cast<OsiUsesBiLinear *> (objects[i]);
     7884      if (obj)
     7885        obj->addBiLinearObjects(solver);
     7886    }
     7887  }
     7888  solver->setNumberObjects(numberObjectsToUse_);
     7889  numberObjects_=numberObjectsToUse_;
     7890  // Use shadow prices
     7891  info->defaultDual_=0.0;
     7892  int numberUnsatisfied=OsiChooseStrong::setupList ( info, initialize);
     7893  solver->setNumberObjects(numberObjects);
     7894  numberObjects_=numberObjects;
     7895  return numberUnsatisfied;
     7896}
     7897/* Choose a variable
     7898   Returns -
     7899   -1 Node is infeasible
     7900   0  Normal termination - we have a candidate
     7901   1  All looks satisfied - no candidate
     7902   2  We can change the bound on a variable - but we also have a strong branching candidate
     7903   3  We can change the bound on a variable - but we have a non-strong branching candidate
     7904   4  We can change the bound on a variable - no other candidates
     7905   We can pick up branch from whichObject() and whichWay()
     7906   We can pick up a forced branch (can change bound) from whichForcedObject() and whichForcedWay()
     7907   If we have a solution then we can pick up from goodObjectiveValue() and goodSolution()
     7908*/
     7909int
     7910OsiChooseStrongSubset::chooseVariable( OsiSolverInterface * solver, OsiBranchingInformation *info, bool fixVariables)
     7911{
     7912  int numberObjects = solver->numberObjects();
     7913  solver->setNumberObjects(numberObjectsToUse_);
     7914  numberObjects_=numberObjectsToUse_;
     7915  // Use shadow prices
     7916  info->defaultDual_=0.0;
     7917  int returnCode=OsiChooseStrong::chooseVariable(solver,info,fixVariables);
     7918  solver->setNumberObjects(numberObjects);
     7919  numberObjects_=numberObjects;
     7920  return returnCode;
     7921}
     7922/** Default Constructor
     7923
     7924  Equivalent to an unspecified binary variable.
     7925*/
     7926OsiUsesBiLinear::OsiUsesBiLinear ()
     7927  : OsiSimpleInteger(),
     7928    numberBiLinear_(0),
     7929    type_(0),
     7930    objects_(NULL)
     7931{
     7932}
     7933
     7934/** Useful constructor
     7935
     7936  Loads actual upper & lower bounds for the specified variable.
     7937*/
     7938OsiUsesBiLinear::OsiUsesBiLinear (const OsiSolverInterface * solver, int iColumn, int type)
     7939  : OsiSimpleInteger(solver,iColumn),
     7940    numberBiLinear_(0),
     7941    type_(type),
     7942    objects_(NULL)
     7943{
     7944  if (type_) {
     7945    assert(originalLower_==floor(originalLower_+0.5));
     7946    assert(originalUpper_==floor(originalUpper_+0.5));
     7947  }
     7948}
     7949
     7950 
     7951// Useful constructor - passed solver index and original bounds
     7952OsiUsesBiLinear::OsiUsesBiLinear ( int iColumn, double lower, double upper, int type)
     7953  : OsiSimpleInteger(iColumn,lower,upper),
     7954    numberBiLinear_(0),
     7955    type_(type),
     7956    objects_(NULL)
     7957{
     7958  if (type_) {
     7959    assert(originalLower_==floor(originalLower_+0.5));
     7960    assert(originalUpper_==floor(originalUpper_+0.5));
     7961  }
     7962}
     7963
     7964// Useful constructor - passed simple integer
     7965OsiUsesBiLinear::OsiUsesBiLinear ( const OsiSimpleInteger &rhs, int type)
     7966  : OsiSimpleInteger(rhs),
     7967    numberBiLinear_(0),
     7968    type_(type),
     7969    objects_(NULL)
     7970{
     7971  if (type_) {
     7972    assert(originalLower_==floor(originalLower_+0.5));
     7973    assert(originalUpper_==floor(originalUpper_+0.5));
     7974  }
     7975}
     7976
     7977// Copy constructor
     7978OsiUsesBiLinear::OsiUsesBiLinear ( const OsiUsesBiLinear & rhs)
     7979  :OsiSimpleInteger(rhs),
     7980    numberBiLinear_(0),
     7981    type_(rhs.type_),
     7982    objects_(NULL)
     7983
     7984{
     7985}
     7986
     7987// Clone
     7988OsiObject *
     7989OsiUsesBiLinear::clone() const
     7990{
     7991  return new OsiUsesBiLinear(*this);
     7992}
     7993
     7994// Assignment operator
     7995OsiUsesBiLinear &
     7996OsiUsesBiLinear::operator=( const OsiUsesBiLinear& rhs)
     7997{
     7998  if (this!=&rhs) {
     7999    OsiSimpleInteger::operator=(rhs);
     8000    delete [] objects_;
     8001    numberBiLinear_ = 0;
     8002    type_ = rhs.type_;
     8003    objects_ = NULL;
     8004  }
     8005  return *this;
     8006}
     8007
     8008// Destructor
     8009OsiUsesBiLinear::~OsiUsesBiLinear ()
     8010{
     8011  delete [] objects_;
     8012}
     8013// Infeasibility - large is 0.5
     8014double
     8015OsiUsesBiLinear::infeasibility(const OsiBranchingInformation * info, int & whichWay) const
     8016{
     8017  assert (type_==0); // just continuous for now
     8018  double value = info->solution_[columnNumber_];
     8019  value = CoinMax(value, info->lower_[columnNumber_]);
     8020  value = CoinMin(value, info->upper_[columnNumber_]);
     8021  infeasibility_ = 0.0;
     8022  for (int i=0;i<numberBiLinear_;i++) {
     8023    OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (objects_[i]);
     8024    assert (obj);
     8025    obj->getPseudoShadow(info);
     8026    infeasibility_ += objects_[i]->infeasibility();
     8027  }
     8028  bool satisfied=false;
     8029  whichWay=-1;
     8030  if (infeasibility_<=info->integerTolerance_) {
     8031    otherInfeasibility_ = 1.0;
     8032    satisfied=true;
     8033    infeasibility_ = 0.0;
     8034  } else {
     8035    otherInfeasibility_ = 10.0*infeasibility_;
     8036    if (value-info->lower_[columnNumber_]>
     8037        info->upper_[columnNumber_]-value)
     8038      whichWay=1;
     8039    else
     8040      whichWay=-1;
     8041  }
     8042  if (preferredWay_>=0&&!satisfied)
     8043    whichWay = preferredWay_;
     8044  whichWay_=whichWay;
     8045  return infeasibility_;
     8046}
     8047// Creates a branching object
     8048OsiBranchingObject *
     8049OsiUsesBiLinear::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) const
     8050{
     8051  double value = info->solution_[columnNumber_];
     8052  value = CoinMax(value, info->lower_[columnNumber_]);
     8053  value = CoinMin(value, info->upper_[columnNumber_]);
     8054  assert (info->upper_[columnNumber_]>info->lower_[columnNumber_]);
     8055  double nearest = floor(value+0.5);
     8056  double integerTolerance = info->integerTolerance_;
     8057  if (fabs(value-nearest)<integerTolerance) {
     8058    // adjust value
     8059    if (nearest!=info->upper_[columnNumber_])
     8060      value = nearest+2.0*integerTolerance;
     8061    else
     8062      value = nearest-2.0*integerTolerance;
     8063  }
     8064  OsiBranchingObject * branch = new OsiIntegerBranchingObject(solver,this,way,
     8065                                                              value,value,value);
     8066  return branch;
     8067}
     8068// This looks at solution and sets bounds to contain solution
     8069/** More precisely: it first forces the variable within the existing
     8070    bounds, and then tightens the bounds to fix the variable at the
     8071    nearest integer value.
     8072*/
     8073double
     8074OsiUsesBiLinear::feasibleRegion(OsiSolverInterface * solver,
     8075                                 const OsiBranchingInformation * info) const
     8076{
     8077  double value = info->solution_[columnNumber_];
     8078  double newValue = CoinMax(value, info->lower_[columnNumber_]);
     8079  newValue = CoinMin(newValue, info->upper_[columnNumber_]);
     8080  solver->setColLower(columnNumber_,newValue);
     8081  solver->setColUpper(columnNumber_,newValue);
     8082  return fabs(value-newValue);
     8083}
     8084// Add all bi-linear objects
     8085void
     8086OsiUsesBiLinear::addBiLinearObjects(OsiSolverLink * solver)
     8087{
     8088  delete [] objects_;
     8089  numberBiLinear_=0;
     8090  OsiObject ** objects = solver->objects();
     8091  int i;
     8092  int numberObjects = solver->numberObjects();
     8093  for (i=0;i<numberObjects;i++) {
     8094    OsiObject * obj = objects[i];
     8095    OsiBiLinear * objB = dynamic_cast<OsiBiLinear *> (obj);
     8096    if (objB) {
     8097      if (objB->xColumn()==columnNumber_||objB->yColumn()==columnNumber_)
     8098        numberBiLinear_++;
     8099    }
     8100  }
     8101  if (numberBiLinear_) {
     8102    objects_ = new OsiObject * [numberBiLinear_];
     8103    numberBiLinear_=0;
     8104    for (i=0;i<numberObjects;i++) {
     8105      OsiObject * obj = objects[i];
     8106      OsiBiLinear * objB = dynamic_cast<OsiBiLinear *> (obj);
     8107      if (objB) {
     8108        if (objB->xColumn()==columnNumber_||objB->yColumn()==columnNumber_)
     8109          objects_[numberBiLinear_++] = obj;;
     8110      }
     8111    }
     8112  } else {
     8113    objects_=NULL;
     8114  }
     8115}
  • trunk/Cbc/src/CbcLinked.hpp

    r640 r687  
    88   CglTemporary
    99*/
    10 #ifndef COIN_HAS_LINK
    11 #ifdef COIN_HAS_ASL
    12 #define COIN_HAS_LINK
    13 #endif
    14 #endif
    1510#include "CoinModel.hpp"
    16 #ifdef COIN_HAS_LINK
    1711#include "OsiClpSolverInterface.hpp"
     12#include "OsiChooseVariable.hpp"
     13#include "CbcFathom.hpp"
    1814class CbcModel;
    1915class CoinPackedMatrix;
     
    2117class OsiObject;
    2218class CglStored;
    23 //#############################################################################
    24 
    2519/**
    2620   
     
    2923*/
    3024
    31 class OsiSolverLink : public OsiClpSolverInterface {
     25class OsiSolverLink : public CbcOsiSolver {
    3226 
    3327public:
     
    6357      heuristic solution
    6458      Returns solution array
     59      mode -
     60      0 just get continuous
     61      1 round and try normal bab
     62      2 use defaultBound_ to bound integer variables near current solution
    6563  */
    6664  double * heuristicSolution(int numberPasses,double deltaTolerance,int mode);
     
    112110  /// Analyze constraints to see which are convex (quadratic)
    113111  void analyzeObjects();
     112  /// Add reformulated bilinear constraints
     113  void addTighterConstraints();
    114114  /// Objective value of best solution found internally
    115115  inline double bestObjectiveValue() const
     
    163163  inline int integerPriority() const
    164164  { return integerPriority_;};
     165  /// Objective transfer variable if one
     166  inline int objectiveVariable() const
     167  { return objectiveVariable_;}
    165168  /// Set biLinear priority
    166169  inline void setBiLinearPriority(int value)
     
    169172  inline int biLinearPriority() const
    170173  { return biLinearPriority_;};
    171   /// Set Cbc Model
    172   inline void setCbcModel(CbcModel * model)
    173   { cbcModel_=model;};
    174174  /// Return CoinModel
    175175  inline const CoinModel * coinModel() const
     
    177177  /// Set all biLinear priorities on x-x variables
    178178  void setBiLinearPriorities(int value, double meshSize=1.0);
     179  /** Set options and priority on all or some biLinear variables
     180      1 - on I-I
     181      2 - on I-x
     182      4 - on x-x
     183      or combinations.
     184      -1 means leave (for priority value and strategy value)
     185  */
     186  void setBranchingStrategyOnVariables(int strategyValue, int priorityValue=-1,
     187                                       int mode=7);
    179188  /// Set all mesh sizes on x-x variables
    180189  void setMeshSizes(double value);
     
    212221  /// Copy of quadratic model if one
    213222  ClpSimplex * quadraticModel_;
    214   /// Pointer back to CbcModel
    215   CbcModel * cbcModel_;
    216223  /// Number of rows with nonLinearities
    217224  int numberNonLinearRows_;
     
    236243     1 bit (2) - quadratic only in objective (add OA cuts)
    237244     2 bit (4) - convex
    238      4 bit (8) - try adding OA cuts
     245     3 bit (8) - try adding OA cuts
     246     4 bit (16) - add linearized constraints
    239247  */
    240248  int specialOptions2_;
     
    681689  */
    682690  OsiBiLinear (OsiSolverInterface * solver, int xColumn,
     691               int yColumn, int xyRow, double coefficient,
     692               double xMesh, double yMesh,
     693               int numberExistingObjects=0,const OsiObject ** objects=NULL );
     694 
     695  /** Useful constructor -
     696      This Adds in rows and variables to construct valid Linked Ordered Set
     697      Adds extra constraints to match other x/y
     698      So note not const model
     699  */
     700  OsiBiLinear (CoinModel * coinModel, int xColumn,
    683701               int yColumn, int xyRow, double coefficient,
    684702               double xMesh, double yMesh,
     
    793811      next bit
    794812      8 set to say don't use in feasible region
     813      next bit
     814      16 set to say - Always satisfied !!
    795815  */
    796816  inline int branchingStrategy() const
     
    820840  /// Compute lambdas (third entry in each .B is current value) (nonzero if bad)
    821841  double computeLambdas(const double xB[3], const double yB[3],const double xybar[4],double lambda[4]) const;
     842  /// Adds in data for extra row with variable coefficients
     843  void addExtraRow(int row, double multiplier);
     844  /// Sets infeasibility and other when pseudo shadow prices
     845  void getPseudoShadow(const OsiBranchingInformation * info);
    822846
    823847protected:
     
    857881      next bit
    858882      8 set to say don't use in feasible region
     883      next bit
     884      16 set to say - Always satisfied !!
    859885  */
    860886  int branchingStrategy_;
     
    875901  /// Convexity row
    876902  int convexity_;
     903  /// Number of extra rows (coefficients to be modified)
     904  int numberExtraRows_;
     905  /// Multiplier for coefficient on row
     906  double * multiplier_;
     907  /// Row number
     908  int * extraRow_;
    877909  /// Which chosen -1 none, 0 x, 1 y
    878910  mutable short chosen_;
     
    9721004  int numberPoints_;
    9731005};
    974 /// Define a single integer class - but one where you kep branching until fixed even if satsified
     1006/// Define a single integer class - but one where you keep branching until fixed even if satisfied
    9751007
    9761008
     
    10131045  /// data
    10141046 
     1047};
     1048/** Define a single variable class which is involved with OsiBiLinear objects.
     1049    This is used so can make better decision on where to branch as it can look at
     1050    all objects.
     1051
     1052    This version sees if it can re-use code from OsiSimpleInteger
     1053    even if not an integer variable.  If not then need to duplicate code.
     1054*/
     1055
     1056
     1057class OsiUsesBiLinear : public OsiSimpleInteger {
     1058
     1059public:
     1060
     1061  /// Default Constructor
     1062  OsiUsesBiLinear ();
     1063
     1064  /// Useful constructor - passed solver index
     1065  OsiUsesBiLinear (const OsiSolverInterface * solver, int iColumn, int type);
     1066 
     1067  /// Useful constructor - passed solver index and original bounds
     1068  OsiUsesBiLinear (int iColumn, double lower, double upper, int type);
     1069 
     1070  /// Useful constructor - passed simple integer
     1071  OsiUsesBiLinear (const OsiSimpleInteger & rhs, int type);
     1072 
     1073  /// Copy constructor
     1074  OsiUsesBiLinear ( const OsiUsesBiLinear & rhs);
     1075   
     1076  /// Clone
     1077  virtual OsiObject * clone() const;
     1078
     1079  /// Assignment operator
     1080  OsiUsesBiLinear & operator=( const OsiUsesBiLinear& rhs);
     1081
     1082  /// Destructor
     1083  virtual ~OsiUsesBiLinear ();
     1084 
     1085  /// Infeasibility - large is 0.5
     1086  virtual double infeasibility(const OsiBranchingInformation * info, int & whichWay) const;
     1087  /** Creates a branching object
     1088
     1089    The preferred direction is set by \p way, 0 for down, 1 for up.
     1090  */
     1091  virtual OsiBranchingObject * createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) const;
     1092  /** Set bounds to fix the variable at the current value.
     1093
     1094    Given an current value, set the lower and upper bounds to fix the
     1095    variable. Returns amount it had to move variable.
     1096  */
     1097  virtual double feasibleRegion(OsiSolverInterface * solver, const OsiBranchingInformation * info) const;
     1098  /// Add all bi-linear objects
     1099  void addBiLinearObjects(OsiSolverLink * solver);
     1100protected:
     1101  /// data
     1102  /// Number of bilinear objects (maybe could be more general)
     1103  int numberBiLinear_;
     1104  /// Type of variable - 0 continuous, 1 integer
     1105  int type_;
     1106  /// Objects
     1107  OsiObject ** objects_;
     1108};
     1109/** This class chooses a variable to branch on
     1110
     1111    This is just as OsiChooseStrong but it fakes it so only
     1112    first so many are looked at in this phase
     1113
     1114*/
     1115
     1116class OsiChooseStrongSubset  : public OsiChooseStrong {
     1117 
     1118public:
     1119   
     1120  /// Default Constructor
     1121  OsiChooseStrongSubset ();
     1122
     1123  /// Constructor from solver (so we can set up arrays etc)
     1124  OsiChooseStrongSubset (const OsiSolverInterface * solver);
     1125
     1126  /// Copy constructor
     1127  OsiChooseStrongSubset (const OsiChooseStrongSubset &);
     1128   
     1129  /// Assignment operator
     1130  OsiChooseStrongSubset & operator= (const OsiChooseStrongSubset& rhs);
     1131
     1132  /// Clone
     1133  virtual OsiChooseVariable * clone() const;
     1134
     1135  /// Destructor
     1136  virtual ~OsiChooseStrongSubset ();
     1137
     1138  /** Sets up strong list and clears all if initialize is true.
     1139      Returns number of infeasibilities.
     1140      If returns -1 then has worked out node is infeasible!
     1141  */
     1142  virtual int setupList ( OsiBranchingInformation *info, bool initialize);
     1143  /** Choose a variable
     1144      Returns -
     1145     -1 Node is infeasible
     1146     0  Normal termination - we have a candidate
     1147     1  All looks satisfied - no candidate
     1148     2  We can change the bound on a variable - but we also have a strong branching candidate
     1149     3  We can change the bound on a variable - but we have a non-strong branching candidate
     1150     4  We can change the bound on a variable - no other candidates
     1151     We can pick up branch from bestObjectIndex() and bestWhichWay()
     1152     We can pick up a forced branch (can change bound) from firstForcedObjectIndex() and firstForcedWhichWay()
     1153     If we have a solution then we can pick up from goodObjectiveValue() and goodSolution()
     1154     If fixVariables is true then 2,3,4 are all really same as problem changed
     1155  */
     1156  virtual int chooseVariable( OsiSolverInterface * solver, OsiBranchingInformation *info, bool fixVariables);
     1157
     1158  /// Number of objects to use
     1159  inline int numberObjectsToUse() const
     1160  { return numberObjectsToUse_;};
     1161  /// Set number of objects to use
     1162  inline void setNumberObjectsToUse(int value)
     1163  { numberObjectsToUse_ = value;};
     1164
     1165protected:
     1166  // Data
     1167  /// Number of objects to be used (and set in solver)
     1168  int numberObjectsToUse_;
    10151169};
    10161170
     
    11471301  //@}
    11481302};
    1149 #endif
    11501303class ClpSimplex;
    11511304/** Return an approximate solution to a CoinModel.
  • trunk/Cbc/src/CbcMessage.cpp

    r640 r687  
    3131  {CBC_STATUS,10,1,"After %d nodes, %d on tree, %g best solution, best possible %g (%.2f seconds)"},
    3232  {CBC_GAP,11,1,"Exiting as integer gap of %g less than %g or %g%%"},
    33   {CBC_ROUNDING,12,1,"Integer solution of %g found by heuristic after %d iterations and %d nodes (%.2f seconds)"},
     33  {CBC_ROUNDING,12,1,"Integer solution of %g found by %s after %d iterations and %d nodes (%.2f seconds)"},
    3434  {CBC_TREE_SOL,24,1,"Integer solution of %g found by subtree after %d iterations and %d nodes (%.2f seconds)"},
    3535  {CBC_ROOT,13,1,"At root node, %d cuts changed objective from %g to %g in %d passes"},
     
    4949  {CBC_START_SUB,28,1,"Starting sub-tree for %s - maximum nodes %d"},
    5050  {CBC_END_SUB,29,1,"Ending sub-tree for %s"},
    51   {CBC_HEURISTIC_SOLUTION,30,1,"solution of %g found by %s after %2.f seconds"},
     51  {CBC_THREAD_STATS,30,1,"%s%? %d used %d times,  waiting to start %g, %?%g cpu time,%? %g waiting for threads, %? %d locks, %g locked, %g waiting for locks"},
    5252  {CBC_CUTS_STATS,31,1,"%d added rows had average density of %g"},
    5353  {CBC_STRONG_STATS,32,1,"Strong branching done %d times (%d iterations), fathomed %d nodes and fixed %d variables"},
     
    5656  {CBC_HEURISTICS_OFF,36,1,"Heuristics switched off as %d branching objects are of wrong type"},
    5757  {CBC_STATUS2,37,1,"%d nodes, %d on tree, best %g - possible %g depth %d unsat %d its %d (%.2f seconds)"},
     58  {CBC_FPUMP1,38,1,"%s"},
     59  {CBC_FPUMP2,39,2,"%s"},
    5860  {CBC_DUMMY_END,999999,0,""}
    5961};
  • trunk/Cbc/src/CbcMessage.hpp

    r640 r687  
    5454  CBC_START_SUB,
    5555  CBC_END_SUB,
    56   CBC_HEURISTIC_SOLUTION,
     56  CBC_THREAD_STATS,
    5757  CBC_CUTS_STATS,
    5858  CBC_STRONG_STATS,
     
    6161  CBC_HEURISTICS_OFF,
    6262  CBC_STATUS2,
     63  CBC_FPUMP1,
     64  CBC_FPUMP2,
    6365  CBC_DUMMY_END
    6466};
  • trunk/Cbc/src/CbcModel.cpp

    r650 r687  
    1414//#define NODE_LOG
    1515//#define GLOBAL_CUTS_JUST_POINTERS
     16#ifndef CLP_FAST_CODE
     17#ifdef NDEBUG
     18#undef NDEBUG
     19#endif
     20#endif
    1621#include <cassert>
    1722#include <cmath>
     
    3641#include "CbcBranchDynamic.hpp"
    3742#include "CbcHeuristic.hpp"
     43#include "CbcHeuristicFPump.hpp"
    3844#include "CbcModel.hpp"
    3945#include "CbcTreeLocal.hpp"
     
    4854#include "CbcCutGenerator.hpp"
    4955#include "CbcFeasibilityBase.hpp"
     56#include "CbcFathom.hpp"
    5057// include Probing
    5158#include "CglProbing.hpp"
     
    6067#include "CbcCompareActual.hpp"
    6168#include "CbcTree.hpp"
     69//#define CBC_THREAD
     70#ifdef CBC_THREAD
     71#include <pthread.h>
     72#ifndef CLP_FAST_CODE
     73//#define CBC_THREAD_DEBUG 1
     74#endif
     75#ifdef CBC_THREAD_DEBUG
     76#ifdef NDEBUG
     77#undef NDEBUG
     78#undef assert
     79#         define assert(expression) {                              \
     80             if (!(expression)) {                                          \
     81                throw CoinError(__STRING(expression), __PRETTY_FUNCTION__, \
     82                                "", __FILE__, __LINE__);                   \
     83             }                                                             \
     84          }
     85#endif
     86#endif
     87// To Pass across to doOneNode
     88typedef struct {
     89  CbcModel * baseModel;
     90  CbcModel * thisModel;
     91  CbcNode * node; // filled in every time
     92  CbcNode * createdNode; // filled in every time on return
     93  pthread_t threadIdOfBase;
     94  pthread_mutex_t * mutex; // for locking data
     95  pthread_mutex_t * mutex2; // for waking up threads
     96  pthread_cond_t * condition2; // for waking up thread
     97  int returnCode; // -1 available, 0 busy, 1 finished , 2??
     98  double timeLocked;
     99  double timeWaitingToLock;
     100  double timeWaitingToStart;
     101  double timeInThread;
     102  int numberTimesLocked;
     103  int numberTimesUnlocked;
     104  int numberTimesWaitingToStart;
     105  int saveStuff[2];
     106  struct timespec absTime;
     107  bool locked;
     108#if CBC_THREAD_DEBUG
     109  int threadNumber;
     110#endif
     111} threadStruct;
     112static void * doNodesThread(void * voidInfo);
     113static void * doCutsThread(void * voidInfo);
     114#endif
    62115/* Various functions local to CbcModel.cpp */
    63116
     
    449502  system held by the solver.
    450503*/
    451 
    452504void CbcModel::branchAndBound(int doStatistics)
    453505
     
    461513  strongInfo_[2]=0;
    462514  numberStrongIterations_ = 0;
     515#ifndef NDEBUG
     516  {
     517    int i;
     518    int n = solver_->getNumCols();
     519    const double *lower = solver_->getColLower() ;
     520    const double *upper = solver_->getColUpper() ;
     521    for (i=0;i<n;i++) {
     522      assert (lower[i]<1.0e10);
     523      assert (upper[i]>-1.0e10);
     524    }
     525    n = solver_->getNumRows();
     526    lower = solver_->getRowLower() ;
     527    upper = solver_->getRowUpper() ;
     528    for (i=0;i<n;i++) {
     529      assert (lower[i]<1.0e10);
     530      assert (upper[i]>-1.0e10);
     531    }
     532  }
     533#endif
    463534  // original solver (only set if pre-processing)
    464535  OsiSolverInterface * originalSolver=NULL;
     
    466537  OsiObject ** originalObject = NULL;
    467538  // Set up strategies
     539#if 0
     540  std::string problemName ;
     541  solver_->getStrParam(OsiProbName,problemName) ;
     542  if (!strcmp(problemName.c_str(),"PP08A")) solver_->activateRowCutDebugger(problemName.c_str()) ;
     543#endif
    468544  if (strategy_) {
    469545    // May do preprocessing
     
    533609        numberObjects_= numberNewIntegers+numberOldIntegers+numberOldOther+nNonInt;
    534610        object_ = new OsiObject * [numberObjects_];
     611        delete [] integerVariable_;
    535612        integerVariable_ = new int [numberNewIntegers+numberOldIntegers];
    536613        /*
     
    9431020  whichGenerator_ = new int[maximumWhich_] ;
    9441021  memset(whichGenerator_,0,maximumWhich_*sizeof(int));
    945   int currentNumberCuts = 0 ;
    9461022  maximumNumberCuts_ = 0 ;
    9471023  currentNumberCuts_ = 0 ;
     
    10051081  maximumDepthActual_=0;
    10061082  numberDJFixed_=0.0;
     1083  // Do heuristics
     1084  {
     1085    double * newSolution = new double [numberColumns] ;
     1086    double heuristicValue = getCutoff() ;
     1087    int found = -1; // no solution found
     1088
     1089    currentPassNumber_ = 1; // so root heuristics will run
     1090    int i;
     1091    for (i = 0;i<numberHeuristics_;i++) {
     1092      // see if heuristic will do anything
     1093      double saveValue = heuristicValue ;
     1094      int ifSol = heuristic_[i]->solution(heuristicValue,
     1095                                          newSolution);
     1096      if (ifSol>0) {
     1097        // better solution found
     1098        found = i ;
     1099        incrementUsed(newSolution);
     1100        // increment number of solutions so other heuristics can test
     1101        numberSolutions_++;
     1102        numberHeuristicSolutions_++;
     1103      } else {
     1104        heuristicValue = saveValue ;
     1105      }
     1106    }
     1107    currentPassNumber_ = 0;
     1108    /*
     1109      Did any of the heuristics turn up a new solution? Record it before we free
     1110      the vector.
     1111    */
     1112    if (found >= 0) {
     1113      // For compiler error on terra cluster!
     1114      if (found<numberHeuristics_)
     1115        lastHeuristic_ = heuristic_[found];
     1116      else
     1117        lastHeuristic_ = heuristic_[0];
     1118      setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;
     1119      CbcTreeLocal * tree
     1120          = dynamic_cast<CbcTreeLocal *> (tree_);
     1121      if (tree)
     1122        tree->passInSolution(bestSolution_,heuristicValue);
     1123    }
     1124    for (i = 0;i<numberHeuristics_;i++) {
     1125      // delete FPump
     1126      CbcHeuristicFPump * pump
     1127        = dynamic_cast<CbcHeuristicFPump *> (heuristic_[i]);
     1128      if (pump) {
     1129        delete pump;
     1130        numberHeuristics_ --;
     1131        for (int j=i;j<numberHeuristics_;j++)
     1132          heuristic_[j] = heuristic_[j+1];
     1133      }
     1134    }
     1135    delete [] newSolution ;
     1136  }
    10071137  statistics_ = NULL;
    10081138  // Do on switch
     
    10371167        feasible = solveWithCuts(cuts, 1,
    10381168                                 NULL);
    1039 #if 0
    1040         currentNumberCuts_ = cuts.sizeRowCuts();
    1041         if (currentNumberCuts_ >= maximumNumberCuts_) {
    1042           maximumNumberCuts_ = currentNumberCuts;
    1043           delete [] addedCuts_;
    1044           addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
    1045         }
    1046 #endif
    10471169      }
    10481170    }
     
    10641186                           *dblParam_[CbcAllowableFractionGap]);
    10651187  if (bestObjective_-bestPossibleObjective_ < testGap && getCutoffIncrement()>=0.0) {
    1066     if (bestPossibleObjective_<getCutoff()) {
     1188    if (bestPossibleObjective_<getCutoff())
    10671189      stoppedOnGap_ = true ;
    1068       messageHandler()->message(CBC_GAP,messages())
    1069         << bestObjective_-bestPossibleObjective_
    1070         << dblParam_[CbcAllowableGap]
    1071         << dblParam_[CbcAllowableFractionGap]*100.0
    1072         << CoinMessageEol ;
    1073       secondaryStatus_ = 2;
    1074     }
    10751190    feasible = false;
    10761191  }
     
    11271242      }
    11281243    }
    1129 #if 0
    1130     while (anyAction == -1)
    1131     {
    1132       // Set objective value (not so obvious if NLP etc)
    1133       setObjectiveValue(newNode,NULL);
    1134       if (numberBeforeTrust_==0 ) {
    1135         anyAction = newNode->chooseBranch(this,NULL,numberPassesLeft) ;
    1136       } else {
    1137         OsiSolverBranch * branches=NULL;
    1138         anyAction = newNode->chooseDynamicBranch(this,NULL,branches,numberPassesLeft) ;
    1139         if (anyAction==-3)
    1140           anyAction = newNode->chooseBranch(this,NULL,numberPassesLeft) ; // dynamic did nothing
    1141       }
    1142       numberPassesLeft--;
    1143       if (anyAction == -1)
    1144       {
    1145         //if (solverCharacteristics_->solverType()!=4)
    1146           feasible = resolve(NULL,10) != 0 ;
    1147           //else
    1148           //feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_,NULL);
    1149       if (problemFeasibility_->feasible(this,0)<0) {
    1150         feasible=false; // pretend infeasible
    1151       }
    1152       reducedCostFix() ;
    1153       resolved = true ;
    1154 #       ifdef CBC_DEBUG
    1155       printf("Resolve (root) as something fixed, Obj value %g %d rows\n",
    1156              solver_->getObjValue(),
    1157              solver_->getNumRows()) ;
    1158 #       endif
    1159       if (!feasible) anyAction = -2 ; }
    1160       if (anyAction == -2||newNode->objectiveValue() >= cutoff)
    1161       { delete newNode ;
    1162         newNode = NULL ;
    1163         feasible = false ; } }
    1164 #else
    1165   OsiSolverBranch * branches = NULL;
    1166   anyAction = chooseBranch(newNode, numberPassesLeft, NULL, cuts,resolved,
    1167                            NULL,NULL,NULL,branches);
    1168   if (anyAction == -2||newNode->objectiveValue() >= cutoff) {
    1169     if (anyAction != -2) {
    1170       // zap parent nodeInfo
     1244    OsiSolverBranch * branches = NULL;
     1245    anyAction = chooseBranch(newNode, numberPassesLeft, NULL, cuts,resolved,
     1246                             NULL,NULL,NULL,branches);
     1247    if (anyAction == -2||newNode->objectiveValue() >= cutoff) {
     1248      if (anyAction != -2) {
     1249        // zap parent nodeInfo
    11711250#ifdef COIN_DEVELOP
    1172       printf("zapping CbcNodeInfo %x\n",newNode->nodeInfo()->parent());
    1173 #endif
    1174       if (newNode->nodeInfo())
    1175         newNode->nodeInfo()->nullParent();
    1176     }
    1177     delete newNode ;
    1178     newNode = NULL ;
    1179     feasible = false ;
    1180   }
    1181 #endif
     1251        printf("zapping CbcNodeInfo %x\n",newNode->nodeInfo()->parent());
     1252#endif
     1253        if (newNode->nodeInfo())
     1254          newNode->nodeInfo()->nullParent();
     1255      }
     1256      delete newNode ;
     1257      newNode = NULL ;
     1258      feasible = false ;
     1259    }
    11821260  }
    11831261/*
     
    11881266  assert (!newNode || newNode->objectiveValue() <= cutoff) ;
    11891267  // Save address of root node as we don't want to delete it
    1190   CbcNode * rootNode = newNode;
    11911268  // initialize for print out
    11921269  int lastDepth=0;
     
    12041281  a valid solution for use by setBestSolution().
    12051282*/
    1206   CoinWarmStartBasis *lastws = 0 ;
     1283  CoinWarmStartBasis *lastws = NULL ;
    12071284  if (feasible && newNode->branchingObject())
    12081285  { if (resolved)
     
    12911368  */
    12921369  numberLongStrong_=0;
     1370  double totalTime = 0.0;
     1371#ifdef CBC_THREAD
     1372  CbcNode * createdNode=NULL;
     1373  CbcModel ** threadModel = NULL;
     1374  pthread_t * threadId = NULL;
     1375  int * threadCount = NULL;
     1376  pthread_mutex_t mutex;
     1377  pthread_cond_t condition_main;
     1378  pthread_mutex_t condition_mutex;
     1379  pthread_mutex_t * mutex2 = NULL;
     1380  pthread_cond_t * condition2 = NULL;
     1381  threadStruct * threadInfo = NULL;
     1382  bool locked=false;
     1383  int threadStats[6];
     1384  memset(threadStats,0,sizeof(threadStats));
     1385  double timeWaiting=0.0;
     1386  // For now just one model
     1387  if (numberThreads_) {
     1388    threadId = new pthread_t [numberThreads_];
     1389    threadCount = new int [numberThreads_];
     1390    CoinZeroN(threadCount,numberThreads_);
     1391    pthread_mutex_init(&mutex,NULL);
     1392    pthread_cond_init(&condition_main,NULL);
     1393    pthread_mutex_init(&condition_mutex,NULL);
     1394    threadModel = new CbcModel * [numberThreads_+1];
     1395    threadInfo = new threadStruct [numberThreads_+1];
     1396    mutex2 = new pthread_mutex_t [numberThreads_];
     1397    condition2 = new pthread_cond_t [numberThreads_];
     1398    // we don't want a strategy object
     1399    CbcStrategy * saveStrategy = strategy_;
     1400    strategy_ = NULL;
     1401    for (int i=0;i<numberThreads_;i++) {
     1402      pthread_mutex_init(mutex2+i,NULL);
     1403      pthread_cond_init(condition2+i,NULL);
     1404      threadId[i]=0;
     1405      threadInfo[i].baseModel=this;
     1406      threadModel[i]=new CbcModel(*this);
     1407#ifdef COIN_HAS_CLP
     1408      // Solver may need to know about model
     1409      CbcModel * thisModel = threadModel[i];
     1410      CbcOsiSolver * solver =
     1411              dynamic_cast<CbcOsiSolver *>(thisModel->solver()) ;
     1412      if (solver)
     1413        solver->setCbcModel(thisModel);
     1414#endif
     1415      mutex_ = (void *) (threadInfo+i);
     1416      threadModel[i]->moveToModel(this,-1);
     1417      threadInfo[i].thisModel=threadModel[i];
     1418      threadInfo[i].node=NULL;
     1419      threadInfo[i].createdNode=NULL;
     1420      threadInfo[i].threadIdOfBase=pthread_self();
     1421      threadInfo[i].mutex=&mutex;
     1422      threadInfo[i].mutex2=mutex2+i;
     1423      threadInfo[i].condition2=condition2+i;
     1424      threadInfo[i].returnCode=-1;
     1425      threadInfo[i].timeLocked=0.0;
     1426      threadInfo[i].timeWaitingToLock=0.0;
     1427      threadInfo[i].timeWaitingToStart=0.0;
     1428      threadInfo[i].timeInThread=0.0;
     1429      threadInfo[i].numberTimesLocked=0;
     1430      threadInfo[i].numberTimesUnlocked=0;
     1431      threadInfo[i].numberTimesWaitingToStart=0;
     1432      threadInfo[i].locked=false;
     1433#if CBC_THREAD_DEBUG
     1434      threadInfo[i].threadNumber=i+2;
     1435#endif
     1436      pthread_create(threadId+i,NULL,doNodesThread,threadInfo+i);
     1437    }
     1438    strategy_ = saveStrategy;
     1439    // Do a partial one for base model
     1440    threadInfo[numberThreads_].baseModel=this;
     1441    threadModel[numberThreads_]=this;
     1442    mutex_ = (void *) (threadInfo+numberThreads_);
     1443    threadInfo[numberThreads_].node=NULL;
     1444    threadInfo[numberThreads_].mutex=&mutex;
     1445    threadInfo[numberThreads_].condition2=&condition_main;
     1446    threadInfo[numberThreads_].mutex2=&condition_mutex;
     1447    threadInfo[numberThreads_].timeLocked=0.0;
     1448    threadInfo[numberThreads_].timeWaitingToLock=0.0;
     1449    threadInfo[numberThreads_].numberTimesLocked=0;
     1450    threadInfo[numberThreads_].numberTimesUnlocked=0;
     1451    threadInfo[numberThreads_].locked=false;
     1452#if CBC_THREAD_DEBUG
     1453    threadInfo[numberThreads_].threadNumber=1;
     1454#endif
     1455  }
     1456#endif
    12931457/*
    12941458  At last, the actual branch-and-cut search loop, which will iterate until
     
    13011465  than the current objective cutoff.
    13021466*/
    1303   while (!tree_->empty()) {
     1467  while (true) {
     1468#ifdef CBC_THREAD
     1469    if (!locked) {
     1470      lockThread();
     1471      locked=true;
     1472    }
     1473#endif
     1474    if (tree_->empty()) {
     1475#ifdef CBC_THREAD
     1476      if (numberThreads_) {
     1477#ifdef COIN_DEVELOP
     1478        printf("empty\n");
     1479#endif
     1480        // may still be outstanding nodes
     1481        int iThread;
     1482        for (iThread=0;iThread<numberThreads_;iThread++) {
     1483          if (threadId[iThread]) {
     1484            if (threadInfo[iThread].returnCode==0)
     1485              break;
     1486          }
     1487        }
     1488        if (iThread<numberThreads_) {
     1489#ifdef COIN_DEVELOP
     1490          printf("waiting for thread %d code 0\n",iThread);
     1491#endif
     1492          unlockThread();
     1493          locked = false;
     1494          pthread_cond_signal(threadInfo[iThread].condition2); // unlock in case
     1495          while (true) {
     1496            pthread_mutex_lock(&condition_mutex);
     1497            struct timespec absTime;
     1498            clock_gettime(CLOCK_REALTIME,&absTime);
     1499            double time = absTime.tv_sec+1.0e-9*absTime.tv_nsec;
     1500            absTime.tv_nsec += 1000000; // millisecond
     1501            if (absTime.tv_nsec>=1000000000) {
     1502              absTime.tv_nsec -= 1000000000;
     1503              absTime.tv_sec++;
     1504            }
     1505            pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime);
     1506            clock_gettime(CLOCK_REALTIME,&absTime);
     1507            double time2 = absTime.tv_sec+1.0e-9*absTime.tv_nsec;
     1508            timeWaiting += time2-time;
     1509            pthread_mutex_unlock(&condition_mutex);
     1510            if (threadInfo[iThread].returnCode!=0)
     1511              break;
     1512            pthread_cond_signal(threadInfo[iThread].condition2); // unlock
     1513          }
     1514          threadModel[iThread]->moveToModel(this,1);
     1515          assert (threadInfo[iThread].returnCode==1);
     1516          // say available
     1517          threadInfo[iThread].returnCode=-1;
     1518          threadStats[4]++;
     1519#ifdef COIN_DEVELOP
     1520          printf("thread %d code now -1\n",iThread);
     1521#endif
     1522          continue;
     1523        } else {
     1524#ifdef COIN_DEVELOP
     1525          printf("no threads at code 0 \n");
     1526#endif
     1527          // now check if any have just finished
     1528          for (iThread=0;iThread<numberThreads_;iThread++) {
     1529            if (threadId[iThread]) {
     1530              if (threadInfo[iThread].returnCode==1)
     1531                break;
     1532            }
     1533          }
     1534          if (iThread<numberThreads_) {
     1535            unlockThread();
     1536            locked = false;
     1537            threadModel[iThread]->moveToModel(this,1);
     1538            assert (threadInfo[iThread].returnCode==1);
     1539            // say available
     1540            threadInfo[iThread].returnCode=-1;
     1541            threadStats[4]++;
     1542#ifdef COIN_DEVELOP
     1543            printf("thread %d code now -1\n",iThread);
     1544#endif
     1545            continue;
     1546          }
     1547        }
     1548        if (!tree_->empty()) {
     1549#ifdef COIN_DEVELOP
     1550          printf("tree not empty!!!!!!\n");
     1551#endif
     1552          continue;
     1553        }
     1554        for (iThread=0;iThread<numberThreads_;iThread++) {
     1555          if (threadId[iThread]) {
     1556            if (threadInfo[iThread].returnCode!=-1) {
     1557              printf("bad end of tree\n");
     1558              abort();
     1559            }
     1560          }
     1561        }
     1562#ifdef COIN_DEVELOP
     1563        printf("finished ************\n");
     1564#endif
     1565      }
     1566      unlockThread();
     1567      locked=false; // not needed as break
     1568#endif
     1569      break;
     1570    }
     1571#ifdef CBC_THREAD
     1572    unlockThread();
     1573    locked = false;
     1574#endif
     1575/*
     1576  Check for abort on limits: node count, solution count, time, integrality gap.
     1577*/
     1578    totalTime = getCurrentSeconds() ;
     1579    if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
     1580          numberSolutions_ < intParam_[CbcMaxNumSol] &&
     1581          totalTime < dblParam_[CbcMaximumSeconds] &&
     1582          !stoppedOnGap_&&!eventHappened_)) {
     1583      // out of loop
     1584      break;
     1585    }
    13041586#ifdef BONMIN
    13051587    assert(!solverCharacteristics_->solutionAddsCuts() || solverCharacteristics_->mipFeasible());
     
    13081590      double newCutoff = getCutoff();
    13091591      if (analyzeResults_) {
    1310         // see if we could fix any (more)
    1311         int n=0;
    1312         double * newLower = analyzeResults_;
    1313         double * objLower = newLower+numberIntegers_;
    1314         double * newUpper = objLower+numberIntegers_;
    1315         double * objUpper = newUpper+numberIntegers_;
    1316         for (int i=0;i<numberIntegers_;i++) {
    1317           if (objLower[i]>newCutoff) {
    1318             n++;
    1319             if (objUpper[i]>newCutoff) {
    1320               newCutoff = -COIN_DBL_MAX;
    1321               break;
    1322             }
    1323           } else if (objUpper[i]>newCutoff) {
    1324             n++;
    1325           }
    1326         }
    1327         if (newCutoff==-COIN_DBL_MAX) {
    1328           printf("Root analysis says finished\n");
    1329         } else if (n>numberFixedNow_) {
    1330           printf("%d more fixed by analysis - now %d\n",n-numberFixedNow_,n);
    1331           numberFixedNow_=n;
    1332         }
     1592        // see if we could fix any (more)
     1593        int n=0;
     1594        double * newLower = analyzeResults_;
     1595        double * objLower = newLower+numberIntegers_;
     1596        double * newUpper = objLower+numberIntegers_;
     1597        double * objUpper = newUpper+numberIntegers_;
     1598        for (int i=0;i<numberIntegers_;i++) {
     1599          if (objLower[i]>newCutoff) {
     1600            n++;
     1601            if (objUpper[i]>newCutoff) {
     1602              newCutoff = -COIN_DBL_MAX;
     1603              break;
     1604            }
     1605          } else if (objUpper[i]>newCutoff) {
     1606            n++;
     1607          }
     1608        }
     1609        if (newCutoff==-COIN_DBL_MAX) {
     1610          printf("Root analysis says finished\n");
     1611        } else if (n>numberFixedNow_) {
     1612          printf("%d more fixed by analysis - now %d\n",n-numberFixedNow_,n);
     1613          numberFixedNow_=n;
     1614        }
    13331615      }
    13341616      if (eventHandler) {
    1335         if (!eventHandler->event(CbcEventHandler::solution)) {
    1336           eventHappened_=true; // exit
    1337         }
    1338       }
     1617        if (!eventHandler->event(CbcEventHandler::solution)) {
     1618          eventHappened_=true; // exit
     1619        }
     1620      }
     1621      lockThread();
    13391622      // Do from deepest
    13401623      tree_->cleanTree(this, newCutoff,bestPossibleObjective_) ;
    13411624      nodeCompare_->newSolution(this) ;
    13421625      nodeCompare_->newSolution(this,continuousObjective_,
    1343                                 continuousInfeasibilities_) ;
     1626                                continuousInfeasibilities_) ;
    13441627      tree_->setComparison(*nodeCompare_) ;
    1345       if (tree_->empty())
    1346         break; // finished
     1628      if (tree_->empty()) {
     1629        unlockThread();
     1630        // For threads we need to check further
     1631        //break; // finished
     1632        continue;
     1633      }
     1634      unlockThread();
    13471635    }
    13481636    cutoff = getCutoff() ;
    13491637/*
    1350   Periodic activities: Opportunities to
     1638    Periodic activities: Opportunities to
    13511639    + tweak the nodeCompare criteria,
    13521640    + check if we've closed the integrality gap enough to quit,
     
    13541642*/
    13551643    if ((numberNodes_%1000) == 0) {
     1644      lockThread();
    13561645      bool redoTree=nodeCompare_->every1000Nodes(this, numberNodes_) ;
    13571646#ifdef CHECK_CUT_SIZE
     
    13611650      if (redoTree)
    13621651        tree_->setComparison(*nodeCompare_) ;
     1652      unlockThread();
    13631653    }
    13641654    if (saveCompare&&!hotstartSolution_) {
     
    13681658      saveCompare=NULL;
    13691659      // redo tree
     1660      lockThread();
    13701661      tree_->setComparison(*nodeCompare_) ;
     1662      unlockThread();
    13711663    }
    13721664    if ((numberNodes_%printFrequency_) == 0) {
     1665      lockThread();
    13731666      int j ;
    13741667      int nNodes = tree_->size() ;
     
    13791672          bestPossibleObjective_ = node->objectiveValue() ;
    13801673      }
     1674      unlockThread();
    13811675      if (!intParam_[CbcPrinting]) {
    13821676        messageHandler()->message(CBC_STATUS,messages())
     
    14051699      stoppedOnGap_ = true ;
    14061700    }
    1407 
     1701   
    14081702#   ifdef CHECK_NODE_FULL
    14091703    verifyTreeNodes(tree_,*this) ;
     
    14121706    verifyCutCounts(tree_,*this) ;
    14131707#   endif
    1414 
    14151708/*
    14161709  Now we come to the meat of the loop. To create the active subproblem, we'll
     
    14231716    if (!node)
    14241717      continue;
     1718#ifndef CBC_THREAD
     1719    int currentNumberCuts = 0 ;
    14251720    currentNode_=node; // so can be accessed elsewhere
    14261721#ifdef CBC_DEBUG
     
    14291724           node->guessedObjectiveValue());
    14301725#endif
     1726#if NEW_UPDATE_OBJECT==0
    14311727    // Save clone in branching decision
    14321728    if(branchingMethod_)
    14331729      branchingMethod_->saveBranchingObject(node->modifiableBranchingObject());
    1434     bool nodeOnTree=false; // Node has been popped
     1730#endif
    14351731    // Say not on optimal path
    14361732    bool onOptimalPath=false;
     
    14601756    CbcNodeInfo * nodeInfo = node->nodeInfo() ;
    14611757    newNode = NULL ;
     1758    int branchesLeft=0;
    14621759    if (!addCuts(node,lastws,numberFixedNow_>numberFixedAtRoot_))
    14631760    { int i ;
     
    14711768        solverCharacteristics_->setBeforeUpper(upperBefore);
    14721769      }
    1473       bool deleteNode ;
    14741770      if (messageHandler()->logLevel()>2)
    14751771        node->modifiableBranchingObject()->print();
    1476       int retCode;
    14771772      if (!useOsiBranching)
    1478         retCode = node->branch(NULL); // old way
     1773        branchesLeft = node->branch(NULL); // old way
    14791774      else
    1480         retCode = node->branch(solver_); // new way
    1481       if (retCode)
    1482       {
     1775        branchesLeft = node->branch(solver_); // new way
     1776      if (branchesLeft) {
    14831777        // set nodenumber correctly
    14841778        node->nodeInfo()->setNodeNumber(numberNodes2_);
     
    14971791        }
    14981792        numberNodes2_++;
    1499         nodeOnTree=true; // back on tree
    1500         deleteNode = false ;
     1793        //nodeOnTree=true; // back on tree
     1794        //deleteNode = false ;
    15011795#       ifdef CHECK_NODE
    15021796        printf("Node %x pushed back on tree - %d left, %d count\n",node,
     
    15041798               nodeInfo->numberPointingToThis()) ;
    15051799#       endif
    1506       }
    1507       else
    1508       { deleteNode = true ;
    1509       if (!nodeInfo->numberBranchesLeft())
    1510         nodeInfo->allBranchesGone(); // can clean up
    1511       }
    1512 
     1800      } else {
     1801        //deleteNode = true ;
     1802        if (!nodeInfo->numberBranchesLeft())
     1803          nodeInfo->allBranchesGone(); // can clean up
     1804      }
    15131805      if ((specialOptions_&1)!=0) {
    15141806        /*
     
    15241816        }
    15251817      }
     1818     
    15261819/*
    15271820  Reoptimize, possibly generating cuts and/or using heuristics to find
     
    15971890      }
    15981891/*
    1599   Check for abort on limits: node count, solution count, time, integrality gap.
    1600 */
    1601       double totalTime = CoinCpuTime()-dblParam_[CbcStartSeconds] ;
    1602       if (numberNodes_ < intParam_[CbcMaxNumNode] &&
    1603           numberSolutions_ < intParam_[CbcMaxNumSol] &&
    1604           totalTime < dblParam_[CbcMaximumSeconds] &&
    1605           !stoppedOnGap_&&!eventHappened_)
    1606       {
    1607 /*
    16081892  Are we still feasible? If so, create a node and do the work to attach a
    16091893  branching object, reoptimising as needed if chooseBranch() identifies
     
    16451929          int numberPassesLeft=5;
    16461930          checkingNode=true;
    1647 #if 0
    1648           while (anyAction == -1)
    1649           {
    1650             // Set objective value (not so obvious if NLP etc)
    1651             setObjectiveValue(newNode,node);
    1652             if (numberBeforeTrust_==0 ) {
    1653               anyAction = newNode->chooseBranch(this,node,numberPassesLeft) ;
    1654             } else {
    1655               OsiSolverBranch * branches=NULL;
    1656               anyAction = newNode->chooseDynamicBranch(this,node,branches,numberPassesLeft) ;
    1657               if (anyAction==-3)
    1658                 anyAction = newNode->chooseBranch(this,node,numberPassesLeft) ; // dynamic did nothing
    1659             }
    1660             if (solverCharacteristics_ &&
    1661                 solverCharacteristics_->solutionAddsCuts() && // we are in some OA based bab
    1662                 feasible && (newNode->numberUnsatisfied()==0) //solution has become integer feasible during strong branching
    1663                 )
    1664               {
    1665                 //in the present case we need to check here integer infeasibility if the node is not fathomed we will have to do the loop
    1666                 // again
    1667                 std::cout<<solver_<<std::endl;
    1668                 resolve(solver_);
    1669                 double objval = solver_->getObjValue();
    1670                 setBestSolution(CBC_SOLUTION, objval,
    1671                                 solver_->getColSolution()) ;
    1672                 lastHeuristic_ = NULL;
    1673                 int easy=2;
    1674                 if (!solverCharacteristics_->mipFeasible())//did we prove that the node could be pruned?
    1675                   feasible = false;
    1676                 // Reset the bound now
    1677                 solverCharacteristics_->setMipBound(-COIN_DBL_MAX);
    1678                
    1679                
    1680                 //numberPassesLeft++;
    1681                 solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
    1682                 feasible &= resolve(node ? node->nodeInfo() : NULL,11) != 0 ;
    1683                 solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
    1684                 resolved = true ;
    1685                 if (problemFeasibility_->feasible(this,0)<0) {
    1686                   feasible=false; // pretend infeasible
    1687                 }
    1688                 if(feasible)
    1689                   anyAction = -1;
    1690                 else
    1691                   anyAction = -2;
    1692               }
    1693 /*
    1694   Yep, false positives for sure. And no easy way to distinguish honest
    1695   infeasibility from `found a solution and tightened objective target.'
    1696 
    1697             if (onOptimalPath)
    1698               assert (anyAction!=-2); // can be useful but gives false positives on strong
    1699 */
    1700             numberPassesLeft--;
    1701             if (numberPassesLeft<=-1) {
    1702               if (!numberLongStrong_)
    1703                 messageHandler()->message(CBC_WARNING_STRONG,
    1704                                           messages()) << CoinMessageEol ;
    1705               numberLongStrong_++;
    1706             }
    1707             if (anyAction == -1)
    1708             {
    1709               // can do quick optimality check
    1710               int easy=2;
    1711               solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
    1712               //if (solverCharacteristics_->solverType()!=4)
    1713                 feasible = resolve(node ? node->nodeInfo() : NULL,11) != 0 ;
    1714                 //else
    1715                 //feasible = solveWithCuts(cuts,maximumCutPasses_,node);
    1716               solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
    1717               resolved = true ;
    1718               if (problemFeasibility_->feasible(this,0)<0) {
    1719                 feasible=false; // pretend infeasible
    1720               }
    1721               if (feasible)
    1722               {
    1723                 // Set objective value (not so obvious if NLP etc)
    1724                 setObjectiveValue(newNode,node);
    1725                 reducedCostFix() ;
    1726                 if (newNode->objectiveValue() >= getCutoff())
    1727                   anyAction=-2;
    1728               }
    1729               else
    1730               { anyAction = -2 ; } } }
    1731           if (anyAction >= 0)
    1732           { if (resolved)
    1733             { bool needValidSolution = (newNode->branchingObject() == NULL) ;
    1734               takeOffCuts(cuts,needValidSolution,NULL) ;
    1735 #             ifdef CHECK_CUT_COUNTS
    1736               { printf("Number of rows after chooseBranch fix (node)"
    1737                        "(active only) %d\n",
    1738                         numberRowsAtContinuous_+numberNewCuts_+
    1739                         numberOldActiveCuts_) ;
    1740                 const CoinWarmStartBasis* debugws =
    1741                   dynamic_cast<const CoinWarmStartBasis*>
    1742                     (solver_->getWarmStart()) ;
    1743                 debugws->print() ;
    1744                 delete debugws ; }
    1745 #             endif
    1746             }
    1747             newNode->createInfo(this,node,lastws,lowerBefore,upperBefore,
    1748                                 numberOldActiveCuts_,numberNewCuts_) ;
    1749             if (newNode->numberUnsatisfied()) {
    1750               maximumDepthActual_ = CoinMax(maximumDepthActual_,newNode->depth());
    1751               newNode->initializeInfo() ;
    1752               newNode->nodeInfo()->addCuts(cuts,newNode->numberBranches(),
    1753                                            whichGenerator_) ; } } }
    1754         else {
    1755           anyAction = -2 ;
    1756           // Reset bound anyway (no harm if not odd)
    1757           solverCharacteristics_->setMipBound(-COIN_DBL_MAX);
    1758           node->nodeInfo()->decrement();
    1759         }
    1760         // May have slipped through i.e. anyAction == 0 and objective above cutoff
    1761         // I think this will screw up cut reference counts if executed.
    1762         // We executed addCuts just above. (lh)
    1763         if ( anyAction >=0 ) {
    1764           assert (newNode);
    1765           if (newNode->objectiveValue() >= getCutoff()) {
    1766             anyAction = -2; // say bad after all
    1767 #ifdef COIN_DEVELOP
    1768             printf("zapping2 CbcNodeInfo %x\n",newNode->nodeInfo()->parent());
    1769 #endif
    1770             // zap parent nodeInfo
    1771             if (newNode->nodeInfo())
    1772               newNode->nodeInfo()->nullParent();
    1773           }
    1774         }
    1775 /*
    1776   If we end up infeasible, we can delete the new node immediately. Since this
    1777   node won't be needing the cuts we collected, decrement the reference counts.
    1778   If we are feasible, then we'll be placing this node into the live set, so
    1779   increment the reference count in the current (parent) nodeInfo.
    1780 */
    1781         if (anyAction == -2)
    1782         { delete newNode ;
    1783           newNode = NULL ;
    1784           // say strong doing well
    1785           if (checkingNode)
    1786             setSpecialOptions(specialOptions_|8);
    1787           for (i = 0 ; i < currentNumberCuts_ ; i++)
    1788           { if (addedCuts_[i])
    1789             { if (!addedCuts_[i]->decrement(1))
    1790                 delete addedCuts_[i] ; } } }
    1791         else
    1792         { nodeInfo->increment() ;
    1793         if ((numberNodes_%20)==0) {
    1794           // say strong not doing as well
    1795           setSpecialOptions(specialOptions_&~8);
    1796         }
    1797 #else
    17981931        OsiSolverBranch * branches=NULL;
    17991932        // point to useful information
     
    18231956          }
    18241957        }
    1825 #endif
    18261958        }
    18271959/*
     
    18661998                generate=false; // only special cuts
    18671999              if (generate) {
    1868                 generator_[i]->generateCuts(theseCuts,true,NULL) ;
     2000                generator_[i]->generateCuts(theseCuts,true,solver_,NULL) ;
    18692001                int numberRowCutsAfter = theseCuts.sizeRowCuts() ;
    18702002                if (numberRowCutsAfter) {
     
    19132045                heurValue = saveValue ; } }
    19142046            if (found >= 0) {
     2047              lastHeuristic_ = heuristic_[found];
    19152048              setBestSolution(CBC_ROUNDING,heurValue,newSolution) ;
    1916               lastHeuristic_ = heuristic_[found];
    19172049            }
    19182050            delete [] newSolution ;
     
    19642096  set.
    19652097*/
    1966         if (deleteNode)
    1967           delete node ; }
    1968 /*
    1969   End of the non-abort actions. The next block of code is executed if we've
    1970   aborted because we hit one of the limits. Clean up by deleting the live set
    1971   and break out of the node processing loop. Note that on an abort, node may
    1972   have been pushed back onto the tree for further processing, in which case
    1973   it'll be deleted in cleanTree. We need to check.
    1974 */
     2098      if (branchesLeft)
     2099      {
     2100      }
    19752101      else
    19762102      {
    1977         if (tree_->size())
    1978           tree_->cleanTree(this,-COIN_DBL_MAX,bestPossibleObjective_) ;
    1979         delete nextRowCut_;
    1980         // We need to get rid of node if is has already been popped from tree
    1981         if (!nodeOnTree&&!stoppedOnGap_&&node!=rootNode)
    1982           delete node;
    1983         if (stoppedOnGap_)
    1984         { messageHandler()->message(CBC_GAP,messages())
    1985             << bestObjective_-bestPossibleObjective_
    1986             << dblParam_[CbcAllowableGap]
    1987             << dblParam_[CbcAllowableFractionGap]*100.0
    1988             << CoinMessageEol ;
    1989         secondaryStatus_ = 2;
    1990           status_ = 0 ; }
    1991         else
    1992         if (isNodeLimitReached())
    1993         { handler_->message(CBC_MAXNODES,messages_) << CoinMessageEol ;
    1994         secondaryStatus_ = 3;
    1995           status_ = 1 ; }
    1996         else
    1997         if (totalTime >= dblParam_[CbcMaximumSeconds])
    1998         { handler_->message(CBC_MAXTIME,messages_) << CoinMessageEol ;
    1999         secondaryStatus_ = 4;
    2000           status_ = 1 ; }
    2001         else
    2002         if (eventHappened_)
    2003         { handler_->message(CBC_EVENT,messages_) << CoinMessageEol ;
    2004         secondaryStatus_ = 5;
    2005           status_ = 5 ; }
    2006         else
    2007         { handler_->message(CBC_MAXSOLS,messages_) << CoinMessageEol ;
    2008         secondaryStatus_ = 6;
    2009           status_ = 1 ; }
    2010         break ; }
     2103        if (!nodeInfo->numberBranchesLeft())
     2104          nodeInfo->allBranchesGone(); // can clean up
     2105        delete node ; }
     2106    } else {
     2107      // add cuts found to be infeasible (on bound)!
     2108      abort();
     2109      delete node;
     2110    }
    20112111/*
    20122112  Delete cuts to get back to the original system.
     
    20222122        { delRows[i] = i+numberRowsAtContinuous_ ; }
    20232123        solver_->deleteRows(numberToDelete,delRows) ;
    2024         delete [] delRows ; } }
    2025 /*
    2026   This node fathomed when addCuts atttempted to revive it. Toss it.
    2027 */
    2028     else
    2029       { delete node ; } }
     2124        delete [] delRows ; }
     2125#else // end of not CBC_THREAD
     2126      if (!numberThreads_) {
     2127        doOneNode(this,node,createdNode);
     2128      } else {
     2129        threadStats[0]++;
     2130        //need to think
     2131        int iThread;
     2132        // Start one off if any available
     2133        for (iThread=0;iThread<numberThreads_;iThread++) {
     2134          if (threadInfo[iThread].returnCode==-1)
     2135            break;
     2136        }
     2137        if (iThread<numberThreads_) {
     2138          threadInfo[iThread].node=node;
     2139          assert (threadInfo[iThread].returnCode==-1);
     2140          // say in use
     2141          threadInfo[iThread].returnCode=0;
     2142          threadModel[iThread]->moveToModel(this,0);
     2143          pthread_cond_signal(threadInfo[iThread].condition2); // unlock
     2144          threadCount[iThread]++;
     2145        }
     2146        lockThread();
     2147        locked=true;
     2148        // see if any finished
     2149        for (iThread=0;iThread<numberThreads_;iThread++) {
     2150          if (threadInfo[iThread].returnCode>0)
     2151            break;
     2152        }
     2153        unlockThread();
     2154        locked=false;
     2155        if (iThread<numberThreads_) {
     2156          threadModel[iThread]->moveToModel(this,1);
     2157          assert (threadInfo[iThread].returnCode==1);
     2158          // say available
     2159          threadInfo[iThread].returnCode=-1;
     2160          // carry on
     2161          threadStats[3]++;
     2162        } else {
     2163          // Start one off if any available
     2164          for (iThread=0;iThread<numberThreads_;iThread++) {
     2165            if (threadInfo[iThread].returnCode==-1)
     2166              break;
     2167          }
     2168          if (iThread<numberThreads_) {
     2169            lockThread();
     2170            locked=true;
     2171            // If any on tree get
     2172            if (!tree_->empty()) {
     2173              //node = tree_->bestNode(cutoff) ;
     2174              //assert (node);
     2175              threadStats[1]++;
     2176              continue; // ** get another node
     2177            }
     2178            unlockThread();
     2179            locked=false;
     2180          }
     2181          // wait (for debug could sleep and use test)
     2182          bool finished=false;
     2183          while (!finished) {
     2184            pthread_mutex_lock(&condition_mutex);
     2185            struct timespec absTime;
     2186            clock_gettime(CLOCK_REALTIME,&absTime);
     2187            double time = absTime.tv_sec+1.0e-9*absTime.tv_nsec;
     2188            absTime.tv_nsec += 1000000; // millisecond
     2189            if (absTime.tv_nsec>=1000000000) {
     2190              absTime.tv_nsec -= 1000000000;
     2191              absTime.tv_sec++;
     2192            }
     2193            pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime);
     2194            clock_gettime(CLOCK_REALTIME,&absTime);
     2195            double time2 = absTime.tv_sec+1.0e-9*absTime.tv_nsec;
     2196            timeWaiting += time2-time;
     2197            pthread_mutex_unlock(&condition_mutex);
     2198            for (iThread=0;iThread<numberThreads_;iThread++) {
     2199              if (threadInfo[iThread].returnCode>0) {
     2200                finished=true;
     2201                break;
     2202              } else if (threadInfo[iThread].returnCode==0) {
     2203                pthread_cond_signal(threadInfo[iThread].condition2); // unlock
     2204              }
     2205            }
     2206          }
     2207          assert (iThread<numberThreads_);
     2208          threadModel[iThread]->moveToModel(this,1);
     2209          node = threadInfo[iThread].node;
     2210          threadInfo[iThread].node=NULL;
     2211          assert (threadInfo[iThread].returnCode==1);
     2212          // say available
     2213          threadInfo[iThread].returnCode=-1;
     2214          // carry on
     2215          threadStats[2]++;
     2216        }
     2217      }
     2218      //lastDepth=node->depth();
     2219      //lastUnsatisfied=node->numberUnsatisfied();
     2220#endif // end of CBC_THREAD
     2221  }
     2222#ifdef CBC_THREAD
     2223  if (numberThreads_) {
     2224    //printf("stats ");
     2225    //for (unsigned int j=0;j<sizeof(threadStats)/sizeof(int);j++)
     2226    //printf("%d ",threadStats[j]);
     2227    //printf("\n");
     2228    int i;
     2229    // Seems to be bug in CoinCpu on Linux - does threads as well despite documentation
     2230    double time=0.0;
     2231    for (i=0;i<numberThreads_;i++)
     2232      time += threadInfo[i].timeInThread;
     2233    bool goodTimer = time<(getCurrentSeconds());
     2234    //bool stopped = (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
     2235    //        numberSolutions_ < intParam_[CbcMaxNumSol] &&
     2236    //        totalTime < dblParam_[CbcMaximumSeconds] &&
     2237    //        !stoppedOnGap_&&!eventHappened_));
     2238    for (i=0;i<numberThreads_;i++) {
     2239      while (threadInfo[i].returnCode==0) {
     2240        pthread_cond_signal(threadInfo[i].condition2); // unlock
     2241        pthread_mutex_lock(&condition_mutex);
     2242        struct timespec absTime;
     2243        clock_gettime(CLOCK_REALTIME,&absTime);
     2244        absTime.tv_nsec += 1000000; // millisecond
     2245        if (absTime.tv_nsec>=1000000000) {
     2246          absTime.tv_nsec -= 1000000000;
     2247          absTime.tv_sec++;
     2248        }
     2249        pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime);
     2250        clock_gettime(CLOCK_REALTIME,&absTime);
     2251        pthread_mutex_unlock(&condition_mutex);
     2252      }
     2253      pthread_cond_signal(threadInfo[i].condition2); // unlock
     2254      pthread_mutex_lock(&condition_mutex); // not sure necessary but have had one hang on interrupt
     2255      threadModel[i]->numberThreads_=0; // say exit
     2256      threadInfo[i].returnCode=0;
     2257      pthread_mutex_unlock(&condition_mutex);
     2258      pthread_cond_signal(threadInfo[i].condition2); // unlock
     2259      //if (!stopped)
     2260        pthread_join(threadId[i],NULL);
     2261        //else
     2262        //pthread_kill(threadId[i]); // kill rather than try and synchronize
     2263      threadModel[i]->moveToModel(this,2);
     2264      pthread_mutex_destroy (threadInfo[i].mutex2);
     2265      pthread_cond_destroy (threadInfo[i].condition2);
     2266      assert (threadInfo[i].numberTimesLocked==threadInfo[i].numberTimesUnlocked);
     2267      handler_->message(CBC_THREAD_STATS,messages_)
     2268        <<"Thread";
     2269      handler_->printing(true)
     2270        <<i<<threadCount[i]<<threadInfo[i].timeWaitingToStart;
     2271      handler_->printing(goodTimer)<<threadInfo[i].timeInThread;
     2272      handler_->printing(false)<<0.0;
     2273      handler_->printing(true)<<threadInfo[i].numberTimesLocked
     2274        <<threadInfo[i].timeLocked<<threadInfo[i].timeWaitingToLock
     2275        <<CoinMessageEol;
     2276    }
     2277    assert (threadInfo[numberThreads_].numberTimesLocked==threadInfo[numberThreads_].numberTimesUnlocked);
     2278    handler_->message(CBC_THREAD_STATS,messages_)
     2279      <<"Main thread";
     2280    handler_->printing(false)<<0<<0<<0.0;
     2281    handler_->printing(false)<<0.0;
     2282    handler_->printing(true)<<timeWaiting;
     2283    handler_->printing(true)<<threadInfo[numberThreads_].numberTimesLocked
     2284      <<threadInfo[numberThreads_].timeLocked<<threadInfo[numberThreads_].timeWaitingToLock
     2285      <<CoinMessageEol;
     2286    pthread_mutex_destroy (&mutex);
     2287    pthread_cond_destroy (&condition_main);
     2288    pthread_mutex_destroy (&condition_mutex);
     2289    // delete models (here in case some point to others)
     2290    for (i=0;i<numberThreads_;i++) {
     2291      delete threadModel[i];
     2292    }
     2293    delete [] mutex2;
     2294    delete [] condition2;
     2295    delete [] threadId;
     2296    delete [] threadInfo;
     2297    delete [] threadModel;
     2298    delete [] threadCount;
     2299    mutex_=NULL;
     2300    // adjust time to allow for children on some systems
     2301    dblParam_[CbcStartSeconds] -= CoinCpuTimeJustChildren();
     2302  }
     2303#endif
     2304/*
     2305  End of the non-abort actions. The next block of code is executed if we've
     2306  aborted because we hit one of the limits. Clean up by deleting the live set
     2307  and break out of the node processing loop. Note that on an abort, node may
     2308  have been pushed back onto the tree for further processing, in which case
     2309  it'll be deleted in cleanTree. We need to check.
     2310*/
     2311    if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
     2312        numberSolutions_ < intParam_[CbcMaxNumSol] &&
     2313        totalTime < dblParam_[CbcMaximumSeconds] &&
     2314        !stoppedOnGap_&&!eventHappened_)) {
     2315      if (tree_->size())
     2316        tree_->cleanTree(this,-COIN_DBL_MAX,bestPossibleObjective_) ;
     2317      delete nextRowCut_;
     2318      if (stoppedOnGap_)
     2319        { messageHandler()->message(CBC_GAP,messages())
     2320          << bestObjective_-bestPossibleObjective_
     2321          << dblParam_[CbcAllowableGap]
     2322          << dblParam_[CbcAllowableFractionGap]*100.0
     2323          << CoinMessageEol ;
     2324        secondaryStatus_ = 2;
     2325        status_ = 0 ; }
     2326        else
     2327          if (isNodeLimitReached())
     2328            { handler_->message(CBC_MAXNODES,messages_) << CoinMessageEol ;
     2329            secondaryStatus_ = 3;
     2330            status_ = 1 ; }
     2331          else
     2332        if (totalTime >= dblParam_[CbcMaximumSeconds])
     2333          { handler_->message(CBC_MAXTIME,messages_) << CoinMessageEol ;
     2334          secondaryStatus_ = 4;
     2335          status_ = 1 ; }
     2336        else
     2337          if (eventHappened_)
     2338            { handler_->message(CBC_EVENT,messages_) << CoinMessageEol ;
     2339            secondaryStatus_ = 5;
     2340            status_ = 5 ; }
     2341          else
     2342            { handler_->message(CBC_MAXSOLS,messages_) << CoinMessageEol ;
     2343            secondaryStatus_ = 6;
     2344            status_ = 1 ; }
     2345    }
    20302346/*
    20312347  That's it, we've exhausted the search tree, or broken out of the loop because
     
    24642780  subTreeModel_(NULL),
    24652781  numberStoppedSubTrees_(0),
     2782  mutex_(NULL),
    24662783  presolve_(0),
    24672784  numberStrong_(5),
     
    25112828  searchStrategy_(-1),
    25122829  numberStrongIterations_(0),
    2513   resolveAfterTakeOffCuts_(true)
     2830  resolveAfterTakeOffCuts_(true),
     2831#if NEW_UPDATE_OBJECT>1
     2832  numberUpdateItems_(0),
     2833  maximumNumberUpdateItems_(0),
     2834  updateItems_(NULL),
     2835#endif
     2836  numberThreads_(0),
     2837  threadMode_(0)
    25142838{
    25152839  memset(intParam_,0,sizeof(intParam_));
     
    25942918  subTreeModel_(NULL),
    25952919  numberStoppedSubTrees_(0),
     2920  mutex_(NULL),
    25962921  presolve_(0),
    25972922  numberStrong_(5),
     
    26412966  searchStrategy_(-1),
    26422967  numberStrongIterations_(0),
    2643   resolveAfterTakeOffCuts_(true)
     2968  resolveAfterTakeOffCuts_(true),
     2969#if NEW_UPDATE_OBJECT>1
     2970  numberUpdateItems_(0),
     2971  maximumNumberUpdateItems_(0),
     2972  updateItems_(NULL),
     2973#endif
     2974  numberThreads_(0),
     2975  threadMode_(0)
    26442976{
    26452977  memset(intParam_,0,sizeof(intParam_));
     
    27803112      numberIntegers_++;
    27813113  }
     3114  delete [] integerVariable_;
    27823115  if (numberIntegers_) {
    2783     delete [] integerVariable_;
    27843116    integerVariable_ = new int [numberIntegers_];
    27853117    numberIntegers_=0;
     
    28193151  subTreeModel_(rhs.subTreeModel_),
    28203152  numberStoppedSubTrees_(rhs.numberStoppedSubTrees_),
     3153  mutex_(NULL),
    28213154  presolve_(rhs.presolve_),
    28223155  numberStrong_(rhs.numberStrong_),
     
    28553188  searchStrategy_(rhs.searchStrategy_),
    28563189  numberStrongIterations_(rhs.numberStrongIterations_),
    2857   resolveAfterTakeOffCuts_(rhs.resolveAfterTakeOffCuts_)
     3190  resolveAfterTakeOffCuts_(rhs.resolveAfterTakeOffCuts_),
     3191#if NEW_UPDATE_OBJECT>1
     3192  numberUpdateItems_(rhs.numberUpdateItems_),
     3193  maximumNumberUpdateItems_(rhs.maximumNumberUpdateItems_),
     3194  updateItems_(NULL),
     3195#endif
     3196  numberThreads_(rhs.numberThreads_),
     3197  threadMode_(rhs.threadMode_)
    28583198{
    28593199  memcpy(intParam_,rhs.intParam_,sizeof(intParam_));
     
    29073247      object_ = new OsiObject * [numberObjects_];
    29083248      int i;
    2909       for (i=0;i<numberObjects_;i++)
     3249      for (i=0;i<numberObjects_;i++) {
    29103250        object_[i]=(rhs.object_[i])->clone();
     3251        CbcObject * obj = dynamic_cast <CbcObject *>(object_[i]) ;
     3252        // Could be OsiObjects
     3253        if (obj)
     3254          obj->setModel(this);
     3255      }
    29113256    } else {
    29123257      object_=NULL;
     
    29323277    originalColumns_=NULL;
    29333278  }
     3279#if NEW_UPDATE_OBJECT>1
     3280  if (maximumNumberUpdateItems_) {
     3281    updateItems_ = new CbcObjectUpdateData [maximumNumberUpdateItems_];
     3282    for (int i=0;i<maximumNumberUpdateItems_;i++)
     3283      updateItems_[i] = rhs.updateItems_[i];
     3284  }
     3285#endif
    29343286  if (maximumWhich_&&rhs.whichGenerator_)
    29353287    whichGenerator_ = CoinCopyOfArray(rhs.whichGenerator_,maximumWhich_);
     
    31053457    subTreeModel_ = rhs.subTreeModel_;
    31063458    numberStoppedSubTrees_ = rhs.numberStoppedSubTrees_;
     3459    mutex_ = NULL;
    31073460    presolve_ = rhs.presolve_;
    31083461    numberStrong_ = rhs.numberStrong_;
     
    31563509    numberNewCuts_ = rhs.numberNewCuts_;
    31573510    resolveAfterTakeOffCuts_=rhs.resolveAfterTakeOffCuts_;
     3511#if NEW_UPDATE_OBJECT>1
     3512    numberUpdateItems_ = rhs.numberUpdateItems_;
     3513    maximumNumberUpdateItems_ = rhs.maximumNumberUpdateItems_;
     3514    delete [] updateItems_;
     3515    if (maximumNumberUpdateItems_) {
     3516      updateItems_ = new CbcObjectUpdateData [maximumNumberUpdateItems_];
     3517      for (i=0;i<maximumNumberUpdateItems_;i++)
     3518        updateItems_[i] = rhs.updateItems_[i];
     3519    } else {
     3520      updateItems_ = NULL;
     3521    }
     3522#endif
     3523    numberThreads_ = rhs.numberThreads_;
     3524    threadMode_ = rhs.threadMode_;
    31583525    sizeMiniTree_ = rhs.sizeMiniTree_;
    31593526    searchStrategy_ = rhs.searchStrategy_;
     
    33293696  originalColumns_=NULL;
    33303697  delete strategy_;
     3698#if NEW_UPDATE_OBJECT>1
     3699  delete [] updateItems_;
     3700  updateItems_=NULL;
     3701  numberUpdateItems_=0;
     3702  maximumNumberUpdateItems_=0;
     3703#endif
    33313704  gutsOfDestructor2();
    33323705}
     
    36864059*/
    36874060  currentNumberCuts_=currentNumberCuts;
    3688   if (currentNumberCuts >= maximumNumberCuts_) {
     4061  if (currentNumberCuts > maximumNumberCuts_) {
    36894062    maximumNumberCuts_ = currentNumberCuts;
    36904063    delete [] addedCuts_;
     
    37844157  reconstructed basis in the solver.
    37854158*/
    3786   if (node->objectiveValue() < cutoff)
     4159  if (node->objectiveValue() < cutoff||numberThreads_)
    37874160  {
    37884161#   ifdef CBC_CHECK_BASIS
     
    38644237  set. Adjust the cut reference counts to reflect that we no longer need to
    38654238  explore the remaining branch arms, hence they will no longer reference any
    3866   cuts. Cuts whose reference count falls to zero are deleted.
     4239  cuts. Cuts whose reference count falls to zero are deleted. 
    38674240*/
    38684241  else
    38694242  { int i;
    3870     int numberLeft = nodeInfo->numberBranchesLeft();
    3871     for (i = 0 ; i < currentNumberCuts ; i++)
    3872     { if (addedCuts_[i])
    3873       { if (!addedCuts_[i]->decrement(numberLeft))
    3874         { delete addedCuts_[i];
    3875           addedCuts_[i] = NULL; } } }
     4243    if (currentNumberCuts) {
     4244      lockThread();
     4245      int numberLeft = nodeInfo->numberBranchesLeft();
     4246      for (i = 0 ; i < currentNumberCuts ; i++)
     4247        { if (addedCuts_[i])
     4248          { if (!addedCuts_[i]->decrement(numberLeft))
     4249            { delete addedCuts_[i];
     4250            addedCuts_[i] = NULL; } } }
     4251      unlockThread();
     4252    }
    38764253    return 1 ; }
    38774254}
     
    39884365    saveClpOptions = clpSolver->specialOptions();
    39894366# endif
     4367#ifdef CBC_THREAD
     4368  CbcModel ** threadModel = NULL;
     4369  pthread_t * threadId = NULL;
     4370  pthread_cond_t condition_main;
     4371  pthread_mutex_t condition_mutex;
     4372  pthread_mutex_t * mutex2 = NULL;
     4373  pthread_cond_t * condition2 = NULL;
     4374  threadStruct * threadInfo = NULL;
     4375  void * saveMutex = NULL;
     4376  if (numberThreads_&&(threadMode_&2)!=0&&!numberNodes_) {
     4377    threadId = new pthread_t [numberThreads_];
     4378    pthread_cond_init(&condition_main,NULL);
     4379    pthread_mutex_init(&condition_mutex,NULL);
     4380    threadModel = new CbcModel * [numberThreads_];
     4381    threadInfo = new threadStruct [numberThreads_+1];
     4382    mutex2 = new pthread_mutex_t [numberThreads_];
     4383    condition2 = new pthread_cond_t [numberThreads_];
     4384    saveMutex = mutex_;
     4385    for (int i=0;i<numberThreads_;i++) {
     4386      pthread_mutex_init(mutex2+i,NULL);
     4387      pthread_cond_init(condition2+i,NULL);
     4388      threadId[i]=0;
     4389      threadModel[i]=new CbcModel;
     4390      threadModel[i]->generator_ = new CbcCutGenerator * [1];
     4391      delete threadModel[i]->solver_;
     4392      threadModel[i]->solver_=NULL;
     4393      threadModel[i]->numberThreads_=numberThreads_;
     4394      mutex_ = (void *) (threadInfo+i);
     4395      threadInfo[i].thisModel=(CbcModel *) threadModel[i];
     4396      threadInfo[i].baseModel=this;
     4397      threadInfo[i].threadIdOfBase=pthread_self();
     4398      threadInfo[i].mutex2=mutex2+i;
     4399      threadInfo[i].condition2=condition2+i;
     4400      threadInfo[i].returnCode=-1;
     4401      pthread_create(threadId+i,NULL,doCutsThread,threadInfo+i);
     4402    }
     4403    // Do a partial one for base model
     4404    threadInfo[numberThreads_].baseModel=this;
     4405    mutex_ = (void *) (threadInfo+numberThreads_);
     4406    threadInfo[numberThreads_].condition2=&condition_main;
     4407    threadInfo[numberThreads_].mutex2=&condition_mutex;
     4408  }
     4409#endif
    39904410  bool feasible = true ;
    39914411  int lastNumberCuts = 0 ;
     
    40264446  //printf("zero iterations on first solve of branch\n");
    40274447#endif
    4028   if (node&&!node->nodeInfo()->numberBranchesLeft())
     4448  if (node&&node->nodeInfo()&&!node->nodeInfo()->numberBranchesLeft())
    40294449    node->nodeInfo()->allBranchesGone(); // can clean up
    40304450  feasible = returnCode  != 0 ;
     
    40354455  }
    40364456 
     4457#if NEW_UPDATE_OBJECT==0
    40374458  // Update branching information if wanted
    40384459  if(node &&branchingMethod_)
    40394460    branchingMethod_->updateInformation(solver_,node);
     4461#elif NEW_UPDATE_OBJECT<2
     4462  // Update branching information if wanted
     4463  if(node &&branchingMethod_) {
     4464    OsiBranchingObject * bobj = node->modifiableBranchingObject();
     4465    CbcBranchingObject * cbcobj = dynamic_cast<CbcBranchingObject *> (bobj);
     4466    if (cbcobj) {
     4467      CbcObject * object = cbcobj->object();
     4468      CbcObjectUpdateData update = object->createUpdateInformation(solver_,node,cbcobj);
     4469      object->updateInformation(update);
     4470    } else {
     4471      branchingMethod_->updateInformation(solver_,node);
     4472    }
     4473  }
     4474#else
     4475  // Update branching information if wanted
     4476  if(node &&branchingMethod_) {
     4477    OsiBranchingObject * bobj = node->modifiableBranchingObject();
     4478    CbcBranchingObject * cbcobj = dynamic_cast<CbcBranchingObject *> (bobj);
     4479    if (cbcobj) {
     4480      CbcObject * object = cbcobj->object();
     4481      CbcObjectUpdateData update = object->createUpdateInformation(solver_,node,cbcobj);
     4482      // have to compute object number as not saved
     4483      CbcSimpleInteger * simpleObject =
     4484          dynamic_cast <CbcSimpleInteger *>(object) ;
     4485      int iObject;
     4486      int iColumn = simpleObject->columnNumber();
     4487      for (iObject = 0 ; iObject < numberObjects_ ; iObject++) {
     4488        simpleObject =
     4489          dynamic_cast <CbcSimpleInteger *>(object_[iObject]) ;
     4490        if (simpleObject->columnNumber()==iColumn)
     4491          break;
     4492      }
     4493      assert (iObject<numberObjects_);
     4494      update.objectNumber_ = iObject;
     4495      addUpdateInformation(update);
     4496    } else {
     4497      branchingMethod_->updateInformation(solver_,node);
     4498    }
     4499  }
     4500#endif
    40404501
    40414502#ifdef CBC_DEBUG
     
    40784539  sumChangeObjective1_ += solver_->getObjValue()*solver_->getObjSense()
    40794540    - objectiveValue ;
    4080   if ( CoinCpuTime()-dblParam_[CbcStartSeconds] > dblParam_[CbcMaximumSeconds] )
     4541  if ( getCurrentSeconds() > dblParam_[CbcMaximumSeconds] )
    40814542    numberTries=0; // exit
    40824543  //if ((numberNodes_%100)==0)
     
    42444705      whichGenerator_[numberViolated++]=-2;
    42454706    }
    4246     double * newSolution = new double [numberColumns] ;
    4247     double heuristicValue = getCutoff() ;
    4248     int found = -1; // no solution found
    42494707
    42504708    // reset probing info
    42514709    if (probingInfo_)
    42524710      probingInfo_->initializeFixing();
    4253     for (int i = 0;i<numberCutGenerators_+numberHeuristics_;i++) {
    4254       int numberRowCutsBefore = theseCuts.sizeRowCuts() ;
    4255       int numberColumnCutsBefore = theseCuts.sizeColCuts() ;
    4256       if (i<numberCutGenerators_) {
    4257         bool generate = generator_[i]->normal();
    4258         // skip if not optimal and should be (maybe a cut generator has fixed variables)
    4259         if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable())
    4260           generate=false;
    4261         if (generator_[i]->switchedOff())
    4262           generate=false;;
     4711    int i;
     4712    if ((threadMode_&2)==0||numberNodes_) {
     4713      for (i = 0;i<numberCutGenerators_;i++) {
     4714        int numberRowCutsBefore = theseCuts.sizeRowCuts() ;
     4715        int numberColumnCutsBefore = theseCuts.sizeColCuts() ;
     4716        bool generate = generator_[i]->normal();
     4717        // skip if not optimal and should be (maybe a cut generator has fixed variables)
     4718        if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable())
     4719          generate=false;
     4720        if (generator_[i]->switchedOff())
     4721          generate=false;;
    42634722        if (generate) {
    42644723          bool mustResolve =
    4265             generator_[i]->generateCuts(theseCuts,fullScan,node) ;
     4724            generator_[i]->generateCuts(theseCuts,fullScan,solver_,node) ;
    42664725          int numberRowCutsAfter = theseCuts.sizeRowCuts() ;
    4267           if(numberRowCutsBefore < numberRowCutsAfter &&
    4268              generator_[i]->mustCallAgain())
    4269             keepGoing=true; // say must go round
     4726          if(numberRowCutsBefore < numberRowCutsAfter &&
     4727             generator_[i]->mustCallAgain())
     4728            keepGoing=true; // say must go round
    42704729          // Check last cut to see if infeasible
    4271           if(numberRowCutsBefore < numberRowCutsAfter) {
     4730          if(numberRowCutsBefore < numberRowCutsAfter) {
    42724731            const OsiRowCut * thisCut = theseCuts.rowCutPtr(numberRowCutsAfter-1) ;
    42734732            if (thisCut->lb()>thisCut->ub()) {
     
    42774736          }
    42784737#ifdef CBC_DEBUG
    4279           {
    4280             int numberRowCutsAfter = theseCuts.sizeRowCuts() ;
    4281             int k ;
    4282             for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
    4283               OsiRowCut thisCut = theseCuts.rowCut(k) ;
    4284               /* check size of elements.
    4285                  We can allow smaller but this helps debug generators as it
    4286                  is unsafe to have small elements */
    4287               int n=thisCut.row().getNumElements();
    4288               const int * column = thisCut.row().getIndices();
    4289               const double * element = thisCut.row().getElements();
    4290               //assert (n);
    4291               for (int i=0;i<n;i++) {
    4292                 double value = element[i];
    4293                 assert(fabs(value)>1.0e-12&&fabs(value)<1.0e20);
    4294               }
    4295             }
    4296           }
     4738          {
     4739            int k ;
     4740            for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
     4741              OsiRowCut thisCut = theseCuts.rowCut(k) ;
     4742              /* check size of elements.
     4743                 We can allow smaller but this helps debug generators as it
     4744                 is unsafe to have small elements */
     4745              int n=thisCut.row().getNumElements();
     4746              const int * column = thisCut.row().getIndices();
     4747              const double * element = thisCut.row().getElements();
     4748              //assert (n);
     4749              for (int i=0;i<n;i++) {
     4750                double value = element[i];
     4751                assert(fabs(value)>1.0e-12&&fabs(value)<1.0e20);
     4752              }
     4753            }
     4754          }
    42974755#endif
    42984756          if (mustResolve) {
    4299             int returncode = resolve(node ? node->nodeInfo() : NULL,2);
    4300             feasible = returnCode  != 0 ;
    4301             if (returncode<0)
    4302               numberTries=0;
    4303             if ((specialOptions_&1)!=0) {
    4304               debugger = solver_->getRowCutDebugger() ;
    4305               if (debugger)
    4306                 onOptimalPath = (debugger->onOptimalPath(*solver_)) ;
    4307               else
    4308                 onOptimalPath=false;
    4309               if (onOptimalPath && !solver_->isDualObjectiveLimitReached())
    4310                 assert(feasible) ;
    4311             }
     4757            int returncode = resolve(node ? node->nodeInfo() : NULL,2);
     4758            feasible = returnCode  != 0 ;
     4759            if (returncode<0)
     4760              numberTries=0;
     4761            if ((specialOptions_&1)!=0) {
     4762              debugger = solver_->getRowCutDebugger() ;
     4763              if (debugger)
     4764                onOptimalPath = (debugger->onOptimalPath(*solver_)) ;
     4765              else
     4766                onOptimalPath=false;
     4767              if (onOptimalPath && !solver_->isDualObjectiveLimitReached())
     4768                assert(feasible) ;
     4769            }
    43124770            if (!feasible)
    43134771              break ;
    43144772          }
    43154773        }
     4774        int numberRowCutsAfter = theseCuts.sizeRowCuts() ;
     4775        int numberColumnCutsAfter = theseCuts.sizeColCuts() ;
     4776       
     4777        if ((specialOptions_&1)!=0) {
     4778          if (onOptimalPath) {
     4779            int k ;
     4780            for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
     4781              OsiRowCut thisCut = theseCuts.rowCut(k) ;
     4782              if(debugger->invalidCut(thisCut)) {
     4783                solver_->writeMps("badCut");
     4784#ifdef NDEBUG
     4785                printf("Cut generator %d (%s) produced invalid cut (%dth in this go)\n",
     4786                       i,generator_[i]->cutGeneratorName(),k-numberRowCutsBefore);
     4787                abort();
     4788#endif
     4789              }
     4790              assert(!debugger->invalidCut(thisCut)) ;
     4791            }
     4792          }
     4793        }
     4794/*
     4795  The cut generator has done its thing, and maybe it generated some
     4796  cuts.  Do a bit of bookkeeping: load
     4797  whichGenerator[i] with the index of the generator responsible for a cut,
     4798  and place cuts flagged as global in the global cut pool for the model.
     4799
     4800  lastNumberCuts is the sum of cuts added in previous iterations; it's the
     4801  offset to the proper starting position in whichGenerator.
     4802*/
     4803        int numberBefore =
     4804          numberRowCutsBefore+numberColumnCutsBefore+lastNumberCuts ;
     4805        int numberAfter =
     4806          numberRowCutsAfter+numberColumnCutsAfter+lastNumberCuts ;
     4807        // possibly extend whichGenerator
     4808        resizeWhichGenerator(numberBefore, numberAfter);
     4809        int j ;
     4810        if (fullScan) {
     4811          // counts
     4812          countColumnCuts[i] += numberColumnCutsAfter-numberColumnCutsBefore ;
     4813        }
     4814        countRowCuts[i] += numberRowCutsAfter-numberRowCutsBefore ;
     4815       
     4816        bool dodgyCuts=false;
     4817        for (j = numberRowCutsBefore;j<numberRowCutsAfter;j++) {
     4818          const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ;
     4819          if (thisCut->lb()>1.0e10||thisCut->ub()<-1.0e10) {
     4820            dodgyCuts=true;
     4821            break;
     4822          }
     4823          whichGenerator_[numberBefore++] = i ;
     4824          if (thisCut->lb()>thisCut->ub())
     4825            violated=-2; // sub-problem is infeasible
     4826          if (thisCut->globallyValid()) {
     4827            // add to global list
     4828            OsiRowCut newCut(*thisCut);
     4829            newCut.setGloballyValid(2);
     4830            newCut.mutableRow().setTestForDuplicateIndex(false);
     4831            globalCuts_.insert(newCut) ;
     4832          }
     4833        }
     4834        if (dodgyCuts) {
     4835          for (int k=numberRowCutsAfter-1;k>=j;k--) {
     4836            const OsiRowCut * thisCut = theseCuts.rowCutPtr(k) ;
     4837            if (thisCut->lb()>thisCut->ub())
     4838              violated=-2; // sub-problem is infeasible
     4839            if (thisCut->lb()>1.0e10||thisCut->ub()<-1.0e10)
     4840              theseCuts.eraseRowCut(k);
     4841          }
     4842          numberRowCutsAfter = theseCuts.sizeRowCuts() ;
     4843          for (;j<numberRowCutsAfter;j++) {
     4844            const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ;
     4845            whichGenerator_[numberBefore++] = i ;
     4846            if (thisCut->globallyValid()) {
     4847              // add to global list
     4848              OsiRowCut newCut(*thisCut);
     4849              newCut.setGloballyValid(2);
     4850              newCut.mutableRow().setTestForDuplicateIndex(false);
     4851              globalCuts_.insert(newCut) ;
     4852            }
     4853          }
     4854        }
     4855        for (j = numberColumnCutsBefore;j<numberColumnCutsAfter;j++) {
     4856          whichGenerator_[numberBefore++] = i ;
     4857          const OsiColCut * thisCut = theseCuts.colCutPtr(j) ;
     4858          if (thisCut->globallyValid()) {
     4859            // add to global list
     4860            OsiColCut newCut(*thisCut);
     4861            newCut.setGloballyValid(2);
     4862            globalCuts_.insert(newCut) ;
     4863          }
     4864        }
     4865      }
     4866      // Add in any violated saved cuts
     4867      if (!theseCuts.sizeRowCuts()&&!theseCuts.sizeColCuts()) {
     4868        int numberOld = theseCuts.sizeRowCuts();
     4869        int numberCuts = slackCuts.sizeRowCuts() ;
     4870        int i;
     4871        // possibly extend whichGenerator
     4872        resizeWhichGenerator(numberOld, numberOld+numberCuts);
     4873        for ( i = 0;i<numberCuts;i++) {
     4874          const OsiRowCut * thisCut = slackCuts.rowCutPtr(i) ;
     4875          if (thisCut->violated(cbcColSolution_)>100.0*primalTolerance) {
     4876            if (messageHandler()->logLevel()>2)
     4877              printf("Old cut added - violation %g\n",
     4878                     thisCut->violated(cbcColSolution_)) ;
     4879            whichGenerator_[numberOld++]=-1;
     4880            theseCuts.insert(*thisCut) ;
     4881          }
     4882        }
     4883      }
     4884    } else {
     4885      // do cuts independently
     4886      OsiCuts * eachCuts = new OsiCuts [numberCutGenerators_];;
     4887#ifdef CBC_THREAD
     4888      if (!threadModel) {
     4889#endif
     4890        // generate cuts
     4891        for (i = 0;i<numberCutGenerators_;i++) {
     4892          bool generate = generator_[i]->normal();
     4893          // skip if not optimal and should be (maybe a cut generator has fixed variables)
     4894          if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable())
     4895            generate=false;
     4896          if (generator_[i]->switchedOff())
     4897            generate=false;;
     4898          if (generate)
     4899            generator_[i]->generateCuts(eachCuts[i],fullScan,solver_,node) ;
     4900        }
     4901#ifdef CBC_THREAD
    43164902      } else {
     4903        for (i=0;i<numberThreads_;i++) {
     4904          // set solver here after cloning
     4905          threadModel[i]->solver_=solver_->clone();
     4906          threadModel[i]->numberNodes_ = (fullScan) ? 1 : 0;
     4907        }
     4908        // generate cuts
     4909        for (i = 0;i<numberCutGenerators_;i++) {
     4910          bool generate = generator_[i]->normal();
     4911          // skip if not optimal and should be (maybe a cut generator has fixed variables)
     4912          if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable())
     4913            generate=false;
     4914          if (generator_[i]->switchedOff())
     4915            generate=false;;
     4916          if (generate) {
     4917            bool finished=false;
     4918            int iThread=-1;
     4919            // see if any available
     4920            for (iThread=0;iThread<numberThreads_;iThread++) {
     4921              if (threadInfo[iThread].returnCode) {
     4922                finished=true;
     4923                break;
     4924              } else if (threadInfo[iThread].returnCode==0) {
     4925                pthread_cond_signal(threadInfo[iThread].condition2); // unlock
     4926              }
     4927            }
     4928            while (!finished) {
     4929              pthread_mutex_lock(&condition_mutex);
     4930              struct timespec absTime;
     4931              clock_gettime(CLOCK_REALTIME,&absTime);
     4932              absTime.tv_nsec += 1000000; // millisecond
     4933              if (absTime.tv_nsec>=1000000000) {
     4934                absTime.tv_nsec -= 1000000000;
     4935                absTime.tv_sec++;
     4936              }
     4937              pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime);
     4938              pthread_mutex_unlock(&condition_mutex);
     4939              for (iThread=0;iThread<numberThreads_;iThread++) {
     4940                if (threadInfo[iThread].returnCode>0) {
     4941                  finished=true;
     4942                  break;
     4943                } else if (threadInfo[iThread].returnCode==0) {
     4944                  pthread_cond_signal(threadInfo[iThread].condition2); // unlock
     4945                }
     4946              }
     4947            }
     4948            assert (iThread<numberThreads_);
     4949            assert (threadInfo[iThread].returnCode);
     4950            threadModel[iThread]->generator_[0]=generator_[i];
     4951            threadModel[iThread]->object_ = (OsiObject **) (eachCuts+i);
     4952            // allow to start
     4953            threadInfo[iThread].returnCode=0;
     4954            pthread_cond_signal(threadInfo[iThread].condition2); // unlock
     4955          }
     4956        }
     4957        // wait
     4958        for (int iThread=0;iThread<numberThreads_;iThread++) {
     4959          if (threadInfo[iThread].returnCode==0) {
     4960            bool finished=false;
     4961            pthread_cond_signal(threadInfo[iThread].condition2); // unlock
     4962            while (!finished) {
     4963              pthread_mutex_lock(&condition_mutex);
     4964              struct timespec absTime;
     4965              clock_gettime(CLOCK_REALTIME,&absTime);
     4966              absTime.tv_nsec += 1000000; // millisecond
     4967              if (absTime.tv_nsec>=1000000000) {
     4968                absTime.tv_nsec -= 1000000000;
     4969                absTime.tv_sec++;
     4970              }
     4971              pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime);
     4972              pthread_mutex_unlock(&condition_mutex);
     4973              if (threadInfo[iThread].returnCode>0) {
     4974                finished=true;
     4975                break;
     4976              } else if (threadInfo[iThread].returnCode==0) {
     4977                pthread_cond_signal(threadInfo[iThread].condition2); // unlock
     4978              }
     4979            }
     4980          }
     4981          assert (threadInfo[iThread].returnCode);
     4982          // say available
     4983          threadInfo[iThread].returnCode=-1;
     4984          delete threadModel[iThread]->solver_;
     4985          threadModel[iThread]->solver_=NULL;
     4986        }
     4987      }
     4988#endif
     4989      // Now put together
     4990      for (i = 0;i<numberCutGenerators_;i++) {
     4991        // add column cuts
     4992        int numberColumnCutsBefore = theseCuts.sizeColCuts() ;
     4993        int numberColumnCuts = eachCuts[i].sizeColCuts();
     4994        int numberColumnCutsAfter = numberColumnCutsBefore
     4995          + numberColumnCuts;
     4996        int j;
     4997        for (j=0;j<numberColumnCuts;j++) {
     4998          theseCuts.insert(eachCuts[i].colCut(j));
     4999        }
     5000        int numberRowCutsBefore = theseCuts.sizeRowCuts() ;
     5001        int numberRowCuts = eachCuts[i].sizeRowCuts();
     5002        int numberRowCutsAfter = numberRowCutsBefore
     5003          + numberRowCuts;
     5004        if (numberRowCuts) {
     5005          for (j=0;j<numberRowCuts;j++) {
     5006            const OsiRowCut * thisCut = eachCuts[i].rowCutPtr(j) ;
     5007            if (thisCut->lb()<=1.0e10&&thisCut->ub()>=-1.0e10)
     5008              theseCuts.insert(eachCuts[i].rowCut(j));
     5009          }
     5010          if (generator_[i]->mustCallAgain())
     5011            keepGoing=true; // say must go round
     5012          // Check last cut to see if infeasible
     5013          const OsiRowCut * thisCut = theseCuts.rowCutPtr(numberRowCutsAfter-1) ;
     5014          if (thisCut->lb()>thisCut->ub()) {
     5015            feasible = false; // sub-problem is infeasible
     5016            break;
     5017          }
     5018        }
     5019#ifdef CBC_DEBUG
     5020        {
     5021          int k ;
     5022          for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
     5023            OsiRowCut thisCut = theseCuts.rowCut(k) ;
     5024            /* check size of elements.
     5025               We can allow smaller but this helps debug generators as it
     5026               is unsafe to have small elements */
     5027            int n=thisCut.row().getNumElements();
     5028            const int * column = thisCut.row().getIndices();
     5029            const double * element = thisCut.row().getElements();
     5030            //assert (n);
     5031            for (int i=0;i<n;i++) {
     5032              double value = element[i];
     5033              assert(fabs(value)>1.0e-12&&fabs(value)<1.0e20);
     5034            }
     5035          }
     5036        }
     5037#endif
     5038        if ((specialOptions_&1)!=0) {
     5039          if (onOptimalPath) {
     5040            int k ;
     5041            for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
     5042              OsiRowCut thisCut = theseCuts.rowCut(k) ;
     5043              if(debugger->invalidCut(thisCut)) {
     5044                solver_->writeMps("badCut");
     5045#ifdef NDEBUG
     5046                printf("Cut generator %d (%s) produced invalid cut (%dth in this go)\n",
     5047                       i,generator_[i]->cutGeneratorName(),k-numberRowCutsBefore);
     5048                abort();
     5049#endif
     5050              }
     5051              assert(!debugger->invalidCut(thisCut)) ;
     5052            }
     5053          }
     5054        }
     5055/*
     5056  The cut generator has done its thing, and maybe it generated some
     5057  cuts.  Do a bit of bookkeeping: load
     5058  whichGenerator[i] with the index of the generator responsible for a cut,
     5059  and place cuts flagged as global in the global cut pool for the model.
     5060
     5061  lastNumberCuts is the sum of cuts added in previous iterations; it's the
     5062  offset to the proper starting position in whichGenerator.
     5063*/
     5064        int numberBefore =
     5065          numberRowCutsBefore+numberColumnCutsBefore+lastNumberCuts ;
     5066        int numberAfter =
     5067          numberRowCutsAfter+numberColumnCutsAfter+lastNumberCuts ;
     5068        // possibly extend whichGenerator
     5069        resizeWhichGenerator(numberBefore, numberAfter);
     5070        if (fullScan) {
     5071          // counts
     5072          countColumnCuts[i] += numberColumnCutsAfter-numberColumnCutsBefore ;
     5073        }
     5074        countRowCuts[i] += numberRowCutsAfter-numberRowCutsBefore ;
     5075       
     5076        for (j = numberRowCutsBefore;j<numberRowCutsAfter;j++) {
     5077          whichGenerator_[numberBefore++] = i ;
     5078          const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ;
     5079          if (thisCut->lb()>thisCut->ub())
     5080            violated=-2; // sub-problem is infeasible
     5081          if (thisCut->globallyValid()) {
     5082            // add to global list
     5083            OsiRowCut newCut(*thisCut);
     5084            newCut.setGloballyValid(2);
     5085            newCut.mutableRow().setTestForDuplicateIndex(false);
     5086            globalCuts_.insert(newCut) ;
     5087          }
     5088        }
     5089        for (j = numberColumnCutsBefore;j<numberColumnCutsAfter;j++) {
     5090          whichGenerator_[numberBefore++] = i ;
     5091          const OsiColCut * thisCut = theseCuts.colCutPtr(j) ;
     5092          if (thisCut->globallyValid()) {
     5093            // add to global list
     5094            OsiColCut newCut(*thisCut);
     5095            newCut.setGloballyValid(2);
     5096            globalCuts_.insert(newCut) ;
     5097          }
     5098        }
     5099      }
     5100      // Add in any violated saved cuts
     5101      if (!theseCuts.sizeRowCuts()&&!theseCuts.sizeColCuts()) {
     5102        int numberOld = theseCuts.sizeRowCuts();
     5103        int numberCuts = slackCuts.sizeRowCuts() ;
     5104        int i;
     5105        // possibly extend whichGenerator
     5106        resizeWhichGenerator(numberOld, numberOld+numberCuts);
     5107        for ( i = 0;i<numberCuts;i++) {
     5108          const OsiRowCut * thisCut = slackCuts.rowCutPtr(i) ;
     5109          if (thisCut->violated(cbcColSolution_)>100.0*primalTolerance) {
     5110            if (messageHandler()->logLevel()>2)
     5111              printf("Old cut added - violation %g\n",
     5112                     thisCut->violated(cbcColSolution_)) ;
     5113            whichGenerator_[numberOld++]=-1;
     5114            theseCuts.insert(*thisCut) ;
     5115          }
     5116        }
     5117      }
     5118      delete [] eachCuts;
     5119    }
     5120    //if (!feasible)
     5121    //break;
     5122/*
     5123  End of the loop to exercise each generator - try heuristics
     5124  - unless at root node and first pass
     5125*/
     5126    if (numberNodes_||currentPassNumber_!=1) {
     5127      double * newSolution = new double [numberColumns] ;
     5128      double heuristicValue = getCutoff() ;
     5129      int found = -1; // no solution found
     5130      for (i = 0;i<numberHeuristics_;i++) {
    43175131        // see if heuristic will do anything
    43185132        double saveValue = heuristicValue ;
    43195133        int ifSol =
    4320           heuristic_[i-numberCutGenerators_]->solution(heuristicValue,
    4321                                                        newSolution,
    4322                                                        theseCuts) ;
     5134          heuristic_[i]->solution(heuristicValue,
     5135                                  newSolution,
     5136                                  theseCuts) ;
    43235137        if (ifSol>0) {
    43245138          // better solution found
    4325           found = i-numberCutGenerators_ ;
     5139          found = i ;
    43265140          incrementUsed(newSolution);
    43275141        } else if (ifSol<0) {
     
    43295143        }
    43305144      }
    4331       int numberRowCutsAfter = theseCuts.sizeRowCuts() ;
    4332       int numberColumnCutsAfter = theseCuts.sizeColCuts() ;
    4333 
    4334       if ((specialOptions_&1)!=0) {
    4335         if (onOptimalPath) {
    4336           int k ;
    4337           for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
    4338             OsiRowCut thisCut = theseCuts.rowCut(k) ;
    4339             if(debugger->invalidCut(thisCut)) {
    4340               solver_->writeMps("badCut");
    4341             }
    4342             assert(!debugger->invalidCut(thisCut)) ;
    4343           }
    4344         }
    4345       }
    4346 /*
    4347   The cut generator/heuristic has done its thing, and maybe it generated some
    4348   cuts and/or a new solution.  Do a bit of bookkeeping: load
    4349   whichGenerator[i] with the index of the generator responsible for a cut,
    4350   and place cuts flagged as global in the global cut pool for the model.
    4351 
    4352   lastNumberCuts is the sum of cuts added in previous iterations; it's the
    4353   offset to the proper starting position in whichGenerator.
    4354 */
    4355       int numberBefore =
    4356             numberRowCutsBefore+numberColumnCutsBefore+lastNumberCuts ;
    4357       int numberAfter =
    4358             numberRowCutsAfter+numberColumnCutsAfter+lastNumberCuts ;
    4359       // possibly extend whichGenerator
    4360       resizeWhichGenerator(numberBefore, numberAfter);
    4361       int j ;
    4362       if (fullScan) {
    4363         // counts
    4364         countColumnCuts[i] += numberColumnCutsAfter-numberColumnCutsBefore ;
    4365       }
    4366       countRowCuts[i] += numberRowCutsAfter-numberRowCutsBefore ;
    4367      
    4368       for (j = numberRowCutsBefore;j<numberRowCutsAfter;j++) {
    4369         whichGenerator_[numberBefore++] = i ;
    4370         const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ;
    4371         if (thisCut->lb()>thisCut->ub())
    4372           violated=-2; // sub-problem is infeasible
    4373         if (thisCut->globallyValid()) {
    4374           // add to global list
    4375           OsiRowCut newCut(*thisCut);
    4376           newCut.setGloballyValid(2);
    4377           newCut.mutableRow().setTestForDuplicateIndex(false);
    4378           globalCuts_.insert(newCut) ;
    4379         }
    4380       }
    4381       for (j = numberColumnCutsBefore;j<numberColumnCutsAfter;j++) {
    4382         whichGenerator_[numberBefore++] = i ;
    4383         const OsiColCut * thisCut = theseCuts.colCutPtr(j) ;
    4384         if (thisCut->globallyValid()) {
    4385           // add to global list
    4386           OsiColCut newCut(*thisCut);
    4387           newCut.setGloballyValid(2);
    4388           globalCuts_.insert(newCut) ;
    4389         }
    4390       }
    4391     }
    4392     // If at root node and first pass do heuristics without cuts
    4393     if (!numberNodes_&&currentPassNumber_==1) {
    4394       // Save number solutions
    4395       int saveNumberSolutions = numberSolutions_;
    4396       int saveNumberHeuristicSolutions = numberHeuristicSolutions_;
    4397       // make sure bounds tight
    4398       {
    4399         for (int i = 0;i < numberObjects_;i++)
    4400           object_[i]->resetBounds(solver_) ;
    4401       }
    4402       for (int i = 0;i<numberHeuristics_;i++) {
    4403         // see if heuristic will do anything
    4404         double saveValue = heuristicValue ;
    4405         int ifSol = heuristic_[i]->solution(heuristicValue,
    4406                                             newSolution);
    4407         if (ifSol>0) {
    4408           // better solution found
    4409           found = i ;
    4410           incrementUsed(newSolution);
    4411           // increment number of solutions so other heuristics can test
    4412           numberSolutions_++;
    4413           numberHeuristicSolutions_++;
    4414         } else {
    4415           heuristicValue = saveValue ;
    4416         }
    4417       }
    4418       // Restore number solutions
    4419       numberSolutions_ = saveNumberSolutions;
    4420       numberHeuristicSolutions_ = saveNumberHeuristicSolutions;
    4421     }
    4422     // Add in any violated saved cuts
    4423     if (!theseCuts.sizeRowCuts()&&!theseCuts.sizeColCuts()) {
    4424       int numberOld = theseCuts.sizeRowCuts();
    4425       int numberCuts = slackCuts.sizeRowCuts() ;
    4426       int i;
    4427       // possibly extend whichGenerator
    4428       resizeWhichGenerator(numberOld, numberOld+numberCuts);
    4429       for ( i = 0;i<numberCuts;i++) {
    4430         const OsiRowCut * thisCut = slackCuts.rowCutPtr(i) ;
    4431         if (thisCut->violated(cbcColSolution_)>100.0*primalTolerance) {
    4432           if (messageHandler()->logLevel()>2)
    4433             printf("Old cut added - violation %g\n",
    4434                    thisCut->violated(cbcColSolution_)) ;
    4435           whichGenerator_[numberOld++]=-1;
    4436           theseCuts.insert(*thisCut) ;
    4437         }
    4438       }
    4439     }
    4440 /*
    4441   End of the loop to exercise each generator/heuristic.
    4442 
     5145/*
    44435146  Did any of the heuristics turn up a new solution? Record it before we free
    44445147  the vector.
    44455148*/
    4446     if (found >= 0) {
    4447       phase_=4;
    4448       incrementUsed(newSolution);
    4449       setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;
    4450       lastHeuristic_ = heuristic_[found];
    4451       CbcTreeLocal * tree
     5149      if (found >= 0) {
     5150        phase_=4;
     5151        incrementUsed(newSolution);
     5152        lastHeuristic_ = heuristic_[found];
     5153        setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;
     5154        CbcTreeLocal * tree
    44525155          = dynamic_cast<CbcTreeLocal *> (tree_);
    4453       if (tree)
    4454         tree->passInSolution(bestSolution_,heuristicValue);
    4455     }
    4456     delete [] newSolution ;
     5156        if (tree)
     5157          tree->passInSolution(bestSolution_,heuristicValue);
     5158      }
     5159      delete [] newSolution ;
     5160    }
    44575161
    44585162#if 0
     
    46165320      }
    46175321      feasible = resolve(node ? node->nodeInfo() : NULL,2) ;
    4618       if ( CoinCpuTime()-dblParam_[CbcStartSeconds] > dblParam_[CbcMaximumSeconds] )
     5322      if ( getCurrentSeconds() > dblParam_[CbcMaximumSeconds] )
    46195323        numberTries=0; // exit
    46205324#     ifdef CBC_DEBUG
     
    46895393    if (!feasible)
    46905394    { int i ;
    4691       for (i = 0;i<currentNumberCuts_;i++) {
    4692         // take off node
    4693         if (addedCuts_[i]) {
    4694           if (!addedCuts_[i]->decrement())
    4695             delete addedCuts_[i] ;
    4696           addedCuts_[i] = NULL ;
     5395      if (currentNumberCuts_) {
     5396        lockThread();
     5397        for (i = 0;i<currentNumberCuts_;i++) {
     5398          // take off node
     5399          if (addedCuts_[i]) {
     5400            if (!addedCuts_[i]->decrement())
     5401              delete addedCuts_[i] ;
     5402            addedCuts_[i] = NULL ;
     5403          }
    46975404        }
    46985405      }
     
    48305537      phase_=4;
    48315538      incrementUsed(newSolution);
     5539      lastHeuristic_ = heuristic_[found];
    48325540      setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;
    4833       lastHeuristic_ = heuristic_[found];
    48345541    }
    48355542    delete [] newSolution ;
     
    51655872    clpSolver->setSpecialOptions(saveClpOptions);
    51665873# endif
     5874#ifdef CBC_THREAD
     5875  if (threadModel) {
     5876    // stop threads
     5877    int i;
     5878    for (i=0;i<numberThreads_;i++) {
     5879      while (threadInfo[i].returnCode==0) {
     5880        pthread_cond_signal(threadInfo[i].condition2); // unlock
     5881        pthread_mutex_lock(&condition_mutex);
     5882        struct timespec absTime;
     5883        cl