source: trunk/CbcBranchActual.cpp @ 122

Last change on this file since 122 was 122, checked in by forrest, 16 years ago

major changes

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 66.5 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#if defined(_MSC_VER)
4// Turn off compiler warning about long names
5#  pragma warning(disable:4786)
6#endif
7#include <cassert>
8#include <cmath>
9#include <cfloat>
10//#define CBC_DEBUG
11
12#include "OsiSolverInterface.hpp"
13#include "CbcModel.hpp"
14#include "CbcMessage.hpp"
15#include "CbcBranchActual.hpp"
16#include "CoinSort.hpp"
17#include "CoinError.hpp"
18
19// Default Constructor
20CbcClique::CbcClique ()
21  : CbcObject(),
22    numberMembers_(0),
23    numberNonSOSMembers_(0),
24    members_(NULL),
25    type_(NULL),
26    cliqueType_(-1),
27    slack_(-1)
28{
29}
30
31// Useful constructor (which are integer indices)
32CbcClique::CbcClique (CbcModel * model, int cliqueType, int numberMembers,
33           const int * which, const char * type, int identifier,int slack)
34  : CbcObject(model)
35{
36  id_=identifier;
37  numberMembers_=numberMembers;
38  if (numberMembers_) {
39    members_ = new int[numberMembers_];
40    memcpy(members_,which,numberMembers_*sizeof(int));
41    type_ = new char[numberMembers_];
42    if (type) {
43      memcpy(type_,type,numberMembers_*sizeof(char));
44    } else {
45      for (int i=0;i<numberMembers_;i++)
46        type_[i]=1;
47    }
48  } else {
49    members_ = NULL;
50    type_ = NULL;
51  }
52  // Find out how many non sos
53  int i;
54  numberNonSOSMembers_=0;
55  for (i=0;i<numberMembers_;i++)
56    if (!type_[i])
57      numberNonSOSMembers_++;
58  cliqueType_ = cliqueType;
59  slack_ = slack;
60}
61
62// Copy constructor
63CbcClique::CbcClique ( const CbcClique & rhs)
64  :CbcObject(rhs)
65{
66  numberMembers_ = rhs.numberMembers_;
67  numberNonSOSMembers_ = rhs.numberNonSOSMembers_;
68  if (numberMembers_) {
69    members_ = new int[numberMembers_];
70    memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
71    type_ = new char[numberMembers_];
72    memcpy(type_,rhs.type_,numberMembers_*sizeof(char));
73  } else {
74    members_ = NULL;
75    type_ = NULL;
76  }
77  cliqueType_ = rhs.cliqueType_;
78  slack_ = rhs.slack_;
79}
80
81// Clone
82CbcObject *
83CbcClique::clone() const
84{
85  return new CbcClique(*this);
86}
87
88// Assignment operator
89CbcClique & 
90CbcClique::operator=( const CbcClique& rhs)
91{
92  if (this!=&rhs) {
93    CbcObject::operator=(rhs);
94    delete [] members_;
95    delete [] type_;
96    numberMembers_ = rhs.numberMembers_;
97    numberNonSOSMembers_ = rhs.numberNonSOSMembers_;
98    if (numberMembers_) {
99      members_ = new int[numberMembers_];
100      memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
101      type_ = new char[numberMembers_];
102      memcpy(type_,rhs.type_,numberMembers_*sizeof(char));
103    } else {
104      members_ = NULL;
105      type_ = NULL;
106    }
107    cliqueType_ = rhs.cliqueType_;
108    slack_ = rhs.slack_;
109  }
110  return *this;
111}
112
113// Destructor
114CbcClique::~CbcClique ()
115{
116  delete [] members_;
117  delete [] type_;
118}
119
120// Infeasibility - large is 0.5
121double 
122CbcClique::infeasibility(int & preferredWay) const
123{
124  int numberUnsatis=0, numberFree=0;
125  int j;
126  const int * integer = model_->integerVariable();
127  OsiSolverInterface * solver = model_->solver();
128  const double * solution = model_->currentSolution();
129  const double * lower = solver->getColLower();
130  const double * upper = solver->getColUpper();
131  double largestValue=0.0;
132  double integerTolerance = 
133    model_->getDblParam(CbcModel::CbcIntegerTolerance);
134  double * sort = new double[numberMembers_];
135
136  double slackValue=0.0;
137  for (j=0;j<numberMembers_;j++) {
138    int sequence = members_[j];
139    int iColumn = integer[sequence];
140    double value = solution[iColumn];
141    value = CoinMax(value, lower[iColumn]);
142    value = CoinMin(value, upper[iColumn]);
143    double nearest = floor(value+0.5);
144    double distance = fabs(value-nearest);
145    if (distance>integerTolerance) {
146      if (!type_[j])
147        value = 1.0-value; // non SOS
148      // if slack then choose that
149      if (j==slack_&&value>0.05)
150        slackValue = value;
151      largestValue = CoinMax(value,largestValue);
152      sort[numberUnsatis++]=-value;
153    } else if (upper[iColumn]>lower[iColumn]) {
154      numberFree++;
155    }
156  }
157  preferredWay=1;
158  double otherWay = 0.0;
159  if (numberUnsatis) {
160    // sort
161    std::sort(sort,sort+numberUnsatis);
162    for (j=0;j<numberUnsatis;j++) {
163      if ((j&1)!=0)
164        otherWay += -sort[j];
165    }
166    // Need to think more
167    double value = 0.2*numberUnsatis+0.01*(numberMembers_-numberFree);
168    if (fabs(largestValue-0.5)<0.1) {
169      // close to half
170      value +=0.1;
171    }
172    if (slackValue) {
173      // branching on slack
174      value += slackValue;
175    }
176    // scale other way
177    otherWay *= value/(1.0-otherWay);
178    delete [] sort;
179    return value;
180  } else {
181    delete [] sort;
182    return 0.0; // satisfied
183  }
184}
185
186// This looks at solution and sets bounds to contain solution
187void 
188CbcClique::feasibleRegion()
189{
190  int j;
191  const int * integer = model_->integerVariable();
192  OsiSolverInterface * solver = model_->solver();
193  const double * solution = model_->currentSolution();
194  const double * lower = solver->getColLower();
195  const double * upper = solver->getColUpper();
196  double integerTolerance = 
197    model_->getDblParam(CbcModel::CbcIntegerTolerance);
198 
199  for (j=0;j<numberMembers_;j++) {
200    int sequence = members_[j];
201    int iColumn = integer[sequence];
202    double value = solution[iColumn];
203    value = CoinMax(value, lower[iColumn]);
204    value = CoinMin(value, upper[iColumn]);
205    double nearest = floor(value+0.5);
206    double distance = fabs(value-nearest);
207    assert(distance<=integerTolerance);
208    solver->setColLower(iColumn,nearest);
209    solver->setColUpper(iColumn,nearest);
210  }
211}
212
213
214// Creates a branching object
215CbcBranchingObject * 
216CbcClique::createBranch(int way) const
217{
218  int numberUnsatis=0;
219  int j;
220  int nUp=0;
221  int nDown=0;
222  int numberFree=numberMembers_;
223  const int * integer = model_->integerVariable();
224  OsiSolverInterface * solver = model_->solver();
225  const double * solution = model_->currentSolution();
226  const double * lower = solver->getColLower();
227  const double * upper = solver->getColUpper();
228  int * upList = new int[numberMembers_];
229  int * downList = new int[numberMembers_];
230  double * sort = new double[numberMembers_];
231  double integerTolerance = 
232      model_->getDblParam(CbcModel::CbcIntegerTolerance);
233
234  double slackValue=0.0;
235  for (j=0;j<numberMembers_;j++) {
236    int sequence = members_[j];
237    int iColumn = integer[sequence];
238    double value = solution[iColumn];
239    value = CoinMax(value, lower[iColumn]);
240    value = CoinMin(value, upper[iColumn]);
241    double nearest = floor(value+0.5);
242    double distance = fabs(value-nearest);
243    if (distance>integerTolerance) {
244      if (!type_[j])
245        value = 1.0-value; // non SOS
246      // if slack then choose that
247      if (j==slack_&&value>0.05)
248        slackValue = value;
249      value = -value; // for sort
250      upList[numberUnsatis]=j;
251      sort[numberUnsatis++]=value;
252    } else if (upper[iColumn]>lower[iColumn]) {
253      upList[--numberFree]=j;
254    }
255  }
256  assert (numberUnsatis);
257  if (!slackValue) {
258    // sort
259    CoinSort_2(sort,sort+numberUnsatis,upList);
260    // put first in up etc
261    int kWay=1;
262    for (j=0;j<numberUnsatis;j++) {
263      if (kWay>0) 
264        upList[nUp++]=upList[j];
265      else
266        downList[nDown++]=upList[j];
267      kWay = -kWay;
268    }
269    for (j=numberFree;j<numberMembers_;j++) {
270      if (kWay>0)
271        upList[nUp++]=upList[j];
272      else
273        downList[nDown++]=upList[j];
274      kWay = -kWay;
275    }
276  } else {
277    // put slack to 0 in first way
278    nUp = 1;
279    upList[0]=slack_;
280    for (j=0;j<numberUnsatis;j++) {
281      downList[nDown++]=upList[j];
282    }
283    for (j=numberFree;j<numberMembers_;j++) {
284      downList[nDown++]=upList[j];
285    }
286  }
287  // create object
288  CbcBranchingObject * branch;
289  if (numberMembers_ <=64)
290     branch = new CbcCliqueBranchingObject(model_,this,way,
291                                         nDown,downList,nUp,upList);
292  else
293    branch = new CbcLongCliqueBranchingObject(model_,this,way,
294                                            nDown,downList,nUp,upList);
295  delete [] upList;
296  delete [] downList;
297  delete [] sort;
298  return branch;
299}
300
301// Default Constructor
302CbcSOS::CbcSOS ()
303  : CbcObject(),
304    members_(NULL),
305    weights_(NULL),
306    numberMembers_(0),
307    sosType_(-1)
308{
309}
310
311// Useful constructor (which are indices)
312CbcSOS::CbcSOS (CbcModel * model,  int numberMembers,
313           const int * which, const double * weights, int identifier,int type)
314  : CbcObject(model),
315    numberMembers_(numberMembers),
316    sosType_(type)
317{
318  id_=identifier;
319  if (numberMembers_) {
320    members_ = new int[numberMembers_];
321    weights_ = new double[numberMembers_];
322    memcpy(members_,which,numberMembers_*sizeof(int));
323    if (weights) {
324      memcpy(weights_,weights,numberMembers_*sizeof(double));
325    } else {
326      for (int i=0;i<numberMembers_;i++)
327        weights_[i]=i;
328    }
329    // sort so weights increasing
330    CoinSort_2(weights_,weights_+numberMembers_,members_);
331  } else {
332    members_ = NULL;
333    weights_ = NULL;
334  }
335  assert (sosType_>0&&sosType_<3);
336}
337
338// Copy constructor
339CbcSOS::CbcSOS ( const CbcSOS & rhs)
340  :CbcObject(rhs)
341{
342  numberMembers_ = rhs.numberMembers_;
343  sosType_ = rhs.sosType_;
344  if (numberMembers_) {
345    members_ = new int[numberMembers_];
346    weights_ = new double[numberMembers_];
347    memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
348    memcpy(weights_,rhs.weights_,numberMembers_*sizeof(double));
349  } else {
350    members_ = NULL;
351    weights_ = NULL;
352  }
353}
354
355// Clone
356CbcObject *
357CbcSOS::clone() const
358{
359  return new CbcSOS(*this);
360}
361
362// Assignment operator
363CbcSOS & 
364CbcSOS::operator=( const CbcSOS& rhs)
365{
366  if (this!=&rhs) {
367    CbcObject::operator=(rhs);
368    delete [] members_;
369    delete [] weights_;
370    numberMembers_ = rhs.numberMembers_;
371    sosType_ = rhs.sosType_;
372    if (numberMembers_) {
373      members_ = new int[numberMembers_];
374      weights_ = new double[numberMembers_];
375      memcpy(members_,rhs.members_,numberMembers_*sizeof(int));
376      memcpy(weights_,rhs.weights_,numberMembers_*sizeof(double));
377    } else {
378      members_ = NULL;
379      weights_ = NULL;
380    }
381  }
382  return *this;
383}
384
385// Destructor
386CbcSOS::~CbcSOS ()
387{
388  delete [] members_;
389  delete [] weights_;
390}
391
392// Infeasibility - large is 0.5
393double 
394CbcSOS::infeasibility(int & preferredWay) const
395{
396  int j;
397  int firstNonZero=-1;
398  int lastNonZero = -1;
399  OsiSolverInterface * solver = model_->solver();
400  const double * solution = model_->currentSolution();
401  const double * lower = solver->getColLower();
402  const double * upper = solver->getColUpper();
403  //double largestValue=0.0;
404  double integerTolerance = 
405    model_->getDblParam(CbcModel::CbcIntegerTolerance);
406  double weight = 0.0;
407  double sum =0.0;
408
409  // check bounds etc
410  double lastWeight=-1.0e100;
411  for (j=0;j<numberMembers_;j++) {
412    int iColumn = members_[j];
413    if (lower[iColumn]&&(lower[iColumn]!=1.0||sosType_!=1))
414      throw CoinError("Non zero lower bound in SOS","infeasibility","CbcSOS");
415    if (lastWeight>=weights_[j]-1.0e-7)
416      throw CoinError("Weights too close together in SOS","infeasibility","CbcSOS");
417    double value = CoinMax(0.0,solution[iColumn]);
418    sum += value;
419    if (value>integerTolerance&&upper[iColumn]) {
420      // Possibly due to scaling a fixed variable might slip through
421      if (value>upper[iColumn]) {
422        // Could change to #ifdef CBC_DEBUG
423#ifndef NDEBUG
424        if (model_->messageHandler()->logLevel()>1)
425          printf("** Variable %d (%d) has value %g and upper bound of %g\n",
426                 iColumn,j,value,upper[iColumn]);
427#endif
428      } 
429      weight += weights_[j]*value;
430      if (firstNonZero<0)
431        firstNonZero=j;
432      lastNonZero=j;
433    }
434  }
435  preferredWay=1;
436  if (lastNonZero-firstNonZero>=sosType_) {
437    // find where to branch
438    assert (sum>0.0);
439    weight /= sum;
440    //int iWhere;
441    //for (iWhere=firstNonZero;iWhere<lastNonZero;iWhere++)
442    //if (weight<weights_[iWhere+1])
443    //break;
444    // probably best to use pseudo duals
445    double value = lastNonZero-firstNonZero+1;
446    value *= 0.5/((double) numberMembers_);
447    return value;
448  } else {
449    return 0.0; // satisfied
450  }
451}
452
453// This looks at solution and sets bounds to contain solution
454void 
455CbcSOS::feasibleRegion()
456{
457  int j;
458  int firstNonZero=-1;
459  int lastNonZero = -1;
460  OsiSolverInterface * solver = model_->solver();
461  const double * solution = model_->currentSolution();
462  //const double * lower = solver->getColLower();
463  const double * upper = solver->getColUpper();
464  double integerTolerance = 
465    model_->getDblParam(CbcModel::CbcIntegerTolerance);
466  double weight = 0.0;
467  double sum =0.0;
468
469  for (j=0;j<numberMembers_;j++) {
470    int iColumn = members_[j];
471    double value = CoinMax(0.0,solution[iColumn]);
472    sum += value;
473    if (value>integerTolerance&&upper[iColumn]) {
474      weight += weights_[j]*value;
475      if (firstNonZero<0)
476        firstNonZero=j;
477      lastNonZero=j;
478    }
479  }
480  assert (lastNonZero-firstNonZero<sosType_) ;
481  for (j=0;j<firstNonZero;j++) {
482    int iColumn = members_[j];
483    solver->setColUpper(iColumn,0.0);
484  }
485  for (j=lastNonZero+1;j<numberMembers_;j++) {
486    int iColumn = members_[j];
487    solver->setColUpper(iColumn,0.0);
488  }
489}
490
491
492// Creates a branching object
493CbcBranchingObject * 
494CbcSOS::createBranch(int way) const
495{
496  int j;
497  const double * solution = model_->currentSolution();
498  double integerTolerance = 
499      model_->getDblParam(CbcModel::CbcIntegerTolerance);
500  OsiSolverInterface * solver = model_->solver();
501  const double * upper = solver->getColUpper();
502  int firstNonFixed=-1;
503  int lastNonFixed=-1;
504  int firstNonZero=-1;
505  int lastNonZero = -1;
506  double weight = 0.0;
507  double sum =0.0;
508  for (j=0;j<numberMembers_;j++) {
509    int iColumn = members_[j];
510    if (upper[iColumn]) {
511      double value = CoinMax(0.0,solution[iColumn]);
512      sum += value;
513      if (firstNonFixed<0)
514        firstNonFixed=j;
515      lastNonFixed=j;
516      if (value>integerTolerance) {
517        weight += weights_[j]*value;
518        if (firstNonZero<0)
519          firstNonZero=j;
520        lastNonZero=j;
521      }
522    }
523  }
524  assert (lastNonZero-firstNonZero>=sosType_) ;
525  // find where to branch
526  assert (sum>0.0);
527  weight /= sum;
528  int iWhere;
529  double separator=0.0;
530  for (iWhere=firstNonZero;iWhere<lastNonZero;iWhere++) 
531    if (weight<weights_[iWhere+1])
532      break;
533  if (sosType_==1) {
534    // SOS 1
535    separator = 0.5 *(weights_[iWhere]+weights_[iWhere+1]);
536  } else {
537    // SOS 2
538    if (iWhere==firstNonFixed)
539      iWhere++;;
540    if (iWhere==lastNonFixed-1)
541      iWhere = lastNonFixed-2;
542    separator = weights_[iWhere+1];
543  }
544  // create object
545  CbcBranchingObject * branch;
546  branch = new CbcSOSBranchingObject(model_,this,way,separator);
547  return branch;
548}
549
550
551
552/** Default Constructor
553
554  Equivalent to an unspecified binary variable.
555*/
556CbcSimpleInteger::CbcSimpleInteger ()
557  : CbcObject(),
558    sequence_(-1),
559    columnNumber_(-1),
560    originalLower_(0.0),
561    originalUpper_(1.0),
562    breakEven_(0.5)
563{
564}
565
566/** Useful constructor
567
568  Loads actual upper & lower bounds for the specified variable.
569*/
570CbcSimpleInteger::CbcSimpleInteger (CbcModel * model, int sequence,
571                                    int iColumn, double breakEven)
572  : CbcObject(model)
573{
574  sequence_ = sequence ;
575  id_ = iColumn; // set as well
576  columnNumber_ = iColumn ;
577  originalLower_ = model_->solver()->getColLower()[columnNumber_] ;
578  originalUpper_ = model_->solver()->getColUpper()[columnNumber_] ;
579  breakEven_ = breakEven;
580  assert (breakEven_>0.0&&breakEven_<1.0);
581}
582
583// Copy constructor
584CbcSimpleInteger::CbcSimpleInteger ( const CbcSimpleInteger & rhs)
585  :CbcObject(rhs)
586
587{
588  sequence_ = rhs.sequence_;
589  columnNumber_ = rhs.columnNumber_;
590  originalLower_ = rhs.originalLower_;
591  originalUpper_ = rhs.originalUpper_;
592  breakEven_ = rhs.breakEven_;
593}
594
595// Clone
596CbcObject *
597CbcSimpleInteger::clone() const
598{
599  return new CbcSimpleInteger(*this);
600}
601
602// Assignment operator
603CbcSimpleInteger & 
604CbcSimpleInteger::operator=( const CbcSimpleInteger& rhs)
605{
606  if (this!=&rhs) {
607    CbcObject::operator=(rhs);
608    columnNumber_ = model_->integerVariable()[sequence_];
609    breakEven_ = rhs.breakEven_;
610  }
611  return *this;
612}
613
614// Destructor
615CbcSimpleInteger::~CbcSimpleInteger ()
616{
617}
618
619// Infeasibility - large is 0.5
620double 
621CbcSimpleInteger::infeasibility(int & preferredWay) const
622{
623  OsiSolverInterface * solver = model_->solver();
624  const double * solution = model_->currentSolution();
625  const double * lower = solver->getColLower();
626  const double * upper = solver->getColUpper();
627  double value = solution[columnNumber_];
628  value = CoinMax(value, lower[columnNumber_]);
629  value = CoinMin(value, upper[columnNumber_]);
630  /*printf("%d %g %g %g %g\n",columnNumber_,value,lower[columnNumber_],
631    solution[columnNumber_],upper[columnNumber_]);*/
632  double nearest = floor(value+(1.0-breakEven_));
633  double integerTolerance = 
634    model_->getDblParam(CbcModel::CbcIntegerTolerance);
635  if (nearest>value) 
636    preferredWay=1;
637  else
638    preferredWay=-1;
639  double weight = fabs(value-nearest);
640  // normalize so weight is 0.5 at break even
641  if (nearest<value)
642    weight = (0.5/breakEven_)*weight;
643  else
644    weight = (0.5/(1.0-breakEven_))*weight;
645  if (fabs(value-nearest)<=integerTolerance) 
646    return 0.0;
647  else
648    return weight;
649}
650
651// This looks at solution and sets bounds to contain solution
652/** More precisely: it first forces the variable within the existing
653    bounds, and then tightens the bounds to fix the variable at the
654    nearest integer value.
655*/
656void 
657CbcSimpleInteger::feasibleRegion()
658{
659  OsiSolverInterface * solver = model_->solver();
660  const double * lower = solver->getColLower();
661  const double * upper = solver->getColUpper();
662  const double * solution = model_->currentSolution();
663  double value = solution[columnNumber_];
664  value = CoinMax(value, lower[columnNumber_]);
665  value = CoinMin(value, upper[columnNumber_]);
666
667  double nearest = floor(value+0.5);
668  //double integerTolerance =
669  //model_->getDblParam(CbcModel::CbcIntegerTolerance);
670  // Scaling may have moved it a bit
671  //assert (fabs(value-nearest)<=100.0*integerTolerance);
672  assert (fabs(value-nearest)<=0.01);
673  solver->setColLower(columnNumber_,nearest);
674  solver->setColUpper(columnNumber_,nearest);
675}
676/* Column number if single column object -1 otherwise,
677   so returns >= 0
678   Used by heuristics
679*/
680int 
681CbcSimpleInteger::columnNumber() const
682{
683  return columnNumber_;
684}
685
686// Creates a branching object
687CbcBranchingObject * 
688CbcSimpleInteger::createBranch(int way) const
689{
690  OsiSolverInterface * solver = model_->solver();
691  const double * solution = model_->currentSolution();
692  const double * lower = solver->getColLower();
693  const double * upper = solver->getColUpper();
694  double value = solution[columnNumber_];
695  value = CoinMax(value, lower[columnNumber_]);
696  value = CoinMin(value, upper[columnNumber_]);
697  double nearest = floor(value+0.5);
698  double integerTolerance = 
699    model_->getDblParam(CbcModel::CbcIntegerTolerance);
700  assert (upper[columnNumber_]>lower[columnNumber_]);
701  int hotstartStrategy=model_->getHotstartStrategy();
702  if (hotstartStrategy<=0) {
703    assert (fabs(value-nearest)>integerTolerance);
704  } else {
705    const double * bestSolution = model_->bestSolution();
706    double targetValue = bestSolution[columnNumber_];
707    if (way>0)
708      value = targetValue-0.1;
709    else
710      value = targetValue+0.1;
711  }
712  return new CbcIntegerBranchingObject(model_,sequence_,way,
713                                             value);
714}
715
716
717/* Given valid solution (i.e. satisfied) and reduced costs etc
718   returns a branching object which would give a new feasible
719   point in direction reduced cost says would be cheaper.
720   If no feasible point returns null
721*/
722CbcBranchingObject * 
723CbcSimpleInteger::preferredNewFeasible() const
724{
725  OsiSolverInterface * solver = model_->solver();
726  double value = model_->currentSolution()[columnNumber_];
727
728  double nearest = floor(value+0.5);
729  double integerTolerance = 
730    model_->getDblParam(CbcModel::CbcIntegerTolerance);
731  assert (fabs(value-nearest)<=integerTolerance);
732  double dj = solver->getObjSense()*solver->getReducedCost()[columnNumber_];
733  CbcIntegerBranchingObject * object = NULL;
734  if (dj>=0.0) {
735    // can we go down
736    if (nearest>originalLower_+0.5) {
737      // yes
738      object = new CbcIntegerBranchingObject(model_,sequence_,-1,
739                                             nearest-1.0,nearest-1.0);
740    }
741  } else {
742    // can we go up
743    if (nearest<originalUpper_-0.5) {
744      // yes
745      object = new CbcIntegerBranchingObject(model_,sequence_,-1,
746                                             nearest+1.0,nearest+1.0);
747    }
748  }
749  return object;
750}
751 
752/* Given valid solution (i.e. satisfied) and reduced costs etc
753   returns a branching object which would give a new feasible
754   point in direction opposite to one reduced cost says would be cheaper.
755   If no feasible point returns null
756*/
757CbcBranchingObject * 
758CbcSimpleInteger::notPreferredNewFeasible() const 
759{
760  OsiSolverInterface * solver = model_->solver();
761  double value = model_->currentSolution()[columnNumber_];
762
763  double nearest = floor(value+0.5);
764  double integerTolerance = 
765    model_->getDblParam(CbcModel::CbcIntegerTolerance);
766  assert (fabs(value-nearest)<=integerTolerance);
767  double dj = solver->getObjSense()*solver->getReducedCost()[columnNumber_];
768  CbcIntegerBranchingObject * object = NULL;
769  if (dj<=0.0) {
770    // can we go down
771    if (nearest>originalLower_+0.5) {
772      // yes
773      object = new CbcIntegerBranchingObject(model_,sequence_,-1,
774                                             nearest-1.0,nearest-1.0);
775    }
776  } else {
777    // can we go up
778    if (nearest<originalUpper_-0.5) {
779      // yes
780      object = new CbcIntegerBranchingObject(model_,sequence_,-1,
781                                             nearest+1.0,nearest+1.0);
782    }
783  }
784  return object;
785}
786 
787/*
788  Bounds may be tightened, so it may be good to be able to refresh the local
789  copy of the original bounds.
790 */
791void 
792CbcSimpleInteger::resetBounds()
793{
794  originalLower_ = model_->solver()->getColLower()[columnNumber_];
795  originalUpper_ = model_->solver()->getColUpper()[columnNumber_];
796}
797
798
799// Default Constructor
800CbcIntegerBranchingObject::CbcIntegerBranchingObject()
801  :CbcBranchingObject()
802{
803  down_[0] = 0.0;
804  down_[1] = 0.0;
805  up_[0] = 0.0;
806  up_[1] = 0.0;
807}
808
809// Useful constructor
810CbcIntegerBranchingObject::CbcIntegerBranchingObject (CbcModel * model, 
811                                                      int variable, int way , double value)
812  :CbcBranchingObject(model,variable,way,value)
813{
814  int iColumn = model_->integerVariable()[variable_];
815  down_[0] = model_->solver()->getColLower()[iColumn];
816  down_[1] = floor(value_);
817  up_[0] = ceil(value_);
818  up_[1] = model->getColUpper()[iColumn];
819}
820// Useful constructor for fixing
821CbcIntegerBranchingObject::CbcIntegerBranchingObject (CbcModel * model, 
822                                                      int variable, int way,
823                                                      double lowerValue, 
824                                                      double upperValue)
825  :CbcBranchingObject(model,variable,way,lowerValue)
826{
827  numberBranchesLeft_=1;
828  down_[0] = lowerValue;
829  down_[1] = upperValue;
830  up_[0] = lowerValue;
831  up_[1] = upperValue;
832}
833 
834
835// Copy constructor
836CbcIntegerBranchingObject::CbcIntegerBranchingObject ( const CbcIntegerBranchingObject & rhs) :CbcBranchingObject(rhs)
837{
838  down_[0] = rhs.down_[0];
839  down_[1] = rhs.down_[1];
840  up_[0] = rhs.up_[0];
841  up_[1] = rhs.up_[1];
842}
843
844// Assignment operator
845CbcIntegerBranchingObject & 
846CbcIntegerBranchingObject::operator=( const CbcIntegerBranchingObject& rhs)
847{
848  if (this != &rhs) {
849    CbcBranchingObject::operator=(rhs);
850    down_[0] = rhs.down_[0];
851    down_[1] = rhs.down_[1];
852    up_[0] = rhs.up_[0];
853    up_[1] = rhs.up_[1];
854  }
855  return *this;
856}
857CbcBranchingObject * 
858CbcIntegerBranchingObject::clone() const
859{ 
860  return (new CbcIntegerBranchingObject(*this));
861}
862
863
864// Destructor
865CbcIntegerBranchingObject::~CbcIntegerBranchingObject ()
866{
867}
868
869/*
870  Perform a branch by adjusting the bounds of the specified variable. Note
871  that each arm of the branch advances the object to the next arm by
872  advancing the value of way_.
873
874  Providing new values for the variable's lower and upper bounds for each
875  branching direction gives a little bit of additional flexibility and will
876  be easily extensible to multi-way branching.
877  Returns change in guessed objective on next branch
878*/
879double
880CbcIntegerBranchingObject::branch(bool normalBranch)
881{
882  if (model_->messageHandler()->logLevel()>2&&normalBranch)
883    print(normalBranch);
884  numberBranchesLeft_--;
885  int iColumn = model_->integerVariable()[variable_];
886  if (way_<0) {
887#ifdef CBC_DEBUG
888  { double olb,oub ;
889    olb = model_->solver()->getColLower()[iColumn] ;
890    oub = model_->solver()->getColUpper()[iColumn] ;
891    printf("branching down on var %d: [%g,%g] => [%g,%g]\n",
892           iColumn,olb,oub,down_[0],down_[1]) ; }
893#endif
894    model_->solver()->setColLower(iColumn,down_[0]);
895    model_->solver()->setColUpper(iColumn,down_[1]);
896    way_=1;
897  } else {
898#ifdef CBC_DEBUG
899  { double olb,oub ;
900    olb = model_->solver()->getColLower()[iColumn] ;
901    oub = model_->solver()->getColUpper()[iColumn] ;
902    printf("branching up on var %d: [%g,%g] => [%g,%g]\n",
903           iColumn,olb,oub,up_[0],up_[1]) ; }
904#endif
905    model_->solver()->setColLower(iColumn,up_[0]);
906    model_->solver()->setColUpper(iColumn,up_[1]);
907    way_=-1;      // Swap direction
908  }
909  return 0.0;
910}
911// Print what would happen 
912void
913CbcIntegerBranchingObject::print(bool normalBranch)
914{
915  int iColumn = model_->integerVariable()[variable_];
916  if (way_<0) {
917  { double olb,oub ;
918    olb = model_->solver()->getColLower()[iColumn] ;
919    oub = model_->solver()->getColUpper()[iColumn] ;
920    printf("CbcInteger would branch down on var %d: [%g,%g] => [%g,%g]\n",
921           iColumn,olb,oub,down_[0],down_[1]) ; }
922  } else {
923  { double olb,oub ;
924    olb = model_->solver()->getColLower()[iColumn] ;
925    oub = model_->solver()->getColUpper()[iColumn] ;
926    printf("CbcInteger would branch up on var %d: [%g,%g] => [%g,%g]\n",
927           iColumn,olb,oub,up_[0],up_[1]) ; }
928  }
929}
930
931
932/** Default Constructor
933
934  Equivalent to an unspecified binary variable.
935*/
936CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost ()
937  : CbcSimpleInteger(),
938    downPseudoCost_(1.0e-5),
939    upPseudoCost_(1.0e-5),
940    method_(0)
941{
942}
943
944/** Useful constructor
945
946  Loads actual upper & lower bounds for the specified variable.
947*/
948CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost (CbcModel * model, int sequence,
949                                    int iColumn, double breakEven)
950  : CbcSimpleInteger(model,sequence,iColumn,breakEven)
951{
952  const double * cost = model->getObjCoefficients();
953  double costValue = CoinMax(1.0e-5,fabs(cost[iColumn]));
954  // treat as if will cost what it says up
955  upPseudoCost_=costValue;
956  // and balance at breakeven
957  downPseudoCost_=((1.0-breakEven_)*upPseudoCost_)/breakEven_;
958  method_=0;
959}
960
961/** Useful constructor
962
963  Loads actual upper & lower bounds for the specified variable.
964*/
965CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost (CbcModel * model, int sequence,
966                                    int iColumn, double downPseudoCost,
967                                                        double upPseudoCost)
968  : CbcSimpleInteger(model,sequence,iColumn)
969{
970  downPseudoCost_ = downPseudoCost;
971  upPseudoCost_ = upPseudoCost;
972  breakEven_ = upPseudoCost_/(upPseudoCost_+downPseudoCost_);
973  method_=0;
974}
975
976// Copy constructor
977CbcSimpleIntegerPseudoCost::CbcSimpleIntegerPseudoCost ( const CbcSimpleIntegerPseudoCost & rhs)
978  :CbcSimpleInteger(rhs),
979   downPseudoCost_(rhs.downPseudoCost_),
980   upPseudoCost_(rhs.upPseudoCost_),
981   method_(rhs.method_)
982
983{
984}
985
986// Clone
987CbcObject *
988CbcSimpleIntegerPseudoCost::clone() const
989{
990  return new CbcSimpleIntegerPseudoCost(*this);
991}
992
993// Assignment operator
994CbcSimpleIntegerPseudoCost & 
995CbcSimpleIntegerPseudoCost::operator=( const CbcSimpleIntegerPseudoCost& rhs)
996{
997  if (this!=&rhs) {
998    CbcSimpleInteger::operator=(rhs);
999    downPseudoCost_=rhs.downPseudoCost_;
1000    upPseudoCost_=rhs.upPseudoCost_;
1001    method_=rhs.method_;
1002  }
1003  return *this;
1004}
1005
1006// Destructor
1007CbcSimpleIntegerPseudoCost::~CbcSimpleIntegerPseudoCost ()
1008{
1009}
1010// Creates a branching object
1011CbcBranchingObject * 
1012CbcSimpleIntegerPseudoCost::createBranch(int way) const
1013{
1014  OsiSolverInterface * solver = model_->solver();
1015  const double * solution = model_->currentSolution();
1016  const double * lower = solver->getColLower();
1017  const double * upper = solver->getColUpper();
1018  double value = solution[columnNumber_];
1019  value = CoinMax(value, lower[columnNumber_]);
1020  value = CoinMin(value, upper[columnNumber_]);
1021  double nearest = floor(value+0.5);
1022  double integerTolerance = 
1023    model_->getDblParam(CbcModel::CbcIntegerTolerance);
1024  assert (upper[columnNumber_]>lower[columnNumber_]);
1025  int hotstartStrategy=model_->getHotstartStrategy();
1026  if (hotstartStrategy<=0) {
1027    assert (fabs(value-nearest)>integerTolerance);
1028  } else {
1029    const double * bestSolution = model_->bestSolution();
1030    double targetValue = bestSolution[columnNumber_];
1031    if (way>0)
1032      value = targetValue-0.1;
1033    else
1034      value = targetValue+0.1;
1035  }
1036  CbcIntegerPseudoCostBranchingObject * newObject = 
1037    new CbcIntegerPseudoCostBranchingObject(model_,sequence_,way,
1038                                            value);
1039  double up =  upPseudoCost_*(ceil(value)-value);
1040  double down =  downPseudoCost_*(value-floor(value));
1041  double changeInGuessed=up-down;
1042  if (way>0)
1043    changeInGuessed = - changeInGuessed;
1044  changeInGuessed=CoinMax(0.0,changeInGuessed);
1045  //if (way>0)
1046  //changeInGuessed += 1.0e8; // bias to stay up
1047  newObject->setChangeInGuessed(changeInGuessed);
1048  return newObject;
1049}
1050// Infeasibility - large is 0.5
1051double 
1052CbcSimpleIntegerPseudoCost::infeasibility(int & preferredWay) const
1053{
1054  OsiSolverInterface * solver = model_->solver();
1055  const double * solution = model_->currentSolution();
1056  const double * lower = solver->getColLower();
1057  const double * upper = solver->getColUpper();
1058  if (upper[columnNumber_]==lower[columnNumber_]) {
1059    // fixed
1060    preferredWay=1;
1061    return 0.0;
1062  }
1063  double value = solution[columnNumber_];
1064  value = CoinMax(value, lower[columnNumber_]);
1065  value = CoinMin(value, upper[columnNumber_]);
1066  /*printf("%d %g %g %g %g\n",columnNumber_,value,lower[columnNumber_],
1067    solution[columnNumber_],upper[columnNumber_]);*/
1068  double nearest = floor(value+0.5);
1069  double integerTolerance = 
1070    model_->getDblParam(CbcModel::CbcIntegerTolerance);
1071  double below = floor(value+integerTolerance);
1072  double above = below+1.0;
1073  if (above>upper[columnNumber_]) {
1074    above=below;
1075    below = above -1;
1076  }
1077  double downCost = CoinMax((value-below)*downPseudoCost_,0.0);
1078  double upCost = CoinMax((above-value)*upPseudoCost_,0.0);
1079  if (downCost>=upCost)
1080    preferredWay=1;
1081  else
1082    preferredWay=-1;
1083  if (fabs(value-nearest)<=integerTolerance) {
1084    return 0.0;
1085  } else {
1086    // can't get at model so 1,2 don't make sense
1087    assert(method_<1||method_>2);
1088    if (!method_)
1089      return CoinMin(downCost,upCost);
1090    else
1091      return CoinMax(downCost,upCost);
1092  }
1093}
1094
1095// Return "up" estimate
1096double 
1097CbcSimpleIntegerPseudoCost::upEstimate() const
1098{
1099  OsiSolverInterface * solver = model_->solver();
1100  const double * solution = model_->currentSolution();
1101  const double * lower = solver->getColLower();
1102  const double * upper = solver->getColUpper();
1103  double value = solution[columnNumber_];
1104  value = CoinMax(value, lower[columnNumber_]);
1105  value = CoinMin(value, upper[columnNumber_]);
1106  if (upper[columnNumber_]==lower[columnNumber_]) {
1107    // fixed
1108    return 0.0;
1109  }
1110  double integerTolerance = 
1111    model_->getDblParam(CbcModel::CbcIntegerTolerance);
1112  double below = floor(value+integerTolerance);
1113  double above = below+1.0;
1114  if (above>upper[columnNumber_]) {
1115    above=below;
1116    below = above -1;
1117  }
1118  double upCost = CoinMax((above-value)*upPseudoCost_,0.0);
1119  return upCost;
1120}
1121// Return "down" estimate
1122double 
1123CbcSimpleIntegerPseudoCost::downEstimate() const
1124{
1125  OsiSolverInterface * solver = model_->solver();
1126  const double * solution = model_->currentSolution();
1127  const double * lower = solver->getColLower();
1128  const double * upper = solver->getColUpper();
1129  double value = solution[columnNumber_];
1130  value = CoinMax(value, lower[columnNumber_]);
1131  value = CoinMin(value, upper[columnNumber_]);
1132  if (upper[columnNumber_]==lower[columnNumber_]) {
1133    // fixed
1134    return 0.0;
1135  }
1136  double integerTolerance = 
1137    model_->getDblParam(CbcModel::CbcIntegerTolerance);
1138  double below = floor(value+integerTolerance);
1139  double above = below+1.0;
1140  if (above>upper[columnNumber_]) {
1141    above=below;
1142    below = above -1;
1143  }
1144  double downCost = CoinMax((value-below)*downPseudoCost_,0.0);
1145  return downCost;
1146}
1147
1148// Default Constructor
1149CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject()
1150  :CbcIntegerBranchingObject()
1151{
1152  changeInGuessed_=1.0e-5;
1153}
1154
1155// Useful constructor
1156CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject (CbcModel * model, 
1157                                                      int variable, int way , double value)
1158  :CbcIntegerBranchingObject(model,variable,way,value)
1159{
1160}
1161// Useful constructor for fixing
1162CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject (CbcModel * model, 
1163                                                      int variable, int way,
1164                                                      double lowerValue, 
1165                                                      double upperValue)
1166  :CbcIntegerBranchingObject(model,variable,way,lowerValue)
1167{
1168  changeInGuessed_=1.0e100;
1169}
1170 
1171
1172// Copy constructor
1173CbcIntegerPseudoCostBranchingObject::CbcIntegerPseudoCostBranchingObject ( 
1174                                 const CbcIntegerPseudoCostBranchingObject & rhs)
1175  :CbcIntegerBranchingObject(rhs)
1176{
1177  changeInGuessed_ = rhs.changeInGuessed_;
1178}
1179
1180// Assignment operator
1181CbcIntegerPseudoCostBranchingObject & 
1182CbcIntegerPseudoCostBranchingObject::operator=( const CbcIntegerPseudoCostBranchingObject& rhs)
1183{
1184  if (this != &rhs) {
1185    CbcIntegerBranchingObject::operator=(rhs);
1186    changeInGuessed_ = rhs.changeInGuessed_;
1187  }
1188  return *this;
1189}
1190CbcBranchingObject * 
1191CbcIntegerPseudoCostBranchingObject::clone() const
1192{ 
1193  return (new CbcIntegerPseudoCostBranchingObject(*this));
1194}
1195
1196
1197// Destructor
1198CbcIntegerPseudoCostBranchingObject::~CbcIntegerPseudoCostBranchingObject ()
1199{
1200}
1201
1202/*
1203  Perform a branch by adjusting the bounds of the specified variable. Note
1204  that each arm of the branch advances the object to the next arm by
1205  advancing the value of way_.
1206
1207  Providing new values for the variable's lower and upper bounds for each
1208  branching direction gives a little bit of additional flexibility and will
1209  be easily extensible to multi-way branching.
1210  Returns change in guessed objective on next branch
1211*/
1212double
1213CbcIntegerPseudoCostBranchingObject::branch(bool normalBranch)
1214{
1215  CbcIntegerBranchingObject::branch(normalBranch);
1216  return changeInGuessed_;
1217}
1218
1219
1220// Default Constructor
1221CbcCliqueBranchingObject::CbcCliqueBranchingObject()
1222  :CbcBranchingObject()
1223{
1224  clique_ = NULL;
1225  downMask_[0]=0;
1226  downMask_[1]=0;
1227  upMask_[0]=0;
1228  upMask_[1]=0;
1229}
1230
1231// Useful constructor
1232CbcCliqueBranchingObject::CbcCliqueBranchingObject (CbcModel * model,
1233                                                    const CbcClique * clique,
1234                                                    int way ,
1235                                                    int numberOnDownSide, const int * down,
1236                                                    int numberOnUpSide, const int * up)
1237  :CbcBranchingObject(model,clique->id(),way,0.5)
1238{
1239  clique_ = clique;
1240  downMask_[0]=0;
1241  downMask_[1]=0;
1242  upMask_[0]=0;
1243  upMask_[1]=0;
1244  int i;
1245  for (i=0;i<numberOnDownSide;i++) {
1246    int sequence = down[i];
1247    int iWord = sequence>>5;
1248    int iBit = sequence - 32*iWord;
1249    unsigned int k = 1<<iBit;
1250    downMask_[iWord] |= k;
1251  }
1252  for (i=0;i<numberOnUpSide;i++) {
1253    int sequence = up[i];
1254    int iWord = sequence>>5;
1255    int iBit = sequence - 32*iWord;
1256    unsigned int k = 1<<iBit;
1257    upMask_[iWord] |= k;
1258  }
1259}
1260
1261// Copy constructor
1262CbcCliqueBranchingObject::CbcCliqueBranchingObject ( const CbcCliqueBranchingObject & rhs) :CbcBranchingObject(rhs)
1263{
1264  clique_=rhs.clique_;
1265  downMask_[0]=rhs.downMask_[0];
1266  downMask_[1]=rhs.downMask_[1];
1267  upMask_[0]=rhs.upMask_[0];
1268  upMask_[1]=rhs.upMask_[1];
1269}
1270
1271// Assignment operator
1272CbcCliqueBranchingObject & 
1273CbcCliqueBranchingObject::operator=( const CbcCliqueBranchingObject& rhs)
1274{
1275  if (this != &rhs) {
1276    CbcBranchingObject::operator=(rhs);
1277    clique_=rhs.clique_;
1278    downMask_[0]=rhs.downMask_[0];
1279    downMask_[1]=rhs.downMask_[1];
1280    upMask_[0]=rhs.upMask_[0];
1281    upMask_[1]=rhs.upMask_[1];
1282  }
1283  return *this;
1284}
1285CbcBranchingObject * 
1286CbcCliqueBranchingObject::clone() const
1287{ 
1288  return (new CbcCliqueBranchingObject(*this));
1289}
1290
1291
1292// Destructor
1293CbcCliqueBranchingObject::~CbcCliqueBranchingObject ()
1294{
1295}
1296double
1297CbcCliqueBranchingObject::branch(bool normalBranch)
1298{
1299  if (model_->messageHandler()->logLevel()>2&&normalBranch)
1300    print(normalBranch);
1301  numberBranchesLeft_--;
1302  int iWord;
1303  int numberMembers = clique_->numberMembers();
1304  const int * which = clique_->members();
1305  const int * integerVariables = model_->integerVariable();
1306  int numberWords=(numberMembers+31)>>5;
1307  // *** for way - up means fix all those in down section
1308  if (way_<0) {
1309#ifdef FULL_PRINT
1310    printf("Down Fix ");
1311#endif
1312    for (iWord=0;iWord<numberWords;iWord++) {
1313      int i;
1314      for (i=0;i<32;i++) {
1315        unsigned int k = 1<<i;
1316        if ((upMask_[iWord]&k)!=0) {
1317          int iColumn = which[i+32*iWord];
1318#ifdef FULL_PRINT
1319          printf("%d ",i+32*iWord);
1320#endif
1321          // fix weak way
1322          if (clique_->type(i+32*iWord))
1323            model_->solver()->setColUpper(integerVariables[iColumn],0.0);
1324          else
1325            model_->solver()->setColLower(integerVariables[iColumn],1.0);
1326        }
1327      }
1328    }
1329    way_=1;       // Swap direction
1330  } else {
1331#ifdef FULL_PRINT
1332    printf("Up Fix ");
1333#endif
1334    for (iWord=0;iWord<numberWords;iWord++) {
1335      int i;
1336      for (i=0;i<32;i++) {
1337        unsigned int k = 1<<i;
1338        if ((downMask_[iWord]&k)!=0) {
1339          int iColumn = which[i+32*iWord];
1340#ifdef FULL_PRINT
1341          printf("%d ",i+32*iWord);
1342#endif
1343          // fix weak way
1344          if (clique_->type(i+32*iWord))
1345            model_->solver()->setColUpper(integerVariables[iColumn],0.0);
1346          else
1347            model_->solver()->setColLower(integerVariables[iColumn],1.0);
1348        }
1349      }
1350    }
1351    way_=-1;      // Swap direction
1352  }
1353#ifdef FULL_PRINT
1354  printf("\n");
1355#endif
1356  return 0.0;
1357}
1358// Print what would happen 
1359void
1360CbcCliqueBranchingObject::print(bool normalBranch)
1361{
1362  int iWord;
1363  int numberMembers = clique_->numberMembers();
1364  const int * which = clique_->members();
1365  const int * integerVariables = model_->integerVariable();
1366  int numberWords=(numberMembers+31)>>5;
1367  // *** for way - up means fix all those in down section
1368  if (way_<0) {
1369    printf("Clique - Down Fix ");
1370    for (iWord=0;iWord<numberWords;iWord++) {
1371      int i;
1372      for (i=0;i<32;i++) {
1373        unsigned int k = 1<<i;
1374        if ((upMask_[iWord]&k)!=0) {
1375          int iColumn = which[i+32*iWord];
1376          printf("%d ",integerVariables[iColumn]);
1377        }
1378      }
1379    }
1380  } else {
1381    printf("Clique - Up Fix ");
1382    for (iWord=0;iWord<numberWords;iWord++) {
1383      int i;
1384      for (i=0;i<32;i++) {
1385        unsigned int k = 1<<i;
1386        if ((downMask_[iWord]&k)!=0) {
1387          int iColumn = which[i+32*iWord];
1388          printf("%d ",integerVariables[iColumn]);
1389        }
1390      }
1391    }
1392  }
1393  printf("\n");
1394}
1395 
1396// Default Constructor
1397CbcLongCliqueBranchingObject::CbcLongCliqueBranchingObject()
1398  :CbcBranchingObject()
1399{
1400  clique_=NULL;
1401  downMask_=NULL;
1402  upMask_=NULL;
1403}
1404
1405// Useful constructor
1406CbcLongCliqueBranchingObject::CbcLongCliqueBranchingObject (CbcModel * model,
1407                                                            const CbcClique * clique, 
1408                                                            int way ,
1409                                                    int numberOnDownSide, const int * down,
1410                                                    int numberOnUpSide, const int * up)
1411  :CbcBranchingObject(model,clique->id(),way,0.5)
1412{
1413  clique_ = clique;
1414  int numberMembers = clique_->numberMembers();
1415  int numberWords=(numberMembers+31)>>5;
1416  downMask_ = new unsigned int [numberWords];
1417  upMask_ = new unsigned int [numberWords];
1418  memset(downMask_,0,numberWords*sizeof(unsigned int));
1419  memset(upMask_,0,numberWords*sizeof(unsigned int));
1420  int i;
1421  for (i=0;i<numberOnDownSide;i++) {
1422    int sequence = down[i];
1423    int iWord = sequence>>5;
1424    int iBit = sequence - 32*iWord;
1425    unsigned int k = 1<<iBit;
1426    downMask_[iWord] |= k;
1427  }
1428  for (i=0;i<numberOnUpSide;i++) {
1429    int sequence = up[i];
1430    int iWord = sequence>>5;
1431    int iBit = sequence - 32*iWord;
1432    unsigned int k = 1<<iBit;
1433    upMask_[iWord] |= k;
1434  }
1435}
1436
1437// Copy constructor
1438CbcLongCliqueBranchingObject::CbcLongCliqueBranchingObject ( const CbcLongCliqueBranchingObject & rhs) :CbcBranchingObject(rhs)
1439{
1440  clique_=rhs.clique_;
1441  if (rhs.downMask_) {
1442    int numberMembers = clique_->numberMembers();
1443    int numberWords=(numberMembers+31)>>5;
1444    downMask_ = new unsigned int [numberWords];
1445    memcpy(downMask_,rhs.downMask_,numberWords*sizeof(unsigned int));
1446    upMask_ = new unsigned int [numberWords];
1447    memcpy(upMask_,rhs.upMask_,numberWords*sizeof(unsigned int));
1448  } else {
1449    downMask_=NULL;
1450    upMask_=NULL;
1451  }   
1452}
1453
1454// Assignment operator
1455CbcLongCliqueBranchingObject & 
1456CbcLongCliqueBranchingObject::operator=( const CbcLongCliqueBranchingObject& rhs)
1457{
1458  if (this != &rhs) {
1459    CbcBranchingObject::operator=(rhs);
1460    clique_=rhs.clique_;
1461    delete [] downMask_;
1462    delete [] upMask_;
1463    if (rhs.downMask_) {
1464      int numberMembers = clique_->numberMembers();
1465      int numberWords=(numberMembers+31)>>5;
1466      downMask_ = new unsigned int [numberWords];
1467      memcpy(downMask_,rhs.downMask_,numberWords*sizeof(unsigned int));
1468      upMask_ = new unsigned int [numberWords];
1469      memcpy(upMask_,rhs.upMask_,numberWords*sizeof(unsigned int));
1470    } else {
1471      downMask_=NULL;
1472      upMask_=NULL;
1473    }   
1474  }
1475  return *this;
1476}
1477CbcBranchingObject * 
1478CbcLongCliqueBranchingObject::clone() const
1479{ 
1480  return (new CbcLongCliqueBranchingObject(*this));
1481}
1482
1483
1484// Destructor
1485CbcLongCliqueBranchingObject::~CbcLongCliqueBranchingObject ()
1486{
1487  delete [] downMask_;
1488  delete [] upMask_;
1489}
1490double
1491CbcLongCliqueBranchingObject::branch(bool normalBranch)
1492{
1493  if (model_->messageHandler()->logLevel()>2&&normalBranch)
1494    print(normalBranch);
1495  numberBranchesLeft_--;
1496  int iWord;
1497  int numberMembers = clique_->numberMembers();
1498  const int * which = clique_->members();
1499  const int * integerVariables = model_->integerVariable();
1500  int numberWords=(numberMembers+31)>>5;
1501  // *** for way - up means fix all those in down section
1502  if (way_<0) {
1503#ifdef FULL_PRINT
1504    printf("Down Fix ");
1505#endif
1506    for (iWord=0;iWord<numberWords;iWord++) {
1507      int i;
1508      for (i=0;i<32;i++) {
1509        unsigned int k = 1<<i;
1510        if ((upMask_[iWord]&k)!=0) {
1511          int iColumn = which[i+32*iWord];
1512#ifdef FULL_PRINT
1513          printf("%d ",i+32*iWord);
1514#endif
1515          // fix weak way
1516          if (clique_->type(i+32*iWord))
1517            model_->solver()->setColUpper(integerVariables[iColumn],0.0);
1518          else
1519            model_->solver()->setColLower(integerVariables[iColumn],1.0);
1520        }
1521      }
1522    }
1523    way_=1;       // Swap direction
1524  } else {
1525#ifdef FULL_PRINT
1526    printf("Up Fix ");
1527#endif
1528    for (iWord=0;iWord<numberWords;iWord++) {
1529      int i;
1530      for (i=0;i<32;i++) {
1531        unsigned int k = 1<<i;
1532        if ((downMask_[iWord]&k)!=0) {
1533          int iColumn = which[i+32*iWord];
1534#ifdef FULL_PRINT
1535          printf("%d ",i+32*iWord);
1536#endif
1537          // fix weak way
1538          if (clique_->type(i+32*iWord))
1539            model_->solver()->setColUpper(integerVariables[iColumn],0.0);
1540          else
1541            model_->solver()->setColLower(integerVariables[iColumn],1.0);
1542        }
1543      }
1544    }
1545    way_=-1;      // Swap direction
1546  }
1547#ifdef FULL_PRINT
1548  printf("\n");
1549#endif
1550  return 0.0;
1551}
1552void
1553CbcLongCliqueBranchingObject::print(bool normalBranch)
1554{
1555  int iWord;
1556  int numberMembers = clique_->numberMembers();
1557  const int * which = clique_->members();
1558  const int * integerVariables = model_->integerVariable();
1559  int numberWords=(numberMembers+31)>>5;
1560  // *** for way - up means fix all those in down section
1561  if (way_<0) {
1562    printf("Clique - Down Fix ");
1563    for (iWord=0;iWord<numberWords;iWord++) {
1564      int i;
1565      for (i=0;i<32;i++) {
1566        unsigned int k = 1<<i;
1567        if ((upMask_[iWord]&k)!=0) {
1568          int iColumn = which[i+32*iWord];
1569          printf("%d ",integerVariables[iColumn]);
1570        }
1571      }
1572    }
1573  } else {
1574    printf("Clique - Up Fix ");
1575    for (iWord=0;iWord<numberWords;iWord++) {
1576      int i;
1577      for (i=0;i<32;i++) {
1578        unsigned int k = 1<<i;
1579        if ((downMask_[iWord]&k)!=0) {
1580          int iColumn = which[i+32*iWord];
1581          printf("%d ",integerVariables[iColumn]);
1582        }
1583      }
1584    }
1585  }
1586  printf("\n");
1587}
1588// Default Constructor
1589CbcSOSBranchingObject::CbcSOSBranchingObject()
1590  :CbcBranchingObject()
1591{
1592  set_ = NULL;
1593  separator_=0.0;
1594}
1595
1596// Useful constructor
1597CbcSOSBranchingObject::CbcSOSBranchingObject (CbcModel * model,
1598                                              const CbcSOS * set,
1599                                              int way ,
1600                                              double separator)
1601  :CbcBranchingObject(model,set->id(),way,0.5)
1602{
1603  set_ = set;
1604  separator_ = separator;
1605}
1606
1607// Copy constructor
1608CbcSOSBranchingObject::CbcSOSBranchingObject ( const CbcSOSBranchingObject & rhs) :CbcBranchingObject(rhs)
1609{
1610  set_=rhs.set_;
1611  separator_ = rhs.separator_;
1612}
1613
1614// Assignment operator
1615CbcSOSBranchingObject & 
1616CbcSOSBranchingObject::operator=( const CbcSOSBranchingObject& rhs)
1617{
1618  if (this != &rhs) {
1619    CbcBranchingObject::operator=(rhs);
1620    set_=rhs.set_;
1621    separator_ = rhs.separator_;
1622  }
1623  return *this;
1624}
1625CbcBranchingObject * 
1626CbcSOSBranchingObject::clone() const
1627{ 
1628  return (new CbcSOSBranchingObject(*this));
1629}
1630
1631
1632// Destructor
1633CbcSOSBranchingObject::~CbcSOSBranchingObject ()
1634{
1635}
1636double
1637CbcSOSBranchingObject::branch(bool normalBranch)
1638{
1639  if (model_->messageHandler()->logLevel()>2&&normalBranch)
1640    print(normalBranch);
1641  numberBranchesLeft_--;
1642  int numberMembers = set_->numberMembers();
1643  const int * which = set_->members();
1644  const double * weights = set_->weights();
1645  OsiSolverInterface * solver = model_->solver();
1646  //const double * lower = solver->getColLower();
1647  //const double * upper = solver->getColUpper();
1648  // *** for way - up means fix all those in down section
1649  if (way_<0) {
1650    int i;
1651    for ( i=0;i<numberMembers;i++) {
1652      if (weights[i] > separator_)
1653        break;
1654    }
1655    assert (i<numberMembers);
1656    for (;i<numberMembers;i++) 
1657      solver->setColUpper(which[i],0.0);
1658    way_=1;       // Swap direction
1659  } else {
1660    int i;
1661    for ( i=0;i<numberMembers;i++) {
1662      if (weights[i] >= separator_)
1663        break;
1664      else
1665        solver->setColUpper(which[i],0.0);
1666    }
1667    assert (i<numberMembers);
1668    way_=-1;      // Swap direction
1669  }
1670  return 0.0;
1671}
1672// Print what would happen 
1673void
1674CbcSOSBranchingObject::print(bool normalBranch)
1675{
1676  int numberMembers = set_->numberMembers();
1677  const int * which = set_->members();
1678  const double * weights = set_->weights();
1679  OsiSolverInterface * solver = model_->solver();
1680  //const double * lower = solver->getColLower();
1681  const double * upper = solver->getColUpper();
1682  int first=numberMembers;
1683  int last=-1;
1684  int numberFixed=0;
1685  int numberOther=0;
1686  int i;
1687  for ( i=0;i<numberMembers;i++) {
1688    double bound = upper[which[i]];
1689    if (bound) {
1690      first = CoinMin(first,i);
1691      last = CoinMax(last,i);
1692    }
1693  }
1694  // *** for way - up means fix all those in down section
1695  if (way_<0) {
1696    printf("SOS Down");
1697    for ( i=0;i<numberMembers;i++) {
1698      double bound = upper[which[i]];
1699      if (weights[i] > separator_)
1700        break;
1701      else if (bound)
1702        numberOther++;
1703    }
1704    assert (i<numberMembers);
1705    for (;i<numberMembers;i++) {
1706      double bound = upper[which[i]];
1707      if (bound)
1708        numberFixed++;
1709    }
1710  } else {
1711    printf("SOS Up");
1712    for ( i=0;i<numberMembers;i++) {
1713      double bound = upper[which[i]];
1714      if (weights[i] >= separator_)
1715        break;
1716      else if (bound)
1717        numberFixed++;
1718    }
1719    assert (i<numberMembers);
1720    for (;i<numberMembers;i++) {
1721      double bound = upper[which[i]];
1722      if (bound)
1723        numberOther++;
1724    }
1725  }
1726  printf(" - at %g, free range %d (%g) => %d (%g), %d would be fixed, %d other way\n",
1727         separator_,which[first],weights[first],which[last],weights[last],numberFixed,numberOther);
1728}
1729 
1730// Default Constructor
1731CbcBranchDefaultDecision::CbcBranchDefaultDecision()
1732  :CbcBranchDecision()
1733{
1734}
1735
1736// Copy constructor
1737CbcBranchDefaultDecision::CbcBranchDefaultDecision (
1738                                    const CbcBranchDefaultDecision & rhs)
1739  :CbcBranchDecision()
1740{
1741}
1742
1743CbcBranchDefaultDecision::~CbcBranchDefaultDecision()
1744{
1745}
1746
1747// Clone
1748CbcBranchDecision * 
1749CbcBranchDefaultDecision::clone() const
1750{
1751  return new CbcBranchDefaultDecision(*this);
1752}
1753
1754// Initialize i.e. before start of choosing at a node
1755void 
1756CbcBranchDefaultDecision::initialize(CbcModel * model)
1757{
1758}
1759
1760
1761/*
1762  Simple default decision algorithm. Compare based on infeasibility (numInfUp,
1763  numInfDn) until a solution is found by search, then switch to change in
1764  objective (changeUp, changeDn). Note that bestSoFar is remembered in
1765  bestObject_, so the parameter bestSoFar is unused.
1766*/
1767
1768int
1769CbcBranchDefaultDecision::betterBranch(CbcBranchingObject * thisOne,
1770                            CbcBranchingObject * bestSoFar,
1771                            double changeUp, int numInfUp,
1772                            double changeDn, int numInfDn)
1773{
1774  printf("Now obsolete CbcBranchDefaultDecision::betterBranch\n");
1775  abort();
1776  return 0;
1777}
1778
1779/* Compare N branching objects. Return index of best
1780   and sets way of branching in chosen object.
1781   
1782   This routine is used only after strong branching.
1783   This is reccommended version as it can be more sophisticated
1784*/
1785
1786int
1787CbcBranchDefaultDecision::bestBranch (CbcBranchingObject ** objects, int numberObjects,
1788                                   int numberUnsatisfied,
1789                                   double * changeUp, int * numberInfeasibilitiesUp,
1790                                   double * changeDown, int * numberInfeasibilitiesDown,
1791                                   double objectiveValue) 
1792{
1793
1794  int bestWay=0;
1795  int whichObject = -1;
1796  if (numberObjects) {
1797    CbcModel * model = objects[0]->model();
1798    // at continuous
1799    //double continuousObjective = model->getContinuousObjective();
1800    //int continuousInfeasibilities = model->getContinuousInfeasibilities();
1801   
1802    // average cost to get rid of infeasibility
1803    //double averageCostPerInfeasibility =
1804    //(objectiveValue-continuousObjective)/
1805    //(double) (abs(continuousInfeasibilities-numberUnsatisfied)+1);
1806    /* beforeSolution is :
1807       0 - before any solution
1808       n - n heuristic solutions but no branched one
1809       -1 - branched solution found
1810    */
1811    int numberSolutions = model->getSolutionCount();
1812    double cutoff = model->getCutoff();
1813    int method=0;
1814    int i;
1815    if (numberSolutions) {
1816      int numberHeuristic = model->getNumberHeuristicSolutions();
1817      if (numberHeuristic<numberSolutions) {
1818        method = 1;
1819      } else {
1820        method = 2;
1821        // look further
1822        for ( i = 0 ; i < numberObjects ; i++) {
1823          int numberNext = numberInfeasibilitiesUp[i];
1824         
1825          if (numberNext<numberUnsatisfied) {
1826            int numberUp = numberUnsatisfied - numberInfeasibilitiesUp[i];
1827            double perUnsatisfied = changeUp[i]/(double) numberUp;
1828            double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
1829            if (estimatedObjective<cutoff) 
1830              method=3;
1831          }
1832          numberNext = numberInfeasibilitiesDown[i];
1833          if (numberNext<numberUnsatisfied) {
1834            int numberDown = numberUnsatisfied - numberInfeasibilitiesDown[i];
1835            double perUnsatisfied = changeDown[i]/(double) numberDown;
1836            double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
1837            if (estimatedObjective<cutoff) 
1838              method=3;
1839          }
1840        }
1841      }
1842      method=2;
1843    } else {
1844      method = 0;
1845    }
1846    // Uncomment next to force method 4
1847    //method=4;
1848    /* Methods :
1849       0 - fewest infeasibilities
1850       1 - largest min change in objective
1851       2 - as 1 but use sum of changes if min close
1852       3 - predicted best solution
1853       4 - take cheapest up branch if infeasibilities same
1854    */
1855    int bestNumber=INT_MAX;
1856    double bestCriterion=-1.0e50;
1857    double alternativeCriterion = -1.0;
1858    double bestEstimate = 1.0e100;
1859    switch (method) {
1860    case 0:
1861      // could add in depth as well
1862      for ( i = 0 ; i < numberObjects ; i++) {
1863        int thisNumber = min(numberInfeasibilitiesUp[i],numberInfeasibilitiesDown[i]);
1864        if (thisNumber<=bestNumber) {
1865          int betterWay=0;
1866          if (numberInfeasibilitiesUp[i]<numberInfeasibilitiesDown[i]) {
1867            if (numberInfeasibilitiesUp[i]<bestNumber) {
1868              betterWay = 1;
1869            } else {
1870              if (changeUp[i]<bestCriterion)
1871                betterWay=1;
1872            }
1873          } else if (numberInfeasibilitiesUp[i]>numberInfeasibilitiesDown[i]) {
1874            if (numberInfeasibilitiesDown[i]<bestNumber) {
1875              betterWay = -1;
1876            } else {
1877              if (changeDown[i]<bestCriterion)
1878                betterWay=-1;
1879            }
1880          } else {
1881            // up and down have same number
1882            bool better=false;
1883            if (numberInfeasibilitiesUp[i]<bestNumber) {
1884              better=true;
1885            } else if (numberInfeasibilitiesUp[i]==bestNumber) {
1886              if (min(changeUp[i],changeDown[i])<bestCriterion)
1887                better=true;;
1888            }
1889            if (better) {
1890              // see which way
1891              if (changeUp[i]<=changeDown[i])
1892                betterWay=1;
1893              else
1894                betterWay=-1;
1895            }
1896          }
1897          if (betterWay) {
1898            bestCriterion = min(changeUp[i],changeDown[i]);
1899            bestNumber = thisNumber;
1900            whichObject = i;
1901            bestWay = betterWay;
1902          }
1903        }
1904      }
1905      break;
1906    case 1:
1907      for ( i = 0 ; i < numberObjects ; i++) {
1908        int betterWay=0;
1909        if (changeUp[i]<=changeDown[i]) {
1910          if (changeUp[i]>bestCriterion)
1911            betterWay=1;
1912        } else {
1913          if (changeDown[i]>bestCriterion)
1914            betterWay=-1;
1915        }
1916        if (betterWay) {
1917          bestCriterion = min(changeUp[i],changeDown[i]);
1918          whichObject = i;
1919          bestWay = betterWay;
1920        }
1921      }
1922      break;
1923    case 2:
1924      for ( i = 0 ; i < numberObjects ; i++) {
1925        double change = min(changeUp[i],changeDown[i]);
1926        double sum = changeUp[i] + changeDown[i];
1927        bool take=false;
1928        if (change>1.1*bestCriterion) 
1929          take=true;
1930        else if (change>0.9*bestCriterion&&sum+change>bestCriterion+alternativeCriterion) 
1931          take=true;
1932        if (take) {
1933          if (changeUp[i]<=changeDown[i]) {
1934            if (changeUp[i]>bestCriterion)
1935              bestWay=1;
1936          } else {
1937            if (changeDown[i]>bestCriterion)
1938              bestWay=-1;
1939          }
1940          bestCriterion = change;
1941          alternativeCriterion = sum;
1942          whichObject = i;
1943        }
1944      }
1945      break;
1946    case 3:
1947      for ( i = 0 ; i < numberObjects ; i++) {
1948        int numberNext = numberInfeasibilitiesUp[i];
1949       
1950        if (numberNext<numberUnsatisfied) {
1951          int numberUp = numberUnsatisfied - numberInfeasibilitiesUp[i];
1952          double perUnsatisfied = changeUp[i]/(double) numberUp;
1953          double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
1954          if (estimatedObjective<bestEstimate) {
1955            bestEstimate = estimatedObjective;
1956            bestWay=1;
1957            whichObject=i;
1958          }
1959        }
1960        numberNext = numberInfeasibilitiesDown[i];
1961        if (numberNext<numberUnsatisfied) {
1962          int numberDown = numberUnsatisfied - numberInfeasibilitiesDown[i];
1963          double perUnsatisfied = changeDown[i]/(double) numberDown;
1964          double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
1965          if (estimatedObjective<bestEstimate) {
1966            bestEstimate = estimatedObjective;
1967            bestWay=-1;
1968            whichObject=i;
1969          }
1970        }
1971      }
1972      break;
1973    case 4:
1974      // if number infeas same then cheapest up
1975      // first get best number or when going down
1976      // now choose smallest change up amongst equal number infeas
1977      for ( i = 0 ; i < numberObjects ; i++) {
1978        int thisNumber = min(numberInfeasibilitiesUp[i],numberInfeasibilitiesDown[i]);
1979        if (thisNumber<=bestNumber) {
1980          int betterWay=0;
1981          if (numberInfeasibilitiesUp[i]<numberInfeasibilitiesDown[i]) {
1982            if (numberInfeasibilitiesUp[i]<bestNumber) {
1983              betterWay = 1;
1984            } else {
1985              if (changeUp[i]<bestCriterion)
1986                betterWay=1;
1987            }
1988          } else if (numberInfeasibilitiesUp[i]>numberInfeasibilitiesDown[i]) {
1989            if (numberInfeasibilitiesDown[i]<bestNumber) {
1990              betterWay = -1;
1991            } else {
1992              if (changeDown[i]<bestCriterion)
1993                betterWay=-1;
1994            }
1995          } else {
1996            // up and down have same number
1997            bool better=false;
1998            if (numberInfeasibilitiesUp[i]<bestNumber) {
1999              better=true;
2000            } else if (numberInfeasibilitiesUp[i]==bestNumber) {
2001              if (min(changeUp[i],changeDown[i])<bestCriterion)
2002                better=true;;
2003            }
2004            if (better) {
2005              // see which way
2006              if (changeUp[i]<=changeDown[i])
2007                betterWay=1;
2008              else
2009                betterWay=-1;
2010            }
2011          }
2012          if (betterWay) {
2013            bestCriterion = min(changeUp[i],changeDown[i]);
2014            bestNumber = thisNumber;
2015            whichObject = i;
2016            bestWay = betterWay;
2017          }
2018        }
2019      }
2020      bestCriterion=1.0e50;
2021      for ( i = 0 ; i < numberObjects ; i++) {
2022        int thisNumber = numberInfeasibilitiesUp[i];
2023        if (thisNumber==bestNumber&&changeUp) {
2024          if (changeUp[i]<bestCriterion) {
2025            bestCriterion = changeUp[i];
2026            whichObject = i;
2027            bestWay = 1;
2028          }
2029        }
2030      }
2031      break;
2032    }
2033    // set way in best
2034    if (whichObject>=0)
2035      objects[whichObject]->way(bestWay);
2036  }
2037  return whichObject;
2038}
2039
2040// Default Constructor
2041CbcFollowOn::CbcFollowOn ()
2042  : CbcObject(),
2043    rhs_(NULL)
2044{
2045}
2046
2047// Useful constructor
2048CbcFollowOn::CbcFollowOn (CbcModel * model)
2049  : CbcObject(model)
2050{
2051  assert (model);
2052  OsiSolverInterface * solver = model_->solver();
2053  matrix_ = *solver->getMatrixByCol();
2054  matrix_.removeGaps();
2055  matrix_.setExtraGap(0.0);
2056  matrixByRow_ = *solver->getMatrixByRow();
2057  int numberRows = matrix_.getNumRows();
2058 
2059  rhs_ = new int[numberRows];
2060  int i;
2061  const double * rowLower = solver->getRowLower();
2062  const double * rowUpper = solver->getRowUpper();
2063  // Row copy
2064  const double * elementByRow = matrixByRow_.getElements();
2065  const int * column = matrixByRow_.getIndices();
2066  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
2067  const int * rowLength = matrixByRow_.getVectorLengths();
2068  for (i=0;i<numberRows;i++) {
2069    rhs_[i]=0;
2070    double value = rowLower[i];
2071    if (value==rowUpper[i]) {
2072      if (floor(value)==value&&value>=1.0&&value<10.0) {
2073        // check elements
2074        bool good=true;
2075        for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
2076          int iColumn = column[j];
2077          if (!solver->isBinary(iColumn))
2078            good=false;
2079          double elValue = elementByRow[j];
2080          if (floor(elValue)!=elValue||value<1.0)
2081            good=false;
2082        }
2083        if (good)
2084          rhs_[i]=(int) value;
2085      }
2086    }
2087  }
2088}
2089
2090// Copy constructor
2091CbcFollowOn::CbcFollowOn ( const CbcFollowOn & rhs)
2092  :CbcObject(rhs),
2093   matrix_(rhs.matrix_),
2094   matrixByRow_(rhs.matrixByRow_)
2095{
2096  int numberRows = matrix_.getNumRows();
2097  rhs_= CoinCopyOfArray(rhs.rhs_,numberRows);
2098}
2099
2100// Clone
2101CbcObject *
2102CbcFollowOn::clone() const
2103{
2104  return new CbcFollowOn(*this);
2105}
2106
2107// Assignment operator
2108CbcFollowOn & 
2109CbcFollowOn::operator=( const CbcFollowOn& rhs)
2110{
2111  if (this!=&rhs) {
2112    CbcObject::operator=(rhs);
2113    delete [] rhs_;
2114    matrix_ = rhs.matrix_;
2115    matrixByRow_ = rhs.matrixByRow_;
2116    int numberRows = matrix_.getNumRows();
2117    rhs_= CoinCopyOfArray(rhs.rhs_,numberRows);
2118  }
2119  return *this;
2120}
2121
2122// Destructor
2123CbcFollowOn::~CbcFollowOn ()
2124{
2125  delete [] rhs_;
2126}
2127// As some computation is needed in more than one place - returns row
2128int 
2129CbcFollowOn::gutsOfFollowOn(int & otherRow, int & preferredWay) const
2130{
2131  int whichRow=-1;
2132  otherRow=-1;
2133  int numberRows = matrix_.getNumRows();
2134 
2135  int i;
2136  // For sorting
2137  int * sort = new int [numberRows];
2138  int * isort = new int [numberRows];
2139  // Column copy
2140  //const double * element = matrix_.getElements();
2141  const int * row = matrix_.getIndices();
2142  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
2143  const int * columnLength = matrix_.getVectorLengths();
2144  // Row copy
2145  const double * elementByRow = matrixByRow_.getElements();
2146  const int * column = matrixByRow_.getIndices();
2147  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
2148  const int * rowLength = matrixByRow_.getVectorLengths();
2149  OsiSolverInterface * solver = model_->solver();
2150  const double * columnLower = solver->getColLower();
2151  const double * columnUpper = solver->getColUpper();
2152  const double * solution = solver->getColSolution();
2153  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
2154  int nSort=0;
2155  for (i=0;i<numberRows;i++) {
2156    if (rhs_[i]) {
2157      // check elements
2158      double smallest=1.0e10;
2159      double largest=0.0;
2160      int rhsValue=rhs_[i];
2161      int number1=0;
2162      int numberUnsatisfied=0;
2163      for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
2164        int iColumn = column[j];
2165        double value = elementByRow[j];
2166        double solValue = solution[iColumn];
2167        if (columnLower[iColumn]!=columnUpper[iColumn]) {
2168          smallest = CoinMin(smallest,value);
2169          largest = CoinMax(largest,value);
2170          if (value==1.0)
2171            number1++;
2172          if (solValue<1.0-integerTolerance&&solValue>integerTolerance)
2173            numberUnsatisfied++;
2174        } else {
2175          rhsValue -= (int)(value*floor(solValue+0.5));
2176        }
2177      }
2178      if (numberUnsatisfied>1) {
2179        if (smallest<largest) {
2180          // probably no good but check a few things
2181          assert (largest<=rhsValue);
2182          if (number1==1&&largest==rhsValue)
2183            printf("could fix\n");
2184        } else if (largest==rhsValue) {
2185          sort[nSort]=i;
2186          isort[nSort++]=-numberUnsatisfied;
2187        }
2188      }
2189    }
2190  }
2191  if (nSort>1) {
2192    CoinSort_2(isort,isort+nSort,sort);
2193    CoinZeroN(isort,numberRows);
2194    double * other = new double[numberRows];
2195    CoinZeroN(other,numberRows);
2196    int * which = new int[numberRows];
2197    //#define COUNT
2198#ifndef COUNT
2199    bool beforeSolution = model_->getSolutionCount()==0;
2200#endif
2201    for (int k=0;k<nSort-1;k++) {
2202      i=sort[k];
2203      int numberUnsatisfied = 0;
2204      int n=0;
2205      int j;
2206      for (j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
2207        int iColumn = column[j];
2208        if (columnLower[iColumn]!=columnUpper[iColumn]) {
2209          double solValue = solution[iColumn]-columnLower[iColumn];
2210          if (solValue<1.0-integerTolerance&&solValue>integerTolerance) {
2211            numberUnsatisfied++;
2212            for (int jj=columnStart[iColumn];jj<columnStart[iColumn]+columnLength[iColumn];jj++) {
2213              int iRow = row[jj];
2214              if (rhs_[iRow]) {
2215                other[iRow]+=solValue;
2216                if (isort[iRow]) {
2217                  isort[iRow]++;
2218                } else {
2219                  isort[iRow]=1;
2220                  which[n++]=iRow;
2221                }
2222              }
2223            }
2224          }
2225        }
2226      }
2227      double total=0.0;
2228      // Take out row
2229      double sumThis=other[i];
2230      other[i]=0.0;
2231      assert (numberUnsatisfied==isort[i]);
2232      // find one nearest half if solution, one if before solution
2233      int iBest=-1;
2234      double dtarget=0.5*total;
2235#ifdef COUNT
2236      int target = (numberUnsatisfied+1)>>1;
2237      int best=numberUnsatisfied;
2238#else
2239      double best;
2240      if (beforeSolution)
2241        best=dtarget;
2242      else
2243        best=1.0e30;
2244#endif
2245      for (j=0;j<n;j++) {
2246        int iRow = which[j];
2247        double dvalue=other[iRow];
2248        other[iRow]=0.0;
2249#ifdef COUNT
2250        int value = isort[iRow];
2251#endif
2252        isort[iRow]=0;
2253        if (fabs(dvalue)<1.0e-8||fabs(sumThis-dvalue)<1.0e-8)
2254          continue;
2255        if (dvalue<integerTolerance||dvalue>1.0-integerTolerance)
2256          continue;
2257#ifdef COUNT
2258        if (abs(value-target)<best&&value!=numberUnsatisfied) {
2259          best=abs(value-target);
2260          iBest=iRow;
2261          if (dvalue<dtarget)
2262            preferredWay=1;
2263          else
2264            preferredWay=-1;
2265        }
2266#else
2267        if (beforeSolution) {
2268          if (fabs(dvalue-dtarget)>best) {
2269            best = fabs(dvalue-dtarget);
2270            iBest=iRow;
2271            if (dvalue<dtarget)
2272              preferredWay=1;
2273            else
2274              preferredWay=-1;
2275          }
2276        } else {
2277          if (fabs(dvalue-dtarget)<best) {
2278            best = fabs(dvalue-dtarget);
2279            iBest=iRow;
2280            if (dvalue<dtarget)
2281              preferredWay=1;
2282            else
2283              preferredWay=-1;
2284          }
2285        }
2286#endif
2287      }
2288      if (iBest>=0) {
2289        whichRow=i;
2290        otherRow=iBest;
2291        break;
2292      }
2293    }
2294    delete [] which;
2295    delete [] other;
2296  }
2297  delete [] sort;
2298  delete [] isort;
2299  return whichRow;
2300}
2301
2302// Infeasibility - large is 0.5
2303double 
2304CbcFollowOn::infeasibility(int & preferredWay) const
2305{
2306  int otherRow=0;
2307  int whichRow = gutsOfFollowOn(otherRow,preferredWay);
2308  if (whichRow<0)
2309    return 0.0;
2310  else
2311  return 2.0* model_->getDblParam(CbcModel::CbcIntegerTolerance);
2312}
2313
2314// This looks at solution and sets bounds to contain solution
2315void 
2316CbcFollowOn::feasibleRegion()
2317{
2318}
2319
2320
2321// Creates a branching object
2322CbcBranchingObject * 
2323CbcFollowOn::createBranch(int way) const
2324{
2325  int otherRow=0;
2326  int preferredWay;
2327  int whichRow = gutsOfFollowOn(otherRow,preferredWay);
2328  assert(way==preferredWay);
2329  assert (whichRow>=0);
2330  int numberColumns = matrix_.getNumCols();
2331 
2332  // Column copy
2333  //const double * element = matrix_.getElements();
2334  const int * row = matrix_.getIndices();
2335  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
2336  const int * columnLength = matrix_.getVectorLengths();
2337  // Row copy
2338  //const double * elementByRow = matrixByRow_.getElements();
2339  const int * column = matrixByRow_.getIndices();
2340  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
2341  const int * rowLength = matrixByRow_.getVectorLengths();
2342  OsiSolverInterface * solver = model_->solver();
2343  const double * columnLower = solver->getColLower();
2344  const double * columnUpper = solver->getColUpper();
2345  //const double * solution = solver->getColSolution();
2346  int nUp=0;
2347  int nDown=0;
2348  int * upList = new int[numberColumns];
2349  int * downList = new int[numberColumns];
2350  int j;
2351  for (j=rowStart[whichRow];j<rowStart[whichRow]+rowLength[whichRow];j++) {
2352    int iColumn = column[j];
2353    if (columnLower[iColumn]!=columnUpper[iColumn]) {
2354      bool up=true;
2355      for (int jj=columnStart[iColumn];jj<columnStart[iColumn]+columnLength[iColumn];jj++) {
2356        int iRow = row[jj];
2357        if (iRow==otherRow) {
2358          up=false;
2359          break;
2360        }
2361      }
2362      if (up)
2363        upList[nUp++]=iColumn;
2364      else
2365        downList[nDown++]=iColumn;
2366    }
2367  }
2368  //printf("way %d\n",way);
2369  // create object
2370  //printf("would fix %d down and %d up\n",nDown,nUp);
2371  CbcBranchingObject * branch
2372     = new CbcFixingBranchingObject(model_,way,
2373                                         nDown,downList,nUp,upList);
2374  delete [] upList;
2375  delete [] downList;
2376  return branch;
2377}
2378// Default Constructor
2379CbcFixingBranchingObject::CbcFixingBranchingObject()
2380  :CbcBranchingObject()
2381{
2382  numberDown_=0;
2383  numberUp_=0;
2384  downList_=NULL;
2385  upList_=NULL;
2386}
2387
2388// Useful constructor
2389CbcFixingBranchingObject::CbcFixingBranchingObject (CbcModel * model,
2390                                                    int way ,
2391                                                    int numberOnDownSide, const int * down,
2392                                                    int numberOnUpSide, const int * up)
2393  :CbcBranchingObject(model,0,way,0.5)
2394{
2395  numberDown_=numberOnDownSide;
2396  numberUp_=numberOnUpSide;
2397  downList_ = CoinCopyOfArray(down,numberDown_);
2398  upList_ = CoinCopyOfArray(up,numberUp_);
2399}
2400
2401// Copy constructor
2402CbcFixingBranchingObject::CbcFixingBranchingObject ( const CbcFixingBranchingObject & rhs) :CbcBranchingObject(rhs)
2403{
2404  numberDown_=rhs.numberDown_;
2405  numberUp_=rhs.numberUp_;
2406  downList_ = CoinCopyOfArray(rhs.downList_,numberDown_);
2407  upList_ = CoinCopyOfArray(rhs.upList_,numberUp_);
2408}
2409
2410// Assignment operator
2411CbcFixingBranchingObject & 
2412CbcFixingBranchingObject::operator=( const CbcFixingBranchingObject& rhs)
2413{
2414  if (this != &rhs) {
2415    CbcBranchingObject::operator=(rhs);
2416    delete [] downList_;
2417    delete [] upList_;
2418    numberDown_=rhs.numberDown_;
2419    numberUp_=rhs.numberUp_;
2420    downList_ = CoinCopyOfArray(rhs.downList_,numberDown_);
2421    upList_ = CoinCopyOfArray(rhs.upList_,numberUp_);
2422  }
2423  return *this;
2424}
2425CbcBranchingObject * 
2426CbcFixingBranchingObject::clone() const
2427{ 
2428  return (new CbcFixingBranchingObject(*this));
2429}
2430
2431
2432// Destructor
2433CbcFixingBranchingObject::~CbcFixingBranchingObject ()
2434{
2435  delete [] downList_;
2436  delete [] upList_;
2437}
2438double
2439CbcFixingBranchingObject::branch(bool normalBranch)
2440{
2441  if (model_->messageHandler()->logLevel()>2&&normalBranch)
2442    print(normalBranch);
2443  numberBranchesLeft_--;
2444  OsiSolverInterface * solver = model_->solver();
2445  const double * columnLower = solver->getColLower();
2446  int i;
2447  // *** for way - up means fix all those in up section
2448  if (way_<0) {
2449#ifdef FULL_PRINT
2450    printf("Down Fix ");
2451#endif
2452    //printf("Down Fix %d\n",numberDown_);
2453    for (i=0;i<numberDown_;i++) {
2454      int iColumn = downList_[i];
2455      model_->solver()->setColUpper(iColumn,columnLower[iColumn]);
2456#ifdef FULL_PRINT
2457      printf("Setting bound on %d to lower bound\n",iColumn);
2458#endif
2459    }
2460    way_=1;       // Swap direction
2461  } else {
2462#ifdef FULL_PRINT
2463    printf("Up Fix ");
2464#endif
2465    //printf("Up Fix %d\n",numberUp_);
2466    for (i=0;i<numberUp_;i++) {
2467      int iColumn = upList_[i];
2468      model_->solver()->setColUpper(iColumn,columnLower[iColumn]);
2469#ifdef FULL_PRINT
2470      printf("Setting bound on %d to lower bound\n",iColumn);
2471#endif
2472    }
2473    way_=-1;      // Swap direction
2474  }
2475#ifdef FULL_PRINT
2476  printf("\n");
2477#endif
2478  return 0.0;
2479}
2480void
2481CbcFixingBranchingObject::print(bool normalBranch)
2482{
2483  int i;
2484  // *** for way - up means fix all those in up section
2485  if (way_<0) {
2486    printf("Down Fix ");
2487    for (i=0;i<numberDown_;i++) {
2488      int iColumn = downList_[i];
2489      printf("%d ",iColumn);
2490    }
2491  } else {
2492    printf("Up Fix ");
2493    for (i=0;i<numberUp_;i++) {
2494      int iColumn = upList_[i];
2495      printf("%d ",iColumn);
2496    }
2497  }
2498  printf("\n");
2499}
Note: See TracBrowser for help on using the repository browser.