source: trunk/CbcBranchActual.cpp @ 129

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

make createBranch non const

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 68.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) 
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) 
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) 
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 (int var %d): [%g,%g] => [%g,%g]\n",
921           iColumn,variable_,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 (int var %d): [%g,%g] => [%g,%g]\n",
927           iColumn,variable_,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) 
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  bestCriterion_ = 0.0;
1735  bestChangeUp_ = 0.0;
1736  bestNumberUp_ = 0;
1737  bestChangeDown_ = 0.0;
1738  bestNumberDown_ = 0;
1739  bestObject_ = NULL;
1740}
1741
1742// Copy constructor
1743CbcBranchDefaultDecision::CbcBranchDefaultDecision (
1744                                    const CbcBranchDefaultDecision & rhs)
1745  :CbcBranchDecision()
1746{
1747  bestCriterion_ = rhs.bestCriterion_;
1748  bestChangeUp_ = rhs.bestChangeUp_;
1749  bestNumberUp_ = rhs.bestNumberUp_;
1750  bestChangeDown_ = rhs.bestChangeDown_;
1751  bestNumberDown_ = rhs.bestNumberDown_;
1752  bestObject_ = rhs.bestObject_;
1753}
1754
1755CbcBranchDefaultDecision::~CbcBranchDefaultDecision()
1756{
1757}
1758
1759// Clone
1760CbcBranchDecision * 
1761CbcBranchDefaultDecision::clone() const
1762{
1763  return new CbcBranchDefaultDecision(*this);
1764}
1765
1766// Initialize i.e. before start of choosing at a node
1767void 
1768CbcBranchDefaultDecision::initialize(CbcModel * model)
1769{
1770  bestCriterion_ = 0.0;
1771  bestChangeUp_ = 0.0;
1772  bestNumberUp_ = 0;
1773  bestChangeDown_ = 0.0;
1774  bestNumberDown_ = 0;
1775  bestObject_ = NULL;
1776}
1777
1778
1779/*
1780  Simple default decision algorithm. Compare based on infeasibility (numInfUp,
1781  numInfDn) until a solution is found by search, then switch to change in
1782  objective (changeUp, changeDn). Note that bestSoFar is remembered in
1783  bestObject_, so the parameter bestSoFar is unused.
1784*/
1785
1786int
1787CbcBranchDefaultDecision::betterBranch(CbcBranchingObject * thisOne,
1788                            CbcBranchingObject * bestSoFar,
1789                            double changeUp, int numInfUp,
1790                            double changeDn, int numInfDn)
1791{
1792  bool beforeSolution = thisOne->model()->getSolutionCount()==
1793    thisOne->model()->getNumberHeuristicSolutions();;
1794  int betterWay=0;
1795  if (beforeSolution) {
1796    if (!bestObject_) {
1797      bestNumberUp_=INT_MAX;
1798      bestNumberDown_=INT_MAX;
1799    }
1800    // before solution - choose smallest number
1801    // could add in depth as well
1802    int bestNumber = CoinMin(bestNumberUp_,bestNumberDown_);
1803    if (numInfUp<numInfDn) {
1804      if (numInfUp<bestNumber) {
1805        betterWay = 1;
1806      } else if (numInfUp==bestNumber) {
1807        if (changeUp<bestCriterion_)
1808          betterWay=1;
1809      }
1810    } else if (numInfUp>numInfDn) {
1811      if (numInfDn<bestNumber) {
1812        betterWay = -1;
1813      } else if (numInfDn==bestNumber) {
1814        if (changeDn<bestCriterion_)
1815          betterWay=-1;
1816      }
1817    } else {
1818      // up and down have same number
1819      bool better=false;
1820      if (numInfUp<bestNumber) {
1821        better=true;
1822      } else if (numInfUp==bestNumber) {
1823        if (min(changeUp,changeDn)<bestCriterion_)
1824          better=true;;
1825      }
1826      if (better) {
1827        // see which way
1828        if (changeUp<=changeDn)
1829          betterWay=1;
1830        else
1831          betterWay=-1;
1832      }
1833    }
1834  } else {
1835    if (!bestObject_) {
1836      bestCriterion_=-1.0;
1837    }
1838    // got a solution
1839    if (changeUp<=changeDn) {
1840      if (changeUp>bestCriterion_)
1841        betterWay=1;
1842    } else {
1843      if (changeDn>bestCriterion_)
1844        betterWay=-1;
1845    }
1846  }
1847  if (betterWay) {
1848    bestCriterion_ = CoinMin(changeUp,changeDn);
1849    bestChangeUp_ = changeUp;
1850    bestNumberUp_ = numInfUp;
1851    bestChangeDown_ = changeDn;
1852    bestNumberDown_ = numInfDn;
1853    bestObject_=thisOne;
1854  }
1855  return betterWay;
1856}
1857
1858/* Compare N branching objects. Return index of best
1859   and sets way of branching in chosen object.
1860   
1861   This routine is used only after strong branching.
1862*/
1863
1864int
1865CbcBranchDefaultDecision::bestBranch (CbcBranchingObject ** objects, int numberObjects,
1866                                   int numberUnsatisfied,
1867                                   double * changeUp, int * numberInfeasibilitiesUp,
1868                                   double * changeDown, int * numberInfeasibilitiesDown,
1869                                   double objectiveValue) 
1870{
1871
1872  int bestWay=0;
1873  int whichObject = -1;
1874  if (numberObjects) {
1875    CbcModel * model = objects[0]->model();
1876    // at continuous
1877    //double continuousObjective = model->getContinuousObjective();
1878    //int continuousInfeasibilities = model->getContinuousInfeasibilities();
1879   
1880    // average cost to get rid of infeasibility
1881    //double averageCostPerInfeasibility =
1882    //(objectiveValue-continuousObjective)/
1883    //(double) (abs(continuousInfeasibilities-numberUnsatisfied)+1);
1884    /* beforeSolution is :
1885       0 - before any solution
1886       n - n heuristic solutions but no branched one
1887       -1 - branched solution found
1888    */
1889    int numberSolutions = model->getSolutionCount();
1890    double cutoff = model->getCutoff();
1891    int method=0;
1892    int i;
1893    if (numberSolutions) {
1894      int numberHeuristic = model->getNumberHeuristicSolutions();
1895      if (numberHeuristic<numberSolutions) {
1896        method = 1;
1897      } else {
1898        method = 2;
1899        // look further
1900        for ( i = 0 ; i < numberObjects ; i++) {
1901          int numberNext = numberInfeasibilitiesUp[i];
1902         
1903          if (numberNext<numberUnsatisfied) {
1904            int numberUp = numberUnsatisfied - numberInfeasibilitiesUp[i];
1905            double perUnsatisfied = changeUp[i]/(double) numberUp;
1906            double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
1907            if (estimatedObjective<cutoff) 
1908              method=3;
1909          }
1910          numberNext = numberInfeasibilitiesDown[i];
1911          if (numberNext<numberUnsatisfied) {
1912            int numberDown = numberUnsatisfied - numberInfeasibilitiesDown[i];
1913            double perUnsatisfied = changeDown[i]/(double) numberDown;
1914            double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
1915            if (estimatedObjective<cutoff) 
1916              method=3;
1917          }
1918        }
1919      }
1920      method=2;
1921    } else {
1922      method = 0;
1923    }
1924    // Uncomment next to force method 4
1925    //method=4;
1926    /* Methods :
1927       0 - fewest infeasibilities
1928       1 - largest min change in objective
1929       2 - as 1 but use sum of changes if min close
1930       3 - predicted best solution
1931       4 - take cheapest up branch if infeasibilities same
1932    */
1933    int bestNumber=INT_MAX;
1934    double bestCriterion=-1.0e50;
1935    double alternativeCriterion = -1.0;
1936    double bestEstimate = 1.0e100;
1937    switch (method) {
1938    case 0:
1939      // could add in depth as well
1940      for ( i = 0 ; i < numberObjects ; i++) {
1941        int thisNumber = min(numberInfeasibilitiesUp[i],numberInfeasibilitiesDown[i]);
1942        if (thisNumber<=bestNumber) {
1943          int betterWay=0;
1944          if (numberInfeasibilitiesUp[i]<numberInfeasibilitiesDown[i]) {
1945            if (numberInfeasibilitiesUp[i]<bestNumber) {
1946              betterWay = 1;
1947            } else {
1948              if (changeUp[i]<bestCriterion)
1949                betterWay=1;
1950            }
1951          } else if (numberInfeasibilitiesUp[i]>numberInfeasibilitiesDown[i]) {
1952            if (numberInfeasibilitiesDown[i]<bestNumber) {
1953              betterWay = -1;
1954            } else {
1955              if (changeDown[i]<bestCriterion)
1956                betterWay=-1;
1957            }
1958          } else {
1959            // up and down have same number
1960            bool better=false;
1961            if (numberInfeasibilitiesUp[i]<bestNumber) {
1962              better=true;
1963            } else if (numberInfeasibilitiesUp[i]==bestNumber) {
1964              if (min(changeUp[i],changeDown[i])<bestCriterion)
1965                better=true;;
1966            }
1967            if (better) {
1968              // see which way
1969              if (changeUp[i]<=changeDown[i])
1970                betterWay=1;
1971              else
1972                betterWay=-1;
1973            }
1974          }
1975          if (betterWay) {
1976            bestCriterion = min(changeUp[i],changeDown[i]);
1977            bestNumber = thisNumber;
1978            whichObject = i;
1979            bestWay = betterWay;
1980          }
1981        }
1982      }
1983      break;
1984    case 1:
1985      for ( i = 0 ; i < numberObjects ; i++) {
1986        int betterWay=0;
1987        if (changeUp[i]<=changeDown[i]) {
1988          if (changeUp[i]>bestCriterion)
1989            betterWay=1;
1990        } else {
1991          if (changeDown[i]>bestCriterion)
1992            betterWay=-1;
1993        }
1994        if (betterWay) {
1995          bestCriterion = min(changeUp[i],changeDown[i]);
1996          whichObject = i;
1997          bestWay = betterWay;
1998        }
1999      }
2000      break;
2001    case 2:
2002      for ( i = 0 ; i < numberObjects ; i++) {
2003        double change = min(changeUp[i],changeDown[i]);
2004        double sum = changeUp[i] + changeDown[i];
2005        bool take=false;
2006        if (change>1.1*bestCriterion) 
2007          take=true;
2008        else if (change>0.9*bestCriterion&&sum+change>bestCriterion+alternativeCriterion) 
2009          take=true;
2010        if (take) {
2011          if (changeUp[i]<=changeDown[i]) {
2012            if (changeUp[i]>bestCriterion)
2013              bestWay=1;
2014          } else {
2015            if (changeDown[i]>bestCriterion)
2016              bestWay=-1;
2017          }
2018          bestCriterion = change;
2019          alternativeCriterion = sum;
2020          whichObject = i;
2021        }
2022      }
2023      break;
2024    case 3:
2025      for ( i = 0 ; i < numberObjects ; i++) {
2026        int numberNext = numberInfeasibilitiesUp[i];
2027       
2028        if (numberNext<numberUnsatisfied) {
2029          int numberUp = numberUnsatisfied - numberInfeasibilitiesUp[i];
2030          double perUnsatisfied = changeUp[i]/(double) numberUp;
2031          double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
2032          if (estimatedObjective<bestEstimate) {
2033            bestEstimate = estimatedObjective;
2034            bestWay=1;
2035            whichObject=i;
2036          }
2037        }
2038        numberNext = numberInfeasibilitiesDown[i];
2039        if (numberNext<numberUnsatisfied) {
2040          int numberDown = numberUnsatisfied - numberInfeasibilitiesDown[i];
2041          double perUnsatisfied = changeDown[i]/(double) numberDown;
2042          double estimatedObjective = objectiveValue + numberUnsatisfied * perUnsatisfied;
2043          if (estimatedObjective<bestEstimate) {
2044            bestEstimate = estimatedObjective;
2045            bestWay=-1;
2046            whichObject=i;
2047          }
2048        }
2049      }
2050      break;
2051    case 4:
2052      // if number infeas same then cheapest up
2053      // first get best number or when going down
2054      // now choose smallest change up amongst equal number infeas
2055      for ( i = 0 ; i < numberObjects ; i++) {
2056        int thisNumber = min(numberInfeasibilitiesUp[i],numberInfeasibilitiesDown[i]);
2057        if (thisNumber<=bestNumber) {
2058          int betterWay=0;
2059          if (numberInfeasibilitiesUp[i]<numberInfeasibilitiesDown[i]) {
2060            if (numberInfeasibilitiesUp[i]<bestNumber) {
2061              betterWay = 1;
2062            } else {
2063              if (changeUp[i]<bestCriterion)
2064                betterWay=1;
2065            }
2066          } else if (numberInfeasibilitiesUp[i]>numberInfeasibilitiesDown[i]) {
2067            if (numberInfeasibilitiesDown[i]<bestNumber) {
2068              betterWay = -1;
2069            } else {
2070              if (changeDown[i]<bestCriterion)
2071                betterWay=-1;
2072            }
2073          } else {
2074            // up and down have same number
2075            bool better=false;
2076            if (numberInfeasibilitiesUp[i]<bestNumber) {
2077              better=true;
2078            } else if (numberInfeasibilitiesUp[i]==bestNumber) {
2079              if (min(changeUp[i],changeDown[i])<bestCriterion)
2080                better=true;;
2081            }
2082            if (better) {
2083              // see which way
2084              if (changeUp[i]<=changeDown[i])
2085                betterWay=1;
2086              else
2087                betterWay=-1;
2088            }
2089          }
2090          if (betterWay) {
2091            bestCriterion = min(changeUp[i],changeDown[i]);
2092            bestNumber = thisNumber;
2093            whichObject = i;
2094            bestWay = betterWay;
2095          }
2096        }
2097      }
2098      bestCriterion=1.0e50;
2099      for ( i = 0 ; i < numberObjects ; i++) {
2100        int thisNumber = numberInfeasibilitiesUp[i];
2101        if (thisNumber==bestNumber&&changeUp) {
2102          if (changeUp[i]<bestCriterion) {
2103            bestCriterion = changeUp[i];
2104            whichObject = i;
2105            bestWay = 1;
2106          }
2107        }
2108      }
2109      break;
2110    }
2111    // set way in best
2112    if (whichObject>=0)
2113      objects[whichObject]->way(bestWay);
2114  }
2115  return whichObject;
2116}
2117
2118// Default Constructor
2119CbcFollowOn::CbcFollowOn ()
2120  : CbcObject(),
2121    rhs_(NULL)
2122{
2123}
2124
2125// Useful constructor
2126CbcFollowOn::CbcFollowOn (CbcModel * model)
2127  : CbcObject(model)
2128{
2129  assert (model);
2130  OsiSolverInterface * solver = model_->solver();
2131  matrix_ = *solver->getMatrixByCol();
2132  matrix_.removeGaps();
2133  matrix_.setExtraGap(0.0);
2134  matrixByRow_ = *solver->getMatrixByRow();
2135  int numberRows = matrix_.getNumRows();
2136 
2137  rhs_ = new int[numberRows];
2138  int i;
2139  const double * rowLower = solver->getRowLower();
2140  const double * rowUpper = solver->getRowUpper();
2141  // Row copy
2142  const double * elementByRow = matrixByRow_.getElements();
2143  const int * column = matrixByRow_.getIndices();
2144  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
2145  const int * rowLength = matrixByRow_.getVectorLengths();
2146  for (i=0;i<numberRows;i++) {
2147    rhs_[i]=0;
2148    double value = rowLower[i];
2149    if (value==rowUpper[i]) {
2150      if (floor(value)==value&&value>=1.0&&value<10.0) {
2151        // check elements
2152        bool good=true;
2153        for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
2154          int iColumn = column[j];
2155          if (!solver->isBinary(iColumn))
2156            good=false;
2157          double elValue = elementByRow[j];
2158          if (floor(elValue)!=elValue||value<1.0)
2159            good=false;
2160        }
2161        if (good)
2162          rhs_[i]=(int) value;
2163      }
2164    }
2165  }
2166}
2167
2168// Copy constructor
2169CbcFollowOn::CbcFollowOn ( const CbcFollowOn & rhs)
2170  :CbcObject(rhs),
2171   matrix_(rhs.matrix_),
2172   matrixByRow_(rhs.matrixByRow_)
2173{
2174  int numberRows = matrix_.getNumRows();
2175  rhs_= CoinCopyOfArray(rhs.rhs_,numberRows);
2176}
2177
2178// Clone
2179CbcObject *
2180CbcFollowOn::clone() const
2181{
2182  return new CbcFollowOn(*this);
2183}
2184
2185// Assignment operator
2186CbcFollowOn & 
2187CbcFollowOn::operator=( const CbcFollowOn& rhs)
2188{
2189  if (this!=&rhs) {
2190    CbcObject::operator=(rhs);
2191    delete [] rhs_;
2192    matrix_ = rhs.matrix_;
2193    matrixByRow_ = rhs.matrixByRow_;
2194    int numberRows = matrix_.getNumRows();
2195    rhs_= CoinCopyOfArray(rhs.rhs_,numberRows);
2196  }
2197  return *this;
2198}
2199
2200// Destructor
2201CbcFollowOn::~CbcFollowOn ()
2202{
2203  delete [] rhs_;
2204}
2205// As some computation is needed in more than one place - returns row
2206int 
2207CbcFollowOn::gutsOfFollowOn(int & otherRow, int & preferredWay) const
2208{
2209  int whichRow=-1;
2210  otherRow=-1;
2211  int numberRows = matrix_.getNumRows();
2212 
2213  int i;
2214  // For sorting
2215  int * sort = new int [numberRows];
2216  int * isort = new int [numberRows];
2217  // Column copy
2218  //const double * element = matrix_.getElements();
2219  const int * row = matrix_.getIndices();
2220  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
2221  const int * columnLength = matrix_.getVectorLengths();
2222  // Row copy
2223  const double * elementByRow = matrixByRow_.getElements();
2224  const int * column = matrixByRow_.getIndices();
2225  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
2226  const int * rowLength = matrixByRow_.getVectorLengths();
2227  OsiSolverInterface * solver = model_->solver();
2228  const double * columnLower = solver->getColLower();
2229  const double * columnUpper = solver->getColUpper();
2230  const double * solution = solver->getColSolution();
2231  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
2232  int nSort=0;
2233  for (i=0;i<numberRows;i++) {
2234    if (rhs_[i]) {
2235      // check elements
2236      double smallest=1.0e10;
2237      double largest=0.0;
2238      int rhsValue=rhs_[i];
2239      int number1=0;
2240      int numberUnsatisfied=0;
2241      for (int j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
2242        int iColumn = column[j];
2243        double value = elementByRow[j];
2244        double solValue = solution[iColumn];
2245        if (columnLower[iColumn]!=columnUpper[iColumn]) {
2246          smallest = CoinMin(smallest,value);
2247          largest = CoinMax(largest,value);
2248          if (value==1.0)
2249            number1++;
2250          if (solValue<1.0-integerTolerance&&solValue>integerTolerance)
2251            numberUnsatisfied++;
2252        } else {
2253          rhsValue -= (int)(value*floor(solValue+0.5));
2254        }
2255      }
2256      if (numberUnsatisfied>1) {
2257        if (smallest<largest) {
2258          // probably no good but check a few things
2259          assert (largest<=rhsValue);
2260          if (number1==1&&largest==rhsValue)
2261            printf("could fix\n");
2262        } else if (largest==rhsValue) {
2263          sort[nSort]=i;
2264          isort[nSort++]=-numberUnsatisfied;
2265        }
2266      }
2267    }
2268  }
2269  if (nSort>1) {
2270    CoinSort_2(isort,isort+nSort,sort);
2271    CoinZeroN(isort,numberRows);
2272    double * other = new double[numberRows];
2273    CoinZeroN(other,numberRows);
2274    int * which = new int[numberRows];
2275    //#define COUNT
2276#ifndef COUNT
2277    bool beforeSolution = model_->getSolutionCount()==0;
2278#endif
2279    for (int k=0;k<nSort-1;k++) {
2280      i=sort[k];
2281      int numberUnsatisfied = 0;
2282      int n=0;
2283      int j;
2284      for (j=rowStart[i];j<rowStart[i]+rowLength[i];j++) {
2285        int iColumn = column[j];
2286        if (columnLower[iColumn]!=columnUpper[iColumn]) {
2287          double solValue = solution[iColumn]-columnLower[iColumn];
2288          if (solValue<1.0-integerTolerance&&solValue>integerTolerance) {
2289            numberUnsatisfied++;
2290            for (int jj=columnStart[iColumn];jj<columnStart[iColumn]+columnLength[iColumn];jj++) {
2291              int iRow = row[jj];
2292              if (rhs_[iRow]) {
2293                other[iRow]+=solValue;
2294                if (isort[iRow]) {
2295                  isort[iRow]++;
2296                } else {
2297                  isort[iRow]=1;
2298                  which[n++]=iRow;
2299                }
2300              }
2301            }
2302          }
2303        }
2304      }
2305      double total=0.0;
2306      // Take out row
2307      double sumThis=other[i];
2308      other[i]=0.0;
2309      assert (numberUnsatisfied==isort[i]);
2310      // find one nearest half if solution, one if before solution
2311      int iBest=-1;
2312      double dtarget=0.5*total;
2313#ifdef COUNT
2314      int target = (numberUnsatisfied+1)>>1;
2315      int best=numberUnsatisfied;
2316#else
2317      double best;
2318      if (beforeSolution)
2319        best=dtarget;
2320      else
2321        best=1.0e30;
2322#endif
2323      for (j=0;j<n;j++) {
2324        int iRow = which[j];
2325        double dvalue=other[iRow];
2326        other[iRow]=0.0;
2327#ifdef COUNT
2328        int value = isort[iRow];
2329#endif
2330        isort[iRow]=0;
2331        if (fabs(dvalue)<1.0e-8||fabs(sumThis-dvalue)<1.0e-8)
2332          continue;
2333        if (dvalue<integerTolerance||dvalue>1.0-integerTolerance)
2334          continue;
2335#ifdef COUNT
2336        if (abs(value-target)<best&&value!=numberUnsatisfied) {
2337          best=abs(value-target);
2338          iBest=iRow;
2339          if (dvalue<dtarget)
2340            preferredWay=1;
2341          else
2342            preferredWay=-1;
2343        }
2344#else
2345        if (beforeSolution) {
2346          if (fabs(dvalue-dtarget)>best) {
2347            best = fabs(dvalue-dtarget);
2348            iBest=iRow;
2349            if (dvalue<dtarget)
2350              preferredWay=1;
2351            else
2352              preferredWay=-1;
2353          }
2354        } else {
2355          if (fabs(dvalue-dtarget)<best) {
2356            best = fabs(dvalue-dtarget);
2357            iBest=iRow;
2358            if (dvalue<dtarget)
2359              preferredWay=1;
2360            else
2361              preferredWay=-1;
2362          }
2363        }
2364#endif
2365      }
2366      if (iBest>=0) {
2367        whichRow=i;
2368        otherRow=iBest;
2369        break;
2370      }
2371    }
2372    delete [] which;
2373    delete [] other;
2374  }
2375  delete [] sort;
2376  delete [] isort;
2377  return whichRow;
2378}
2379
2380// Infeasibility - large is 0.5
2381double 
2382CbcFollowOn::infeasibility(int & preferredWay) const
2383{
2384  int otherRow=0;
2385  int whichRow = gutsOfFollowOn(otherRow,preferredWay);
2386  if (whichRow<0)
2387    return 0.0;
2388  else
2389  return 2.0* model_->getDblParam(CbcModel::CbcIntegerTolerance);
2390}
2391
2392// This looks at solution and sets bounds to contain solution
2393void 
2394CbcFollowOn::feasibleRegion()
2395{
2396}
2397
2398
2399// Creates a branching object
2400CbcBranchingObject * 
2401CbcFollowOn::createBranch(int way) 
2402{
2403  int otherRow=0;
2404  int preferredWay;
2405  int whichRow = gutsOfFollowOn(otherRow,preferredWay);
2406  assert(way==preferredWay);
2407  assert (whichRow>=0);
2408  int numberColumns = matrix_.getNumCols();
2409 
2410  // Column copy
2411  //const double * element = matrix_.getElements();
2412  const int * row = matrix_.getIndices();
2413  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
2414  const int * columnLength = matrix_.getVectorLengths();
2415  // Row copy
2416  //const double * elementByRow = matrixByRow_.getElements();
2417  const int * column = matrixByRow_.getIndices();
2418  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
2419  const int * rowLength = matrixByRow_.getVectorLengths();
2420  OsiSolverInterface * solver = model_->solver();
2421  const double * columnLower = solver->getColLower();
2422  const double * columnUpper = solver->getColUpper();
2423  //const double * solution = solver->getColSolution();
2424  int nUp=0;
2425  int nDown=0;
2426  int * upList = new int[numberColumns];
2427  int * downList = new int[numberColumns];
2428  int j;
2429  for (j=rowStart[whichRow];j<rowStart[whichRow]+rowLength[whichRow];j++) {
2430    int iColumn = column[j];
2431    if (columnLower[iColumn]!=columnUpper[iColumn]) {
2432      bool up=true;
2433      for (int jj=columnStart[iColumn];jj<columnStart[iColumn]+columnLength[iColumn];jj++) {
2434        int iRow = row[jj];
2435        if (iRow==otherRow) {
2436          up=false;
2437          break;
2438        }
2439      }
2440      if (up)
2441        upList[nUp++]=iColumn;
2442      else
2443        downList[nDown++]=iColumn;
2444    }
2445  }
2446  //printf("way %d\n",way);
2447  // create object
2448  //printf("would fix %d down and %d up\n",nDown,nUp);
2449  CbcBranchingObject * branch
2450     = new CbcFixingBranchingObject(model_,way,
2451                                         nDown,downList,nUp,upList);
2452  delete [] upList;
2453  delete [] downList;
2454  return branch;
2455}
2456// Default Constructor
2457CbcFixingBranchingObject::CbcFixingBranchingObject()
2458  :CbcBranchingObject()
2459{
2460  numberDown_=0;
2461  numberUp_=0;
2462  downList_=NULL;
2463  upList_=NULL;
2464}
2465
2466// Useful constructor
2467CbcFixingBranchingObject::CbcFixingBranchingObject (CbcModel * model,
2468                                                    int way ,
2469                                                    int numberOnDownSide, const int * down,
2470                                                    int numberOnUpSide, const int * up)
2471  :CbcBranchingObject(model,0,way,0.5)
2472{
2473  numberDown_=numberOnDownSide;
2474  numberUp_=numberOnUpSide;
2475  downList_ = CoinCopyOfArray(down,numberDown_);
2476  upList_ = CoinCopyOfArray(up,numberUp_);
2477}
2478
2479// Copy constructor
2480CbcFixingBranchingObject::CbcFixingBranchingObject ( const CbcFixingBranchingObject & rhs) :CbcBranchingObject(rhs)
2481{
2482  numberDown_=rhs.numberDown_;
2483  numberUp_=rhs.numberUp_;
2484  downList_ = CoinCopyOfArray(rhs.downList_,numberDown_);
2485  upList_ = CoinCopyOfArray(rhs.upList_,numberUp_);
2486}
2487
2488// Assignment operator
2489CbcFixingBranchingObject & 
2490CbcFixingBranchingObject::operator=( const CbcFixingBranchingObject& rhs)
2491{
2492  if (this != &rhs) {
2493    CbcBranchingObject::operator=(rhs);
2494    delete [] downList_;
2495    delete [] upList_;
2496    numberDown_=rhs.numberDown_;
2497    numberUp_=rhs.numberUp_;
2498    downList_ = CoinCopyOfArray(rhs.downList_,numberDown_);
2499    upList_ = CoinCopyOfArray(rhs.upList_,numberUp_);
2500  }
2501  return *this;
2502}
2503CbcBranchingObject * 
2504CbcFixingBranchingObject::clone() const
2505{ 
2506  return (new CbcFixingBranchingObject(*this));
2507}
2508
2509
2510// Destructor
2511CbcFixingBranchingObject::~CbcFixingBranchingObject ()
2512{
2513  delete [] downList_;
2514  delete [] upList_;
2515}
2516double
2517CbcFixingBranchingObject::branch(bool normalBranch)
2518{
2519  if (model_->messageHandler()->logLevel()>2&&normalBranch)
2520    print(normalBranch);
2521  numberBranchesLeft_--;
2522  OsiSolverInterface * solver = model_->solver();
2523  const double * columnLower = solver->getColLower();
2524  int i;
2525  // *** for way - up means fix all those in up section
2526  if (way_<0) {
2527#ifdef FULL_PRINT
2528    printf("Down Fix ");
2529#endif
2530    //printf("Down Fix %d\n",numberDown_);
2531    for (i=0;i<numberDown_;i++) {
2532      int iColumn = downList_[i];
2533      model_->solver()->setColUpper(iColumn,columnLower[iColumn]);
2534#ifdef FULL_PRINT
2535      printf("Setting bound on %d to lower bound\n",iColumn);
2536#endif
2537    }
2538    way_=1;       // Swap direction
2539  } else {
2540#ifdef FULL_PRINT
2541    printf("Up Fix ");
2542#endif
2543    //printf("Up Fix %d\n",numberUp_);
2544    for (i=0;i<numberUp_;i++) {
2545      int iColumn = upList_[i];
2546      model_->solver()->setColUpper(iColumn,columnLower[iColumn]);
2547#ifdef FULL_PRINT
2548      printf("Setting bound on %d to lower bound\n",iColumn);
2549#endif
2550    }
2551    way_=-1;      // Swap direction
2552  }
2553#ifdef FULL_PRINT
2554  printf("\n");
2555#endif
2556  return 0.0;
2557}
2558void
2559CbcFixingBranchingObject::print(bool normalBranch)
2560{
2561  int i;
2562  // *** for way - up means fix all those in up section
2563  if (way_<0) {
2564    printf("Down Fix ");
2565    for (i=0;i<numberDown_;i++) {
2566      int iColumn = downList_[i];
2567      printf("%d ",iColumn);
2568    }
2569  } else {
2570    printf("Up Fix ");
2571    for (i=0;i<numberUp_;i++) {
2572      int iColumn = upList_[i];
2573      printf("%d ",iColumn);
2574    }
2575  }
2576  printf("\n");
2577}
Note: See TracBrowser for help on using the repository browser.