Ignore:
Timestamp:
Jan 3, 2019 2:03:23 PM (2 years ago)
Author:
unxusr
Message:

.clang-format with proposal for formatting code

File:
1 edited

Legend:

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

    r2382 r2464  
    66#if defined(_MSC_VER)
    77// Turn off compiler warning about long names
    8 #  pragma warning(disable:4786)
     8#pragma warning(disable : 4786)
    99#endif
    1010
     
    3838//==============================================================================
    3939
    40 CbcHeuristicNode::CbcHeuristicNode(const CbcHeuristicNode& rhs)
    41 {
    42     numObjects_ = rhs.numObjects_;
    43     brObj_ = new CbcBranchingObject*[numObjects_];
    44     for (int i = 0; i < numObjects_; ++i) {
    45         brObj_[i] = rhs.brObj_[i]->clone();
    46     }
    47 }
    48 
    49 void
    50 CbcHeuristicNodeList::gutsOfDelete()
    51 {
    52     for (int i = (static_cast<int>(nodes_.size())) - 1; i >= 0; --i) {
    53         delete nodes_[i];
    54     }
    55 }
    56 
    57 void
    58 CbcHeuristicNodeList::gutsOfCopy(const CbcHeuristicNodeList& rhs)
    59 {
    60     append(rhs);
    61 }
    62 
    63 CbcHeuristicNodeList::CbcHeuristicNodeList(const CbcHeuristicNodeList& rhs)
    64 {
     40CbcHeuristicNode::CbcHeuristicNode(const CbcHeuristicNode &rhs)
     41{
     42  numObjects_ = rhs.numObjects_;
     43  brObj_ = new CbcBranchingObject *[numObjects_];
     44  for (int i = 0; i < numObjects_; ++i) {
     45    brObj_[i] = rhs.brObj_[i]->clone();
     46  }
     47}
     48
     49void CbcHeuristicNodeList::gutsOfDelete()
     50{
     51  for (int i = (static_cast<int>(nodes_.size())) - 1; i >= 0; --i) {
     52    delete nodes_[i];
     53  }
     54}
     55
     56void CbcHeuristicNodeList::gutsOfCopy(const CbcHeuristicNodeList &rhs)
     57{
     58  append(rhs);
     59}
     60
     61CbcHeuristicNodeList::CbcHeuristicNodeList(const CbcHeuristicNodeList &rhs)
     62{
     63  gutsOfCopy(rhs);
     64}
     65
     66CbcHeuristicNodeList &CbcHeuristicNodeList::operator=(const CbcHeuristicNodeList &rhs)
     67{
     68  if (this != &rhs) {
     69    gutsOfDelete();
    6570    gutsOfCopy(rhs);
    66 }
    67 
    68 CbcHeuristicNodeList& CbcHeuristicNodeList::operator=
    69 (const CbcHeuristicNodeList & rhs)
    70 {
    71     if (this != &rhs) {
    72         gutsOfDelete();
    73         gutsOfCopy(rhs);
    74     }
    75     return *this;
     71  }
     72  return *this;
    7673}
    7774
    7875CbcHeuristicNodeList::~CbcHeuristicNodeList()
    7976{
    80     gutsOfDelete();
    81 }
    82 
    83 void
    84 CbcHeuristicNodeList::append(CbcHeuristicNode*& node)
    85 {
    86     nodes_.push_back(node);
    87     node = NULL;
    88 }
    89 
    90 void
    91 CbcHeuristicNodeList::append(const CbcHeuristicNodeList& nodes)
    92 {
    93     nodes_.reserve(nodes_.size() + nodes.size());
    94     for (int i = 0; i < nodes.size(); ++i) {
    95         CbcHeuristicNode* node = new CbcHeuristicNode(*nodes.node(i));
    96         append(node);
    97     }
     77  gutsOfDelete();
     78}
     79
     80void CbcHeuristicNodeList::append(CbcHeuristicNode *&node)
     81{
     82  nodes_.push_back(node);
     83  node = NULL;
     84}
     85
     86void CbcHeuristicNodeList::append(const CbcHeuristicNodeList &nodes)
     87{
     88  nodes_.reserve(nodes_.size() + nodes.size());
     89  for (int i = 0; i < nodes.size(); ++i) {
     90    CbcHeuristicNode *node = new CbcHeuristicNode(*nodes.node(i));
     91    append(node);
     92  }
    9893}
    9994
    10095//==============================================================================
    101 #define DEFAULT_WHERE ((255-2-16)*(1+256))
     96#define DEFAULT_WHERE ((255 - 2 - 16) * (1 + 256))
    10297// Default Constructor
    103 CbcHeuristic::CbcHeuristic() :
    104         model_(NULL),
    105         when_(2),
    106         numberNodes_(200),
    107         feasibilityPumpOptions_(-1),
    108         fractionSmall_(1.0),
    109         heuristicName_("Unknown"),
    110         howOften_(1),
    111         decayFactor_(0.0),
    112         switches_(0),
    113         whereFrom_(DEFAULT_WHERE),
    114         shallowDepth_(1),
    115         howOftenShallow_(1),
    116         numInvocationsInShallow_(0),
    117         numInvocationsInDeep_(0),
    118         lastRunDeep_(0),
    119         numRuns_(0),
    120         minDistanceToRun_(1),
    121         runNodes_(),
    122         numCouldRun_(0),
    123         numberSolutionsFound_(0),
    124         numberNodesDone_(0),
    125         inputSolution_(NULL)
    126 {
    127     // As CbcHeuristic virtual need to modify .cpp if above change
     98CbcHeuristic::CbcHeuristic()
     99  : model_(NULL)
     100  , when_(2)
     101  , numberNodes_(200)
     102  , feasibilityPumpOptions_(-1)
     103  , fractionSmall_(1.0)
     104  , heuristicName_("Unknown")
     105  , howOften_(1)
     106  , decayFactor_(0.0)
     107  , switches_(0)
     108  , whereFrom_(DEFAULT_WHERE)
     109  , shallowDepth_(1)
     110  , howOftenShallow_(1)
     111  , numInvocationsInShallow_(0)
     112  , numInvocationsInDeep_(0)
     113  , lastRunDeep_(0)
     114  , numRuns_(0)
     115  , minDistanceToRun_(1)
     116  , runNodes_()
     117  , numCouldRun_(0)
     118  , numberSolutionsFound_(0)
     119  , numberNodesDone_(0)
     120  , inputSolution_(NULL)
     121{
     122  // As CbcHeuristic virtual need to modify .cpp if above change
    128123}
    129124
    130125// Constructor from model
    131 CbcHeuristic::CbcHeuristic(CbcModel & model) :
    132         model_(&model),
    133         when_(2),
    134         numberNodes_(200),
    135         feasibilityPumpOptions_(-1),
    136         fractionSmall_(1.0),
    137         heuristicName_("Unknown"),
    138         howOften_(1),
    139         decayFactor_(0.0),
    140         switches_(0),
    141         whereFrom_(DEFAULT_WHERE),
    142         shallowDepth_(1),
    143         howOftenShallow_(1),
    144         numInvocationsInShallow_(0),
    145         numInvocationsInDeep_(0),
    146         lastRunDeep_(0),
    147         numRuns_(0),
    148         minDistanceToRun_(1),
    149         runNodes_(),
    150         numCouldRun_(0),
    151         numberSolutionsFound_(0),
    152         numberNodesDone_(0),
    153         inputSolution_(NULL)
    154 {
    155 }
    156 
    157 void
    158 CbcHeuristic::gutsOfCopy(const CbcHeuristic & rhs)
    159 {
    160     model_ = rhs.model_;
    161     when_ = rhs.when_;
    162     numberNodes_ = rhs.numberNodes_;
    163     feasibilityPumpOptions_ = rhs.feasibilityPumpOptions_;
    164     fractionSmall_ = rhs.fractionSmall_;
    165     randomNumberGenerator_ = rhs.randomNumberGenerator_;
    166     heuristicName_ = rhs.heuristicName_;
    167     howOften_ = rhs.howOften_;
    168     decayFactor_ = rhs.decayFactor_;
    169     switches_ = rhs.switches_;
    170     whereFrom_ = rhs.whereFrom_;
    171     shallowDepth_ = rhs.shallowDepth_;
    172     howOftenShallow_ = rhs.howOftenShallow_;
    173     numInvocationsInShallow_ = rhs.numInvocationsInShallow_;
    174     numInvocationsInDeep_ = rhs.numInvocationsInDeep_;
    175     lastRunDeep_ = rhs.lastRunDeep_;
    176     numRuns_ = rhs.numRuns_;
    177     numCouldRun_ = rhs.numCouldRun_;
    178     minDistanceToRun_ = rhs.minDistanceToRun_;
    179     runNodes_ = rhs.runNodes_;
    180     numberSolutionsFound_ = rhs.numberSolutionsFound_;
    181     numberNodesDone_ = rhs.numberNodesDone_;
    182     if (rhs.inputSolution_) {
    183         int numberColumns = model_->getNumCols();
    184         setInputSolution(rhs.inputSolution_, rhs.inputSolution_[numberColumns]);
    185     }
     126CbcHeuristic::CbcHeuristic(CbcModel &model)
     127  : model_(&model)
     128  , when_(2)
     129  , numberNodes_(200)
     130  , feasibilityPumpOptions_(-1)
     131  , fractionSmall_(1.0)
     132  , heuristicName_("Unknown")
     133  , howOften_(1)
     134  , decayFactor_(0.0)
     135  , switches_(0)
     136  , whereFrom_(DEFAULT_WHERE)
     137  , shallowDepth_(1)
     138  , howOftenShallow_(1)
     139  , numInvocationsInShallow_(0)
     140  , numInvocationsInDeep_(0)
     141  , lastRunDeep_(0)
     142  , numRuns_(0)
     143  , minDistanceToRun_(1)
     144  , runNodes_()
     145  , numCouldRun_(0)
     146  , numberSolutionsFound_(0)
     147  , numberNodesDone_(0)
     148  , inputSolution_(NULL)
     149{
     150}
     151
     152void CbcHeuristic::gutsOfCopy(const CbcHeuristic &rhs)
     153{
     154  model_ = rhs.model_;
     155  when_ = rhs.when_;
     156  numberNodes_ = rhs.numberNodes_;
     157  feasibilityPumpOptions_ = rhs.feasibilityPumpOptions_;
     158  fractionSmall_ = rhs.fractionSmall_;
     159  randomNumberGenerator_ = rhs.randomNumberGenerator_;
     160  heuristicName_ = rhs.heuristicName_;
     161  howOften_ = rhs.howOften_;
     162  decayFactor_ = rhs.decayFactor_;
     163  switches_ = rhs.switches_;
     164  whereFrom_ = rhs.whereFrom_;
     165  shallowDepth_ = rhs.shallowDepth_;
     166  howOftenShallow_ = rhs.howOftenShallow_;
     167  numInvocationsInShallow_ = rhs.numInvocationsInShallow_;
     168  numInvocationsInDeep_ = rhs.numInvocationsInDeep_;
     169  lastRunDeep_ = rhs.lastRunDeep_;
     170  numRuns_ = rhs.numRuns_;
     171  numCouldRun_ = rhs.numCouldRun_;
     172  minDistanceToRun_ = rhs.minDistanceToRun_;
     173  runNodes_ = rhs.runNodes_;
     174  numberSolutionsFound_ = rhs.numberSolutionsFound_;
     175  numberNodesDone_ = rhs.numberNodesDone_;
     176  if (rhs.inputSolution_) {
     177    int numberColumns = model_->getNumCols();
     178    setInputSolution(rhs.inputSolution_, rhs.inputSolution_[numberColumns]);
     179  }
    186180}
    187181// Copy constructor
    188 CbcHeuristic::CbcHeuristic(const CbcHeuristic & rhs)
    189 {
    190     inputSolution_ = NULL;
    191     gutsOfCopy(rhs);
     182CbcHeuristic::CbcHeuristic(const CbcHeuristic &rhs)
     183{
     184  inputSolution_ = NULL;
     185  gutsOfCopy(rhs);
    192186}
    193187
    194188// Assignment operator
    195189CbcHeuristic &
    196 CbcHeuristic::operator=( const CbcHeuristic & rhs)
    197 {
    198     if (this != &rhs) {
    199         gutsOfDelete();
    200         gutsOfCopy(rhs);
    201     }
    202     return *this;
    203 }
    204 
    205 void CbcHeurDebugNodes(CbcModel* model_)
    206 {
    207     CbcNode* node = model_->currentNode();
    208     CbcNodeInfo* nodeInfo = node->nodeInfo();
    209     std::cout << "===============================================================\n";
    210     while (nodeInfo) {
    211         const CbcNode* node = nodeInfo->owner();
    212         printf("nodeinfo: node %i\n", nodeInfo->nodeNumber());
    213         {
    214             const CbcIntegerBranchingObject* brPrint =
    215                 dynamic_cast<const CbcIntegerBranchingObject*>(nodeInfo->parentBranch());
    216             if (!brPrint) {
    217                 printf("    parentBranch: NULL\n");
    218             } else {
    219                 const double* downBounds = brPrint->downBounds();
    220                 const double* upBounds = brPrint->upBounds();
    221                 int variable = brPrint->variable();
    222                 int way = brPrint->way();
    223                 printf("   parentBranch: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
    224                        variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
    225                        static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
    226             }
    227         }
    228         if (! node) {
    229             printf("    owner: NULL\n");
    230         } else {
    231             printf("    owner: node %i depth %i onTree %i active %i",
    232                    node->nodeNumber(), node->depth(), node->onTree(), node->active());
    233             const OsiBranchingObject* osibr =
    234                 nodeInfo->owner()->branchingObject();
    235             const CbcBranchingObject* cbcbr =
    236                 dynamic_cast<const CbcBranchingObject*>(osibr);
    237             const CbcIntegerBranchingObject* brPrint =
    238                 dynamic_cast<const CbcIntegerBranchingObject*>(cbcbr);
    239             if (!brPrint) {
    240                 printf("        ownerBranch: NULL\n");
    241             } else {
    242                 const double* downBounds = brPrint->downBounds();
    243                 const double* upBounds = brPrint->upBounds();
    244                 int variable = brPrint->variable();
    245                 int way = brPrint->way();
    246                 printf("        ownerbranch: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
    247                        variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
    248                        static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
    249             }
    250         }
    251         nodeInfo = nodeInfo->parent();
    252     }
    253 }
    254 
    255 void
    256 CbcHeuristic::debugNodes()
    257 {
    258     CbcHeurDebugNodes(model_);
    259 }
    260 
    261 void
    262 CbcHeuristic::printDistanceToNodes()
    263 {
    264     const CbcNode* currentNode = model_->currentNode();
    265     if (currentNode != NULL) {
    266         CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_);
    267         for (int i = runNodes_.size() - 1; i >= 0; --i) {
    268             nodeDesc->distance(runNodes_.node(i));
    269         }
    270         runNodes_.append(nodeDesc);
    271     }
    272 }
    273 
    274 bool
    275 CbcHeuristic::shouldHeurRun(int whereFrom)
    276 {
    277     assert (whereFrom >= 0 && whereFrom < 16);
    278     // take off 8 (code - likes new solution)
    279     whereFrom &= 7;
    280     if ((whereFrom_&(1 << whereFrom)) == 0)
    281         return false;
     190CbcHeuristic::operator=(const CbcHeuristic &rhs)
     191{
     192  if (this != &rhs) {
     193    gutsOfDelete();
     194    gutsOfCopy(rhs);
     195  }
     196  return *this;
     197}
     198
     199void CbcHeurDebugNodes(CbcModel *model_)
     200{
     201  CbcNode *node = model_->currentNode();
     202  CbcNodeInfo *nodeInfo = node->nodeInfo();
     203  std::cout << "===============================================================\n";
     204  while (nodeInfo) {
     205    const CbcNode *node = nodeInfo->owner();
     206    printf("nodeinfo: node %i\n", nodeInfo->nodeNumber());
     207    {
     208      const CbcIntegerBranchingObject *brPrint = dynamic_cast<const CbcIntegerBranchingObject *>(nodeInfo->parentBranch());
     209      if (!brPrint) {
     210        printf("    parentBranch: NULL\n");
     211      } else {
     212        const double *downBounds = brPrint->downBounds();
     213        const double *upBounds = brPrint->upBounds();
     214        int variable = brPrint->variable();
     215        int way = brPrint->way();
     216        printf("   parentBranch: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
     217          variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
     218          static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
     219      }
     220    }
     221    if (!node) {
     222      printf("    owner: NULL\n");
     223    } else {
     224      printf("    owner: node %i depth %i onTree %i active %i",
     225        node->nodeNumber(), node->depth(), node->onTree(), node->active());
     226      const OsiBranchingObject *osibr = nodeInfo->owner()->branchingObject();
     227      const CbcBranchingObject *cbcbr = dynamic_cast<const CbcBranchingObject *>(osibr);
     228      const CbcIntegerBranchingObject *brPrint = dynamic_cast<const CbcIntegerBranchingObject *>(cbcbr);
     229      if (!brPrint) {
     230        printf("        ownerBranch: NULL\n");
     231      } else {
     232        const double *downBounds = brPrint->downBounds();
     233        const double *upBounds = brPrint->upBounds();
     234        int variable = brPrint->variable();
     235        int way = brPrint->way();
     236        printf("        ownerbranch: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
     237          variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
     238          static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
     239      }
     240    }
     241    nodeInfo = nodeInfo->parent();
     242  }
     243}
     244
     245void CbcHeuristic::debugNodes()
     246{
     247  CbcHeurDebugNodes(model_);
     248}
     249
     250void CbcHeuristic::printDistanceToNodes()
     251{
     252  const CbcNode *currentNode = model_->currentNode();
     253  if (currentNode != NULL) {
     254    CbcHeuristicNode *nodeDesc = new CbcHeuristicNode(*model_);
     255    for (int i = runNodes_.size() - 1; i >= 0; --i) {
     256      nodeDesc->distance(runNodes_.node(i));
     257    }
     258    runNodes_.append(nodeDesc);
     259  }
     260}
     261
     262bool CbcHeuristic::shouldHeurRun(int whereFrom)
     263{
     264  assert(whereFrom >= 0 && whereFrom < 16);
     265  // take off 8 (code - likes new solution)
     266  whereFrom &= 7;
     267  if ((whereFrom_ & (1 << whereFrom)) == 0)
     268    return false;
    282269    // No longer used for original purpose - so use for ever run at all JJF
    283270#ifndef JJF_ONE
    284     // Don't run if hot start or no rows!
    285     if (model_ && (model_->hotstartSolution()||!model_->getNumRows()))
    286         return false;
    287     else
    288         return true;
     271  // Don't run if hot start or no rows!
     272  if (model_ && (model_->hotstartSolution() || !model_->getNumRows()))
     273    return false;
     274  else
     275    return true;
    289276#else
    290277#ifdef JJF_ZERO
    291     const CbcNode* currentNode = model_->currentNode();
    292     if (currentNode == NULL) {
    293         return false;
    294     }
    295 
    296     debugNodes();
    297 //   return false;
    298 
    299     const int depth = currentNode->depth();
     278  const CbcNode *currentNode = model_->currentNode();
     279  if (currentNode == NULL) {
     280    return false;
     281  }
     282
     283  debugNodes();
     284  //   return false;
     285
     286  const int depth = currentNode->depth();
    300287#else
    301     int depth = model_->currentDepth();
    302 #endif
    303 
    304     const int nodeCount = model_->getNodeCount(); // FIXME: check that this is
    305     // correct in parallel
    306 
    307     if (nodeCount == 0 || depth <= shallowDepth_) {
    308         // what to do when we are in the shallow part of the tree
    309         if (model_->getCurrentPassNumber() == 1) {
    310             // first time in the node...
    311             numInvocationsInShallow_ = 0;
    312         }
    313         ++numInvocationsInShallow_;
    314         // Very large howOftenShallow_ will give the original test:
    315         // (model_->getCurrentPassNumber() != 1)
    316         //    if ((numInvocationsInShallow_ % howOftenShallow_) != 1) {
    317         if ((numInvocationsInShallow_ % howOftenShallow_) != 0) {
    318             return false;
    319         }
    320         // LL: should we save these nodes in the list of nodes where the heur was
    321         // LL: run?
     288  int depth = model_->currentDepth();
     289#endif
     290
     291  const int nodeCount = model_->getNodeCount(); // FIXME: check that this is
     292  // correct in parallel
     293
     294  if (nodeCount == 0 || depth <= shallowDepth_) {
     295    // what to do when we are in the shallow part of the tree
     296    if (model_->getCurrentPassNumber() == 1) {
     297      // first time in the node...
     298      numInvocationsInShallow_ = 0;
     299    }
     300    ++numInvocationsInShallow_;
     301    // Very large howOftenShallow_ will give the original test:
     302    // (model_->getCurrentPassNumber() != 1)
     303    //    if ((numInvocationsInShallow_ % howOftenShallow_) != 1) {
     304    if ((numInvocationsInShallow_ % howOftenShallow_) != 0) {
     305      return false;
     306    }
     307    // LL: should we save these nodes in the list of nodes where the heur was
     308    // LL: run?
    322309#ifndef JJF_ONE
    323         if (currentNode != NULL) {
    324             // Get where we are and create the appropriate CbcHeuristicNode object
    325             CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_);
    326             runNodes_.append(nodeDesc);
    327         }
    328 #endif
    329     } else {
    330         // deeper in the tree
    331         if (model_->getCurrentPassNumber() == 1) {
    332             // first time in the node...
    333             ++numInvocationsInDeep_;
    334         }
    335         if (numInvocationsInDeep_ - lastRunDeep_ < howOften_) {
    336             return false;
    337         }
    338         if (model_->getCurrentPassNumber() > 1) {
    339             // Run the heuristic only when first entering the node.
    340             // LL: I don't think this is right. It should run just before strong
    341             // LL: branching, I believe.
    342             return false;
    343         }
    344         // Get where we are and create the appropriate CbcHeuristicNode object
    345         CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_);
    346         //#ifdef PRINT_DEBUG
     310    if (currentNode != NULL) {
     311      // Get where we are and create the appropriate CbcHeuristicNode object
     312      CbcHeuristicNode *nodeDesc = new CbcHeuristicNode(*model_);
     313      runNodes_.append(nodeDesc);
     314    }
     315#endif
     316  } else {
     317    // deeper in the tree
     318    if (model_->getCurrentPassNumber() == 1) {
     319      // first time in the node...
     320      ++numInvocationsInDeep_;
     321    }
     322    if (numInvocationsInDeep_ - lastRunDeep_ < howOften_) {
     323      return false;
     324    }
     325    if (model_->getCurrentPassNumber() > 1) {
     326      // Run the heuristic only when first entering the node.
     327      // LL: I don't think this is right. It should run just before strong
     328      // LL: branching, I believe.
     329      return false;
     330    }
     331    // Get where we are and create the appropriate CbcHeuristicNode object
     332    CbcHeuristicNode *nodeDesc = new CbcHeuristicNode(*model_);
     333    //#ifdef PRINT_DEBUG
    347334#ifndef JJF_ONE
    348         const double minDistanceToRun = 1.5 * log((double)depth) / log((double)2);
     335    const double minDistanceToRun = 1.5 * log((double)depth) / log((double)2);
    349336#else
    350337    const double minDistanceToRun = minDistanceToRun_;
    351338#endif
    352339#ifdef PRINT_DEBUG
    353         double minDistance = nodeDesc->minDistance(runNodes_);
    354         std::cout << "minDistance = " << minDistance
    355                   << ", minDistanceToRun = " << minDistanceToRun << std::endl;
    356 #endif
    357         if (nodeDesc->minDistanceIsSmall(runNodes_, minDistanceToRun)) {
    358             delete nodeDesc;
    359             return false;
    360         }
    361         runNodes_.append(nodeDesc);
    362         lastRunDeep_ = numInvocationsInDeep_;
    363         //    ++lastRunDeep_;
    364     }
    365     ++numRuns_;
    366     return true;
    367 #endif
    368 }
    369 
    370 bool
    371 CbcHeuristic::shouldHeurRun_randomChoice()
    372 {
    373     if (!when_)
    374         return false;
    375     int depth = model_->currentDepth();
    376     // when_ -999 is special marker to force to run
    377     if (depth != 0 && when_ != -999) {
    378         const double numerator = depth * depth;
    379         const double denominator = exp(depth * log(2.0));
    380         double probability = numerator / denominator;
    381         double randomNumber = randomNumberGenerator_.randomDouble();
    382         int when = when_ % 100;
    383         if (when > 2 && when < 8) {
    384             /* JJF adjustments
     340    double minDistance = nodeDesc->minDistance(runNodes_);
     341    std::cout << "minDistance = " << minDistance
     342              << ", minDistanceToRun = " << minDistanceToRun << std::endl;
     343#endif
     344    if (nodeDesc->minDistanceIsSmall(runNodes_, minDistanceToRun)) {
     345      delete nodeDesc;
     346      return false;
     347    }
     348    runNodes_.append(nodeDesc);
     349    lastRunDeep_ = numInvocationsInDeep_;
     350    //    ++lastRunDeep_;
     351  }
     352  ++numRuns_;
     353  return true;
     354#endif
     355}
     356
     357bool CbcHeuristic::shouldHeurRun_randomChoice()
     358{
     359  if (!when_)
     360    return false;
     361  int depth = model_->currentDepth();
     362  // when_ -999 is special marker to force to run
     363  if (depth != 0 && when_ != -999) {
     364    const double numerator = depth * depth;
     365    const double denominator = exp(depth * log(2.0));
     366    double probability = numerator / denominator;
     367    double randomNumber = randomNumberGenerator_.randomDouble();
     368    int when = when_ % 100;
     369    if (when > 2 && when < 8) {
     370      /* JJF adjustments
    385371            3 only at root and if no solution
    386372            4 only at root and if this heuristic has not got solution
     
    389375            7 run up to 2 times if solution found 4 otherwise
    390376            */
    391             switch (when) {
    392             case 3:
    393             default:
    394                 if (model_->bestSolution())
    395                     probability = -1.0;
    396                 break;
    397             case 4:
    398                 if (numberSolutionsFound_)
    399                     probability = -1.0;
    400                 break;
    401             case 5:
    402                 assert (decayFactor_);
    403                 if (model_->bestSolution()) {
    404                     probability = -1.0;
    405                 } else if (numCouldRun_ > 1000) {
    406                     decayFactor_ *= 0.99;
    407                     probability *= decayFactor_;
    408                 }
    409                 break;
    410             case 6:
    411                 if (depth >= 3) {
    412                     if ((numCouldRun_ % howOften_) == 0 &&
    413                             numberSolutionsFound_*howOften_ < numCouldRun_) {
    414                       //#define COIN_DEVELOP
     377      switch (when) {
     378      case 3:
     379      default:
     380        if (model_->bestSolution())
     381          probability = -1.0;
     382        break;
     383      case 4:
     384        if (numberSolutionsFound_)
     385          probability = -1.0;
     386        break;
     387      case 5:
     388        assert(decayFactor_);
     389        if (model_->bestSolution()) {
     390          probability = -1.0;
     391        } else if (numCouldRun_ > 1000) {
     392          decayFactor_ *= 0.99;
     393          probability *= decayFactor_;
     394        }
     395        break;
     396      case 6:
     397        if (depth >= 3) {
     398          if ((numCouldRun_ % howOften_) == 0 && numberSolutionsFound_ * howOften_ < numCouldRun_) {
     399            //#define COIN_DEVELOP
    415400#ifdef COIN_DEVELOP
    416                         int old = howOften_;
    417 #endif
    418                         howOften_ = CoinMin(CoinMax(static_cast<int> (howOften_ * 1.1), howOften_ + 1), 1000000);
     401            int old = howOften_;
     402#endif
     403            howOften_ = CoinMin(CoinMax(static_cast<int>(howOften_ * 1.1), howOften_ + 1), 1000000);
    419404#ifdef COIN_DEVELOP
    420                         printf("Howoften changed from %d to %d for %s\n",
    421                                old, howOften_, heuristicName_.c_str());
    422 #endif
    423                     }
    424                     probability = 1.0 / howOften_;
    425                     if (model_->bestSolution())
    426                         probability *= 0.5;
    427                 } else {
    428                     probability = 1.1;
    429                 }
    430                 break;
    431             case 7:
    432                 if ((model_->bestSolution() && numRuns_ >= 2) || numRuns_ >= 4)
    433                     probability = -1.0;
    434                 break;
    435             }
    436         }
    437         if (randomNumber > probability)
    438             return false;
    439 
    440         if (model_->getCurrentPassNumber() > 1)
    441             return false;
     405            printf("Howoften changed from %d to %d for %s\n",
     406              old, howOften_, heuristicName_.c_str());
     407#endif
     408          }
     409          probability = 1.0 / howOften_;
     410          if (model_->bestSolution())
     411            probability *= 0.5;
     412        } else {
     413          probability = 1.1;
     414        }
     415        break;
     416      case 7:
     417        if ((model_->bestSolution() && numRuns_ >= 2) || numRuns_ >= 4)
     418          probability = -1.0;
     419        break;
     420      }
     421    }
     422    if (randomNumber > probability)
     423      return false;
     424
     425    if (model_->getCurrentPassNumber() > 1)
     426      return false;
    442427#ifdef COIN_DEVELOP
    443         printf("Running %s, random %g probability %g\n",
    444                heuristicName_.c_str(), randomNumber, probability);
    445 #endif
    446     } else {
     428    printf("Running %s, random %g probability %g\n",
     429      heuristicName_.c_str(), randomNumber, probability);
     430#endif
     431  } else {
    447432#ifdef COIN_DEVELOP
    448         printf("Running %s, depth %d when %d\n",
    449                heuristicName_.c_str(), depth, when_);
    450 #endif
    451     }
    452     ++numRuns_;
    453     return true;
     433    printf("Running %s, depth %d when %d\n",
     434      heuristicName_.c_str(), depth, when_);
     435#endif
     436  }
     437  ++numRuns_;
     438  return true;
    454439}
    455440
    456441// Resets stuff if model changes
    457 void
    458 CbcHeuristic::resetModel(CbcModel * model)
    459 {
    460     model_ = model;
     442void CbcHeuristic::resetModel(CbcModel *model)
     443{
     444  model_ = model;
    461445}
    462446// Set seed
    463 void
    464 CbcHeuristic::setSeed(int value)
    465 {
    466     if (value==0) {
    467       double time = fabs(CoinGetTimeOfDay());
    468       while (time>=COIN_INT_MAX)
    469         time *= 0.5;
    470       value = static_cast<int>(time);
    471       char printArray[100];
    472       sprintf(printArray, "using time of day seed was changed from %d to %d",
    473               randomNumberGenerator_.getSeed(), value);
    474       if (model_)
    475         model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
    476           << printArray
    477           << CoinMessageEol;
    478     }
    479     randomNumberGenerator_.setSeed(value);
     447void CbcHeuristic::setSeed(int value)
     448{
     449  if (value == 0) {
     450    double time = fabs(CoinGetTimeOfDay());
     451    while (time >= COIN_INT_MAX)
     452      time *= 0.5;
     453    value = static_cast<int>(time);
     454    char printArray[100];
     455    sprintf(printArray, "using time of day seed was changed from %d to %d",
     456      randomNumberGenerator_.getSeed(), value);
     457    if (model_)
     458      model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
     459        << printArray
     460        << CoinMessageEol;
     461  }
     462  randomNumberGenerator_.setSeed(value);
    480463}
    481464// Get seed
    482 int
    483 CbcHeuristic::getSeed() const
     465int CbcHeuristic::getSeed() const
    484466{
    485467  return randomNumberGenerator_.getSeed();
     
    487469
    488470// Create C++ lines to get to current state
    489 void
    490 CbcHeuristic::generateCpp( FILE * fp, const char * heuristic)
    491 {
    492     // hard coded as CbcHeuristic virtual
    493     if (when_ != 2)
    494         fprintf(fp, "3  %s.setWhen(%d);\n", heuristic, when_);
    495     else
    496         fprintf(fp, "4  %s.setWhen(%d);\n", heuristic, when_);
    497     if (numberNodes_ != 200)
    498         fprintf(fp, "3  %s.setNumberNodes(%d);\n", heuristic, numberNodes_);
    499     else
    500         fprintf(fp, "4  %s.setNumberNodes(%d);\n", heuristic, numberNodes_);
    501     if (feasibilityPumpOptions_ != -1)
    502         fprintf(fp, "3  %s.setFeasibilityPumpOptions(%d);\n", heuristic, feasibilityPumpOptions_);
    503     else
    504         fprintf(fp, "4  %s.setFeasibilityPumpOptions(%d);\n", heuristic, feasibilityPumpOptions_);
    505     if (fractionSmall_ != 1.0)
    506         fprintf(fp, "3  %s.setFractionSmall(%g);\n", heuristic, fractionSmall_);
    507     else
    508         fprintf(fp, "4  %s.setFractionSmall(%g);\n", heuristic, fractionSmall_);
    509     if (heuristicName_ != "Unknown")
    510         fprintf(fp, "3  %s.setHeuristicName(\"%s\");\n",
    511                 heuristic, heuristicName_.c_str()) ;
    512     else
    513         fprintf(fp, "4  %s.setHeuristicName(\"%s\");\n",
    514                 heuristic, heuristicName_.c_str()) ;
    515     if (decayFactor_ != 0.0)
    516         fprintf(fp, "3  %s.setDecayFactor(%g);\n", heuristic, decayFactor_);
    517     else
    518         fprintf(fp, "4  %s.setDecayFactor(%g);\n", heuristic, decayFactor_);
    519     if (switches_ != 0)
    520         fprintf(fp, "3  %s.setSwitches(%d);\n", heuristic, switches_);
    521     else
    522         fprintf(fp, "4  %s.setSwitches(%d);\n", heuristic, switches_);
    523     if (whereFrom_ != DEFAULT_WHERE)
    524         fprintf(fp, "3  %s.setWhereFrom(%d);\n", heuristic, whereFrom_);
    525     else
    526         fprintf(fp, "4  %s.setWhereFrom(%d);\n", heuristic, whereFrom_);
    527     if (shallowDepth_ != 1)
    528         fprintf(fp, "3  %s.setShallowDepth(%d);\n", heuristic, shallowDepth_);
    529     else
    530         fprintf(fp, "4  %s.setShallowDepth(%d);\n", heuristic, shallowDepth_);
    531     if (howOftenShallow_ != 1)
    532         fprintf(fp, "3  %s.setHowOftenShallow(%d);\n", heuristic, howOftenShallow_);
    533     else
    534         fprintf(fp, "4  %s.setHowOftenShallow(%d);\n", heuristic, howOftenShallow_);
    535     if (minDistanceToRun_ != 1)
    536         fprintf(fp, "3  %s.setMinDistanceToRun(%d);\n", heuristic, minDistanceToRun_);
    537     else
    538         fprintf(fp, "4  %s.setMinDistanceToRun(%d);\n", heuristic, minDistanceToRun_);
     471void CbcHeuristic::generateCpp(FILE *fp, const char *heuristic)
     472{
     473  // hard coded as CbcHeuristic virtual
     474  if (when_ != 2)
     475    fprintf(fp, "3  %s.setWhen(%d);\n", heuristic, when_);
     476  else
     477    fprintf(fp, "4  %s.setWhen(%d);\n", heuristic, when_);
     478  if (numberNodes_ != 200)
     479    fprintf(fp, "3  %s.setNumberNodes(%d);\n", heuristic, numberNodes_);
     480  else
     481    fprintf(fp, "4  %s.setNumberNodes(%d);\n", heuristic, numberNodes_);
     482  if (feasibilityPumpOptions_ != -1)
     483    fprintf(fp, "3  %s.setFeasibilityPumpOptions(%d);\n", heuristic, feasibilityPumpOptions_);
     484  else
     485    fprintf(fp, "4  %s.setFeasibilityPumpOptions(%d);\n", heuristic, feasibilityPumpOptions_);
     486  if (fractionSmall_ != 1.0)
     487    fprintf(fp, "3  %s.setFractionSmall(%g);\n", heuristic, fractionSmall_);
     488  else
     489    fprintf(fp, "4  %s.setFractionSmall(%g);\n", heuristic, fractionSmall_);
     490  if (heuristicName_ != "Unknown")
     491    fprintf(fp, "3  %s.setHeuristicName(\"%s\");\n",
     492      heuristic, heuristicName_.c_str());
     493  else
     494    fprintf(fp, "4  %s.setHeuristicName(\"%s\");\n",
     495      heuristic, heuristicName_.c_str());
     496  if (decayFactor_ != 0.0)
     497    fprintf(fp, "3  %s.setDecayFactor(%g);\n", heuristic, decayFactor_);
     498  else
     499    fprintf(fp, "4  %s.setDecayFactor(%g);\n", heuristic, decayFactor_);
     500  if (switches_ != 0)
     501    fprintf(fp, "3  %s.setSwitches(%d);\n", heuristic, switches_);
     502  else
     503    fprintf(fp, "4  %s.setSwitches(%d);\n", heuristic, switches_);
     504  if (whereFrom_ != DEFAULT_WHERE)
     505    fprintf(fp, "3  %s.setWhereFrom(%d);\n", heuristic, whereFrom_);
     506  else
     507    fprintf(fp, "4  %s.setWhereFrom(%d);\n", heuristic, whereFrom_);
     508  if (shallowDepth_ != 1)
     509    fprintf(fp, "3  %s.setShallowDepth(%d);\n", heuristic, shallowDepth_);
     510  else
     511    fprintf(fp, "4  %s.setShallowDepth(%d);\n", heuristic, shallowDepth_);
     512  if (howOftenShallow_ != 1)
     513    fprintf(fp, "3  %s.setHowOftenShallow(%d);\n", heuristic, howOftenShallow_);
     514  else
     515    fprintf(fp, "4  %s.setHowOftenShallow(%d);\n", heuristic, howOftenShallow_);
     516  if (minDistanceToRun_ != 1)
     517    fprintf(fp, "3  %s.setMinDistanceToRun(%d);\n", heuristic, minDistanceToRun_);
     518  else
     519    fprintf(fp, "4  %s.setMinDistanceToRun(%d);\n", heuristic, minDistanceToRun_);
    539520}
    540521// Destructor
    541 CbcHeuristic::~CbcHeuristic ()
    542 {
    543     delete [] inputSolution_;
     522CbcHeuristic::~CbcHeuristic()
     523{
     524  delete[] inputSolution_;
    544525}
    545526
    546527// update model
    547 void CbcHeuristic::setModel(CbcModel * model)
    548 {
    549     model_ = model;
     528void CbcHeuristic::setModel(CbcModel *model)
     529{
     530  model_ = model;
    550531}
    551532/* Clone but ..
     
    556537CbcHeuristic::cloneBut(int type)
    557538{
    558     OsiSolverInterface * solver;
    559     if ((type&1) == 0 || !model_->continuousSolver())
    560         solver = model_->solver()->clone();
    561     else
    562         solver = model_->continuousSolver()->clone();
     539  OsiSolverInterface *solver;
     540  if ((type & 1) == 0 || !model_->continuousSolver())
     541    solver = model_->solver()->clone();
     542  else
     543    solver = model_->continuousSolver()->clone();
    563544#ifdef COIN_HAS_CLP
    564     OsiClpSolverInterface * clpSolver
    565     = dynamic_cast<OsiClpSolverInterface *> (solver);
    566 #endif
    567     if ((type&2) != 0) {
    568         int n = model_->numberObjects();
    569         int priority = model_->continuousPriority();
    570         if (priority < COIN_INT_MAX) {
    571             for (int i = 0; i < n; i++) {
    572                 const OsiObject * obj = model_->object(i);
    573                 const CbcSimpleInteger * thisOne =
    574                     dynamic_cast <const CbcSimpleInteger *> (obj);
    575                 if (thisOne) {
    576                     int iColumn = thisOne->columnNumber();
    577                     if (thisOne->priority() >= priority)
    578                         solver->setContinuous(iColumn);
    579                 }
    580             }
    581         }
     545  OsiClpSolverInterface *clpSolver
     546    = dynamic_cast<OsiClpSolverInterface *>(solver);
     547#endif
     548  if ((type & 2) != 0) {
     549    int n = model_->numberObjects();
     550    int priority = model_->continuousPriority();
     551    if (priority < COIN_INT_MAX) {
     552      for (int i = 0; i < n; i++) {
     553        const OsiObject *obj = model_->object(i);
     554        const CbcSimpleInteger *thisOne = dynamic_cast<const CbcSimpleInteger *>(obj);
     555        if (thisOne) {
     556          int iColumn = thisOne->columnNumber();
     557          if (thisOne->priority() >= priority)
     558            solver->setContinuous(iColumn);
     559        }
     560      }
     561    }
    582562#ifdef COIN_HAS_CLP
    583         if (clpSolver) {
    584             for (int i = 0; i < n; i++) {
    585                 const OsiObject * obj = model_->object(i);
    586                 const CbcSimpleInteger * thisOne =
    587                     dynamic_cast <const CbcSimpleInteger *> (obj);
    588                 if (thisOne) {
    589                     int iColumn = thisOne->columnNumber();
    590                     if (clpSolver->isOptionalInteger(iColumn))
    591                         clpSolver->setContinuous(iColumn);
    592                 }
    593             }
    594         }
    595 #endif
    596     }
     563    if (clpSolver) {
     564      for (int i = 0; i < n; i++) {
     565        const OsiObject *obj = model_->object(i);
     566        const CbcSimpleInteger *thisOne = dynamic_cast<const CbcSimpleInteger *>(obj);
     567        if (thisOne) {
     568          int iColumn = thisOne->columnNumber();
     569          if (clpSolver->isOptionalInteger(iColumn))
     570            clpSolver->setContinuous(iColumn);
     571        }
     572      }
     573    }
     574#endif
     575  }
    597576#ifdef COIN_HAS_CLP
    598     if ((type&4) != 0 && clpSolver) {
    599         int options = clpSolver->getModelPtr()->moreSpecialOptions();
    600         clpSolver->getModelPtr()->setMoreSpecialOptions(options | 64);
    601     }
    602 #endif
    603     return solver;
     577  if ((type & 4) != 0 && clpSolver) {
     578    int options = clpSolver->getModelPtr()->moreSpecialOptions();
     579    clpSolver->getModelPtr()->setMoreSpecialOptions(options | 64);
     580  }
     581#endif
     582  return solver;
    604583}
    605584// Whether to exit at once on gap
    606 bool
    607 CbcHeuristic::exitNow(double bestObjective) const
    608 {
    609     if ((switches_&2048) != 0) {
    610         // exit may be forced - but unset for next time
    611         switches_ &= ~2048;
    612         if ((switches_&1024) != 0)
    613             return true;
    614     } else if ((switches_&1) == 0) {
    615         return false;
    616     }
    617     // See if can stop on gap
    618     OsiSolverInterface * solver = model_->solver();
    619     double bestPossibleObjective = solver->getObjValue() * solver->getObjSense();
    620     double absGap = CoinMax(model_->getAllowableGap(),
    621                             model_->getHeuristicGap());
    622     double fracGap = CoinMax(model_->getAllowableFractionGap(),
    623                              model_->getHeuristicFractionGap());
    624     double testGap = CoinMax(absGap, fracGap *
    625                              CoinMax(fabs(bestObjective),
    626                                      fabs(bestPossibleObjective)));
    627 
    628     if (bestObjective - bestPossibleObjective < testGap
    629             && model_->getCutoffIncrement() >= 0.0) {
    630         return true;
    631     } else {
    632         return false;
    633     }
     585bool CbcHeuristic::exitNow(double bestObjective) const
     586{
     587  if ((switches_ & 2048) != 0) {
     588    // exit may be forced - but unset for next time
     589    switches_ &= ~2048;
     590    if ((switches_ & 1024) != 0)
     591      return true;
     592  } else if ((switches_ & 1) == 0) {
     593    return false;
     594  }
     595  // See if can stop on gap
     596  OsiSolverInterface *solver = model_->solver();
     597  double bestPossibleObjective = solver->getObjValue() * solver->getObjSense();
     598  double absGap = CoinMax(model_->getAllowableGap(),
     599    model_->getHeuristicGap());
     600  double fracGap = CoinMax(model_->getAllowableFractionGap(),
     601    model_->getHeuristicFractionGap());
     602  double testGap = CoinMax(absGap, fracGap * CoinMax(fabs(bestObjective), fabs(bestPossibleObjective)));
     603
     604  if (bestObjective - bestPossibleObjective < testGap
     605    && model_->getCutoffIncrement() >= 0.0) {
     606    return true;
     607  } else {
     608    return false;
     609  }
    634610}
    635611#ifdef HISTORY_STATISTICS
     
    637613#endif
    638614static double sizeRatio(int numberRowsNow, int numberColumnsNow,
    639                         int numberRowsStart, int numberColumnsStart)
    640 {
    641     double valueNow;
    642     if (numberRowsNow*10 > numberColumnsNow || numberColumnsNow < 200) {
    643         valueNow = 2 * numberRowsNow + numberColumnsNow;
    644     } else {
    645         // long and thin - rows are more important
    646         if (numberRowsNow*40 > numberColumnsNow)
    647             valueNow = 10 * numberRowsNow + numberColumnsNow;
    648         else
    649             valueNow = 200 * numberRowsNow + numberColumnsNow;
    650     }
    651     double valueStart;
    652     if (numberRowsStart*10 > numberColumnsStart || numberColumnsStart < 200) {
    653         valueStart = 2 * numberRowsStart + numberColumnsStart;
    654     } else {
    655         // long and thin - rows are more important
    656         if (numberRowsStart*40 > numberColumnsStart)
    657             valueStart = 10 * numberRowsStart + numberColumnsStart;
    658         else
    659             valueStart = 200 * numberRowsStart + numberColumnsStart;
    660     }
    661     //printf("sizeProblem Now %g, %d rows, %d columns\nsizeProblem Start %g, %d rows, %d columns\n",
    662     // valueNow,numberRowsNow,numberColumnsNow,
    663     // valueStart,numberRowsStart,numberColumnsStart);
    664     if (10*numberRowsNow < 8*numberRowsStart || 10*numberColumnsNow < 7*numberColumnsStart)
    665         return valueNow / valueStart;
    666     else if (10*numberRowsNow < 9*numberRowsStart)
    667         return 1.1*(valueNow / valueStart);
    668     else if (numberRowsNow < numberRowsStart)
    669         return 1.5*(valueNow / valueStart);
     615  int numberRowsStart, int numberColumnsStart)
     616{
     617  double valueNow;
     618  if (numberRowsNow * 10 > numberColumnsNow || numberColumnsNow < 200) {
     619    valueNow = 2 * numberRowsNow + numberColumnsNow;
     620  } else {
     621    // long and thin - rows are more important
     622    if (numberRowsNow * 40 > numberColumnsNow)
     623      valueNow = 10 * numberRowsNow + numberColumnsNow;
    670624    else
    671         return 2.0*(valueNow / valueStart);
     625      valueNow = 200 * numberRowsNow + numberColumnsNow;
     626  }
     627  double valueStart;
     628  if (numberRowsStart * 10 > numberColumnsStart || numberColumnsStart < 200) {
     629    valueStart = 2 * numberRowsStart + numberColumnsStart;
     630  } else {
     631    // long and thin - rows are more important
     632    if (numberRowsStart * 40 > numberColumnsStart)
     633      valueStart = 10 * numberRowsStart + numberColumnsStart;
     634    else
     635      valueStart = 200 * numberRowsStart + numberColumnsStart;
     636  }
     637  //printf("sizeProblem Now %g, %d rows, %d columns\nsizeProblem Start %g, %d rows, %d columns\n",
     638  // valueNow,numberRowsNow,numberColumnsNow,
     639  // valueStart,numberRowsStart,numberColumnsStart);
     640  if (10 * numberRowsNow < 8 * numberRowsStart || 10 * numberColumnsNow < 7 * numberColumnsStart)
     641    return valueNow / valueStart;
     642  else if (10 * numberRowsNow < 9 * numberRowsStart)
     643    return 1.1 * (valueNow / valueStart);
     644  else if (numberRowsNow < numberRowsStart)
     645    return 1.5 * (valueNow / valueStart);
     646  else
     647    return 2.0 * (valueNow / valueStart);
    672648}
    673649
    674650//static int saveModel=0;
    675651// Do mini branch and bound (return 1 if solution)
    676 int
    677 CbcHeuristic::smallBranchAndBound(OsiSolverInterface * solver, int numberNodes,
    678                                   double * newSolution, double & newSolutionValue,
    679                                   double cutoff, std::string name) const
    680 {
    681   CbcEventHandler *eventHandler = model_->getEventHandler() ;
     652int CbcHeuristic::smallBranchAndBound(OsiSolverInterface *solver, int numberNodes,
     653  double *newSolution, double &newSolutionValue,
     654  double cutoff, std::string name) const
     655{
     656  CbcEventHandler *eventHandler = model_->getEventHandler();
    682657  // Use this fraction
    683658  double fractionSmall = fractionSmall_;
    684   int maximumSolutions =  model_->getMaximumSolutions();
     659  int maximumSolutions = model_->getMaximumSolutions();
    685660  int iterationMultiplier = 100;
    686661  if (eventHandler) {
     
    688663      double fractionSmall;
    689664      double spareDouble[3];
    690       OsiSolverInterface * solver;
    691       void * sparePointer[2];
     665      OsiSolverInterface *solver;
     666      void *sparePointer[2];
    692667      int numberNodes;
    693668      int maximumSolutions;
     
    697672    } SmallMod;
    698673    SmallMod temp;
    699     temp.solver=solver;
    700     temp.fractionSmall=fractionSmall;
    701     temp.numberNodes=numberNodes;
    702     temp.iterationMultiplier=iterationMultiplier;
    703     temp.howOften=howOften_;
    704     temp.maximumSolutions=maximumSolutions;
    705     CbcEventHandler::CbcAction status =
    706       eventHandler->event(CbcEventHandler::smallBranchAndBound,
    707                           &temp);
    708     if (status==CbcEventHandler::killSolution)
     674    temp.solver = solver;
     675    temp.fractionSmall = fractionSmall;
     676    temp.numberNodes = numberNodes;
     677    temp.iterationMultiplier = iterationMultiplier;
     678    temp.howOften = howOften_;
     679    temp.maximumSolutions = maximumSolutions;
     680    CbcEventHandler::CbcAction status = eventHandler->event(CbcEventHandler::smallBranchAndBound,
     681      &temp);
     682    if (status == CbcEventHandler::killSolution)
    709683      return -1;
    710     if (status==CbcEventHandler::takeAction) {
    711       fractionSmall=temp.fractionSmall;
    712       numberNodes=temp.numberNodes;
    713       iterationMultiplier=temp.iterationMultiplier;
    714       howOften_=temp.howOften;
    715       maximumSolutions=temp.maximumSolutions;
     684    if (status == CbcEventHandler::takeAction) {
     685      fractionSmall = temp.fractionSmall;
     686      numberNodes = temp.numberNodes;
     687      iterationMultiplier = temp.iterationMultiplier;
     688      howOften_ = temp.howOften;
     689      maximumSolutions = temp.maximumSolutions;
    716690    }
    717691  }
     
    722696  }
    723697#endif
    724     // size before
    725     int shiftRows = 0;
    726     if (numberNodes < 0)
    727         shiftRows = solver->getNumRows() - numberNodes_;
    728     int numberRowsStart = solver->getNumRows() - shiftRows;
    729     int numberColumnsStart = solver->getNumCols();
     698  // size before
     699  int shiftRows = 0;
     700  if (numberNodes < 0)
     701    shiftRows = solver->getNumRows() - numberNodes_;
     702  int numberRowsStart = solver->getNumRows() - shiftRows;
     703  int numberColumnsStart = solver->getNumCols();
    730704#ifdef CLP_INVESTIGATE
    731     printf("%s has %d rows, %d columns\n",
    732            name.c_str(), solver->getNumRows(), solver->getNumCols());
    733 #endif
    734     double before = 2 * numberRowsStart + numberColumnsStart;
    735     if (before > 40000.0) {
    736         // fairly large - be more conservative
    737         double multiplier = 1.0 - 0.3 * CoinMin(100000.0, before - 40000.0) / 100000.0;
    738         if (multiplier < 1.0) {
    739             fractionSmall *= multiplier;
     705  printf("%s has %d rows, %d columns\n",
     706    name.c_str(), solver->getNumRows(), solver->getNumCols());
     707#endif
     708  double before = 2 * numberRowsStart + numberColumnsStart;
     709  if (before > 40000.0) {
     710    // fairly large - be more conservative
     711    double multiplier = 1.0 - 0.3 * CoinMin(100000.0, before - 40000.0) / 100000.0;
     712    if (multiplier < 1.0) {
     713      fractionSmall *= multiplier;
    740714#ifdef CLP_INVESTIGATE
    741             printf("changing fractionSmall from %g to %g for %s\n",
    742                    fractionSmall_, fractionSmall, name.c_str());
    743 #endif
    744         }
    745     }
     715      printf("changing fractionSmall from %g to %g for %s\n",
     716        fractionSmall_, fractionSmall, name.c_str());
     717#endif
     718    }
     719  }
    746720#ifdef COIN_HAS_CLP
    747     OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
    748     if (clpSolver && (clpSolver->specialOptions()&65536) == 0) {
    749         // go faster stripes
    750         if (clpSolver->getNumRows() < 300 && clpSolver->getNumCols() < 500) {
    751             clpSolver->setupForRepeatedUse(2, 0);
    752         } else {
    753             clpSolver->setupForRepeatedUse(0, 0);
    754         }
    755         // Turn this off if you get problems
    756         // Used to be automatically set
    757         clpSolver->setSpecialOptions(clpSolver->specialOptions() | (128 + 64 - 128));
    758         ClpSimplex * lpSolver = clpSolver->getModelPtr();
    759         lpSolver->setSpecialOptions(lpSolver->specialOptions() | 0x01000000); // say is Cbc (and in branch and bound)
    760         lpSolver->setSpecialOptions(lpSolver->specialOptions() |
    761                                     (/*16384+*/4096 + 512 + 128));
    762     }
     721  OsiClpSolverInterface *clpSolver = dynamic_cast<OsiClpSolverInterface *>(solver);
     722  if (clpSolver && (clpSolver->specialOptions() & 65536) == 0) {
     723    // go faster stripes
     724    if (clpSolver->getNumRows() < 300 && clpSolver->getNumCols() < 500) {
     725      clpSolver->setupForRepeatedUse(2, 0);
     726    } else {
     727      clpSolver->setupForRepeatedUse(0, 0);
     728    }
     729    // Turn this off if you get problems
     730    // Used to be automatically set
     731    clpSolver->setSpecialOptions(clpSolver->specialOptions() | (128 + 64 - 128));
     732    ClpSimplex *lpSolver = clpSolver->getModelPtr();
     733    lpSolver->setSpecialOptions(lpSolver->specialOptions() | 0x01000000); // say is Cbc (and in branch and bound)
     734    lpSolver->setSpecialOptions(lpSolver->specialOptions() | (/*16384+*/ 4096 + 512 + 128));
     735  }
    763736#endif
    764737#ifdef HISTORY_STATISTICS
    765     getHistoryStatistics_ = false;
     738  getHistoryStatistics_ = false;
    766739#endif
    767740#ifdef COIN_DEVELOP
    768     int status = 0;
    769 #endif
    770     int logLevel = model_->logLevel();
     741  int status = 0;
     742#endif
     743  int logLevel = model_->logLevel();
    771744#define LEN_PRINT 250
    772     char generalPrint[LEN_PRINT];
    773     // Do presolve to see if possible
    774     int numberColumns = solver->getNumCols();
    775     char * reset = NULL;
    776     int returnCode = 1;
    777     int saveModelOptions = model_->specialOptions();
    778     //assert ((saveModelOptions&2048) == 0);
    779     model_->setSpecialOptions(saveModelOptions | 2048);
    780     if (fractionSmall<1.0) {
    781         int saveLogLevel = solver->messageHandler()->logLevel();
    782         if (saveLogLevel == 1)
    783             solver->messageHandler()->setLogLevel(0);
    784         OsiPresolve * pinfo = new OsiPresolve();
    785         int presolveActions = 0;
    786         // Allow dual stuff on integers
    787         presolveActions = 1;
    788         // Do not allow all +1 to be tampered with
    789         //if (allPlusOnes)
    790         //presolveActions |= 2;
    791         // allow transfer of costs
    792         // presolveActions |= 4;
    793         pinfo->setPresolveActions(presolveActions);
    794         OsiSolverInterface * presolvedModel = pinfo->presolvedModel(*solver, 1.0e-8, true, 2);
    795         delete pinfo;
    796         // see if too big
    797 
    798         if (presolvedModel) {
    799             int afterRows = presolvedModel->getNumRows();
    800             int afterCols = presolvedModel->getNumCols();
    801             //#define COIN_DEVELOP
     745  char generalPrint[LEN_PRINT];
     746  // Do presolve to see if possible
     747  int numberColumns = solver->getNumCols();
     748  char *reset = NULL;
     749  int returnCode = 1;
     750  int saveModelOptions = model_->specialOptions();
     751  //assert ((saveModelOptions&2048) == 0);
     752  model_->setSpecialOptions(saveModelOptions | 2048);
     753  if (fractionSmall < 1.0) {
     754    int saveLogLevel = solver->messageHandler()->logLevel();
     755    if (saveLogLevel == 1)
     756      solver->messageHandler()->setLogLevel(0);
     757    OsiPresolve *pinfo = new OsiPresolve();
     758    int presolveActions = 0;
     759    // Allow dual stuff on integers
     760    presolveActions = 1;
     761    // Do not allow all +1 to be tampered with
     762    //if (allPlusOnes)
     763    //presolveActions |= 2;
     764    // allow transfer of costs
     765    // presolveActions |= 4;
     766    pinfo->setPresolveActions(presolveActions);
     767    OsiSolverInterface *presolvedModel = pinfo->presolvedModel(*solver, 1.0e-8, true, 2);
     768    delete pinfo;
     769    // see if too big
     770
     771    if (presolvedModel) {
     772      int afterRows = presolvedModel->getNumRows();
     773      int afterCols = presolvedModel->getNumCols();
     774      //#define COIN_DEVELOP
    802775#ifdef COIN_DEVELOP_z
    803             if (numberNodes < 0) {
    804                 solver->writeMpsNative("before.mps", NULL, NULL, 2, 1);
    805                 presolvedModel->writeMpsNative("after1.mps", NULL, NULL, 2, 1);
     776      if (numberNodes < 0) {
     777        solver->writeMpsNative("before.mps", NULL, NULL, 2, 1);
     778        presolvedModel->writeMpsNative("after1.mps", NULL, NULL, 2, 1);
     779      }
     780#endif
     781      delete presolvedModel;
     782      double ratio = sizeRatio(afterRows - shiftRows, afterCols,
     783        numberRowsStart, numberColumnsStart);
     784      double after = 2 * afterRows + afterCols;
     785      if (ratio > fractionSmall && after > 300 && numberNodes >= 0) {
     786        // Need code to try again to compress further using used
     787        const int *used = model_->usedInSolution();
     788        int maxUsed = 0;
     789        int iColumn;
     790        const double *lower = solver->getColLower();
     791        const double *upper = solver->getColUpper();
     792        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     793          if (upper[iColumn] > lower[iColumn]) {
     794            if (solver->isBinary(iColumn))
     795              maxUsed = CoinMax(maxUsed, used[iColumn]);
     796          }
     797        }
     798        if (maxUsed) {
     799          reset = new char[numberColumns];
     800          int nFix = 0;
     801          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     802            reset[iColumn] = 0;
     803            if (upper[iColumn] > lower[iColumn]) {
     804              if (solver->isBinary(iColumn) && used[iColumn] == maxUsed) {
     805                bool setValue = true;
     806                if (maxUsed == 1) {
     807                  double randomNumber = randomNumberGenerator_.randomDouble();
     808                  if (randomNumber > 0.3)
     809                    setValue = false;
     810                }
     811                if (setValue) {
     812                  reset[iColumn] = 1;
     813                  solver->setColLower(iColumn, 1.0);
     814                  nFix++;
     815                }
     816              }
    806817            }
    807 #endif
     818          }
     819          pinfo = new OsiPresolve();
     820          presolveActions = 0;
     821          // Allow dual stuff on integers
     822          presolveActions = 1;
     823          // Do not allow all +1 to be tampered with
     824          //if (allPlusOnes)
     825          //presolveActions |= 2;
     826          // allow transfer of costs
     827          // presolveActions |= 4;
     828          pinfo->setPresolveActions(presolveActions);
     829          presolvedModel = pinfo->presolvedModel(*solver, 1.0e-8, true, 2);
     830          delete pinfo;
     831          if (presolvedModel) {
     832            // see if too big
     833            int afterRows2 = presolvedModel->getNumRows();
     834            int afterCols2 = presolvedModel->getNumCols();
    808835            delete presolvedModel;
    809             double ratio = sizeRatio(afterRows - shiftRows, afterCols,
    810                                      numberRowsStart, numberColumnsStart);
    811             double after = 2 * afterRows + afterCols;
    812             if (ratio > fractionSmall && after > 300 && numberNodes >= 0) {
    813                 // Need code to try again to compress further using used
    814                 const int * used =  model_->usedInSolution();
    815                 int maxUsed = 0;
    816                 int iColumn;
    817                 const double * lower = solver->getColLower();
    818                 const double * upper = solver->getColUpper();
    819                 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    820                     if (upper[iColumn] > lower[iColumn]) {
    821                         if (solver->isBinary(iColumn))
    822                             maxUsed = CoinMax(maxUsed, used[iColumn]);
    823                     }
    824                 }
    825                 if (maxUsed) {
    826                     reset = new char [numberColumns];
    827                     int nFix = 0;
    828                     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    829                         reset[iColumn] = 0;
    830                         if (upper[iColumn] > lower[iColumn]) {
    831                             if (solver->isBinary(iColumn) && used[iColumn] == maxUsed) {
    832                                 bool setValue = true;
    833                                 if (maxUsed == 1) {
    834                                     double randomNumber = randomNumberGenerator_.randomDouble();
    835                                     if (randomNumber > 0.3)
    836                                         setValue = false;
    837                                 }
    838                                 if (setValue) {
    839                                     reset[iColumn] = 1;
    840                                     solver->setColLower(iColumn, 1.0);
    841                                     nFix++;
    842                                 }
    843                             }
    844                         }
    845                     }
    846                     pinfo = new OsiPresolve();
    847                     presolveActions = 0;
    848                     // Allow dual stuff on integers
    849                     presolveActions = 1;
    850                     // Do not allow all +1 to be tampered with
    851                     //if (allPlusOnes)
    852                     //presolveActions |= 2;
    853                     // allow transfer of costs
    854                     // presolveActions |= 4;
    855                     pinfo->setPresolveActions(presolveActions);
    856                     presolvedModel = pinfo->presolvedModel(*solver, 1.0e-8, true, 2);
    857                     delete pinfo;
    858                     if (presolvedModel) {
    859                         // see if too big
    860                         int afterRows2 = presolvedModel->getNumRows();
    861                         int afterCols2 = presolvedModel->getNumCols();
    862                         delete presolvedModel;
    863                         double ratio = sizeRatio(afterRows2 - shiftRows, afterCols2,
    864                                                  numberRowsStart, numberColumnsStart);
    865                         double after = 2 * afterRows2 + afterCols2;
    866                         if (ratio > fractionSmall && (after > 300 || numberNodes < 0)) {
    867                             sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - %d fixed gives %d, %d - still too large",
    868                                     solver->getNumRows(), solver->getNumCols(),
    869                                     afterRows, afterCols, nFix, afterRows2, afterCols2);
    870                             // If much too big - give up
    871                             if (ratio > 0.75)
    872                                 returnCode = -1;
    873                         } else {
    874                             sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - %d fixed gives %d, %d - ok now",
    875                                     solver->getNumRows(), solver->getNumCols(),
    876                                     afterRows, afterCols, nFix, afterRows2, afterCols2);
    877                         }
    878                         model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
    879                         << generalPrint
    880                         << CoinMessageEol;
    881                     } else {
    882                         returnCode = 2; // infeasible
    883                     }
    884                 }
    885             } else if (ratio > fractionSmall && after > 300 && numberNodes >=0) {
     836            double ratio = sizeRatio(afterRows2 - shiftRows, afterCols2,
     837              numberRowsStart, numberColumnsStart);
     838            double after = 2 * afterRows2 + afterCols2;
     839            if (ratio > fractionSmall && (after > 300 || numberNodes < 0)) {
     840              sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - %d fixed gives %d, %d - still too large",
     841                solver->getNumRows(), solver->getNumCols(),
     842                afterRows, afterCols, nFix, afterRows2, afterCols2);
     843              // If much too big - give up
     844              if (ratio > 0.75)
    886845                returnCode = -1;
     846            } else {
     847              sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - %d fixed gives %d, %d - ok now",
     848                solver->getNumRows(), solver->getNumCols(),
     849                afterRows, afterCols, nFix, afterRows2, afterCols2);
    887850            }
    888         } else {
     851            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
     852              << generalPrint
     853              << CoinMessageEol;
     854          } else {
    889855            returnCode = 2; // infeasible
    890         }
    891         solver->messageHandler()->setLogLevel(saveLogLevel);
    892     }
    893     if (returnCode == 2 || returnCode == -1) {
    894         model_->setSpecialOptions(saveModelOptions);
    895         delete [] reset;
     856          }
     857        }
     858      } else if (ratio > fractionSmall && after > 300 && numberNodes >= 0) {
     859        returnCode = -1;
     860      }
     861    } else {
     862      returnCode = 2; // infeasible
     863    }
     864    solver->messageHandler()->setLogLevel(saveLogLevel);
     865  }
     866  if (returnCode == 2 || returnCode == -1) {
     867    model_->setSpecialOptions(saveModelOptions);
     868    delete[] reset;
    896869#ifdef HISTORY_STATISTICS
    897         getHistoryStatistics_ = true;
    898 #endif
    899         //printf("small no good\n");
    900         return returnCode;
    901     }
    902     // Reduce printout
    903     bool takeHint;
    904     OsiHintStrength strength;
    905     solver->getHintParam(OsiDoReducePrint, takeHint, strength);
    906     solver->setHintParam(OsiDoReducePrint, true, OsiHintTry);
    907     solver->setHintParam(OsiDoPresolveInInitial, false, OsiHintTry);
    908     double signedCutoff = cutoff*solver->getObjSense();
    909     solver->setDblParam(OsiDualObjectiveLimit, signedCutoff);
    910     solver->initialSolve();
    911     if (solver->isProvenOptimal()) {
    912         CglPreProcess process;
    913         OsiSolverInterface * solver2 = NULL;
    914         if ((model_->moreSpecialOptions()&65536)!=0)
    915           process.setOptions(2+4+8+16); // no cuts
    916         else
    917           process.setOptions(16); // no complicated dupcol stuff
    918         /* Do not try and produce equality cliques and
     870    getHistoryStatistics_ = true;
     871#endif
     872    //printf("small no good\n");
     873    return returnCode;
     874  }
     875  // Reduce printout
     876  bool takeHint;
     877  OsiHintStrength strength;
     878  solver->getHintParam(OsiDoReducePrint, takeHint, strength);
     879  solver->setHintParam(OsiDoReducePrint, true, OsiHintTry);
     880  solver->setHintParam(OsiDoPresolveInInitial, false, OsiHintTry);
     881  double signedCutoff = cutoff * solver->getObjSense();
     882  solver->setDblParam(OsiDualObjectiveLimit, signedCutoff);
     883  solver->initialSolve();
     884  if (solver->isProvenOptimal()) {
     885    CglPreProcess process;
     886    OsiSolverInterface *solver2 = NULL;
     887    if ((model_->moreSpecialOptions() & 65536) != 0)
     888      process.setOptions(2 + 4 + 8 + 16); // no cuts
     889    else
     890      process.setOptions(16); // no complicated dupcol stuff
     891    /* Do not try and produce equality cliques and
    919892           do up to 2 passes (normally) 5 if restart */
    920         int numberPasses = 2;
    921         if ((model_->moreSpecialOptions2()&16)!=0) {
    922           // quick
    923           process.setOptions(2+4+8+16); // no cuts
    924           numberPasses = 1;
    925         }
    926         if (numberNodes < 0) {
    927           numberPasses = 5;
    928           // Say some rows cuts
    929           int numberRows = solver->getNumRows();
    930           if (numberNodes_ < numberRows && true /* think */) {
    931             char * type = new char[numberRows];
    932             memset(type, 0, numberNodes_);
    933             memset(type + numberNodes_, 1, numberRows - numberNodes_);
    934             process.passInRowTypes(type, numberRows);
    935             delete [] type;
    936           }
    937         }
    938         if (logLevel <= 1)
    939           process.messageHandler()->setLogLevel(0);
    940         if (!solver->defaultHandler()&&
    941             solver->messageHandler()->logLevel(0)!=-1000)
    942           process.passInMessageHandler(solver->messageHandler());
     893    int numberPasses = 2;
     894    if ((model_->moreSpecialOptions2() & 16) != 0) {
     895      // quick
     896      process.setOptions(2 + 4 + 8 + 16); // no cuts
     897      numberPasses = 1;
     898    }
     899    if (numberNodes < 0) {
     900      numberPasses = 5;
     901      // Say some rows cuts
     902      int numberRows = solver->getNumRows();
     903      if (numberNodes_ < numberRows && true /* think */) {
     904        char *type = new char[numberRows];
     905        memset(type, 0, numberNodes_);
     906        memset(type + numberNodes_, 1, numberRows - numberNodes_);
     907        process.passInRowTypes(type, numberRows);
     908        delete[] type;
     909      }
     910    }
     911    if (logLevel <= 1)
     912      process.messageHandler()->setLogLevel(0);
     913    if (!solver->defaultHandler() && solver->messageHandler()->logLevel(0) != -1000)
     914      process.passInMessageHandler(solver->messageHandler());
    943915#ifdef CGL_DEBUG
    944         /*
     916    /*
    945917          We're debugging. (specialOptions 1)
    946918        */
    947         if ((model_->specialOptions()&1) != 0) {
    948           const OsiRowCutDebugger *debugger = solver->getRowCutDebugger() ;
    949           if (debugger) {
    950             process.setApplicationData(const_cast<double *>(debugger->optimalSolution()));
    951           }
    952         }
     919    if ((model_->specialOptions() & 1) != 0) {
     920      const OsiRowCutDebugger *debugger = solver->getRowCutDebugger();
     921      if (debugger) {
     922        process.setApplicationData(const_cast<double *>(debugger->optimalSolution()));
     923      }
     924    }
    953925#endif
    954926#ifdef COIN_HAS_CLP
    955         OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
    956         // See if SOS
    957         if (clpSolver&&clpSolver->numberSOS()) {
    958           // SOS
    959           int numberSOS = clpSolver->numberSOS();
    960           const CoinSet * setInfo = clpSolver->setInfo();
    961           int *sosStart = new int [numberSOS+1];
    962           char *sosType = new char [numberSOS];
    963           int i;
    964           int nTotal = 0;
    965           sosStart[0] = 0;
    966           for ( i = 0; i < numberSOS; i++) {
    967             int type = setInfo[i].setType();
    968             int n = setInfo[i].numberEntries();
    969             sosType[i] = static_cast<char>(type);
    970             nTotal += n;
    971             sosStart[i+1] = nTotal;
    972           }
    973           int * sosIndices = new int[nTotal];
    974           double * sosReference = new double [nTotal];
    975           for (i = 0; i < numberSOS; i++) {
    976             int n = setInfo[i].numberEntries();
    977             const int * which = setInfo[i].which();
    978             const double * weights = setInfo[i].weights();
    979             int base = sosStart[i];
    980             for (int j = 0; j < n; j++) {
    981               int k = which[j];
    982               sosIndices[j+base] = k;
    983               sosReference[j+base] = weights ? weights[j] : static_cast<double> (j);
    984             }
    985           }
    986           int numberColumns = solver->getNumCols();
    987           char * prohibited = new char[numberColumns];
    988           memset(prohibited, 0, numberColumns);
    989           int n = sosStart[numberSOS];
    990           for (int i = 0; i < n; i++) {
    991             int iColumn = sosIndices[i];
    992             prohibited[iColumn] = 1;
    993           }
    994           delete [] sosIndices;
    995           delete [] sosReference;
    996           delete [] sosStart;
    997           delete [] sosType;
    998           process.passInProhibited(prohibited, numberColumns);
    999           delete [] prohibited;
    1000         }
    1001 #endif
    1002         solver2 = process.preProcessNonDefault(*solver, false,
    1003                                                numberPasses);
    1004           if (!solver2) {
    1005             if (logLevel > 1)
    1006               printf("Pre-processing says infeasible\n");
    1007             returnCode = 2; // so will be infeasible
    1008           } else {
     927    OsiClpSolverInterface *clpSolver = dynamic_cast<OsiClpSolverInterface *>(solver);
     928    // See if SOS
     929    if (clpSolver && clpSolver->numberSOS()) {
     930      // SOS
     931      int numberSOS = clpSolver->numberSOS();
     932      const CoinSet *setInfo = clpSolver->setInfo();
     933      int *sosStart = new int[numberSOS + 1];
     934      char *sosType = new char[numberSOS];
     935      int i;
     936      int nTotal = 0;
     937      sosStart[0] = 0;
     938      for (i = 0; i < numberSOS; i++) {
     939        int type = setInfo[i].setType();
     940        int n = setInfo[i].numberEntries();
     941        sosType[i] = static_cast<char>(type);
     942        nTotal += n;
     943        sosStart[i + 1] = nTotal;
     944      }
     945      int *sosIndices = new int[nTotal];
     946      double *sosReference = new double[nTotal];
     947      for (i = 0; i < numberSOS; i++) {
     948        int n = setInfo[i].numberEntries();
     949        const int *which = setInfo[i].which();
     950        const double *weights = setInfo[i].weights();
     951        int base = sosStart[i];
     952        for (int j = 0; j < n; j++) {
     953          int k = which[j];
     954          sosIndices[j + base] = k;
     955          sosReference[j + base] = weights ? weights[j] : static_cast<double>(j);
     956        }
     957      }
     958      int numberColumns = solver->getNumCols();
     959      char *prohibited = new char[numberColumns];
     960      memset(prohibited, 0, numberColumns);
     961      int n = sosStart[numberSOS];
     962      for (int i = 0; i < n; i++) {
     963        int iColumn = sosIndices[i];
     964        prohibited[iColumn] = 1;
     965      }
     966      delete[] sosIndices;
     967      delete[] sosReference;
     968      delete[] sosStart;
     969      delete[] sosType;
     970      process.passInProhibited(prohibited, numberColumns);
     971      delete[] prohibited;
     972    }
     973#endif
     974    solver2 = process.preProcessNonDefault(*solver, false,
     975      numberPasses);
     976    if (!solver2) {
     977      if (logLevel > 1)
     978        printf("Pre-processing says infeasible\n");
     979      returnCode = 2; // so will be infeasible
     980    } else {
    1009981#ifdef COIN_DEVELOP_z
    1010             if (numberNodes < 0) {
    1011               solver2->writeMpsNative("after2.mps", NULL, NULL, 2, 1);
     982      if (numberNodes < 0) {
     983        solver2->writeMpsNative("after2.mps", NULL, NULL, 2, 1);
     984      }
     985#endif
     986      // see if too big
     987      double ratio = sizeRatio(solver2->getNumRows() - shiftRows, solver2->getNumCols(),
     988        numberRowsStart, numberColumnsStart);
     989      double after = 2 * solver2->getNumRows() + solver2->getNumCols();
     990      if (ratio > fractionSmall && (after > 300 || numberNodes < 0)) {
     991        sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - too large",
     992          solver->getNumRows(), solver->getNumCols(),
     993          solver2->getNumRows(), solver2->getNumCols());
     994        model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
     995          << generalPrint
     996          << CoinMessageEol;
     997        returnCode = -1;
     998        //printf("small no good2\n");
     999      } else {
     1000        sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns",
     1001          solver->getNumRows(), solver->getNumCols(),
     1002          solver2->getNumRows(), solver2->getNumCols());
     1003        model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
     1004          << generalPrint
     1005          << CoinMessageEol;
     1006      }
     1007#ifdef CGL_DEBUG
     1008      if ((model_->specialOptions() & 1) != 0) {
     1009        const OsiRowCutDebugger *debugger = solver2->getRowCutDebugger();
     1010        if (debugger) {
     1011          printf("On optimal path after preprocessing\n");
     1012        }
     1013      }
     1014#endif
     1015      if (returnCode == 1) {
     1016        solver2->resolve();
     1017        CbcModel model(*solver2);
     1018        double startTime = model_->getDblParam(CbcModel::CbcStartSeconds);
     1019        model.setDblParam(CbcModel::CbcStartSeconds, startTime);
     1020        // move seed across
     1021        model.randomNumberGenerator()->setSeed(model_->randomNumberGenerator()->getSeed());
     1022#ifdef COIN_HAS_CLP
     1023        // redo SOS
     1024        OsiClpSolverInterface *clpSolver
     1025          = dynamic_cast<OsiClpSolverInterface *>(model.solver());
     1026        if (clpSolver && clpSolver->numberSOS()) {
     1027          int numberColumns = clpSolver->getNumCols();
     1028          const int *originalColumns = process.originalColumns();
     1029          CoinSet *setInfo = const_cast<CoinSet *>(clpSolver->setInfo());
     1030          int numberSOS = clpSolver->numberSOS();
     1031          for (int iSOS = 0; iSOS < numberSOS; iSOS++) {
     1032            //int type = setInfo[iSOS].setType();
     1033            int n = setInfo[iSOS].numberEntries();
     1034            int *which = setInfo[iSOS].modifiableWhich();
     1035            double *weights = setInfo[iSOS].modifiableWeights();
     1036            int n2 = 0;
     1037            for (int j = 0; j < n; j++) {
     1038              int iColumn = which[j];
     1039              int i;
     1040              for (i = 0; i < numberColumns; i++) {
     1041                if (originalColumns[i] == iColumn)
     1042                  break;
     1043              }
     1044              if (i < numberColumns) {
     1045                which[n2] = i;
     1046                weights[n2++] = weights[j];
     1047              }
    10121048            }
    1013 #endif
    1014             // see if too big
    1015             double ratio = sizeRatio(solver2->getNumRows() - shiftRows, solver2->getNumCols(),
    1016                                      numberRowsStart, numberColumnsStart);
    1017             double after = 2 * solver2->getNumRows() + solver2->getNumCols();
    1018             if (ratio > fractionSmall && (after > 300 || numberNodes < 0)) {
    1019                 sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - too large",
    1020                         solver->getNumRows(), solver->getNumCols(),
    1021                         solver2->getNumRows(), solver2->getNumCols());
    1022                 model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
    1023                 << generalPrint
    1024                 << CoinMessageEol;
    1025                 returnCode = -1;
    1026                 //printf("small no good2\n");
    1027             } else {
    1028                 sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns",
    1029                         solver->getNumRows(), solver->getNumCols(),
    1030                         solver2->getNumRows(), solver2->getNumCols());
    1031                 model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
    1032                 << generalPrint
    1033                 << CoinMessageEol;
    1034             }
    1035 #ifdef CGL_DEBUG
    1036             if ((model_->specialOptions()&1) != 0) {
    1037               const OsiRowCutDebugger *debugger = solver2->getRowCutDebugger() ;
    1038               if (debugger) {
    1039                 printf("On optimal path after preprocessing\n");
    1040               }
    1041             }
    1042 #endif
    1043             if (returnCode == 1) {
    1044                 solver2->resolve();
    1045                 CbcModel model(*solver2);
    1046                 double startTime=model_->getDblParam(CbcModel::CbcStartSeconds);
    1047                 model.setDblParam(CbcModel::CbcStartSeconds,startTime);
    1048                 // move seed across
    1049                 model.randomNumberGenerator()->setSeed(model_->randomNumberGenerator()->getSeed());
    1050 #ifdef COIN_HAS_CLP
    1051                 // redo SOS
    1052                 OsiClpSolverInterface * clpSolver
    1053                   = dynamic_cast<OsiClpSolverInterface *> (model.solver());
    1054                 if (clpSolver && clpSolver->numberSOS()) {
    1055                   int numberColumns = clpSolver->getNumCols();
    1056                   const int * originalColumns = process.originalColumns();
    1057                   CoinSet * setInfo =
    1058                     const_cast<CoinSet *>(clpSolver->setInfo());
    1059                   int numberSOS = clpSolver->numberSOS();
    1060                   for (int iSOS = 0; iSOS < numberSOS; iSOS++) {
    1061                     //int type = setInfo[iSOS].setType();
    1062                     int n = setInfo[iSOS].numberEntries();
    1063                     int * which = setInfo[iSOS].modifiableWhich();
    1064                     double * weights = setInfo[iSOS].modifiableWeights();
    1065                     int n2=0;
    1066                     for (int j=0;j<n;j++) {
    1067                       int iColumn=which[j];
    1068                       int i;
    1069                       for (i = 0; i < numberColumns; i++) {
    1070                         if (originalColumns[i] == iColumn)
    1071                           break;
    1072                       }
    1073                       if (i < numberColumns) {
    1074                         which[n2] = i;
    1075                         weights[n2++] = weights[j];
    1076                       }
    1077                     }
    1078                     setInfo[iSOS].setNumberEntries(n2);
    1079                   }
    1080                 }
    1081 #endif
    1082                 if (numberNodes >= 0) {
    1083                     // normal
    1084                     model.setSpecialOptions(saveModelOptions | 2048);
    1085                     if (logLevel <= 1 && feasibilityPumpOptions_ != -3)
    1086                         model.setLogLevel(0);
    1087                     else
    1088                         model.setLogLevel(logLevel);
    1089                     // No small fathoming
    1090                     model.setFastNodeDepth(-1);
    1091                     model.setCutoff(signedCutoff);
    1092                     model.setStrongStrategy(0);
    1093                     // Don't do if original fraction > 1.0 and too large
    1094                     if (fractionSmall_>1.0 && fractionSmall_ < 1000000.0) {
    1095                       /* 1.4 means -1 nodes if >.4
     1049            setInfo[iSOS].setNumberEntries(n2);
     1050          }
     1051        }
     1052#endif
     1053        if (numberNodes >= 0) {
     1054          // normal
     1055          model.setSpecialOptions(saveModelOptions | 2048);
     1056          if (logLevel <= 1 && feasibilityPumpOptions_ != -3)
     1057            model.setLogLevel(0);
     1058          else
     1059            model.setLogLevel(logLevel);
     1060          // No small fathoming
     1061          model.setFastNodeDepth(-1);
     1062          model.setCutoff(signedCutoff);
     1063          model.setStrongStrategy(0);
     1064          // Don't do if original fraction > 1.0 and too large
     1065          if (fractionSmall_ > 1.0 && fractionSmall_ < 1000000.0) {
     1066            /* 1.4 means -1 nodes if >.4
    10961067                         2.4 means -1 nodes if >.5 and 0 otherwise
    10971068                         3.4 means -1 nodes if >.6 and 0 or 5
    10981069                         4.4 means -1 nodes if >.7 and 0, 5 or 10
    10991070                      */
    1100                       double fraction = fractionSmall_-floor(fractionSmall_);
    1101                       if (ratio>fraction) {
    1102                         int type = static_cast<int>(floor(fractionSmall_*0.1));
    1103                         int over = static_cast<int>(ceil(ratio-fraction));
    1104                         int maxNodes[]={-1,0,5,10};
    1105                         if (type>over)
    1106                           numberNodes=maxNodes[type-over];
    1107                         else
    1108                           numberNodes=-1;
    1109                       }
    1110                     }
    1111                     model.setMaximumNodes(numberNodes);
    1112                     model.solver()->setHintParam(OsiDoReducePrint, true, OsiHintTry);
    1113                     if ((saveModelOptions&2048) == 0)
    1114                       model.setMoreSpecialOptions(model_->moreSpecialOptions());
    1115                     model.setMoreSpecialOptions2(model_->moreSpecialOptions2());
    1116                     // off conflict analysis
    1117                     model.setMoreSpecialOptions(model.moreSpecialOptions()&~4194304);
    1118                    
    1119                     // Lightweight
    1120                     CbcStrategyDefaultSubTree strategy(model_, 1, 5, 1, 0);
    1121                     model.setStrategy(strategy);
    1122                     model.solver()->setIntParam(OsiMaxNumIterationHotStart, 10);
    1123                     model.setMaximumCutPassesAtRoot(CoinMin(20, CoinAbs(model_->getMaximumCutPassesAtRoot())));
    1124                     model.setMaximumCutPasses(CoinMin(10, model_->getMaximumCutPasses()));
    1125                     // Set best solution (even if bad for this submodel)
    1126                     if (model_->bestSolution()) {
    1127                       const double * bestSolution = model_->bestSolution();
    1128                       int numberColumns2 = model.solver()->getNumCols();
    1129                       double * bestSolution2 = new double [numberColumns2];
    1130                       const int * originalColumns = process.originalColumns();
    1131                       for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
    1132                         int jColumn = originalColumns[iColumn];
    1133                         bestSolution2[iColumn] = bestSolution[jColumn];
    1134                       }
    1135                       model.setBestSolution(bestSolution2,numberColumns2,
    1136                                             1.0e50,
    1137                                             false);
    1138                       model.setSolutionCount(1);
    1139                       maximumSolutions++;
    1140                       delete [] bestSolution2;
    1141                     }
    1142                 } else {
    1143                     // modify for event handler
    1144                     model.setSpecialOptions(saveModelOptions);
    1145                     model_->messageHandler()->message(CBC_RESTART, model_->messages())
    1146                     << solver2->getNumRows() << solver2->getNumCols()
    1147                     << CoinMessageEol;
    1148                     // going for full search and copy across more stuff
    1149                     model.gutsOfCopy(*model_, 2);
     1071            double fraction = fractionSmall_ - floor(fractionSmall_);
     1072            if (ratio > fraction) {
     1073              int type = static_cast<int>(floor(fractionSmall_ * 0.1));
     1074              int over = static_cast<int>(ceil(ratio - fraction));
     1075              int maxNodes[] = { -1, 0, 5, 10 };
     1076              if (type > over)
     1077                numberNodes = maxNodes[type - over];
     1078              else
     1079                numberNodes = -1;
     1080            }
     1081          }
     1082          model.setMaximumNodes(numberNodes);
     1083          model.solver()->setHintParam(OsiDoReducePrint, true, OsiHintTry);
     1084          if ((saveModelOptions & 2048) == 0)
     1085            model.setMoreSpecialOptions(model_->moreSpecialOptions());
     1086          model.setMoreSpecialOptions2(model_->moreSpecialOptions2());
     1087          // off conflict analysis
     1088          model.setMoreSpecialOptions(model.moreSpecialOptions() & ~4194304);
     1089
     1090          // Lightweight
     1091          CbcStrategyDefaultSubTree strategy(model_, 1, 5, 1, 0);
     1092          model.setStrategy(strategy);
     1093          model.solver()->setIntParam(OsiMaxNumIterationHotStart, 10);
     1094          model.setMaximumCutPassesAtRoot(CoinMin(20, CoinAbs(model_->getMaximumCutPassesAtRoot())));
     1095          model.setMaximumCutPasses(CoinMin(10, model_->getMaximumCutPasses()));
     1096          // Set best solution (even if bad for this submodel)
     1097          if (model_->bestSolution()) {
     1098            const double *bestSolution = model_->bestSolution();
     1099            int numberColumns2 = model.solver()->getNumCols();
     1100            double *bestSolution2 = new double[numberColumns2];
     1101            const int *originalColumns = process.originalColumns();
     1102            for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
     1103              int jColumn = originalColumns[iColumn];
     1104              bestSolution2[iColumn] = bestSolution[jColumn];
     1105            }
     1106            model.setBestSolution(bestSolution2, numberColumns2,
     1107              1.0e50,
     1108              false);
     1109            model.setSolutionCount(1);
     1110            maximumSolutions++;
     1111            delete[] bestSolution2;
     1112          }
     1113        } else {
     1114          // modify for event handler
     1115          model.setSpecialOptions(saveModelOptions);
     1116          model_->messageHandler()->message(CBC_RESTART, model_->messages())
     1117            << solver2->getNumRows() << solver2->getNumCols()
     1118            << CoinMessageEol;
     1119          // going for full search and copy across more stuff
     1120          model.gutsOfCopy(*model_, 2);
    11501121#ifdef CGL_DEBUG
    1151                     if ((model_->specialOptions()&1) != 0) {
    1152                       const OsiRowCutDebugger *debugger = model.solver()->getRowCutDebugger() ;
    1153                       if (debugger) {
    1154                         printf("On optimal path BB\n");
    1155                       }
    1156                     }
    1157 #endif
    1158                     assert (!model_->heuristicModel());
    1159                     model_->setHeuristicModel(&model);
    1160                     for (int i = 0; i < model.numberCutGenerators(); i++) {
    1161                         CbcCutGenerator * generator = model.cutGenerator(i);
    1162                         CglGomory * gomory = dynamic_cast<CglGomory *>
    1163                           (generator->generator());
    1164                         if (gomory&&gomory->originalSolver())
    1165                           gomory->passInOriginalSolver(model.solver());
    1166                         generator->setTiming(true);
    1167                         // Turn on if was turned on
    1168                         int iOften = model_->cutGenerator(i)->howOften();
     1122          if ((model_->specialOptions() & 1) != 0) {
     1123            const OsiRowCutDebugger *debugger = model.solver()->getRowCutDebugger();
     1124            if (debugger) {
     1125              printf("On optimal path BB\n");
     1126            }
     1127          }
     1128#endif
     1129          assert(!model_->heuristicModel());
     1130          model_->setHeuristicModel(&model);
     1131          for (int i = 0; i < model.numberCutGenerators(); i++) {
     1132            CbcCutGenerator *generator = model.cutGenerator(i);
     1133            CglGomory *gomory = dynamic_cast<CglGomory *>(generator->generator());
     1134            if (gomory && gomory->originalSolver())
     1135              gomory->passInOriginalSolver(model.solver());
     1136            generator->setTiming(true);
     1137            // Turn on if was turned on
     1138            int iOften = model_->cutGenerator(i)->howOften();
    11691139#ifdef CLP_INVESTIGATE
    1170                         printf("Gen %d often %d %d\n",
    1171                                i, generator->howOften(),
    1172                                iOften);
    1173 #endif
    1174                         if (iOften > 0)
    1175                             generator->setHowOften(iOften % 1000000);
    1176                         if (model_->cutGenerator(i)->howOftenInSub() == -200)
    1177                             generator->setHowOften(-100);
    1178                     }
    1179                     model.setCutoff(signedCutoff);
    1180                     // make sure can't do nested search! but allow heuristics
    1181                     model.setSpecialOptions((model.specialOptions()&(~(512 + 2048))) | 1024);
    1182                     // but say we are doing full search
    1183                     model.setSpecialOptions(model.specialOptions()|67108864);
    1184                     bool takeHint;
    1185                     OsiHintStrength strength;
    1186                     // Switch off printing if asked to
    1187                     model_->solver()->getHintParam(OsiDoReducePrint, takeHint, strength);
    1188                     model.solver()->setHintParam(OsiDoReducePrint, takeHint, strength);
    1189                     // no cut generators if none in parent
    1190                     CbcStrategyDefault
    1191                       strategy(model_->numberCutGenerators() ? 1 : -1,
    1192                                model_->numberStrong(),
    1193                                model_->numberBeforeTrust());
    1194                     // Set up pre-processing - no
    1195                     strategy.setupPreProcessing(0); // was (4);
    1196                     model.setStrategy(strategy);
    1197                     //model.solver()->writeMps("crunched");
    1198                     int numberCuts = process.cuts().sizeRowCuts();
    1199                     if (numberCuts) {
    1200                         // add in cuts
    1201                         CglStored cuts = process.cuts();
    1202                         model.addCutGenerator(&cuts, 1, "Stored from first");
    1203                         model.cutGenerator(model.numberCutGenerators()-1)->setGlobalCuts(true);
    1204                     }
    1205                 }
    1206                 // Do search
    1207                 if (logLevel > 1)
    1208                     model_->messageHandler()->message(CBC_START_SUB, model_->messages())
    1209                     << name
    1210                     << model.getMaximumNodes()
    1211                     << CoinMessageEol;
    1212                 // probably faster to use a basis to get integer solutions
    1213                 model.setSpecialOptions(model.specialOptions() | 2);
     1140            printf("Gen %d often %d %d\n",
     1141              i, generator->howOften(),
     1142              iOften);
     1143#endif
     1144            if (iOften > 0)
     1145              generator->setHowOften(iOften % 1000000);
     1146            if (model_->cutGenerator(i)->howOftenInSub() == -200)
     1147              generator->setHowOften(-100);
     1148          }
     1149          model.setCutoff(signedCutoff);
     1150          // make sure can't do nested search! but allow heuristics
     1151          model.setSpecialOptions((model.specialOptions() & (~(512 + 2048))) | 1024);
     1152          // but say we are doing full search
     1153          model.setSpecialOptions(model.specialOptions() | 67108864);
     1154          bool takeHint;
     1155          OsiHintStrength strength;
     1156          // Switch off printing if asked to
     1157          model_->solver()->getHintParam(OsiDoReducePrint, takeHint, strength);
     1158          model.solver()->setHintParam(OsiDoReducePrint, takeHint, strength);
     1159          // no cut generators if none in parent
     1160          CbcStrategyDefault
     1161            strategy(model_->numberCutGenerators() ? 1 : -1,
     1162              model_->numberStrong(),
     1163              model_->numberBeforeTrust());
     1164          // Set up pre-processing - no
     1165          strategy.setupPreProcessing(0); // was (4);
     1166          model.setStrategy(strategy);
     1167          //model.solver()->writeMps("crunched");
     1168          int numberCuts = process.cuts().sizeRowCuts();
     1169          if (numberCuts) {
     1170            // add in cuts
     1171            CglStored cuts = process.cuts();
     1172            model.addCutGenerator(&cuts, 1, "Stored from first");
     1173            model.cutGenerator(model.numberCutGenerators() - 1)->setGlobalCuts(true);
     1174          }
     1175        }
     1176        // Do search
     1177        if (logLevel > 1)
     1178          model_->messageHandler()->message(CBC_START_SUB, model_->messages())
     1179            << name
     1180            << model.getMaximumNodes()
     1181            << CoinMessageEol;
     1182        // probably faster to use a basis to get integer solutions
     1183        model.setSpecialOptions(model.specialOptions() | 2);
    12141184#ifdef CBC_THREAD
    1215                 if (model_->getNumberThreads() > 0 && (model_->getThreadMode()&4) != 0) {
    1216                     // See if at root node
    1217                     bool atRoot = model_->getNodeCount() == 0;
    1218                     int passNumber = model_->getCurrentPassNumber();
    1219                     if (atRoot && passNumber == 1)
    1220                         model.setNumberThreads(model_->getNumberThreads());
    1221                 }
    1222 #endif
    1223                 model.setParentModel(*model_);
    1224                 model.setMaximumSolutions(maximumSolutions);
    1225                 model.setOriginalColumns(process.originalColumns());
    1226                 model.setSearchStrategy(-1);
    1227                 // If no feasibility pump then insert a lightweight one
    1228                 if (feasibilityPumpOptions_ >= 0 || feasibilityPumpOptions_ == -2) {
    1229                     CbcHeuristicFPump * fpump = NULL;
    1230                     for (int i = 0; i < model.numberHeuristics(); i++) {
    1231                         CbcHeuristicFPump* pump =
    1232                             dynamic_cast<CbcHeuristicFPump*>(model.heuristic(i));
    1233                         if (pump) {
    1234                             fpump = pump;
    1235                             break;
    1236                         }
    1237                     }
    1238                     if (!fpump) {
    1239                         CbcHeuristicFPump heuristic4;
    1240                         // use any cutoff
    1241                         heuristic4.setFakeCutoff(0.5*COIN_DBL_MAX);
    1242                         if (fractionSmall_<=1.0)
    1243                           heuristic4.setMaximumPasses(10);
    1244                         int pumpTune = feasibilityPumpOptions_;
    1245                         if (pumpTune==-2)
    1246                           pumpTune = 4; // proximity
    1247                         if (pumpTune > 0) {
    1248                             /*
     1185        if (model_->getNumberThreads() > 0 && (model_->getThreadMode() & 4) != 0) {
     1186          // See if at root node
     1187          bool atRoot = model_->getNodeCount() == 0;
     1188          int passNumber = model_->getCurrentPassNumber();
     1189          if (atRoot && passNumber == 1)
     1190            model.setNumberThreads(model_->getNumberThreads());
     1191        }
     1192#endif
     1193        model.setParentModel(*model_);
     1194        model.setMaximumSolutions(maximumSolutions);
     1195        model.setOriginalColumns(process.originalColumns());
     1196        model.setSearchStrategy(-1);
     1197        // If no feasibility pump then insert a lightweight one
     1198        if (feasibilityPumpOptions_ >= 0 || feasibilityPumpOptions_ == -2) {
     1199          CbcHeuristicFPump *fpump = NULL;
     1200          for (int i = 0; i < model.numberHeuristics(); i++) {
     1201            CbcHeuristicFPump *pump = dynamic_cast<CbcHeuristicFPump *>(model.heuristic(i));
     1202            if (pump) {
     1203              fpump = pump;
     1204              break;
     1205            }
     1206          }
     1207          if (!fpump) {
     1208            CbcHeuristicFPump heuristic4;
     1209            // use any cutoff
     1210            heuristic4.setFakeCutoff(0.5 * COIN_DBL_MAX);
     1211            if (fractionSmall_ <= 1.0)
     1212              heuristic4.setMaximumPasses(10);
     1213            int pumpTune = feasibilityPumpOptions_;
     1214            if (pumpTune == -2)
     1215              pumpTune = 4; // proximity
     1216            if (pumpTune > 0) {
     1217              /*
    12491218                            >=10000000 for using obj
    12501219                            >=1000000 use as accumulate switch
     
    12561225                            6 as 3 but all slack basis!
    12571226                            */
    1258                             double value = solver2->getObjSense() * solver2->getObjValue();
    1259                             int w = pumpTune / 10;
    1260                             int ix = w % 10;
    1261                             w /= 10;
    1262                             int c = w % 10;
    1263                             w /= 10;
    1264                             int r = w;
    1265                             int accumulate = r / 1000;
    1266                             r -= 1000 * accumulate;
    1267                             if (accumulate >= 10) {
    1268                                 int which = accumulate / 10;
    1269                                 accumulate -= 10 * which;
    1270                                 which--;
    1271                                 // weights and factors
    1272                                 double weight[] = {0.1, 0.1, 0.5, 0.5, 1.0, 1.0, 5.0, 5.0};
    1273                                 double factor[] = {0.1, 0.5, 0.1, 0.5, 0.1, 0.5, 0.1, 0.5};
    1274                                 heuristic4.setInitialWeight(weight[which]);
    1275                                 heuristic4.setWeightFactor(factor[which]);
    1276                             }
    1277                             // fake cutoff
    1278                             if (c) {
    1279                                 double cutoff;
    1280                                 solver2->getDblParam(OsiDualObjectiveLimit, cutoff);
    1281                                 cutoff = CoinMin(cutoff, value + 0.1 * fabs(value) * c);
    1282                                 heuristic4.setFakeCutoff(cutoff);
    1283                             }
    1284                             if (r) {
    1285                                 // also set increment
    1286                                 //double increment = (0.01*i+0.005)*(fabs(value)+1.0e-12);
    1287                                 double increment = 0.0;
    1288                                 heuristic4.setAbsoluteIncrement(increment);
    1289                                 heuristic4.setAccumulate(accumulate);
    1290                                 heuristic4.setMaximumRetries(r + 1);
    1291                             }
    1292                             pumpTune = pumpTune % 100;
    1293                             if (pumpTune == 6)
    1294                                 pumpTune = 13;
    1295                             if (pumpTune != 13)
    1296                                 pumpTune = pumpTune % 10;
    1297                             heuristic4.setWhen(pumpTune);
    1298                             if (ix) {
    1299                                 heuristic4.setFeasibilityPumpOptions(ix*10);
    1300                             }
    1301                         }
    1302                         model.addHeuristic(&heuristic4, "feasibility pump", 0);
    1303        
    1304                     }
    1305                 } else if (feasibilityPumpOptions_==-3) {
    1306                   // add all (except this)
    1307                   for (int i = 0; i < model_->numberHeuristics(); i++) {
    1308                     if (strcmp(heuristicName(),model_->heuristic(i)->heuristicName()))
    1309                       model.addHeuristic(model_->heuristic(i));
    1310                   }
    1311                 }
    1312                 // modify heuristics
    1313                 for (int i = 0; i < model.numberHeuristics(); i++) {
    1314                   // reset lastNode
    1315                   CbcHeuristicRINS * rins =
    1316                     dynamic_cast<CbcHeuristicRINS*>(model.heuristic(i));
    1317                   if (rins) {
    1318                     rins->setLastNode(-1000);
    1319                     rins->setSolutionCount(0);
    1320                   }
    1321                 }
    1322                 //printf("sol %x\n",inputSolution_);
     1227              double value = solver2->getObjSense() * solver2->getObjValue();
     1228              int w = pumpTune / 10;
     1229              int ix = w % 10;
     1230              w /= 10;
     1231              int c = w % 10;
     1232              w /= 10;
     1233              int r = w;
     1234              int accumulate = r / 1000;
     1235              r -= 1000 * accumulate;
     1236              if (accumulate >= 10) {
     1237                int which = accumulate / 10;
     1238                accumulate -= 10 * which;
     1239                which--;
     1240                // weights and factors
     1241                double weight[] = { 0.1, 0.1, 0.5, 0.5, 1.0, 1.0, 5.0, 5.0 };
     1242                double factor[] = { 0.1, 0.5, 0.1, 0.5, 0.1, 0.5, 0.1, 0.5 };
     1243                heuristic4.setInitialWeight(weight[which]);
     1244                heuristic4.setWeightFactor(factor[which]);
     1245              }
     1246              // fake cutoff
     1247              if (c) {
     1248                double cutoff;
     1249                solver2->getDblParam(OsiDualObjectiveLimit, cutoff);
     1250                cutoff = CoinMin(cutoff, value + 0.1 * fabs(value) * c);
     1251                heuristic4.setFakeCutoff(cutoff);
     1252              }
     1253              if (r) {
     1254                // also set increment
     1255                //double increment = (0.01*i+0.005)*(fabs(value)+1.0e-12);
     1256                double increment = 0.0;
     1257                heuristic4.setAbsoluteIncrement(increment);
     1258                heuristic4.setAccumulate(accumulate);
     1259                heuristic4.setMaximumRetries(r + 1);
     1260              }
     1261              pumpTune = pumpTune % 100;
     1262              if (pumpTune == 6)
     1263                pumpTune = 13;
     1264              if (pumpTune != 13)
     1265                pumpTune = pumpTune % 10;
     1266              heuristic4.setWhen(pumpTune);
     1267              if (ix) {
     1268                heuristic4.setFeasibilityPumpOptions(ix * 10);
     1269              }
     1270            }
     1271            model.addHeuristic(&heuristic4, "feasibility pump", 0);
     1272          }
     1273        } else if (feasibilityPumpOptions_ == -3) {
     1274          // add all (except this)
     1275          for (int i = 0; i < model_->numberHeuristics(); i++) {
     1276            if (strcmp(heuristicName(), model_->heuristic(i)->heuristicName()))
     1277              model.addHeuristic(model_->heuristic(i));
     1278          }
     1279        }
     1280        // modify heuristics
     1281        for (int i = 0; i < model.numberHeuristics(); i++) {
     1282          // reset lastNode
     1283          CbcHeuristicRINS *rins = dynamic_cast<CbcHeuristicRINS *>(model.heuristic(i));
     1284          if (rins) {
     1285            rins->setLastNode(-1000);
     1286            rins->setSolutionCount(0);
     1287          }
     1288        }
     1289        //printf("sol %x\n",inputSolution_);
    13231290#ifdef CGL_DEBUG
    1324                 if ((model_->specialOptions()&1) != 0) {
    1325                   const OsiRowCutDebugger *debugger = model.solver()->getRowCutDebugger() ;
    1326                   if (debugger) {
    1327                     printf("On optimal path CC\n");
    1328                   }
    1329                 }
    1330 #endif
    1331                 if (inputSolution_) {
    1332                     // translate and add a serendipity heuristic
    1333                     int numberColumns = solver2->getNumCols();
    1334                     const int * which = process.originalColumns();
    1335                     OsiSolverInterface * solver3 = solver2->clone();
    1336                     for (int i = 0; i < numberColumns; i++) {
    1337                         if (isHeuristicInteger(solver3,i)) {
    1338                             int k = which[i];
    1339                             double value = inputSolution_[k];
    1340                             //if (value)
    1341                             //printf("orig col %d now %d val %g\n",
    1342                             //       k,i,value);
    1343                             solver3->setColLower(i, value);
    1344                             solver3->setColUpper(i, value);
    1345                         }
    1346                     }
    1347                     solver3->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX);
    1348                     solver3->resolve();
    1349                     if (!solver3->isProvenOptimal()) {
    1350                         // Try just setting nonzeros
    1351                         OsiSolverInterface * solver4 = solver2->clone();
    1352                         for (int i = 0; i < numberColumns; i++) {
    1353                             if (isHeuristicInteger(solver4,i)) {
    1354                                 int k = which[i];
    1355                                 double value = floor(inputSolution_[k] + 0.5);
    1356                                 if (value) {
    1357                                     solver3->setColLower(i, value);
    1358                                     solver3->setColUpper(i, value);
    1359                                 }
    1360                             }
    1361                         }
    1362                         solver4->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX);
    1363                         solver4->resolve();
    1364                         int nBad = -1;
    1365                         if (solver4->isProvenOptimal()) {
    1366                             nBad = 0;
    1367                             const double * solution = solver4->getColSolution();
    1368                             for (int i = 0; i < numberColumns; i++) {
    1369                                 if (isHeuristicInteger(solver4,i)) {
    1370                                     double value = floor(solution[i] + 0.5);
    1371                                     if (fabs(value - solution[i]) > 1.0e-6)
    1372                                         nBad++;
    1373                                 }
    1374                             }
    1375                         }
    1376                         if (nBad) {
    1377                             delete solver4;
    1378                         } else {
    1379                             delete solver3;
    1380                             solver3 = solver4;
    1381                         }
    1382                     }
    1383                     if (solver3->isProvenOptimal()) {
    1384                         // good
    1385                         CbcSerendipity heuristic(model);
    1386                         double value = solver3->getObjSense() * solver3->getObjValue();
    1387                         heuristic.setInputSolution(solver3->getColSolution(), value);
    1388                         value = value + 1.0e-7*(1.0 + fabs(value));
    1389                         value *= solver3->getObjSense();
    1390                         model.setCutoff(value);
    1391                         model.addHeuristic(&heuristic, "Previous solution", 0);
    1392                         //printf("added seren\n");
    1393                     } else {
    1394                         double value = model_->getMinimizationObjValue();
    1395                         value = value + 1.0e-7*(1.0 + fabs(value));
    1396                         value *= solver3->getObjSense();
    1397                         model.setCutoff(value);
    1398                         sprintf(generalPrint, "Unable to insert previous solution - using cutoff of %g",
    1399                                 value);
    1400                         model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
    1401                         << generalPrint
    1402                         << CoinMessageEol;
     1291        if ((model_->specialOptions() & 1) != 0) {
     1292          const OsiRowCutDebugger *debugger = model.solver()->getRowCutDebugger();
     1293          if (debugger) {
     1294            printf("On optimal path CC\n");
     1295          }
     1296        }
     1297#endif
     1298        if (inputSolution_) {
     1299          // translate and add a serendipity heuristic
     1300          int numberColumns = solver2->getNumCols();
     1301          const int *which = process.originalColumns();
     1302          OsiSolverInterface *solver3 = solver2->clone();
     1303          for (int i = 0; i < numberColumns; i++) {
     1304            if (isHeuristicInteger(solver3, i)) {
     1305              int k = which[i];
     1306              double value = inputSolution_[k];
     1307              //if (value)
     1308              //printf("orig col %d now %d val %g\n",
     1309              //       k,i,value);
     1310              solver3->setColLower(i, value);
     1311              solver3->setColUpper(i, value);
     1312            }
     1313          }
     1314          solver3->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX);
     1315          solver3->resolve();
     1316          if (!solver3->isProvenOptimal()) {
     1317            // Try just setting nonzeros
     1318            OsiSolverInterface *solver4 = solver2->clone();
     1319            for (int i = 0; i < numberColumns; i++) {
     1320              if (isHeuristicInteger(solver4, i)) {
     1321                int k = which[i];
     1322                double value = floor(inputSolution_[k] + 0.5);
     1323                if (value) {
     1324                  solver3->setColLower(i, value);
     1325                  solver3->setColUpper(i, value);
     1326                }
     1327              }
     1328            }
     1329            solver4->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX);
     1330            solver4->resolve();
     1331            int nBad = -1;
     1332            if (solver4->isProvenOptimal()) {
     1333              nBad = 0;
     1334              const double *solution = solver4->getColSolution();
     1335              for (int i = 0; i < numberColumns; i++) {
     1336                if (isHeuristicInteger(solver4, i)) {
     1337                  double value = floor(solution[i] + 0.5);
     1338                  if (fabs(value - solution[i]) > 1.0e-6)
     1339                    nBad++;
     1340                }
     1341              }
     1342            }
     1343            if (nBad) {
     1344              delete solver4;
     1345            } else {
     1346              delete solver3;
     1347              solver3 = solver4;
     1348            }
     1349          }
     1350          if (solver3->isProvenOptimal()) {
     1351            // good
     1352            CbcSerendipity heuristic(model);
     1353            double value = solver3->getObjSense() * solver3->getObjValue();
     1354            heuristic.setInputSolution(solver3->getColSolution(), value);
     1355            value = value + 1.0e-7 * (1.0 + fabs(value));
     1356            value *= solver3->getObjSense();
     1357            model.setCutoff(value);
     1358            model.addHeuristic(&heuristic, "Previous solution", 0);
     1359            //printf("added seren\n");
     1360          } else {
     1361            double value = model_->getMinimizationObjValue();
     1362            value = value + 1.0e-7 * (1.0 + fabs(value));
     1363            value *= solver3->getObjSense();
     1364            model.setCutoff(value);
     1365            sprintf(generalPrint, "Unable to insert previous solution - using cutoff of %g",
     1366              value);
     1367            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
     1368              << generalPrint
     1369              << CoinMessageEol;
    14031370#ifdef CLP_INVESTIGATE
    1404                         printf("NOT added seren\n");
    1405                         solver3->writeMps("bad_seren");
    1406                         solver->writeMps("orig_seren");
    1407 #endif
    1408                     }
    1409                     delete solver3;
    1410                 }
    1411                 if (model_->searchStrategy() == 2) {
    1412                     model.setNumberStrong(5);
    1413                     model.setNumberBeforeTrust(5);
    1414                 }
    1415                 if (model.getNumCols()) {
    1416                     if (numberNodes >= 0) {
    1417                         setCutAndHeuristicOptions(model);
    1418                         // not too many iterations
    1419                         model.setMaximumNumberIterations(iterationMultiplier*(numberNodes + 10));
    1420                         // Not fast stuff
    1421                         model.setFastNodeDepth(-1);
    1422                         //model.solver()->writeMps("before");
    1423                     } else if (model.fastNodeDepth() >= 1000000) {
    1424                         // already set
    1425                         model.setFastNodeDepth(model.fastNodeDepth() - 1000000);
    1426                     }
    1427                     model.setWhenCuts(999998);
     1371            printf("NOT added seren\n");
     1372            solver3->writeMps("bad_seren");
     1373            solver->writeMps("orig_seren");
     1374#endif
     1375          }
     1376          delete solver3;
     1377        }
     1378        if (model_->searchStrategy() == 2) {
     1379          model.setNumberStrong(5);
     1380          model.setNumberBeforeTrust(5);
     1381        }
     1382        if (model.getNumCols()) {
     1383          if (numberNodes >= 0) {
     1384            setCutAndHeuristicOptions(model);
     1385            // not too many iterations
     1386            model.setMaximumNumberIterations(iterationMultiplier * (numberNodes + 10));
     1387            // Not fast stuff
     1388            model.setFastNodeDepth(-1);
     1389            //model.solver()->writeMps("before");
     1390          } else if (model.fastNodeDepth() >= 1000000) {
     1391            // already set
     1392            model.setFastNodeDepth(model.fastNodeDepth() - 1000000);
     1393          }
     1394          model.setWhenCuts(999998);
    14281395#define ALWAYS_DUAL
    14291396#ifdef ALWAYS_DUAL
    1430                     OsiSolverInterface * solverD = model.solver();
    1431                     bool takeHint;
    1432                     OsiHintStrength strength;
    1433                     solverD->getHintParam(OsiDoDualInResolve, takeHint, strength);
    1434                     solverD->setHintParam(OsiDoDualInResolve, true, OsiHintDo);
    1435 #endif
    1436                     model.passInEventHandler(model_->getEventHandler());
    1437                     // say model_ is sitting there
    1438                     int saveOptions = model_->specialOptions();
    1439                     model_->setSpecialOptions(saveOptions|1048576);
    1440                     // and switch off debugger
    1441                     model.setSpecialOptions(model.specialOptions()&(~1));
     1397          OsiSolverInterface *solverD = model.solver();
     1398          bool takeHint;
     1399          OsiHintStrength strength;
     1400          solverD->getHintParam(OsiDoDualInResolve, takeHint, strength);
     1401          solverD->setHintParam(OsiDoDualInResolve, true, OsiHintDo);
     1402#endif
     1403          model.passInEventHandler(model_->getEventHandler());
     1404          // say model_ is sitting there
     1405          int saveOptions = model_->specialOptions();
     1406          model_->setSpecialOptions(saveOptions | 1048576);
     1407          // and switch off debugger
     1408          model.setSpecialOptions(model.specialOptions() & (~1));
    14421409#if 0 //def COIN_HAS_CLP
    14431410                    OsiClpSolverInterface * clpSolver
     
    14471414#endif
    14481415#ifdef CONFLICT_CUTS
    1449                     if ((model_->moreSpecialOptions()&4194304)!=0)
    1450                       model.zapGlobalCuts();
     1416          if ((model_->moreSpecialOptions() & 4194304) != 0)
     1417            model.zapGlobalCuts();
    14511418#endif
    14521419#ifdef CGL_DEBUG
    1453                     if ((model_->specialOptions()&1) != 0) {
    1454                       const OsiRowCutDebugger *debugger = model.solver()->getRowCutDebugger() ;
    1455                       if (debugger) {
    1456                         printf("On optimal path DD\n");
    1457                       }
    1458                     }
    1459 #endif
    1460                     model.setPreProcess(&process);
    1461                     model.branchAndBound();
    1462                     model_->setHeuristicModel(NULL);
    1463                     model_->setSpecialOptions(saveOptions);
     1420          if ((model_->specialOptions() & 1) != 0) {
     1421            const OsiRowCutDebugger *debugger = model.solver()->getRowCutDebugger();
     1422            if (debugger) {
     1423              printf("On optimal path DD\n");
     1424            }
     1425          }
     1426#endif
     1427          model.setPreProcess(&process);
     1428          model.branchAndBound();
     1429          model_->setHeuristicModel(NULL);
     1430          model_->setSpecialOptions(saveOptions);
    14641431#ifdef ALWAYS_DUAL
    1465                     solverD = model.solver();
    1466                     solverD->setHintParam(OsiDoDualInResolve, takeHint, strength);
    1467 #endif
    1468                     numberNodesDone_ = model.getNodeCount();
     1432          solverD = model.solver();
     1433          solverD->setHintParam(OsiDoDualInResolve, takeHint, strength);
     1434#endif
     1435          numberNodesDone_ = model.getNodeCount();
    14691436#ifdef COIN_DEVELOP
    1470                     printf("sub branch %d nodes, %d iterations - max %d\n",
    1471                            model.getNodeCount(), model.getIterationCount(),
    1472                            100*(numberNodes + 10));
    1473 #endif
    1474                     if (numberNodes < 0) {
    1475                         model_->incrementIterationCount(model.getIterationCount());
    1476                         model_->incrementNodeCount(model.getNodeCount());
    1477                         // update best solution (in case ctrl-c)
    1478                         // !!! not a good idea - think a bit harder
    1479                         //model_->setMinimizationObjValue(model.getMinimizationObjValue());
    1480                         for (int iGenerator = 0; iGenerator < model.numberCutGenerators(); iGenerator++) {
    1481                             CbcCutGenerator * generator = model.cutGenerator(iGenerator);
    1482                             sprintf(generalPrint,
    1483                                     "%s was tried %d times and created %d cuts of which %d were active after adding rounds of cuts (%.3f seconds)",
    1484                                     generator->cutGeneratorName(),
    1485                                     generator->numberTimesEntered(),
    1486                                     generator->numberCutsInTotal() +
    1487                                     generator->numberColumnCuts(),
    1488                                     generator->numberCutsActive(),
    1489                                     generator->timeInCutGenerator());
    1490                             CglStored * stored = dynamic_cast<CglStored*>(generator->generator());
    1491                             if (stored && !generator->numberCutsInTotal())
    1492                                 continue;
     1437          printf("sub branch %d nodes, %d iterations - max %d\n",
     1438            model.getNodeCount(), model.getIterationCount(),
     1439            100 * (numberNodes + 10));
     1440#endif
     1441          if (numberNodes < 0) {
     1442            model_->incrementIterationCount(model.getIterationCount());
     1443            model_->incrementNodeCount(model.getNodeCount());
     1444            // update best solution (in case ctrl-c)
     1445            // !!! not a good idea - think a bit harder
     1446            //model_->setMinimizationObjValue(model.getMinimizationObjValue());
     1447            for (int iGenerator = 0; iGenerator < model.numberCutGenerators(); iGenerator++) {
     1448              CbcCutGenerator *generator = model.cutGenerator(iGenerator);
     1449              sprintf(generalPrint,
     1450                "%s was tried %d times and created %d cuts of which %d were active after adding rounds of cuts (%.3f seconds)",
     1451                generator->cutGeneratorName(),
     1452                generator->numberTimesEntered(),
     1453                generator->numberCutsInTotal() + generator->numberColumnCuts(),
     1454                generator->numberCutsActive(),
     1455                generator->timeInCutGenerator());
     1456              CglStored *stored = dynamic_cast<CglStored *>(generator->generator());
     1457              if (stored && !generator->numberCutsInTotal())
     1458                continue;
    14931459#ifndef CLP_INVESTIGATE
    1494                             CglImplication * implication = dynamic_cast<CglImplication*>(generator->generator());
    1495                             if (implication)
    1496                                 continue;
    1497 #endif
    1498                             model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
    1499                             << generalPrint
    1500                             << CoinMessageEol;
    1501                         }
     1460              CglImplication *implication = dynamic_cast<CglImplication *>(generator->generator());
     1461              if (implication)
     1462                continue;
     1463#endif
     1464              model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
     1465                << generalPrint
     1466                << CoinMessageEol;
     1467            }
     1468          }
     1469        } else {
     1470          // empty model
     1471          model.setMinimizationObjValue(model.solver()->getObjSense() * model.solver()->getObjValue());
     1472        }
     1473        if (logLevel > 1)
     1474          model_->messageHandler()->message(CBC_END_SUB, model_->messages())
     1475            << name
     1476            << CoinMessageEol;
     1477        if (model.getMinimizationObjValue() < CoinMin(cutoff, 1.0e30)) {
     1478          // solution
     1479          if (model.getNumCols())
     1480            returnCode = model.isProvenOptimal() ? 3 : 1;
     1481          else
     1482            returnCode = 3;
     1483            // post process
     1484#ifdef COIN_HAS_CLP
     1485          OsiClpSolverInterface *clpSolver = dynamic_cast<OsiClpSolverInterface *>(model.solver());
     1486          if (clpSolver) {
     1487            ClpSimplex *lpSolver = clpSolver->getModelPtr();
     1488            lpSolver->setSpecialOptions(lpSolver->specialOptions() | 0x01000000); // say is Cbc (and in branch and bound)
     1489          }
     1490#endif
     1491          //if (fractionSmall_ < 1000000.0)
     1492          process.postProcess(*model.solver());
     1493          if (solver->isProvenOptimal() && solver->getObjValue() * solver->getObjSense() < cutoff) {
     1494            // Solution now back in solver
     1495            int numberColumns = solver->getNumCols();
     1496            memcpy(newSolution, solver->getColSolution(),
     1497              numberColumns * sizeof(double));
     1498            newSolutionValue = model.getMinimizationObjValue();
     1499#ifdef COIN_HAS_CLP
     1500            if (clpSolver) {
     1501              if (clpSolver && clpSolver->numberSOS()) {
     1502                // SOS
     1503                int numberSOS = clpSolver->numberSOS();
     1504                const CoinSet *setInfo = clpSolver->setInfo();
     1505                int i;
     1506                for (i = 0; i < numberSOS; i++) {
     1507                  int type = setInfo[i].setType();
     1508                  int n = setInfo[i].numberEntries();
     1509                  const int *which = setInfo[i].which();
     1510                  int first = -1;
     1511                  int last = -1;
     1512                  for (int j = 0; j < n; j++) {
     1513                    int iColumn = which[j];
     1514                    if (fabs(newSolution[iColumn]) > 1.0e-7) {
     1515                      last = j;
     1516                      if (first < 0)
     1517                        first = j;
    15021518                    }
    1503                 } else {
    1504                     // empty model
    1505                     model.setMinimizationObjValue(model.solver()->getObjSense()*model.solver()->getObjValue());
     1519                  }
     1520                  assert(last - first < type);
     1521                  for (int j = 0; j < n; j++) {
     1522                    if (j < first || j > last) {
     1523                      int iColumn = which[j];
     1524                      // do I want to fix??
     1525                      solver->setColLower(iColumn, 0.0);
     1526                      solver->setColUpper(iColumn, 0.0);
     1527                      newSolution[iColumn] = 0.0;
     1528                    }
     1529                  }
    15061530                }
    1507                 if (logLevel > 1)
    1508                     model_->messageHandler()->message(CBC_END_SUB, model_->messages())
    1509                     << name
    1510                     << CoinMessageEol;
    1511                 if (model.getMinimizationObjValue() < CoinMin(cutoff, 1.0e30)) {
    1512                     // solution
    1513                     if (model.getNumCols())
    1514                         returnCode = model.isProvenOptimal() ? 3 : 1;
    1515                     else
    1516                         returnCode = 3;
    1517                     // post process
    1518 #ifdef COIN_HAS_CLP
    1519                     OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (model.solver());
    1520                     if (clpSolver) {
    1521                         ClpSimplex * lpSolver = clpSolver->getModelPtr();
    1522                         lpSolver->setSpecialOptions(lpSolver->specialOptions() | 0x01000000); // say is Cbc (and in branch and bound)
    1523                     }
    1524 #endif
    1525                     //if (fractionSmall_ < 1000000.0)
    1526                     process.postProcess(*model.solver());
    1527                     if (solver->isProvenOptimal() && solver->getObjValue()*solver->getObjSense() < cutoff) {
    1528                         // Solution now back in solver
    1529                         int numberColumns = solver->getNumCols();
    1530                         memcpy(newSolution, solver->getColSolution(),
    1531                                numberColumns*sizeof(double));
    1532                         newSolutionValue = model.getMinimizationObjValue();
    1533 #ifdef COIN_HAS_CLP
    1534                       if (clpSolver) {
    1535                         if (clpSolver && clpSolver->numberSOS()) {
    1536                           // SOS
    1537                           int numberSOS = clpSolver->numberSOS();
    1538                           const CoinSet * setInfo = clpSolver->setInfo();
    1539                           int i;
    1540                           for ( i = 0; i < numberSOS; i++) {
    1541                             int type = setInfo[i].setType();
    1542                             int n = setInfo[i].numberEntries();
    1543                             const int * which = setInfo[i].which();
    1544                             int first = -1;
    1545                             int last = -1;
    1546                             for (int j = 0; j < n; j++) {
    1547                               int iColumn = which[j];
    1548                               if (fabs(newSolution[iColumn]) > 1.0e-7) {
    1549                                 last = j;
    1550                                 if (first < 0)
    1551                                   first = j;
    1552                               }
    1553                             }
    1554                             assert (last - first < type);
    1555                             for (int j = 0; j < n; j++) {
    1556                               if (j < first || j > last) {
    1557                                 int iColumn = which[j];
    1558                                 // do I want to fix??
    1559                                 solver->setColLower(iColumn, 0.0);
    1560                                 solver->setColUpper(iColumn, 0.0);
    1561                                 newSolution[iColumn]=0.0;
    1562                               }
    1563                             }
    1564                           }
    1565                         }
    1566                       }
    1567 #endif
    1568                     } else {
    1569                         // odd - but no good
    1570                         returnCode = 0; // so will be infeasible
    1571                     }
    1572                 } else {
    1573                     // no good
    1574                     returnCode = model.isProvenInfeasible() ? 2 : 0; // so will be infeasible
    1575                 }
    1576                 int totalNumberIterations = model.getIterationCount() +
    1577                                             process.numberIterationsPre() +
    1578                                             process.numberIterationsPost();
    1579                 if (totalNumberIterations > 100*(numberNodes + 10)
    1580                     && fractionSmall_ < 1000000.0) {
    1581                     // only allow smaller problems
    1582                     fractionSmall = fractionSmall_;
    1583                     fractionSmall_ *= 0.9;
     1531              }
     1532            }
     1533#endif
     1534          } else {
     1535            // odd - but no good
     1536            returnCode = 0; // so will be infeasible
     1537          }
     1538        } else {
     1539          // no good
     1540          returnCode = model.isProvenInfeasible() ? 2 : 0; // so will be infeasible
     1541        }
     1542        int totalNumberIterations = model.getIterationCount() + process.numberIterationsPre() + process.numberIterationsPost();
     1543        if (totalNumberIterations > 100 * (numberNodes + 10)
     1544          && fractionSmall_ < 1000000.0) {
     1545          // only allow smaller problems
     1546          fractionSmall = fractionSmall_;
     1547          fractionSmall_ *= 0.9;
    15841548#ifdef CLP_INVESTIGATE
    1585                     printf("changing fractionSmall from %g to %g for %s as %d iterations\n",
    1586                            fractionSmall, fractionSmall_, name.c_str(), totalNumberIterations);
    1587 #endif
    1588                 }
    1589                 if (model.status() == 5)
    1590                    model_->sayEventHappened();
     1549          printf("changing fractionSmall from %g to %g for %s as %d iterations\n",
     1550            fractionSmall, fractionSmall_, name.c_str(), totalNumberIterations);
     1551#endif
     1552        }
     1553        if (model.status() == 5)
     1554          model_->sayEventHappened();
    15911555#ifdef COIN_DEVELOP
    1592                 if (model.isProvenInfeasible())
    1593                     status = 1;
    1594                 else if (model.isProvenOptimal())
    1595                     status = 2;
    1596 #endif
     1556        if (model.isProvenInfeasible())
     1557          status = 1;
     1558        else if (model.isProvenOptimal())
     1559          status = 2;
     1560#endif
     1561      }
     1562    }
     1563  } else {
     1564    returnCode = 2; // infeasible finished
     1565    if (logLevel > 1) {
     1566      printf("Infeasible on initial solve\n");
     1567    }
     1568  }
     1569  model_->setSpecialOptions(saveModelOptions);
     1570  model_->setLogLevel(logLevel);
     1571  if (returnCode == 1 || returnCode == 2) {
     1572    OsiSolverInterface *solverC = model_->continuousSolver();
     1573    if (false && solverC) {
     1574      const double *lower = solver->getColLower();
     1575      const double *upper = solver->getColUpper();
     1576      const double *lowerC = solverC->getColLower();
     1577      const double *upperC = solverC->getColUpper();
     1578      bool good = true;
     1579      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     1580        if (isHeuristicInteger(solverC, iColumn)) {
     1581          if (lower[iColumn] > lowerC[iColumn] && upper[iColumn] < upperC[iColumn]) {
     1582            good = false;
     1583            printf("CUT - can't add\n");
     1584            break;
     1585          }
     1586        }
     1587      }
     1588      if (good) {
     1589        double *cut = new double[numberColumns];
     1590        int *which = new int[numberColumns];
     1591        double rhs = -1.0;
     1592        int n = 0;
     1593        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     1594          if (isHeuristicInteger(solverC, iColumn)) {
     1595            if (lower[iColumn] == upperC[iColumn]) {
     1596              rhs += lower[iColumn];
     1597              cut[n] = 1.0;
     1598              which[n++] = iColumn;
     1599            } else if (upper[iColumn] == lowerC[iColumn]) {
     1600              rhs -= upper[iColumn];
     1601              cut[n] = -1.0;
     1602              which[n++] = iColumn;
    15971603            }
    1598         }
     1604          }
     1605        }
     1606        printf("CUT has %d entries\n", n);
     1607        OsiRowCut newCut;
     1608        newCut.setLb(-COIN_DBL_MAX);
     1609        newCut.setUb(rhs);
     1610        newCut.setRow(n, which, cut, false);
     1611        model_->makeGlobalCut(newCut);
     1612        delete[] cut;
     1613        delete[] which;
     1614      }
     1615    }
     1616#ifdef COIN_DEVELOP
     1617    if (status == 1)
     1618      printf("heuristic could add cut because infeasible (%s)\n", heuristicName_.c_str());
     1619    else if (status == 2)
     1620      printf("heuristic could add cut because optimal (%s)\n", heuristicName_.c_str());
     1621#endif
     1622  }
     1623  if (reset) {
     1624    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     1625      if (reset[iColumn])
     1626        solver->setColLower(iColumn, 0.0);
     1627    }
     1628    delete[] reset;
     1629  }
     1630#ifdef HISTORY_STATISTICS
     1631  getHistoryStatistics_ = true;
     1632#endif
     1633  solver->setHintParam(OsiDoReducePrint, takeHint, strength);
     1634  return returnCode;
     1635}
     1636// Set input solution
     1637void CbcHeuristic::setInputSolution(const double *solution, double objValue)
     1638{
     1639  delete[] inputSolution_;
     1640  inputSolution_ = NULL;
     1641  if (model_ && solution) {
     1642    int numberColumns = model_->getNumCols();
     1643    inputSolution_ = new double[numberColumns + 1];
     1644    memcpy(inputSolution_, solution, numberColumns * sizeof(double));
     1645    inputSolution_[numberColumns] = objValue;
     1646  }
     1647}
     1648
     1649//##############################################################################
     1650
     1651inline int compare3BranchingObjects(const CbcBranchingObject *br0,
     1652  const CbcBranchingObject *br1)
     1653{
     1654  const int t0 = br0->type();
     1655  const int t1 = br1->type();
     1656  if (t0 < t1) {
     1657    return -1;
     1658  }
     1659  if (t0 > t1) {
     1660    return 1;
     1661  }
     1662  return br0->compareOriginalObject(br1);
     1663}
     1664
     1665//==============================================================================
     1666
     1667inline bool compareBranchingObjects(const CbcBranchingObject *br0,
     1668  const CbcBranchingObject *br1)
     1669{
     1670  return compare3BranchingObjects(br0, br1) < 0;
     1671}
     1672
     1673//==============================================================================
     1674
     1675void CbcHeuristicNode::gutsOfConstructor(CbcModel &model)
     1676{
     1677  //  CbcHeurDebugNodes(&model);
     1678  CbcNode *node = model.currentNode();
     1679  brObj_ = new CbcBranchingObject *[node->depth()];
     1680  CbcNodeInfo *nodeInfo = node->nodeInfo();
     1681  int cnt = 0;
     1682  while (nodeInfo->parentBranch() != NULL) {
     1683    const OsiBranchingObject *br = nodeInfo->parentBranch();
     1684    const CbcBranchingObject *cbcbr = dynamic_cast<const CbcBranchingObject *>(br);
     1685    if (!cbcbr) {
     1686      throw CoinError("CbcHeuristicNode can be used only with CbcBranchingObjects.\n",
     1687        "gutsOfConstructor",
     1688        "CbcHeuristicNode",
     1689        __FILE__, __LINE__);
     1690    }
     1691    brObj_[cnt] = cbcbr->clone();
     1692    brObj_[cnt]->previousBranch();
     1693    ++cnt;
     1694    nodeInfo = nodeInfo->parent();
     1695  }
     1696  std::sort(brObj_, brObj_ + cnt, compareBranchingObjects);
     1697  if (cnt <= 1) {
     1698    numObjects_ = cnt;
     1699  } else {
     1700    numObjects_ = 0;
     1701    CbcBranchingObject *br = NULL; // What should this be?
     1702    for (int i = 1; i < cnt; ++i) {
     1703      if (compare3BranchingObjects(brObj_[numObjects_], brObj_[i]) == 0) {
     1704        int comp = brObj_[numObjects_]->compareBranchingObject(brObj_[i], br != 0);
     1705        switch (comp) {
     1706        case CbcRangeSame: // the same range
     1707        case CbcRangeDisjoint: // disjoint decisions
     1708          // should not happen! we are on a chain!
     1709          abort();
     1710        case CbcRangeSubset: // brObj_[numObjects_] is a subset of brObj_[i]
     1711          delete brObj_[i];
     1712          break;
     1713        case CbcRangeSuperset: // brObj_[i] is a subset of brObj_[numObjects_]
     1714          delete brObj_[numObjects_];
     1715          brObj_[numObjects_] = brObj_[i];
     1716          break;
     1717        case CbcRangeOverlap: // overlap
     1718          delete brObj_[i];
     1719          delete brObj_[numObjects_];
     1720          brObj_[numObjects_] = br;
     1721          break;
     1722        }
     1723        continue;
     1724      } else {
     1725        brObj_[++numObjects_] = brObj_[i];
     1726      }
     1727    }
     1728    ++numObjects_;
     1729  }
     1730}
     1731
     1732//==============================================================================
     1733
     1734CbcHeuristicNode::CbcHeuristicNode(CbcModel &model)
     1735{
     1736  gutsOfConstructor(model);
     1737}
     1738
     1739//==============================================================================
     1740
     1741double
     1742CbcHeuristicNode::distance(const CbcHeuristicNode *node) const
     1743{
     1744
     1745  const double disjointWeight = 1;
     1746  const double overlapWeight = 0.4;
     1747  const double subsetWeight = 0.2;
     1748  int countDisjointWeight = 0;
     1749  int countOverlapWeight = 0;
     1750  int countSubsetWeight = 0;
     1751  int i = 0;
     1752  int j = 0;
     1753  double dist = 0.0;
     1754#ifdef PRINT_DEBUG
     1755  printf(" numObjects_ = %i, node->numObjects_ = %i\n",
     1756    numObjects_, node->numObjects_);
     1757#endif
     1758  while (i < numObjects_ && j < node->numObjects_) {
     1759    CbcBranchingObject *br0 = brObj_[i];
     1760    const CbcBranchingObject *br1 = node->brObj_[j];
     1761#ifdef PRINT_DEBUG
     1762    const CbcIntegerBranchingObject *brPrint0 = dynamic_cast<const CbcIntegerBranchingObject *>(br0);
     1763    const double *downBounds = brPrint0->downBounds();
     1764    const double *upBounds = brPrint0->upBounds();
     1765    int variable = brPrint0->variable();
     1766    int way = brPrint0->way();
     1767    printf("   br0: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
     1768      variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
     1769      static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
     1770    const CbcIntegerBranchingObject *brPrint1 = dynamic_cast<const CbcIntegerBranchingObject *>(br1);
     1771    downBounds = brPrint1->downBounds();
     1772    upBounds = brPrint1->upBounds();
     1773    variable = brPrint1->variable();
     1774    way = brPrint1->way();
     1775    printf("   br1: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
     1776      variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
     1777      static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
     1778#endif
     1779    const int brComp = compare3BranchingObjects(br0, br1);
     1780    if (brComp < 0) {
     1781      dist += subsetWeight;
     1782      countSubsetWeight++;
     1783      ++i;
     1784    } else if (brComp > 0) {
     1785      dist += subsetWeight;
     1786      countSubsetWeight++;
     1787      ++j;
    15991788    } else {
    1600         returnCode = 2; // infeasible finished
    1601         if (logLevel > 1){
    1602            printf("Infeasible on initial solve\n");
    1603         }
    1604     }
    1605     model_->setSpecialOptions(saveModelOptions);
    1606     model_->setLogLevel(logLevel);
    1607     if (returnCode == 1 || returnCode == 2) {
    1608         OsiSolverInterface * solverC = model_->continuousSolver();
    1609         if (false && solverC) {
    1610             const double * lower = solver->getColLower();
    1611             const double * upper = solver->getColUpper();
    1612             const double * lowerC = solverC->getColLower();
    1613             const double * upperC = solverC->getColUpper();
    1614             bool good = true;
    1615             for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    1616                 if (isHeuristicInteger(solverC,iColumn)) {
    1617                     if (lower[iColumn] > lowerC[iColumn] &&
    1618                             upper[iColumn] < upperC[iColumn]) {
    1619                         good = false;
    1620                         printf("CUT - can't add\n");
    1621                         break;
    1622                     }
    1623                 }
    1624             }
    1625             if (good) {
    1626                 double * cut = new double [numberColumns];
    1627                 int * which = new int [numberColumns];
    1628                 double rhs = -1.0;
    1629                 int n = 0;
    1630                 for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    1631                     if (isHeuristicInteger(solverC,iColumn)) {
    1632                         if (lower[iColumn] == upperC[iColumn]) {
    1633                             rhs += lower[iColumn];
    1634                             cut[n] = 1.0;
    1635                             which[n++] = iColumn;
    1636                         } else if (upper[iColumn] == lowerC[iColumn]) {
    1637                             rhs -= upper[iColumn];
    1638                             cut[n] = -1.0;
    1639                             which[n++] = iColumn;
    1640                         }
    1641                     }
    1642                 }
    1643                 printf("CUT has %d entries\n", n);
    1644                 OsiRowCut newCut;
    1645                 newCut.setLb(-COIN_DBL_MAX);
    1646                 newCut.setUb(rhs);
    1647                 newCut.setRow(n, which, cut, false);
    1648                 model_->makeGlobalCut(newCut);
    1649                 delete [] cut;
    1650                 delete [] which;
    1651             }
    1652         }
    1653 #ifdef COIN_DEVELOP
    1654         if (status == 1)
    1655             printf("heuristic could add cut because infeasible (%s)\n", heuristicName_.c_str());
    1656         else if (status == 2)
    1657             printf("heuristic could add cut because optimal (%s)\n", heuristicName_.c_str());
    1658 #endif
    1659     }
    1660     if (reset) {
    1661         for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    1662             if (reset[iColumn])
    1663                 solver->setColLower(iColumn, 0.0);
    1664         }
    1665         delete [] reset;
    1666     }
    1667 #ifdef HISTORY_STATISTICS
    1668     getHistoryStatistics_ = true;
    1669 #endif
    1670     solver->setHintParam(OsiDoReducePrint, takeHint, strength);
    1671     return returnCode;
    1672 }
    1673 // Set input solution
    1674 void
    1675 CbcHeuristic::setInputSolution(const double * solution, double objValue)
    1676 {
    1677     delete [] inputSolution_;
    1678     inputSolution_ = NULL;
    1679     if (model_ && solution) {
    1680         int numberColumns = model_->getNumCols();
    1681         inputSolution_ = new double [numberColumns+1];
    1682         memcpy(inputSolution_, solution, numberColumns*sizeof(double));
    1683         inputSolution_[numberColumns] = objValue;
    1684     }
    1685 }
    1686 
    1687 //##############################################################################
    1688 
    1689 inline int compare3BranchingObjects(const CbcBranchingObject* br0,
    1690                                     const CbcBranchingObject* br1)
    1691 {
    1692     const int t0 = br0->type();
    1693     const int t1 = br1->type();
    1694     if (t0 < t1) {
    1695         return -1;
    1696     }
    1697     if (t0 > t1) {
    1698         return 1;
    1699     }
    1700     return br0->compareOriginalObject(br1);
     1789      const int comp = br0->compareBranchingObject(br1, false);
     1790      switch (comp) {
     1791      case CbcRangeSame:
     1792        // do nothing
     1793        break;
     1794      case CbcRangeDisjoint: // disjoint decisions
     1795        dist += disjointWeight;
     1796        countDisjointWeight++;
     1797        break;
     1798      case CbcRangeSubset: // subset one way or another
     1799      case CbcRangeSuperset:
     1800        dist += subsetWeight;
     1801        countSubsetWeight++;
     1802        break;
     1803      case CbcRangeOverlap: // overlap
     1804        dist += overlapWeight;
     1805        countOverlapWeight++;
     1806        break;
     1807      }
     1808      ++i;
     1809      ++j;
     1810    }
     1811  }
     1812  dist += subsetWeight * (numObjects_ - i + node->numObjects_ - j);
     1813  countSubsetWeight += (numObjects_ - i + node->numObjects_ - j);
     1814  COIN_DETAIL_PRINT(printf("subset = %i, overlap = %i, disjoint = %i\n", countSubsetWeight,
     1815    countOverlapWeight, countDisjointWeight));
     1816  return dist;
    17011817}
    17021818
    17031819//==============================================================================
    17041820
    1705 inline bool compareBranchingObjects(const CbcBranchingObject* br0,
    1706                                     const CbcBranchingObject* br1)
    1707 {
    1708     return compare3BranchingObjects(br0, br1) < 0;
     1821CbcHeuristicNode::~CbcHeuristicNode()
     1822{
     1823  for (int i = 0; i < numObjects_; ++i) {
     1824    delete brObj_[i];
     1825  }
     1826  delete[] brObj_;
    17091827}
    17101828
    17111829//==============================================================================
    17121830
    1713 void
    1714 CbcHeuristicNode::gutsOfConstructor(CbcModel& model)
    1715 {
    1716     //  CbcHeurDebugNodes(&model);
    1717     CbcNode* node = model.currentNode();
    1718     brObj_ = new CbcBranchingObject*[node->depth()];
    1719     CbcNodeInfo* nodeInfo = node->nodeInfo();
    1720     int cnt = 0;
    1721     while (nodeInfo->parentBranch() != NULL) {
    1722         const OsiBranchingObject* br = nodeInfo->parentBranch();
    1723         const CbcBranchingObject* cbcbr = dynamic_cast<const CbcBranchingObject*>(br);
    1724         if (! cbcbr) {
    1725             throw CoinError("CbcHeuristicNode can be used only with CbcBranchingObjects.\n",
    1726                             "gutsOfConstructor",
    1727                             "CbcHeuristicNode",
    1728                             __FILE__, __LINE__);
    1729         }
    1730         brObj_[cnt] = cbcbr->clone();
    1731         brObj_[cnt]->previousBranch();
    1732         ++cnt;
    1733         nodeInfo = nodeInfo->parent();
    1734     }
    1735     std::sort(brObj_, brObj_ + cnt, compareBranchingObjects);
    1736     if (cnt <= 1) {
    1737         numObjects_ = cnt;
     1831double
     1832CbcHeuristicNode::minDistance(const CbcHeuristicNodeList &nodeList) const
     1833{
     1834  double minDist = COIN_DBL_MAX;
     1835  for (int i = nodeList.size() - 1; i >= 0; --i) {
     1836    minDist = CoinMin(minDist, distance(nodeList.node(i)));
     1837  }
     1838  return minDist;
     1839}
     1840
     1841//==============================================================================
     1842
     1843bool CbcHeuristicNode::minDistanceIsSmall(const CbcHeuristicNodeList &nodeList,
     1844  const double threshold) const
     1845{
     1846  for (int i = nodeList.size() - 1; i >= 0; --i) {
     1847    if (distance(nodeList.node(i)) >= threshold) {
     1848      continue;
    17381849    } else {
    1739         numObjects_ = 0;
    1740         CbcBranchingObject* br = NULL; // What should this be?
    1741         for (int i = 1; i < cnt; ++i) {
    1742             if (compare3BranchingObjects(brObj_[numObjects_], brObj_[i]) == 0) {
    1743                 int comp = brObj_[numObjects_]->compareBranchingObject(brObj_[i], br != 0);
    1744                 switch (comp) {
    1745                 case CbcRangeSame: // the same range
    1746                 case CbcRangeDisjoint: // disjoint decisions
    1747                     // should not happen! we are on a chain!
    1748                     abort();
    1749                 case CbcRangeSubset: // brObj_[numObjects_] is a subset of brObj_[i]
    1750                     delete brObj_[i];
    1751                     break;
    1752                 case CbcRangeSuperset: // brObj_[i] is a subset of brObj_[numObjects_]
    1753                     delete brObj_[numObjects_];
    1754                     brObj_[numObjects_] = brObj_[i];
    1755                     break;
    1756                 case CbcRangeOverlap: // overlap
    1757                     delete brObj_[i];
    1758                     delete brObj_[numObjects_];
    1759                     brObj_[numObjects_] = br;
    1760                     break;
    1761                 }
    1762                 continue;
    1763             } else {
    1764                 brObj_[++numObjects_] = brObj_[i];
    1765             }
    1766         }
    1767         ++numObjects_;
    1768     }
     1850      return true;
     1851    }
     1852  }
     1853  return false;
    17691854}
    17701855
    17711856//==============================================================================
    17721857
    1773 CbcHeuristicNode::CbcHeuristicNode(CbcModel& model)
    1774 {
    1775     gutsOfConstructor(model);
    1776 }
    1777 
    1778 //==============================================================================
    1779 
    17801858double
    1781 CbcHeuristicNode::distance(const CbcHeuristicNode* node) const
    1782 {
    1783 
    1784     const double disjointWeight = 1;
    1785     const double overlapWeight = 0.4;
    1786     const double subsetWeight = 0.2;
    1787     int countDisjointWeight = 0;
    1788     int countOverlapWeight = 0;
    1789     int countSubsetWeight = 0;
    1790     int i = 0;
    1791     int j = 0;
    1792     double dist = 0.0;
    1793 #ifdef PRINT_DEBUG
    1794     printf(" numObjects_ = %i, node->numObjects_ = %i\n",
    1795            numObjects_, node->numObjects_);
    1796 #endif
    1797     while ( i < numObjects_ && j < node->numObjects_) {
    1798         CbcBranchingObject* br0 = brObj_[i];
    1799         const CbcBranchingObject* br1 = node->brObj_[j];
    1800 #ifdef PRINT_DEBUG
    1801         const CbcIntegerBranchingObject* brPrint0 =
    1802             dynamic_cast<const CbcIntegerBranchingObject*>(br0);
    1803         const double* downBounds = brPrint0->downBounds();
    1804         const double* upBounds = brPrint0->upBounds();
    1805         int variable = brPrint0->variable();
    1806         int way = brPrint0->way();
    1807         printf("   br0: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
    1808                variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
    1809                static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
    1810         const CbcIntegerBranchingObject* brPrint1 =
    1811             dynamic_cast<const CbcIntegerBranchingObject*>(br1);
    1812         downBounds = brPrint1->downBounds();
    1813         upBounds = brPrint1->upBounds();
    1814         variable = brPrint1->variable();
    1815         way = brPrint1->way();
    1816         printf("   br1: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
    1817                variable, static_cast<int>(downBounds[0]), static_cast<int>(downBounds[1]),
    1818                static_cast<int>(upBounds[0]), static_cast<int>(upBounds[1]), way);
    1819 #endif
    1820         const int brComp = compare3BranchingObjects(br0, br1);
    1821         if (brComp < 0) {
    1822             dist += subsetWeight;
    1823             countSubsetWeight++;
    1824             ++i;
    1825         } else if (brComp > 0) {
    1826             dist += subsetWeight;
    1827             countSubsetWeight++;
    1828             ++j;
    1829         } else {
    1830             const int comp = br0->compareBranchingObject(br1, false);
    1831             switch (comp) {
    1832             case CbcRangeSame:
    1833                 // do nothing
    1834                 break;
    1835             case CbcRangeDisjoint: // disjoint decisions
    1836                 dist += disjointWeight;
    1837                 countDisjointWeight++;
    1838                 break;
    1839             case CbcRangeSubset: // subset one way or another
    1840             case CbcRangeSuperset:
    1841                 dist += subsetWeight;
    1842                 countSubsetWeight++;
    1843                 break;
    1844             case CbcRangeOverlap: // overlap
    1845                 dist += overlapWeight;
    1846                 countOverlapWeight++;
    1847                 break;
    1848             }
    1849             ++i;
    1850             ++j;
    1851         }
    1852     }
    1853     dist += subsetWeight * (numObjects_ - i + node->numObjects_ - j);
    1854     countSubsetWeight += (numObjects_ - i + node->numObjects_ - j);
    1855     COIN_DETAIL_PRINT(printf("subset = %i, overlap = %i, disjoint = %i\n", countSubsetWeight,
    1856                              countOverlapWeight, countDisjointWeight));
    1857     return dist;
    1858 }
    1859 
    1860 //==============================================================================
    1861 
    1862 CbcHeuristicNode::~CbcHeuristicNode()
    1863 {
    1864     for (int i = 0; i < numObjects_; ++i) {
    1865         delete brObj_[i];
    1866     }
    1867     delete [] brObj_;
    1868 }
    1869 
    1870 //==============================================================================
    1871 
    1872 double
    1873 CbcHeuristicNode::minDistance(const CbcHeuristicNodeList& nodeList) const
    1874 {
    1875     double minDist = COIN_DBL_MAX;
    1876     for (int i = nodeList.size() - 1; i >= 0; --i) {
    1877         minDist = CoinMin(minDist, distance(nodeList.node(i)));
    1878     }
    1879     return minDist;
    1880 }
    1881 
    1882 //==============================================================================
    1883 
    1884 bool
    1885 CbcHeuristicNode::minDistanceIsSmall(const CbcHeuristicNodeList& nodeList,
    1886                                      const double threshold) const
    1887 {
    1888     for (int i = nodeList.size() - 1; i >= 0; --i) {
    1889         if (distance(nodeList.node(i)) >= threshold) {
    1890             continue;
    1891         } else {
    1892             return true;
    1893         }
    1894     }
    1895     return false;
    1896 }
    1897 
    1898 //==============================================================================
    1899 
    1900 double
    1901 CbcHeuristicNode::avgDistance(const CbcHeuristicNodeList& nodeList) const
    1902 {
    1903     if (nodeList.size() == 0) {
    1904         return COIN_DBL_MAX;
    1905     }
    1906     double sumDist = 0;
    1907     for (int i = nodeList.size() - 1; i >= 0; --i) {
    1908         sumDist += distance(nodeList.node(i));
    1909     }
    1910     return sumDist / nodeList.size();
     1859CbcHeuristicNode::avgDistance(const CbcHeuristicNodeList &nodeList) const
     1860{
     1861  if (nodeList.size() == 0) {
     1862    return COIN_DBL_MAX;
     1863  }
     1864  double sumDist = 0;
     1865  for (int i = nodeList.size() - 1; i >= 0; --i) {
     1866    sumDist += distance(nodeList.node(i));
     1867  }
     1868  return sumDist / nodeList.size();
    19111869}
    19121870
     
    19151873// Default Constructor
    19161874CbcRounding::CbcRounding()
    1917         : CbcHeuristic()
    1918 {
    1919     // matrix and row copy will automatically be empty
    1920     seed_ = 7654321;
    1921     down_ = NULL;
    1922     up_ = NULL;
    1923     equal_ = NULL;
    1924     //whereFrom_ |= 16*(1+256); // allow more often
     1875  : CbcHeuristic()
     1876{
     1877  // matrix and row copy will automatically be empty
     1878  seed_ = 7654321;
     1879  down_ = NULL;
     1880  up_ = NULL;
     1881  equal_ = NULL;
     1882  //whereFrom_ |= 16*(1+256); // allow more often
    19251883}
    19261884
    19271885// Constructor from model
    1928 CbcRounding::CbcRounding(CbcModel & model)
    1929         : CbcHeuristic(model)
    1930 {
    1931     // Get a copy of original matrix (and by row for rounding);
    1932     assert(model.solver());
    1933     if (model.solver()->getNumRows()) {
    1934         matrix_ = *model.solver()->getMatrixByCol();
    1935         matrixByRow_ = *model.solver()->getMatrixByRow();
    1936         validate();
    1937     }
    1938     down_ = NULL;
    1939     up_ = NULL;
    1940     equal_ = NULL;
    1941     seed_ = 7654321;
    1942     //whereFrom_ |= 16*(1+256); // allow more often
     1886CbcRounding::CbcRounding(CbcModel &model)
     1887  : CbcHeuristic(model)
     1888{
     1889  // Get a copy of original matrix (and by row for rounding);
     1890  assert(model.solver());
     1891  if (model.solver()->getNumRows()) {
     1892    matrix_ = *model.solver()->getMatrixByCol();
     1893    matrixByRow_ = *model.solver()->getMatrixByRow();
     1894    validate();
     1895  }
     1896  down_ = NULL;
     1897  up_ = NULL;
     1898  equal_ = NULL;
     1899  seed_ = 7654321;
     1900  //whereFrom_ |= 16*(1+256); // allow more often
    19431901}
    19441902
    19451903// Destructor
    1946 CbcRounding::~CbcRounding ()
    1947 {
    1948     delete [] down_;
    1949     delete [] up_;
    1950     delete [] equal_;
     1904CbcRounding::~CbcRounding()
     1905{
     1906  delete[] down_;
     1907  delete[] up_;
     1908  delete[] equal_;
    19511909}
    19521910
     
    19551913CbcRounding::clone() const
    19561914{
    1957     return new CbcRounding(*this);
     1915  return new CbcRounding(*this);
    19581916}
    19591917// Create C++ lines to get to current state
    1960 void
    1961 CbcRounding::generateCpp( FILE * fp)
    1962 {
    1963     CbcRounding other;
    1964     fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n");
    1965     fprintf(fp, "3  CbcRounding rounding(*cbcModel);\n");
    1966     CbcHeuristic::generateCpp(fp, "rounding");
    1967     if (seed_ != other.seed_)
    1968         fprintf(fp, "3  rounding.setSeed(%d);\n", seed_);
    1969     else
    1970         fprintf(fp, "4  rounding.setSeed(%d);\n", seed_);
    1971     fprintf(fp, "3  cbcModel->addHeuristic(&rounding);\n");
     1918void CbcRounding::generateCpp(FILE *fp)
     1919{
     1920  CbcRounding other;
     1921  fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n");
     1922  fprintf(fp, "3  CbcRounding rounding(*cbcModel);\n");
     1923  CbcHeuristic::generateCpp(fp, "rounding");
     1924  if (seed_ != other.seed_)
     1925    fprintf(fp, "3  rounding.setSeed(%d);\n", seed_);
     1926  else
     1927    fprintf(fp, "4  rounding.setSeed(%d);\n", seed_);
     1928  fprintf(fp, "3  cbcModel->addHeuristic(&rounding);\n");
    19721929}
    19731930//#define NEW_ROUNDING
    19741931// Copy constructor
    1975 CbcRounding::CbcRounding(const CbcRounding & rhs)
    1976         :
    1977         CbcHeuristic(rhs),
    1978         matrix_(rhs.matrix_),
    1979         matrixByRow_(rhs.matrixByRow_),
    1980         seed_(rhs.seed_)
     1932CbcRounding::CbcRounding(const CbcRounding &rhs)
     1933  : CbcHeuristic(rhs)
     1934  , matrix_(rhs.matrix_)
     1935  , matrixByRow_(rhs.matrixByRow_)
     1936  , seed_(rhs.seed_)
    19811937{
    19821938#ifdef NEW_ROUNDING
     1939  int numberColumns = matrix_.getNumCols();
     1940  down_ = CoinCopyOfArray(rhs.down_, numberColumns);
     1941  up_ = CoinCopyOfArray(rhs.up_, numberColumns);
     1942  equal_ = CoinCopyOfArray(rhs.equal_, numberColumns);
     1943#else
     1944  down_ = NULL;
     1945  up_ = NULL;
     1946  equal_ = NULL;
     1947#endif
     1948}
     1949
     1950// Assignment operator
     1951CbcRounding &
     1952CbcRounding::operator=(const CbcRounding &rhs)
     1953{
     1954  if (this != &rhs) {
     1955    CbcHeuristic::operator=(rhs);
     1956    matrix_ = rhs.matrix_;
     1957    matrixByRow_ = rhs.matrixByRow_;
     1958#ifdef NEW_ROUNDING
     1959    delete[] down_;
     1960    delete[] up_;
     1961    delete[] equal_;
    19831962    int numberColumns = matrix_.getNumCols();
    19841963    down_ = CoinCopyOfArray(rhs.down_, numberColumns);
     
    19901969    equal_ = NULL;
    19911970#endif
    1992 }
    1993 
    1994 // Assignment operator
    1995 CbcRounding &
    1996 CbcRounding::operator=( const CbcRounding & rhs)
    1997 {
    1998     if (this != &rhs) {
    1999         CbcHeuristic::operator=(rhs);
    2000         matrix_ = rhs.matrix_;
    2001         matrixByRow_ = rhs.matrixByRow_;
    2002 #ifdef NEW_ROUNDING
    2003         delete [] down_;
    2004         delete [] up_;
    2005         delete [] equal_;
    2006         int numberColumns = matrix_.getNumCols();
    2007         down_ = CoinCopyOfArray(rhs.down_, numberColumns);
    2008         up_ = CoinCopyOfArray(rhs.up_, numberColumns);
    2009         equal_ = CoinCopyOfArray(rhs.equal_, numberColumns);
    2010 #else
    2011         down_ = NULL;
    2012         up_ = NULL;
    2013         equal_ = NULL;
    2014 #endif
    2015         seed_ = rhs.seed_;
    2016     }
    2017     return *this;
     1971    seed_ = rhs.seed_;
     1972  }
     1973  return *this;
    20181974}
    20191975
    20201976// Resets stuff if model changes
    2021 void
    2022 CbcRounding::resetModel(CbcModel * model)
    2023 {
    2024     model_ = model;
    2025     // Get a copy of original matrix (and by row for rounding);
    2026     assert(model_->solver());
    2027     matrix_ = *model_->solver()->getMatrixByCol();
    2028     matrixByRow_ = *model_->solver()->getMatrixByRow();
    2029     validate();
     1977void CbcRounding::resetModel(CbcModel *model)
     1978{
     1979  model_ = model;
     1980  // Get a copy of original matrix (and by row for rounding);
     1981  assert(model_->solver());
     1982  matrix_ = *model_->solver()->getMatrixByCol();
     1983  matrixByRow_ = *model_->solver()->getMatrixByRow();
     1984  validate();
    20301985}
    20311986/* Check whether the heuristic should run at all
     
    20371992   8 added if previous heuristic in loop found solution
    20381993*/
    2039 bool
    2040 CbcRounding::shouldHeurRun(int whereFrom)
    2041 {
    2042   if (whereFrom!=4) {
     1994bool CbcRounding::shouldHeurRun(int whereFrom)
     1995{
     1996  if (whereFrom != 4) {
    20431997    return CbcHeuristic::shouldHeurRun(whereFrom);
    20441998  } else {
     
    20532007// Fix values if asked for
    20542008// Returns 1 if solution, 0 if not
    2055 int
    2056 CbcRounding::solution(double & solutionValue,
    2057                       double * betterSolution)
    2058 {
    2059 
    2060     numCouldRun_++;
    2061     // See if to do
    2062     if (!when() || (when() % 10 == 1 && model_->phase() != 1) ||
    2063             (when() % 10 == 2 && (model_->phase() != 2 && model_->phase() != 3)))
    2064         return 0; // switched off
    2065     numRuns_++;
     2009int CbcRounding::solution(double &solutionValue,
     2010  double *betterSolution)
     2011{
     2012
     2013  numCouldRun_++;
     2014  // See if to do
     2015  if (!when() || (when() % 10 == 1 && model_->phase() != 1) || (when() % 10 == 2 && (model_->phase() != 2 && model_->phase() != 3)))
     2016    return 0; // switched off
     2017  numRuns_++;
    20662018#ifdef HEURISTIC_INFORM
    2067     printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
    2068            heuristicName(),numRuns_,numCouldRun_,when_);
    2069 #endif
    2070     OsiSolverInterface * solver = model_->solver();
    2071     double direction = solver->getObjSense();
    2072     double newSolutionValue = direction * solver->getObjValue();
    2073     return solution(solutionValue, betterSolution, newSolutionValue);
     2019  printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
     2020    heuristicName(), numRuns_, numCouldRun_, when_);
     2021#endif
     2022  OsiSolverInterface *solver = model_->solver();
     2023  double direction = solver->getObjSense();
     2024  double newSolutionValue = direction * solver->getObjValue();
     2025  return solution(solutionValue, betterSolution, newSolutionValue);
    20742026}
    20752027// See if rounding will give solution
     
    20792031// Fix values if asked for
    20802032// Returns 1 if solution, 0 if not
    2081 int
    2082 CbcRounding::solution(double & solutionValue,
    2083                       double * betterSolution,
    2084                       double newSolutionValue)
    2085 {
    2086 
    2087     // See if to do
    2088     if (!when() || (when() % 10 == 1 && model_->phase() != 1) ||
    2089             (when() % 10 == 2 && (model_->phase() != 2 && model_->phase() != 3)))
    2090         return 0; // switched off
    2091     OsiSolverInterface * solver = model_->solver();
    2092     const double * lower = solver->getColLower();
    2093     const double * upper = solver->getColUpper();
    2094     const double * rowLower = solver->getRowLower();
    2095     const double * rowUpper = solver->getRowUpper();
    2096     const double * solution = solver->getColSolution();
    2097     const double * objective = solver->getObjCoefficients();
    2098     double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
    2099     double primalTolerance;
    2100     solver->getDblParam(OsiPrimalTolerance, primalTolerance);
    2101     double useTolerance = primalTolerance;
    2102 
    2103     int numberRows = matrix_.getNumRows();
    2104     assert (numberRows <= solver->getNumRows());
    2105     if (numberRows == 0){
    2106        return 0;
    2107     }
    2108     int numberIntegers = model_->numberIntegers();
    2109     const int * integerVariable = model_->integerVariable();
    2110     int i;
    2111     double direction = solver->getObjSense();
    2112     //double newSolutionValue = direction*solver->getObjValue();
    2113     int returnCode = 0;
    2114     // Column copy
    2115     const double * element = matrix_.getElements();
    2116     const int * row = matrix_.getIndices();
    2117     const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    2118     const int * columnLength = matrix_.getVectorLengths();
    2119     // Row copy
    2120     const double * elementByRow = matrixByRow_.getElements();
    2121     const int * column = matrixByRow_.getIndices();
    2122     const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
    2123     const int * rowLength = matrixByRow_.getVectorLengths();
    2124 
    2125     // Get solution array for heuristic solution
    2126     int numberColumns = solver->getNumCols();
    2127     double * newSolution = new double [numberColumns];
    2128     memcpy(newSolution, solution, numberColumns*sizeof(double));
    2129 
    2130     double * rowActivity = new double[numberRows];
    2131     memset(rowActivity, 0, numberRows*sizeof(double));
    2132     for (i = 0; i < numberColumns; i++) {
     2033int CbcRounding::solution(double &solutionValue,
     2034  double *betterSolution,
     2035  double newSolutionValue)
     2036{
     2037
     2038  // See if to do
     2039  if (!when() || (when() % 10 == 1 && model_->phase() != 1) || (when() % 10 == 2 && (model_->phase() != 2 && model_->phase() != 3)))
     2040    return 0; // switched off
     2041  OsiSolverInterface *solver = model_->solver();
     2042  const double *lower = solver->getColLower();
     2043  const double *upper = solver->getColUpper();
     2044  const double *rowLower = solver->getRowLower();
     2045  const double *rowUpper = solver->getRowUpper();
     2046  const double *solution = solver->getColSolution();
     2047  const double *objective = solver->getObjCoefficients();
     2048  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
     2049  double primalTolerance;
     2050  solver->getDblParam(OsiPrimalTolerance, primalTolerance);
     2051  double useTolerance = primalTolerance;
     2052
     2053  int numberRows = matrix_.getNumRows();
     2054  assert(numberRows <= solver->getNumRows());
     2055  if (numberRows == 0) {
     2056    return 0;
     2057  }
     2058  int numberIntegers = model_->numberIntegers();
     2059  const int *integerVariable = model_->integerVariable();
     2060  int i;
     2061  double direction = solver->getObjSense();
     2062  //double newSolutionValue = direction*solver->getObjValue();
     2063  int returnCode = 0;
     2064  // Column copy
     2065  const double *element = matrix_.getElements();
     2066  const int *row = matrix_.getIndices();
     2067  const CoinBigIndex *columnStart = matrix_.getVectorStarts();
     2068  const int *columnLength = matrix_.getVectorLengths();
     2069  // Row copy
     2070  const double *elementByRow = matrixByRow_.getElements();
     2071  const int *column = matrixByRow_.getIndices();
     2072  const CoinBigIndex *rowStart = matrixByRow_.getVectorStarts();
     2073  const int *rowLength = matrixByRow_.getVectorLengths();
     2074
     2075  // Get solution array for heuristic solution
     2076  int numberColumns = solver->getNumCols();
     2077  double *newSolution = new double[numberColumns];
     2078  memcpy(newSolution, solution, numberColumns * sizeof(double));
     2079
     2080  double *rowActivity = new double[numberRows];
     2081  memset(rowActivity, 0, numberRows * sizeof(double));
     2082  for (i = 0; i < numberColumns; i++) {
     2083    CoinBigIndex j;
     2084    double value = newSolution[i];
     2085    if (value < lower[i]) {
     2086      value = lower[i];
     2087      newSolution[i] = value;
     2088    } else if (value > upper[i]) {
     2089      value = upper[i];
     2090      newSolution[i] = value;
     2091    }
     2092    if (value) {
     2093      for (j = columnStart[i];
     2094           j < columnStart[i] + columnLength[i]; j++) {
     2095        int iRow = row[j];
     2096        rowActivity[iRow] += value * element[j];
     2097      }
     2098    }
     2099  }
     2100  // check was feasible - if not adjust (cleaning may move)
     2101  for (i = 0; i < numberRows; i++) {
     2102    if (rowActivity[i] < rowLower[i]) {
     2103      //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
     2104      rowActivity[i] = rowLower[i];
     2105    } else if (rowActivity[i] > rowUpper[i]) {
     2106      //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
     2107      rowActivity[i] = rowUpper[i];
     2108    }
     2109  }
     2110  for (i = 0; i < numberIntegers; i++) {
     2111    int iColumn = integerVariable[i];
     2112    double value = newSolution[iColumn];
     2113    //double thisTolerance = integerTolerance;
     2114    if (fabs(floor(value + 0.5) - value) > integerTolerance) {
     2115      double below = floor(value);
     2116      double newValue = newSolution[iColumn];
     2117      double cost = direction * objective[iColumn];
     2118      double move;
     2119      if (cost > 0.0) {
     2120        // try up
     2121        move = 1.0 - (value - below);
     2122      } else if (cost < 0.0) {
     2123        // try down
     2124        move = below - value;
     2125      } else {
     2126        // won't be able to move unless we can grab another variable
     2127        double randomNumber = randomNumberGenerator_.randomDouble();
     2128        // which way?
     2129        if (randomNumber < 0.5)
     2130          move = below - value;
     2131        else
     2132          move = 1.0 - (value - below);
     2133      }
     2134      newValue += move;
     2135      newSolution[iColumn] = newValue;
     2136      newSolutionValue += move * cost;
     2137      CoinBigIndex j;
     2138      for (j = columnStart[iColumn];
     2139           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     2140        int iRow = row[j];
     2141        rowActivity[iRow] += move * element[j];
     2142      }
     2143    }
     2144  }
     2145
     2146  double penalty = 0.0;
     2147  // see if feasible - just using singletons
     2148  for (i = 0; i < numberRows; i++) {
     2149    double value = rowActivity[i];
     2150    double thisInfeasibility = 0.0;
     2151    if (value < rowLower[i] - primalTolerance)
     2152      thisInfeasibility = value - rowLower[i];
     2153    else if (value > rowUpper[i] + primalTolerance)
     2154      thisInfeasibility = value - rowUpper[i];
     2155    if (thisInfeasibility) {
     2156      // See if there are any slacks I can use to fix up
     2157      // maybe put in coding for multiple slacks?
     2158      double bestCost = 1.0e50;
     2159      CoinBigIndex k;
     2160      int iBest = -1;
     2161      double addCost = 0.0;
     2162      double newValue = 0.0;
     2163      double changeRowActivity = 0.0;
     2164      double absInfeasibility = fabs(thisInfeasibility);
     2165      for (k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
     2166        int iColumn = column[k];
     2167        // See if all elements help
     2168        if (columnLength[iColumn] == 1) {
     2169          double currentValue = newSolution[iColumn];
     2170          double elementValue = elementByRow[k];
     2171          double lowerValue = lower[iColumn];
     2172          double upperValue = upper[iColumn];
     2173          double gap = rowUpper[i] - rowLower[i];
     2174          double absElement = fabs(elementValue);
     2175          if (thisInfeasibility * elementValue > 0.0) {
     2176            // we want to reduce
     2177            if ((currentValue - lowerValue) * absElement >= absInfeasibility) {
     2178              // possible - check if integer
     2179              double distance = absInfeasibility / absElement;
     2180              double thisCost = -direction * objective[iColumn] * distance;
     2181              if (isHeuristicInteger(solver, iColumn)) {
     2182                distance = ceil(distance - useTolerance);
     2183                if (currentValue - distance >= lowerValue - useTolerance) {
     2184                  if (absInfeasibility - distance * absElement < -gap - useTolerance)
     2185                    thisCost = 1.0e100; // no good
     2186                  else
     2187                    thisCost = -direction * objective[iColumn] * distance;
     2188                } else {
     2189                  thisCost = 1.0e100; // no good
     2190                }
     2191              }
     2192              if (thisCost < bestCost) {
     2193                bestCost = thisCost;
     2194                iBest = iColumn;
     2195                addCost = thisCost;
     2196                newValue = currentValue - distance;
     2197                changeRowActivity = -distance * elementValue;
     2198              }
     2199            }
     2200          } else {
     2201            // we want to increase
     2202            if ((upperValue - currentValue) * absElement >= absInfeasibility) {
     2203              // possible - check if integer
     2204              double distance = absInfeasibility / absElement;
     2205              double thisCost = direction * objective[iColumn] * distance;
     2206              if (isHeuristicInteger(solver, iColumn)) {
     2207                distance = ceil(distance - 1.0e-7);
     2208                assert(currentValue - distance <= upperValue + useTolerance);
     2209                if (absInfeasibility - distance * absElement < -gap - useTolerance)
     2210                  thisCost = 1.0e100; // no good
     2211                else
     2212                  thisCost = direction * objective[iColumn] * distance;
     2213              }
     2214              if (thisCost < bestCost) {
     2215                bestCost = thisCost;
     2216                iBest = iColumn;
     2217                addCost = thisCost;
     2218                newValue = currentValue + distance;
     2219                changeRowActivity = distance * elementValue;
     2220              }
     2221            }
     2222          }
     2223        }
     2224      }
     2225      if (iBest >= 0) {
     2226        /*printf("Infeasibility of %g on row %d cost %g\n",
     2227                  thisInfeasibility,i,addCost);*/
     2228        newSolution[iBest] = newValue;
     2229        thisInfeasibility = 0.0;
     2230        newSolutionValue += addCost;
     2231        rowActivity[i] += changeRowActivity;
     2232      }
     2233      penalty += fabs(thisInfeasibility);
     2234    }
     2235  }
     2236  if (penalty) {
     2237    // see if feasible using any
     2238    // first continuous
     2239    double penaltyChange = 0.0;
     2240    int iColumn;
     2241    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     2242      if (isHeuristicInteger(solver, iColumn))
     2243        continue;
     2244      double currentValue = newSolution[iColumn];
     2245      double lowerValue = lower[iColumn];
     2246      double upperValue = upper[iColumn];
     2247      CoinBigIndex j;
     2248      int anyBadDown = 0;
     2249      int anyBadUp = 0;
     2250      double upImprovement = 0.0;
     2251      double downImprovement = 0.0;
     2252      for (j = columnStart[iColumn];
     2253           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     2254        int iRow = row[j];
     2255        if (rowUpper[iRow] > rowLower[iRow]) {
     2256          double value = element[j];
     2257          if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
     2258            // infeasible above
     2259            downImprovement += value;
     2260            upImprovement -= value;
     2261            if (value > 0.0)
     2262              anyBadUp++;
     2263            else
     2264              anyBadDown++;
     2265          } else if (rowActivity[iRow] > rowUpper[iRow] - primalTolerance) {
     2266            // feasible at ub
     2267            if (value > 0.0) {
     2268              upImprovement -= value;
     2269              anyBadUp++;
     2270            } else {
     2271              downImprovement += value;
     2272              anyBadDown++;
     2273            }
     2274          } else if (rowActivity[iRow] > rowLower[iRow] + primalTolerance) {
     2275            // feasible in interior
     2276          } else if (rowActivity[iRow] > rowLower[iRow] - primalTolerance) {
     2277            // feasible at lb
     2278            if (value < 0.0) {
     2279              upImprovement += value;
     2280              anyBadUp++;
     2281            } else {
     2282              downImprovement -= value;
     2283              anyBadDown++;
     2284            }
     2285          } else {
     2286            // infeasible below
     2287            downImprovement -= value;
     2288            upImprovement += value;
     2289            if (value < 0.0)
     2290              anyBadUp++;
     2291            else
     2292              anyBadDown++;
     2293          }
     2294        } else {
     2295          // equality row
     2296          double value = element[j];
     2297          if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
     2298            // infeasible above
     2299            downImprovement += value;
     2300            upImprovement -= value;
     2301            if (value > 0.0)
     2302              anyBadUp++;
     2303            else
     2304              anyBadDown++;
     2305          } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
     2306            // infeasible below
     2307            downImprovement -= value;
     2308            upImprovement += value;
     2309            if (value < 0.0)
     2310              anyBadUp++;
     2311            else
     2312              anyBadDown++;
     2313          } else {
     2314            // feasible - no good
     2315            anyBadUp = -1;
     2316            anyBadDown = -1;
     2317            break;
     2318          }
     2319        }
     2320      }
     2321      // could change tests for anyBad
     2322      if (anyBadUp)
     2323        upImprovement = 0.0;
     2324      if (anyBadDown)
     2325        downImprovement = 0.0;
     2326      double way = 0.0;
     2327      double improvement = 0.0;
     2328      if (downImprovement > 0.0 && currentValue > lowerValue) {
     2329        way = -1.0;
     2330        improvement = downImprovement;
     2331      } else if (upImprovement > 0.0 && currentValue < upperValue) {
     2332        way = 1.0;
     2333        improvement = upImprovement;
     2334      }
     2335      if (way) {
     2336        // can improve
     2337        double distance;
     2338        if (way > 0.0)
     2339          distance = upperValue - currentValue;
     2340        else
     2341          distance = currentValue - lowerValue;
     2342        for (j = columnStart[iColumn];
     2343             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     2344          int iRow = row[j];
     2345          double value = element[j] * way;
     2346          if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
     2347            // infeasible above
     2348            assert(value < 0.0);
     2349            double gap = rowActivity[iRow] - rowUpper[iRow];
     2350            if (gap + value * distance < 0.0)
     2351              distance = -gap / value;
     2352          } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
     2353            // infeasible below
     2354            assert(value > 0.0);
     2355            double gap = rowActivity[iRow] - rowLower[iRow];
     2356            if (gap + value * distance > 0.0)
     2357              distance = -gap / value;
     2358          } else {
     2359            // feasible
     2360            if (value > 0) {
     2361              double gap = rowActivity[iRow] - rowUpper[iRow];
     2362              if (gap + value * distance > 0.0)
     2363                distance = -gap / value;
     2364            } else {
     2365              double gap = rowActivity[iRow] - rowLower[iRow];
     2366              if (gap + value * distance < 0.0)
     2367                distance = -gap / value;
     2368            }
     2369          }
     2370        }
     2371        //move
     2372        penaltyChange += improvement * distance;
     2373        distance *= way;
     2374        newSolution[iColumn] += distance;
     2375        newSolutionValue += direction * objective[iColumn] * distance;
     2376        for (j = columnStart[iColumn];
     2377             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     2378          int iRow = row[j];
     2379          double value = element[j];
     2380          rowActivity[iRow] += distance * value;
     2381        }
     2382      }
     2383    }
     2384    // and now all if improving
     2385    double lastChange = penaltyChange ? 1.0 : 0.0;
     2386    int numberPasses = 0;
     2387    while (lastChange > 1.0e-2 && numberPasses < 1000) {
     2388      lastChange = 0;
     2389      numberPasses++;
     2390      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     2391        bool isInteger = isHeuristicInteger(solver, iColumn);
     2392        double currentValue = newSolution[iColumn];
     2393        double lowerValue = lower[iColumn];
     2394        double upperValue = upper[iColumn];
     2395        CoinBigIndex j;
     2396        int anyBadDown = 0;
     2397        int anyBadUp = 0;
     2398        double upImprovement = 0.0;
     2399        double downImprovement = 0.0;
     2400        for (j = columnStart[iColumn];
     2401             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     2402          int iRow = row[j];
     2403          double value = element[j];
     2404          if (isInteger) {
     2405            if (value > 0.0) {
     2406              if (rowActivity[iRow] + value > rowUpper[iRow] + primalTolerance)
     2407                anyBadUp++;
     2408              if (rowActivity[iRow] - value < rowLower[iRow] - primalTolerance)
     2409                anyBadDown++;
     2410            } else {
     2411              if (rowActivity[iRow] - value > rowUpper[iRow] + primalTolerance)
     2412                anyBadDown++;
     2413              if (rowActivity[iRow] + value < rowLower[iRow] - primalTolerance)
     2414                anyBadUp++;
     2415            }
     2416          }
     2417          if (rowUpper[iRow] > rowLower[iRow]) {
     2418            if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
     2419              // infeasible above
     2420              downImprovement += value;
     2421              upImprovement -= value;
     2422              if (value > 0.0)
     2423                anyBadUp++;
     2424              else
     2425                anyBadDown++;
     2426            } else if (rowActivity[iRow] > rowUpper[iRow] - primalTolerance) {
     2427              // feasible at ub
     2428              if (value > 0.0) {
     2429                upImprovement -= value;
     2430                anyBadUp++;
     2431              } else {
     2432                downImprovement += value;
     2433                anyBadDown++;
     2434              }
     2435            } else if (rowActivity[iRow] > rowLower[iRow] + primalTolerance) {
     2436              // feasible in interior
     2437            } else if (rowActivity[iRow] > rowLower[iRow] - primalTolerance) {
     2438              // feasible at lb
     2439              if (value < 0.0) {
     2440                upImprovement += value;
     2441                anyBadUp++;
     2442              } else {
     2443                downImprovement -= value;
     2444                anyBadDown++;
     2445              }
     2446            } else {
     2447              // infeasible below
     2448              downImprovement -= value;
     2449              upImprovement += value;
     2450              if (value < 0.0)
     2451                anyBadUp++;
     2452              else
     2453                anyBadDown++;
     2454            }
     2455          } else {
     2456            // equality row
     2457            if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
     2458              // infeasible above
     2459              downImprovement += value;
     2460              upImprovement -= value;
     2461              if (value > 0.0)
     2462                anyBadUp++;
     2463              else
     2464                anyBadDown++;
     2465            } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
     2466              // infeasible below
     2467              downImprovement -= value;
     2468              upImprovement += value;
     2469              if (value < 0.0)
     2470                anyBadUp++;
     2471              else
     2472                anyBadDown++;
     2473            } else {
     2474              // feasible - no good
     2475              anyBadUp = -1;
     2476              anyBadDown = -1;
     2477              break;
     2478            }
     2479          }
     2480        }
     2481        // could change tests for anyBad
     2482        if (anyBadUp)
     2483          upImprovement = 0.0;
     2484        if (anyBadDown)
     2485          downImprovement = 0.0;
     2486        double way = 0.0;
     2487        double improvement = 0.0;
     2488        if (downImprovement > 0.0 && currentValue > lowerValue) {
     2489          way = -1.0;
     2490          improvement = downImprovement;
     2491          if (isInteger && currentValue < lowerValue + 0.99)
     2492            continue; // no good
     2493        } else if (upImprovement > 0.0 && currentValue < upperValue) {
     2494          way = 1.0;
     2495          improvement = upImprovement;
     2496          if (isInteger && currentValue > upperValue - 0.99)
     2497            continue; // no good
     2498        }
     2499        if (way) {
     2500          // can improve
     2501          double distance = COIN_DBL_MAX;
     2502          for (j = columnStart[iColumn];
     2503               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     2504            int iRow = row[j];
     2505            double value = element[j] * way;
     2506            if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
     2507              // infeasible above
     2508              assert(value < 0.0);
     2509              double gap = rowActivity[iRow] - rowUpper[iRow];
     2510              if (gap + value * distance < 0.0) {
     2511                // If integer then has to move by 1
     2512                if (!isInteger)
     2513                  distance = -gap / value;
     2514                else
     2515                  distance = CoinMax(-gap / value, 1.0);
     2516              }
     2517            } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
     2518              // infeasible below
     2519              assert(value > 0.0);
     2520              double gap = rowActivity[iRow] - rowLower[iRow];
     2521              if (gap + value * distance > 0.0) {
     2522                // If integer then has to move by 1
     2523                if (!isInteger)
     2524                  distance = -gap / value;
     2525                else
     2526                  distance = CoinMax(-gap / value, 1.0);
     2527              }
     2528            } else {
     2529              // feasible
     2530              if (value > 0) {
     2531                double gap = rowActivity[iRow] - rowUpper[iRow];
     2532                if (gap + value * distance > 0.0)
     2533                  distance = -gap / value;
     2534              } else {
     2535                double gap = rowActivity[iRow] - rowLower[iRow];
     2536                if (gap + value * distance < 0.0)
     2537                  distance = -gap / value;
     2538              }
     2539            }
     2540          }
     2541          if (isInteger)
     2542            distance = floor(distance + 1.05e-8);
     2543          if (!distance) {
     2544            // should never happen
     2545            //printf("zero distance in CbcRounding - debug\n");
     2546          }
     2547          //move
     2548          lastChange += improvement * distance;
     2549          distance *= way;
     2550          newSolution[iColumn] += distance;
     2551          newSolutionValue += direction * objective[iColumn] * distance;
     2552          for (j = columnStart[iColumn];
     2553               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     2554            int iRow = row[j];
     2555            double value = element[j];
     2556            rowActivity[iRow] += distance * value;
     2557          }
     2558        }
     2559      }
     2560      penaltyChange += lastChange;
     2561    }
     2562    penalty -= penaltyChange;
     2563    if (penalty < 1.0e-5 * fabs(penaltyChange)) {
     2564      // recompute
     2565      penalty = 0.0;
     2566      for (i = 0; i < numberRows; i++) {
     2567        double value = rowActivity[i];
     2568        if (value < rowLower[i] - primalTolerance)
     2569          penalty += rowLower[i] - value;
     2570        else if (value > rowUpper[i] + primalTolerance)
     2571          penalty += value - rowUpper[i];
     2572      }
     2573    }
     2574  }
     2575
     2576  // Could also set SOS (using random) and repeat
     2577  if (!penalty) {
     2578    // See if we can do better
     2579    //seed_++;
     2580    //CoinSeedRandom(seed_);
     2581    // Random number between 0 and 1.
     2582    double randomNumber = randomNumberGenerator_.randomDouble();
     2583    int iPass;
     2584    int start[2];
     2585    int end[2];
     2586    int iRandom = static_cast<int>(randomNumber * (static_cast<double>(numberIntegers)));
     2587    start[0] = iRandom;
     2588    end[0] = numberIntegers;
     2589    start[1] = 0;
     2590    end[1] = iRandom;
     2591    for (iPass = 0; iPass < 2; iPass++) {
     2592      int i;
     2593      for (i = start[iPass]; i < end[iPass]; i++) {
     2594        int iColumn = integerVariable[i];
     2595#ifndef NDEBUG
     2596        double value = newSolution[iColumn];
     2597        assert(fabs(floor(value + 0.5) - value) <= integerTolerance);
     2598#endif
     2599        double cost = direction * objective[iColumn];
     2600        double move = 0.0;
     2601        if (cost > 0.0)
     2602          move = -1.0;
     2603        else if (cost < 0.0)
     2604          move = 1.0;
     2605        while (move) {
     2606          bool good = true;
     2607          double newValue = newSolution[iColumn] + move;
     2608          if (newValue < lower[iColumn] - useTolerance || newValue > upper[iColumn] + useTolerance) {
     2609            move = 0.0;
     2610          } else {
     2611            // see if we can move
     2612            CoinBigIndex j;
     2613            for (j = columnStart[iColumn];
     2614                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     2615              int iRow = row[j];
     2616              double newActivity = rowActivity[iRow] + move * element[j];
     2617              if (newActivity < rowLower[iRow] - primalTolerance || newActivity > rowUpper[iRow] + primalTolerance) {
     2618                good = false;
     2619                break;
     2620              }
     2621            }
     2622            if (good) {
     2623              newSolution[iColumn] = newValue;
     2624              newSolutionValue += move * cost;
     2625              CoinBigIndex j;
     2626              for (j = columnStart[iColumn];
     2627                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     2628                int iRow = row[j];
     2629                rowActivity[iRow] += move * element[j];
     2630              }
     2631            } else {
     2632              move = 0.0;
     2633            }
     2634          }
     2635        }
     2636      }
     2637    }
     2638    // Just in case of some stupidity
     2639    double objOffset = 0.0;
     2640    solver->getDblParam(OsiObjOffset, objOffset);
     2641    newSolutionValue = -objOffset;
     2642    for (i = 0; i < numberColumns; i++)
     2643      newSolutionValue += objective[i] * newSolution[i];
     2644    newSolutionValue *= direction;
     2645    //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
     2646    if (newSolutionValue < solutionValue) {
     2647      // paranoid check
     2648      memset(rowActivity, 0, numberRows * sizeof(double));
     2649      for (i = 0; i < numberColumns; i++) {
    21332650        CoinBigIndex j;
    21342651        double value = newSolution[i];
    2135         if (value < lower[i]) {
    2136             value = lower[i];
    2137             newSolution[i] = value;
    2138         } else if (value > upper[i]) {
    2139             value = upper[i];
    2140             newSolution[i] = value;
    2141         }
    21422652        if (value) {
    2143             for (j = columnStart[i];
    2144                     j < columnStart[i] + columnLength[i]; j++) {
    2145                 int iRow = row[j];
    2146                 rowActivity[iRow] += value * element[j];
    2147             }
    2148         }
     2653          for (j = columnStart[i];
     2654               j < columnStart[i] + columnLength[i]; j++) {
     2655            int iRow = row[j];
     2656            rowActivity[iRow] += value * element[j];
     2657          }
     2658        }
     2659      }
     2660      // check was approximately feasible
     2661      bool feasible = true;
     2662      for (i = 0; i < numberRows; i++) {
     2663        if (rowActivity[i] < rowLower[i]) {
     2664          if (rowActivity[i] < rowLower[i] - 1000.0 * primalTolerance)
     2665            feasible = false;
     2666        } else if (rowActivity[i] > rowUpper[i]) {
     2667          if (rowActivity[i] > rowUpper[i] + 1000.0 * primalTolerance)
     2668            feasible = false;
     2669        }
     2670      }
     2671      if (feasible) {
     2672        // new solution
     2673        memcpy(betterSolution, newSolution, numberColumns * sizeof(double));
     2674        solutionValue = newSolutionValue;
     2675        //printf("** Solution of %g found by rounding\n",newSolutionValue);
     2676        returnCode = 1;
     2677      } else {
     2678        // Can easily happen
     2679        //printf("Debug CbcRounding giving bad solution\n");
     2680      }
     2681    }
     2682  }
     2683#ifdef NEW_ROUNDING
     2684  if (!returnCode) {
     2685#ifdef JJF_ZERO
     2686    // back to starting point
     2687    memcpy(newSolution, solution, numberColumns * sizeof(double));
     2688    memset(rowActivity, 0, numberRows * sizeof(double));
     2689    for (i = 0; i < numberColumns; i++) {
     2690      int j;
     2691      double value = newSolution[i];
     2692      if (value < lower[i]) {
     2693        value = lower[i];
     2694        newSolution[i] = value;
     2695      } else if (value > upper[i]) {
     2696        value = upper[i];
     2697        newSolution[i] = value;
     2698      }
     2699      if (value) {
     2700        for (j = columnStart[i];
     2701             j < columnStart[i] + columnLength[i]; j++) {
     2702          int iRow = row[j];
     2703          rowActivity[iRow] += value * element[j];
     2704        }
     2705      }
    21492706    }
    21502707    // check was feasible - if not adjust (cleaning may move)
    21512708    for (i = 0; i < numberRows; i++) {
    2152         if (rowActivity[i] < rowLower[i]) {
    2153             //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
    2154             rowActivity[i] = rowLower[i];
    2155         } else if (rowActivity[i] > rowUpper[i]) {
    2156             //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
    2157             rowActivity[i] = rowUpper[i];
    2158         }
    2159     }
    2160     for (i = 0; i < numberIntegers; i++) {
    2161         int iColumn = integerVariable[i];
    2162         double value = newSolution[iColumn];
    2163         //double thisTolerance = integerTolerance;
    2164         if (fabs(floor(value + 0.5) - value) > integerTolerance) {
    2165             double below = floor(value);
    2166             double newValue = newSolution[iColumn];
    2167             double cost = direction * objective[iColumn];
    2168             double move;
    2169             if (cost > 0.0) {
    2170                 // try up
    2171                 move = 1.0 - (value - below);
    2172             } else if (cost < 0.0) {
    2173                 // try down
    2174                 move = below - value;
    2175             } else {
    2176                 // won't be able to move unless we can grab another variable
    2177                 double randomNumber = randomNumberGenerator_.randomDouble();
    2178                 // which way?
    2179                 if (randomNumber < 0.5)
    2180                     move = below - value;
    2181                 else
    2182                     move = 1.0 - (value - below);
    2183             }
    2184             newValue += move;
    2185             newSolution[iColumn] = newValue;
    2186             newSolutionValue += move * cost;
    2187             CoinBigIndex j;
    2188             for (j = columnStart[iColumn];
    2189                     j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    2190                 int iRow = row[j];
    2191                 rowActivity[iRow] += move * element[j];
    2192             }
    2193         }
    2194     }
    2195 
    2196     double penalty = 0.0;
    2197     // see if feasible - just using singletons
    2198     for (i = 0; i < numberRows; i++) {
    2199         double value = rowActivity[i];
    2200         double thisInfeasibility = 0.0;
    2201         if (value < rowLower[i] - primalTolerance)
    2202             thisInfeasibility = value - rowLower[i];
    2203         else if (value > rowUpper[i] + primalTolerance)
    2204             thisInfeasibility = value - rowUpper[i];
    2205         if (thisInfeasibility) {
    2206             // See if there are any slacks I can use to fix up
    2207             // maybe put in coding for multiple slacks?
    2208             double bestCost = 1.0e50;
    2209             CoinBigIndex k;
    2210             int iBest = -1;
    2211             double addCost = 0.0;
    2212             double newValue = 0.0;
    2213             double changeRowActivity = 0.0;
    2214             double absInfeasibility = fabs(thisInfeasibility);
    2215             for (k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
    2216                 int iColumn = column[k];
    2217                 // See if all elements help
    2218                 if (columnLength[iColumn] == 1) {
    2219                     double currentValue = newSolution[iColumn];
    2220                     double elementValue = elementByRow[k];
    2221                     double lowerValue = lower[iColumn];
    2222                     double upperValue = upper[iColumn];
    2223                     double gap = rowUpper[i] - rowLower[i];
    2224                     double absElement = fabs(elementValue);
    2225                     if (thisInfeasibility*elementValue > 0.0) {
    2226                         // we want to reduce
    2227                         if ((currentValue - lowerValue)*absElement >= absInfeasibility) {
    2228                             // possible - check if integer
    2229                             double distance = absInfeasibility / absElement;
    2230                             double thisCost = -direction * objective[iColumn] * distance;
    2231                             if (isHeuristicInteger(solver,iColumn)) {
    2232                                 distance = ceil(distance - useTolerance);
    2233                                 if (currentValue - distance >= lowerValue - useTolerance) {
    2234                                     if (absInfeasibility - distance*absElement < -gap - useTolerance)
    2235                                         thisCost = 1.0e100; // no good
    2236                                     else
    2237                                         thisCost = -direction * objective[iColumn] * distance;
    2238                                 } else {
    2239                                     thisCost = 1.0e100; // no good
    2240                                 }
    2241                             }
    2242                             if (thisCost < bestCost) {
    2243                                 bestCost = thisCost;
    2244                                 iBest = iColumn;
    2245                                 addCost = thisCost;
    2246                                 newValue = currentValue - distance;
    2247                                 changeRowActivity = -distance * elementValue;
    2248                             }
    2249                         }
    2250                     } else {
    2251                         // we want to increase
    2252                         if ((upperValue - currentValue)*absElement >= absInfeasibility) {
    2253                             // possible - check if integer
    2254                             double distance = absInfeasibility / absElement;
    2255                             double thisCost = direction * objective[iColumn] * distance;
    2256                             if (isHeuristicInteger(solver,iColumn)) {
    2257                                 distance = ceil(distance - 1.0e-7);
    2258                                 assert (currentValue - distance <= upperValue + useTolerance);
    2259                                 if (absInfeasibility - distance*absElement < -gap - useTolerance)
    2260                                     thisCost = 1.0e100; // no good
    2261                                 else
    2262                                     thisCost = direction * objective[iColumn] * distance;
    2263                             }
    2264                             if (thisCost < bestCost) {
    2265                                 bestCost = thisCost;
    2266                                 iBest = iColumn;
    2267                                 addCost = thisCost;
    2268                                 newValue = currentValue + distance;
    2269                                 changeRowActivity = distance * elementValue;
    2270                             }
    2271                         }
    2272                     }
    2273                 }
    2274             }
    2275             if (iBest >= 0) {
    2276                 /*printf("Infeasibility of %g on row %d cost %g\n",
    2277                   thisInfeasibility,i,addCost);*/
    2278                 newSolution[iBest] = newValue;
    2279                 thisInfeasibility = 0.0;
    2280                 newSolutionValue += addCost;
    2281                 rowActivity[i] += changeRowActivity;
    2282             }
    2283             penalty += fabs(thisInfeasibility);
    2284         }
    2285     }
    2286     if (penalty) {
    2287         // see if feasible using any
    2288         // first continuous
    2289         double penaltyChange = 0.0;
    2290         int iColumn;
    2291         for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    2292             if (isHeuristicInteger(solver,iColumn))
    2293                 continue;
    2294             double currentValue = newSolution[iColumn];
    2295             double lowerValue = lower[iColumn];
    2296             double upperValue = upper[iColumn];
    2297             CoinBigIndex j;
    2298             int anyBadDown = 0;
    2299             int anyBadUp = 0;
    2300             double upImprovement = 0.0;
    2301             double downImprovement = 0.0;
    2302             for (j = columnStart[iColumn];
    2303                     j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    2304                 int iRow = row[j];
    2305                 if (rowUpper[iRow] > rowLower[iRow]) {
    2306                     double value = element[j];
    2307                     if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
    2308                         // infeasible above
    2309                         downImprovement += value;
    2310                         upImprovement -= value;
    2311                         if (value > 0.0)
    2312                             anyBadUp++;
    2313                         else
    2314                             anyBadDown++;
    2315                     } else if (rowActivity[iRow] > rowUpper[iRow] - primalTolerance) {
    2316                         // feasible at ub
    2317                         if (value > 0.0) {
    2318                             upImprovement -= value;
    2319                             anyBadUp++;
    2320                         } else {
    2321                             downImprovement += value;
    2322                             anyBadDown++;
    2323                         }
    2324                     } else if (rowActivity[iRow] > rowLower[iRow] + primalTolerance) {
    2325                         // feasible in interior
    2326                     } else if (rowActivity[iRow] > rowLower[iRow] - primalTolerance) {
    2327                         // feasible at lb
    2328                         if (value < 0.0) {
    2329                             upImprovement += value;
    2330                             anyBadUp++;
    2331                         } else {
    2332                             downImprovement -= value;
    2333                             anyBadDown++;
    2334                         }
    2335                     } else {
    2336                         // infeasible below
    2337                         downImprovement -= value;
    2338                         upImprovement += value;
    2339                         if (value < 0.0)
    2340                             anyBadUp++;
    2341                         else
    2342                             anyBadDown++;
    2343                     }
    2344                 } else {
    2345                     // equality row
    2346                     double value = element[j];
    2347                     if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
    2348                         // infeasible above
    2349                         downImprovement += value;
    2350                         upImprovement -= value;
    2351                         if (value > 0.0)
    2352                             anyBadUp++;
    2353                         else
    2354                             anyBadDown++;
    2355                     } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
    2356                         // infeasible below
    2357                         downImprovement -= value;
    2358                         upImprovement += value;
    2359                         if (value < 0.0)
    2360                             anyBadUp++;
    2361                         else
    2362                             anyBadDown++;
    2363                     } else {
    2364                         // feasible - no good
    2365                         anyBadUp = -1;
    2366                         anyBadDown = -1;
    2367                         break;
    2368                     }
    2369                 }
    2370             }
    2371             // could change tests for anyBad
    2372             if (anyBadUp)
    2373                 upImprovement = 0.0;
    2374             if (anyBadDown)
    2375                 downImprovement = 0.0;
    2376             double way = 0.0;
    2377             double improvement = 0.0;
    2378             if (downImprovement > 0.0 && currentValue > lowerValue) {
    2379                 way = -1.0;
    2380                 improvement = downImprovement;
    2381             } else if (upImprovement > 0.0 && currentValue < upperValue) {
    2382                 way = 1.0;
    2383                 improvement = upImprovement;
    2384             }
    2385             if (way) {
    2386                 // can improve
    2387                 double distance;
    2388                 if (way > 0.0)
    2389                     distance = upperValue - currentValue;
    2390                 else
    2391                     distance = currentValue - lowerValue;
    2392                 for (j = columnStart[iColumn];
    2393                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    2394                     int iRow = row[j];
    2395                     double value = element[j] * way;
    2396                     if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
    2397                         // infeasible above
    2398                         assert (value < 0.0);
    2399                         double gap = rowActivity[iRow] - rowUpper[iRow];
    2400                         if (gap + value*distance < 0.0)
    2401                             distance = -gap / value;
    2402                     } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
    2403                         // infeasible below
    2404                         assert (value > 0.0);
    2405                         double gap = rowActivity[iRow] - rowLower[iRow];
    2406                         if (gap + value*distance > 0.0)
    2407                             distance = -gap / value;
    2408                     } else {
    2409                         // feasible
    2410                         if (value > 0) {
    2411                             double gap = rowActivity[iRow] - rowUpper[iRow];
    2412                             if (gap + value*distance > 0.0)
    2413                                 distance = -gap / value;
    2414                         } else {
    2415                             double gap = rowActivity[iRow] - rowLower[iRow];
    2416                             if (gap + value*distance < 0.0)
    2417                                 distance = -gap / value;
    2418                         }
    2419                     }
    2420                 }
    2421                 //move
    2422                 penaltyChange += improvement * distance;
    2423                 distance *= way;
    2424                 newSolution[iColumn] += distance;
    2425                 newSolutionValue += direction * objective[iColumn] * distance;
    2426                 for (j = columnStart[iColumn];
    2427                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    2428                     int iRow = row[j];
    2429                     double value = element[j];
    2430                     rowActivity[iRow] += distance * value;
    2431                 }
    2432             }
    2433         }
    2434         // and now all if improving
    2435         double lastChange = penaltyChange ? 1.0 : 0.0;
    2436         int numberPasses=0;
    2437         while (lastChange > 1.0e-2 && numberPasses < 1000) {
    2438             lastChange = 0;
    2439             numberPasses++;
    2440             for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    2441                 bool isInteger = isHeuristicInteger(solver,iColumn);
    2442                 double currentValue = newSolution[iColumn];
    2443                 double lowerValue = lower[iColumn];
    2444                 double upperValue = upper[iColumn];
    2445                 CoinBigIndex j;
    2446                 int anyBadDown = 0;
    2447                 int anyBadUp = 0;
    2448                 double upImprovement = 0.0;
    2449                 double downImprovement = 0.0;
    2450                 for (j = columnStart[iColumn];
    2451                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    2452                     int iRow = row[j];
    2453                     double value = element[j];
    2454                     if (isInteger) {
    2455                         if (value > 0.0) {
    2456                             if (rowActivity[iRow] + value > rowUpper[iRow] + primalTolerance)
    2457                                 anyBadUp++;
    2458                             if (rowActivity[iRow] - value < rowLower[iRow] - primalTolerance)
    2459                                 anyBadDown++;
    2460                         } else {
    2461                             if (rowActivity[iRow] - value > rowUpper[iRow] + primalTolerance)
    2462                                 anyBadDown++;
    2463                             if (rowActivity[iRow] + value < rowLower[iRow] - primalTolerance)
    2464                                 anyBadUp++;
    2465                         }
    2466                     }
    2467                     if (rowUpper[iRow] > rowLower[iRow]) {
    2468                         if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
    2469                             // infeasible above
    2470                             downImprovement += value;
    2471                             upImprovement -= value;
    2472                             if (value > 0.0)
    2473                                 anyBadUp++;
    2474                             else
    2475                                 anyBadDown++;
    2476                         } else if (rowActivity[iRow] > rowUpper[iRow] - primalTolerance) {
    2477                             // feasible at ub
    2478                             if (value > 0.0) {
    2479                                 upImprovement -= value;
    2480                                 anyBadUp++;
    2481                             } else {
    2482                                 downImprovement += value;
    2483                                 anyBadDown++;
    2484                             }
    2485                         } else if (rowActivity[iRow] > rowLower[iRow] + primalTolerance) {
    2486                             // feasible in interior
    2487                         } else if (rowActivity[iRow] > rowLower[iRow] - primalTolerance) {
    2488                             // feasible at lb
    2489                             if (value < 0.0) {
    2490                                 upImprovement += value;
    2491                                 anyBadUp++;
    2492                             } else {
    2493                                 downImprovement -= value;
    2494                                 anyBadDown++;
    2495                             }
    2496                         } else {
    2497                             // infeasible below
    2498                             downImprovement -= value;
    2499                             upImprovement += value;
    2500                             if (value < 0.0)
    2501                                 anyBadUp++;
    2502                             else
    2503                                 anyBadDown++;
    2504                         }
    2505                     } else {
    2506                         // equality row
    2507                         if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
    2508                             // infeasible above
    2509                             downImprovement += value;
    2510                             upImprovement -= value;
    2511                             if (value > 0.0)
    2512                                 anyBadUp++;
    2513                             else
    2514                                 anyBadDown++;
    2515                         } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
    2516                             // infeasible below
    2517                             downImprovement -= value;
    2518                             upImprovement += value;
    2519                             if (value < 0.0)
    2520                                 anyBadUp++;
    2521                             else
    2522                                 anyBadDown++;
    2523                         } else {
    2524                             // feasible - no good
    2525                             anyBadUp = -1;
    2526                             anyBadDown = -1;
    2527                             break;
    2528                         }
    2529                     }
    2530                 }
    2531                 // could change tests for anyBad
    2532                 if (anyBadUp)
    2533                     upImprovement = 0.0;
    2534                 if (anyBadDown)
    2535                     downImprovement = 0.0;
    2536                 double way = 0.0;
    2537                 double improvement = 0.0;
    2538                 if (downImprovement > 0.0 && currentValue > lowerValue) {
    2539                     way = -1.0;
    2540                     improvement = downImprovement;
    2541                     if (isInteger&&currentValue<lowerValue+0.99)
    2542                       continue; // no good
    2543                 } else if (upImprovement > 0.0 && currentValue < upperValue) {
    2544                     way = 1.0;
    2545                     improvement = upImprovement;
    2546                     if (isInteger&&currentValue>upperValue-0.99)
    2547                       continue; // no good
    2548                 }
    2549                 if (way) {
    2550                     // can improve
    2551                     double distance = COIN_DBL_MAX;
    2552                     for (j = columnStart[iColumn];
    2553                             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    2554                         int iRow = row[j];
    2555                         double value = element[j] * way;
    2556                         if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
    2557                             // infeasible above
    2558                             assert (value < 0.0);
    2559                             double gap = rowActivity[iRow] - rowUpper[iRow];
    2560                             if (gap + value*distance < 0.0) {
    2561                                 // If integer then has to move by 1
    2562                                 if (!isInteger)
    2563                                     distance = -gap / value;
    2564                                 else
    2565                                     distance = CoinMax(-gap / value, 1.0);
    2566                             }
    2567                         } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
    2568                             // infeasible below
    2569                             assert (value > 0.0);
    2570                             double gap = rowActivity[iRow] - rowLower[iRow];
    2571                             if (gap + value*distance > 0.0) {
    2572                                 // If integer then has to move by 1
    2573                                 if (!isInteger)
    2574                                     distance = -gap / value;
    2575                                 else
    2576                                     distance = CoinMax(-gap / value, 1.0);
    2577                             }
    2578                         } else {
    2579                             // feasible
    2580                             if (value > 0) {
    2581                                 double gap = rowActivity[iRow] - rowUpper[iRow];
    2582                                 if (gap + value*distance > 0.0)
    2583                                     distance = -gap / value;
    2584                             } else {
    2585                                 double gap = rowActivity[iRow] - rowLower[iRow];
    2586                                 if (gap + value*distance < 0.0)
    2587                                     distance = -gap / value;
    2588                             }
    2589                         }
    2590                     }
    2591                     if (isInteger)
    2592                         distance = floor(distance + 1.05e-8);
    2593                     if (!distance) {
    2594                         // should never happen
    2595                         //printf("zero distance in CbcRounding - debug\n");
    2596                     }
    2597                     //move
    2598                     lastChange += improvement * distance;
    2599                     distance *= way;
    2600                     newSolution[iColumn] += distance;
    2601                     newSolutionValue += direction * objective[iColumn] * distance;
    2602                     for (j = columnStart[iColumn];
    2603                             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    2604                         int iRow = row[j];
    2605                         double value = element[j];
    2606                         rowActivity[iRow] += distance * value;
    2607                     }
    2608                 }
    2609             }
    2610             penaltyChange += lastChange;
    2611         }
    2612         penalty -= penaltyChange;
    2613         if (penalty < 1.0e-5*fabs(penaltyChange)) {
    2614             // recompute
    2615             penalty = 0.0;
    2616             for (i = 0; i < numberRows; i++) {
    2617                 double value = rowActivity[i];
    2618                 if (value < rowLower[i] - primalTolerance)
    2619                     penalty += rowLower[i] - value;
    2620                 else if (value > rowUpper[i] + primalTolerance)
    2621                     penalty += value - rowUpper[i];
    2622             }
    2623         }
    2624     }
    2625 
    2626     // Could also set SOS (using random) and repeat
    2627     if (!penalty) {
    2628         // See if we can do better
    2629         //seed_++;
    2630         //CoinSeedRandom(seed_);
    2631         // Random number between 0 and 1.
    2632         double randomNumber = randomNumberGenerator_.randomDouble();
    2633         int iPass;
    2634         int start[2];
    2635         int end[2];
    2636         int iRandom = static_cast<int> (randomNumber * (static_cast<double> (numberIntegers)));
    2637         start[0] = iRandom;
    2638         end[0] = numberIntegers;
    2639         start[1] = 0;
    2640         end[1] = iRandom;
    2641         for (iPass = 0; iPass < 2; iPass++) {
    2642             int i;
    2643             for (i = start[iPass]; i < end[iPass]; i++) {
    2644                 int iColumn = integerVariable[i];
    2645 #ifndef NDEBUG
    2646                 double value = newSolution[iColumn];
    2647                 assert (fabs(floor(value + 0.5) - value) <= integerTolerance);
    2648 #endif
    2649                 double cost = direction * objective[iColumn];
    2650                 double move = 0.0;
    2651                 if (cost > 0.0)
    2652                     move = -1.0;
    2653                 else if (cost < 0.0)
    2654                     move = 1.0;
    2655                 while (move) {
    2656                     bool good = true;
    2657                     double newValue = newSolution[iColumn] + move;
    2658                     if (newValue < lower[iColumn] - useTolerance ||
    2659                             newValue > upper[iColumn] + useTolerance) {
    2660                         move = 0.0;
    2661                     } else {
    2662                         // see if we can move
    2663                         CoinBigIndex j;
    2664                         for (j = columnStart[iColumn];
    2665                                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    2666                             int iRow = row[j];
    2667                             double newActivity = rowActivity[iRow] + move * element[j];
    2668                             if (newActivity < rowLower[iRow] - primalTolerance ||
    2669                                     newActivity > rowUpper[iRow] + primalTolerance) {
    2670                                 good = false;
    2671                                 break;
    2672                             }
    2673                         }
    2674                         if (good) {
    2675                             newSolution[iColumn] = newValue;
    2676                             newSolutionValue += move * cost;
    2677                             CoinBigIndex j;
    2678                             for (j = columnStart[iColumn];
    2679                                     j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    2680                                 int iRow = row[j];
    2681                                 rowActivity[iRow] += move * element[j];
    2682                             }
    2683                         } else {
    2684                             move = 0.0;
    2685                         }
    2686                     }
    2687                 }
    2688             }
    2689         }
    2690         // Just in case of some stupidity
    2691         double objOffset = 0.0;
    2692         solver->getDblParam(OsiObjOffset, objOffset);
    2693         newSolutionValue = -objOffset;
    2694         for ( i = 0 ; i < numberColumns ; i++ )
    2695             newSolutionValue += objective[i] * newSolution[i];
    2696         newSolutionValue *= direction;
    2697         //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
    2698         if (newSolutionValue < solutionValue) {
    2699             // paranoid check
    2700             memset(rowActivity, 0, numberRows*sizeof(double));
    2701             for (i = 0; i < numberColumns; i++) {
    2702                 CoinBigIndex j;
    2703                 double value = newSolution[i];
    2704                 if (value) {
    2705                     for (j = columnStart[i];
    2706                             j < columnStart[i] + columnLength[i]; j++) {
    2707                         int iRow = row[j];
    2708                         rowActivity[iRow] += value * element[j];
    2709                     }
    2710                 }
    2711             }
    2712             // check was approximately feasible
    2713             bool feasible = true;
    2714             for (i = 0; i < numberRows; i++) {
    2715                 if (rowActivity[i] < rowLower[i]) {
    2716                     if (rowActivity[i] < rowLower[i] - 1000.0*primalTolerance)
    2717                         feasible = false;
    2718                 } else if (rowActivity[i] > rowUpper[i]) {
    2719                     if (rowActivity[i] > rowUpper[i] + 1000.0*primalTolerance)
    2720                         feasible = false;
    2721                 }
    2722             }
    2723             if (feasible) {
    2724                 // new solution
    2725                 memcpy(betterSolution, newSolution, numberColumns*sizeof(double));
    2726                 solutionValue = newSolutionValue;
    2727                 //printf("** Solution of %g found by rounding\n",newSolutionValue);
    2728                 returnCode = 1;
    2729             } else {
    2730                 // Can easily happen
    2731                 //printf("Debug CbcRounding giving bad solution\n");
    2732             }
    2733         }
    2734     }
     2709      if (rowActivity[i] < rowLower[i]) {
     2710        //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
     2711        rowActivity[i] = rowLower[i];
     2712      } else if (rowActivity[i] > rowUpper[i]) {
     2713        //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
     2714        rowActivity[i] = rowUpper[i];
     2715      }
     2716    }
     2717#endif
     2718    int *candidate = new int[numberColumns];
     2719    int nCandidate = 0;
     2720    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     2721      bool isInteger = isHeuristicInteger(solver, iColumn);
     2722      if (isInteger) {
     2723        double currentValue = newSolution[iColumn];
     2724        if (fabs(currentValue - floor(currentValue + 0.5)) > 1.0e-8)
     2725          candidate[nCandidate++] = iColumn;
     2726      }
     2727    }
     2728    if (true) {
     2729      // Rounding as in Berthold
     2730      while (nCandidate) {
     2731        double infeasibility = 1.0e-7;
     2732        int iRow = -1;
     2733        for (i = 0; i < numberRows; i++) {
     2734          double value = 0.0;
     2735          if (rowActivity[i] < rowLower[i]) {
     2736            value = rowLower[i] - rowActivity[i];
     2737          } else if (rowActivity[i] > rowUpper[i]) {
     2738            value = rowActivity[i] - rowUpper[i];
     2739          }
     2740          if (value > infeasibility) {
     2741            infeasibility = value;
     2742            iRow = i;
     2743          }
     2744        }
     2745        if (iRow >= 0) {
     2746          // infeasible
     2747        } else {
     2748          // feasible
     2749        }
     2750      }
     2751    } else {
     2752      // Shifting as in Berthold
     2753    }
     2754    delete[] candidate;
     2755  }
     2756#endif
     2757  delete[] newSolution;
     2758  delete[] rowActivity;
     2759  return returnCode;
     2760}
     2761// update model
     2762void CbcRounding::setModel(CbcModel *model)
     2763{
     2764  model_ = model;
     2765  // Get a copy of original matrix (and by row for rounding);
     2766  assert(model_->solver());
     2767  if (model_->solver()->getNumRows()) {
     2768    matrix_ = *model_->solver()->getMatrixByCol();
     2769    matrixByRow_ = *model_->solver()->getMatrixByRow();
     2770    // make sure model okay for heuristic
     2771    validate();
     2772  }
     2773}
     2774// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
     2775void CbcRounding::validate()
     2776{
     2777  if (model_ && (when() % 100) < 10) {
     2778    if (model_->numberIntegers() != model_->numberObjects() && (model_->numberObjects() || (model_->specialOptions() & 1024) == 0)) {
     2779      int numberOdd = 0;
     2780      for (int i = 0; i < model_->numberObjects(); i++) {
     2781        if (!model_->object(i)->canDoHeuristics())
     2782          numberOdd++;
     2783      }
     2784      if (numberOdd)
     2785        setWhen(0);
     2786    }
     2787  }
    27352788#ifdef NEW_ROUNDING
    2736     if (!returnCode) {
    2737 #ifdef JJF_ZERO
    2738         // back to starting point
    2739         memcpy(newSolution, solution, numberColumns*sizeof(double));
    2740         memset(rowActivity, 0, numberRows*sizeof(double));
    2741         for (i = 0; i < numberColumns; i++) {
    2742             int j;
    2743             double value = newSolution[i];
    2744             if (value < lower[i]) {
    2745                 value = lower[i];
    2746                 newSolution[i] = value;
    2747             } else if (value > upper[i]) {
    2748                 value = upper[i];
    2749                 newSolution[i] = value;
    2750             }
    2751             if (value) {
    2752                 for (j = columnStart[i];
    2753                         j < columnStart[i] + columnLength[i]; j++) {
    2754                     int iRow = row[j];
    2755                     rowActivity[iRow] += value * element[j];
    2756                 }
    2757             }
    2758         }
    2759         // check was feasible - if not adjust (cleaning may move)
    2760         for (i = 0; i < numberRows; i++) {
    2761             if (rowActivity[i] < rowLower[i]) {
    2762                 //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
    2763                 rowActivity[i] = rowLower[i];
    2764             } else if (rowActivity[i] > rowUpper[i]) {
    2765                 //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
    2766                 rowActivity[i] = rowUpper[i];
    2767             }
    2768         }
    2769 #endif
    2770         int * candidate = new int [numberColumns];
    2771         int nCandidate = 0;
    2772         for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    2773             bool isInteger = isHeuristicInteger(solver,iColumn);
    2774             if (isInteger) {
    2775                 double currentValue = newSolution[iColumn];
    2776                 if (fabs(currentValue - floor(currentValue + 0.5)) > 1.0e-8)
    2777                     candidate[nCandidate++] = iColumn;
    2778             }
    2779         }
    2780         if (true) {
    2781             // Rounding as in Berthold
    2782             while (nCandidate) {
    2783                 double infeasibility = 1.0e-7;
    2784                 int iRow = -1;
    2785                 for (i = 0; i < numberRows; i++) {
    2786                     double value = 0.0;
    2787                     if (rowActivity[i] < rowLower[i]) {
    2788                         value = rowLower[i] - rowActivity[i];
    2789                     } else if (rowActivity[i] > rowUpper[i]) {
    2790                         value = rowActivity[i] - rowUpper[i];
    2791                     }
    2792                     if (value > infeasibility) {
    2793                         infeasibility = value;
    2794                         iRow = i;
    2795                     }
    2796                 }
    2797                 if (iRow >= 0) {
    2798                     // infeasible
    2799                 } else {
    2800                     // feasible
    2801                 }
    2802             }
    2803         } else {
    2804             // Shifting as in Berthold
    2805         }
    2806         delete [] candidate;
    2807     }
    2808 #endif
    2809     delete [] newSolution;
    2810     delete [] rowActivity;
    2811     return returnCode;
    2812 }
    2813 // update model
    2814 void CbcRounding::setModel(CbcModel * model)
    2815 {
    2816     model_ = model;
    2817     // Get a copy of original matrix (and by row for rounding);
    2818     assert(model_->solver());
    2819     if (model_->solver()->getNumRows()) {
    2820         matrix_ = *model_->solver()->getMatrixByCol();
    2821         matrixByRow_ = *model_->solver()->getMatrixByRow();
    2822         // make sure model okay for heuristic
    2823         validate();
    2824     }
    2825 }
    2826 // Validate model i.e. sets when_ to 0 if necessary (may be NULL)
    2827 void
    2828 CbcRounding::validate()
    2829 {
    2830     if (model_ && (when() % 100) < 10) {
    2831         if (model_->numberIntegers() !=
    2832                 model_->numberObjects() && (model_->numberObjects() ||
    2833                                             (model_->specialOptions()&1024) == 0)) {
    2834             int numberOdd = 0;
    2835             for (int i = 0; i < model_->numberObjects(); i++) {
    2836                 if (!model_->object(i)->canDoHeuristics())
    2837                     numberOdd++;
    2838             }
    2839             if (numberOdd)
    2840                 setWhen(0);
    2841         }
    2842     }
    2843 #ifdef NEW_ROUNDING
    2844     int numberColumns = matrix_.getNumCols();
    2845     down_ = new unsigned short [numberColumns];
    2846     up_ = new unsigned short [numberColumns];
    2847     equal_ = new unsigned short [numberColumns];
    2848     // Column copy
    2849     const double * element = matrix_.getElements();
    2850     const int * row = matrix_.getIndices();
    2851     const CoinBigIndex * columnStart = matrix_.getVectorStarts();
    2852     const int * columnLength = matrix_.getVectorLengths();
    2853     const double * rowLower = model.solver()->getRowLower();
    2854     const double * rowUpper = model.solver()->getRowUpper();
    2855     for (int i = 0; i < numberColumns; i++) {
    2856         int down = 0;
    2857         int up = 0;
    2858         int equal = 0;
    2859         if (columnLength[i] > 65535) {
    2860             equal[0] = 65535;
    2861             break; // unlikely to work
    2862         }
    2863         for (CoinBigIndex j = columnStart[i];
    2864                 j < columnStart[i] + columnLength[i]; j++) {
    2865             int iRow = row[j];
    2866             if (rowLower[iRow] > -1.0e20 && rowUpper[iRow] < 1.0e20) {
    2867                 equal++;
    2868             } else if (element[j] > 0.0) {
    2869                 if (rowUpper[iRow] < 1.0e20)
    2870                     up++;
    2871                 else
    2872                     down--;
    2873             } else {
    2874                 if (rowLower[iRow] > -1.0e20)
    2875                     up++;
    2876                 else
    2877                     down--;
    2878             }
    2879         }
    2880         down_[i] = (unsigned short) down;
    2881         up_[i] = (unsigned short) up;
    2882         equal_[i] = (unsigned short) equal;
    2883     }
     2789  int numberColumns = matrix_.getNumCols();
     2790  down_ = new unsigned short[numberColumns];
     2791  up_ = new unsigned short[numberColumns];
     2792  equal_ = new unsigned short[numberColumns];
     2793  // Column copy
     2794  const double *element = matrix_.getElements();
     2795  const int *row = matrix_.getIndices();
     2796  const CoinBigIndex *columnStart = matrix_.getVectorStarts();
     2797  const int *columnLength = matrix_.getVectorLengths();
     2798  const double *rowLower = model.solver()->getRowLower();
     2799  const double *rowUpper = model.solver()->getRowUpper();
     2800  for (int i = 0; i < numberColumns; i++) {
     2801    int down = 0;
     2802    int up = 0;
     2803    int equal = 0;
     2804    if (columnLength[i] > 65535) {
     2805      equal[0] = 65535;
     2806      break; // unlikely to work
     2807    }
     2808    for (CoinBigIndex j = columnStart[i];
     2809         j < columnStart[i] + columnLength[i]; j++) {
     2810      int iRow = row[j];
     2811      if (rowLower[iRow] > -1.0e20 && rowUpper[iRow] < 1.0e20) {
     2812        equal++;
     2813      } else if (element[j] > 0.0) {
     2814        if (rowUpper[iRow] < 1.0e20)
     2815          up++;
     2816        else
     2817          down--;
     2818      } else {
     2819        if (rowLower[iRow] > -1.0e20)
     2820          up++;
     2821        else
     2822          down--;
     2823      }
     2824    }
     2825    down_[i] = (unsigned short)down;
     2826    up_[i] = (unsigned short)up;
     2827    equal_[i] = (unsigned short)equal;
     2828  }
    28842829#else
    2885     down_ = NULL;
    2886     up_ = NULL;
    2887     equal_ = NULL;
     2830  down_ = NULL;
     2831  up_ = NULL;
     2832  equal_ = NULL;
    28882833#endif
    28892834}
     
    28912836// Default Constructor
    28922837CbcHeuristicPartial::CbcHeuristicPartial()
    2893         : CbcHeuristic()
    2894 {
    2895     fixPriority_ = 10000;
     2838  : CbcHeuristic()
     2839{
     2840  fixPriority_ = 10000;
    28962841}
    28972842
    28982843// Constructor from model
    2899 CbcHeuristicPartial::CbcHeuristicPartial(CbcModel & model, int fixPriority, int numberNodes)
    2900         : CbcHeuristic(model)
    2901 {
    2902     fixPriority_ = fixPriority;
    2903     setNumberNodes(numberNodes);
    2904     validate();
     2844CbcHeuristicPartial::CbcHeuristicPartial(CbcModel &model, int fixPriority, int numberNodes)
     2845  : CbcHeuristic(model)
     2846{
     2847  fixPriority_ = fixPriority;
     2848  setNumberNodes(numberNodes);
     2849  validate();
    29052850}
    29062851
    29072852// Destructor
    2908 CbcHeuristicPartial::~CbcHeuristicPartial ()
     2853CbcHeuristicPartial::~CbcHeuristicPartial()
    29092854{
    29102855}
     
    29142859CbcHeuristicPartial::clone() const
    29152860{
    2916     return new CbcHeuristicPartial(*this);
     2861  return new CbcHeuristicPartial(*this);
    29172862}
    29182863// Create C++ lines to get to current state
    2919 void
    2920 CbcHeuristicPartial::generateCpp( FILE * fp)
    2921 {
    2922     CbcHeuristicPartial other;
    2923     fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n");
    2924     fprintf(fp, "3  CbcHeuristicPartial partial(*cbcModel);\n");
    2925     CbcHeuristic::generateCpp(fp, "partial");
    2926     if (fixPriority_ != other.fixPriority_)
    2927         fprintf(fp, "3  partial.setFixPriority(%d);\n", fixPriority_);
    2928     else
    2929         fprintf(fp, "4  partial.setFixPriority(%d);\n", fixPriority_);
    2930     fprintf(fp, "3  cbcModel->addHeuristic(&partial);\n");
     2864void CbcHeuristicPartial::generateCpp(FILE *fp)
     2865{
     2866  CbcHeuristicPartial other;
     2867  fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n");
     2868  fprintf(fp, "3  CbcHeuristicPartial partial(*cbcModel);\n");
     2869  CbcHeuristic::generateCpp(fp, "partial");
     2870  if (fixPriority_ != other.fixPriority_)
     2871    fprintf(fp, "3  partial.setFixPriority(%d);\n", fixPriority_);
     2872  else
     2873    fprintf(fp, "4  partial.setFixPriority(%d);\n", fixPriority_);
     2874  fprintf(fp, "3  cbcModel->addHeuristic(&partial);\n");
    29312875}
    29322876//#define NEW_PARTIAL
    29332877// Copy constructor
    2934 CbcHeuristicPartial::CbcHeuristicPartial(const CbcHeuristicPartial & rhs)
    2935         :
    2936         CbcHeuristic(rhs),
    2937         fixPriority_(rhs.fixPriority_)
     2878CbcHeuristicPartial::CbcHeuristicPartial(const CbcHeuristicPartial &rhs)
     2879  : CbcHeuristic(rhs)
     2880  , fixPriority_(rhs.fixPriority_)
    29382881{
    29392882}
     
    29412884// Assignment operator
    29422885CbcHeuristicPartial &
    2943 CbcHeuristicPartial::operator=( const CbcHeuristicPartial & rhs)
    2944 {
    2945     if (this != &rhs) {
    2946         CbcHeuristic::operator=(rhs);
    2947         fixPriority_ = rhs.fixPriority_;
    2948     }
    2949     return *this;
     2886CbcHeuristicPartial::operator=(const CbcHeuristicPartial &rhs)
     2887{
     2888  if (this != &rhs) {
     2889    CbcHeuristic::operator=(rhs);
     2890    fixPriority_ = rhs.fixPriority_;
     2891  }
     2892  return *this;
    29502893}
    29512894
    29522895// Resets stuff if model changes
    2953 void
    2954 CbcHeuristicPartial::resetModel(CbcModel * model)
    2955 {
    2956     model_ = model;
    2957     // Get a copy of original matrix (and by row for partial);
    2958     assert(model_->solver());
    2959     validate();
     2896void CbcHeuristicPartial::resetModel(CbcModel *model)
     2897{
     2898  model_ = model;
     2899  // Get a copy of original matrix (and by row for partial);
     2900  assert(model_->solver());
     2901  validate();
    29602902}
    29612903// See if partial will give solution
     
    29652907// Fix values if asked for
    29662908// Returns 1 if solution, 0 if not
    2967 int
    2968 CbcHeuristicPartial::solution(double & solutionValue,
    2969                               double * betterSolution)
    2970 {
    2971     // Return if already done
    2972     if (fixPriority_ < 0)
    2973         return 0; // switched off
     2909int CbcHeuristicPartial::solution(double &solutionValue,
     2910  double *betterSolution)
     2911{
     2912  // Return if already done
     2913  if (fixPriority_ < 0)
     2914    return 0; // switched off
    29742915#ifdef HEURISTIC_INFORM
    2975     printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
    2976            heuristicName(),numRuns_,numCouldRun_,when_);
    2977 #endif
    2978     const double * hotstartSolution = model_->hotstartSolution();
    2979     const int * hotstartPriorities = model_->hotstartPriorities();
    2980     if (!hotstartSolution)
    2981         return 0;
    2982     OsiSolverInterface * solver = model_->solver();
    2983 
    2984     int numberIntegers = model_->numberIntegers();
    2985     const int * integerVariable = model_->integerVariable();
    2986 
    2987     OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
    2988     const double * colLower = newSolver->getColLower();
    2989     const double * colUpper = newSolver->getColUpper();
    2990 
    2991     double primalTolerance;
    2992     solver->getDblParam(OsiPrimalTolerance, primalTolerance);
    2993 
    2994     int i;
    2995     int numberFixed = 0;
    2996     int returnCode = 0;
    2997 
    2998     for (i = 0; i < numberIntegers; i++) {
    2999         int iColumn = integerVariable[i];
    3000         if (abs(hotstartPriorities[iColumn]) <= fixPriority_) {
    3001             double value = hotstartSolution[iColumn];
    3002             double lower = colLower[iColumn];
    3003             double upper = colUpper[iColumn];
    3004             value = CoinMax(value, lower);
    3005             value = CoinMin(value, upper);
    3006             if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
    3007                 value = floor(value + 0.5);
    3008                 newSolver->setColLower(iColumn, value);
    3009                 newSolver->setColUpper(iColumn, value);
    3010                 numberFixed++;
    3011             }
    3012         }
    3013     }
    3014     if (numberFixed > numberIntegers / 5 - 100000000) {
     2916  printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
     2917    heuristicName(), numRuns_, numCouldRun_, when_);
     2918#endif
     2919  const double *hotstartSolution = model_->hotstartSolution();
     2920  const int *hotstartPriorities = model_->hotstartPriorities();
     2921  if (!hotstartSolution)
     2922    return 0;
     2923  OsiSolverInterface *solver = model_->solver();
     2924
     2925  int numberIntegers = model_->numberIntegers();
     2926  const int *integerVariable = model_->integerVariable();
     2927
     2928  OsiSolverInterface *newSolver = model_->continuousSolver()->clone();
     2929  const double *colLower = newSolver->getColLower();
     2930  const double *colUpper = newSolver->getColUpper();
     2931
     2932  double primalTolerance;
     2933  solver->getDblParam(OsiPrimalTolerance, primalTolerance);
     2934
     2935  int i;
     2936  int numberFixed = 0;
     2937  int returnCode = 0;
     2938
     2939  for (i = 0; i < numberIntegers; i++) {
     2940    int iColumn = integerVariable[i];
     2941    if (abs(hotstartPriorities[iColumn]) <= fixPriority_) {
     2942      double value = hotstartSolution[iColumn];
     2943      double lower = colLower[iColumn];
     2944      double upper = colUpper[iColumn];
     2945      value = CoinMax(value, lower);
     2946      value = CoinMin(value, upper);
     2947      if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
     2948        value = floor(value + 0.5);
     2949        newSolver->setColLower(iColumn, value);
     2950        newSolver->setColUpper(iColumn, value);
     2951        numberFixed++;
     2952      }
     2953    }
     2954  }
     2955  if (numberFixed > numberIntegers / 5 - 100000000) {
    30152956#ifdef COIN_DEVELOP
    3016         printf("%d integers fixed\n", numberFixed);
    3017 #endif
    3018         returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
    3019                                          model_->getCutoff(), "CbcHeuristicPartial");
    3020         if (returnCode < 0)
    3021             returnCode = 0; // returned on size
    3022         //printf("return code %d",returnCode);
    3023         if ((returnCode&2) != 0) {
    3024             // could add cut
    3025             returnCode &= ~2;
    3026             //printf("could add cut with %d elements (if all 0-1)\n",nFix);
    3027         } else {
    3028             //printf("\n");
    3029         }
    3030     }
    3031     fixPriority_ = -1; // switch off
    3032 
    3033     delete newSolver;
    3034     return returnCode;
     2957    printf("%d integers fixed\n", numberFixed);
     2958#endif
     2959    returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
     2960      model_->getCutoff(), "CbcHeuristicPartial");
     2961    if (returnCode < 0)
     2962      returnCode = 0; // returned on size
     2963    //printf("return code %d",returnCode);
     2964    if ((returnCode & 2) != 0) {
     2965      // could add cut
     2966      returnCode &= ~2;
     2967      //printf("could add cut with %d elements (if all 0-1)\n",nFix);
     2968    } else {
     2969      //printf("\n");
     2970    }
     2971  }
     2972  fixPriority_ = -1; // switch off
     2973
     2974  delete newSolver;
     2975  return returnCode;
    30352976}
    30362977// update model
    3037 void CbcHeuristicPartial::setModel(CbcModel * model)
    3038 {
    3039     model_ = model;
    3040     assert(model_->solver());
    3041     // make sure model okay for heuristic
    3042     validate();
     2978void CbcHeuristicPartial::setModel(CbcModel *model)
     2979{
     2980  model_ = model;
     2981  assert(model_->solver());
     2982  // make sure model okay for heuristic
     2983  validate();
    30432984}
    30442985// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
    3045 void
    3046 CbcHeuristicPartial::validate()
    3047 {
    3048     if (model_ && (when() % 100) < 10) {
    3049         if (model_->numberIntegers() !=
    3050                 model_->numberObjects())
    3051             setWhen(0);
    3052     }
    3053 }
    3054 bool
    3055 CbcHeuristicPartial::shouldHeurRun(int /*whereFrom*/)
    3056 {
    3057     return true;
     2986void CbcHeuristicPartial::validate()
     2987{
     2988  if (model_ && (when() % 100) < 10) {
     2989    if (model_->numberIntegers() != model_->numberObjects())
     2990      setWhen(0);
     2991  }
     2992}
     2993bool CbcHeuristicPartial::shouldHeurRun(int /*whereFrom*/)
     2994{
     2995  return true;
    30582996}
    30592997
    30602998// Default Constructor
    30612999CbcSerendipity::CbcSerendipity()
    3062         : CbcHeuristic()
     3000  : CbcHeuristic()
    30633001{
    30643002}
    30653003
    30663004// Constructor from model
    3067 CbcSerendipity::CbcSerendipity(CbcModel & model)
    3068         : CbcHeuristic(model)
     3005CbcSerendipity::CbcSerendipity(CbcModel &model)
     3006  : CbcHeuristic(model)
    30693007{
    30703008}
    30713009
    30723010// Destructor
    3073 CbcSerendipity::~CbcSerendipity ()
     3011CbcSerendipity::~CbcSerendipity()
    3